1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
26 Remove a layer of cells and prepare addressing data
28 \*---------------------------------------------------------------------------*/
30 #include "layerAdditionRemoval.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"
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 bool Foam::layerAdditionRemoval::validCollapse() const
45 // Check for valid layer collapse
46 // - no boundary-to-boundary collapse
50 Pout << "Checking layer collapse for object " << name() << endl;
53 // Grab the face collapse mapping
54 const polyMesh& mesh = topoChanger().mesh();
56 const labelList& ftc = facesPairing();
57 const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
59 label nBoundaryHits = 0;
65 !mesh.isInternalFace(mf[faceI])
66 && !mesh.isInternalFace(ftc[faceI])
76 Pout<< "Finished checking layer collapse for object "
77 << name() <<". Number of boundary-on-boundary hits: "
78 << nBoundaryHits << endl;
81 if (returnReduce(nBoundaryHits, sumOp<label>()) > 0)
91 void Foam::layerAdditionRemoval::removeCellLayer
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.
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.
106 Pout << "Removing the cell layer for object " << name() << endl;
109 const polyMesh& mesh = topoChanger().mesh();
111 const labelList& ptc = pointsPairing();
112 const labelList& ftc = facesPairing();
114 // Remove all the cells from the master layer
115 const labelList& mc =
116 topoChanger().mesh().faceZones()[faceZoneID_.index()].masterCells();
120 ref.setAction(polyRemoveCell(mc[cellI]));
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_);
127 const cellList& cells = mesh.cells();
131 const cell& curCell = cells[mc[cellI]];
133 forAll (curCell, faceI)
135 // Check if the face is in the master zone. If not, remove it
138 mesh.faceZones().whichZone(curCell[faceI])
139 != faceZoneID_.index()
142 facesToRemoveMap.insert(curCell[faceI]);
147 forAllConstIter(labelHashSet, facesToRemoveMap, iter)
149 ref.setAction(polyRemoveFace(iter.key()));
152 // Remove all points that will be collapsed
155 ref.setAction(polyRemovePoint(ptc[pointI]));
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.
162 // Make a map of all point to be removed, giving the master point label
165 Map<label> removedPointMap(2*ptc.size());
167 const labelList& meshPoints =
168 mesh.faceZones()[faceZoneID_.index()]().meshPoints();
172 removedPointMap.insert(ptc[pointI], meshPoints[pointI]);
175 const labelListList& pf = mesh.pointFaces();
177 const faceList& faces = mesh.faces();
178 const labelList& own = mesh.faceOwner();
179 const labelList& nei = mesh.faceNeighbour();
181 // Make a list of faces to be modified using the map to avoid duplicates
182 labelHashSet facesToModify(ptc.size()*primitiveMesh::facesPerPoint_);
186 const labelList& curFaces = pf[ptc[pointI]];
188 forAll (curFaces, faceI)
190 if (!facesToRemoveMap.found(curFaces[faceI]))
192 facesToModify.insert(curFaces[faceI]);
197 labelList ftm = facesToModify.toc();
199 //Pout << "faces to modify: " << ftm << endl;
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];
208 face newFace(faces[curFaceID]);
210 forAll (newFace, pointI)
212 Map<label>::iterator rpmIter =
213 removedPointMap.find(newFace[pointI]);
215 if (rpmIter != removedPointMap.end())
217 // Point mapped. Replace it
218 newFace[pointI] = rpmIter();
222 //Pout<< "face label: " << curFaceID
223 // << " old face: " << faces[curFaceID]
224 // << " new face: " << newFace << endl;
226 // Get face zone and its flip
227 label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
228 bool modifiedFaceZoneFlip = false;
230 if (modifiedFaceZone >= 0)
232 modifiedFaceZoneFlip =
233 mesh.faceZones()[modifiedFaceZone].flipMap()
235 mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
241 if (curFaceID < mesh.nInternalFaces())
243 newNei = nei[curFaceID];
255 newFace, // modified face
256 curFaceID, // label of face being modified
257 own[curFaceID], // owner
260 mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
261 false, // remove from zone
262 modifiedFaceZone, // zone for face
263 modifiedFaceZoneFlip // face flip in zone
268 // Modify the faces in the master layer to point past the removed cells
270 const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
271 const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
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]];
279 if (masterSideCell == mc[faceI])
281 // Owner cell of the face is being removed.
282 // Grab the neighbour instead
283 masterSideCell = nei[mf[faceI]];
286 label slaveSideCell = own[ftc[faceI]];
288 if (slaveSideCell == mc[faceI])
290 // Owner cell of the face is being removed.
291 // Grab the neighbour instead
292 slaveSideCell = nei[ftc[faceI]];
295 // Find out if the face needs to be flipped
297 label newNeighbour = -1;
298 bool flipFace = false;
299 label newPatchID = -1;
300 label newZoneID = -1;
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.
306 if (!mesh.isInternalFace(mf[faceI]))
308 // Master is the boundary face: it gets a new owner but no flip
309 newOwner = slaveSideCell;
312 newPatchID = mesh.boundaryMesh().whichPatch(mf[faceI]);
313 newZoneID = mesh.faceZones().whichZone(mf[faceI]);
315 else if (!mesh.isInternalFace(ftc[faceI]))
317 // Slave is the boundary face: grab its patch
318 newOwner = slaveSideCell;
321 // Find out if the face flip is necessary
322 if (own[mf[faceI]] == slaveSideCell)
331 newPatchID = mesh.boundaryMesh().whichPatch(ftc[faceI]);
333 // The zone of the master face is preserved
334 newZoneID = mesh.faceZones().whichZone(mf[faceI]);
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);
343 if (newOwner == own[mf[faceI]] || newNeighbour == nei[mf[faceI]])
354 // The zone of the master face is preserved
355 newZoneID = mesh.faceZones().whichZone(mf[faceI]);
358 // Modify the face and flip if necessary
359 face newFace = faces[mf[faceI]];
360 bool zoneFlip = mfFlip[faceI];
364 newFace = newFace.reverseFace();
365 zoneFlip = !zoneFlip;
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;
382 newFace, // modified face
383 mf[faceI], // label of face being modified
385 newNeighbour, // neighbour
387 newPatchID, // patch for face
388 false, // remove from zone
389 newZoneID, // new zone
390 zoneFlip // face flip in zone
397 // ************************************************************************* //