src/finiteVolume/fvMesh/fvMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "fvMesh.H"
  27 #include "Time.H"
  28 #include "volFields.H"
  29 #include "surfaceFields.H"
  30 #include "slicedVolFields.H"
  31 #include "slicedSurfaceFields.H"
  32 #include "SubField.H"
  33 #include "cyclicFvPatchFields.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             IOobject::NO_READ,
  69             IOobject::NO_WRITE,
  70             false
  71         ),
  72         *this,
  73         dimArea,
  74         faceAreas()
  75     );
  76 }
  77
  78
  79 void fvMesh::makeMagSf() const
  80 {
  81     if (debug)
  82     {
  83         Info<< "void fvMesh::makeMagSf() : "
  84             << "assembling mag face areas"
  85             << endl;
  86     }
  87
  88     // It is an error to attempt to recalculate
  89     // if the pointer is already set
  90     if (magSfPtr_)
  91     {
  92         FatalErrorIn("void fvMesh::makeMagSf()")
  93             << "mag face areas already exist"
  94             << abort(FatalError);
  95     }
  96
  97     // Note: Added stabilisation for faces with exactly zero area.
  98     // These should be caught on mesh checking but at least this stops
  99     // the code from producing Nans.
 100     magSfPtr_ = new surfaceScalarField
 101     (
 102         IOobject
 103         (
 104             "magSf",
 105             pointsInstance(),
 106             meshSubDir,
 107             *this,
 108             IOobject::NO_READ,
 109             IOobject::NO_WRITE,
 110             false
 111         ),
 112         mag(Sf()) + dimensionedScalar("vs", dimArea, VSMALL)
 113     );
 114 }
 115
 116
 117 void fvMesh::makeC() const
 118 {
 119     if (debug)
 120     {
 121         Info<< "void fvMesh::makeC() : "
 122             << "assembling cell centres"
 123             << endl;
 124     }
 125
 126     // It is an error to attempt to recalculate
 127     // if the pointer is already set
 128     if (CPtr_)
 129     {
 130         FatalErrorIn("fvMesh::makeC()")
 131             << "cell centres already exist"
 132             << abort(FatalError);
 133     }
 134
 135     CPtr_ = new slicedVolVectorField
 136     (
 137         IOobject
 138         (
 139             "C",
 140             pointsInstance(),
 141             meshSubDir,
 142             *this,
 143             IOobject::NO_READ,
 144             IOobject::NO_WRITE,
 145             false
 146         ),
 147         *this,
 148         dimLength,
 149         cellCentres(),
 150         faceCentres()
 151     );
 152
 153
 154     // Need to correct for cyclics transformation since absolute quantity.
 155     // Ok on processor patches since hold opposite cell centre (no
 156     // transformation)
 157     slicedVolVectorField& C = *CPtr_;
 158
 159     forAll(C.boundaryField(), patchi)
 160     {
 161         if (isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi]))
 162         {
 163             // Note: cyclic is not slice but proper field
 164             C.boundaryField()[patchi] == static_cast<const vectorField&>
 165             (
 166                 static_cast<const List<vector>&>
 167                 (
 168                     boundary_[patchi].patchSlice(faceCentres())
 169                 )
 170             );
 171         }
 172     }
 173 }
 174
 175
 176 void fvMesh::makeCf() const
 177 {
 178     if (debug)
 179     {
 180         Info<< "void fvMesh::makeCf() : "
 181             << "assembling face centres"
 182             << endl;
 183     }
 184
 185     // It is an error to attempt to recalculate
 186     // if the pointer is already set
 187     if (CfPtr_)
 188     {
 189         FatalErrorIn("fvMesh::makeCf()")
 190             << "face centres already exist"
 191             << abort(FatalError);
 192     }
 193
 194     CfPtr_ = new slicedSurfaceVectorField
 195     (
 196         IOobject
 197         (
 198             "Cf",
 199             pointsInstance(),
 200             meshSubDir,
 201             *this,
 202             IOobject::NO_READ,
 203             IOobject::NO_WRITE,
 204             false
 205         ),
 206         *this,
 207         dimLength,
 208         faceCentres()
 209     );
 210 }
 211
 212
 213 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 214
 215 const volScalarField::DimensionedInternalField& fvMesh::V() const
 216 {
 217     if (!VPtr_)
 218     {
 219         VPtr_ = new slicedVolScalarField::DimensionedInternalField
 220         (
 221             IOobject
 222             (
 223                 "V",
 224                 time().timeName(),
 225                 *this,
 226                 IOobject::NO_READ,
 227                 IOobject::NO_WRITE
 228             ),
 229             *this,
 230             dimVolume,
 231             cellVolumes()
 232         );
 233     }
 234
 235     return *static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
 236 }
 237
 238
 239 const volScalarField::DimensionedInternalField& fvMesh::V0() const
 240 {
 241     if (!V0Ptr_)
 242     {
 243         FatalErrorIn("fvMesh::V0() const")
 244             << "V0 is not available"
 245             << abort(FatalError);
 246     }
 247
 248     return *V0Ptr_;
 249 }
 250
 251
 252 volScalarField::DimensionedInternalField& fvMesh::setV0()
 253 {
 254     if (!V0Ptr_)
 255     {
 256         FatalErrorIn("fvMesh::setV0()")
 257             << "V0 is not available"
 258             << abort(FatalError);
 259     }
 260
 261     return *V0Ptr_;
 262 }
 263
 264
 265 const volScalarField::DimensionedInternalField& fvMesh::V00() const
 266 {
 267     if (!V00Ptr_)
 268     {
 269         V00Ptr_ = new DimensionedField<scalar, volMesh>
 270         (
 271             IOobject
 272             (
 273                 "V00",
 274                 time().timeName(),
 275                 *this,
 276                 IOobject::NO_READ,
 277                 IOobject::NO_WRITE
 278             ),
 279             V0()
 280         );
 281
 282         // If V00 is used then V0 should be stored for restart
 283         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
 284     }
 285
 286     return *V00Ptr_;
 287 }
 288
 289
 290 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
 291 {
 292     if (moving() && time().subCycling())
 293     {
 294         const TimeState& ts = time();
 295         const TimeState& ts0 = time().prevTimeState();
 296
 297         scalar tFrac =
 298         (
 299             ts.value() - (ts0.value() - ts0.deltaTValue())
 300         )/ts0.deltaTValue();
 301
 302         if (tFrac < (1 - SMALL))
 303         {
 304             return V0() + tFrac*(V() - V0());
 305         }
 306         else
 307         {
 308             return V();
 309         }
 310     }
 311     else
 312     {
 313         return V();
 314     }
 315 }
 316
 317
 318 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
 319 {
 320     if (moving() && time().subCycling())
 321     {
 322         const TimeState& ts = time();
 323         const TimeState& ts0 = time().prevTimeState();
 324
 325         scalar t0Frac =
 326         (
 327             (ts.value() - ts.deltaTValue())
 328           - (ts0.value() - ts0.deltaTValue())
 329         )/ts0.deltaTValue();
 330
 331         if (t0Frac > SMALL)
 332         {
 333             return V0() + t0Frac*(V() - V0());
 334         }
 335         else
 336         {
 337             return V0();
 338         }
 339     }
 340     else
 341     {
 342         return V0();
 343     }
 344 }
 345
 346
 347 const surfaceVectorField& fvMesh::Sf() const
 348 {
 349     if (!SfPtr_)
 350     {
 351         makeSf();
 352     }
 353
 354     return *SfPtr_;
 355 }
 356
 357
 358 const surfaceScalarField& fvMesh::magSf() const
 359 {
 360     if (!magSfPtr_)
 361     {
 362         makeMagSf();
 363     }
 364
 365     return *magSfPtr_;
 366 }
 367
 368
 369 const volVectorField& fvMesh::C() const
 370 {
 371     if (!CPtr_)
 372     {
 373         makeC();
 374     }
 375
 376     return *CPtr_;
 377 }
 378
 379
 380 const surfaceVectorField& fvMesh::Cf() const
 381 {
 382     if (!CfPtr_)
 383     {
 384         makeCf();
 385     }
 386
 387     return *CfPtr_;
 388 }
 389
 390
 391 const surfaceScalarField& fvMesh::phi() const
 392 {
 393     if (!phiPtr_)
 394     {
 395         FatalErrorIn("fvMesh::phi()")
 396             << "mesh flux field does not exists, is the mesh actually moving?"
 397             << exit(FatalError);
 398     }
 399
 400     // Set zero current time
 401     // mesh motion fluxes if the time has been incremented
 402     if (phiPtr_->timeIndex() != time().timeIndex())
 403     {
 404         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 405     }
 406
 407     return *phiPtr_;
 408 }
 409
 410
 411 surfaceScalarField& fvMesh::setPhi()
 412 {
 413     if (!phiPtr_)
 414     {
 415         FatalErrorIn("fvMesh::setPhi()")
 416             << "mesh flux field does not exists, is the mesh actually moving?"
 417             << exit(FatalError);
 418     }
 419
 420     return *phiPtr_;
 421 }
 422
 423
 424 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 425
 426 } // End namespace Foam
 427
 428 // ************************************************************************* //