src/finiteVolume/fvMesh/fvMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "fvMesh.H"
  28 #include "Time.H"
  29 #include "volFields.H"
  30 #include "surfaceFields.H"
  31 #include "slicedVolFields.H"
  32 #include "slicedSurfaceFields.H"
  33 #include "SubField.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void fvMesh::makeSf() const
  43 {
  44     if (debug)
  45     {
  46         Info<< "void fvMesh::makeSf() : "
  47             << "assembling face areas"
  48             << endl;
  49     }
  50
  51     // It is an error to attempt to recalculate
  52     // if the pointer is already set
  53     if (SfPtr_)
  54     {
  55         FatalErrorIn("fvMesh::makeSf()")
  56             << "face areas already exist"
  57             << abort(FatalError);
  58     }
  59
  60     SfPtr_ = new slicedSurfaceVectorField
  61     (
  62         IOobject
  63         (
  64             "S",
  65             pointsInstance(),
  66             meshSubDir,
  67             *this
  68         ),
  69         *this,
  70         dimArea,
  71         faceAreas()
  72     );
  73 }
  74
  75
  76 void fvMesh::makeMagSf() const
  77 {
  78     if (debug)
  79     {
  80         Info<< "void fvMesh::makeMagSf() : "
  81             << "assembling mag face areas"
  82             << endl;
  83     }
  84
  85     // It is an error to attempt to recalculate
  86     // if the pointer is already set
  87     if (magSfPtr_)
  88     {
  89         FatalErrorIn("void fvMesh::makeMagSf()")
  90             << "mag face areas already exist"
  91             << abort(FatalError);
  92     }
  93
  94     // Note: Added stabilisation for faces with exactly zero area.
  95     // These should be caught on mesh checking but at least this stops
  96     // the code from producing Nans.
  97     magSfPtr_ = new surfaceScalarField
  98     (
  99         IOobject
 100         (
 101             "magSf",
 102             pointsInstance(),
 103             meshSubDir,
 104             *this,
 105             IOobject::NO_READ,
 106             IOobject::NO_WRITE,
 107             false
 108         ),
 109         mag(Sf()) + dimensionedScalar("vs", dimArea, VSMALL)
 110     );
 111 }
 112
 113
 114 void fvMesh::makeC() const
 115 {
 116     if (debug)
 117     {
 118         Info<< "void fvMesh::makeC() : "
 119             << "assembling cell centres"
 120             << endl;
 121     }
 122
 123     // It is an error to attempt to recalculate
 124     // if the pointer is already set
 125     if (CPtr_)
 126     {
 127         FatalErrorIn("fvMesh::makeC()")
 128             << "cell centres already exist"
 129             << abort(FatalError);
 130     }
 131
 132     CPtr_ = new slicedVolVectorField
 133     (
 134         IOobject
 135         (
 136             "C",
 137             pointsInstance(),
 138             meshSubDir,
 139             *this,
 140             IOobject::NO_READ,
 141             IOobject::NO_WRITE,
 142             false
 143         ),
 144         *this,
 145         dimLength,
 146         cellCentres(),
 147         faceCentres()
 148     );
 149 }
 150
 151
 152 void fvMesh::makeCf() const
 153 {
 154     if (debug)
 155     {
 156         Info<< "void fvMesh::makeCf() : "
 157             << "assembling face centres"
 158             << endl;
 159     }
 160
 161     // It is an error to attempt to recalculate
 162     // if the pointer is already set
 163     if (CfPtr_)
 164     {
 165         FatalErrorIn("fvMesh::makeCf()")
 166             << "face centres already exist"
 167             << abort(FatalError);
 168     }
 169
 170     CfPtr_ = new slicedSurfaceVectorField
 171     (
 172         IOobject
 173         (
 174             "Cf",
 175             pointsInstance(),
 176             meshSubDir,
 177             *this,
 178             IOobject::NO_READ,
 179             IOobject::NO_WRITE,
 180             false
 181         ),
 182         *this,
 183         dimLength,
 184         faceCentres()
 185     );
 186 }
 187
 188
 189 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 190
 191 const volScalarField::DimensionedInternalField& fvMesh::V() const
 192 {
 193     if (!VPtr_)
 194     {
 195         VPtr_ = new slicedVolScalarField::DimensionedInternalField
 196         (
 197             IOobject
 198             (
 199                 "V",
 200                 time().timeName(),
 201                 *this,
 202                 IOobject::NO_READ,
 203                 IOobject::NO_WRITE
 204             ),
 205             *this,
 206             dimVolume,
 207             cellVolumes()
 208         );
 209     }
 210
 211     return *static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
 212 }
 213
 214
 215 const volScalarField::DimensionedInternalField& fvMesh::V0() const
 216 {
 217     if (!V0Ptr_)
 218     {
 219         FatalErrorIn("fvMesh::V0() const")
 220             << "V0 is not available"
 221             << abort(FatalError);
 222     }
 223
 224     return *V0Ptr_;
 225 }
 226
 227
 228 volScalarField::DimensionedInternalField& fvMesh::setV0()
 229 {
 230     if (!V0Ptr_)
 231     {
 232         FatalErrorIn("fvMesh::setV0()")
 233             << "V0 is not available"
 234             << abort(FatalError);
 235     }
 236
 237     return *V0Ptr_;
 238 }
 239
 240
 241 const volScalarField::DimensionedInternalField& fvMesh::V00() const
 242 {
 243     if (!V00Ptr_)
 244     {
 245         V00Ptr_ = new DimensionedField<scalar, volMesh>
 246         (
 247             IOobject
 248             (
 249                 "V00",
 250                 time().timeName(),
 251                 *this,
 252                 IOobject::NO_READ,
 253                 IOobject::NO_WRITE
 254             ),
 255             V0()
 256         );
 257
 258         // If V00 is used then V0 should be stored for restart
 259         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
 260     }
 261
 262     return *V00Ptr_;
 263 }
 264
 265
 266 const surfaceVectorField& fvMesh::Sf() const
 267 {
 268     if (!SfPtr_)
 269     {
 270         makeSf();
 271     }
 272
 273     return *SfPtr_;
 274 }
 275
 276
 277 const surfaceScalarField& fvMesh::magSf() const
 278 {
 279     if (!magSfPtr_)
 280     {
 281         makeMagSf();
 282     }
 283
 284     return *magSfPtr_;
 285 }
 286
 287
 288 const volVectorField& fvMesh::C() const
 289 {
 290     if (!CPtr_)
 291     {
 292         makeC();
 293     }
 294
 295     return *CPtr_;
 296 }
 297
 298
 299 const surfaceVectorField& fvMesh::Cf() const
 300 {
 301     if (!CfPtr_)
 302     {
 303         makeCf();
 304     }
 305
 306     return *CfPtr_;
 307 }
 308
 309
 310 const surfaceScalarField& fvMesh::phi() const
 311 {
 312     if (!phiPtr_)
 313     {
 314         FatalErrorIn("fvMesh::phi()")
 315             << "mesh flux field does not exists, is the mesh actually moving?"
 316             << exit(FatalError);
 317     }
 318
 319     // Set zero current time
 320     // mesh motion fluxes if the time has been incremented
 321     if (phiPtr_->timeIndex() != time().timeIndex())
 322     {
 323         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 324     }
 325
 326     return *phiPtr_;
 327 }
 328
 329
 330 surfaceScalarField& fvMesh::setPhi()
 331 {
 332     if (!phiPtr_)
 333     {
 334         FatalErrorIn("fvMesh::setPhi()")
 335             << "mesh flux field does not exists, is the mesh actually moving?"
 336             << exit(FatalError);
 337     }
 338
 339     return *phiPtr_;
 340 }
 341
 342
 343 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 344
 345 } // End namespace Foam
 346
 347 // ************************************************************************* //