src/dynamicMesh/slidingInterface/slidingInterfaceAttachedAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "slidingInterface.H"
  27 #include "polyMesh.H"
  28 #include "mapPolyMesh.H"
  29 #include "polyTopoChanger.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 void Foam::slidingInterface::calcAttachedAddressing() const
  34 {
  35     if (debug)
  36     {
  37         Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
  38             << " for object " << name() << " : "
  39             << "Calculating zone face-cell addressing."
  40             << endl;
  41     }
  42
  43     if (!attached_)
  44     {
  45         // Clear existing addressing
  46         clearAttachedAddressing();
  47
  48         const polyMesh& mesh = topoChanger().mesh();
  49         const labelList& own = mesh.faceOwner();
  50         const labelList& nei = mesh.faceNeighbour();
  51         const faceZoneMesh& faceZones = mesh.faceZones();
  52
  53         // Master side
  54
  55         const primitiveFacePatch& masterPatch =
  56             faceZones[masterFaceZoneID_.index()]();
  57
  58         const labelList& masterPatchFaces =
  59             faceZones[masterFaceZoneID_.index()];
  60
  61         const boolList& masterFlip =
  62             faceZones[masterFaceZoneID_.index()].flipMap();
  63
  64         masterFaceCellsPtr_ = new labelList(masterPatchFaces.size());
  65         labelList& mfc = *masterFaceCellsPtr_;
  66
  67         forAll(masterPatchFaces, faceI)
  68         {
  69             if (masterFlip[faceI])
  70             {
  71                 mfc[faceI] = nei[masterPatchFaces[faceI]];
  72             }
  73             else
  74             {
  75                 mfc[faceI] = own[masterPatchFaces[faceI]];
  76             }
  77         }
  78
  79         // Slave side
  80
  81         const primitiveFacePatch& slavePatch =
  82             faceZones[slaveFaceZoneID_.index()]();
  83
  84         const labelList& slavePatchFaces =
  85             faceZones[slaveFaceZoneID_.index()];
  86
  87         const boolList& slaveFlip =
  88             faceZones[slaveFaceZoneID_.index()].flipMap();
  89
  90         slaveFaceCellsPtr_ = new labelList(slavePatchFaces.size());
  91         labelList& sfc = *slaveFaceCellsPtr_;
  92
  93         forAll(slavePatchFaces, faceI)
  94         {
  95             if (slaveFlip[faceI])
  96             {
  97                 sfc[faceI] = nei[slavePatchFaces[faceI]];
  98             }
  99             else
 100             {
 101                 sfc[faceI] = own[slavePatchFaces[faceI]];
 102             }
 103         }
 104
 105         // Check that the addressing is valid
 106         if (min(mfc) < 0 || min(sfc) < 0)
 107         {
 108             if (debug)
 109             {
 110                 forAll(mfc, faceI)
 111                 {
 112                     if (mfc[faceI] < 0)
 113                     {
 114                         Pout<< "No cell next to master patch face " << faceI
 115                             << ".  Global face no: " << mfc[faceI]
 116                             << " own: " << own[masterPatchFaces[faceI]]
 117                             << " nei: " << nei[masterPatchFaces[faceI]]
 118                             << " flip: " << masterFlip[faceI] << endl;
 119                     }
 120                 }
 121
 122                 forAll(sfc, faceI)
 123                 {
 124                     if (sfc[faceI] < 0)
 125                     {
 126                         Pout<< "No cell next to slave patch face " << faceI
 127                             << ".  Global face no: " << sfc[faceI]
 128                             << " own: " << own[slavePatchFaces[faceI]]
 129                             << " nei: " << nei[slavePatchFaces[faceI]]
 130                             << " flip: " << slaveFlip[faceI] << endl;
 131                     }
 132                 }
 133             }
 134
 135             FatalErrorIn
 136             (
 137                 "void slidingInterface::calcAttachedAddressing()"
 138                 "const"
 139             )   << "Error is zone face-cell addressing.  Probable error in "
 140                 << "decoupled mesh or sliding interface definition."
 141                 << abort(FatalError);
 142         }
 143
 144         // Calculate stick-out faces
 145         const labelListList& pointFaces = mesh.pointFaces();
 146
 147         // Master side
 148         labelHashSet masterStickOutFaceMap
 149         (
 150             primitiveMesh::facesPerCell_*(masterPatch.size())
 151         );
 152
 153         const labelList& masterMeshPoints = masterPatch.meshPoints();
 154
 155         forAll(masterMeshPoints, pointI)
 156         {
 157             const labelList& curFaces = pointFaces[masterMeshPoints[pointI]];
 158
 159             forAll(curFaces, faceI)
 160             {
 161                 // Check if the face belongs to the master face zone;
 162                 // if not add it
 163                 if
 164                 (
 165                     faceZones.whichZone(curFaces[faceI])
 166                  != masterFaceZoneID_.index()
 167                 )
 168                 {
 169                     masterStickOutFaceMap.insert(curFaces[faceI]);
 170                 }
 171             }
 172         }
 173
 174         masterStickOutFacesPtr_ = new labelList(masterStickOutFaceMap.toc());
 175
 176         // Slave side
 177         labelHashSet slaveStickOutFaceMap
 178         (
 179             primitiveMesh::facesPerCell_*(slavePatch.size())
 180         );
 181
 182         const labelList& slaveMeshPoints = slavePatch.meshPoints();
 183
 184         forAll(slaveMeshPoints, pointI)
 185         {
 186             const labelList& curFaces = pointFaces[slaveMeshPoints[pointI]];
 187
 188             forAll(curFaces, faceI)
 189             {
 190                 // Check if the face belongs to the slave face zone;
 191                 // if not add it
 192                 if
 193                 (
 194                     faceZones.whichZone(curFaces[faceI])
 195                  != slaveFaceZoneID_.index()
 196                 )
 197                 {
 198                     slaveStickOutFaceMap.insert(curFaces[faceI]);
 199                 }
 200             }
 201         }
 202
 203         slaveStickOutFacesPtr_ = new labelList(slaveStickOutFaceMap.toc());
 204
 205
 206         // Retired point addressing does not exist at this stage.
 207         // It will be filled when the interface is coupled.
 208         retiredPointMapPtr_ =
 209             new Map<label>
 210             (
 211                 2*faceZones[slaveFaceZoneID_.index()]().nPoints()
 212             );
 213
 214         // Ditto for cut point edge map.  This is a rough guess of its size
 215         //
 216         cutPointEdgePairMapPtr_ =
 217             new Map<Pair<edge> >
 218             (
 219                 faceZones[slaveFaceZoneID_.index()]().nEdges()
 220             );
 221     }
 222     else
 223     {
 224         FatalErrorIn
 225         (
 226             "void slidingInterface::calcAttachedAddressing() const"
 227         )   << "The interface is attached.  The zone face-cell addressing "
 228             << "cannot be assembled for object " << name()
 229             << abort(FatalError);
 230     }
 231
 232     if (debug)
 233     {
 234         Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
 235             << " for object " << name() << " : "
 236             << "Finished calculating zone face-cell addressing."
 237             << endl;
 238     }
 239 }
 240
 241
 242 void Foam::slidingInterface::clearAttachedAddressing() const
 243 {
 244     deleteDemandDrivenData(masterFaceCellsPtr_);
 245     deleteDemandDrivenData(slaveFaceCellsPtr_);
 246
 247     deleteDemandDrivenData(masterStickOutFacesPtr_);
 248     deleteDemandDrivenData(slaveStickOutFacesPtr_);
 249
 250     deleteDemandDrivenData(retiredPointMapPtr_);
 251     deleteDemandDrivenData(cutPointEdgePairMapPtr_);
 252 }
 253
 254
 255 void Foam::slidingInterface::renumberAttachedAddressing
 256 (
 257     const mapPolyMesh& m
 258 ) const
 259 {
 260     // Get reference to reverse cell renumbering
 261     // The renumbering map is needed the other way around, i.e. giving
 262     // the new cell number for every old cell next to the interface.
 263     const labelList& reverseCellMap = m.reverseCellMap();
 264
 265     const labelList& mfc = masterFaceCells();
 266     const labelList& sfc = slaveFaceCells();
 267
 268     // Master side
 269     labelList* newMfcPtr = new labelList(mfc.size(), -1);
 270     labelList& newMfc = *newMfcPtr;
 271
 272     const labelList& mfzRenumber =
 273         m.faceZoneFaceMap()[masterFaceZoneID_.index()];
 274
 275     forAll(mfc, faceI)
 276     {
 277         label newCellI = reverseCellMap[mfc[mfzRenumber[faceI]]];
 278
 279         if (newCellI >= 0)
 280         {
 281             newMfc[faceI] = newCellI;
 282         }
 283     }
 284
 285     // Slave side
 286     labelList* newSfcPtr = new labelList(sfc.size(), -1);
 287     labelList& newSfc = *newSfcPtr;
 288
 289     const labelList& sfzRenumber =
 290         m.faceZoneFaceMap()[slaveFaceZoneID_.index()];
 291
 292     forAll(sfc, faceI)
 293     {
 294         label newCellI = reverseCellMap[sfc[sfzRenumber[faceI]]];
 295
 296         if (newCellI >= 0)
 297         {
 298             newSfc[faceI] = newCellI;
 299         }
 300     }
 301
 302     if (debug)
 303     {
 304         // Check if all the mapped cells are live
 305         if (min(newMfc) < 0 || min(newSfc) < 0)
 306         {
 307             FatalErrorIn
 308             (
 309                 "void slidingInterface::renumberAttachedAddressing("
 310                 "const mapPolyMesh& m) const"
 311             )   << "Error in cell renumbering for object " << name()
 312                 << ".  Some of master cells next "
 313                 << "to the interface have been removed."
 314                 << abort(FatalError);
 315         }
 316     }
 317
 318     // Renumber stick-out faces
 319
 320     const labelList& reverseFaceMap = m.reverseFaceMap();
 321
 322     // Master side
 323     const labelList& msof = masterStickOutFaces();
 324
 325     labelList* newMsofPtr = new labelList(msof.size(), -1);
 326     labelList& newMsof = *newMsofPtr;
 327
 328     forAll(msof, faceI)
 329     {
 330         label newFaceI = reverseFaceMap[msof[faceI]];
 331
 332         if (newFaceI >= 0)
 333         {
 334             newMsof[faceI] = newFaceI;
 335         }
 336     }
 337 //     Pout<< "newMsof: " << newMsof << endl;
 338     // Slave side
 339     const labelList& ssof = slaveStickOutFaces();
 340
 341     labelList* newSsofPtr = new labelList(ssof.size(), -1);
 342     labelList& newSsof = *newSsofPtr;
 343
 344     forAll(ssof, faceI)
 345     {
 346         label newFaceI = reverseFaceMap[ssof[faceI]];
 347
 348         if (newFaceI >= 0)
 349         {
 350             newSsof[faceI] = newFaceI;
 351         }
 352     }
 353 //     Pout<< "newSsof: " << newSsof << endl;
 354     if (debug)
 355     {
 356         // Check if all the mapped cells are live
 357         if (min(newMsof) < 0 || min(newSsof) < 0)
 358         {
 359             FatalErrorIn
 360             (
 361                 "void slidingInterface::renumberAttachedAddressing("
 362                 "const mapPolyMesh& m) const"
 363             )   << "Error in face renumbering for object " << name()
 364                 << ".  Some of stick-out next "
 365                 << "to the interface have been removed."
 366                 << abort(FatalError);
 367         }
 368     }
 369
 370     // Renumber the retired point map. Need to take a copy!
 371     const Map<label> rpm = retiredPointMap();
 372
 373     Map<label>* newRpmPtr = new Map<label>(rpm.size());
 374     Map<label>& newRpm = *newRpmPtr;
 375
 376     const labelList rpmToc = rpm.toc();
 377
 378     // Get reference to point renumbering
 379     const labelList& reversePointMap = m.reversePointMap();
 380
 381     label key, value;
 382
 383     forAll(rpmToc, rpmTocI)
 384     {
 385         key = reversePointMap[rpmToc[rpmTocI]];
 386
 387         value = reversePointMap[rpm.find(rpmToc[rpmTocI])()];
 388
 389         if (debug)
 390         {
 391             // Check if all the mapped cells are live
 392             if (key < 0 || value < 0)
 393             {
 394                 FatalErrorIn
 395                 (
 396                     "void slidingInterface::renumberAttachedAddressing("
 397                     "const mapPolyMesh& m) const"
 398                 )   << "Error in retired point numbering for object "
 399                     << name() << ".  Some of master "
 400                     << "points have been removed."
 401                     << abort(FatalError);
 402             }
 403         }
 404
 405         newRpm.insert(key, value);
 406     }
 407
 408     // Renumber the cut point edge pair map. Need to take a copy!
 409     const Map<Pair<edge> > cpepm = cutPointEdgePairMap();
 410
 411     Map<Pair<edge> >* newCpepmPtr = new Map<Pair<edge> >(cpepm.size());
 412     Map<Pair<edge> >& newCpepm = *newCpepmPtr;
 413
 414     const labelList cpepmToc = cpepm.toc();
 415
 416     forAll(cpepmToc, cpepmTocI)
 417     {
 418         key = reversePointMap[cpepmToc[cpepmTocI]];
 419
 420         const Pair<edge>& oldPe = cpepm.find(cpepmToc[cpepmTocI])();
 421
 422         // Re-do the edges in global addressing
 423         const label ms = reversePointMap[oldPe.first().start()];
 424         const label me = reversePointMap[oldPe.first().end()];
 425
 426         const label ss = reversePointMap[oldPe.second().start()];
 427         const label se = reversePointMap[oldPe.second().end()];
 428
 429         if (debug)
 430         {
 431             // Check if all the mapped cells are live
 432             if (key < 0 || ms < 0 || me < 0 || ss < 0 || se < 0)
 433             {
 434                 FatalErrorIn
 435                 (
 436                     "void slidingInterface::renumberAttachedAddressing("
 437                     "const mapPolyMesh& m) const"
 438                 )   << "Error in cut point edge pair map numbering for object "
 439                     << name() << ".  Some of master points have been removed."
 440                     << abort(FatalError);
 441             }
 442         }
 443
 444         newCpepm.insert(key, Pair<edge>(edge(ms, me), edge(ss, se)));
 445     }
 446
 447     if (!projectedSlavePointsPtr_)
 448     {
 449         FatalErrorIn
 450         (
 451             "void slidingInterface::renumberAttachedAddressing("
 452             "const mapPolyMesh& m) const"
 453         )   << "Error in projected point numbering for object " << name()
 454             << abort(FatalError);
 455     }
 456
 457     // Renumber the projected slave zone points
 458     const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
 459
 460     pointField* newProjectedSlavePointsPtr
 461     (
 462         new pointField(projectedSlavePoints.size())
 463     );
 464     pointField& newProjectedSlavePoints = *newProjectedSlavePointsPtr;
 465
 466     const labelList& sfzPointRenumber =
 467         m.faceZonePointMap()[slaveFaceZoneID_.index()];
 468
 469     forAll(newProjectedSlavePoints, pointI)
 470     {
 471         if (sfzPointRenumber[pointI] > -1)
 472         {
 473             newProjectedSlavePoints[pointI] =
 474                 projectedSlavePoints[sfzPointRenumber[pointI]];
 475         }
 476     }
 477
 478     // Re-set the lists
 479     clearAttachedAddressing();
 480
 481     deleteDemandDrivenData(projectedSlavePointsPtr_);
 482
 483     masterFaceCellsPtr_ = newMfcPtr;
 484     slaveFaceCellsPtr_ = newSfcPtr;
 485
 486     masterStickOutFacesPtr_ = newMsofPtr;
 487     slaveStickOutFacesPtr_ = newSsofPtr;
 488
 489     retiredPointMapPtr_ = newRpmPtr;
 490     cutPointEdgePairMapPtr_ = newCpepmPtr;
 491     projectedSlavePointsPtr_ = newProjectedSlavePointsPtr;
 492 }
 493
 494
 495 const Foam::labelList& Foam::slidingInterface::masterFaceCells() const
 496 {
 497     if (!masterFaceCellsPtr_)
 498     {
 499         FatalErrorIn
 500         (
 501             "const labelList& slidingInterface::masterFaceCells() const"
 502         )   << "Master zone face-cell addressing not available for object "
 503             << name()
 504             << abort(FatalError);
 505     }
 506
 507     return *masterFaceCellsPtr_;
 508 }
 509
 510
 511 const Foam::labelList& Foam::slidingInterface::slaveFaceCells() const
 512 {
 513     if (!slaveFaceCellsPtr_)
 514     {
 515         FatalErrorIn
 516         (
 517             "const labelList& slidingInterface::slaveFaceCells() const"
 518         )   << "Slave zone face-cell addressing not available for object "
 519             << name()
 520             << abort(FatalError);
 521     }
 522
 523     return *slaveFaceCellsPtr_;
 524 }
 525
 526
 527 const Foam::labelList& Foam::slidingInterface::masterStickOutFaces() const
 528 {
 529     if (!masterStickOutFacesPtr_)
 530     {
 531         FatalErrorIn
 532         (
 533             "const labelList& slidingInterface::masterStickOutFaces() const"
 534         )   << "Master zone stick-out face addressing not available for object "
 535             << name()
 536             << abort(FatalError);
 537     }
 538
 539     return *masterStickOutFacesPtr_;
 540 }
 541
 542
 543 const Foam::labelList& Foam::slidingInterface::slaveStickOutFaces() const
 544 {
 545     if (!slaveStickOutFacesPtr_)
 546     {
 547         FatalErrorIn
 548         (
 549             "const labelList& slidingInterface::slaveStickOutFaces() const"
 550         )   << "Slave zone stick-out face addressing not available for object "
 551             << name()
 552             << abort(FatalError);
 553     }
 554
 555     return *slaveStickOutFacesPtr_;
 556 }
 557
 558
 559 const Foam::Map<Foam::label>& Foam::slidingInterface::retiredPointMap() const
 560 {
 561     if (!retiredPointMapPtr_)
 562     {
 563         FatalErrorIn
 564         (
 565             "const Map<label>& slidingInterface::retiredPointMap() const"
 566         )   << "Retired point map not available for object " << name()
 567             << abort(FatalError);
 568     }
 569
 570     return *retiredPointMapPtr_;
 571 }
 572
 573
 574 const Foam::Map<Foam::Pair<Foam::edge> >&
 575 Foam::slidingInterface::cutPointEdgePairMap() const
 576 {
 577     if (!cutPointEdgePairMapPtr_)
 578     {
 579         FatalErrorIn
 580         (
 581             "const Map<Pair<edge> >& slidingInterface::"
 582             "cutPointEdgePairMap() const"
 583         )   << "Retired point map not available for object " << name()
 584             << abort(FatalError);
 585     }
 586
 587     return *cutPointEdgePairMapPtr_;
 588 }
 589
 590
 591 // ************************************************************************* //