src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.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 "leastSquaresGrad.H"
  28 #include "leastSquaresVectors.H"
  29 #include "gaussGrad.H"
  30 #include "fvMesh.H"
  31 #include "volMesh.H"
  32 #include "surfaceMesh.H"
  33 #include "GeometricField.H"
  34 #include "zeroGradientFvPatchField.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace fv
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 template<class Type>
  49 tmp
  50 <
  51     GeometricField
  52     <
  53         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
  54     >
  55 >
  56 leastSquaresGrad<Type>::grad
  57 (
  58     const GeometricField<Type, fvPatchField, volMesh>& vsf
  59 ) const
  60 {
  61     typedef typename outerProduct<vector, Type>::type GradType;
  62
  63     const fvMesh& mesh = vsf.mesh();
  64
  65     tmp<GeometricField<GradType, fvPatchField, volMesh> > tlsGrad
  66     (
  67         new GeometricField<GradType, fvPatchField, volMesh>
  68         (
  69             IOobject
  70             (
  71                 "grad("+vsf.name()+')',
  72                 vsf.instance(),
  73                 mesh,
  74                 IOobject::NO_READ,
  75                 IOobject::NO_WRITE
  76             ),
  77             mesh,
  78             dimensioned<GradType>
  79             (
  80                 "zero",
  81                 vsf.dimensions()/dimLength,
  82                 pTraits<GradType>::zero
  83             ),
  84             zeroGradientFvPatchField<GradType>::typeName
  85         )
  86     );
  87     GeometricField<GradType, fvPatchField, volMesh>& lsGrad = tlsGrad();
  88
  89     // Get reference to least square vectors
  90     const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
  91
  92     const surfaceVectorField& ownLs = lsv.pVectors();
  93     const surfaceVectorField& neiLs = lsv.nVectors();
  94
  95     const unallocLabelList& own = mesh.owner();
  96     const unallocLabelList& nei = mesh.neighbour();
  97
  98     forAll(own, facei)
  99     {
 100         register label ownFaceI = own[facei];
 101         register label neiFaceI = nei[facei];
 102
 103         Type deltaVsf = vsf[neiFaceI] - vsf[ownFaceI];
 104
 105         lsGrad[ownFaceI] += ownLs[facei]*deltaVsf;
 106         lsGrad[neiFaceI] -= neiLs[facei]*deltaVsf;
 107     }
 108
 109     // Boundary faces
 110     forAll(vsf.boundaryField(), patchi)
 111     {
 112         const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
 113
 114         const unallocLabelList& faceCells =
 115             lsGrad.boundaryField()[patchi].patch().faceCells();
 116
 117         if (vsf.boundaryField()[patchi].coupled())
 118         {
 119             Field<Type> neiVsf =
 120                 vsf.boundaryField()[patchi].patchNeighbourField();
 121
 122             forAll(neiVsf, patchFaceI)
 123             {
 124                 lsGrad[faceCells[patchFaceI]] +=
 125                     patchOwnLs[patchFaceI]
 126                    *(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
 127             }
 128         }
 129         else
 130         {
 131             const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
 132
 133             forAll(patchVsf, patchFaceI)
 134             {
 135                 lsGrad[faceCells[patchFaceI]] +=
 136                      patchOwnLs[patchFaceI]
 137                     *(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
 138             }
 139         }
 140     }
 141
 142
 143     lsGrad.correctBoundaryConditions();
 144     gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);
 145
 146     return tlsGrad;
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 151
 152 } // End namespace fv
 153
 154 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 155
 156 } // End namespace Foam
 157
 158 // ************************************************************************* //