src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/reconCentral/reconCentral.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 "reconCentral.H"
  28 #include "fvMesh.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 template<class Type>
  33 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
  34 Foam::reconCentral<Type>::interpolate
  35 (
  36     const GeometricField<Type, fvPatchField, volMesh>& vf
  37 ) const
  38 {
  39     const fvMesh& mesh = this->mesh();
  40
  41     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsf
  42     (
  43         new GeometricField<Type, fvsPatchField, surfaceMesh>
  44         (
  45             IOobject
  46             (
  47                 "interpolate("+vf.name()+')',
  48                 vf.instance(),
  49                 vf.db()
  50             ),
  51             mesh,
  52             dimensioned<Type>(vf.name(), vf.dimensions(), pTraits<Type>::zero)
  53         )
  54     );
  55
  56     GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf();
  57
  58     const labelList& owner = mesh.owner();
  59     const labelList& neighbour = mesh.neighbour();
  60
  61     const volVectorField& C = mesh.C();
  62     const surfaceVectorField& Cf = mesh.Cf();
  63
  64     GeometricField
  65         <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
  66         gradVf = gradScheme_().grad(vf);
  67
  68     // Note: in order for the patchNeighbourField to be correct on coupled
  69     // boundaries, correctBoundaryConditions needs to be called.
  70     // The call shall be moved into the library fvc operators
  71     gradVf.correctBoundaryConditions();
  72
  73     Field<Type>& sfIn = sf.internalField();
  74
  75     forAll(sfIn, facei)
  76     {
  77         // Owner contribution
  78         label own = owner[facei];
  79         sfIn[facei] += 0.5*(vf[own] + ((Cf[facei] - C[own]) & gradVf[own]));
  80
  81         // Neighbour contribution
  82         label nei = neighbour[facei];
  83         sfIn[facei] += 0.5*(vf[nei] + ((Cf[facei] - C[nei]) & gradVf[nei]));
  84     }
  85
  86
  87     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
  88         GeometricBoundaryField& bSf = sf.boundaryField();
  89
  90     forAll(bSf, patchi)
  91     {
  92         const fvPatch& p = mesh.boundary()[patchi];
  93
  94         fvsPatchField<Type>& pSf = bSf[patchi];
  95
  96         const unallocLabelList& pOwner = p.faceCells();
  97
  98         const vectorField& pCf = Cf.boundaryField()[patchi];
  99
 100         if (pSf.coupled())
 101         {
 102             Field<Type> vfNei =
 103                 vf.boundaryField()[patchi].patchNeighbourField();
 104
 105             Field<typename outerProduct<vector, Type>::type> pGradVfNei =
 106                 gradVf.boundaryField()[patchi].patchNeighbourField();
 107
 108             // Build the d-vectors.  Used to calculate neighbour face centre
 109             // HJ, 19/Apr/2010
 110             // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
 111             vectorField pd = p.delta();
 112
 113             forAll(pOwner, facei)
 114             {
 115                 label own = pOwner[facei];
 116
 117                 // Owner contribution
 118                 pSf[facei] +=
 119                     0.5*(vf[own] + ((pCf[facei] - C[own]) & gradVf[own]));
 120
 121                 // Neighbour contribution
 122                 pSf[facei] +=
 123                     0.5*
 124                     (
 125                         vfNei[facei]
 126                       + (
 127                             (pCf[facei] - pd[facei] - C[own])
 128                           & pGradVfNei[facei]
 129                         )
 130                     );
 131             }
 132         }
 133         else if (vf.boundaryField()[patchi].fixesValue())
 134         {
 135             // For fixed boundary patches copy the value
 136             pSf = vf.boundaryField()[patchi];
 137         }
 138         else
 139         {
 140             // For patches that do not fix the value, calculate
 141             // extrapolated field
 142             forAll(pOwner, facei)
 143             {
 144                 label own = pOwner[facei];
 145
 146                 pSf[facei] =
 147                     (vf[own] + ((pCf[facei] - C[own]) & gradVf[own]));
 148             }
 149         }
 150     }
 151
 152     return tsf;
 153 }
 154
 155
 156 template<class Type>
 157 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 158 Foam::reconCentral<Type>::correction
 159 (
 160     const GeometricField<Type, fvPatchField, volMesh>& vf
 161 ) const
 162 {
 163     // Note: Correction is calculated by assembling the complete interpolation
 164     // including extrapolated gradient contribution and subtracting the
 165     // implicit contribution.  HJ, 27/Mar/2010
 166     const fvMesh& mesh = this->mesh();
 167
 168     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
 169     (
 170         new GeometricField<Type, fvsPatchField, surfaceMesh>
 171         (
 172             IOobject
 173             (
 174                 "reconCentralCorrection(" + vf.name() + ')',
 175                 mesh.time().timeName(),
 176                 mesh,
 177                 IOobject::NO_READ,
 178                 IOobject::NO_WRITE,
 179                 false
 180             ),
 181             reconCentral<Type>::interpolate(vf)
 182           - surfaceInterpolationScheme<Type>::interpolate
 183             (
 184                 vf,
 185                 this->weights()
 186             )
 187         )
 188     );
 189
 190     return tsfCorr;
 191 }
 192
 193
 194 namespace Foam
 195 {
 196     //makelimitedSurfaceInterpolationScheme(reconCentral)
 197     makelimitedSurfaceInterpolationTypeScheme(reconCentral, scalar)
 198     makelimitedSurfaceInterpolationTypeScheme(reconCentral, vector)
 199 }
 200
 201 // ************************************************************************* //