src/finiteVolume/fvMesh/fvMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "fvMesh.H"
  27 #include "foamTime.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 #include "cyclicGgiFvPatchFields.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() const : "
  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() const")
  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() const : "
  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() const")
  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.  HJ, date deleted
  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() const : "
 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() const")
 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     // This piece of code is necessary for cyclic and cyclicGgi interfaces
 152     // using a separationOffset transform.
 153     // Those two interfaces will be using the method ::patchNeighbourField()
 154     // to evaluate the field C on the shadow patch. For cyclic and cyclicGgi
 155     // translational-only interfaces, the separationOffset transform is never
 156     // applied directly in ::patchNeighbourField() because this transform is
 157     // only pertinent for 3D coordinates, and the method ::patchNeighbourField()
 158     // does not discriminate the type of field it is operating on.
 159     // So, because the separationOffset transform is not applied, the evaluation
 160     // of a 3D position field like 'C' will always be wrong on the shadow patches
 161     // of translational cyclic and cyclicGgi interfaces.
 162     // For cyclic and cyclicGgi interfaces using a rotational transform, the
 163     // evaluation of the field C will be valid, but since we are only
 164     // interested in the patch face centers for these interfaces, we can override
 165     // those values as well.
 166     // See also:
 167     //    https://sourceforge.net/apps/mantisbt/openfoam-extend/view.php?id=42
 168     // MB, 12/Dec/2010
 169     //
 170     // Need to correct for cyclics transformation since absolute quantity.
 171     // Ok on processor patches since hold opposite cell centre (no
 172     // transformation)
 173     slicedVolVectorField& C = *CPtr_;
 174
 175     forAll(C.boundaryField(), patchi)
 176     {
 177         if
 178         (
 179             isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi])
 180          || isA<cyclicGgiFvPatchVectorField>(C.boundaryField()[patchi])
 181         )
 182         {
 183             // Note: cyclic is not slice but proper field
 184             C.boundaryField()[patchi] == static_cast<const vectorField&>
 185             (
 186                 static_cast<const List<vector>&>
 187                 (
 188                     boundary_[patchi].patchSlice(faceCentres())
 189                 )
 190             );
 191         }
 192     }
 193 }
 194
 195
 196 void fvMesh::makeCf() const
 197 {
 198     if (debug)
 199     {
 200         Info<< "void fvMesh::makeCf() const : "
 201             << "assembling face centres"
 202             << endl;
 203     }
 204
 205     // It is an error to attempt to recalculate
 206     // if the pointer is already set
 207     if (CfPtr_)
 208     {
 209         FatalErrorIn("fvMesh::makeCf() const")
 210             << "face centres already exist"
 211             << abort(FatalError);
 212     }
 213
 214     CfPtr_ = new slicedSurfaceVectorField
 215     (
 216         IOobject
 217         (
 218             "Cf",
 219             pointsInstance(),
 220             meshSubDir,
 221             *this,
 222             IOobject::NO_READ,
 223             IOobject::NO_WRITE,
 224             false
 225         ),
 226         *this,
 227         dimLength,
 228         faceCentres()
 229     );
 230 }
 231
 232
 233 void fvMesh::makePhi() const
 234 {
 235     if (debug)
 236     {
 237         Info<< "void fvMesh::makePhi() const : "
 238             << "reading old time flux field if present and creating "
 239             << "zero current time flux field"
 240             << endl;
 241     }
 242
 243     // It is an error to attempt to recalculate
 244     // if the pointer is already set
 245     if (phiPtr_)
 246     {
 247         FatalErrorIn("fvMesh::makePhi() const")
 248             << "flux field already exists"
 249             << abort(FatalError);
 250     }
 251
 252     // Reading old time mesh motion flux if it exists and
 253     // creating zero current time mesh motion flux
 254
 255     scalar t0 = this->time().value() - this->time().deltaT().value();
 256
 257     IOobject meshPhiHeader
 258     (
 259         "meshPhi",
 260         this->time().timeName(t0),
 261         *this,
 262         IOobject::NO_READ
 263     );
 264
 265     if (meshPhiHeader.headerOk())
 266     {
 267         if (debug)
 268         {
 269             InfoIn("void fvMesh::makePhi()")
 270                 << "Reading mesh fluxes" << endl;
 271         }
 272
 273         phiPtr_ = new surfaceScalarField
 274         (
 275             IOobject
 276             (
 277                 "meshPhi",
 278                 this->time().timeName(t0),
 279                 *this,
 280                 IOobject::MUST_READ,
 281                 IOobject::AUTO_WRITE
 282             ),
 283             *this
 284         );
 285
 286         phiPtr_->oldTime();
 287
 288         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 289
 290         // This mesh is moving: set the motion to true
 291     }
 292     else
 293     {
 294         if (debug)
 295         {
 296             InfoIn("void fvMesh::makePhi()")
 297                 << "Creating null mesh motion fluxes" << endl;
 298         }
 299
 300         phiPtr_ = new surfaceScalarField
 301         (
 302             IOobject
 303             (
 304                 "meshPhi",
 305                 this->time().timeName(),
 306                 *this,
 307                 IOobject::NO_READ,
 308                 IOobject::AUTO_WRITE
 309             ),
 310             *this,
 311             dimensionedScalar("0", dimVolume/dimTime, 0.0)
 312         );
 313     }
 314 }
 315
 316
 317 void fvMesh::updatePhi(const scalarField& sweptVols) const
 318 {
 319     // Fill in mesh motion fluxes given swept volumes for all faces
 320
 321     if (!phiPtr_)
 322     {
 323         makePhi();
 324     }
 325
 326     surfaceScalarField& phi = *phiPtr_;
 327
 328     scalar rDeltaT = 1.0/time().deltaT().value();
 329
 330     phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
 331     phi.internalField() *= rDeltaT;
 332
 333     const fvPatchList& patches = boundary();
 334
 335     forAll (patches, patchI)
 336     {
 337         phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
 338         phi.boundaryField()[patchI] *= rDeltaT;
 339     }
 340 }
 341
 342
 343 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 344
 345 const volScalarField::DimensionedInternalField& fvMesh::V() const
 346 {
 347     if (!VPtr_)
 348     {
 349         if (debug)
 350         {
 351             InfoIn
 352             (
 353                 "const volScalarField::DimensionedInternalField& "
 354                 "fvMesh::V() const"
 355             )   << "Calculating cell volumes." << endl;
 356         }
 357
 358         VPtr_ = new DimensionedField<scalar, volMesh>
 359         (
 360             IOobject
 361             (
 362                 "V",
 363                 time().timeName(),
 364                 *this,
 365                 IOobject::NO_READ,
 366                 IOobject::NO_WRITE
 367             ),
 368             *this,
 369             dimVolume,
 370             cellVolumes()
 371         );
 372     }
 373
 374     return *VPtr_;
 375 }
 376
 377
 378 const volScalarField::DimensionedInternalField& fvMesh::V0() const
 379 {
 380     if (!V0Ptr_)
 381     {
 382         FatalErrorIn("fvMesh::V0() const")
 383             << "V0 is not available"
 384             << abort(FatalError);
 385     }
 386
 387     return *V0Ptr_;
 388 }
 389
 390
 391 DimensionedField<scalar, volMesh>& fvMesh::setV0()
 392 {
 393     // Delete old volume and mesh motion fluxes.  setV0() must be followed by
 394     // another mesh motion.  HJ, 25/Feb/2009
 395     deleteDemandDrivenData(phiPtr_);
 396     deleteDemandDrivenData(V0Ptr_);
 397
 398     if (debug)
 399     {
 400         InfoIn("DimensionedField<scalar, volMesh>& fvMesh::setV0()")
 401             << "Setting old cell volumes" << endl;
 402     }
 403
 404     V0Ptr_ = new DimensionedField<scalar, volMesh>
 405     (
 406         IOobject
 407         (
 408             "V0",
 409             time().timeName(),
 410             *this,
 411             IOobject::NO_READ,
 412             IOobject::NO_WRITE
 413         ),
 414         V()
 415     );
 416
 417     return *V0Ptr_;
 418 }
 419
 420
 421 const volScalarField::DimensionedInternalField& fvMesh::V00() const
 422 {
 423     if (!V00Ptr_)
 424     {
 425         V00Ptr_ = new DimensionedField<scalar, volMesh>
 426         (
 427             IOobject
 428             (
 429                 "V00",
 430                 time().timeName(),
 431                 *this,
 432                 IOobject::NO_READ,
 433                 IOobject::NO_WRITE
 434             ),
 435             V0()
 436         );
 437
 438         // If V00 is used then V0 should be stored for restart
 439         V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
 440     }
 441
 442     return *V00Ptr_;
 443 }
 444
 445
 446 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc() const
 447 {
 448     if (moving() && time().subCycling())
 449     {
 450         const TimeState& ts = time();
 451         const TimeState& ts0 = time().prevTimeState();
 452
 453         scalar tFrac =
 454         (
 455             ts.value() - (ts0.value() - ts0.deltaTValue())
 456         )/ts0.deltaTValue();
 457
 458         if (tFrac < (1 - SMALL))
 459         {
 460             return V0() + tFrac*(V() - V0());
 461         }
 462         else
 463         {
 464             return V();
 465         }
 466     }
 467     else
 468     {
 469         return V();
 470     }
 471 }
 472
 473
 474 tmp<volScalarField::DimensionedInternalField> fvMesh::Vsc0() const
 475 {
 476     if (moving() && time().subCycling())
 477     {
 478         const TimeState& ts = time();
 479         const TimeState& ts0 = time().prevTimeState();
 480
 481         scalar t0Frac =
 482         (
 483             (ts.value() - ts.deltaTValue())
 484           - (ts0.value() - ts0.deltaTValue())
 485         )/ts0.deltaTValue();
 486
 487         if (t0Frac > SMALL)
 488         {
 489             return V0() + t0Frac*(V() - V0());
 490         }
 491         else
 492         {
 493             return V0();
 494         }
 495     }
 496     else
 497     {
 498         return V0();
 499     }
 500 }
 501
 502
 503 const surfaceVectorField& fvMesh::Sf() const
 504 {
 505     if (!SfPtr_)
 506     {
 507         makeSf();
 508     }
 509
 510     return *SfPtr_;
 511 }
 512
 513
 514 const surfaceScalarField& fvMesh::magSf() const
 515 {
 516     if (!magSfPtr_)
 517     {
 518         makeMagSf();
 519     }
 520
 521     return *magSfPtr_;
 522 }
 523
 524
 525 const volVectorField& fvMesh::C() const
 526 {
 527     if (!CPtr_)
 528     {
 529         makeC();
 530     }
 531
 532     return *CPtr_;
 533 }
 534
 535
 536 const surfaceVectorField& fvMesh::Cf() const
 537 {
 538     if (!CfPtr_)
 539     {
 540         makeCf();
 541     }
 542
 543     return *CfPtr_;
 544 }
 545
 546
 547 const surfaceScalarField& fvMesh::phi() const
 548 {
 549     if (!phiPtr_)
 550     {
 551         makePhi();
 552     }
 553
 554     // Set zero current time
 555     // mesh motion fluxes if the time has been incremented
 556     if (phiPtr_->timeIndex() != time().timeIndex())
 557     {
 558         phiPtr_->oldTime();
 559
 560         if (debug)
 561         {
 562             InfoIn("const surfaceScalarField& fvMesh::phi() const")
 563                 << "Resetting mesh motion fluxes to zero" << endl;
 564         }
 565
 566         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
 567     }
 568
 569     return *phiPtr_;
 570 }
 571
 572
 573 surfaceScalarField& fvMesh::setPhi()
 574 {
 575     if (!phiPtr_)
 576     {
 577         makePhi();
 578     }
 579
 580     return *phiPtr_;
 581 }
 582
 583
 584 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 585
 586 } // End namespace Foam
 587
 588 // ************************************************************************* //