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