src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 "leastSquaresVectors.H"
  28 #include "surfaceFields.H"
  29 #include "volFields.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35     defineTypeNameAndDebug(leastSquaresVectors, 0);
  36 }
  37
  38
  39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  40
  41 Foam::leastSquaresVectors::leastSquaresVectors(const fvMesh& mesh)
  42 :
  43     MeshObject<fvMesh, leastSquaresVectors>(mesh),
  44     pVectorsPtr_(NULL),
  45     nVectorsPtr_(NULL)
  46 {}
  47
  48
  49 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
  50
  51 Foam::leastSquaresVectors::~leastSquaresVectors()
  52 {
  53     deleteDemandDrivenData(pVectorsPtr_);
  54     deleteDemandDrivenData(nVectorsPtr_);
  55 }
  56
  57
  58 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  59
  60 void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
  61 {
  62     if (debug)
  63     {
  64         Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
  65             << "Constructing least square gradient vectors"
  66             << endl;
  67     }
  68
  69     const fvMesh& mesh = mesh_;
  70
  71     pVectorsPtr_ = new surfaceVectorField
  72     (
  73         IOobject
  74         (
  75             "LeastSquaresP",
  76             mesh_.pointsInstance(),
  77             mesh_,
  78             IOobject::NO_READ,
  79             IOobject::NO_WRITE,
  80             false
  81         ),
  82         mesh_,
  83         dimensionedVector("zero", dimless/dimLength, vector::zero)
  84     );
  85     surfaceVectorField& lsP = *pVectorsPtr_;
  86
  87     nVectorsPtr_ = new surfaceVectorField
  88     (
  89         IOobject
  90         (
  91             "LeastSquaresN",
  92             mesh_.pointsInstance(),
  93             mesh_,
  94             IOobject::NO_READ,
  95             IOobject::NO_WRITE,
  96             false
  97         ),
  98         mesh_,
  99         dimensionedVector("zero", dimless/dimLength, vector::zero)
 100     );
 101     surfaceVectorField& lsN = *nVectorsPtr_;
 102
 103     // Set local references to mesh data
 104     const unallocLabelList& owner = mesh_.owner();
 105     const unallocLabelList& neighbour = mesh_.neighbour();
 106
 107     const volVectorField& C = mesh.C();
 108     const surfaceScalarField& w = mesh.weights();
 109     const surfaceScalarField& magSf = mesh.magSf();
 110
 111
 112     // Set up temporary storage for the dd tensor (before inversion)
 113     symmTensorField dd(mesh_.nCells(), symmTensor::zero);
 114
 115     forAll(owner, facei)
 116     {
 117         label own = owner[facei];
 118         label nei = neighbour[facei];
 119
 120         vector d = C[nei] - C[own];
 121         symmTensor wdd = (magSf[facei]/magSqr(d))*sqr(d);
 122
 123         dd[own] += (1 - w[facei])*wdd;
 124         dd[nei] += w[facei]*wdd;
 125     }
 126
 127
 128     forAll(lsP.boundaryField(), patchi)
 129     {
 130         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
 131         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
 132
 133         const fvPatch& p = pw.patch();
 134         const unallocLabelList& faceCells = p.patch().faceCells();
 135
 136         // Build the d-vectors
 137         vectorField pd =
 138             mesh.Sf().boundaryField()[patchi]
 139            /(
 140                mesh.magSf().boundaryField()[patchi]
 141               *mesh.deltaCoeffs().boundaryField()[patchi]
 142            );
 143
 144         if (!mesh.orthogonal())
 145         {
 146             pd -= mesh.correctionVectors().boundaryField()[patchi]
 147                 /mesh.deltaCoeffs().boundaryField()[patchi];
 148         }
 149
 150         if (p.coupled())
 151         {
 152             forAll(pd, patchFacei)
 153             {
 154                 const vector& d = pd[patchFacei];
 155
 156                 dd[faceCells[patchFacei]] +=
 157                     ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))*sqr(d);
 158             }
 159         }
 160         else
 161         {
 162             forAll(pd, patchFacei)
 163             {
 164                 const vector& d = pd[patchFacei];
 165
 166                 dd[faceCells[patchFacei]] +=
 167                     (pMagSf[patchFacei]/magSqr(d))*sqr(d);
 168             }
 169         }
 170     }
 171
 172
 173     // Invert the dd tensor
 174     symmTensorField invDd = inv(dd);
 175
 176
 177     // Revisit all faces and calculate the lsP and lsN vectors
 178     forAll(owner, facei)
 179     {
 180         label own = owner[facei];
 181         label nei = neighbour[facei];
 182
 183         vector d = C[nei] - C[own];
 184         scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
 185
 186         lsP[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
 187         lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
 188     }
 189
 190     forAll(lsP.boundaryField(), patchi)
 191     {
 192         fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi];
 193
 194         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
 195         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
 196
 197         const fvPatch& p = pw.patch();
 198         const unallocLabelList& faceCells = p.faceCells();
 199
 200         // Build the d-vectors
 201         vectorField pd =
 202             mesh.Sf().boundaryField()[patchi]
 203            /(
 204                mesh.magSf().boundaryField()[patchi]
 205               *mesh.deltaCoeffs().boundaryField()[patchi]
 206            );
 207
 208         if (!mesh.orthogonal())
 209         {
 210             pd -= mesh.correctionVectors().boundaryField()[patchi]
 211                 /mesh.deltaCoeffs().boundaryField()[patchi];
 212         }
 213
 214
 215         if (p.coupled())
 216         {
 217             forAll(pd, patchFacei)
 218             {
 219                 const vector& d = pd[patchFacei];
 220
 221                 patchLsP[patchFacei] =
 222                     ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))
 223                    *(invDd[faceCells[patchFacei]] & d);
 224             }
 225         }
 226         else
 227         {
 228             forAll(pd, patchFacei)
 229             {
 230                 const vector& d = pd[patchFacei];
 231
 232                 patchLsP[patchFacei] =
 233                     pMagSf[patchFacei]*(1.0/magSqr(d))
 234                    *(invDd[faceCells[patchFacei]] & d);
 235             }
 236         }
 237     }
 238
 239
 240     // For 3D meshes check the determinant of the dd tensor and switch to
 241     // Gauss if it is less than 3
 242     /* Currently the det(dd[celli]) criterion is incorrect: dd is weighted by Sf
 243     if (mesh.nGeometricD() == 3)
 244     {
 245         label nBadCells = 0;
 246
 247         const cellList& cells = mesh.cells();
 248         const scalarField& V = mesh.V();
 249         const surfaceVectorField& Sf = mesh.Sf();
 250         const surfaceScalarField& w = mesh.weights();
 251
 252         forAll (dd, celli)
 253         {
 254             if (det(dd[celli]) < 3)
 255             {
 256                 nBadCells++;
 257
 258                 const cell& c = cells[celli];
 259
 260                 forAll(c, cellFacei)
 261                 {
 262                     label facei = c[cellFacei];
 263
 264                     if (mesh.isInternalFace(facei))
 265                     {
 266                         scalar wf = max(min(w[facei], 0.8), 0.2);
 267
 268                         if (celli == owner[facei])
 269                         {
 270                             lsP[facei] = (1 - wf)*Sf[facei]/V[celli];
 271                         }
 272                         else
 273                         {
 274                             lsN[facei] = -wf*Sf[facei]/V[celli];
 275                         }
 276                     }
 277                     else
 278                     {
 279                         label patchi = mesh.boundaryMesh().whichPatch(facei);
 280
 281                         if (mesh.boundary()[patchi].size())
 282                         {
 283                             label patchFacei =
 284                                 facei - mesh.boundaryMesh()[patchi].start();
 285
 286                             if (mesh.boundary()[patchi].coupled())
 287                             {
 288                                 scalar wf = max
 289                                 (
 290                                     min
 291                                     (
 292                                         w.boundaryField()[patchi][patchFacei],
 293                                         0.8
 294                                     ),
 295                                     0.2
 296                                 );
 297
 298                                 lsP.boundaryField()[patchi][patchFacei] =
 299                                     (1 - wf)
 300                                    *Sf.boundaryField()[patchi][patchFacei]
 301                                    /V[celli];
 302                             }
 303                             else
 304                             {
 305                                 lsP.boundaryField()[patchi][patchFacei] =
 306                                     Sf.boundaryField()[patchi][patchFacei]
 307                                    /V[celli];
 308                             }
 309                         }
 310                     }
 311                 }
 312             }
 313         }
 314
 315         if (debug)
 316         {
 317             InfoIn("leastSquaresVectors::makeLeastSquaresVectors()")
 318                 << "number of bad cells switched to Gauss = " << nBadCells
 319                 << endl;
 320         }
 321     }
 322     */
 323
 324     if (debug)
 325     {
 326         Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
 327             << "Finished constructing least square gradient vectors"
 328             << endl;
 329     }
 330 }
 331
 332
 333 const Foam::surfaceVectorField& Foam::leastSquaresVectors::pVectors() const
 334 {
 335     if (!pVectorsPtr_)
 336     {
 337         makeLeastSquaresVectors();
 338     }
 339
 340     return *pVectorsPtr_;
 341 }
 342
 343
 344 const Foam::surfaceVectorField& Foam::leastSquaresVectors::nVectors() const
 345 {
 346     if (!nVectorsPtr_)
 347     {
 348         makeLeastSquaresVectors();
 349     }
 350
 351     return *nVectorsPtr_;
 352 }
 353
 354
 355 bool Foam::leastSquaresVectors::movePoints()
 356 {
 357     deleteDemandDrivenData(pVectorsPtr_);
 358     deleteDemandDrivenData(nVectorsPtr_);
 359
 360     return true;
 361 }
 362
 363
 364 // ************************************************************************* //