src/finiteVolume/fvMesh/fvMesh.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 "volFields.H"
  28 #include "surfaceFields.H"
  29 #include "slicedVolFields.H"
  30 #include "slicedSurfaceFields.H"
  31 #include "SubField.H"
  32 #include "demandDrivenData.H"
  33 #include "fvMeshLduAddressing.H"
  34 #include "emptyPolyPatch.H"
  35 #include "mapPolyMesh.H"
  36 #include "MapFvFields.H"
  37 #include "fvMeshMapper.H"
  38 #include "mapClouds.H"
  39
  40 #include "volPointInterpolation.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(Foam::fvMesh, 0);
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 void Foam::fvMesh::clearGeomNotOldVol()
  50 {
  51     deleteDemandDrivenData(VPtr_);
  52
  53     deleteDemandDrivenData(SfPtr_);
  54     deleteDemandDrivenData(magSfPtr_);
  55     deleteDemandDrivenData(CPtr_);
  56     deleteDemandDrivenData(CfPtr_);
  57 }
  58
  59
  60 void Foam::fvMesh::clearGeom()
  61 {
  62     clearGeomNotOldVol();
  63
  64     deleteDemandDrivenData(V0Ptr_);
  65     deleteDemandDrivenData(V00Ptr_);
  66
  67     // Mesh motion flux cannot be deleted here because the old-time flux
  68     // needs to be saved.
  69
  70     // Geometry dependent object updated through call-back
  71     // and handled by polyMesh
  72     // HJ, 29/Aug/2010
  73 }
  74
  75
  76 void Foam::fvMesh::clearAddressing()
  77 {
  78     deleteDemandDrivenData(lduPtr_);
  79
  80     // Geometry dependent object updated through call-back
  81     // and handled by polyMesh
  82     // HJ, 29/Aug/2010
  83 }
  84
  85
  86 void Foam::fvMesh::clearOut()
  87 {
  88     clearGeom();
  89     surfaceInterpolation::clearOut();
  90
  91     clearAddressing();
  92
  93     // Clear mesh motion flux
  94     deleteDemandDrivenData(phiPtr_);
  95
  96     polyMesh::clearOut();
  97 }
  98
  99
 100 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 101
 102 Foam::fvMesh::fvMesh(const IOobject& io)
 103 :
 104     polyMesh(io),
 105     surfaceInterpolation(*this),
 106     boundary_(*this, boundaryMesh()),
 107     lduPtr_(NULL),
 108     curTimeIndex_(time().timeIndex()),
 109     VPtr_(NULL),
 110     V0Ptr_(NULL),
 111     V00Ptr_(NULL),
 112     SfPtr_(NULL),
 113     magSfPtr_(NULL),
 114     CPtr_(NULL),
 115     CfPtr_(NULL),
 116     phiPtr_(NULL)
 117 {
 118     if (debug)
 119     {
 120         Info<< "Constructing fvMesh from IOobject"
 121             << endl;
 122     }
 123
 124     // Check the existance of the cell volumes and read if present
 125     // and set the storage of V00
 126     if (isFile(time().timePath()/"V0"))
 127     {
 128         if (debug)
 129         {
 130             Info<< "Reading old cell volumes" << endl;
 131         }
 132
 133         V0Ptr_ = new DimensionedField<scalar, volMesh>
 134         (
 135             IOobject
 136             (
 137                 "V0",
 138                 time().timeName(),
 139                 *this,
 140                 IOobject::MUST_READ,
 141                 IOobject::NO_WRITE
 142             ),
 143             *this
 144         );
 145
 146         V00();
 147     }
 148
 149     // Check the existance of the mesh fluxes, read if present and set the
 150     // mesh to be moving
 151     if (isFile(time().timePath()/"meshPhi"))
 152     {
 153         if (debug)
 154         {
 155             Info<< "Reading motion fluxes" << endl;
 156         }
 157
 158         phiPtr_ = new surfaceScalarField
 159         (
 160             IOobject
 161             (
 162                 "meshPhi",
 163                 time().timeName(),
 164                 *this,
 165                 IOobject::MUST_READ,
 166                 IOobject::AUTO_WRITE
 167             ),
 168             *this
 169         );
 170
 171         // The mesh is now considered moving so the old-time cell volumes
 172         // will be required for the time derivatives so if they haven't been
 173         // read initialise to the current cell volumes
 174         if (!V0Ptr_)
 175         {
 176             V0Ptr_ = new DimensionedField<scalar, volMesh>
 177             (
 178                 IOobject
 179                 (
 180                     "V0",
 181                     time().timeName(),
 182                     *this,
 183                     IOobject::NO_READ,
 184                     IOobject::NO_WRITE,
 185                     false
 186                 ),
 187                 V()
 188             );
 189         }
 190
 191         moving(true);
 192     }
 193 }
 194
 195
 196 Foam::fvMesh::fvMesh
 197 (
 198     const IOobject& io,
 199     const Xfer<pointField>& points,
 200     const Xfer<faceList>& faces,
 201     const Xfer<labelList>& allOwner,
 202     const Xfer<labelList>& allNeighbour,
 203     const bool syncPar
 204 )
 205 :
 206     polyMesh(io, points, faces, allOwner, allNeighbour, syncPar),
 207     surfaceInterpolation(*this),
 208     boundary_(*this),
 209     lduPtr_(NULL),
 210     curTimeIndex_(time().timeIndex()),
 211     VPtr_(NULL),
 212     V0Ptr_(NULL),
 213     V00Ptr_(NULL),
 214     SfPtr_(NULL),
 215     magSfPtr_(NULL),
 216     CPtr_(NULL),
 217     CfPtr_(NULL),
 218     phiPtr_(NULL)
 219 {
 220     if (debug)
 221     {
 222         Info<< "Constructing fvMesh from components" << endl;
 223     }
 224 }
 225
 226
 227 Foam::fvMesh::fvMesh
 228 (
 229     const IOobject& io,
 230     const Xfer<pointField>& points,
 231     const Xfer<faceList>& faces,
 232     const Xfer<cellList>& cells,
 233     const bool syncPar
 234 )
 235 :
 236     polyMesh(io, points, faces, cells, syncPar),
 237     surfaceInterpolation(*this),
 238     boundary_(*this),
 239     lduPtr_(NULL),
 240     curTimeIndex_(time().timeIndex()),
 241     VPtr_(NULL),
 242     V0Ptr_(NULL),
 243     V00Ptr_(NULL),
 244     SfPtr_(NULL),
 245     magSfPtr_(NULL),
 246     CPtr_(NULL),
 247     CfPtr_(NULL),
 248     phiPtr_(NULL)
 249 {
 250     if (debug)
 251     {
 252         Info<< "Constructing fvMesh from components" << endl;
 253     }
 254 }
 255
 256
 257 Foam::fvMesh::fvMesh
 258 (
 259     const IOobject& io,
 260     const Xfer<pointField>& points,
 261     const cellShapeList& shapes,
 262     const faceListList& boundaryFaces,
 263     const wordList& boundaryPatchNames,
 264     const wordList& boundaryPatchTypes,
 265     const word& defaultBoundaryPatchName,
 266     const word& defaultBoundaryPatchType,
 267     const wordList& boundaryPatchPhysicalTypes,
 268     const bool syncPar
 269 )
 270 :
 271     polyMesh
 272     (
 273         io,
 274         points,
 275         shapes,
 276         boundaryFaces,
 277         boundaryPatchNames,
 278         boundaryPatchTypes,
 279         defaultBoundaryPatchName,
 280         defaultBoundaryPatchType,
 281         boundaryPatchPhysicalTypes,
 282         syncPar
 283     ),
 284     surfaceInterpolation(*this),
 285     boundary_(*this),
 286     lduPtr_(NULL),
 287     curTimeIndex_(time().timeIndex()),
 288     VPtr_(NULL),
 289     V0Ptr_(NULL),
 290     V00Ptr_(NULL),
 291     SfPtr_(NULL),
 292     magSfPtr_(NULL),
 293     CPtr_(NULL),
 294     CfPtr_(NULL),
 295     phiPtr_(NULL)
 296 {
 297     if (debug)
 298     {
 299         Info<< "Constructing fvMesh from components" << endl;
 300     }
 301 }
 302
 303
 304 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 305
 306 Foam::fvMesh::~fvMesh()
 307 {
 308     clearOut();
 309 }
 310
 311
 312 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 313
 314 void Foam::fvMesh::addFvPatches
 315 (
 316     const List<polyPatch*> & p,
 317     const bool validBoundary
 318 )
 319 {
 320     if (boundary().size())
 321     {
 322         FatalErrorIn
 323         (
 324             "fvMesh::addFvPatches(const List<polyPatch*>&, const bool)"
 325         )   << " boundary already exists"
 326             << abort(FatalError);
 327     }
 328
 329     // first add polyPatches
 330     addPatches(p, validBoundary);
 331     boundary_.addPatches(boundaryMesh());
 332 }
 333
 334
 335 void Foam::fvMesh::removeFvBoundary()
 336 {
 337     if (debug)
 338     {
 339         Info<< "void fvMesh::removeFvBoundary(): "
 340             << "Removing boundary patches."
 341             << endl;
 342     }
 343
 344     // Remove fvBoundaryMesh data first.
 345     boundary_.clear();
 346     boundary_.setSize(0);
 347     polyMesh::removeBoundary();
 348
 349     clearOut();
 350 }
 351
 352
 353 Foam::polyMesh::readUpdateState Foam::fvMesh::readUpdate()
 354 {
 355     if (debug)
 356     {
 357         Info<< "polyMesh::readUpdateState fvMesh::readUpdate() : "
 358             << "Updating fvMesh.  ";
 359     }
 360
 361     // Note: issues with update: should meshObject update happen
 362     // in polyMesh or fvMesh?  HJ, 18/Feb/2011
 363     polyMesh::readUpdateState state = polyMesh::readUpdate();
 364
 365     if (state == polyMesh::TOPO_PATCH_CHANGE)
 366     {
 367         if (debug)
 368         {
 369             Info << "Boundary and topological update" << endl;
 370         }
 371
 372         boundary_.readUpdate(boundaryMesh());
 373
 374         clearOut();
 375
 376     }
 377     else if (state == polyMesh::TOPO_CHANGE)
 378     {
 379         if (debug)
 380         {
 381             Info << "Topological update" << endl;
 382         }
 383
 384         clearOut();
 385     }
 386     else if (state == polyMesh::POINTS_MOVED)
 387     {
 388         if (debug)
 389         {
 390             Info << "Point motion update" << endl;
 391         }
 392
 393         clearGeom();
 394     }
 395     else
 396     {
 397         if (debug)
 398         {
 399             Info << "No update" << endl;
 400         }
 401     }
 402
 403     return state;
 404 }
 405
 406
 407 const Foam::fvBoundaryMesh& Foam::fvMesh::boundary() const
 408 {
 409     return boundary_;
 410 }
 411
 412
 413 const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
 414 {
 415     if (!lduPtr_)
 416     {
 417         lduPtr_ = new fvMeshLduAddressing(*this);
 418     }
 419
 420     return *lduPtr_;
 421 }
 422
 423
 424 void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap) const
 425 {
 426     if (debug)
 427     {
 428         Info<< "void fvMesh::mapFields(const mapPolyMesh& meshMap) const: "
 429             << "Mapping fv fields."
 430             << endl;
 431     }
 432
 433     // Create a mapper
 434     const fvMeshMapper mapper(*this, meshMap);
 435
 436     // Map all the volFields in the objectRegistry
 437     MapGeometricFields<scalar, fvPatchField, fvMeshMapper, volMesh>(mapper);
 438     MapGeometricFields<vector, fvPatchField, fvMeshMapper, volMesh>(mapper);
 439     MapGeometricFields<sphericalTensor, fvPatchField, fvMeshMapper, volMesh>
 440         (mapper);
 441
 442     MapGeometricFields<symmTensor, fvPatchField, fvMeshMapper, volMesh>
 443         (mapper);
 444
 445     MapGeometricFields<symmTensor4thOrder, fvPatchField, fvMeshMapper, volMesh>
 446       (mapper);
 447
 448     MapGeometricFields<diagTensor, fvPatchField, fvMeshMapper, volMesh>
 449       (mapper);
 450
 451     MapGeometricFields<tensor, fvPatchField, fvMeshMapper, volMesh>(mapper);
 452
 453     // Map all the surfaceFields in the objectRegistry
 454     MapGeometricFields<scalar, fvsPatchField, fvMeshMapper, surfaceMesh>
 455         (mapper);
 456
 457     MapGeometricFields<vector, fvsPatchField, fvMeshMapper, surfaceMesh>
 458         (mapper);
 459
 460     MapGeometricFields
 461         <sphericalTensor, fvsPatchField, fvMeshMapper, surfaceMesh>(mapper);
 462     MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 463         (mapper);
 464
 465     MapGeometricFields<symmTensor4thOrder, fvsPatchField, fvMeshMapper, surfaceMesh>
 466       (mapper);
 467
 468     MapGeometricFields<diagTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 469       (mapper);
 470
 471     MapGeometricFields<tensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 472         (mapper);
 473
 474     // Map all the clouds in the objectRegistry
 475     mapClouds(*this, meshMap);
 476 }
 477
 478
 479 void Foam::fvMesh::mapOldVolumes(const mapPolyMesh& meshMap)
 480 {
 481     const labelList& cellMap = meshMap.cellMap();
 482
 483     // Map the old volume. Just map to new cell labels.
 484     if (V0Ptr_)
 485     {
 486         if (debug)
 487         {
 488             InfoIn("void fvMesh::mapOldVolumes(const mapPolyMesh& meshMap)")
 489                 << "Mapping old cell volumes." << endl;
 490         }
 491
 492         scalarField& V0 = *V0Ptr_;
 493
 494         scalarField savedV0(V0);
 495         V0.setSize(nCells());
 496
 497         forAll (V0, i)
 498         {
 499             if (cellMap[i] > -1)
 500             {
 501                 V0[i] = savedV0[cellMap[i]];
 502             }
 503             else
 504             {
 505                 V0[i] = 0.0;
 506             }
 507         }
 508     }
 509
 510     // Map the old-old volume. Just map to new cell labels.
 511     if (V00Ptr_)
 512     {
 513         if (debug)
 514         {
 515             InfoIn("void fvMesh::mapOldVolumes(const mapPolyMesh& meshMap)")
 516                 << "Mapping old-old cell volumes." << endl;
 517         }
 518
 519         scalarField& V00 = *V00Ptr_;
 520
 521         scalarField savedV00(V00);
 522         V00.setSize(nCells());
 523
 524         forAll (V00, i)
 525         {
 526             if (cellMap[i] > -1)
 527             {
 528                 V00[i] = savedV00[cellMap[i]];
 529             }
 530             else
 531             {
 532                 V00[i] = 0.0;
 533             }
 534         }
 535     }
 536 }
 537
 538
 539 void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
 540 {
 541     // Update polyMesh. This needs to keep volume existent!
 542     polyMesh::updateMesh(mpm);
 543
 544     surfaceInterpolation::clearOut();
 545     clearGeomNotOldVol();
 546
 547     // Map all fields
 548     mapFields(mpm);
 549
 550     // Map old-volumes
 551     mapOldVolumes(mpm);
 552
 553     clearAddressing();
 554
 555     // Mesh morphing should also clear out the surfaceInterpolation.
 556     // This is a temporary solution
 557     surfaceInterpolation::movePoints();
 558
 559     // Function object update moved to polyMesh
 560     // HJ, 29/Aug/2010
 561 }
 562
 563
 564 void Foam::fvMesh::syncUpdateMesh()
 565 {
 566     // Update polyMesh. This needs to keep volume existent!
 567     polyMesh::syncUpdateMesh();
 568
 569     // Not sure how much clean-up is needed here.  HJ, 27/Nov/2009
 570
 571     surfaceInterpolation::clearOut();
 572     clearGeomNotOldVol();
 573
 574     clearAddressing();
 575
 576     // handleMorph() should also clear out the surfaceInterpolation.
 577     // This is a temporary solution
 578     surfaceInterpolation::movePoints();
 579
 580     // Function object update moved to polyMesh
 581     // HJ, 29/Aug/2010
 582 }
 583
 584
 585 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 586 {
 587     // Grab old time volumes if the time has been incremented
 588     if (curTimeIndex_ < time().timeIndex())
 589     {
 590         if (V00Ptr_ && V0Ptr_)
 591         {
 592             if (debug)
 593             {
 594                 InfoIn("void fvMesh::movePoints(const pointField& p)")
 595                     << "Grabbing old-old cell volumes." << endl;
 596             }
 597
 598             *V00Ptr_ = *V0Ptr_;
 599         }
 600
 601         if (V0Ptr_)
 602         {
 603             if (debug)
 604             {
 605                 InfoIn("void fvMesh::movePoints(const pointField& p)")
 606                     << "Grabbing old cell volumes." << endl;
 607             }
 608
 609             *V0Ptr_ = V();
 610         }
 611         else
 612         {
 613             if (debug)
 614             {
 615                 InfoIn("void fvMesh::movePoints(const pointField& p)")
 616                     << "Creating old cell volumes." << endl;
 617             }
 618
 619             V0Ptr_ = new DimensionedField<scalar, volMesh>
 620             (
 621                 IOobject
 622                 (
 623                     "V0",
 624                     time().timeName(),
 625                     *this,
 626                     IOobject::NO_READ,
 627                     IOobject::NO_WRITE
 628                 ),
 629                 V()
 630             );
 631         }
 632
 633         curTimeIndex_ = time().timeIndex();
 634     }
 635
 636
 637     // Delete out of date geometrical information
 638     clearGeomNotOldVol();
 639
 640
 641     if (!phiPtr_)
 642     {
 643         // Create mesh motion flux
 644         if (debug)
 645         {
 646             InfoIn("tmp<scalarField> fvMesh::movePoints(const pointField& p)")
 647                 << "Creating new mesh motion fluxes" << endl;
 648         }
 649
 650         phiPtr_ = new surfaceScalarField
 651         (
 652             IOobject
 653             (
 654                 "meshPhi",
 655                 this->time().timeName(),
 656                 *this,
 657                 IOobject::NO_READ,
 658                 IOobject::AUTO_WRITE
 659             ),
 660             *this,
 661             dimVolume/dimTime
 662         );
 663     }
 664     else
 665     {
 666         // Grab old time mesh motion fluxes if the time has been incremented
 667         if (phiPtr_->timeIndex() < time().timeIndex())
 668         {
 669             phiPtr_->oldTime();
 670         }
 671     }
 672
 673     // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
 674     tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
 675
 676     updatePhi(tsweptVols());
 677
 678     boundary_.movePoints();
 679     surfaceInterpolation::movePoints();
 680
 681     // Function object update moved to polyMesh
 682     // HJ, 29/Aug/2010
 683
 684     return tsweptVols;
 685 }
 686
 687
 688 bool Foam::fvMesh::writeObjects
 689 (
 690     IOstream::streamFormat fmt,
 691     IOstream::versionNumber ver,
 692     IOstream::compressionType cmp
 693 ) const
 694 {
 695     return polyMesh::writeObject(fmt, ver, cmp);
 696 }
 697
 698
 699 //- Write mesh using IO settings from the time
 700 bool Foam::fvMesh::write() const
 701 {
 702     return polyMesh::write();
 703 }
 704
 705
 706 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 707
 708 bool Foam::fvMesh::operator!=(const fvMesh& bm) const
 709 {
 710     return &bm != this;
 711 }
 712
 713
 714 bool Foam::fvMesh::operator==(const fvMesh& bm) const
 715 {
 716     return &bm == this;
 717 }
 718
 719
 720 // ************************************************************************* //