src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "linearUpwind.H"
  28 #include "fvMesh.H"
  29 #include "zeroGradientFvPatchField.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 template<class Type>
  34 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
  35 Foam::linearUpwind<Type>::correction
  36 (
  37     const GeometricField<Type, fvPatchField, volMesh>& vf
  38 ) const
  39 {
  40     const fvMesh& mesh = this->mesh();
  41
  42     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
  43     (
  44         new GeometricField<Type, fvsPatchField, surfaceMesh>
  45         (
  46             IOobject
  47             (
  48                 "linearUpwindCorrection(" + vf.name() + ')',
  49                 mesh.time().timeName(),
  50                 mesh,
  51                 IOobject::NO_READ,
  52                 IOobject::NO_WRITE,
  53                 false
  54             ),
  55             mesh,
  56             dimensioned<Type>(vf.name(), vf.dimensions(), pTraits<Type>::zero)
  57         )
  58     );
  59
  60     GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr();
  61
  62     const surfaceScalarField& faceFlux = this->faceFlux_;
  63
  64     const labelList& owner = mesh.owner();
  65     const labelList& neighbour = mesh.neighbour();
  66
  67     const volVectorField& C = mesh.C();
  68     const surfaceVectorField& Cf = mesh.Cf();
  69
  70     GeometricField
  71         <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
  72         gradVf = gradScheme_().grad(vf);
  73
  74     // Note: in order for the patchNeighbourField to be correct on coupled
  75     // boundaries, correctBoundaryConditions needs to be called.
  76     // The call shall be moved into the library fvc operators
  77     gradVf.correctBoundaryConditions();
  78
  79     forAll(faceFlux, facei)
  80     {
  81         if (faceFlux[facei] > 0)
  82         {
  83             label own = owner[facei];
  84             sfCorr[facei] = (Cf[facei] - C[own]) & gradVf[own];
  85         }
  86         else
  87         {
  88             label nei = neighbour[facei];
  89             sfCorr[facei] = (Cf[facei] - C[nei]) & gradVf[nei];
  90         }
  91     }
  92
  93
  94     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
  95         GeometricBoundaryField& bSfCorr = sfCorr.boundaryField();
  96
  97     forAll(bSfCorr, patchi)
  98     {
  99         fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
 100
 101         if (pSfCorr.coupled())
 102         {
 103             const fvPatch& p = mesh.boundary()[patchi];
 104
 105             const unallocLabelList& pOwner = p.faceCells();
 106
 107             const vectorField& pCf = Cf.boundaryField()[patchi];
 108
 109             const scalarField& pFaceFlux = faceFlux.boundaryField()[patchi];
 110
 111             Field<typename outerProduct<vector, Type>::type> pGradVfNei =
 112                 gradVf.boundaryField()[patchi].patchNeighbourField();
 113
 114             // Build the d-vectors
 115             // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
 116             vectorField pd = p.delta();
 117
 118             forAll(pOwner, facei)
 119             {
 120                 label own = pOwner[facei];
 121
 122                 if (pFaceFlux[facei] > 0)
 123                 {
 124                     pSfCorr[facei] = (pCf[facei] - C[own]) & gradVf[own];
 125                 }
 126                 else
 127                 {
 128                     pSfCorr[facei] =
 129                         (pCf[facei] - pd[facei] - C[own]) & pGradVfNei[facei];
 130                 }
 131             }
 132         }
 133     }
 134
 135     return tsfCorr;
 136 }
 137
 138
 139 namespace Foam
 140 {
 141     //makelimitedSurfaceInterpolationScheme(linearUpwind)
 142     makelimitedSurfaceInterpolationTypeScheme(linearUpwind, scalar)
 143     makelimitedSurfaceInterpolationTypeScheme(linearUpwind, vector)
 144 }
 145
 146 // ************************************************************************* //