src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "gaussGrad.H"
  28 #include "zeroGradientFvPatchField.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace fv
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 template<class Type>
  43 tmp
  44 <
  45     GeometricField
  46     <
  47         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
  48     >
  49 >
  50 gaussGrad<Type>::grad
  51 (
  52     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
  53 )
  54 {
  55     typedef typename outerProduct<vector, Type>::type GradType;
  56
  57     const fvMesh& mesh = ssf.mesh();
  58
  59     tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad
  60     (
  61         new GeometricField<GradType, fvPatchField, volMesh>
  62         (
  63             IOobject
  64             (
  65                 "grad("+ssf.name()+')',
  66                 ssf.instance(),
  67                 mesh,
  68                 IOobject::NO_READ,
  69                 IOobject::NO_WRITE
  70             ),
  71             mesh,
  72             dimensioned<GradType>
  73             (
  74                 "0",
  75                 ssf.dimensions()/dimLength,
  76                 pTraits<GradType>::zero
  77             ),
  78             zeroGradientFvPatchField<GradType>::typeName
  79         )
  80     );
  81     GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();
  82
  83     const unallocLabelList& owner = mesh.owner();
  84     const unallocLabelList& neighbour = mesh.neighbour();
  85     const vectorField& Sf = mesh.Sf();
  86
  87     Field<GradType>& igGrad = gGrad;
  88     const Field<Type>& issf = ssf;
  89
  90     forAll(owner, facei)
  91     {
  92         GradType Sfssf = Sf[facei]*issf[facei];
  93
  94         igGrad[owner[facei]] += Sfssf;
  95         igGrad[neighbour[facei]] -= Sfssf;
  96     }
  97
  98     forAll(mesh.boundary(), patchi)
  99     {
 100         const unallocLabelList& pFaceCells =
 101             mesh.boundary()[patchi].faceCells();
 102
 103         const vectorField& pSf = mesh.Sf().boundaryField()[patchi];
 104
 105         const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
 106
 107         forAll(mesh.boundary()[patchi], facei)
 108         {
 109             igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
 110         }
 111     }
 112
 113     igGrad /= mesh.V();
 114
 115     gGrad.correctBoundaryConditions();
 116
 117     return tgGrad;
 118 }
 119
 120
 121 template<class Type>
 122 tmp
 123 <
 124     GeometricField
 125     <
 126         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
 127     >
 128 >
 129 gaussGrad<Type>::grad
 130 (
 131     const GeometricField<Type, fvPatchField, volMesh>& vsf
 132 ) const
 133 {
 134     typedef typename outerProduct<vector, Type>::type GradType;
 135
 136     tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad
 137     (
 138         grad(tinterpScheme_().interpolate(vsf))
 139     );
 140     GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();
 141
 142     gGrad.rename("grad(" + vsf.name() + ')');
 143     correctBoundaryConditions(vsf, gGrad);
 144
 145     return tgGrad;
 146 }
 147
 148
 149 template<class Type>
 150 void gaussGrad<Type>::correctBoundaryConditions
 151 (
 152     const GeometricField<Type, fvPatchField, volMesh>& vsf,
 153     GeometricField
 154     <
 155         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
 156     >& gGrad
 157 )
 158 {
 159     forAll(vsf.boundaryField(), patchi)
 160     {
 161         if (!vsf.boundaryField()[patchi].coupled())
 162         {
 163             vectorField n =
 164                 vsf.mesh().Sf().boundaryField()[patchi]
 165                /vsf.mesh().magSf().boundaryField()[patchi];
 166
 167             gGrad.boundaryField()[patchi] += n *
 168             (
 169                 vsf.boundaryField()[patchi].snGrad()
 170               - (n & gGrad.boundaryField()[patchi])
 171             );
 172         }
 173      }
 174 }
 175
 176
 177 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 178
 179 } // End namespace fv
 180
 181 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 182
 183 } // End namespace Foam
 184
 185 // ************************************************************************* //