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