src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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 \*---------------------------------------------------------------------------*/
  25
  26 #include "basicSymmetryFvPatchField.H"
  27 #include "volFields.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  35
  36 template<>
  37 tmp<scalarField > basicSymmetryFvPatchField<scalar>::snGrad() const
  38 {
  39     return tmp<scalarField>(new scalarField(size(), scalar(0)));
  40 }
  41
  42
  43 // Note:
  44 // Skew correction and second order is currently handled in template
  45 // specialisations because
  46 // Note 2: Currently repeated code in evaluate after typedefs.
  47 // Refactor.  HJ, 18/Mar/2015
  48 template<>
  49 void basicSymmetryFvPatchField<scalar>::evaluate(const Pstream::commsTypes)
  50 {
  51     // Local typedefs
  52     typedef scalar Type;
  53     typedef outerProduct<vector, Type>::type gradType;
  54     typedef GeometricField<gradType, fvPatchField, volMesh> gradFieldType;
  55
  56     if (!updated())
  57     {
  58         updateCoeffs();
  59     }
  60
  61     vectorField nHat = this->patch().nf();
  62
  63     // Get patch internal field
  64     Field<Type> pif = this->patchInternalField();
  65
  66     // Skew corrected treatment
  67     if (skewCorrected_)
  68     {
  69         // Access the gradient
  70         const word DName = this->dimensionedInternalField().name();
  71         const word gradDName("grad(" + DName + ")");
  72
  73         if (!this->db().foundObject<gradFieldType>(gradDName))
  74         {
  75             InfoIn
  76             (
  77                 "void basicSymmetryFvPatchField<Type>::"
  78                 "evaluate(const Pstream::commsTypes)"
  79             )   << "Cannot access " << gradDName << " field for field "
  80                 << DName << " for patch "<< this->patch().name()
  81                 << ".  Evaluating without skew correction"
  82                 << endl;
  83         }
  84         else
  85         {
  86             // Calculate skew correction vector k
  87             vectorField delta = patch().delta();
  88             vectorField k = delta - nHat*(nHat&delta);
  89
  90             // Get gradient
  91             const fvPatchField<gradType>& gradD =
  92                 patch().lookupPatchField<gradFieldType, gradType>(gradDName);
  93
  94             // Correct cell value for k
  95             pif += (k & gradD.patchInternalField());
  96
  97             if (secondOrder_)
  98             {
  99                 Field<Type> nGradD = (nHat & gradD.patchInternalField());
 100
 101                 Field<Type>::operator=
 102                 (
 103                     transform
 104                     (
 105                         I - sqr(nHat),
 106                         pif + 0.5*nGradD/this->patch().deltaCoeffs()
 107                     )
 108                 );
 109
 110                 // Return, skew corrected, second order
 111                 transformFvPatchField<Type>::evaluate();
 112                 return;
 113             }
 114             else
 115             {
 116                 Field<Type>::operator=
 117                 (
 118                     0.5*(pif + transform(I - 2.0*sqr(nHat), pif))
 119                 );
 120
 121                 // Return, skew corrected without second order correction
 122                 transformFvPatchField<Type>::evaluate();
 123                 return;
 124             }
 125         }
 126     }
 127
 128     // Without skew correction
 129     Field<Type>::operator=
 130     (
 131         0.5*(pif + transform(I - 2.0*sqr(nHat), pif))
 132     );
 133
 134     transformFvPatchField<Type>::evaluate();
 135 }
 136
 137
 138 template<>
 139 tmp<vectorField> basicSymmetryFvPatchField<vector>::snGrad() const
 140 {
 141     // Local typedefs
 142     typedef vector Type;
 143     typedef outerProduct<vector, Type>::type gradType;
 144     typedef GeometricField<gradType, fvPatchField, volMesh> gradFieldType;
 145
 146     vectorField nHat = this->patch().nf();
 147
 148     // Get patch internal field
 149     Field<Type> pif = this->patchInternalField();
 150
 151     // Skew corrected treatment
 152     if (skewCorrected_)
 153     {
 154         // Access the gradient
 155         const word DName = this->dimensionedInternalField().name();
 156         const word gradDName("grad(" + DName + ")");
 157
 158         if (!this->db().foundObject<gradFieldType>(gradDName))
 159         {
 160             InfoIn
 161             (
 162                 "void basicSymmetryFvPatchField<Type>::snGrad() const"
 163             )   << "Cannot access " << gradDName << " field for field "
 164                 << DName << " for patch "<< this->patch().name()
 165                 << ".  Evaluating without skew correction"
 166                 << endl;
 167         }
 168         else
 169         {
 170             // Calculate skew correction vector k
 171             vectorField delta = patch().delta();
 172             vectorField k = delta - nHat*(nHat&delta);
 173
 174             // Get gradient
 175             const fvPatchField<gradType>& gradD =
 176                 patch().lookupPatchField<gradFieldType, gradType>(gradDName);
 177
 178             // Correct cell value for k
 179             pif += (k & gradD.patchInternalField());
 180
 181             if (secondOrder_)
 182             {
 183                 Field<Type> nGradD = (nHat & gradD.patchInternalField());
 184
 185                 // Return, skew corrected, second order
 186                 return 2*
 187                     (
 188                         transform(I - 2.0*sqr(nHat), pif) - pif
 189                     )*(this->patch().deltaCoeffs()/2.0)
 190                   - transform(sqr(nHat), nGradD);
 191             }
 192             else
 193             {
 194                 // Return, skew corrected, without second order correction
 195                 return
 196                 (
 197                     transform(I - 2.0*sqr(nHat), pif) - pif
 198                 )*(this->patch().deltaCoeffs()/2.0);
 199             }
 200         }
 201     }
 202
 203     // Without skew correction
 204     return
 205     (
 206         transform(I - 2.0*sqr(nHat), pif) - pif
 207     )*(this->patch().deltaCoeffs()/2.0);
 208 }
 209
 210
 211 template<>
 212 void basicSymmetryFvPatchField<vector>::evaluate(const Pstream::commsTypes)
 213 {
 214     // Local typedefs
 215     typedef vector Type;
 216     typedef outerProduct<vector, Type>::type gradType;
 217     typedef GeometricField<gradType, fvPatchField, volMesh> gradFieldType;
 218
 219     if (!updated())
 220     {
 221         updateCoeffs();
 222     }
 223
 224     if (!updated())
 225     {
 226         updateCoeffs();
 227     }
 228
 229     vectorField nHat = this->patch().nf();
 230
 231     // Get patch internal field
 232     Field<Type> pif = this->patchInternalField();
 233
 234     // Skew corrected treatment
 235     if (skewCorrected_)
 236     {
 237         // Access the gradient
 238         const word DName = this->dimensionedInternalField().name();
 239         const word gradDName("grad(" + DName + ")");
 240
 241         if (!this->db().foundObject<gradFieldType>(gradDName))
 242         {
 243             InfoIn
 244             (
 245                 "void basicSymmetryFvPatchField<Type>::"
 246                 "evaluate(const Pstream::commsTypes)"
 247             )   << "Cannot access " << gradDName << " field for field "
 248                 << DName << " for patch "<< this->patch().name()
 249                 << ".  Evaluating without skew correction"
 250                 << endl;
 251         }
 252         else
 253         {
 254             // Calculate skew correction vector k
 255             vectorField delta = patch().delta();
 256             vectorField k = delta - nHat*(nHat&delta);
 257
 258             // Get gradient
 259             const fvPatchField<gradType>& gradD =
 260                 patch().lookupPatchField<gradFieldType, gradType>(gradDName);
 261
 262             // Correct cell value for k
 263             pif += (k & gradD.patchInternalField());
 264
 265             if (secondOrder_)
 266             {
 267                 Field<Type> nGradD = (nHat & gradD.patchInternalField());
 268
 269                 Field<Type>::operator=
 270                 (
 271                     transform
 272                     (
 273                         I - sqr(nHat),
 274                         pif + 0.5*nGradD/this->patch().deltaCoeffs()
 275                     )
 276                 );
 277
 278                 // Return, skew corrected, second order
 279                 transformFvPatchField<Type>::evaluate();
 280                 return;
 281             }
 282             else
 283             {
 284                 Field<Type>::operator=
 285                 (
 286                     0.5*(pif + transform(I - 2.0*sqr(nHat), pif))
 287                 );
 288
 289                 // Return, skew corrected without second order correction
 290                 transformFvPatchField<Type>::evaluate();
 291                 return;
 292             }
 293         }
 294     }
 295
 296     // Without skew correction
 297     Field<Type>::operator=
 298     (
 299         0.5*(pif + transform(I - 2.0*sqr(nHat), pif))
 300     );
 301
 302     transformFvPatchField<Type>::evaluate();
 303 }
 304
 305
 306 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 307
 308 } // End namespace Foam
 309
 310 // ************************************************************************* //