src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C

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