edgeCollapser.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 "edgeCollapser.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "globalMeshData.H"
30 #include "syncTools.H"
31 #include "PointEdgeWave.H"
32 #include "globalIndex.H"
33 #include "removePoints.H"
34 #include "motionSmoother.H"
35 #include "OFstream.H"
36 
37 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 defineTypeNameAndDebug(edgeCollapser, 0);
42 }
43 
44 
46 (
47  const polyMesh& mesh,
48  const dictionary& meshQualityDict
49 )
50 {
51  labelHashSet badFaces(mesh.nFaces()/100);
52  DynamicList<label> checkFaces(mesh.nFaces());
53 
54  const vectorField& fAreas = mesh.faceAreas();
55 
56  scalar faceAreaLimit = SMALL;
57 
58  forAll(fAreas, fI)
59  {
60  if (mag(fAreas[fI]) > faceAreaLimit)
61  {
62  checkFaces.append(fI);
63  }
64  }
65 
66  Info<< endl;
67 
69  (
70  false,
71  mesh,
72  meshQualityDict,
73  checkFaces,
74  badFaces
75  );
76 
77  return badFaces;
78 }
79 
80 
82 (
83  const polyMesh& mesh,
84  const dictionary& meshQualityDict,
85  PackedBoolList& isErrorPoint
86 )
87 {
89  (
90  mesh,
91  meshQualityDict
92  );
93 
94  label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
95 
96  forAllConstIter(labelHashSet, badFaces, iter)
97  {
98  const face& f = mesh.faces()[iter.key()];
99 
100  forAll(f, pI)
101  {
102  isErrorPoint[f[pI]] = true;
103  }
104  }
105 
107  (
108  mesh,
109  isErrorPoint,
111  0
112  );
113 
114  return nBadFaces;
115 }
116 
117 
118 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
119 
120 Foam::labelList Foam::edgeCollapser::edgesFromPoints
121 (
122  const label& facei,
123  const labelList& pointLabels
124 ) const
125 {
126  labelList edgeLabels(pointLabels.size() - 1, -1);
127 
128  const labelList& faceEdges = mesh_.faceEdges()[facei];
129  const edgeList& edges = mesh_.edges();
130 
131  label count = 0;
132 
133  forAll(faceEdges, eI)
134  {
135  const label edgeI = faceEdges[eI];
136  const edge& e = edges[edgeI];
137 
138  label match = 0;
139 
140  forAll(pointLabels, pI)
141  {
142  if (e[0] == pointLabels[pI])
143  {
144  match++;
145  }
146 
147  if (e[1] == pointLabels[pI])
148  {
149  match++;
150  }
151  }
152 
153  if (match == 2)
154  {
155  // Edge found
156  edgeLabels[count++] = edgeI;
157  }
158  }
159 
160  if (count != edgeLabels.size())
161  {
162  edgeLabels.setSize(count);
163  }
164 
165  return edgeLabels;
166 }
167 
168 
169 void Foam::edgeCollapser::collapseToEdge
170 (
171  const label facei,
172  const pointField& pts,
173  const labelList& pointPriority,
174  const vector& collapseAxis,
175  const point& fC,
176  const labelList& facePtsNeg,
177  const labelList& facePtsPos,
178  const scalarList& dNeg,
179  const scalarList& dPos,
180  const scalar dShift,
181  PackedBoolList& collapseEdge,
182  Map<point>& collapsePointToLocation
183 ) const
184 {
185  // Negative half
186 
187  Foam::point collapseToPtA(GREAT, GREAT, GREAT);
188  //collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
189 
190  label maxPriority = labelMin;
191  DynamicList<label> maxPriorityPts(max(dNeg.size(), dPos.size()));
192 
193  forAll(facePtsNeg, fPtI)
194  {
195  const label facePointi = facePtsNeg[fPtI];
196  const label facePtPriority = pointPriority[facePointi];
197 
198  if (facePtPriority > maxPriority)
199  {
200  maxPriority = facePtPriority;
201  maxPriorityPts.clear();
202  maxPriorityPts.append(facePointi);
203  }
204  else if (facePtPriority == maxPriority)
205  {
206  maxPriorityPts.append(facePointi);
207  }
208  }
209 
210  if (!maxPriorityPts.empty())
211  {
212  Foam::point averagePt(Zero);
213 
214  forAll(maxPriorityPts, ptI)
215  {
216  averagePt += pts[maxPriorityPts[ptI]];
217  }
218 
219  collapseToPtA = averagePt/maxPriorityPts.size();
220 // collapseToPtA = pts[maxPriorityPts.first()];
221  }
222 
223  maxPriority = labelMin;
224  maxPriorityPts.clear();
225 
226  labelList faceEdgesNeg = edgesFromPoints(facei, facePtsNeg);
227 
228  forAll(faceEdgesNeg, edgeI)
229  {
230  collapseEdge[faceEdgesNeg[edgeI]] = true;
231  }
232 
233  forAll(facePtsNeg, pI)
234  {
235  collapsePointToLocation.set(facePtsNeg[pI], collapseToPtA);
236  }
237 
238 
239  // Positive half
240  Foam::point collapseToPtB(GREAT, GREAT, GREAT);
241 // = collapseAxis*(sum(dPos)/dPos.size() - dShift) + fC;
242 
243  forAll(facePtsPos, fPtI)
244  {
245  const label facePointi = facePtsPos[fPtI];
246  const label facePtPriority = pointPriority[facePointi];
247 
248  if (facePtPriority > maxPriority)
249  {
250  maxPriority = facePtPriority;
251  maxPriorityPts.clear();
252  maxPriorityPts.append(facePointi);
253  }
254  else if (facePtPriority == maxPriority)
255  {
256  maxPriorityPts.append(facePointi);
257  }
258  }
259 
260  if (!maxPriorityPts.empty())
261  {
262  Foam::point averagePt(Zero);
263 
264  forAll(maxPriorityPts, ptI)
265  {
266  averagePt += pts[maxPriorityPts[ptI]];
267  }
268 
269  collapseToPtB = averagePt/maxPriorityPts.size();
270 // collapseToPtB = pts[maxPriorityPts.first()];
271  }
272 
273  labelList faceEdgesPos = edgesFromPoints(facei, facePtsPos);
274 
275  forAll(faceEdgesPos, edgeI)
276  {
277  collapseEdge[faceEdgesPos[edgeI]] = true;
278  }
279 
280  forAll(facePtsPos, pI)
281  {
282  collapsePointToLocation.set(facePtsPos[pI], collapseToPtB);
283  }
284 }
285 
286 
287 void Foam::edgeCollapser::collapseToPoint
288 (
289  const label& facei,
290  const pointField& pts,
291  const labelList& pointPriority,
292  const point& fC,
293  const labelList& facePts,
294  PackedBoolList& collapseEdge,
295  Map<point>& collapsePointToLocation
296 ) const
297 {
298  const face& f = mesh_.faces()[facei];
299 
300  Foam::point collapseToPt = fC;
301 
302  label maxPriority = labelMin;
303  DynamicList<label> maxPriorityPts(f.size());
304 
305  forAll(facePts, fPtI)
306  {
307  const label facePointi = facePts[fPtI];
308  const label facePtPriority = pointPriority[facePointi];
309 
310  if (facePtPriority > maxPriority)
311  {
312  maxPriority = facePtPriority;
313  maxPriorityPts.clear();
314  maxPriorityPts.append(facePointi);
315  }
316  else if (facePtPriority == maxPriority)
317  {
318  maxPriorityPts.append(facePointi);
319  }
320  }
321 
322  if (!maxPriorityPts.empty())
323  {
324  Foam::point averagePt(Zero);
325 
326  forAll(maxPriorityPts, ptI)
327  {
328  averagePt += pts[maxPriorityPts[ptI]];
329  }
330 
331  collapseToPt = averagePt/maxPriorityPts.size();
332 
333 // collapseToPt = pts[maxPriorityPts.first()];
334  }
335 
336 // DynamicList<label> faceBoundaryPts(f.size());
337 // DynamicList<label> faceFeaturePts(f.size());
338 //
339 // forAll(facePts, fPtI)
340 // {
341 // if (pointPriority[facePts[fPtI]] == 1)
342 // {
343 // faceFeaturePts.append(facePts[fPtI]);
344 // }
345 // else if (pointPriority[facePts[fPtI]] == 0)
346 // {
347 // faceBoundaryPts.append(facePts[fPtI]);
348 // }
349 // }
350 //
351 // if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
352 // {
353 // if (!faceFeaturePts.empty())
354 // {
355 // collapseToPt = pts[faceFeaturePts.first()];
356 // }
357 // else if (faceBoundaryPts.size() == 2)
358 // {
359 // collapseToPt =
360 // 0.5
361 // *(
362 // pts[faceBoundaryPts[0]]
363 // + pts[faceBoundaryPts[1]]
364 // );
365 // }
366 // else if (faceBoundaryPts.size() <= f.size())
367 // {
368 // face bFace(faceBoundaryPts);
369 //
370 // collapseToPt = bFace.centre(pts);
371 // }
372 // }
373 
374  const labelList& faceEdges = mesh_.faceEdges()[facei];
375 
376  forAll(faceEdges, eI)
377  {
378  const label edgeI = faceEdges[eI];
379  collapseEdge[edgeI] = true;
380  }
381 
382  forAll(f, pI)
383  {
384  collapsePointToLocation.set(f[pI], collapseToPt);
385  }
386 }
387 
388 
389 void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
390 (
391  const face& f,
392  const point& fC,
393  vector& collapseAxis,
394  scalar& aspectRatio
395 ) const
396 {
397  const pointField& pts = mesh_.points();
398 
399  tensor J = f.inertia(pts, fC);
400 
401  // Find the dominant collapse direction by finding the eigenvector
402  // that corresponds to the normal direction, discarding it. The
403  // eigenvector corresponding to the smaller of the two remaining
404  // eigenvalues is the dominant axis in a high aspect ratio face.
405 
406  scalar magJ = mag(J);
407 
408  scalar detJ = SMALL;
409 
410  if (magJ > VSMALL)
411  {
412  // Normalise inertia tensor to remove problems with small values
413 
414  J /= mag(J);
415  // J /= cmptMax(J);
416  // J /= max(eigenValues(J).x(), SMALL);
417 
418  // Calculating determinant, including stabilisation for zero or
419  // small negative values
420 
421  detJ = max(det(J), SMALL);
422  }
423 
424  if (detJ < 1e-5)
425  {
426  collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts);
427 
428  // It is possible that all the points of a face are the same
429  if (magSqr(collapseAxis) > VSMALL)
430  {
431  collapseAxis /= mag(collapseAxis);
432  }
433 
434  // Empirical correlation for high aspect ratio faces
435 
436  aspectRatio = Foam::sqrt(0.35/detJ);
437  }
438  else
439  {
440  vector eVals = eigenValues(J);
441 
442  if (mag(eVals.y() - eVals.x()) < 100*SMALL)
443  {
444  // First two eigenvalues are the same: i.e. a square face
445 
446  // Cannot necessarily determine linearly independent
447  // eigenvectors, or any at all, use longest edge direction.
448 
449  collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts);
450 
451  collapseAxis /= mag(collapseAxis);
452 
453  aspectRatio = 1.0;
454  }
455  else
456  {
457  // The maximum eigenvalue (z()) must be the direction of the
458  // normal, as it has the greatest value. The minimum eigenvalue
459  // is the dominant collapse axis for high aspect ratio faces.
460 
461  collapseAxis = eigenVector(J, eVals.x());
462 
463  // The inertia calculation describes the mass distribution as a
464  // function of distance squared to the axis, so the square root of
465  // the ratio of face-plane moments gives a good indication of the
466  // aspect ratio.
467 
468  aspectRatio = Foam::sqrt(eVals.y()/max(eVals.x(), SMALL));
469  }
470  }
471 }
472 
473 
474 Foam::scalarField Foam::edgeCollapser::calcTargetFaceSizes() const
475 {
476  scalarField targetFaceSizes(mesh_.nFaces(), -1);
477 
478  const scalarField& V = mesh_.cellVolumes();
479  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
480 
481  const labelList& cellOwner = mesh_.faceOwner();
482  const labelList& cellNeighbour = mesh_.faceNeighbour();
483 
484  const label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
485 
486  // Calculate face size from cell volumes for internal faces
487  for (label intFacei = 0; intFacei < mesh_.nInternalFaces(); ++intFacei)
488  {
489  const scalar cellOwnerVol = max(0.0, V[cellOwner[intFacei]]);
490  const scalar cellNeighbourVol = max(0.0, V[cellNeighbour[intFacei]]);
491 
492  scalar targetFaceSizeA = Foam::pow(cellOwnerVol, 1.0/3.0);
493  scalar targetFaceSizeB = Foam::pow(cellNeighbourVol, 1.0/3.0);
494 
495  targetFaceSizes[intFacei] = 0.5*(targetFaceSizeA + targetFaceSizeB);
496  }
497 
498  scalarField neiCellVolumes(nBoundaryFaces, -1);
499 
500  // Now do boundary faces
501  forAll(patches, patchi)
502  {
503  const polyPatch& patch = patches[patchi];
504 
505  label bFacei = patch.start() - mesh_.nInternalFaces();
506 
507  if (patch.coupled())
508  {
509  // Processor boundary face: Need to get the cell volume on the other
510  // processor
511  const labelUList& faceCells = patch.faceCells();
512 
513  forAll(faceCells, facei)
514  {
515  neiCellVolumes[bFacei++] = max(0.0, V[faceCells[facei]]);
516  }
517  }
518  else
519  {
520  // Normal boundary face: Just use owner cell volume to calculate
521  // the target face size
522  forAll(patch, patchFacei)
523  {
524  const label extFacei = patchFacei + patch.start();
525  const scalar cellOwnerVol = max(0.0, V[cellOwner[extFacei]]);
526 
527  targetFaceSizes[extFacei] = Foam::pow(cellOwnerVol, 1.0/3.0);
528  }
529  }
530  }
531 
532  syncTools::swapBoundaryFaceList(mesh_, neiCellVolumes);
533 
534  forAll(patches, patchi)
535  {
536  const polyPatch& patch = patches[patchi];
537 
538  label bFacei = patch.start() - mesh_.nInternalFaces();
539 
540  if (patch.coupled())
541  {
542  forAll(patch, patchFacei)
543  {
544  const label localFacei = patchFacei + patch.start();
545  const scalar cellOwnerVol = max(0.0, V[cellOwner[localFacei]]);
546  const scalar cellNeighbourVol = neiCellVolumes[bFacei++];
547 
548  scalar targetFaceSizeA = Foam::pow(cellOwnerVol, 1.0/3.0);
549  scalar targetFaceSizeB = Foam::pow(cellNeighbourVol, 1.0/3.0);
550 
551  targetFaceSizes[localFacei]
552  = 0.5*(targetFaceSizeA + targetFaceSizeB);
553  }
554  }
555  }
556 
557  // Returns a characteristic length, not an area
558  return targetFaceSizes;
559 }
560 
561 
562 Foam::edgeCollapser::collapseType Foam::edgeCollapser::collapseFace
563 (
564  const labelList& pointPriority,
565  const face& f,
566  const label facei,
567  const scalar targetFaceSize,
568  PackedBoolList& collapseEdge,
569  Map<point>& collapsePointToLocation,
570  const scalarField& faceFilterFactor
571 ) const
572 {
573  const scalar collapseSizeLimitCoeff = faceFilterFactor[facei];
574 
575  const pointField& pts = mesh_.points();
576 
577  labelList facePts(f);
578 
579  const Foam::point fC = f.centre(pts);
580 
581  const scalar fA = f.mag(pts);
582 
583  vector collapseAxis = Zero;
584  scalar aspectRatio = 1.0;
585 
586  faceCollapseAxisAndAspectRatio(f, fC, collapseAxis, aspectRatio);
587 
588  // The signed distance along the collapse axis passing through the
589  // face centre that each vertex projects to.
590 
591  scalarField d(f.size());
592 
593  forAll(f, fPtI)
594  {
595  const Foam::point& pt = pts[f[fPtI]];
596 
597  d[fPtI] = (collapseAxis & (pt - fC));
598  }
599 
600  // Sort the projected distances and the corresponding vertex
601  // indices along the collapse axis
602 
603  labelList oldToNew;
604 
605  sortedOrder(d, oldToNew);
606 
607  oldToNew = invert(oldToNew.size(), oldToNew);
608 
609  inplaceReorder(oldToNew, d);
610 
611  inplaceReorder(oldToNew, facePts);
612 
613  // Shift the points so that they are relative to the centre of the
614  // collapse line.
615 
616  scalar dShift = -0.5*(d.first() + d.last());
617 
618  d += dShift;
619 
620  // Form two lists, one for each half of the set of points
621  // projected along the collapse axis.
622 
623  // Middle value, index of first entry in the second half
624  label middle = -1;
625 
626  forAll(d, dI)
627  {
628  if (d[dI] > 0)
629  {
630  middle = dI;
631 
632  break;
633  }
634  }
635 
636  if (middle == -1)
637  {
638 // SeriousErrorInFunction
639 // << "middle == -1, " << f << " " << d
640 // << endl;//abort(FatalError);
641 
642  return noCollapse;
643  }
644 
645  // Negative half
646  SubList<scalar> dNeg(d, middle, 0);
647  SubList<label> facePtsNeg(facePts, middle, 0);
648 
649  // Positive half
650  SubList<scalar> dPos(d, d.size() - middle, middle);
651  SubList<label> facePtsPos(facePts, d.size() - middle, middle);
652 
653  // Defining how close to the midpoint (M) of the projected
654  // vertices line a projected vertex (X) can be before making this
655  // an invalid edge collapse
656  //
657  // X---X-g----------------M----X-----------g----X--X
658  //
659  // Only allow a collapse if all projected vertices are outwith
660  // guardFraction (g) of the distance form the face centre to the
661  // furthest vertex in the considered direction
662 
663  if (dNeg.size() == 0 || dPos.size() == 0)
664  {
666  << "All points on one side of face centre, not collapsing."
667  << endl;
668  }
669 
670 // Info<< "Face : " << f << nl
671 // << " Collapse Axis: " << collapseAxis << nl
672 // << " Aspect Ratio : " << aspectRatio << endl;
673 
674  collapseType typeOfCollapse = noCollapse;
675 
676  if (magSqr(collapseAxis) < VSMALL)
677  {
678  typeOfCollapse = toPoint;
679  }
680  else if (fA < aspectRatio*sqr(targetFaceSize*collapseSizeLimitCoeff))
681  {
682  if
683  (
684  allowEarlyCollapseToPoint_
685  && (d.last() - d.first())
686  < targetFaceSize
687  *allowEarlyCollapseCoeff_*maxCollapseFaceToPointSideLengthCoeff_
688  )
689  {
690  typeOfCollapse = toPoint;
691  }
692  else if
693  (
694  (dNeg.last() < guardFraction_*dNeg.first())
695  && (dPos.first() > guardFraction_*dPos.last())
696  )
697  {
698  typeOfCollapse = toEdge;
699  }
700  else if
701  (
702  (d.last() - d.first())
703  < targetFaceSize
704  *maxCollapseFaceToPointSideLengthCoeff_
705  )
706  {
707  // If the face can't be collapsed to an edge, and it has a
708  // small enough span, collapse it to a point.
709  typeOfCollapse = toPoint;
710  }
711  }
712 
713  if (typeOfCollapse == toPoint)
714  {
715  collapseToPoint
716  (
717  facei,
718  pts,
719  pointPriority,
720  fC,
721  facePts,
722  collapseEdge,
723  collapsePointToLocation
724  );
725  }
726  else if (typeOfCollapse == toEdge)
727  {
728  collapseToEdge
729  (
730  facei,
731  pts,
732  pointPriority,
733  collapseAxis,
734  fC,
735  facePtsNeg,
736  facePtsPos,
737  dNeg,
738  dPos,
739  dShift,
740  collapseEdge,
741  collapsePointToLocation
742  );
743  }
744 
745  return typeOfCollapse;
746 }
747 
748 
749 Foam::label Foam::edgeCollapser::edgeMaster
750 (
751  const labelList& pointPriority,
752  const edge& e
753 ) const
754 {
755  label masterPoint = -1;
756 
757  const label e0 = e.start();
758  const label e1 = e.end();
759 
760  const label e0Priority = pointPriority[e0];
761  const label e1Priority = pointPriority[e1];
762 
763  if (e0Priority > e1Priority)
764  {
765  masterPoint = e0;
766  }
767  else if (e0Priority < e1Priority)
768  {
769  masterPoint = e1;
770  }
771  else if (e0Priority == e1Priority)
772  {
773  masterPoint = e0;
774  }
775 
776 // // Collapse edge to point with higher priority.
777 // if (pointPriority[e0] >= 0)
778 // {
779 // if (pointPriority[e1] >= 0)
780 // {
781 // // Both points have high priority. Choose one to collapse to.
782 // // Note: should look at feature edges/points!
783 // masterPoint = e0;
784 // }
785 // else
786 // {
787 // masterPoint = e0;
788 // }
789 // }
790 // else
791 // {
792 // if (pointPriority[e1] >= 0)
793 // {
794 // masterPoint = e1;
795 // }
796 // else
797 // {
798 // // None on boundary. Neither is a master.
799 // return -1;
800 // }
801 // }
802 
803  return masterPoint;
804 }
805 
806 
807 void Foam::edgeCollapser::checkBoundaryPointMergeEdges
808 (
809  const label pointi,
810  const label otherPointi,
811  const labelList& pointPriority,
812  Map<point>& collapsePointToLocation
813 ) const
814 {
815  const pointField& points = mesh_.points();
816 
817  const label e0Priority = pointPriority[pointi];
818  const label e1Priority = pointPriority[otherPointi];
819 
820  if (e0Priority > e1Priority)
821  {
822  collapsePointToLocation.set
823  (
824  otherPointi,
825  points[pointi]
826  );
827  }
828  else if (e0Priority < e1Priority)
829  {
830  collapsePointToLocation.set
831  (
832  pointi,
833  points[otherPointi]
834  );
835  }
836  else // e0Priority == e1Priority
837  {
838  collapsePointToLocation.set
839  (
840  pointi,
841  points[otherPointi]
842  );
843 
844 // Foam::point averagePt
845 // (
846 // 0.5*(points[otherPointi] + points[pointi])
847 // );
848 //
849 // collapsePointToLocation.set(pointi, averagePt);
850 // collapsePointToLocation.set(otherPointi, averagePt);
851  }
852 }
853 
854 
855 Foam::label Foam::edgeCollapser::breakStringsAtEdges
856 (
857  const PackedBoolList& markedEdges,
858  PackedBoolList& collapseEdge,
859  List<pointEdgeCollapse>& allPointInfo
860 ) const
861 {
862  const edgeList& edges = mesh_.edges();
863  const labelListList& pointEdges = mesh_.pointEdges();
864 
865  label nUncollapsed = 0;
866 
867  forAll(edges, eI)
868  {
869  if (markedEdges[eI])
870  {
871  const edge& e = edges[eI];
872 
873  const label startCollapseIndex
874  = allPointInfo[e.start()].collapseIndex();
875 
876  if (startCollapseIndex != -1 && startCollapseIndex != -2)
877  {
878  const label endCollapseIndex
879  = allPointInfo[e.end()].collapseIndex();
880 
881  if
882  (
883  !collapseEdge[eI]
884  && startCollapseIndex == endCollapseIndex
885  )
886  {
887  const labelList& ptEdgesStart = pointEdges[e.start()];
888 
889  forAll(ptEdgesStart, ptEdgeI)
890  {
891  const label edgeI = ptEdgesStart[ptEdgeI];
892 
893  const label nbrPointi
894  = edges[edgeI].otherVertex(e.start());
895  const label nbrIndex
896  = allPointInfo[nbrPointi].collapseIndex();
897 
898  if
899  (
900  collapseEdge[edgeI]
901  && nbrIndex == startCollapseIndex
902  )
903  {
904  collapseEdge[edgeI] = false;
905  nUncollapsed++;
906  }
907  }
908  }
909  }
910  }
911  }
912 
913  return nUncollapsed;
914 }
915 
916 
917 void Foam::edgeCollapser::determineDuplicatePointsOnFace
918 (
919  const face& f,
920  PackedBoolList& markedPoints,
921  labelHashSet& uniqueCollapses,
922  labelHashSet& duplicateCollapses,
923  List<pointEdgeCollapse>& allPointInfo
924 ) const
925 {
926  uniqueCollapses.clear();
927  duplicateCollapses.clear();
928 
929  forAll(f, fpI)
930  {
931  label index = allPointInfo[f[fpI]].collapseIndex();
932 
933  // Check for consecutive duplicate
934  if (index != allPointInfo[f.prevLabel(fpI)].collapseIndex())
935  {
936  if (!uniqueCollapses.insert(index))
937  {
938  // Failed inserting so duplicate
939  duplicateCollapses.insert(index);
940  }
941  }
942  }
943 
944  // Now duplicateCollapses contains duplicate collapse indices.
945  // Convert to points.
946  forAll(f, fpI)
947  {
948  label index = allPointInfo[f[fpI]].collapseIndex();
949  if (duplicateCollapses.found(index))
950  {
951  markedPoints[f[fpI]] = true;
952  }
953  }
954 }
955 
956 
957 Foam::label Foam::edgeCollapser::countEdgesOnFace
958 (
959  const face& f,
960  List<pointEdgeCollapse>& allPointInfo
961 ) const
962 {
963  label nEdges = 0;
964 
965  forAll(f, fpI)
966  {
967  const label pointi = f[fpI];
968  const label newPointi = allPointInfo[pointi].collapseIndex();
969 
970  if (newPointi == -2)
971  {
972  nEdges++;
973  }
974  else
975  {
976  const label prevPointi = f[f.fcIndex(fpI)];
977  const label prevNewPointi
978  = allPointInfo[prevPointi].collapseIndex();
979 
980  if (newPointi != prevNewPointi)
981  {
982  nEdges++;
983  }
984  }
985  }
986 
987  return nEdges;
988 }
989 
990 
991 bool Foam::edgeCollapser::isFaceCollapsed
992 (
993  const face& f,
994  List<pointEdgeCollapse>& allPointInfo
995 ) const
996 {
997  label nEdges = countEdgesOnFace(f, allPointInfo);
998 
999  // Polygons must have 3 or more edges to be valid
1000  if (nEdges < 3)
1001  {
1002  return true;
1003  }
1004 
1005  return false;
1006 }
1007 
1008 
1009 // Create consistent set of collapses.
1010 // collapseEdge : per edge:
1011 // -1 : do not collapse
1012 // 0 : collapse to start
1013 // 1 : collapse to end
1014 // Note: collapseEdge has to be parallel consistent (in orientation)
1015 Foam::label Foam::edgeCollapser::syncCollapse
1016 (
1017  const globalIndex& globalPoints,
1018  const labelList& pointPriority,
1019  const PackedBoolList& collapseEdge,
1020  const Map<point>& collapsePointToLocation,
1021  List<pointEdgeCollapse>& allPointInfo
1022 ) const
1023 {
1024  const edgeList& edges = mesh_.edges();
1025 
1026  label nCollapsed = 0;
1027 
1028  DynamicList<label> initPoints(mesh_.nPoints());
1029  DynamicList<pointEdgeCollapse> initPointInfo(mesh_.nPoints());
1030 
1031  allPointInfo.clear();
1032  allPointInfo.setSize(mesh_.nPoints());
1033 
1034  // Initialise edges to no collapse
1035  List<pointEdgeCollapse> allEdgeInfo
1036  (
1037  mesh_.nEdges(),
1038  pointEdgeCollapse(Zero, -1, -1)
1039  );
1040 
1041  // Mark selected edges for collapse
1042  forAll(edges, edgeI)
1043  {
1044  if (collapseEdge[edgeI])
1045  {
1046  const edge& e = edges[edgeI];
1047 
1048  label masterPointi = e.start();
1049 
1050  // Choose the point on the edge with the highest priority.
1051  if (pointPriority[e.end()] > pointPriority[e.start()])
1052  {
1053  masterPointi = e.end();
1054  }
1055 
1056  label masterPointPriority = pointPriority[masterPointi];
1057 
1058  label index = globalPoints.toGlobal(masterPointi);
1059 
1060  if (!collapsePointToLocation.found(masterPointi))
1061  {
1062  const label otherVertex = e.otherVertex(masterPointi);
1063 
1064  if (!collapsePointToLocation.found(otherVertex))
1065  {
1067  << masterPointi << " on edge " << edgeI << " " << e
1068  << " is not marked for collapse."
1069  << abort(FatalError);
1070  }
1071  else
1072  {
1073  masterPointi = otherVertex;
1074  masterPointPriority = pointPriority[masterPointi];
1075  index = globalPoints.toGlobal(masterPointi);
1076  }
1077  }
1078 
1079  const point& collapsePoint = collapsePointToLocation[masterPointi];
1080 
1081  const pointEdgeCollapse pec
1082  (
1083  collapsePoint,
1084  index,
1085  masterPointPriority
1086  );
1087 
1088  // Mark as collapsable but with nonsense master so it gets
1089  // overwritten and starts an update wave
1090  allEdgeInfo[edgeI] = pointEdgeCollapse
1091  (
1092  collapsePoint,
1093  labelMax,
1094  labelMin
1095  );
1096 
1097  initPointInfo.append(pec);
1098  initPoints.append(e.start());
1099 
1100  initPointInfo.append(pec);
1101  initPoints.append(e.end());
1102 
1103  nCollapsed++;
1104  }
1105  }
1106 
1107  PointEdgeWave<pointEdgeCollapse> collapsePropagator
1108  (
1109  mesh_,
1110  initPoints,
1111  initPointInfo,
1112  allPointInfo,
1113  allEdgeInfo,
1114  mesh_.globalData().nTotalPoints() // Maximum number of iterations
1115  );
1116 
1117  return nCollapsed;
1118 }
1119 
1120 
1121 void Foam::edgeCollapser::filterFace
1122 (
1123  const Map<DynamicList<label>>& collapseStrings,
1124  const List<pointEdgeCollapse>& allPointInfo,
1125  face& f
1126 ) const
1127 {
1128  label newFp = 0;
1129 
1130  face oldFace = f;
1131 
1132  forAll(f, fp)
1133  {
1134  label pointi = f[fp];
1135 
1136  label collapseIndex = allPointInfo[pointi].collapseIndex();
1137 
1138  // Do we have a local point for this index?
1139  if (collapseStrings.found(collapseIndex))
1140  {
1141  label localPointi = collapseStrings[collapseIndex][0];
1142 
1143  if (findIndex(SubList<label>(f, newFp), localPointi) == -1)
1144  {
1145  f[newFp++] = localPointi;
1146  }
1147  }
1148  else if (collapseIndex == -1)
1149  {
1151  << "Point " << pointi << " was not visited by PointEdgeWave"
1152  << endl;
1153  }
1154  else
1155  {
1156  f[newFp++] = pointi;
1157  }
1158  }
1159 
1160 
1161  // Check for pinched face. Tries to correct
1162  // - consecutive duplicate vertex. Removes duplicate vertex.
1163  // - duplicate vertex with one other vertex in between (spike).
1164  // Both of these should not really occur! and should be checked before
1165  // collapsing edges.
1166 
1167  const label size = newFp;
1168 
1169  newFp = 2;
1170 
1171  for (label fp = 2; fp < size; fp++)
1172  {
1173  label fp1 = fp-1;
1174  label fp2 = fp-2;
1175 
1176  label pointi = f[fp];
1177 
1178  // Search for previous occurrence.
1179  label index = findIndex(SubList<label>(f, fp), pointi);
1180 
1181  if (index == fp1)
1182  {
1184  << "Removing consecutive duplicate vertex in face "
1185  << f << endl;
1186  // Don't store current pointi
1187  }
1188  else if (index == fp2)
1189  {
1191  << "Removing non-consecutive duplicate vertex in face "
1192  << f << endl;
1193  // Don't store current pointi and remove previous
1194  newFp--;
1195  }
1196  else if (index != -1)
1197  {
1199  << "Pinched face " << f << endl;
1200  f[newFp++] = pointi;
1201  }
1202  else
1203  {
1204  f[newFp++] = pointi;
1205  }
1206  }
1207 
1208  f.setSize(newFp);
1209 }
1210 
1211 
1212 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1213 
1214 Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh)
1215 :
1216  mesh_(mesh),
1217  guardFraction_(0),
1218  maxCollapseFaceToPointSideLengthCoeff_(0),
1219  allowEarlyCollapseToPoint_(false),
1220  allowEarlyCollapseCoeff_(0)
1221 {}
1222 
1223 
1224 Foam::edgeCollapser::edgeCollapser
1226  const polyMesh& mesh,
1227  const dictionary& dict
1228 )
1229 :
1230  mesh_(mesh),
1231  guardFraction_
1232  (
1233  dict.lookupOrDefault<scalar>("guardFraction", 0)
1234  ),
1235  maxCollapseFaceToPointSideLengthCoeff_
1236  (
1237  dict.lookupOrDefault<scalar>("maxCollapseFaceToPointSideLengthCoeff", 0)
1238  ),
1239  allowEarlyCollapseToPoint_
1240  (
1241  dict.lookupOrDefault<Switch>("allowEarlyCollapseToPoint", true)
1242  ),
1243  allowEarlyCollapseCoeff_
1244  (
1245  dict.lookupOrDefault<scalar>("allowEarlyCollapseCoeff", 0)
1246  )
1247 {
1248  if (debug)
1249  {
1250  Info<< "Edge Collapser Settings:" << nl
1251  << " Guard Fraction = " << guardFraction_ << nl
1252  << " Max collapse face to point side length = "
1253  << maxCollapseFaceToPointSideLengthCoeff_ << nl
1254  << " " << (allowEarlyCollapseToPoint_ ? "Allow" : "Do not allow")
1255  << " early collapse to point" << nl
1256  << " Early collapse coeff = " << allowEarlyCollapseCoeff_
1257  << endl;
1258  }
1259 }
1260 
1261 
1262 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1263 
1266  const List<pointEdgeCollapse>& allPointInfo,
1267  polyTopoChange& meshMod
1268 ) const
1269 {
1270  const cellList& cells = mesh_.cells();
1271  const labelList& faceOwner = mesh_.faceOwner();
1272  const labelList& faceNeighbour = mesh_.faceNeighbour();
1273  const labelListList& pointFaces = mesh_.pointFaces();
1274  const pointZoneMesh& pointZones = mesh_.pointZones();
1275 
1276 
1277 
1278 
1279 // // Dump point collapses
1280 // label count = 0;
1281 // forAll(allPointInfo, ptI)
1282 // {
1283 // const pointEdgeCollapse& pec = allPointInfo[ptI];
1284 //
1285 // if (mesh_.points()[ptI] != pec.collapsePoint())
1286 // {
1287 // count++;
1288 // }
1289 // }
1290 //
1291 // OFstream str("collapses_" + name(count) + ".obj");
1292 // // Dump point collapses
1293 // forAll(allPointInfo, ptI)
1294 // {
1295 // const pointEdgeCollapse& pec = allPointInfo[ptI];
1296 //
1297 // if
1298 // (
1299 // mesh_.points()[ptI] != pec.collapsePoint()
1300 // && pec.collapsePoint() != vector(GREAT, GREAT, GREAT)
1301 // )
1302 // {
1303 // meshTools::writeOBJ
1304 // (
1305 // str,
1306 // mesh_.points()[ptI],
1307 // pec.collapsePoint()
1308 // );
1309 // }
1310 // }
1311 
1312 
1313 
1314  bool meshChanged = false;
1315 
1316  PackedBoolList removedPoints(mesh_.nPoints());
1317 
1318  // Create strings of edges.
1319  // Map from collapseIndex(=global master point) to set of points
1320  Map<DynamicList<label>> collapseStrings;
1321  {
1322  // 1. Count elements per collapseIndex
1323  Map<label> nPerIndex(mesh_.nPoints()/10);
1324  forAll(allPointInfo, pointi)
1325  {
1326  label collapseIndex = allPointInfo[pointi].collapseIndex();
1327 
1328  if (collapseIndex != -1 && collapseIndex != -2)
1329  {
1330  Map<label>::iterator fnd = nPerIndex.find(collapseIndex);
1331  if (fnd != nPerIndex.end())
1332  {
1333  fnd()++;
1334  }
1335  else
1336  {
1337  nPerIndex.insert(collapseIndex, 1);
1338  }
1339  }
1340  }
1341 
1342  // 2. Size
1343  collapseStrings.resize(2*nPerIndex.size());
1344  forAllConstIter(Map<label>, nPerIndex, iter)
1345  {
1346  collapseStrings.insert(iter.key(), DynamicList<label>(iter()));
1347  }
1348 
1349  // 3. Fill
1350  forAll(allPointInfo, pointi)
1351  {
1352  const label collapseIndex = allPointInfo[pointi].collapseIndex();
1353 
1354  if (collapseIndex != -1 && collapseIndex != -2)
1355  {
1356  collapseStrings[collapseIndex].append(pointi);
1357  }
1358  }
1359  }
1360 
1361 
1362 
1363 
1364 // OFstream str2("collapseStrings_" + name(count) + ".obj");
1365 // // Dump point collapses
1366 // forAllConstIter(Map<DynamicList<label>>, collapseStrings, iter)
1367 // {
1368 // const label masterPoint = iter.key();
1369 // const DynamicList<label>& edgeCollapses = iter();
1370 //
1371 // forAll(edgeCollapses, eI)
1372 // {
1373 // meshTools::writeOBJ
1374 // (
1375 // str2,
1376 // mesh_.points()[edgeCollapses[eI]],
1377 // mesh_.points()[masterPoint]
1378 // );
1379 // }
1380 // }
1381 
1382 
1383 
1384  // Current faces (is also collapseStatus: f.size() < 3)
1385  faceList newFaces(mesh_.faces());
1386 
1387  // Current cellCollapse status
1388  boolList cellRemoved(mesh_.nCells(), false);
1389 
1390  label nUnvisited = 0;
1391  label nUncollapsed = 0;
1392  label nCollapsed = 0;
1393 
1394  forAll(allPointInfo, pI)
1395  {
1396  const pointEdgeCollapse& pec = allPointInfo[pI];
1397 
1398  if (pec.collapseIndex() == -1)
1399  {
1400  nUnvisited++;
1401  }
1402  else if (pec.collapseIndex() == -2)
1403  {
1404  nUncollapsed++;
1405  }
1406  else
1407  {
1408  nCollapsed++;
1409  }
1410  }
1411 
1412  label nPoints = allPointInfo.size();
1413 
1414  reduce(nPoints, sumOp<label>());
1415  reduce(nUnvisited, sumOp<label>());
1416  reduce(nUncollapsed, sumOp<label>());
1417  reduce(nCollapsed, sumOp<label>());
1418 
1419  Info<< incrIndent;
1420  Info<< indent << "Number of points : " << nPoints << nl
1421  << indent << "Not visited : " << nUnvisited << nl
1422  << indent << "Not collapsed : " << nUncollapsed << nl
1423  << indent << "Collapsed : " << nCollapsed << nl
1424  << endl;
1425  Info<< decrIndent;
1426 
1427  do
1428  {
1429  forAll(newFaces, facei)
1430  {
1431  filterFace(collapseStrings, allPointInfo, newFaces[facei]);
1432  }
1433 
1434  // Check if faces to be collapsed cause cells to become collapsed.
1435  label nCellCollapsed = 0;
1436 
1437  forAll(cells, celli)
1438  {
1439  if (!cellRemoved[celli])
1440  {
1441  const cell& cFaces = cells[celli];
1442 
1443  label nFaces = cFaces.size();
1444 
1445  forAll(cFaces, i)
1446  {
1447  label facei = cFaces[i];
1448 
1449  if (newFaces[facei].size() < 3)
1450  {
1451  --nFaces;
1452 
1453  if (nFaces < 4)
1454  {
1455  Pout<< "Cell:" << celli
1456  << " uses faces:" << cFaces
1457  << " of which too many are marked for removal:"
1458  << endl
1459  << " ";
1460 
1461 
1462  forAll(cFaces, j)
1463  {
1464  if (newFaces[cFaces[j]].size() < 3)
1465  {
1466  Pout<< ' '<< cFaces[j];
1467  }
1468  }
1469  Pout<< endl;
1470 
1471  cellRemoved[celli] = true;
1472 
1473  // Collapse all edges of cell to nothing
1474 // collapseEdges(cellEdges[celli]);
1475 
1476  nCellCollapsed++;
1477 
1478  break;
1479  }
1480  }
1481  }
1482  }
1483  }
1484 
1485  reduce(nCellCollapsed, sumOp<label>());
1486  Info<< indent << "Collapsing " << nCellCollapsed << " cells" << endl;
1487 
1488  if (nCellCollapsed == 0)
1489  {
1490  break;
1491  }
1492 
1493  } while (true);
1494 
1495 
1496  // Keep track of faces that have been done already.
1497  boolList doneFace(mesh_.nFaces(), false);
1498 
1499  {
1500  // Mark points used.
1501  boolList usedPoint(mesh_.nPoints(), false);
1502 
1503  forAll(cellRemoved, celli)
1504  {
1505  if (cellRemoved[celli])
1506  {
1507  meshMod.removeCell(celli, -1);
1508  }
1509  }
1510 
1511  // Remove faces
1512  forAll(newFaces, facei)
1513  {
1514  const face& f = newFaces[facei];
1515 
1516  if (f.size() < 3)
1517  {
1518  meshMod.removeFace(facei, -1);
1519  meshChanged = true;
1520 
1521  // Mark face as been done.
1522  doneFace[facei] = true;
1523  }
1524  else
1525  {
1526  // Kept face. Mark vertices
1527  forAll(f, fp)
1528  {
1529  usedPoint[f[fp]] = true;
1530  }
1531  }
1532  }
1533 
1534  // Remove unused vertices that have not been marked for removal already
1535  forAll(usedPoint, pointi)
1536  {
1537  if (!usedPoint[pointi])
1538  {
1539  removedPoints[pointi] = true;
1540  meshMod.removePoint(pointi, -1);
1541  meshChanged = true;
1542  }
1543  }
1544  }
1545 
1546  // Modify the point location of the remaining points
1547  forAll(allPointInfo, pointi)
1548  {
1549  const label collapseIndex = allPointInfo[pointi].collapseIndex();
1550  const point& collapsePoint = allPointInfo[pointi].collapsePoint();
1551 
1552  if
1553  (
1554  removedPoints[pointi] == false
1555  && collapseIndex != -1
1556  && collapseIndex != -2
1557  )
1558  {
1559  meshMod.modifyPoint
1560  (
1561  pointi,
1562  collapsePoint,
1563  pointZones.whichZone(pointi),
1564  false
1565  );
1566  }
1567  }
1568 
1569 
1570  const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
1571  const faceZoneMesh& faceZones = mesh_.faceZones();
1572 
1573  // Renumber faces that use points
1574  forAll(allPointInfo, pointi)
1575  {
1576  if (removedPoints[pointi] == true)
1577  {
1578  const labelList& changedFaces = pointFaces[pointi];
1579 
1580  forAll(changedFaces, changedFacei)
1581  {
1582  label facei = changedFaces[changedFacei];
1583 
1584  if (!doneFace[facei])
1585  {
1586  doneFace[facei] = true;
1587 
1588  // Get current zone info
1589  label zoneID = faceZones.whichZone(facei);
1590 
1591  bool zoneFlip = false;
1592 
1593  if (zoneID >= 0)
1594  {
1595  const faceZone& fZone = faceZones[zoneID];
1596 
1597  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
1598  }
1599 
1600  // Get current connectivity
1601  label own = faceOwner[facei];
1602  label nei = -1;
1603  label patchID = -1;
1604 
1605  if (mesh_.isInternalFace(facei))
1606  {
1607  nei = faceNeighbour[facei];
1608  }
1609  else
1610  {
1611  patchID = boundaryMesh.whichPatch(facei);
1612  }
1613 
1614  meshMod.modifyFace
1615  (
1616  newFaces[facei], // face
1617  facei, // facei to change
1618  own, // owner
1619  nei, // neighbour
1620  false, // flipFaceFlux
1621  patchID, // patch
1622  zoneID,
1623  zoneFlip
1624  );
1625 
1626  meshChanged = true;
1627  }
1628  }
1629  }
1630  }
1631 
1632  return meshChanged;
1633 }
1634 
1635 
1638  const globalIndex& globalPoints,
1639  const labelList& pointPriority,
1640  const Map<point>& collapsePointToLocation,
1641  PackedBoolList& collapseEdge,
1642  List<pointEdgeCollapse>& allPointInfo,
1643  const bool allowCellCollapse
1644 ) const
1645 {
1646  // Make sure we don't collapse cells
1647  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1648  const faceList faces = mesh_.faces();
1649  const edgeList& edges = mesh_.edges();
1650  const labelListList& faceEdges = mesh_.faceEdges();
1651  const labelListList& pointEdges = mesh_.pointEdges();
1652  const cellList& cells = mesh_.cells();
1653 
1654  labelHashSet uniqueCollapses;
1655  labelHashSet duplicateCollapses;
1656 
1657  while (true)
1658  {
1659  label nUncollapsed = 0;
1660 
1662  (
1663  mesh_,
1664  collapseEdge,
1666  0
1667  );
1668 
1669  // Create consistent set of collapses
1670  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1671  // Note: requires collapseEdge to be synchronised.
1672  syncCollapse
1673  (
1674  globalPoints,
1675  pointPriority,
1676  collapseEdge,
1677  collapsePointToLocation,
1678  allPointInfo
1679  );
1680 
1681  // Get collapsed faces
1682 
1683  PackedBoolList isCollapsedFace(mesh_.nFaces());
1684  PackedBoolList markedPoints(mesh_.nPoints());
1685 
1686  forAll(faces, facei)
1687  {
1688  const face& f = faces[facei];
1689 
1690  isCollapsedFace[facei] = isFaceCollapsed(f, allPointInfo);
1691 
1692  if (isCollapsedFace[facei] < 1)
1693  {
1694  determineDuplicatePointsOnFace
1695  (
1696  f,
1697  markedPoints,
1698  uniqueCollapses,
1699  duplicateCollapses,
1700  allPointInfo
1701  );
1702  }
1703  }
1704 
1705  // Synchronise the marked points
1707  (
1708  mesh_,
1709  markedPoints,
1711  0
1712  );
1713 
1714  // Mark all edges attached to the point for collapse
1715  forAll(markedPoints, pointi)
1716  {
1717  if (markedPoints[pointi])
1718  {
1719  const label index = allPointInfo[pointi].collapseIndex();
1720 
1721  const labelList& ptEdges = pointEdges[pointi];
1722 
1723  forAll(ptEdges, ptEdgeI)
1724  {
1725  const label edgeI = ptEdges[ptEdgeI];
1726  const label nbrPointi = edges[edgeI].otherVertex(pointi);
1727  const label nbrIndex
1728  = allPointInfo[nbrPointi].collapseIndex();
1729 
1730  if (collapseEdge[edgeI] && nbrIndex == index)
1731  {
1732  collapseEdge[edgeI] = false;
1733  nUncollapsed++;
1734  }
1735  }
1736  }
1737  }
1738 
1739  PackedBoolList markedEdges(mesh_.nEdges());
1740 
1741  if (!allowCellCollapse)
1742  {
1743  // Check collapsed cells
1744  forAll(cells, celli)
1745  {
1746  const cell& cFaces = cells[celli];
1747 
1748  label nFaces = cFaces.size();
1749 
1750  forAll(cFaces, fI)
1751  {
1752  label facei = cFaces[fI];
1753 
1754  if (isCollapsedFace[facei])
1755  {
1756  nFaces--;
1757  }
1758  }
1759 
1760  if (nFaces < 4)
1761  {
1762  forAll(cFaces, fI)
1763  {
1764  label facei = cFaces[fI];
1765 
1766  const labelList& fEdges = faceEdges[facei];
1767 
1768  // Unmark this face for collapse
1769  forAll(fEdges, fEdgeI)
1770  {
1771  label edgeI = fEdges[fEdgeI];
1772 
1773  if (collapseEdge[edgeI])
1774  {
1775  collapseEdge[edgeI] = false;
1776  nUncollapsed++;
1777  }
1778 
1779  markedEdges[edgeI] = true;
1780  }
1781 
1782  // Uncollapsed this face.
1783  isCollapsedFace[facei] = false;
1784  nFaces++;
1785  }
1786  }
1787 
1788  if (nFaces < 4)
1789  {
1791  << "Cell " << celli << " " << cFaces << nl
1792  << "is " << nFaces << ", "
1793  << "but cell collapse has been disabled."
1794  << abort(FatalError);
1795  }
1796  }
1797  }
1798 
1800  (
1801  mesh_,
1802  markedEdges,
1804  0
1805  );
1806 
1807  nUncollapsed += breakStringsAtEdges
1808  (
1809  markedEdges,
1810  collapseEdge,
1811  allPointInfo
1812  );
1813 
1814  reduce(nUncollapsed, sumOp<label>());
1815 
1816  Info<< " Uncollapsed edges = " << nUncollapsed << " / "
1817  << returnReduce(mesh_.nEdges(), sumOp<label>()) << endl;
1818 
1819  if (nUncollapsed == 0)
1820  {
1821  break;
1822  }
1823  }
1824 }
1825 
1826 
1829  const scalarField& minEdgeLen,
1830  const labelList& pointPriority,
1831  PackedBoolList& collapseEdge,
1832  Map<point>& collapsePointToLocation
1833 ) const
1834 {
1835  // Work out which edges to collapse
1836  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1837 
1838  const pointField& points = mesh_.points();
1839  const edgeList& edges = mesh_.edges();
1840 
1841  label nCollapsed = 0;
1842 
1843  forAll(edges, edgeI)
1844  {
1845  const edge& e = edges[edgeI];
1846 
1847  if (!collapseEdge[edgeI])
1848  {
1849  if (e.mag(points) < minEdgeLen[edgeI])
1850  {
1851  collapseEdge[edgeI] = true;
1852 
1853  label masterPointi = edgeMaster(pointPriority, e);
1854 
1855  if (masterPointi == -1)
1856  {
1857  const point average
1858  = 0.5*(points[e.start()] + points[e.end()]);
1859 
1860  collapsePointToLocation.set(e.start(), average);
1861  }
1862  else
1863  {
1864  const point& collapsePt = points[masterPointi];
1865 
1866  collapsePointToLocation.set(masterPointi, collapsePt);
1867  }
1868 
1869 
1870  nCollapsed++;
1871  }
1872  }
1873  }
1874 
1875  return nCollapsed;
1876 }
1877 
1878 
1881  const scalar maxCos,
1882  const labelList& pointPriority,
1883  PackedBoolList& collapseEdge,
1884  Map<point>& collapsePointToLocation
1885 ) const
1886 {
1887  const edgeList& edges = mesh_.edges();
1888  const pointField& points = mesh_.points();
1889  const labelListList& pointEdges = mesh_.pointEdges();
1890 
1891  // Point removal engine
1892  removePoints pointRemover(mesh_, false);
1893 
1894  // Find out points that can be deleted
1895  boolList pointCanBeDeleted;
1896  label nTotRemove = pointRemover.countPointUsage(maxCos, pointCanBeDeleted);
1897 
1898  // Rework point-to-remove into edge-to-collapse.
1899 
1900  label nCollapsed = 0;
1901 
1902  if (nTotRemove > 0)
1903  {
1904  forAll(pointEdges, pointi)
1905  {
1906  if (pointCanBeDeleted[pointi])
1907  {
1908  const labelList& pEdges = pointEdges[pointi];
1909 
1910  if (pEdges.size() == 2)
1911  {
1912  // Always the case?
1913 
1914  label e0 = pEdges[0];
1915  label e1 = pEdges[1];
1916 
1917  if (!collapseEdge[e0] && !collapseEdge[e1])
1918  {
1919  // Get lengths of both edges and choose the smallest
1920  scalar e0length = mag
1921  (
1922  points[edges[e0][0]] - points[edges[e0][1]]
1923  );
1924 
1925  scalar e1length = mag
1926  (
1927  points[edges[e1][0]] - points[edges[e1][1]]
1928  );
1929 
1930  if (e0length <= e1length)
1931  {
1932  collapseEdge[e0] = true;
1933 
1934  checkBoundaryPointMergeEdges
1935  (
1936  pointi,
1937  edges[e0].otherVertex(pointi),
1938  pointPriority,
1939  collapsePointToLocation
1940  );
1941  }
1942  else
1943  {
1944  collapseEdge[e1] = true;
1945 
1946  checkBoundaryPointMergeEdges
1947  (
1948  pointi,
1949  edges[e1].otherVertex(pointi),
1950  pointPriority,
1951  collapsePointToLocation
1952  );
1953  }
1954 
1955  nCollapsed++;
1956  }
1957  }
1958  }
1959  }
1960  }
1961 
1962  return nCollapsed;
1963 }
1964 
1965 
1968  const scalarField& faceFilterFactor,
1969  const labelList& pointPriority,
1970  PackedBoolList& collapseEdge,
1971  Map<point>& collapsePointToLocation
1972 ) const
1973 {
1974  const faceList& faces = mesh_.faces();
1975 
1976  const scalarField targetFaceSizes = calcTargetFaceSizes();
1977 
1978  // Calculate number of faces that will be collapsed to a point or an edge
1979  label nCollapseToPoint = 0;
1980  label nCollapseToEdge = 0;
1981 
1982  forAll(faces, fI)
1983  {
1984  const face& f = faces[fI];
1985 
1986  if (faceFilterFactor[fI] <= 0)
1987  {
1988  continue;
1989  }
1990 
1991  collapseType flagCollapseFace = collapseFace
1992  (
1993  pointPriority,
1994  f,
1995  fI,
1996  targetFaceSizes[fI],
1997  collapseEdge,
1998  collapsePointToLocation,
1999  faceFilterFactor
2000  );
2001 
2002  if (flagCollapseFace == noCollapse)
2003  {
2004  continue;
2005  }
2006  else if (flagCollapseFace == toPoint)
2007  {
2008  nCollapseToPoint++;
2009  }
2010  else if (flagCollapseFace == toEdge)
2011  {
2012  nCollapseToEdge++;
2013  }
2014  else
2015  {
2017  << "Face is marked to be collapsed to " << flagCollapseFace
2018  << ". Currently can only collapse to point/edge."
2019  << abort(FatalError);
2020  }
2021  }
2022 
2023  return labelPair(nCollapseToPoint, nCollapseToEdge);
2024 }
2025 
2026 
2029  const faceZone& fZone,
2030  const scalarField& faceFilterFactor,
2031  const labelList& pointPriority,
2032  PackedBoolList& collapseEdge,
2033  Map<point>& collapsePointToLocation
2034 ) const
2035 {
2036  const faceList& faces = mesh_.faces();
2037 
2038  const scalarField targetFaceSizes = calcTargetFaceSizes();
2039 
2040  // Calculate number of faces that will be collapsed to a point or an edge
2041  label nCollapseToPoint = 0;
2042  label nCollapseToEdge = 0;
2043 
2044  forAll(faces, fI)
2045  {
2046  if (fZone.whichFace(fI) == -1)
2047  {
2048  continue;
2049  }
2050 
2051  const face& f = faces[fI];
2052 
2053  if (faceFilterFactor[fI] <= 0)
2054  {
2055  continue;
2056  }
2057 
2058  collapseType flagCollapseFace = collapseFace
2059  (
2060  pointPriority,
2061  f,
2062  fI,
2063  targetFaceSizes[fI],
2064  collapseEdge,
2065  collapsePointToLocation,
2066  faceFilterFactor
2067  );
2068 
2069  if (flagCollapseFace == noCollapse)
2070  {
2071  continue;
2072  }
2073  else if (flagCollapseFace == toPoint)
2074  {
2075  nCollapseToPoint++;
2076  }
2077  else if (flagCollapseFace == toEdge)
2078  {
2079  nCollapseToEdge++;
2080  }
2081  else
2082  {
2084  << "Face is marked to be collapsed to " << flagCollapseFace
2085  << ". Currently can only collapse to point/edge."
2086  << abort(FatalError);
2087  }
2088  }
2089 
2090  return labelPair(nCollapseToPoint, nCollapseToEdge);
2091 
2092 // const edgeList& edges = mesh_.edges();
2093 // const pointField& points = mesh_.points();
2094 // const labelListList& edgeFaces = mesh_.edgeFaces();
2095 // const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
2096 //
2097 // forAll(edges, eI)
2098 // {
2099 // const edge& e = edges[eI];
2100 //
2101 // const labelList& eFaces = edgeFaces[eI];
2102 //
2103 // bool keepEdge = false;
2104 //
2105 // label nInternalFaces = 0;
2106 // label nPatchFaces = 0;
2107 // label nIndirectFaces = 0;
2108 //
2109 // bool coupled = false;
2110 //
2111 // forAll(eFaces, eFacei)
2112 // {
2113 // const label eFaceIndex = eFaces[eFacei];
2114 //
2115 // if (mesh_.isInternalFace(eFaceIndex))
2116 // {
2117 // nInternalFaces++;
2118 // }
2119 // else
2120 // {
2121 // const label patchIndex = bMesh.whichPatch(eFaceIndex);
2122 // const polyPatch& pPatch = bMesh[patchIndex];
2123 //
2124 // if (pPatch.coupled())
2125 // {
2126 // coupled = true;
2127 // nInternalFaces++;
2128 // }
2129 // else
2130 // {
2131 // // Keep the edge if an attached face is not in the zone
2132 // if (fZone.whichFace(eFaceIndex) == -1)
2133 // {
2134 // nPatchFaces++;
2135 // }
2136 // else
2137 // {
2138 // nIndirectFaces++;
2139 // }
2140 // }
2141 // }
2142 // }
2143 //
2144 // if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
2145 // {
2146 // Pout<< eFaces.size() << " ("
2147 // << nInternalFaces << "/" << nPatchFaces << "/"
2148 // << nIndirectFaces << ")" << endl;
2149 // }
2150 //
2151 // if
2152 // (
2153 // eFaces.size() == nInternalFaces
2154 // || nIndirectFaces < (coupled ? 1 : 2)
2155 // )
2156 // {
2157 // keepEdge = true;
2158 // }
2159 //
2160 // if (!keepEdge)
2161 // {
2162 // collapseEdge[eI] = true;
2163 //
2164 // const Foam::point collapsePoint =
2165 // 0.5*(points[e.end()] + points[e.start()]);
2166 //
2167 // collapsePointToLocation.insert(e.start(), collapsePoint);
2168 // collapsePointToLocation.insert(e.end(), collapsePoint);
2169 // }
2170 // }
2171 
2172 // OFstream str
2173 // (
2174 // mesh_.time().path()
2175 // /"markedEdges_" + name(collapseEdge.count()) + ".obj"
2176 // );
2177 // label count = 0;
2178 //
2179 // forAll(collapseEdge, eI)
2180 // {
2181 // if (collapseEdge[eI])
2182 // {
2183 // const edge& e = edges[eI];
2184 //
2185 // meshTools::writeOBJ
2186 // (
2187 // str,
2188 // points[e.start()],
2189 // points[e.end()],
2190 // count
2191 // );
2192 // }
2193 // }
2194 }
2195 
2196 
2197 // ************************************************************************* //
const labelListList & pointFaces() const
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
label end() const
Return end vertex label.
Definition: edgeI.H:92
labelPair markSmallSliverFaces(const scalarField &faceFilterFactor, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Find small faces and sliver faces in the mesh and mark the.
point centre(const pointField &) const
Centre point of face.
Definition: face.C:488
A HashTable with keys but without contents.
Definition: HashSet.H:59
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#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
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:253
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:106
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Cmpt & x() const
Definition: VectorI.H:75
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
static HashSet< label > checkBadFaces(const polyMesh &mesh, const dictionary &meshQualityDict)
Calls motionSmoother::checkMesh and returns a set of bad faces.
Definition: edgeCollapser.C:46
tensor inertia(const pointField &, const point &refPt=vector::zero, scalar density=1.0) const
Return the inertia tensor, with optional reference.
Definition: face.C:726
label countPointUsage(const scalar minCos, boolList &pointCanBeDeleted) const
Mark in pointCanBeDeleted the points that can be deleted.
Definition: removePoints.C:145
const vectorField & faceAreas() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
label markMergeEdges(const scalar maxCos, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to merge.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Definition: boundaryMesh.H:59
dimensionedVector eigenValues(const dimensionedTensor &dt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
label prevLabel(const label i) const
Previous vertex on face.
Definition: faceI.H:124
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:103
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
void removePoint(const label, const label)
Remove/merge point.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
patches[0]
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
T & first()
Return the first element of the list.
Definition: UListI.H:114
const cellList & cells() const
void removeFace(const label, const label)
Remove/merge face.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
bool insert(const label &, const T &newElmt)
Insert a new hashedEntry.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:138
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:305
labelPair markFaceZoneEdges(const faceZone &fZone, const scalarField &faceFilterFactor, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Marks edges in the faceZone indirectPatchFaces for collapse.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
A List obtained as a section of another List.
Definition: SubList.H:53
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
const cellShapeList & cells
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:356
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
label nPoints
void clear()
Clear all entries from table.
Definition: HashTable.C:464
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
const Cmpt & y() const
Definition: VectorI.H:81
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
PointFrompoint toPoint(const Foam::point &p)
static const label labelMax
Definition: label.H:62
const labelListList & pointEdges() const
static const zero Zero
Definition: zero.H:91
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
edgeList edges() const
Return edges in face point ordering,.
Definition: face.C:761
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
Foam::polyBoundaryMesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
label start() const
Return start vertex label.
Definition: edgeI.H:81
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpaceI.H:93
static const char nl
Definition: Ostream.H:262
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:163
static label checkMeshQuality(const polyMesh &mesh, const dictionary &meshQualityDict, PackedBoolList &isErrorPoint)
Check mesh and mark points on faces in error.
Definition: edgeCollapser.C:82
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
label markSmallEdges(const scalarField &minEdgeLen, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to collapse.
label nEdges() const
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:457
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label longestEdge(const face &f, const pointField &pts)
Find the longest edge on a face. Face point labels index into pts.
Definition: face.C:869
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label nFaces() const
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void setSize(const label)
Reset size of List.
Definition: List.C:295
bool set(const label &, const T &newElmt)
Assign a new hashedEntry, overwriting existing entries.
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Modify vertices or cell of face.
A bit-packed bool list.
label patchi
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
#define WarningInFunction
Report a warning using Foam::Warning.
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
label newPointi
Definition: readKivaGrid.H:501
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:221
Direct mesh changes based on v1.3 polyTopoChange syntax.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
void modifyPoint(const label, const point &, const label newZoneID, const bool inCell)
Modify coordinate.
label nPoints() const
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
Determines length of string of edges walked to point.
void consistentCollapse(const globalIndex &globalPoints, const labelList &pointPriority, const Map< point > &collapsePointToLocation, PackedBoolList &collapseEdge, List< pointEdgeCollapse > &allPointInfo, const bool allowCellCollapse=false) const
Ensure that the collapse is parallel consistent and update.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Removes selected points from mesh and updates faces using these points.
Definition: removePoints.H:58
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
bool setRefinement(const List< pointEdgeCollapse > &allPointInfo, polyTopoChange &meshMod) const
Play commands into polyTopoChange to create mesh.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
T & last()
Return the last element of the list.
Definition: UListI.H:128
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:224
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
Wave propagation of information through grid. Every iteration information goes through one layer of e...
Definition: PointEdgeWave.H:85
const labelListList & faceEdges() const
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
vector eigenVector(const tensor &, const scalar lambda)
Definition: tensor.C:190
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
scalar mag(const pointField &) const
Magnitude of face area.
Definition: faceI.H:97
static const label labelMin
Definition: label.H:61
void removeCell(const label, const label)
Remove/merge cell.
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49