src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Class to create the weighting-factors based on the NVD
  26     (Normalised Variable Diagram).
  27     The particular differencing scheme class is supplied as a template argument,
  28     the weight function of which is called by the weight function of this class
  29     for the internal edges as well as edges of coupled patches
  30     (e.g. processor-processor patches). The weight function is supplied the
  31     central-differencing weighting factor, the edge-flux, the cell and edge
  32     gradients (from which the normalised variable distribution may be created)
  33     and the cell centre distance.
  34
  35     This code organisation is both neat and efficient, allowing for convenient
  36     implementation of new schemes to run on parallelised cases.
  37
  38 \*---------------------------------------------------------------------------*/
  39
  40 #include "areaFields.H"
  41 #include "edgeFields.H"
  42 #include "facGrad.H"
  43 #include "coupledFaPatchFields.H"
  44 #include "upwindEdgeInterpolation.H"
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 namespace Foam
  49 {
  50
  51 inline tmp<areaScalarField> limiter(const areaScalarField& phi)
  52 {
  53     return phi;
  54 }
  55
  56 inline tmp<areaScalarField> limiter(const areaVectorField& phi)
  57 {
  58     return magSqr(phi);
  59 }
  60
  61 inline tmp<areaScalarField> limiter(const areaTensorField& phi)
  62 {
  63     return magSqr(phi);
  64 }
  65
  66
  67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  68 //- Return the interpolation weighting factors
  69
  70 template<class Type, class NVDweight>
  71 tmp<edgeScalarField> faNVDscheme<Type,NVDweight>::weights
  72 (
  73     const GeometricField<Type, faPatchField, areaMesh>& phi
  74 ) const
  75 {
  76     const faMesh& mesh = this->mesh();
  77
  78     tmp<edgeScalarField> tWeightingFactors
  79     (
  80         new edgeScalarField(mesh.edgeInterpolation::weights())
  81     );
  82     edgeScalarField& weightingFactors = tWeightingFactors();
  83
  84     scalarField& weights = weightingFactors.internalField();
  85
  86     tmp<areaScalarField> tvf = limiter(phi);
  87     const areaScalarField& vf = tvf();
  88
  89     areaVectorField gradc(fac::grad(vf));
  90
  91 //     edgeVectorField d =
  92 //         mesh.Le()
  93 //        /(mesh.magLe()*mesh.edgeInterpolation::deltaCoeffs());
  94
  95 //     if (!mesh.orthogonal())
  96 //     {
  97 //         d -=
  98 //             mesh.edgeInterpolation::correctionVectors()
  99 //            /mesh.edgeInterpolation::deltaCoeffs();
 100 //     }
 101
 102     const unallocLabelList& owner = mesh.owner();
 103     const unallocLabelList& neighbour = mesh.neighbour();
 104     const vectorField& n = mesh.faceAreaNormals().internalField();
 105     const vectorField& c = mesh.areaCentres().internalField();
 106
 107     forAll(weights, edge)
 108     {
 109         vector d = vector::zero;
 110
 111         if(edgeFlux_[edge] > 0)
 112         {
 113             d = c[neighbour[edge]] - c[owner[edge]];
 114             d -= n[owner[edge]]*(n[owner[edge]]&d);
 115             d /= mag(d)/mesh.edgeInterpolation::lPN().internalField()[edge];
 116         }
 117         else
 118         {
 119             d = c[neighbour[edge]] - c[owner[edge]];
 120             d -= n[neighbour[edge]]*(n[neighbour[edge]]&d);
 121             d /= mag(d)/mesh.edgeInterpolation::lPN().internalField()[edge];
 122         }
 123
 124         weights[edge] =
 125             this->weight
 126             (
 127                 weights[edge],
 128                 edgeFlux_[edge],
 129                 vf[owner[edge]],
 130                 vf[neighbour[edge]],
 131                 gradc[owner[edge]],
 132                 gradc[neighbour[edge]],
 133                 d
 134             );
 135     }
 136
 137
 138     GeometricField<scalar, faePatchField, edgeMesh>::GeometricBoundaryField&
 139         bWeights = weightingFactors.boundaryField();
 140
 141     forAll(bWeights, patchI)
 142     {
 143         if (bWeights[patchI].coupled())
 144         {
 145             scalarField& pWeights = bWeights[patchI];
 146
 147             const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchI];
 148
 149             scalarField pVfP =
 150                 vf.boundaryField()[patchI].patchInternalField();
 151
 152             scalarField pVfN =
 153                 vf.boundaryField()[patchI].patchNeighbourField();
 154
 155             vectorField pGradcP =
 156                 gradc.boundaryField()[patchI].patchInternalField();
 157
 158             vectorField pGradcN =
 159                 gradc.boundaryField()[patchI].patchNeighbourField();
 160
 161             vectorField CP = mesh.areaCentres().boundaryField()[patchI]
 162                 .patchInternalField();
 163
 164             vectorField CN =
 165                 mesh.areaCentres().boundaryField()[patchI]
 166                 .patchNeighbourField();
 167
 168             vectorField nP =
 169                 mesh.faceAreaNormals().boundaryField()[patchI]
 170                .patchInternalField();
 171
 172             vectorField nN =
 173                 mesh.faceAreaNormals().boundaryField()[patchI]
 174                .patchNeighbourField();
 175
 176             scalarField pLPN =
 177                 mesh.edgeInterpolation::lPN().boundaryField()[patchI];
 178
 179             forAll(pWeights, edgeI)
 180             {
 181                 vector d = vector::zero;
 182
 183                 if(pEdgeFlux[edgeI] > 0)
 184                 {
 185                     d = CN[edgeI] - CP[edgeI];
 186                     d -= nP[edgeI]*(nP[edgeI]&d);
 187                     d /= mag(d)/pLPN[edgeI];
 188                 }
 189                 else
 190                 {
 191                     d = CN[edgeI] - CP[edgeI];
 192                     d -= nN[edgeI]*(nN[edgeI]&d);
 193                     d /= mag(d)/pLPN[edgeI];
 194                 }
 195
 196                 pWeights[edgeI] =
 197                     this->weight
 198                     (
 199                         pWeights[edgeI],
 200                         pEdgeFlux[edgeI],
 201                         pVfP[edgeI],
 202                         pVfN[edgeI],
 203                         pGradcP[edgeI],
 204                         pGradcN[edgeI],
 205                         d
 206                     );
 207             }
 208         }
 209     }
 210
 211     return tWeightingFactors;
 212 }
 213
 214
 215 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 216
 217 } // End namespace Foam
 218
 219 // ************************************************************************* //