removeCells.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 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 "removeCells.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "polyRemoveCell.H"
30 #include "polyRemoveFace.H"
31 #include "polyModifyFace.H"
32 #include "polyRemovePoint.H"
33 #include "syncTools.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(removeCells, 0);
40 }
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 void Foam::removeCells::uncount
45 (
46  const labelList& f,
47  labelList& nUsage
48 )
49 {
50  forAll(f, fp)
51  {
52  nUsage[f[fp]]--;
53  }
54 }
55 
56 
57 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
58 
60 (
61  const polyMesh& mesh,
62  const bool syncPar
63 )
64 :
65  mesh_(mesh),
66  syncPar_(syncPar)
67 {}
68 
69 
70 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71 
73 (
74  const labelList& cellLabels
75 ) const
76 {
77  // Create list of cells to be removed
78  boolList removedCell(mesh_.nCells(), false);
79 
80  // Go from labelList of cells-to-remove to a boolList.
81  forAll(cellLabels, i)
82  {
83  removedCell[cellLabels[i]] = true;
84  }
85 
86 
87  const labelList& faceOwner = mesh_.faceOwner();
88  const labelList& faceNeighbour = mesh_.faceNeighbour();
89 
90  // Count cells using face.
91  labelList nCellsUsingFace(mesh_.nFaces(), 0);
92 
93  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
94  {
95  label own = faceOwner[faceI];
96  label nei = faceNeighbour[faceI];
97 
98  if (!removedCell[own])
99  {
100  nCellsUsingFace[faceI]++;
101  }
102  if (!removedCell[nei])
103  {
104  nCellsUsingFace[faceI]++;
105  }
106  }
107 
108  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
109  {
110  label own = faceOwner[faceI];
111 
112  if (!removedCell[own])
113  {
114  nCellsUsingFace[faceI]++;
115  }
116  }
117 
118  // Coupled faces: add number of cells using face across couple.
119  if (syncPar_)
120  {
122  (
123  mesh_,
124  nCellsUsingFace,
126  );
127  }
128 
129  // Now nCellsUsingFace:
130  // 0 : internal face whose both cells get deleted
131  // boundary face whose all cells get deleted
132  // 1 : internal face that gets exposed
133  // unaffected (uncoupled) boundary face
134  // coupled boundary face that gets exposed ('uncoupled')
135  // 2 : unaffected internal face
136  // unaffected coupled boundary face
137 
138  DynamicList<label> exposedFaces(mesh_.nFaces()/10);
139 
140  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
141  {
142  if (nCellsUsingFace[faceI] == 1)
143  {
144  exposedFaces.append(faceI);
145  }
146  }
147 
148  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
149 
150  forAll(patches, patchI)
151  {
152  const polyPatch& pp = patches[patchI];
153 
154  if (pp.coupled())
155  {
156  label faceI = pp.start();
157 
158  forAll(pp, i)
159  {
160  label own = faceOwner[faceI];
161 
162  if (nCellsUsingFace[faceI] == 1 && !removedCell[own])
163  {
164  // My owner not removed but other side is so has to become
165  // normal, uncoupled, boundary face
166  exposedFaces.append(faceI);
167  }
168 
169  faceI++;
170  }
171  }
172  }
173 
174  return exposedFaces.shrink();
175 }
176 
177 
179 (
180  const labelList& cellLabels,
181  const labelList& exposedFaceLabels,
182  const labelList& exposedPatchIDs,
183  polyTopoChange& meshMod
184 ) const
185 {
186  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
187 
188  if (exposedFaceLabels.size() != exposedPatchIDs.size())
189  {
191  (
192  "removeCells::setRefinement(const labelList&"
193  ", const labelList&, const labelList&, polyTopoChange&)"
194  ) << "Size of exposedFaceLabels " << exposedFaceLabels.size()
195  << " differs from size of exposedPatchIDs "
196  << exposedPatchIDs.size()
197  << abort(FatalError);
198  }
199 
200  // List of new patchIDs
201  labelList newPatchID(mesh_.nFaces(), -1);
202 
203  forAll(exposedFaceLabels, i)
204  {
205  label patchI = exposedPatchIDs[i];
206 
207  if (patchI < 0 || patchI >= patches.size())
208  {
210  (
211  "removeCells::setRefinement(const labelList&"
212  ", const labelList&, const labelList&, polyTopoChange&)"
213  ) << "Invalid patch " << patchI
214  << " for exposed face " << exposedFaceLabels[i] << endl
215  << "Valid patches 0.." << patches.size()-1
216  << abort(FatalError);
217  }
218 
219  if (patches[patchI].coupled())
220  {
222  (
223  "removeCells::setRefinement(const labelList&"
224  ", const labelList&, const labelList&, polyTopoChange&)"
225  ) << "Trying to put exposed face " << exposedFaceLabels[i]
226  << " into a coupled patch : " << patches[patchI].name()
227  << endl
228  << "This is illegal."
229  << abort(FatalError);
230  }
231 
232  newPatchID[exposedFaceLabels[i]] = patchI;
233  }
234 
235 
236  // Create list of cells to be removed
237  boolList removedCell(mesh_.nCells(), false);
238 
239  // Go from labelList of cells-to-remove to a boolList and remove all
240  // cells mentioned.
241  forAll(cellLabels, i)
242  {
243  label cellI = cellLabels[i];
244 
245  removedCell[cellI] = true;
246 
247  //Pout<< "Removing cell " << cellI
248  // << " cc:" << mesh_.cellCentres()[cellI] << endl;
249 
250  meshMod.setAction(polyRemoveCell(cellI));
251  }
252 
253 
254  // Remove faces that are no longer used. Modify faces that
255  // are used by one cell only.
256 
257  const faceList& faces = mesh_.faces();
258  const labelList& faceOwner = mesh_.faceOwner();
259  const labelList& faceNeighbour = mesh_.faceNeighbour();
260  const faceZoneMesh& faceZones = mesh_.faceZones();
261 
262  // Count starting number of faces using each point. Keep up to date whenever
263  // removing a face.
264  labelList nFacesUsingPoint(mesh_.nPoints(), 0);
265 
266  forAll(faces, faceI)
267  {
268  const face& f = faces[faceI];
269 
270  forAll(f, fp)
271  {
272  nFacesUsingPoint[f[fp]]++;
273  }
274  }
275 
276 
277  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
278  {
279  const face& f = faces[faceI];
280  label own = faceOwner[faceI];
281  label nei = faceNeighbour[faceI];
282 
283  if (removedCell[own])
284  {
285  if (removedCell[nei])
286  {
287  // Face no longer used
288  //Pout<< "Removing internal face " << faceI
289  // << " fc:" << mesh_.faceCentres()[faceI] << endl;
290 
291  meshMod.setAction(polyRemoveFace(faceI));
292  uncount(f, nFacesUsingPoint);
293  }
294  else
295  {
296  if (newPatchID[faceI] == -1)
297  {
299  (
300  "removeCells::setRefinement(const labelList&"
301  ", const labelList&, const labelList&"
302  ", polyTopoChange&)"
303  ) << "No patchID provided for exposed face " << faceI
304  << " on cell " << nei << nl
305  << "Did you provide patch IDs for all exposed faces?"
306  << abort(FatalError);
307  }
308 
309  // nei is remaining cell. FaceI becomes external cell
310 
311  label zoneID = faceZones.whichZone(faceI);
312  bool zoneFlip = false;
313 
314  if (zoneID >= 0)
315  {
316  const faceZone& fZone = faceZones[zoneID];
317  // Note: we reverse the owner/neighbour of the face
318  // so should also select the other side of the zone
319  zoneFlip = !fZone.flipMap()[fZone.whichFace(faceI)];
320  }
321 
322  //Pout<< "Putting exposed internal face " << faceI
323  // << " fc:" << mesh_.faceCentres()[faceI]
324  // << " into patch " << newPatchID[faceI] << endl;
325 
326  meshMod.setAction
327  (
329  (
330  f.reverseFace(), // modified face
331  faceI, // label of face being modified
332  nei, // owner
333  -1, // neighbour
334  true, // face flip
335  newPatchID[faceI], // patch for face
336  false, // remove from zone
337  zoneID, // zone for face
338  zoneFlip // face flip in zone
339  )
340  );
341  }
342  }
343  else if (removedCell[nei])
344  {
345  if (newPatchID[faceI] == -1)
346  {
348  (
349  "removeCells::setRefinement(const labelList&"
350  ", const labelList&, const labelList&"
351  ", polyTopoChange&)"
352  ) << "No patchID provided for exposed face " << faceI
353  << " on cell " << own << nl
354  << "Did you provide patch IDs for all exposed faces?"
355  << abort(FatalError);
356  }
357 
358  //Pout<< "Putting exposed internal face " << faceI
359  // << " fc:" << mesh_.faceCentres()[faceI]
360  // << " into patch " << newPatchID[faceI] << endl;
361 
362  // own is remaining cell. FaceI becomes external cell.
363  label zoneID = faceZones.whichZone(faceI);
364  bool zoneFlip = false;
365 
366  if (zoneID >= 0)
367  {
368  const faceZone& fZone = faceZones[zoneID];
369  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
370  }
371 
372  meshMod.setAction
373  (
375  (
376  f, // modified face
377  faceI, // label of face being modified
378  own, // owner
379  -1, // neighbour
380  false, // face flip
381  newPatchID[faceI], // patch for face
382  false, // remove from zone
383  zoneID, // zone for face
384  zoneFlip // face flip in zone
385  )
386  );
387  }
388  }
389 
390  forAll(patches, patchI)
391  {
392  const polyPatch& pp = patches[patchI];
393 
394  if (pp.coupled())
395  {
396  label faceI = pp.start();
397 
398  forAll(pp, i)
399  {
400  if (newPatchID[faceI] != -1)
401  {
402  //Pout<< "Putting uncoupled coupled face " << faceI
403  // << " fc:" << mesh_.faceCentres()[faceI]
404  // << " into patch " << newPatchID[faceI] << endl;
405 
406  label zoneID = faceZones.whichZone(faceI);
407  bool zoneFlip = false;
408 
409  if (zoneID >= 0)
410  {
411  const faceZone& fZone = faceZones[zoneID];
412  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
413  }
414 
415  meshMod.setAction
416  (
418  (
419  faces[faceI], // modified face
420  faceI, // label of face
421  faceOwner[faceI], // owner
422  -1, // neighbour
423  false, // face flip
424  newPatchID[faceI], // patch for face
425  false, // remove from zone
426  zoneID, // zone for face
427  zoneFlip // face flip in zone
428  )
429  );
430  }
431  else if (removedCell[faceOwner[faceI]])
432  {
433  // Face no longer used
434  //Pout<< "Removing boundary face " << faceI
435  // << " fc:" << mesh_.faceCentres()[faceI]
436  // << endl;
437 
438  meshMod.setAction(polyRemoveFace(faceI));
439  uncount(faces[faceI], nFacesUsingPoint);
440  }
441 
442  faceI++;
443  }
444  }
445  else
446  {
447  label faceI = pp.start();
448 
449  forAll(pp, i)
450  {
451  if (newPatchID[faceI] != -1)
452  {
454  (
455  "removeCells::setRefinement(const labelList&"
456  ", const labelList&, const labelList&"
457  ", polyTopoChange&)"
458  ) << "new patchID provided for boundary face " << faceI
459  << " even though it is not on a coupled face."
460  << abort(FatalError);
461  }
462 
463  if (removedCell[faceOwner[faceI]])
464  {
465  // Face no longer used
466  //Pout<< "Removing boundary face " << faceI
467  // << " fc:" << mesh_.faceCentres()[faceI]
468  // << endl;
469 
470  meshMod.setAction(polyRemoveFace(faceI));
471  uncount(faces[faceI], nFacesUsingPoint);
472  }
473 
474  faceI++;
475  }
476  }
477  }
478 
479 
480  // Remove points that are no longer used.
481  // Loop rewritten to not use pointFaces.
482 
483  forAll(nFacesUsingPoint, pointI)
484  {
485  if (nFacesUsingPoint[pointI] == 0)
486  {
487  //Pout<< "Removing unused point " << pointI
488  // << " at:" << mesh_.points()[pointI] << endl;
489 
490  meshMod.setAction(polyRemovePoint(pointI));
491  }
492  else if (nFacesUsingPoint[pointI] == 1)
493  {
494  WarningIn
495  (
496  "removeCells::setRefinement(const labelList&"
497  ", const labelList&, const labelList&"
498  ", polyTopoChange&)"
499  ) << "point " << pointI << " at coordinate "
500  << mesh_.points()[pointI]
501  << " is only used by 1 face after removing cells."
502  << " This probably results in an illegal mesh."
503  << endl;
504  }
505  }
506 }
507 
508 
509 // ************************************************************************* //
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
Definition: removeCells.C:73
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
label size() const
Return the number of elements in the PtrList.
Definition: PtrListI.H:32
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:232
face reverseFace() const
Return face with reverse direction.
Definition: face.C:616
labelList f(nPoints)
Class containing data for point removal.
Class containing data for cell removal.
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
void setRefinement(const labelList &cellsToRemove, const labelList &facesToExpose, const labelList &patchIDs, polyTopoChange &) const
Play commands into polyTopoChange to remove cells.
Definition: removeCells.C:179
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:380
Namespace for OpenFOAM.
Class containing data for face removal.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:317
patches[0]
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
Foam::polyBoundaryMesh.
#define forAll(list, i)
Definition: UList.H:421
Class describing modification of a face.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & name() const
Return name.
Definition: IOobject.H:260
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
removeCells(const polyMesh &mesh, const bool syncPar=true)
Construct from mesh. syncPar: do parallel synchronization.
Definition: removeCells.C:60
List< label > labelList
A List of labels.
Definition: labelList.H:56
Direct mesh changes based on v1.3 polyTopoChange syntax.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
defineTypeNameAndDebug(combustionModel, 0)