undoableMeshCutter.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-2022 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 "undoableMeshCutter.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "DynamicList.H"
30 #include "meshCutter.H"
31 #include "cellCuts.H"
32 #include "splitCell.H"
33 #include "polyTopoChangeMap.H"
34 #include "unitConversion.H"
35 #include "meshTools.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 defineTypeNameAndDebug(undoableMeshCutter, 0);
42 }
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 // For debugging
48 void Foam::undoableMeshCutter::printCellRefTree
49 (
50  Ostream& os,
51  const word& indent,
52  const splitCell* splitCellPtr
53 ) const
54 {
55  if (splitCellPtr)
56  {
57  os << indent << splitCellPtr->cellLabel() << endl;
58 
59  word subIndent = indent + "--";
60 
61  printCellRefTree(os, subIndent, splitCellPtr->master());
62 
63  printCellRefTree(os, subIndent, splitCellPtr->slave());
64  }
65 }
66 
67 
68 // For debugging
69 void Foam::undoableMeshCutter::printRefTree(Ostream& os) const
70 {
71  forAllConstIter(Map<splitCell*>, liveSplitCells_, iter)
72  {
73  const splitCell* splitPtr = iter();
74 
75  // Walk to top (master path only)
76  while (splitPtr->parent())
77  {
78  if (!splitPtr->isMaster())
79  {
80  splitPtr = nullptr;
81 
82  break;
83  }
84  else
85  {
86  splitPtr = splitPtr->parent();
87  }
88  }
89 
90  // If we have reached top along master path start printing.
91  if (splitPtr)
92  {
93  // Print from top down
94  printCellRefTree(os, word(""), splitPtr);
95  }
96  }
97 }
98 
99 
100 // Update all (cell) labels on splitCell structure.
101 // Do in two passes to prevent allocation if nothing changed.
102 void Foam::undoableMeshCutter::updateLabels
103 (
104  const labelList& map,
105  Map<splitCell*>& liveSplitCells
106 )
107 {
108  // Pass1 : check if changed
109 
110  bool changed = false;
111 
112  forAllConstIter(Map<splitCell*>, liveSplitCells, iter)
113  {
114  const splitCell* splitPtr = iter();
115 
116  if (!splitPtr)
117  {
119  << "Problem: null pointer on liveSplitCells list"
120  << abort(FatalError);
121  }
122 
123  label celli = splitPtr->cellLabel();
124 
125  if (celli != map[celli])
126  {
127  changed = true;
128 
129  break;
130  }
131  }
132 
133 
134  // Pass2: relabel
135 
136  if (changed)
137  {
138  // Build new liveSplitCells
139  // since new labels (= keys in Map) might clash with existing ones.
140  Map<splitCell*> newLiveSplitCells(2*liveSplitCells.size());
141 
142  forAllIter(Map<splitCell*>, liveSplitCells, iter)
143  {
144  splitCell* splitPtr = iter();
145 
146  label celli = splitPtr->cellLabel();
147 
148  label newCelli = map[celli];
149 
150  if (debug && (celli != newCelli))
151  {
152  Pout<< "undoableMeshCutter::updateLabels :"
153  << " Updating live (split)cell from " << celli
154  << " to " << newCelli << endl;
155  }
156 
157  if (newCelli >= 0)
158  {
159  // Update splitCell. Can do inplace since only one celli will
160  // refer to this structure.
161  splitPtr->cellLabel() = newCelli;
162 
163  // Update liveSplitCells
164  newLiveSplitCells.insert(newCelli, splitPtr);
165  }
166  }
167  liveSplitCells = newLiveSplitCells;
168  }
169 }
170 
171 
172 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
173 
174 // Construct from components
176 (
177  const polyMesh& mesh,
178  const bool undoable
179 )
180 :
181  meshCutter(mesh),
182  undoable_(undoable),
183  liveSplitCells_(mesh.nCells()/100 + 100),
184  faceRemover_
185  (
186  mesh,
187  Foam::cos(degToRad(30.0))
188  )
189 {}
190 
191 
192 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
193 
195 {
196  // Clean split cell tree.
197 
198  forAllIter(Map<splitCell*>, liveSplitCells_, iter)
199  {
200  splitCell* splitPtr = iter();
201 
202  while (splitPtr)
203  {
204  splitCell* parentPtr = splitPtr->parent();
205 
206  // Sever ties with parent. Also of other side of refinement since
207  // we are handling rest of tree so other side will not have to.
208  if (parentPtr)
209  {
210  splitCell* otherSidePtr = splitPtr->getOther();
211 
212  otherSidePtr->parent() = nullptr;
213 
214  splitPtr->parent() = nullptr;
215  }
216 
217  // Delete splitCell (updates pointer on parent to itself)
218  delete splitPtr;
219 
220  splitPtr = parentPtr;
221  }
222  }
223 }
224 
225 
226 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
227 
229 (
230  const cellCuts& cuts,
231  polyTopoChange& meshMod
232 )
233 {
234  // Insert commands to actually cut cells
235  meshCutter::setRefinement(cuts, meshMod);
236 
237  if (undoable_)
238  {
239  // Use cells cut in this iteration to update splitCell tree.
240  forAllConstIter(Map<label>, addedCells(), iter)
241  {
242  label celli = iter.key();
243 
244  label addedCelli = iter();
245 
246 
247  // Newly created split cell. (celli -> celli + addedCelli)
248 
249  // Check if celli already part of split.
250  Map<splitCell*>::iterator findCell =
251  liveSplitCells_.find(celli);
252 
253  if (findCell == liveSplitCells_.end())
254  {
255  // Celli not yet split. It cannot be unlive split cell
256  // since that would be illegal to split in the first
257  // place.
258 
259  // Create 0th level. Null parent to denote this.
260  splitCell* parentPtr = new splitCell(celli, nullptr);
261 
262  splitCell* masterPtr = new splitCell(celli, parentPtr);
263 
264  splitCell* slavePtr = new splitCell(addedCelli, parentPtr);
265 
266  // Store newly created cells on parent together with face
267  // that splits them
268  parentPtr->master() = masterPtr;
269  parentPtr->slave() = slavePtr;
270 
271  // Insert master and slave into live splitcell list
272 
273  if (liveSplitCells_.found(addedCelli))
274  {
276  << "problem addedCell:" << addedCelli
277  << abort(FatalError);
278  }
279 
280  liveSplitCells_.insert(celli, masterPtr);
281  liveSplitCells_.insert(addedCelli, slavePtr);
282  }
283  else
284  {
285  // Cell that was split has been split again.
286  splitCell* parentPtr = findCell();
287 
288  // It is no longer live
289  liveSplitCells_.erase(findCell);
290 
291  splitCell* masterPtr = new splitCell(celli, parentPtr);
292 
293  splitCell* slavePtr = new splitCell(addedCelli, parentPtr);
294 
295  // Store newly created cells on parent together with face
296  // that splits them
297  parentPtr->master() = masterPtr;
298  parentPtr->slave() = slavePtr;
299 
300  // Insert master and slave into live splitcell list
301 
302  if (liveSplitCells_.found(addedCelli))
303  {
305  << "problem addedCell:" << addedCelli
306  << abort(FatalError);
307  }
308 
309  liveSplitCells_.insert(celli, masterPtr);
310  liveSplitCells_.insert(addedCelli, slavePtr);
311  }
312  }
313 
314  if (debug & 2)
315  {
316  Pout<< "** After refinement: liveSplitCells_:" << endl;
317 
318  printRefTree(Pout);
319  }
320  }
321 }
322 
323 
325 {
326  // Update mesh cutter for new labels.
328 
329  // No need to update cell walker for new labels since does not store any.
330 
331  // Update faceRemover for new labels
332  faceRemover_.topoChange(map);
333 
334  if (undoable_)
335  {
336  // Update all live split cells for mesh mapper.
337  updateLabels(map.reverseCellMap(), liveSplitCells_);
338  }
339 }
340 
341 
343 {
344  if (!undoable_)
345  {
347  << "Only call if constructed with unrefinement capability"
348  << abort(FatalError);
349  }
350 
351  DynamicList<label> liveSplitFaces(liveSplitCells_.size());
352 
353  forAllConstIter(Map<splitCell*>, liveSplitCells_, iter)
354  {
355  const splitCell* splitPtr = iter();
356 
357  if (!splitPtr->parent())
358  {
360  << "Live split cell without parent" << endl
361  << "splitCell:" << splitPtr->cellLabel()
362  << abort(FatalError);
363  }
364 
365  // Check if not top of refinement and whether it is the master side
366  if (splitPtr->isMaster())
367  {
368  splitCell* slavePtr = splitPtr->getOther();
369 
370  if
371  (
372  liveSplitCells_.found(slavePtr->cellLabel())
373  && splitPtr->isUnrefined()
374  && slavePtr->isUnrefined()
375  )
376  {
377  // Both master and slave are live and are not refined.
378  // Find common face.
379 
380  label celli = splitPtr->cellLabel();
381 
382  label slaveCelli = slavePtr->cellLabel();
383 
384  label commonFacei =
386  (
387  mesh(),
388  celli,
389  slaveCelli
390  );
391 
392  liveSplitFaces.append(commonFacei);
393  }
394  }
395  }
396 
397  return liveSplitFaces.shrink();
398 }
399 
400 
402 {
403  // (code copied from getSplitFaces)
404 
405  if (!undoable_)
406  {
408  << "Only call if constructed with unrefinement capability"
409  << abort(FatalError);
410  }
411 
412  Map<label> addedCells(liveSplitCells_.size());
413 
414  forAllConstIter(Map<splitCell*>, liveSplitCells_, iter)
415  {
416  const splitCell* splitPtr = iter();
417 
418  if (!splitPtr->parent())
419  {
421  << "Live split cell without parent" << endl
422  << "splitCell:" << splitPtr->cellLabel()
423  << abort(FatalError);
424  }
425 
426  // Check if not top of refinement and whether it is the master side
427  if (splitPtr->isMaster())
428  {
429  splitCell* slavePtr = splitPtr->getOther();
430 
431  if
432  (
433  liveSplitCells_.found(slavePtr->cellLabel())
434  && splitPtr->isUnrefined()
435  && slavePtr->isUnrefined()
436  )
437  {
438  // Both master and slave are live and are not refined.
439  addedCells.insert(splitPtr->cellLabel(), slavePtr->cellLabel());
440  }
441  }
442  }
443  return addedCells;
444 }
445 
446 
448 (
449  const labelList& splitFaces,
450  polyTopoChange& meshMod
451 )
452 {
453  if (!undoable_)
454  {
456  << "Only call if constructed with unrefinement capability"
457  << abort(FatalError);
458  }
459 
460  // Check with faceRemover what faces will get removed. Note that this can
461  // be more (but never less) than splitFaces provided.
462  labelList cellRegion;
463  labelList cellRegionMaster;
464  labelList facesToRemove;
465 
466  faceRemover().compatibleRemoves
467  (
468  splitFaces, // pierced faces
469  cellRegion, // per cell -1 or region it is merged into
470  cellRegionMaster, // per region the master cell
471  facesToRemove // new faces to be removed.
472  );
473 
474  if (facesToRemove.size() != splitFaces.size())
475  {
476  Pout<< "cellRegion:" << cellRegion << endl;
477  Pout<< "cellRegionMaster:" << cellRegionMaster << endl;
478 
480  << "Faces to remove:" << splitFaces << endl
481  << "to be removed:" << facesToRemove
482  << abort(FatalError);
483  }
484 
485 
486  // Every face removed will result in neighbour and owner being merged
487  // into owner.
488  forAll(facesToRemove, facesToRemoveI)
489  {
490  label facei = facesToRemove[facesToRemoveI];
491 
492  if (!mesh().isInternalFace(facei))
493  {
495  << "Trying to remove face that is not internal"
496  << abort(FatalError);
497  }
498 
499  label own = mesh().faceOwner()[facei];
500 
501  label nbr = mesh().faceNeighbour()[facei];
502 
503  Map<splitCell*>::iterator ownFind = liveSplitCells_.find(own);
504 
505  Map<splitCell*>::iterator nbrFind = liveSplitCells_.find(nbr);
506 
507  if
508  (
509  (ownFind == liveSplitCells_.end())
510  || (nbrFind == liveSplitCells_.end())
511  )
512  {
513  // Can happen because of removeFaces adding extra faces to
514  // original splitFaces
515  }
516  else
517  {
518  // Face is original splitFace.
519 
520  splitCell* ownPtr = ownFind();
521 
522  splitCell* nbrPtr = nbrFind();
523 
524  splitCell* parentPtr = ownPtr->parent();
525 
526  // Various checks on sanity.
527 
528  if (debug)
529  {
530  Pout<< "Updating for removed splitFace " << facei
531  << " own:" << own << " nbr:" << nbr
532  << " ownPtr:" << ownPtr->cellLabel()
533  << " nbrPtr:" << nbrPtr->cellLabel()
534  << endl;
535  }
536  if (!parentPtr)
537  {
539  << "No parent for owner " << ownPtr->cellLabel()
540  << abort(FatalError);
541  }
542 
543  if (!nbrPtr->parent())
544  {
546  << "No parent for neighbour " << nbrPtr->cellLabel()
547  << abort(FatalError);
548  }
549 
550  if (parentPtr != nbrPtr->parent())
551  {
553  << "Owner and neighbour liveSplitCell entries do not have"
554  << " same parent. facei:" << facei << " owner:" << own
555  << " ownparent:" << parentPtr->cellLabel()
556  << " neighbour:" << nbr
557  << " nbrparent:" << nbrPtr->parent()->cellLabel()
558  << abort(FatalError);
559  }
560 
561  if
562  (
563  !ownPtr->isUnrefined()
564  || !nbrPtr->isUnrefined()
565  || parentPtr->isUnrefined()
566  )
567  {
568  // Live owner and neighbour are refined themselves.
570  << "Owner and neighbour liveSplitCell entries are"
571  << " refined themselves or the parent is not refined"
572  << endl
573  << "owner unrefined:" << ownPtr->isUnrefined()
574  << " neighbour unrefined:" << nbrPtr->isUnrefined()
575  << " master unrefined:" << parentPtr->isUnrefined()
576  << abort(FatalError);
577  }
578 
579  // Delete from liveSplitCell
580  liveSplitCells_.erase(ownFind);
581 
583  liveSplitCells_.erase(liveSplitCells_.find(nbr));
584 
585  // Delete entries themselves
586  delete ownPtr;
587  delete nbrPtr;
588 
589  //
590  // Update parent:
591  // - has parent itself: is part of split cell. Update cellLabel
592  // with merged cell one.
593  // - has no parent: is start of tree. Completely remove.
594 
595  if (parentPtr->parent())
596  {
597  // Update parent cell label to be new merged cell label
598  // (will be owner)
599  parentPtr->cellLabel() = own;
600 
601  // And insert into live cells (is ok since old entry with
602  // own as key has been removed above)
603  liveSplitCells_.insert(own, parentPtr);
604  }
605  else
606  {
607  // No parent so is start of tree. No need to keep splitCell
608  // tree.
609  delete parentPtr;
610  }
611  }
612  }
613 
614  // Insert all commands to combine cells. Never fails so don't have to
615  // test for success.
616  faceRemover().setRefinement
617  (
618  facesToRemove,
619  cellRegion,
620  cellRegionMaster,
621  meshMod
622  );
623 
624  return facesToRemove;
625 }
626 
627 
628 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
bool isMaster() const
Check if this is master cell of split.
Definition: splitCell.C:70
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
splitCell * getOther() const
Returns other half of split cell. I.e. slave if this is master.
Definition: splitCell.C:106
Unit conversion functions.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label cellLabel() const
Definition: splitCell.H:87
fvMesh & mesh
Description of cuts across cells.
Definition: cellCuts.H:108
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
undoableMeshCutter(const polyMesh &mesh, const bool undoable=true)
Construct from mesh and flag whether refinement pattern needs.
splitCell * slave() const
Definition: splitCell.H:117
An STL-conforming iterator.
Definition: HashTable.H:426
label getSharedFace(const primitiveMesh &, const label cell0, const label cell1)
Return face shared by two cells. Throws error if none found.
Definition: meshTools.C:419
dimensionedScalar cos(const dimensionedScalar &ds)
bool isUnrefined() const
Check if this is unrefined (i.e. has no master or slave)
Definition: splitCell.C:100
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
splitCell * parent() const
Definition: splitCell.H:97
Description of cell after splitting. Contains cellLabel and pointers to cells it it split in...
Definition: splitCell.H:51
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Cuts (splits) cells.
Definition: meshCutter.H:134
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
void setRefinement(const cellCuts &cuts, polyTopoChange &)
Refine cells acc. to cellCuts. Plays topology changes.
Map< label > getAddedCells() const
Like getSplitFaces but returns map from original to added cell.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
Definition: meshCutter.C:525
defineTypeNameAndDebug(combustionModel, 0)
const labelList & reverseCellMap() const
Reverse cell map.
splitCell * master() const
Definition: splitCell.H:107
labelList removeSplitFaces(const labelList &splitFaces, polyTopoChange &)
Remove some refinement. Needs to be supplied subset of.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Direct mesh changes based on v1.3 polyTopoChange syntax.
void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
Definition: meshCutter.C:998
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
labelList getSplitFaces() const
Calculate split faces from current liveCells. Only.
void topoChange(const polyTopoChangeMap &map)
Update stored refinement pattern for changes to mesh. Only.
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49