removeCellLayer.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 \*---------------------------------------------------------------------------*/
25 
26 #include "layerAdditionRemoval.H"
27 #include "polyMesh.H"
28 #include "primitiveMesh.H"
29 #include "polyTopoChange.H"
30 #include "oppositeFace.H"
31 #include "polyTopoChanger.H"
32 #include "polyRemoveCell.H"
33 #include "polyRemoveFace.H"
34 #include "polyRemovePoint.H"
35 #include "polyModifyFace.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 
88 void Foam::layerAdditionRemoval::removeCellLayer
89 (
90  polyTopoChange& ref
91 ) const
92 {
93  // Algorithm for layer removal. Second phase: topological change
94  // 2) Go through all the faces of the master cell layer and remove
95  // the ones that are not in the master face zone.
96  //
97  // 3) Grab all the faces coming out of points that are collapsed
98  // and select the ones that are not marked for removal. Go
99  // through the remaining faces and replace the point to remove by
100  // the equivalent point in the master face zone.
101  if (debug)
102  {
103  Pout<< "Removing the cell layer for object " << name() << endl;
104  }
105 
106  const polyMesh& mesh = topoChanger().mesh();
107 
108  const labelList& own = mesh.faceOwner();
109  const labelList& nei = mesh.faceNeighbour();
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  mesh.faceZones()[faceZoneID_.index()].masterCells();
117 
118  forAll(mc, facei)
119  {
120  label slaveSideCell = own[ftc[facei]];
121 
122  if (mesh.isInternalFace(ftc[facei]) && slaveSideCell == mc[facei])
123  {
124  // Owner cell of the face is being removed.
125  // Grab the neighbour instead
126  slaveSideCell = nei[ftc[facei]];
127  }
128 
129  ref.setAction(polyRemoveCell(mc[facei], slaveSideCell));
130  }
131 
132  // Remove all the faces from the master layer cells which are not in
133  // the master face layer
134  labelHashSet facesToRemoveMap(mc.size()*primitiveMesh::facesPerCell_);
135 
136  const cellList& cells = mesh.cells();
137 
138  forAll(mc, celli)
139  {
140  const cell& curCell = cells[mc[celli]];
141 
142  forAll(curCell, facei)
143  {
144  // Check if the face is in the master zone. If not, remove it
145  if
146  (
147  mesh.faceZones().whichZone(curCell[facei])
148  != faceZoneID_.index()
149  )
150  {
151  facesToRemoveMap.insert(curCell[facei]);
152  }
153  }
154  }
155 
156  forAllConstIter(labelHashSet, facesToRemoveMap, iter)
157  {
158  ref.setAction(polyRemoveFace(iter.key()));
159  }
160 
161  // Remove all points that will be collapsed
162  forAll(ptc, pointi)
163  {
164  ref.setAction(polyRemovePoint(ptc[pointi]));
165  }
166 
167  // Grab all faces coming off points to be deleted. If the face
168  // has not been deleted, replace the removed point with the
169  // equivalent from the master layer.
170 
171  // Make a map of all point to be removed, giving the master point label
172  // for each of them
173 
174  Map<label> removedPointMap(2*ptc.size());
175 
176  const labelList& meshPoints =
177  mesh.faceZones()[faceZoneID_.index()]().meshPoints();
178 
179  forAll(ptc, pointi)
180  {
181  removedPointMap.insert(ptc[pointi], meshPoints[pointi]);
182  }
183 
184  const labelListList& pf = mesh.pointFaces();
185 
186  const faceList& faces = mesh.faces();
187 
188  // Make a list of faces to be modified using the map to avoid duplicates
189  labelHashSet facesToModify(ptc.size()*primitiveMesh::facesPerPoint_);
190 
191  forAll(ptc, pointi)
192  {
193  const labelList& curFaces = pf[ptc[pointi]];
194 
195  forAll(curFaces, facei)
196  {
197  if (!facesToRemoveMap.found(curFaces[facei]))
198  {
199  facesToModify.insert(curFaces[facei]);
200  }
201  }
202  }
203 
204  labelList ftm = facesToModify.toc();
205 
206  if (debug > 1)
207  {
208  Pout<< "faces to modify: " << ftm << endl;
209  }
210 
211  forAll(ftm, facei)
212  {
213  // For every face to modify, copy the face and re-map the vertices.
214  // It is known all the faces will be changed since they hang off
215  // re-mapped vertices
216  label curFaceID = ftm[facei];
217 
218  face newFace(faces[curFaceID]);
219 
220  forAll(newFace, pointi)
221  {
222  Map<label>::iterator rpmIter =
223  removedPointMap.find(newFace[pointi]);
224 
225  if (rpmIter != removedPointMap.end())
226  {
227  // Point mapped. Replace it
228  newFace[pointi] = rpmIter();
229  }
230  }
231 
232  if (debug > 1)
233  {
234  Pout<< "face label: " << curFaceID
235  << " old face: " << faces[curFaceID]
236  << " new face: " << newFace << endl;
237  }
238 
239  // Get face zone and its flip
240  label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
241  bool modifiedFaceZoneFlip = false;
242 
243  if (modifiedFaceZone >= 0)
244  {
245  const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
246  modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
247  }
248 
249  label newNeighbour = -1;
250 
251  if (curFaceID < mesh.nInternalFaces())
252  {
253  newNeighbour = nei[curFaceID];
254  }
255 
256  // Modify the face
257  ref.setAction
258  (
259  polyModifyFace
260  (
261  newFace, // modified face
262  curFaceID, // label of face being modified
263  own[curFaceID], // owner
264  newNeighbour, // neighbour
265  false, // face flip
266  mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
267  false, // remove from zone
268  modifiedFaceZone, // zone for face
269  modifiedFaceZoneFlip // face flip in zone
270  )
271  );
272  }
273 
274  // Modify the faces in the master layer to point past the removed cells
275 
276  const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
277  const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
278 
279  forAll(mf, facei)
280  {
281  // Grab the owner and neighbour of the faces to be collapsed and get rid
282  // of the cell to be removed
283  label masterSideCell = own[mf[facei]];
284 
285  if (masterSideCell == mc[facei])
286  {
287  if (mesh.isInternalFace(mf[facei]))
288  {
289  // Owner cell of the face is being removed.
290  // Grab the neighbour instead
291  masterSideCell = nei[mf[facei]];
292  }
293  else
294  {
295  masterSideCell = -1;
296  }
297  }
298 
299  label slaveSideCell = own[ftc[facei]];
300 
301  if (slaveSideCell == mc[facei])
302  {
303  if (mesh.isInternalFace(ftc[facei]))
304  {
305  // Owner cell of the face is being removed.
306  // Grab the neighbour instead
307  slaveSideCell = nei[ftc[facei]];
308  }
309  else
310  {
311  slaveSideCell = -1;
312  }
313  }
314 
315  // Find out if the face needs to be flipped
316  label newOwner = -1;
317  label newNeighbour = -1;
318  bool flipFace = false;
319  label newPatchID = -1;
320  label newZoneID = -1;
321 
322  // Is any of the faces a boundary face? If so, grab the patch
323  // A boundary-to-boundary collapse is checked for in validCollapse()
324  // and cannot happen here.
325 
326  if (!mesh.isInternalFace(mf[facei]))
327  {
328  // Master is the boundary face: it gets a new owner but no flip
329  newOwner = slaveSideCell;
330  newNeighbour = -1;
331  flipFace = false;
332  newPatchID = mesh.boundaryMesh().whichPatch(mf[facei]);
333  newZoneID = mesh.faceZones().whichZone(mf[facei]);
334  }
335  else if (!mesh.isInternalFace(ftc[facei]))
336  {
337  // Slave is the boundary face: grab its patch
338  newOwner = slaveSideCell;
339  newNeighbour = -1;
340 
341  // Find out if the face flip is necessary
342  if (own[mf[facei]] == slaveSideCell)
343  {
344  flipFace = false;
345  }
346  else
347  {
348  flipFace = true;
349  }
350 
351  newPatchID = mesh.boundaryMesh().whichPatch(ftc[facei]);
352 
353  // The zone of the master face is preserved
354  newZoneID = mesh.faceZones().whichZone(mf[facei]);
355  }
356  else
357  {
358  // Both faces are internal: flip is decided based on which of the
359  // new cells around it has a lower label
360  newOwner = min(masterSideCell, slaveSideCell);
361  newNeighbour = max(masterSideCell, slaveSideCell);
362 
363  if (newOwner == own[mf[facei]] || newNeighbour == nei[mf[facei]])
364  {
365  flipFace = false;
366  }
367  else
368  {
369  flipFace = true;
370  }
371 
372  newPatchID = -1;
373 
374  // The zone of the master face is preserved
375  newZoneID = mesh.faceZones().whichZone(mf[facei]);
376  }
377 
378  // Modify the face and flip if necessary
379  face newFace = faces[mf[facei]];
380  bool zoneFlip = mfFlip[facei];
381 
382  if (flipFace)
383  {
384  newFace.flip();
385  zoneFlip = !zoneFlip;
386  }
387 
388  if (debug > 1)
389  {
390  Pout<< "Modifying face " << mf[facei]
391  << " newFace: " << newFace << nl
392  << " newOwner: " << newOwner
393  << " newNeighbour: " << newNeighbour
394  << " flipFace: " << flipFace
395  << " newPatchID: " << newPatchID
396  << " newZoneID: " << newZoneID << nl
397  << " oldOwn: " << own[mf[facei]]
398  << " oldNei: " << nei[mf[facei]] << endl;
399  }
400 
401  ref.setAction
402  (
403  polyModifyFace
404  (
405  newFace, // modified face
406  mf[facei], // label of face being modified
407  newOwner, // owner
408  newNeighbour, // neighbour
409  flipFace, // flip
410  newPatchID, // patch for face
411  false, // remove from zone
412  newZoneID, // new zone
413  zoneFlip // face flip in zone
414  )
415  );
416  }
417 }
418 
419 
420 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
HashTable< label, label, Hash< label > >::iterator iterator
Definition: Map.H:56
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static const unsigned facesPerPoint_
Estimated number of faces per point.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
static const unsigned facesPerCell_
Estimated number of faces per cell.
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
static const char nl
Definition: Ostream.H:262
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const word & name() const
Return name of this modifier.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const polyMesh & mesh() const
Return the mesh reference.
List< cell > cellList
list of cells
Definition: cellList.H:42