src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.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 Description
  26     A subset of mesh faces.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "faceZone.H"
  31 #include "addToRunTimeSelectionTable.H"
  32 #include "faceZoneMesh.H"
  33 #include "polyMesh.H"
  34 #include "primitiveMesh.H"
  35 #include "demandDrivenData.H"
  36 #include "mapPolyMesh.H"
  37 #include "syncTools.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43     defineTypeNameAndDebug(faceZone, 0);
  44     defineRunTimeSelectionTable(faceZone, dictionary);
  45     addToRunTimeSelectionTable(faceZone, faceZone, dictionary);
  46 }
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 void Foam::faceZone::calcFaceZonePatch() const
  51 {
  52     if (debug)
  53     {
  54         Info<< "void faceZone::calcFaceZonePatch() const : "
  55             << "Calculating primitive patch"
  56             << endl;
  57     }
  58
  59     if (patchPtr_)
  60     {
  61         FatalErrorIn
  62         (
  63             "void faceZone::calcFaceZonePatch() const"
  64         )   << "primitive face zone patch already calculated"
  65             << abort(FatalError);
  66     }
  67
  68     patchPtr_ =
  69         new primitiveFacePatch
  70         (
  71             faceList(size()),
  72             zoneMesh().mesh().points()
  73         );
  74
  75     primitiveFacePatch& patch = *patchPtr_;
  76
  77     const faceList& f = zoneMesh().mesh().faces();
  78
  79     const labelList& addr = *this;
  80     const boolList& flip = flipMap();
  81
  82     forAll (addr, faceI)
  83     {
  84         if (flip[faceI])
  85         {
  86             patch[faceI] = f[addr[faceI]].reverseFace();
  87         }
  88         else
  89         {
  90             patch[faceI] = f[addr[faceI]];
  91         }
  92     }
  93
  94     if (debug)
  95     {
  96         Info<< "void faceZone::calcFaceZonePatch() const : "
  97             << "Finished calculating primitive patch"
  98             << endl;
  99     }
 100 }
 101
 102
 103 const Foam::Map<Foam::label>& Foam::faceZone::faceLookupMap() const
 104 {
 105     if (!faceLookupMapPtr_)
 106     {
 107         calcFaceLookupMap();
 108     }
 109
 110     return *faceLookupMapPtr_;
 111 }
 112
 113
 114 void Foam::faceZone::calcFaceLookupMap() const
 115 {
 116     if (debug)
 117     {
 118         Info<< "void faceZone::calcFaceLookupMap() const : "
 119             << "Calculating face lookup map"
 120             << endl;
 121     }
 122
 123     if (faceLookupMapPtr_)
 124     {
 125         FatalErrorIn
 126         (
 127             "void faceZone::calcFaceLookupMap() const"
 128         )   << "face lookup map already calculated"
 129             << abort(FatalError);
 130     }
 131
 132     const labelList& addr = *this;
 133
 134     faceLookupMapPtr_ = new Map<label>(2*addr.size());
 135     Map<label>& flm = *faceLookupMapPtr_;
 136
 137     forAll (addr, faceI)
 138     {
 139         flm.insert(addr[faceI], faceI);
 140     }
 141
 142     if (debug)
 143     {
 144         Info<< "void faceZone::calcFaceLookupMap() const : "
 145             << "Finished calculating face lookup map"
 146             << endl;
 147     }
 148 }
 149
 150
 151 void Foam::faceZone::calcCellLayers() const
 152 {
 153     if (debug)
 154     {
 155         Info<< "void Foam::faceZone::calcCellLayers() const : "
 156             << "calculating master cells"
 157             << endl;
 158     }
 159
 160     // It is an error to attempt to recalculate edgeCells
 161     // if the pointer is already set
 162     if (masterCellsPtr_ || slaveCellsPtr_)
 163     {
 164         FatalErrorIn("void faceZone::calcCellLayers() const")
 165             << "cell layers already calculated"
 166             << abort(FatalError);
 167     }
 168     else
 169     {
 170         // Go through all the faces in the master zone.  Choose the
 171         // master or slave cell based on the face flip
 172
 173         const labelList& own = zoneMesh().mesh().faceOwner();
 174         const labelList& nei = zoneMesh().mesh().faceNeighbour();
 175
 176         const labelList& mf = *this;
 177
 178         const boolList& faceFlip = flipMap();
 179
 180         masterCellsPtr_ = new labelList(mf.size());
 181         labelList& mc = *masterCellsPtr_;
 182
 183         slaveCellsPtr_ = new labelList(mf.size());
 184         labelList& sc = *slaveCellsPtr_;
 185
 186         forAll (mf, faceI)
 187         {
 188             if (!faceFlip[faceI])
 189             {
 190                 // Face is oriented correctly, no flip needed
 191                 mc[faceI] = nei[mf[faceI]];
 192                 sc[faceI] = own[mf[faceI]];
 193             }
 194             else
 195             {
 196                 mc[faceI] = own[mf[faceI]];
 197                 sc[faceI] = nei[mf[faceI]];
 198             }
 199         }
 200         //Info << "masterCells: " << mc << endl;
 201         //Info << "slaveCells: " << sc << endl;
 202     }
 203 }
 204
 205
 206 void Foam::faceZone::checkAddressing() const
 207 {
 208     if (size() != flipMap_.size())
 209     {
 210         FatalErrorIn("void Foam::faceZone::checkAddressing() const")
 211             << "Different sizes of the addressing and flipMap arrays.  "
 212             << "Size of addressing: " << size()
 213             << " size of flip map: " << flipMap_.size()
 214             << abort(FatalError);
 215     }
 216 }
 217
 218
 219 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 220
 221 // Construct from components
 222 Foam::faceZone::faceZone
 223 (
 224     const word& name,
 225     const labelList& addr,
 226     const boolList& fm,
 227     const label index,
 228     const faceZoneMesh& zm
 229 )
 230 :
 231     labelList(addr),
 232     name_(name),
 233     flipMap_(fm),
 234     index_(index),
 235     zoneMesh_(zm),
 236     patchPtr_(NULL),
 237     masterCellsPtr_(NULL),
 238     slaveCellsPtr_(NULL),
 239     mePtr_(NULL),
 240     faceLookupMapPtr_(NULL)
 241 {
 242     checkAddressing();
 243 }
 244
 245
 246 Foam::faceZone::faceZone
 247 (
 248     const word& name,
 249     const Xfer<labelList>& addr,
 250     const Xfer<boolList>& fm,
 251     const label index,
 252     const faceZoneMesh& zm
 253 )
 254 :
 255     labelList(addr),
 256     name_(name),
 257     flipMap_(fm),
 258     index_(index),
 259     zoneMesh_(zm),
 260     patchPtr_(NULL),
 261     masterCellsPtr_(NULL),
 262     slaveCellsPtr_(NULL),
 263     mePtr_(NULL),
 264     faceLookupMapPtr_(NULL)
 265 {
 266     checkAddressing();
 267 }
 268
 269
 270 // Construct from dictionary
 271 Foam::faceZone::faceZone
 272 (
 273     const word& name,
 274     const dictionary& dict,
 275     const label index,
 276     const faceZoneMesh& zm
 277 )
 278 :
 279     labelList(dict.lookup("faceLabels")),
 280     name_(name),
 281     flipMap_(dict.lookup("flipMap")),
 282     index_(index),
 283     zoneMesh_(zm),
 284     patchPtr_(NULL),
 285     masterCellsPtr_(NULL),
 286     slaveCellsPtr_(NULL),
 287     mePtr_(NULL),
 288     faceLookupMapPtr_(NULL)
 289 {
 290     checkAddressing();
 291 }
 292
 293
 294 // Construct given the original zone and resetting the
 295 // face list and zone mesh information
 296 Foam::faceZone::faceZone
 297 (
 298     const faceZone& fz,
 299     const labelList& addr,
 300     const boolList& fm,
 301     const label index,
 302     const faceZoneMesh& zm
 303 )
 304 :
 305     labelList(addr),
 306     name_(fz.name()),
 307     flipMap_(fm),
 308     index_(index),
 309     zoneMesh_(zm),
 310     patchPtr_(NULL),
 311     masterCellsPtr_(NULL),
 312     slaveCellsPtr_(NULL),
 313     mePtr_(NULL),
 314     faceLookupMapPtr_(NULL)
 315 {
 316     checkAddressing();
 317 }
 318
 319
 320 Foam::faceZone::faceZone
 321 (
 322     const faceZone& fz,
 323     const Xfer<labelList>& addr,
 324     const Xfer<boolList>& fm,
 325     const label index,
 326     const faceZoneMesh& zm
 327 )
 328 :
 329     labelList(addr),
 330     name_(fz.name()),
 331     flipMap_(fm),
 332     index_(index),
 333     zoneMesh_(zm),
 334     patchPtr_(NULL),
 335     masterCellsPtr_(NULL),
 336     slaveCellsPtr_(NULL),
 337     mePtr_(NULL),
 338     faceLookupMapPtr_(NULL)
 339 {
 340     checkAddressing();
 341 }
 342
 343
 344 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 345
 346 Foam::faceZone::~faceZone()
 347 {
 348     clearAddressing();
 349 }
 350
 351
 352 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 353
 354 Foam::label Foam::faceZone::whichFace(const label globalFaceID) const
 355 {
 356     const Map<label>& flm = faceLookupMap();
 357
 358     Map<label>::const_iterator flmIter = flm.find(globalFaceID);
 359
 360     if (flmIter == flm.end())
 361     {
 362         return -1;
 363     }
 364     else
 365     {
 366         return flmIter();
 367     }
 368 }
 369
 370
 371 const Foam::faceZoneMesh& Foam::faceZone::zoneMesh() const
 372 {
 373     return zoneMesh_;
 374 }
 375
 376
 377 const Foam::primitiveFacePatch& Foam::faceZone::operator()() const
 378 {
 379     if (!patchPtr_)
 380     {
 381         calcFaceZonePatch();
 382     }
 383
 384     return *patchPtr_;
 385 }
 386
 387
 388 const Foam::labelList& Foam::faceZone::masterCells() const
 389 {
 390     if (!masterCellsPtr_)
 391     {
 392         calcCellLayers();
 393     }
 394
 395     return *masterCellsPtr_;
 396 }
 397
 398
 399 const Foam::labelList& Foam::faceZone::slaveCells() const
 400 {
 401     if (!slaveCellsPtr_)
 402     {
 403         calcCellLayers();
 404     }
 405
 406     return *slaveCellsPtr_;
 407 }
 408
 409
 410 const Foam::labelList& Foam::faceZone::meshEdges() const
 411 {
 412     if (!mePtr_)
 413     {
 414         //labelList faceCells(size());
 415         //
 416         //const labelList& own = zoneMesh().mesh().faceOwner();
 417         //
 418         //const labelList& faceLabels = *this;
 419         //
 420         //forAll (faceCells, faceI)
 421         //{
 422         //    faceCells[faceI] = own[faceLabels[faceI]];
 423         //}
 424         //
 425         //mePtr_ =
 426         //    new labelList
 427         //    (
 428         //        operator()().meshEdges
 429         //        (
 430         //            zoneMesh().mesh().edges(),
 431         //            zoneMesh().mesh().cellEdges(),
 432         //            faceCells
 433         //        )
 434         //    );
 435
 436         mePtr_ =
 437             new labelList
 438             (
 439                 operator()().meshEdges
 440                 (
 441                     zoneMesh().mesh().edges(),
 442                     zoneMesh().mesh().pointEdges()
 443                 )
 444             );
 445     }
 446
 447     return *mePtr_;
 448 }
 449
 450
 451 void Foam::faceZone::clearAddressing()
 452 {
 453     deleteDemandDrivenData(patchPtr_);
 454
 455     deleteDemandDrivenData(masterCellsPtr_);
 456     deleteDemandDrivenData(slaveCellsPtr_);
 457
 458     deleteDemandDrivenData(mePtr_);
 459     deleteDemandDrivenData(faceLookupMapPtr_);
 460 }
 461
 462
 463 void Foam::faceZone::resetAddressing
 464 (
 465     const labelList& addr,
 466     const boolList& flipMap
 467 )
 468 {
 469     clearAddressing();
 470     labelList::operator=(addr);
 471     flipMap_ = flipMap;
 472 }
 473
 474
 475 void Foam::faceZone::updateMesh(const mapPolyMesh& mpm)
 476 {
 477     clearAddressing();
 478
 479     labelList newAddressing(size());
 480     boolList newFlipMap(flipMap_.size());
 481     label nFaces = 0;
 482
 483     const labelList& faceMap = mpm.reverseFaceMap();
 484
 485     forAll(*this, i)
 486     {
 487         label faceI = operator[](i);
 488
 489         if (faceMap[faceI] >= 0)
 490         {
 491             newAddressing[nFaces] = faceMap[faceI];
 492             newFlipMap[nFaces] = flipMap_[i];       // Keep flip map.
 493             nFaces++;
 494         }
 495     }
 496
 497     newAddressing.setSize(nFaces);
 498     newFlipMap.setSize(nFaces);
 499
 500     transfer(newAddressing);
 501     flipMap_.transfer(newFlipMap);
 502 }
 503
 504
 505 bool Foam::faceZone::checkDefinition(const bool report) const
 506 {
 507     const labelList& addr = *this;
 508
 509     bool boundaryError = false;
 510
 511     forAll(addr, i)
 512     {
 513         if (addr[i] < 0 || addr[i] >= zoneMesh().mesh().faces().size())
 514         {
 515             boundaryError = true;
 516
 517             if (report)
 518             {
 519                 SeriousErrorIn
 520                 (
 521                     "bool faceZone::checkDefinition("
 522                     "const bool report) const"
 523                 )   << "Zone " << name()
 524                     << " contains invalid face label " << addr[i] << nl
 525                     << "Valid face labels are 0.."
 526                     << zoneMesh().mesh().faces().size()-1 << endl;
 527             }
 528         }
 529     }
 530     return boundaryError;
 531 }
 532
 533
 534 bool Foam::faceZone::checkParallelSync(const bool report) const
 535 {
 536     const polyMesh& mesh = zoneMesh().mesh();
 537     const polyBoundaryMesh& bm = mesh.boundaryMesh();
 538
 539     bool boundaryError = false;
 540
 541
 542     // Check that zone faces are synced
 543     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 544
 545     {
 546         boolList neiZoneFace(mesh.nFaces()-mesh.nInternalFaces(), false);
 547         boolList neiZoneFlip(mesh.nFaces()-mesh.nInternalFaces(), false);
 548         forAll(*this, i)
 549         {
 550             label faceI = operator[](i);
 551
 552             if (!mesh.isInternalFace(faceI))
 553             {
 554                 neiZoneFace[faceI-mesh.nInternalFaces()] = true;
 555                 neiZoneFlip[faceI-mesh.nInternalFaces()] = flipMap()[i];
 556             }
 557         }
 558         boolList myZoneFace(neiZoneFace);
 559         syncTools::swapBoundaryFaceList(mesh, neiZoneFace, false);
 560         boolList myZoneFlip(neiZoneFlip);
 561         syncTools::swapBoundaryFaceList(mesh, neiZoneFlip, false);
 562
 563         forAll(*this, i)
 564         {
 565             label faceI = operator[](i);
 566
 567             label patchI = bm.whichPatch(faceI);
 568
 569             if (patchI != -1 && bm[patchI].coupled())
 570             {
 571                 label bFaceI = faceI-mesh.nInternalFaces();
 572
 573                 // Check face in zone on both sides
 574                 if (myZoneFace[bFaceI] != neiZoneFace[bFaceI])
 575                 {
 576                     boundaryError = true;
 577
 578                     if (report)
 579                     {
 580                         Pout<< " ***Problem with faceZone " << index()
 581                             << " named " << name()
 582                             << ". Face " << faceI
 583                             << " on coupled patch "
 584                             << bm[patchI].name()
 585                             << " is not consistent with its coupled neighbour."
 586                             << endl;
 587                     }
 588                 }
 589
 590                 // Flip state should be opposite.
 591                 if (myZoneFlip[bFaceI] == neiZoneFlip[bFaceI])
 592                 {
 593                     boundaryError = true;
 594
 595                     if (report)
 596                     {
 597                         Pout<< " ***Problem with faceZone " << index()
 598                             << " named " << name()
 599                             << ". Face " << faceI
 600                             << " on coupled patch "
 601                             << bm[patchI].name()
 602                             << " does not have consistent flipMap"
 603                             << " across coupled faces."
 604                             << endl;
 605                     }
 606                 }
 607             }
 608         }
 609     }
 610
 611     return returnReduce(boundaryError, orOp<bool>());
 612 }
 613
 614
 615 void Foam::faceZone::movePoints(const pointField& p)
 616 {
 617     if (patchPtr_)
 618     {
 619         patchPtr_->movePoints(p);
 620     }
 621 }
 622
 623 void Foam::faceZone::write(Ostream& os) const
 624 {
 625     os  << nl << name()
 626         << nl << static_cast<const labelList&>(*this)
 627         << nl << flipMap();
 628 }
 629
 630
 631 void Foam::faceZone::writeDict(Ostream& os) const
 632 {
 633     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
 634         << "    type " << type() << token::END_STATEMENT << nl;
 635
 636     writeEntry("faceLabels", os);
 637     flipMap().writeEntry("flipMap", os);
 638
 639     os  << token::END_BLOCK << endl;
 640 }
 641
 642
 643 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 644
 645 Foam::Ostream& Foam::operator<<(Ostream& os, const faceZone& p)
 646 {
 647     p.write(os);
 648     os.check("Ostream& operator<<(Ostream& f, const faceZone& p");
 649     return os;
 650 }
 651
 652
 653 // ************************************************************************* //