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