src/finiteVolume/fvMesh/fvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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 #include "extendedLeastSquaresVectors.H"
  42 #include "extendedLeastSquaresVectors.H"
  43 #include "leastSquaresVectors.H"
  44 #include "CentredFitData.H"
  45 #include "linearFitPolynomial.H"
  46 #include "quadraticFitPolynomial.H"
  47 #include "quadraticLinearFitPolynomial.H"
  48 //#include "quadraticFitSnGradData.H"
  49 #include "skewCorrectionVectors.H"
  50
  51
  52 #include "centredCECCellToFaceStencilObject.H"
  53 #include "centredCFCCellToFaceStencilObject.H"
  54 #include "centredCPCCellToFaceStencilObject.H"
  55 #include "centredFECCellToFaceStencilObject.H"
  56 #include "upwindCECCellToFaceStencilObject.H"
  57 #include "upwindCFCCellToFaceStencilObject.H"
  58 #include "upwindCPCCellToFaceStencilObject.H"
  59 #include "upwindFECCellToFaceStencilObject.H"
  60
  61 #include "centredCFCFaceToCellStencilObject.H"
  62
  63 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  64
  65 defineTypeNameAndDebug(Foam::fvMesh, 0);
  66
  67 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  68
  69 void Foam::fvMesh::clearGeomNotOldVol()
  70 {
  71     slicedVolScalarField::DimensionedInternalField* VPtr =
  72         static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
  73     deleteDemandDrivenData(VPtr);
  74     VPtr_ = NULL;
  75
  76     deleteDemandDrivenData(SfPtr_);
  77     deleteDemandDrivenData(magSfPtr_);
  78     deleteDemandDrivenData(CPtr_);
  79     deleteDemandDrivenData(CfPtr_);
  80 }
  81
  82
  83 void Foam::fvMesh::clearGeom()
  84 {
  85     clearGeomNotOldVol();
  86
  87     deleteDemandDrivenData(V0Ptr_);
  88     deleteDemandDrivenData(V00Ptr_);
  89
  90     // Mesh motion flux cannot be deleted here because the old-time flux
  91     // needs to be saved.
  92
  93     // Things geometry dependent that are not updated.
  94     volPointInterpolation::Delete(*this);
  95     extendedLeastSquaresVectors::Delete(*this);
  96     leastSquaresVectors::Delete(*this);
  97     CentredFitData<linearFitPolynomial>::Delete(*this);
  98     CentredFitData<quadraticFitPolynomial>::Delete(*this);
  99     CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
 100     skewCorrectionVectors::Delete(*this);
 101     //quadraticFitSnGradData::Delete(*this);
 102 }
 103
 104
 105 void Foam::fvMesh::clearAddressing()
 106 {
 107     deleteDemandDrivenData(lduPtr_);
 108
 109     // Hack until proper callbacks. Below are all the fvMesh-MeshObjects.
 110
 111     volPointInterpolation::Delete(*this);
 112     extendedLeastSquaresVectors::Delete(*this);
 113     leastSquaresVectors::Delete(*this);
 114     CentredFitData<linearFitPolynomial>::Delete(*this);
 115     CentredFitData<quadraticFitPolynomial>::Delete(*this);
 116     CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
 117     skewCorrectionVectors::Delete(*this);
 118     //quadraticFitSnGradData::Delete(*this);
 119
 120     centredCECCellToFaceStencilObject::Delete(*this);
 121     centredCFCCellToFaceStencilObject::Delete(*this);
 122     centredCPCCellToFaceStencilObject::Delete(*this);
 123     centredFECCellToFaceStencilObject::Delete(*this);
 124     // Is this geometry related - cells distorting to upwind direction?
 125     upwindCECCellToFaceStencilObject::Delete(*this);
 126     upwindCFCCellToFaceStencilObject::Delete(*this);
 127     upwindCPCCellToFaceStencilObject::Delete(*this);
 128     upwindFECCellToFaceStencilObject::Delete(*this);
 129
 130     centredCFCFaceToCellStencilObject::Delete(*this);
 131 }
 132
 133
 134 void Foam::fvMesh::clearOut()
 135 {
 136     clearGeom();
 137     surfaceInterpolation::clearOut();
 138
 139     clearAddressing();
 140
 141     // Clear mesh motion flux
 142     deleteDemandDrivenData(phiPtr_);
 143
 144     polyMesh::clearOut();
 145 }
 146
 147
 148 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 149
 150 Foam::fvMesh::fvMesh(const IOobject& io)
 151 :
 152     polyMesh(io),
 153     surfaceInterpolation(*this),
 154     fvSchemes(static_cast<const objectRegistry&>(*this)),
 155     fvSolution(static_cast<const objectRegistry&>(*this)),
 156     data(static_cast<const objectRegistry&>(*this)),
 157     boundary_(*this, boundaryMesh()),
 158     lduPtr_(NULL),
 159     curTimeIndex_(time().timeIndex()),
 160     VPtr_(NULL),
 161     V0Ptr_(NULL),
 162     V00Ptr_(NULL),
 163     SfPtr_(NULL),
 164     magSfPtr_(NULL),
 165     CPtr_(NULL),
 166     CfPtr_(NULL),
 167     phiPtr_(NULL)
 168 {
 169     if (debug)
 170     {
 171         Info<< "Constructing fvMesh from IOobject"
 172             << endl;
 173     }
 174
 175     // Check the existance of the cell volumes and read if present
 176     // and set the storage of V00
 177     if (isFile(time().timePath()/"V0"))
 178     {
 179         V0Ptr_ = new DimensionedField<scalar, volMesh>
 180         (
 181             IOobject
 182             (
 183                 "V0",
 184                 time().timeName(),
 185                 *this,
 186                 IOobject::MUST_READ,
 187                 IOobject::NO_WRITE
 188             ),
 189             *this
 190         );
 191
 192         V00();
 193     }
 194
 195     // Check the existance of the mesh fluxes, read if present and set the
 196     // mesh to be moving
 197     if (isFile(time().timePath()/"meshPhi"))
 198     {
 199         phiPtr_ = new surfaceScalarField
 200         (
 201             IOobject
 202             (
 203                 "meshPhi",
 204                 time().timeName(),
 205                 *this,
 206                 IOobject::MUST_READ,
 207                 IOobject::AUTO_WRITE
 208             ),
 209             *this
 210         );
 211
 212         // The mesh is now considered moving so the old-time cell volumes
 213         // will be required for the time derivatives so if they haven't been
 214         // read initialise to the current cell volumes
 215         if (!V0Ptr_)
 216         {
 217             V0Ptr_ = new DimensionedField<scalar, volMesh>
 218             (
 219                 IOobject
 220                 (
 221                     "V0",
 222                     time().timeName(),
 223                     *this,
 224                     IOobject::NO_READ,
 225                     IOobject::NO_WRITE,
 226                     false
 227                 ),
 228                 V()
 229             );
 230         }
 231
 232         moving(true);
 233     }
 234 }
 235
 236
 237 Foam::fvMesh::fvMesh
 238 (
 239     const IOobject& io,
 240     const Xfer<pointField>& points,
 241     const Xfer<faceList>& faces,
 242     const Xfer<labelList>& allOwner,
 243     const Xfer<labelList>& allNeighbour,
 244     const bool syncPar
 245 )
 246 :
 247     polyMesh(io, points, faces, allOwner, allNeighbour, syncPar),
 248     surfaceInterpolation(*this),
 249     fvSchemes(static_cast<const objectRegistry&>(*this)),
 250     fvSolution(static_cast<const objectRegistry&>(*this)),
 251     data(static_cast<const objectRegistry&>(*this)),
 252     boundary_(*this),
 253     lduPtr_(NULL),
 254     curTimeIndex_(time().timeIndex()),
 255     VPtr_(NULL),
 256     V0Ptr_(NULL),
 257     V00Ptr_(NULL),
 258     SfPtr_(NULL),
 259     magSfPtr_(NULL),
 260     CPtr_(NULL),
 261     CfPtr_(NULL),
 262     phiPtr_(NULL)
 263 {
 264     if (debug)
 265     {
 266         Info<< "Constructing fvMesh from components" << endl;
 267     }
 268 }
 269
 270
 271 Foam::fvMesh::fvMesh
 272 (
 273     const IOobject& io,
 274     const Xfer<pointField>& points,
 275     const Xfer<faceList>& faces,
 276     const Xfer<cellList>& cells,
 277     const bool syncPar
 278 )
 279 :
 280     polyMesh(io, points, faces, cells, syncPar),
 281     surfaceInterpolation(*this),
 282     fvSchemes(static_cast<const objectRegistry&>(*this)),
 283     fvSolution(static_cast<const objectRegistry&>(*this)),
 284     data(static_cast<const objectRegistry&>(*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     polyMesh::readUpdateState state = polyMesh::readUpdate();
 362
 363     if (state == polyMesh::TOPO_PATCH_CHANGE)
 364     {
 365         if (debug)
 366         {
 367             Info<< "Boundary and topological update" << endl;
 368         }
 369
 370         boundary_.readUpdate(boundaryMesh());
 371
 372         clearOut();
 373
 374     }
 375     else if (state == polyMesh::TOPO_CHANGE)
 376     {
 377         if (debug)
 378         {
 379             Info<< "Topological update" << endl;
 380         }
 381
 382         clearOut();
 383     }
 384     else if (state == polyMesh::POINTS_MOVED)
 385     {
 386         if (debug)
 387         {
 388             Info<< "Point motion update" << endl;
 389         }
 390
 391         clearGeom();
 392     }
 393     else
 394     {
 395         if (debug)
 396         {
 397             Info<< "No update" << endl;
 398         }
 399     }
 400
 401     return state;
 402 }
 403
 404
 405 const Foam::fvBoundaryMesh& Foam::fvMesh::boundary() const
 406 {
 407     return boundary_;
 408 }
 409
 410
 411 const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
 412 {
 413     if (!lduPtr_)
 414     {
 415         lduPtr_ = new fvMeshLduAddressing(*this);
 416     }
 417
 418     return *lduPtr_;
 419 }
 420
 421
 422 void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 423 {
 424     // Create a mapper
 425     const fvMeshMapper mapper(*this, meshMap);
 426
 427     // Map all the volFields in the objectRegistry
 428     MapGeometricFields<scalar, fvPatchField, fvMeshMapper, volMesh>
 429     (mapper);
 430     MapGeometricFields<vector, fvPatchField, fvMeshMapper, volMesh>
 431     (mapper);
 432     MapGeometricFields<sphericalTensor, fvPatchField, fvMeshMapper, volMesh>
 433     (mapper);
 434     MapGeometricFields<symmTensor, fvPatchField, fvMeshMapper, volMesh>
 435     (mapper);
 436     MapGeometricFields<tensor, fvPatchField, fvMeshMapper, volMesh>
 437     (mapper);
 438
 439     // Map all the surfaceFields in the objectRegistry
 440     MapGeometricFields<scalar, fvsPatchField, fvMeshMapper, surfaceMesh>
 441     (mapper);
 442     MapGeometricFields<vector, fvsPatchField, fvMeshMapper, surfaceMesh>
 443     (mapper);
 444     MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 445     (mapper);
 446     MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 447     (mapper);
 448     MapGeometricFields<tensor, fvsPatchField, fvMeshMapper, surfaceMesh>
 449     (mapper);
 450
 451     // Map all the clouds in the objectRegistry
 452     mapClouds(*this, meshMap);
 453
 454
 455     const labelList& cellMap = meshMap.cellMap();
 456
 457     // Map the old volume. Just map to new cell labels.
 458     if (V0Ptr_)
 459     {
 460         scalarField& V0 = *V0Ptr_;
 461
 462         scalarField savedV0(V0);
 463         V0.setSize(nCells());
 464
 465         forAll(V0, i)
 466         {
 467             if (cellMap[i] > -1)
 468             {
 469                 V0[i] = savedV0[cellMap[i]];
 470             }
 471             else
 472             {
 473                 V0[i] = 0.0;
 474             }
 475         }
 476     }
 477
 478     // Map the old-old volume. Just map to new cell labels.
 479     if (V00Ptr_)
 480     {
 481         scalarField& V00 = *V00Ptr_;
 482
 483         scalarField savedV00(V00);
 484         V00.setSize(nCells());
 485
 486         forAll(V00, i)
 487         {
 488             if (cellMap[i] > -1)
 489             {
 490                 V00[i] = savedV00[cellMap[i]];
 491             }
 492             else
 493             {
 494                 V00[i] = 0.0;
 495             }
 496         }
 497     }
 498 }
 499
 500
 501 // Temporary helper function to call move points on
 502 // MeshObjects
 503 template<class Type>
 504 void MeshObjectMovePoints(const Foam::fvMesh& mesh)
 505 {
 506     if (mesh.thisDb().foundObject<Type>(Type::typeName))
 507     {
 508         const_cast<Type&>
 509         (
 510             mesh.thisDb().lookupObject<Type>
 511             (
 512                 Type::typeName
 513             )
 514         ).movePoints();
 515     }
 516 }
 517
 518
 519 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 520 {
 521     // Grab old time volumes if the time has been incremented
 522     if (curTimeIndex_ < time().timeIndex())
 523     {
 524         if (V00Ptr_ && V0Ptr_)
 525         {
 526             *V00Ptr_ = *V0Ptr_;
 527         }
 528
 529         if (V0Ptr_)
 530         {
 531             *V0Ptr_ = V();
 532         }
 533         else
 534         {
 535             V0Ptr_ = new DimensionedField<scalar, volMesh>
 536             (
 537                 IOobject
 538                 (
 539                     "V0",
 540                     time().timeName(),
 541                     *this,
 542                     IOobject::NO_READ,
 543                     IOobject::NO_WRITE
 544                 ),
 545                 V()
 546             );
 547         }
 548
 549         curTimeIndex_ = time().timeIndex();
 550     }
 551
 552
 553     // delete out of date geometrical information
 554     clearGeomNotOldVol();
 555
 556
 557     if (!phiPtr_)
 558     {
 559         // Create mesh motion flux
 560         phiPtr_ = new surfaceScalarField
 561         (
 562             IOobject
 563             (
 564                 "meshPhi",
 565                 this->time().timeName(),
 566                 *this,
 567                 IOobject::NO_READ,
 568                 IOobject::AUTO_WRITE
 569             ),
 570             *this,
 571             dimVolume/dimTime
 572         );
 573     }
 574     else
 575     {
 576         // Grab old time mesh motion fluxes if the time has been incremented
 577         if (phiPtr_->timeIndex() != time().timeIndex())
 578         {
 579             phiPtr_->oldTime();
 580         }
 581     }
 582
 583     surfaceScalarField& phi = *phiPtr_;
 584
 585     // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
 586
 587     scalar rDeltaT = 1.0/time().deltaTValue();
 588
 589     tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
 590     scalarField& sweptVols = tsweptVols();
 591
 592     phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
 593     phi.internalField() *= rDeltaT;
 594
 595     const fvPatchList& patches = boundary();
 596
 597     forAll(patches, patchI)
 598     {
 599         phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
 600         phi.boundaryField()[patchI] *= rDeltaT;
 601     }
 602
 603     boundary_.movePoints();
 604     surfaceInterpolation::movePoints();
 605
 606
 607     // Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
 608     // movePoints function.
 609     MeshObjectMovePoints<volPointInterpolation>(*this);
 610     MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
 611     MeshObjectMovePoints<leastSquaresVectors>(*this);
 612     MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
 613     MeshObjectMovePoints<CentredFitData<quadraticFitPolynomial> >(*this);
 614     MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
 615     MeshObjectMovePoints<skewCorrectionVectors>(*this);
 616     //MeshObjectMovePoints<quadraticFitSnGradData>(*this);
 617
 618     return tsweptVols;
 619 }
 620
 621
 622 void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
 623 {
 624     // Update polyMesh. This needs to keep volume existent!
 625     polyMesh::updateMesh(mpm);
 626
 627     // Clear the sliced fields
 628     clearGeomNotOldVol();
 629
 630     // Map all fields using current (i.e. not yet mapped) volume
 631     mapFields(mpm);
 632
 633     // Clear the current volume and other geometry factors
 634     surfaceInterpolation::clearOut();
 635
 636     clearAddressing();
 637
 638     // handleMorph() should also clear out the surfaceInterpolation.
 639     // This is a temporary solution
 640     surfaceInterpolation::movePoints();
 641 }
 642
 643
 644 bool Foam::fvMesh::writeObjects
 645 (
 646     IOstream::streamFormat fmt,
 647     IOstream::versionNumber ver,
 648     IOstream::compressionType cmp
 649 ) const
 650 {
 651     return polyMesh::writeObject(fmt, ver, cmp);
 652 }
 653
 654
 655 //- Write mesh using IO settings from the time
 656 bool Foam::fvMesh::write() const
 657 {
 658     return polyMesh::write();
 659 }
 660
 661
 662 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 663
 664 bool Foam::fvMesh::operator!=(const fvMesh& bm) const
 665 {
 666     return &bm != this;
 667 }
 668
 669
 670 bool Foam::fvMesh::operator==(const fvMesh& bm) const
 671 {
 672     return &bm == this;
 673 }
 674
 675
 676 // ************************************************************************* //