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