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