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-2017 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 = eigenVectors(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 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
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
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:313
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
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:466
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:110
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
edgeList edges() const
Return edges in face point ordering,.
Definition: face.C:761
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 > &)
const labelListList & faceEdges() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1055
static HashSet< label > checkBadFaces(const polyMesh &mesh, const dictionary &meshQualityDict)
Calls motionSmoother::checkMesh and returns a set of bad faces.
Definition: edgeCollapser.C:46
const labelListList & pointEdges() const
label nFaces() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:253
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label nCells() const
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
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 markMergeEdges(const scalar maxCos, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to merge.
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
point centre(const pointField &) const
Centre point of face.
Definition: face.C:488
const cellList & cells() const
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:306
void removePoint(const label, const label)
Remove/merge point.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
patches[0]
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:103
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
void removeFace(const label, const label)
Remove/merge face.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const Cmpt & y() const
Definition: VectorI.H:81
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1011
label prevLabel(const label i) const
Previous vertex on face.
Definition: faceI.H:124
bool insert(const label &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
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
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
const cellShapeList & cells
const pointField & points
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.
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
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: ListI.H:124
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:460
label nPoints
void clear()
Clear all entries from table.
Definition: HashTable.C:468
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1049
scalar mag(const pointField &) const
Magnitude of face area.
Definition: faceI.H:97
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
PointFrompoint toPoint(const Foam::point &p)
static const label labelMax
Definition: label.H:62
static const zero Zero
Definition: zero.H:91
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:163
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1036
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:61
const Cmpt & x() const
Definition: VectorI.H:75
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpaceI.H:93
static const char nl
Definition: Ostream.H:262
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,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
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.
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
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
label nEdges() const
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:221
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label countPointUsage(const scalar minCos, boolList &pointCanBeDeleted) const
Mark in pointCanBeDeleted the points that can be deleted.
Definition: removePoints.C:145
void setSize(const label)
Reset size of List.
Definition: List.C:281
bool set(const label &, const T &newElmt)
Assign a new hashedEntry, overwriting existing entries.
Definition: HashTableI.H:91
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.
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
A bit-packed bool list.
label patchi
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.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
label newPointi
Definition: readKivaGrid.H:501
label markSmallEdges(const scalarField &minEdgeLen, const labelList &pointPriority, PackedBoolList &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to collapse.
label end() const
Return end vertex label.
Definition: edgeI.H:92
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
bool setRefinement(const List< pointEdgeCollapse > &allPointInfo, polyTopoChange &meshMod) const
Play commands into polyTopoChange to create mesh.
const vectorField & faceAreas() const
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
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.
messageStream Info
const labelListList & pointFaces() const
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
Determines length of string of edges walked to point.
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
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
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:224
Wave propagation of information through grid. Every iteration information goes through one layer of e...
Definition: PointEdgeWave.H:85
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
static const label labelMin
Definition: label.H:61
void removeCell(const label, const label)
Remove/merge cell.
label start() const
Return start vertex label.
Definition: edgeI.H:81
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49