removeCells.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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  << "Size of exposedFaceLabels " << exposedFaceLabels.size()
192  << " differs from size of exposedPatchIDs "
193  << exposedPatchIDs.size()
194  << abort(FatalError);
195  }
196 
197  // List of new patchIDs
198  labelList newPatchID(mesh_.nFaces(), -1);
199 
200  forAll(exposedFaceLabels, i)
201  {
202  label patchi = exposedPatchIDs[i];
203 
204  if (patchi < 0 || patchi >= patches.size())
205  {
207  << "Invalid patch " << patchi
208  << " for exposed face " << exposedFaceLabels[i] << endl
209  << "Valid patches 0.." << patches.size()-1
210  << abort(FatalError);
211  }
212 
213  if (patches[patchi].coupled())
214  {
216  << "Trying to put exposed face " << exposedFaceLabels[i]
217  << " into a coupled patch : " << patches[patchi].name()
218  << endl
219  << "This is illegal."
220  << abort(FatalError);
221  }
222 
223  newPatchID[exposedFaceLabels[i]] = patchi;
224  }
225 
226 
227  // Create list of cells to be removed
228  boolList removedCell(mesh_.nCells(), false);
229 
230  // Go from labelList of cells-to-remove to a boolList and remove all
231  // cells mentioned.
232  forAll(cellLabels, i)
233  {
234  label celli = cellLabels[i];
235 
236  removedCell[celli] = true;
237 
238  // Pout<< "Removing cell " << celli
239  // << " cc:" << mesh_.cellCentres()[celli] << endl;
240 
241  meshMod.setAction(polyRemoveCell(celli));
242  }
243 
244 
245  // Remove faces that are no longer used. Modify faces that
246  // are used by one cell only.
247 
248  const faceList& faces = mesh_.faces();
249  const labelList& faceOwner = mesh_.faceOwner();
250  const labelList& faceNeighbour = mesh_.faceNeighbour();
251  const meshFaceZones& faceZones = mesh_.faceZones();
252 
253  // Count starting number of faces using each point. Keep up to date whenever
254  // removing a face.
255  labelList nFacesUsingPoint(mesh_.nPoints(), 0);
256 
257  forAll(faces, facei)
258  {
259  const face& f = faces[facei];
260 
261  forAll(f, fp)
262  {
263  nFacesUsingPoint[f[fp]]++;
264  }
265  }
266 
267 
268  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
269  {
270  const face& f = faces[facei];
271  label own = faceOwner[facei];
272  label nei = faceNeighbour[facei];
273 
274  if (removedCell[own])
275  {
276  if (removedCell[nei])
277  {
278  // Face no longer used
279  // Pout<< "Removing internal face " << facei
280  // << " fc:" << mesh_.faceCentres()[facei] << endl;
281 
282  meshMod.setAction(polyRemoveFace(facei));
283  uncount(f, nFacesUsingPoint);
284  }
285  else
286  {
287  if (newPatchID[facei] == -1)
288  {
290  << "No patchID provided for exposed face " << facei
291  << " on cell " << nei << nl
292  << "Did you provide patch IDs for all exposed faces?"
293  << abort(FatalError);
294  }
295 
296  // nei is remaining cell. Facei becomes external cell
297 
298  label zoneID = faceZones.whichZone(facei);
299  bool zoneFlip = false;
300 
301  if (zoneID >= 0)
302  {
303  const faceZone& fZone = faceZones[zoneID];
304  // Note: we reverse the owner/neighbour of the face
305  // so should also select the other side of the zone
306  zoneFlip = !fZone.flipMap()[fZone.whichFace(facei)];
307  }
308 
309  // Pout<< "Putting exposed internal face " << facei
310  // << " fc:" << mesh_.faceCentres()[facei]
311  // << " into patch " << newPatchID[facei] << endl;
312 
313  meshMod.setAction
314  (
316  (
317  f.reverseFace(), // modified face
318  facei, // label of face being modified
319  nei, // owner
320  -1, // neighbour
321  true, // face flip
322  newPatchID[facei], // patch for face
323  false, // remove from zone
324  zoneID, // zone for face
325  zoneFlip // face flip in zone
326  )
327  );
328  }
329  }
330  else if (removedCell[nei])
331  {
332  if (newPatchID[facei] == -1)
333  {
335  << "No patchID provided for exposed face " << facei
336  << " on cell " << own << nl
337  << "Did you provide patch IDs for all exposed faces?"
338  << abort(FatalError);
339  }
340 
341  // Pout<< "Putting exposed internal face " << facei
342  // << " fc:" << mesh_.faceCentres()[facei]
343  // << " into patch " << newPatchID[facei] << endl;
344 
345  // own is remaining cell. Facei becomes external cell.
346  label zoneID = faceZones.whichZone(facei);
347  bool zoneFlip = false;
348 
349  if (zoneID >= 0)
350  {
351  const faceZone& fZone = faceZones[zoneID];
352  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
353  }
354 
355  meshMod.setAction
356  (
358  (
359  f, // modified face
360  facei, // label of face being modified
361  own, // owner
362  -1, // neighbour
363  false, // face flip
364  newPatchID[facei], // patch for face
365  false, // remove from zone
366  zoneID, // zone for face
367  zoneFlip // face flip in zone
368  )
369  );
370  }
371  }
372 
373  forAll(patches, patchi)
374  {
375  const polyPatch& pp = patches[patchi];
376 
377  if (pp.coupled())
378  {
379  label facei = pp.start();
380 
381  forAll(pp, i)
382  {
383  if (newPatchID[facei] != -1)
384  {
385  // Pout<< "Putting uncoupled coupled face " << facei
386  // << " fc:" << mesh_.faceCentres()[facei]
387  // << " into patch " << newPatchID[facei] << endl;
388 
389  label zoneID = faceZones.whichZone(facei);
390  bool zoneFlip = false;
391 
392  if (zoneID >= 0)
393  {
394  const faceZone& fZone = faceZones[zoneID];
395  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
396  }
397 
398  meshMod.setAction
399  (
401  (
402  faces[facei], // modified face
403  facei, // label of face
404  faceOwner[facei], // owner
405  -1, // neighbour
406  false, // face flip
407  newPatchID[facei], // patch for face
408  false, // remove from zone
409  zoneID, // zone for face
410  zoneFlip // face flip in zone
411  )
412  );
413  }
414  else if (removedCell[faceOwner[facei]])
415  {
416  // Face no longer used
417  // Pout<< "Removing boundary face " << facei
418  // << " fc:" << mesh_.faceCentres()[facei]
419  // << endl;
420 
421  meshMod.setAction(polyRemoveFace(facei));
422  uncount(faces[facei], nFacesUsingPoint);
423  }
424 
425  facei++;
426  }
427  }
428  else
429  {
430  label facei = pp.start();
431 
432  forAll(pp, i)
433  {
434  if (newPatchID[facei] != -1)
435  {
437  << "new patchID provided for boundary face " << facei
438  << " even though it is not on a coupled face."
439  << abort(FatalError);
440  }
441 
442  if (removedCell[faceOwner[facei]])
443  {
444  // Face no longer used
445  // Pout<< "Removing boundary face " << facei
446  // << " fc:" << mesh_.faceCentres()[facei]
447  // << endl;
448 
449  meshMod.setAction(polyRemoveFace(facei));
450  uncount(faces[facei], nFacesUsingPoint);
451  }
452 
453  facei++;
454  }
455  }
456  }
457 
458 
459  // Remove points that are no longer used.
460  // Loop rewritten to not use pointFaces.
461 
462  forAll(nFacesUsingPoint, pointi)
463  {
464  if (nFacesUsingPoint[pointi] == 0)
465  {
466  // Pout<< "Removing unused point " << pointi
467  // << " at:" << mesh_.points()[pointi] << endl;
468 
469  meshMod.setAction(polyRemovePoint(pointi));
470  }
471  else if (nFacesUsingPoint[pointi] == 1)
472  {
474  << "point " << pointi << " at coordinate "
475  << mesh_.points()[pointi]
476  << " is only used by 1 face after removing cells."
477  << " This probably results in an illegal mesh."
478  << endl;
479  }
480  }
481 }
482 
483 
484 // ************************************************************************* //
const fvPatchList & patches
Class containing data for face removal.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const word & name() const
Return name.
Definition: IOobject.H:315
Class describing modification of a face.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void setRefinement(const labelList &cellsToRemove, const labelList &facesToExpose, const labelList &patchIDs, polyTopoChange &) const
Play commands into polyTopoChange to remove cells.
Definition: removeCells.C:179
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:307
face reverseFace() const
Return face with reverse direction.
Definition: face.C:256
Class containing data for cell removal.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:316
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
removeCells(const polyMesh &mesh, const bool syncPar=true)
Construct from mesh. syncPar: do parallel synchronisation.
Definition: removeCells.C:60
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronise values on all mesh faces.
Definition: syncTools.H:387
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: MeshZones.C:221
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
Definition: removeCells.C:73
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Class containing data for point removal.
Namespace for OpenFOAM.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.