src/finiteVolume/fvMesh/fvMesh.C

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