src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.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 "extendedLeastSquaresGrad.H"
  28 #include "extendedLeastSquaresVectors.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 extendedLeastSquaresGrad<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 extendedLeastSquaresVectors& lsv = extendedLeastSquaresVectors::New
  91     (
  92         mesh,
  93         minDet_
  94     );
  95
  96     const surfaceVectorField& ownLs = lsv.pVectors();
  97     const surfaceVectorField& neiLs = lsv.nVectors();
  98
  99     const unallocLabelList& owner = mesh.owner();
 100     const unallocLabelList& neighbour = mesh.neighbour();
 101
 102     forAll(owner, facei)
 103     {
 104         label own = owner[facei];
 105         label nei = neighbour[facei];
 106
 107         Type deltaVsf = vsf[nei] - vsf[own];
 108
 109         lsGrad[own] += ownLs[facei]*deltaVsf;
 110         lsGrad[nei] -= neiLs[facei]*deltaVsf;
 111     }
 112
 113     // Boundary faces
 114     forAll(vsf.boundaryField(), patchi)
 115     {
 116         const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
 117
 118         const unallocLabelList& faceCells =
 119             lsGrad.boundaryField()[patchi].patch().faceCells();
 120
 121         if (vsf.boundaryField()[patchi].coupled())
 122         {
 123             Field<Type> neiVsf =
 124                 vsf.boundaryField()[patchi].patchNeighbourField();
 125
 126             forAll(neiVsf, patchFaceI)
 127             {
 128                 lsGrad[faceCells[patchFaceI]] +=
 129                     patchOwnLs[patchFaceI]
 130                    *(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
 131             }
 132         }
 133         else
 134         {
 135             const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
 136
 137             forAll(patchVsf, patchFaceI)
 138             {
 139                 lsGrad[faceCells[patchFaceI]] +=
 140                      patchOwnLs[patchFaceI]
 141                     *(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
 142             }
 143         }
 144     }
 145
 146
 147     const List<labelPair>& additionalCells = lsv.additionalCells();
 148     const vectorField& additionalVectors = lsv.additionalVectors();
 149
 150     forAll(additionalCells, i)
 151     {
 152         lsGrad[additionalCells[i][0]] +=
 153             additionalVectors[i]
 154            *(vsf[additionalCells[i][1]] - vsf[additionalCells[i][0]]);
 155     }
 156
 157
 158     lsGrad.correctBoundaryConditions();
 159     gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);
 160
 161     return tlsGrad;
 162 }
 163
 164
 165 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 166
 167 } // End namespace fv
 168
 169 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 170
 171 } // End namespace Foam
 172
 173 // ************************************************************************* //