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