src/finiteVolume/fvMesh/fvMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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() const : "
  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() const")
  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() const : "
  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() const")
  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() const : "
 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() const")
 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 /* HJ, I think this is wrong.  HJ, 6/Jul/2010
 151     // Need to correct for cyclics transformation since absolute quantity.
 152     // Ok on processor patches since hold opposite cell centre (no
 153     // transformation)
 154     slicedVolVectorField& C = *CPtr_;
 155
 156     forAll(C.boundaryField(), patchi)
 157     {
 158         if (isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi]))
 159         {
 160             // Note: cyclic is not slice but proper field
 161             C.boundaryField()[patchi] == static_cast<const vectorField&>
 162             (
 163                 static_cast<const List<vector>&>
 164                 (
 165                     boundary_[patchi].patchSlice(faceCentres())
 166                 )
 167             );
 168         }
 169     }
 170 */
 171 }
 172
 173
 174 void fvMesh::makeCf() const
 175 {
 176     if (debug)
 177     {
 178         Info<< "void fvMesh::makeCf() const : "
 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() const")
 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 void fvMesh::makePhi() const
 212 {
 213     if (debug)
 214     {
 215         Info<< "void fvMesh::makePhi() const : "
 216             << "reading old time flux field if present and creating "
 217             << "zero current time flux field"
 218             << endl;
 219     }
 220
 221     // It is an error to attempt to recalculate
 222     // if the pointer is already set
 223     if (phiPtr_)
 224     {
 225         FatalErrorIn("fvMesh::makePhi() const")
 226             << "flux field already exists"
 227             << abort(FatalError);
 228     }
 229
 230     // Reading old time mesh motion flux if it exists and
 231     // creating zero current time mesh motion flux
 232
 233     scalar t0 = this->time().value() - this->time().deltaT().value();
 234
 235     IOobject meshPhiHeader
 236     (
 237         "meshPhi",
 238         this->time().timeName(t0),
 239         *this,
 240         IOobject::NO_READ
 241     );
 242
 243     if (meshPhiHeader.headerOk())
 244     {
 245         if (debug)
 246         {
 247             InfoIn("void fvMesh::makePhi()")
 248                 << "Reading mesh fluxes" << endl;
 249         }
 250
 251         phiPtr_ = new surfaceScalarField
 252         (
 253             IOobject
 254             (
 255                 "meshPhi",
 256                 this->time().timeName(t0),
 257                 *this,
 258                 IOobject::MUST_READ,
 259                 IOobject::AUTO_WRITE
 260             ),
 261             *this
 262         );
 263
 264         phiPtr_->oldTime();
 265
 266         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 267
 268         // This mesh is moving: set the motion to true
 269     }
 270     else
 271     {
 272         if (debug)
 273         {
 274             InfoIn("void fvMesh::makePhi()")
 275                 << "Creating null mesh motion fluxes" << endl;
 276         }
 277
 278         phiPtr_ = new surfaceScalarField
 279         (
 280             IOobject
 281             (
 282                 "meshPhi",
 283                 this->time().timeName(),
 284                 *this,
 285                 IOobject::NO_READ,
 286                 IOobject::AUTO_WRITE
 287             ),
 288             *this,
 289             dimensionedScalar("0", dimVolume/dimTime, 0.0)
 290         );
 291     }
 292 }
 293
 294
 295 void fvMesh::updatePhi(const scalarField& sweptVols) const
 296 {
 297     // Fill in mesh motion fluxes given swept volumes for all faces
 298
 299     if (!phiPtr_)
 300     {
 301         makePhi();
 302     }
 303
 304     surfaceScalarField& phi = *phiPtr_;
 305
 306     scalar rDeltaT = 1.0/time().deltaT().value();
 307
 308     phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
 309     phi.internalField() *= rDeltaT;
 310
 311     const fvPatchList& patches = boundary();
 312
 313     forAll (patches, patchI)
 314     {
 315         phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
 316         phi.boundaryField()[patchI] *= rDeltaT;
 317     }
 318 }
 319
 320
 321 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 322
 323 const volScalarField::DimensionedInternalField& fvMesh::V() const
 324 {
 325     if (!VPtr_)
 326     {
 327         if (debug)
 328         {
 329             InfoIn
 330             (
 331                 "const volScalarField::DimensionedInternalField& "
 332                 "fvMesh::V() const"
 333             )   << "Calculating cell volumes." << endl;
 334         }
 335
 336         VPtr_ = new DimensionedField<scalar, volMesh>
 337         (
 338             IOobject
 339             (
 340                 "V",
 341                 time().timeName(),
 342                 *this,
 343                 IOobject::NO_READ,
 344                 IOobject::NO_WRITE
 345             ),
 346             *this,
 347             dimVolume,
 348             cellVolumes()
 349         );
 350     }
 351
 352     return *VPtr_;
 353 }
 354
 355
 356 const volScalarField::DimensionedInternalField& fvMesh::V0() const
 357 {
 358     if (!V0Ptr_)
 359     {
 360         FatalErrorIn("fvMesh::V0() const")
 361             << "V0 is not available"
 362             << abort(FatalError);
 363     }
 364
 365     return *V0Ptr_;
 366 }
 367
 368
 369 DimensionedField<scalar, volMesh>& fvMesh::setV0()
 370 {
 371     if (!V0Ptr_)
 372     {
 373         FatalErrorIn("fvMesh::setV0()")
 374             << "V0 is not available"
 375             << abort(FatalError);
 376     }
 377
 378     // Delete old volume and mesh motion fluxes.  setV0() must be followed by
 379     // another mesh motion.  HJ, 25/Feb/2009
 380     deleteDemandDrivenData(phiPtr_);
 381     deleteDemandDrivenData(V0Ptr_);
 382
 383     if (debug)
 384     {
 385         InfoIn("DimensionedField<scalar, volMesh>& fvMesh::setV0()")
 386             << "Setting old cell volumes" << endl;
 387     }
 388
 389     V0Ptr_ = new DimensionedField<scalar, volMesh>
 390     (
 391         IOobject
 392         (
 393             "V0",
 394             time().timeName(),
 395             *this,
 396             IOobject::NO_READ,
 397             IOobject::NO_WRITE
 398         ),
 399         V()
 400     );
 401
 402     return *V0Ptr_;
 403 }
 404
 405
 406 const volScalarField::DimensionedInternalField& fvMesh::V00() const
 407 {
 408     if (!V00Ptr_)
 409     {
 410         V00Ptr_ = new DimensionedField<scalar, volMesh>
 411         (
 412             IOobject
 413             (
 414                 "V00",
 415                 time().timeName(),
 416                 *this,
 417                 IOobject::NO_READ,
 418                 IOobject::NO_WRITE
 419             ),
 420             V0()
 421         );
 422
 423         // If V00 is used then V0 should be stored for restart
 424         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
 425     }
 426
 427     return *V00Ptr_;
 428 }
 429
 430
 431 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
 432 {
 433     if (moving() && time().subCycling())
 434     {
 435         const TimeState& ts = time();
 436         const TimeState& ts0 = time().prevTimeState();
 437
 438         scalar tFrac =
 439         (
 440             ts.value() - (ts0.value() - ts0.deltaTValue())
 441         )/ts0.deltaTValue();
 442
 443         if (tFrac < (1 - SMALL))
 444         {
 445             return V0() + tFrac*(V() - V0());
 446         }
 447         else
 448         {
 449             return V();
 450         }
 451     }
 452     else
 453     {
 454         return V();
 455     }
 456 }
 457
 458
 459 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
 460 {
 461     if (moving() && time().subCycling())
 462     {
 463         const TimeState& ts = time();
 464         const TimeState& ts0 = time().prevTimeState();
 465
 466         scalar t0Frac =
 467         (
 468             (ts.value() - ts.deltaTValue())
 469           - (ts0.value() - ts0.deltaTValue())
 470         )/ts0.deltaTValue();
 471
 472         if (t0Frac > SMALL)
 473         {
 474             return V0() + t0Frac*(V() - V0());
 475         }
 476         else
 477         {
 478             return V0();
 479         }
 480     }
 481     else
 482     {
 483         return V0();
 484     }
 485 }
 486
 487
 488 const surfaceVectorField& fvMesh::Sf() const
 489 {
 490     if (!SfPtr_)
 491     {
 492         makeSf();
 493     }
 494
 495     return *SfPtr_;
 496 }
 497
 498
 499 const surfaceScalarField& fvMesh::magSf() const
 500 {
 501     if (!magSfPtr_)
 502     {
 503         makeMagSf();
 504     }
 505
 506     return *magSfPtr_;
 507 }
 508
 509
 510 const volVectorField& fvMesh::C() const
 511 {
 512     if (!CPtr_)
 513     {
 514         makeC();
 515     }
 516
 517     return *CPtr_;
 518 }
 519
 520
 521 const surfaceVectorField& fvMesh::Cf() const
 522 {
 523     if (!CfPtr_)
 524     {
 525         makeCf();
 526     }
 527
 528     return *CfPtr_;
 529 }
 530
 531
 532 const surfaceScalarField& fvMesh::phi() const
 533 {
 534     if (!phiPtr_)
 535     {
 536         makePhi();
 537     }
 538
 539     // Set zero current time
 540     // mesh motion fluxes if the time has been incremented
 541     if (phiPtr_->timeIndex() != time().timeIndex())
 542     {
 543         phiPtr_->oldTime();
 544
 545         if (debug)
 546         {
 547             InfoIn("const surfaceScalarField& fvMesh::phi() const")
 548                 << "Resetting mesh motion fluxes to zero" << endl;
 549         }
 550
 551         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 552     }
 553
 554     return *phiPtr_;
 555 }
 556
 557
 558 surfaceScalarField& fvMesh::setPhi()
 559 {
 560     if (!phiPtr_)
 561     {
 562         makePhi();
 563     }
 564
 565     return *phiPtr_;
 566 }
 567
 568
 569 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 570
 571 } // End namespace Foam
 572
 573 // ************************************************************************* //