src/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/removeCellLayer.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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     Remove a layer of cells and prepare addressing data
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "layerAdditionRemoval.H"
  31 #include "polyMesh.H"
  32 #include "primitiveMesh.H"
  33 #include "polyTopoChange.H"
  34 #include "oppositeFace.H"
  35 #include "polyTopoChanger.H"
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 bool Foam::layerAdditionRemoval::validCollapse() const
  40 {
  41     // Check for valid layer collapse
  42     // - no boundary-to-boundary collapse
  43
  44     if (debug)
  45     {
  46         Pout << "Checking layer collapse for object " << name() << endl;
  47     }
  48
  49     // Grab the face collapse mapping
  50     const polyMesh& mesh = topoChanger().mesh();
  51
  52     const labelList& ftc = facesPairing();
  53     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
  54
  55     label nBoundaryHits = 0;
  56
  57     forAll (mf, faceI)
  58     {
  59         if
  60         (
  61             !mesh.isInternalFace(mf[faceI])
  62          && !mesh.isInternalFace(ftc[faceI])
  63         )
  64         {
  65             nBoundaryHits++;
  66         }
  67     }
  68
  69
  70     if (debug)
  71     {
  72         Pout<< "Finished checking layer collapse for object "
  73             << name() <<".  Number of boundary-on-boundary hits: "
  74             << nBoundaryHits << endl;
  75     }
  76
  77     if (returnReduce(nBoundaryHits, sumOp<label>()) > 0)
  78     {
  79         return false;
  80     }
  81     else
  82     {
  83         return true;
  84     }
  85 }
  86
  87 void Foam::layerAdditionRemoval::removeCellLayer
  88 (
  89     polyTopoChange& ref
  90 ) const
  91 {
  92     // Algorithm for layer removal.  Second phase: topological change
  93     // 2)  Go through all the faces of the master cell layer and remove
  94     //     the ones that are not in the master face zone.
  95     //
  96     // 3)  Grab all the faces coming out of points that are collapsed
  97     //     and select the ones that are not marked for removal.  Go
  98     //     through the remaining faces and replace the point to remove by
  99     //     the equivalent point in the master face zone.
 100     if (debug)
 101     {
 102         Pout << "Removing the cell layer for object " << name() << endl;
 103     }
 104
 105     const polyMesh& mesh = topoChanger().mesh();
 106
 107     const labelList& ptc = pointsPairing();
 108     const labelList& ftc = facesPairing();
 109
 110     // Remove all the cells from the master layer
 111     const labelList& mc =
 112         topoChanger().mesh().faceZones()[faceZoneID_.index()].masterCells();
 113
 114     forAll (mc, cellI)
 115     {
 116         ref.setAction(polyRemoveCell(mc[cellI]));
 117     }
 118
 119     // Remove all the faces from the master layer cells which are not in
 120     // the master face layer
 121     labelHashSet facesToRemoveMap(mc.size()*primitiveMesh::facesPerCell_);
 122
 123     const cellList& cells = mesh.cells();
 124
 125     forAll (mc, cellI)
 126     {
 127         const cell& curCell = cells[mc[cellI]];
 128
 129         forAll (curCell, faceI)
 130         {
 131             // Check if the face is in the master zone.  If not, remove it
 132             if
 133             (
 134                 mesh.faceZones().whichZone(curCell[faceI])
 135              != faceZoneID_.index()
 136             )
 137             {
 138                 facesToRemoveMap.insert(curCell[faceI]);
 139             }
 140         }
 141     }
 142
 143     forAllConstIter(labelHashSet, facesToRemoveMap, iter)
 144     {
 145         ref.setAction(polyRemoveFace(iter.key()));
 146     }
 147
 148     // Remove all points that will be collapsed
 149     forAll (ptc, pointI)
 150     {
 151         ref.setAction(polyRemovePoint(ptc[pointI]));
 152     }
 153
 154     // Grab all faces coming off points to be deleted.  If the face
 155     // has not been deleted, replace the removed point with the
 156     // equivalent from the master layer.
 157
 158     // Make a map of all point to be removed, giving the master point label
 159     // for each of them
 160
 161     Map<label> removedPointMap(2*ptc.size());
 162
 163     const labelList& meshPoints =
 164         mesh.faceZones()[faceZoneID_.index()]().meshPoints();
 165
 166     forAll (ptc, pointI)
 167     {
 168         removedPointMap.insert(ptc[pointI], meshPoints[pointI]);
 169     }
 170
 171     const labelListList& pf = mesh.pointFaces();
 172
 173     const faceList& faces = mesh.faces();
 174     const labelList& own = mesh.faceOwner();
 175     const labelList& nei = mesh.faceNeighbour();
 176
 177     // Make a list of faces to be modified using the map to avoid duplicates
 178     labelHashSet facesToModify(ptc.size()*primitiveMesh::facesPerPoint_);
 179
 180     forAll (ptc, pointI)
 181     {
 182         const labelList& curFaces = pf[ptc[pointI]];
 183
 184         forAll (curFaces, faceI)
 185         {
 186             if (!facesToRemoveMap.found(curFaces[faceI]))
 187             {
 188                 facesToModify.insert(curFaces[faceI]);
 189             }
 190         }
 191     }
 192
 193     labelList ftm = facesToModify.toc();
 194
 195 //Pout << "faces to modify: " << ftm << endl;
 196
 197     forAll (ftm, faceI)
 198     {
 199         // For every face to modify, copy the face and re-map the vertices.
 200         // It is known all the faces will be changed since they hang off
 201         // re-mapped vertices
 202         label curFaceID = ftm[faceI];
 203
 204         face newFace(faces[curFaceID]);
 205
 206         forAll (newFace, pointI)
 207         {
 208             Map<label>::iterator rpmIter =
 209                 removedPointMap.find(newFace[pointI]);
 210
 211             if (rpmIter != removedPointMap.end())
 212             {
 213                 // Point mapped. Replace it
 214                 newFace[pointI] = rpmIter();
 215             }
 216         }
 217
 218 //Pout<< "face label: " << curFaceID
 219 //    << " old face: " << faces[curFaceID]
 220 //    << " new face: " << newFace << endl;
 221
 222         // Get face zone and its flip
 223         label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
 224         bool modifiedFaceZoneFlip = false;
 225
 226         if (modifiedFaceZone >= 0)
 227         {
 228             modifiedFaceZoneFlip =
 229                 mesh.faceZones()[modifiedFaceZone].flipMap()
 230                 [
 231                     mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
 232                 ];
 233         }
 234
 235         label newNeighbour = -1;
 236
 237         if (curFaceID < mesh.nInternalFaces())
 238         {
 239             newNeighbour = nei[curFaceID];
 240         }
 241         else
 242         {
 243             newNeighbour = -1;
 244         }
 245
 246         // Modify the face
 247         ref.setAction
 248         (
 249             polyModifyFace
 250             (
 251                 newFace,                // modified face
 252                 curFaceID,              // label of face being modified
 253                 own[curFaceID],         // owner
 254                 newNeighbour,           // neighbour
 255                 false,                  // face flip
 256                 mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
 257                 false,                  // remove from zone
 258                 modifiedFaceZone,       // zone for face
 259                 modifiedFaceZoneFlip    // face flip in zone
 260             )
 261         );
 262     }
 263
 264     // Modify the faces in the master layer to point past the removed cells
 265
 266     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
 267     const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
 268
 269     forAll (mf, faceI)
 270     {
 271         // Grab the owner and neighbour of the faces to be collapsed and get rid
 272         // of the cell to be removed
 273         label masterSideCell = own[mf[faceI]];
 274
 275         if (masterSideCell == mc[faceI])
 276         {
 277             // Owner cell of the face is being removed.
 278             // Grab the neighbour instead
 279             if (mesh.isInternalFace(mf[faceI]))
 280             {
 281                 masterSideCell = nei[mf[faceI]];
 282             }
 283             else
 284             {
 285                 masterSideCell = -1;
 286             }
 287         }
 288
 289         label slaveSideCell = own[ftc[faceI]];
 290
 291         if (slaveSideCell == mc[faceI])
 292         {
 293             // Owner cell of the face is being removed.
 294             // Grab the neighbour instead
 295             if (mesh.isInternalFace(ftc[faceI]))
 296             {
 297                 slaveSideCell = nei[ftc[faceI]];
 298             }
 299             else
 300             {
 301                 slaveSideCell = -1;
 302             }
 303         }
 304
 305         // Find out if the face needs to be flipped
 306         label newOwner = -1;
 307         label newNeighbour = -1;
 308         bool flipFace = false;
 309         label newPatchID = -1;
 310         label newZoneID = -1;
 311
 312         // Is any of the faces a boundary face?  If so, grab the patch
 313         // A boundary-to-boundary collapse is checked for in validCollapse()
 314         // and cannnot happen here.
 315
 316         if (!mesh.isInternalFace(mf[faceI]))
 317         {
 318             // Master is the boundary face: it gets a new owner but no flip
 319             newOwner = slaveSideCell;
 320             newNeighbour = -1;
 321             flipFace = false;
 322             newPatchID = mesh.boundaryMesh().whichPatch(mf[faceI]);
 323             newZoneID = mesh.faceZones().whichZone(mf[faceI]);
 324         }
 325         else if (!mesh.isInternalFace(ftc[faceI]))
 326         {
 327             // Slave is the boundary face: grab its patch
 328             newOwner = slaveSideCell;
 329             newNeighbour = -1;
 330
 331             // Find out if the face flip is necessary
 332             if (own[mf[faceI]] == slaveSideCell)
 333             {
 334                 flipFace = false;
 335             }
 336             else
 337             {
 338                 flipFace = true;
 339             }
 340
 341             newPatchID = mesh.boundaryMesh().whichPatch(ftc[faceI]);
 342
 343             // The zone of the master face is preserved
 344             newZoneID = mesh.faceZones().whichZone(mf[faceI]);
 345         }
 346         else
 347         {
 348             // Both faces are internal: flip is decided based on which of the
 349             // new cells around it has a lower label
 350             newOwner = min(masterSideCell, slaveSideCell);
 351             newNeighbour = max(masterSideCell, slaveSideCell);
 352
 353             if (newOwner == own[mf[faceI]] || newNeighbour == nei[mf[faceI]])
 354             {
 355                 flipFace = false;
 356             }
 357             else
 358             {
 359                 flipFace = true;
 360             }
 361
 362             newPatchID = -1;
 363
 364             // The zone of the master face is preserved
 365             newZoneID = mesh.faceZones().whichZone(mf[faceI]);
 366         }
 367
 368         // Modify the face and flip if necessary
 369         face newFace = faces[mf[faceI]];
 370         bool zoneFlip = mfFlip[faceI];
 371
 372         if (flipFace)
 373         {
 374             newFace = newFace.reverseFace();
 375             zoneFlip = !zoneFlip;
 376         }
 377
 378 // Pout<< "Modifying face " << mf[faceI]
 379 //     << " newFace: " << newFace << nl
 380 //     << " newOwner: " << newOwner
 381 //     << " newNeighbour: " << newNeighbour
 382 //     << " flipFace: " << flipFace
 383 //     << " newPatchID: " << newPatchID
 384 //     << " newZoneID: " << newZoneID << nl
 385 //     << " oldOwn: " << own[mf[faceI]]
 386 //     << " oldNei: " << nei[mf[faceI]] << endl;
 387
 388         ref.setAction
 389         (
 390             polyModifyFace
 391             (
 392                 newFace,       // modified face
 393                 mf[faceI],     // label of face being modified
 394                 newOwner,      // owner
 395                 newNeighbour,  // neighbour
 396                 flipFace,      // flip
 397                 newPatchID,    // patch for face
 398                 false,         // remove from zone
 399                 newZoneID,     // new zone
 400                 zoneFlip       // face flip in zone
 401             )
 402         );
 403     }
 404 }
 405
 406
 407 // ************************************************************************* //