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