src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Class
  25     Foam::cubic
  26
  27 Description
  28     Cubic interpolation scheme class derived from linear and returns
  29     linear weighting factors but also applies an explicit correction.
  30
  31 SourceFiles
  32     cubic.C
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #ifndef cubic_H
  37 #define cubic_H
  38
  39 #include "linear.H"
  40 #include "gaussGrad.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 /*---------------------------------------------------------------------------*\
  48                            Class cubic Declaration
  49 \*---------------------------------------------------------------------------*/
  50
  51 template<class Type>
  52 class cubic
  53 :
  54     public linear<Type>
  55 {
  56     // Private Member Functions
  57
  58         //- Disallow default bitwise copy construct
  59         cubic(const cubic&);
  60
  61         //- Disallow default bitwise assignment
  62         void operator=(const cubic&);
  63
  64
  65 public:
  66
  67     //- Runtime type information
  68     TypeName("cubic");
  69
  70
  71     // Constructors
  72
  73         //- Construct from mesh
  74         cubic(const fvMesh& mesh)
  75         :
  76             linear<Type>(mesh)
  77         {}
  78
  79         //- Construct from mesh and Istream
  80         cubic
  81         (
  82             const fvMesh& mesh,
  83             Istream&
  84         )
  85         :
  86             linear<Type>(mesh)
  87         {}
  88
  89         //- Construct from mesh, faceFlux and Istream
  90         cubic
  91         (
  92             const fvMesh& mesh,
  93             const surfaceScalarField&,
  94             Istream&
  95         )
  96         :
  97             linear<Type>(mesh)
  98         {}
  99
 100
 101     // Member Functions
 102
 103         //- Return true if this scheme uses an explicit correction
 104         virtual bool corrected() const
 105         {
 106             return true;
 107         }
 108
 109         //- Return the explicit correction to the face-interpolate
 110         virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
 111         correction
 112         (
 113             const GeometricField<Type, fvPatchField, volMesh>& vf
 114         ) const
 115         {
 116             const fvMesh& mesh = this->mesh();
 117
 118             // calculate the appropriate interpolation factors
 119             const surfaceScalarField& lambda = mesh.weights();
 120
 121             const surfaceScalarField kSc
 122             (
 123                 lambda*(scalar(1) - lambda*(scalar(3) - scalar(2)*lambda))
 124             );
 125
 126             const surfaceScalarField kVecP(sqr(scalar(1) - lambda)*lambda);
 127             const surfaceScalarField kVecN(sqr(lambda)*(lambda - scalar(1)));
 128
 129             tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
 130             (
 131                 new GeometricField<Type, fvsPatchField, surfaceMesh>
 132                 (
 133                     IOobject
 134                     (
 135                         "cubic::correction(" + vf.name() +')',
 136                         mesh.time().timeName(),
 137                         mesh,
 138                         IOobject::NO_READ,
 139                         IOobject::NO_WRITE,
 140                         false
 141                     ),
 142                     surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc)
 143                 )
 144             );
 145
 146             GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr =
 147                 tsfCorr();
 148
 149             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
 150             {
 151                 sfCorr.replace
 152                 (
 153                     cmpt,
 154                     sfCorr.component(cmpt)
 155                   + (
 156                         surfaceInterpolationScheme
 157                         <
 158                             typename outerProduct
 159                             <
 160                                 vector,
 161                                 typename pTraits<Type>::cmptType
 162                             >::type
 163                         >::interpolate
 164                         (
 165                             fv::gaussGrad
 166                             <typename pTraits<Type>::cmptType>(mesh)
 167                            .grad(vf.component(cmpt)),
 168                             kVecP,
 169                             kVecN
 170                         ) & mesh.Sf()
 171                     )/mesh.magSf()/mesh.surfaceInterpolation::deltaCoeffs()
 172                 );
 173             }
 174
 175             forAll(sfCorr.boundaryField(), pi)
 176             {
 177                 if (!sfCorr.boundaryField()[pi].coupled())
 178                 {
 179                     sfCorr.boundaryField()[pi] = pTraits<Type>::zero;
 180                 }
 181             }
 182
 183             return tsfCorr;
 184         }
 185 };
 186
 187
 188 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 189
 190 } // End namespace Foam
 191
 192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 193
 194 #endif
 195
 196 // ************************************************************************* //