edgeCollapser.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "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  PackedBoolList& collapseEdge,
858  List<pointEdgeCollapse>& allPointInfo
859 ) const
860 {
861  const edgeList& edges = mesh_.edges();
862  const labelListList& pointEdges = mesh_.pointEdges();
863 
864  label nUncollapsed = 0;
865 
866  forAll(edges, eI)
867  {
868  const edge& e = edges[eI];
869 
870  const label startCollapseIndex
871  = allPointInfo[e.start()].collapseIndex();
872 
873  if (startCollapseIndex != -1 && startCollapseIndex != -2)
874  {
875  const label endCollapseIndex
876  = allPointInfo[e.end()].collapseIndex();
877 
878  if (!collapseEdge[eI] && startCollapseIndex == endCollapseIndex)
879  {
880  const labelList& ptEdgesStart = pointEdges[e.start()];
881 
882  forAll(ptEdgesStart, ptEdgeI)
883  {
884  const label edgeI = ptEdgesStart[ptEdgeI];
885 
886  const label nbrPointi = edges[edgeI].otherVertex(e.start());
887  const label nbrIndex =
888  allPointInfo[nbrPointi].collapseIndex();
889 
890  if (collapseEdge[edgeI] && nbrIndex == startCollapseIndex)
891  {
892  collapseEdge[edgeI] = false;
893  nUncollapsed++;
894  }
895  }
896  }
897  }
898  }
899 
900  return nUncollapsed;
901 }
902 
903 
904 void Foam::edgeCollapser::determineDuplicatePointsOnFace
905 (
906  const face& f,
907  PackedBoolList& markedPoints,
908  labelHashSet& uniqueCollapses,
909  labelHashSet& duplicateCollapses,
910  List<pointEdgeCollapse>& allPointInfo
911 ) const
912 {
913  uniqueCollapses.clear();
914  duplicateCollapses.clear();
915 
916  forAll(f, fpI)
917  {
918  label index = allPointInfo[f[fpI]].collapseIndex();
919 
920  // Check for consecutive duplicate
921  if (index != allPointInfo[f.prevLabel(fpI)].collapseIndex())
922  {
923  if (!uniqueCollapses.insert(index))
924  {
925  // Failed inserting so duplicate
926  duplicateCollapses.insert(index);
927  }
928  }
929  }
930 
931  // Now duplicateCollapses contains duplicate collapse indices.
932  // Convert to points.
933  forAll(f, fpI)
934  {
935  label index = allPointInfo[f[fpI]].collapseIndex();
936  if (duplicateCollapses.found(index))
937  {
938  markedPoints[f[fpI]] = true;
939  }
940  }
941 }
942 
943 
944 Foam::label Foam::edgeCollapser::countEdgesOnFace
945 (
946  const face& f,
947  List<pointEdgeCollapse>& allPointInfo
948 ) const
949 {
950  label nEdges = 0;
951 
952  forAll(f, fpI)
953  {
954  const label pointi = f[fpI];
955  const label newPointi = allPointInfo[pointi].collapseIndex();
956 
957  if (newPointi == -2)
958  {
959  nEdges++;
960  }
961  else
962  {
963  const label prevPointi = f[f.fcIndex(fpI)];
964  const label prevNewPointi
965  = allPointInfo[prevPointi].collapseIndex();
966 
967  if (newPointi != prevNewPointi)
968  {
969  nEdges++;
970  }
971  }
972  }
973 
974  return nEdges;
975 }
976 
977 
978 bool Foam::edgeCollapser::isFaceCollapsed
979 (
980  const face& f,
981  List<pointEdgeCollapse>& allPointInfo
982 ) const
983 {
984  label nEdges = countEdgesOnFace(f, allPointInfo);
985 
986  // Polygons must have 3 or more edges to be valid
987  if (nEdges < 3)
988  {
989  return true;
990  }
991 
992  return false;
993 }
994 
995 
996 // Create consistent set of collapses.
997 // collapseEdge : per edge:
998 // -1 : do not collapse
999 // 0 : collapse to start
1000 // 1 : collapse to end
1001 // Note: collapseEdge has to be parallel consistent (in orientation)
1002 Foam::label Foam::edgeCollapser::syncCollapse
1003 (
1004  const globalIndex& globalPoints,
1005  const labelList& pointPriority,
1006  const PackedBoolList& collapseEdge,
1007  const Map<point>& collapsePointToLocation,
1008  List<pointEdgeCollapse>& allPointInfo
1009 ) const
1010 {
1011  const edgeList& edges = mesh_.edges();
1012 
1013  label nCollapsed = 0;
1014 
1015  DynamicList<label> initPoints(mesh_.nPoints());
1016  DynamicList<pointEdgeCollapse> initPointInfo(mesh_.nPoints());
1017 
1018  allPointInfo.clear();
1019  allPointInfo.setSize(mesh_.nPoints());
1020 
1021  // Initialise edges to no collapse
1022  List<pointEdgeCollapse> allEdgeInfo
1023  (
1024  mesh_.nEdges(),
1025  pointEdgeCollapse(Zero, -1, -1)
1026  );
1027 
1028  // Mark selected edges for collapse
1029  forAll(edges, edgeI)
1030  {
1031  if (collapseEdge[edgeI])
1032  {
1033  const edge& e = edges[edgeI];
1034 
1035  label masterPointi = e.start();
1036 
1037  // Choose the point on the edge with the highest priority.
1038  if (pointPriority[e.end()] > pointPriority[e.start()])
1039  {
1040  masterPointi = e.end();
1041  }
1042 
1043  label masterPointPriority = pointPriority[masterPointi];
1044 
1045  label index = globalPoints.toGlobal(masterPointi);
1046 
1047  if (!collapsePointToLocation.found(masterPointi))
1048  {
1049  const label otherVertex = e.otherVertex(masterPointi);
1050 
1051  if (!collapsePointToLocation.found(otherVertex))
1052  {
1054  << masterPointi << " on edge " << edgeI << " " << e
1055  << " is not marked for collapse."
1056  << abort(FatalError);
1057  }
1058  else
1059  {
1060  masterPointi = otherVertex;
1061  masterPointPriority = pointPriority[masterPointi];
1062  index = globalPoints.toGlobal(masterPointi);
1063  }
1064  }
1065 
1066  const point& collapsePoint = collapsePointToLocation[masterPointi];
1067 
1068  const pointEdgeCollapse pec
1069  (
1070  collapsePoint,
1071  index,
1072  masterPointPriority
1073  );
1074 
1075  // Mark as collapsible but with nonsense master so it gets
1076  // overwritten and starts an update wave
1077  allEdgeInfo[edgeI] = pointEdgeCollapse
1078  (
1079  collapsePoint,
1080  labelMax,
1081  labelMin
1082  );
1083 
1084  initPointInfo.append(pec);
1085  initPoints.append(e.start());
1086 
1087  initPointInfo.append(pec);
1088  initPoints.append(e.end());
1089 
1090  nCollapsed++;
1091  }
1092  }
1093 
1094  PointEdgeWave<pointEdgeCollapse> collapsePropagator
1095  (
1096  mesh_,
1097  initPoints,
1098  initPointInfo,
1099  allPointInfo,
1100  allEdgeInfo,
1101  mesh_.globalData().nTotalPoints() // Maximum number of iterations
1102  );
1103 
1104  return nCollapsed;
1105 }
1106 
1107 
1108 void Foam::edgeCollapser::filterFace
1109 (
1110  const Map<DynamicList<label>>& collapseStrings,
1111  const List<pointEdgeCollapse>& allPointInfo,
1112  face& f
1113 ) const
1114 {
1115  label newFp = 0;
1116 
1117  face oldFace = f;
1118 
1119  forAll(f, fp)
1120  {
1121  label pointi = f[fp];
1122 
1123  label collapseIndex = allPointInfo[pointi].collapseIndex();
1124 
1125  // Do we have a local point for this index?
1126  if (collapseStrings.found(collapseIndex))
1127  {
1128  label localPointi = collapseStrings[collapseIndex][0];
1129 
1130  if (findIndex(SubList<label>(f, newFp), localPointi) == -1)
1131  {
1132  f[newFp++] = localPointi;
1133  }
1134  }
1135  else if (collapseIndex == -1)
1136  {
1138  << "Point " << pointi << " was not visited by PointEdgeWave"
1139  << endl;
1140  }
1141  else
1142  {
1143  f[newFp++] = pointi;
1144  }
1145  }
1146 
1147 
1148  // Check for pinched face. Tries to correct
1149  // - consecutive duplicate vertex. Removes duplicate vertex.
1150  // - duplicate vertex with one other vertex in between (spike).
1151  // Both of these should not really occur! and should be checked before
1152  // collapsing edges.
1153 
1154  const label size = newFp;
1155 
1156  newFp = 2;
1157 
1158  for (label fp = 2; fp < size; fp++)
1159  {
1160  label fp1 = fp-1;
1161  label fp2 = fp-2;
1162 
1163  label pointi = f[fp];
1164 
1165  // Search for previous occurrence.
1166  label index = findIndex(SubList<label>(f, fp), pointi);
1167 
1168  if (index == fp1)
1169  {
1171  << "Removing consecutive duplicate vertex in face "
1172  << f << endl;
1173  // Don't store current pointi
1174  }
1175  else if (index == fp2)
1176  {
1178  << "Removing non-consecutive duplicate vertex in face "
1179  << f << endl;
1180  // Don't store current pointi and remove previous
1181  newFp--;
1182  }
1183  else if (index != -1)
1184  {
1186  << "Pinched face " << f << endl;
1187  f[newFp++] = pointi;
1188  }
1189  else
1190  {
1191  f[newFp++] = pointi;
1192  }
1193  }
1194 
1195  f.setSize(newFp);
1196 }
1197 
1198 
1199 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1200 
1202 :
1203  mesh_(mesh),
1204  guardFraction_(0),
1205  maxCollapseFaceToPointSideLengthCoeff_(0),
1206  allowEarlyCollapseToPoint_(false),
1207  allowEarlyCollapseCoeff_(0)
1208 {}
1209 
1210 
1213  const polyMesh& mesh,
1214  const dictionary& dict
1215 )
1216 :
1217  mesh_(mesh),
1218  guardFraction_
1219  (
1220  dict.lookupOrDefault<scalar>("guardFraction", 0)
1221  ),
1222  maxCollapseFaceToPointSideLengthCoeff_
1223  (
1224  dict.lookupOrDefault<scalar>("maxCollapseFaceToPointSideLengthCoeff", 0)
1225  ),
1226  allowEarlyCollapseToPoint_
1227  (
1228  dict.lookupOrDefault<Switch>("allowEarlyCollapseToPoint", true)
1229  ),
1230  allowEarlyCollapseCoeff_
1231  (
1232  dict.lookupOrDefault<scalar>("allowEarlyCollapseCoeff", 0)
1233  )
1234 {
1235  if (debug)
1236  {
1237  Info<< "Edge Collapser Settings:" << nl
1238  << " Guard Fraction = " << guardFraction_ << nl
1239  << " Max collapse face to point side length = "
1240  << maxCollapseFaceToPointSideLengthCoeff_ << nl
1241  << " " << (allowEarlyCollapseToPoint_ ? "Allow" : "Do not allow")
1242  << " early collapse to point" << nl
1243  << " Early collapse coeff = " << allowEarlyCollapseCoeff_
1244  << endl;
1245  }
1246 }
1247 
1248 
1249 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1250 
1253  const List<pointEdgeCollapse>& allPointInfo,
1254  polyTopoChange& meshMod
1255 ) const
1256 {
1257  const cellList& cells = mesh_.cells();
1258  const labelList& faceOwner = mesh_.faceOwner();
1259  const labelList& faceNeighbour = mesh_.faceNeighbour();
1260  const labelListList& pointFaces = mesh_.pointFaces();
1261  const meshPointZones& pointZones = mesh_.pointZones();
1262 
1263 
1264 
1265 
1266 // // Dump point collapses
1267 // label count = 0;
1268 // forAll(allPointInfo, ptI)
1269 // {
1270 // const pointEdgeCollapse& pec = allPointInfo[ptI];
1271 //
1272 // if (mesh_.points()[ptI] != pec.collapsePoint())
1273 // {
1274 // count++;
1275 // }
1276 // }
1277 //
1278 // OFstream str("collapses_" + name(count) + ".obj");
1279 // // Dump point collapses
1280 // forAll(allPointInfo, ptI)
1281 // {
1282 // const pointEdgeCollapse& pec = allPointInfo[ptI];
1283 //
1284 // if
1285 // (
1286 // mesh_.points()[ptI] != pec.collapsePoint()
1287 // && pec.collapsePoint() != vector(great, great, great)
1288 // )
1289 // {
1290 // meshTools::writeOBJ
1291 // (
1292 // str,
1293 // mesh_.points()[ptI],
1294 // pec.collapsePoint()
1295 // );
1296 // }
1297 // }
1298 
1299 
1300 
1301  bool meshChanged = false;
1302 
1303  PackedBoolList removedPoints(mesh_.nPoints());
1304 
1305  // Create strings of edges.
1306  // Map from collapseIndex(=global master point) to set of points
1307  Map<DynamicList<label>> collapseStrings;
1308  {
1309  // 1. Count elements per collapseIndex
1310  Map<label> nPerIndex(mesh_.nPoints()/10);
1311  forAll(allPointInfo, pointi)
1312  {
1313  label collapseIndex = allPointInfo[pointi].collapseIndex();
1314 
1315  if (collapseIndex != -1 && collapseIndex != -2)
1316  {
1317  Map<label>::iterator fnd = nPerIndex.find(collapseIndex);
1318  if (fnd != nPerIndex.end())
1319  {
1320  fnd()++;
1321  }
1322  else
1323  {
1324  nPerIndex.insert(collapseIndex, 1);
1325  }
1326  }
1327  }
1328 
1329  // 2. Size
1330  collapseStrings.resize(2*nPerIndex.size());
1331  forAllConstIter(Map<label>, nPerIndex, iter)
1332  {
1333  collapseStrings.insert(iter.key(), DynamicList<label>(iter()));
1334  }
1335 
1336  // 3. Fill
1337  forAll(allPointInfo, pointi)
1338  {
1339  const label collapseIndex = allPointInfo[pointi].collapseIndex();
1340 
1341  if (collapseIndex != -1 && collapseIndex != -2)
1342  {
1343  collapseStrings[collapseIndex].append(pointi);
1344  }
1345  }
1346  }
1347 
1348 
1349 
1350 
1351 // OFstream str2("collapseStrings_" + name(count) + ".obj");
1352 // // Dump point collapses
1353 // forAllConstIter(Map<DynamicList<label>>, collapseStrings, iter)
1354 // {
1355 // const label masterPoint = iter.key();
1356 // const DynamicList<label>& edgeCollapses = iter();
1357 //
1358 // forAll(edgeCollapses, eI)
1359 // {
1360 // meshTools::writeOBJ
1361 // (
1362 // str2,
1363 // mesh_.points()[edgeCollapses[eI]],
1364 // mesh_.points()[masterPoint]
1365 // );
1366 // }
1367 // }
1368 
1369 
1370 
1371  // Current faces (is also collapseStatus: f.size() < 3)
1372  faceList newFaces(mesh_.faces());
1373 
1374  // Current cellCollapse status
1375  boolList cellRemoved(mesh_.nCells(), false);
1376 
1377  label nUnvisited = 0;
1378  label nUncollapsed = 0;
1379  label nCollapsed = 0;
1380 
1381  forAll(allPointInfo, pI)
1382  {
1383  const pointEdgeCollapse& pec = allPointInfo[pI];
1384 
1385  if (pec.collapseIndex() == -1)
1386  {
1387  nUnvisited++;
1388  }
1389  else if (pec.collapseIndex() == -2)
1390  {
1391  nUncollapsed++;
1392  }
1393  else
1394  {
1395  nCollapsed++;
1396  }
1397  }
1398 
1399  label nPoints = allPointInfo.size();
1400 
1401  reduce(nPoints, sumOp<label>());
1402  reduce(nUnvisited, sumOp<label>());
1403  reduce(nUncollapsed, sumOp<label>());
1404  reduce(nCollapsed, sumOp<label>());
1405 
1406  Info<< incrIndent;
1407  Info<< indent << "Number of points : " << nPoints << nl
1408  << indent << "Not visited : " << nUnvisited << nl
1409  << indent << "Not collapsed : " << nUncollapsed << nl
1410  << indent << "Collapsed : " << nCollapsed << nl
1411  << endl;
1412  Info<< decrIndent;
1413 
1414  do
1415  {
1416  forAll(newFaces, facei)
1417  {
1418  filterFace(collapseStrings, allPointInfo, newFaces[facei]);
1419  }
1420 
1421  // Check if faces to be collapsed cause cells to become collapsed.
1422  label nCellCollapsed = 0;
1423 
1424  forAll(cells, celli)
1425  {
1426  if (!cellRemoved[celli])
1427  {
1428  const cell& cFaces = cells[celli];
1429 
1430  label nFaces = cFaces.size();
1431 
1432  forAll(cFaces, i)
1433  {
1434  label facei = cFaces[i];
1435 
1436  if (newFaces[facei].size() < 3)
1437  {
1438  --nFaces;
1439 
1440  if (nFaces < 4)
1441  {
1442  Pout<< "Cell:" << celli
1443  << " uses faces:" << cFaces
1444  << " of which too many are marked for removal:"
1445  << endl
1446  << " ";
1447 
1448 
1449  forAll(cFaces, j)
1450  {
1451  if (newFaces[cFaces[j]].size() < 3)
1452  {
1453  Pout<< ' '<< cFaces[j];
1454  }
1455  }
1456  Pout<< endl;
1457 
1458  cellRemoved[celli] = true;
1459 
1460  // Collapse all edges of cell to nothing
1461 // collapseEdges(cellEdges[celli]);
1462 
1463  nCellCollapsed++;
1464 
1465  break;
1466  }
1467  }
1468  }
1469  }
1470  }
1471 
1472  reduce(nCellCollapsed, sumOp<label>());
1473  Info<< indent << "Collapsing " << nCellCollapsed << " cells" << endl;
1474 
1475  if (nCellCollapsed == 0)
1476  {
1477  break;
1478  }
1479 
1480  } while (true);
1481 
1482 
1483  // Keep track of faces that have been done already.
1484  boolList doneFace(mesh_.nFaces(), false);
1485 
1486  {
1487  // Mark points used.
1488  boolList usedPoint(mesh_.nPoints(), false);
1489 
1490  forAll(cellRemoved, celli)
1491  {
1492  if (cellRemoved[celli])
1493  {
1494  meshMod.removeCell(celli, -1);
1495  }
1496  }
1497 
1498  // Remove faces
1499  forAll(newFaces, facei)
1500  {
1501  const face& f = newFaces[facei];
1502 
1503  if (f.size() < 3)
1504  {
1505  meshMod.removeFace(facei, -1);
1506  meshChanged = true;
1507 
1508  // Mark face as been done.
1509  doneFace[facei] = true;
1510  }
1511  else
1512  {
1513  // Kept face. Mark vertices
1514  forAll(f, fp)
1515  {
1516  usedPoint[f[fp]] = true;
1517  }
1518  }
1519  }
1520 
1521  // Remove unused vertices that have not been marked for removal already
1522  forAll(usedPoint, pointi)
1523  {
1524  if (!usedPoint[pointi])
1525  {
1526  removedPoints[pointi] = true;
1527  meshMod.removePoint(pointi, -1);
1528  meshChanged = true;
1529  }
1530  }
1531  }
1532 
1533  // Modify the point location of the remaining points
1534  forAll(allPointInfo, pointi)
1535  {
1536  const label collapseIndex = allPointInfo[pointi].collapseIndex();
1537  const point& collapsePoint = allPointInfo[pointi].collapsePoint();
1538 
1539  if
1540  (
1541  removedPoints[pointi] == false
1542  && collapseIndex != -1
1543  && collapseIndex != -2
1544  )
1545  {
1546  meshMod.modifyPoint
1547  (
1548  pointi,
1549  collapsePoint,
1550  pointZones.whichZone(pointi),
1551  false
1552  );
1553  }
1554  }
1555 
1556 
1557  const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
1558  const meshFaceZones& faceZones = mesh_.faceZones();
1559 
1560  // Renumber faces that use points
1561  forAll(allPointInfo, pointi)
1562  {
1563  if (removedPoints[pointi] == true)
1564  {
1565  const labelList& changedFaces = pointFaces[pointi];
1566 
1567  forAll(changedFaces, changedFacei)
1568  {
1569  label facei = changedFaces[changedFacei];
1570 
1571  if (!doneFace[facei])
1572  {
1573  doneFace[facei] = true;
1574 
1575  // Get current zone info
1576  label zoneID = faceZones.whichZone(facei);
1577 
1578  bool zoneFlip = false;
1579 
1580  if (zoneID >= 0)
1581  {
1582  const faceZone& fZone = faceZones[zoneID];
1583 
1584  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
1585  }
1586 
1587  // Get current connectivity
1588  label own = faceOwner[facei];
1589  label nei = -1;
1590  label patchID = -1;
1591 
1592  if (mesh_.isInternalFace(facei))
1593  {
1594  nei = faceNeighbour[facei];
1595  }
1596  else
1597  {
1598  patchID = boundaryMesh.whichPatch(facei);
1599  }
1600 
1601  meshMod.modifyFace
1602  (
1603  newFaces[facei], // face
1604  facei, // facei to change
1605  own, // owner
1606  nei, // neighbour
1607  false, // flipFaceFlux
1608  patchID, // patch
1609  zoneID,
1610  zoneFlip
1611  );
1612 
1613  meshChanged = true;
1614  }
1615  }
1616  }
1617  }
1618 
1619  return meshChanged;
1620 }
1621 
1622 
1625  const globalIndex& globalPoints,
1626  const labelList& pointPriority,
1627  const Map<point>& collapsePointToLocation,
1628  PackedBoolList& collapseEdge,
1629  List<pointEdgeCollapse>& allPointInfo,
1630  const bool allowCellCollapse
1631 ) const
1632 {
1633  // Make sure we don't collapse cells
1634  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1635  const faceList faces = mesh_.faces();
1636  const edgeList& edges = mesh_.edges();
1637  const labelListList& faceEdges = mesh_.faceEdges();
1638  const labelListList& pointEdges = mesh_.pointEdges();
1639  const cellList& cells = mesh_.cells();
1640 
1641  labelHashSet uniqueCollapses;
1642  labelHashSet duplicateCollapses;
1643 
1644  while (true)
1645  {
1646  label nUncollapsed = 0;
1647 
1649  (
1650  mesh_,
1651  collapseEdge,
1653  0
1654  );
1655 
1656  // Create consistent set of collapses
1657  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1658  // Note: requires collapseEdge to be synchronised.
1659  syncCollapse
1660  (
1661  globalPoints,
1662  pointPriority,
1663  collapseEdge,
1664  collapsePointToLocation,
1665  allPointInfo
1666  );
1667 
1668  // Get collapsed faces
1669 
1670  PackedBoolList isCollapsedFace(mesh_.nFaces());
1671  PackedBoolList markedPoints(mesh_.nPoints());
1672 
1673  forAll(faces, facei)
1674  {
1675  const face& f = faces[facei];
1676 
1677  isCollapsedFace[facei] = isFaceCollapsed(f, allPointInfo);
1678 
1679  if (isCollapsedFace[facei] < 1)
1680  {
1681  determineDuplicatePointsOnFace
1682  (
1683  f,
1684  markedPoints,
1685  uniqueCollapses,
1686  duplicateCollapses,
1687  allPointInfo
1688  );
1689  }
1690  }
1691 
1692  // Synchronise the marked points
1694  (
1695  mesh_,
1696  markedPoints,
1698  0
1699  );
1700 
1701  // Mark all edges attached to the point for collapse
1702  forAll(markedPoints, pointi)
1703  {
1704  if (markedPoints[pointi])
1705  {
1706  const label index = allPointInfo[pointi].collapseIndex();
1707 
1708  const labelList& ptEdges = pointEdges[pointi];
1709 
1710  forAll(ptEdges, ptEdgeI)
1711  {
1712  const label edgeI = ptEdges[ptEdgeI];
1713  const label nbrPointi = edges[edgeI].otherVertex(pointi);
1714  const label nbrIndex
1715  = allPointInfo[nbrPointi].collapseIndex();
1716 
1717  if (collapseEdge[edgeI] && nbrIndex == index)
1718  {
1719  collapseEdge[edgeI] = false;
1720  nUncollapsed++;
1721  }
1722  }
1723  }
1724  }
1725 
1726  if (!allowCellCollapse)
1727  {
1728  // Check collapsed cells
1729  forAll(cells, celli)
1730  {
1731  const cell& cFaces = cells[celli];
1732 
1733  label nFaces = cFaces.size();
1734 
1735  forAll(cFaces, fI)
1736  {
1737  label facei = cFaces[fI];
1738 
1739  if (isCollapsedFace[facei])
1740  {
1741  nFaces--;
1742  }
1743  }
1744 
1745  if (nFaces < 4)
1746  {
1747  forAll(cFaces, fI)
1748  {
1749  label facei = cFaces[fI];
1750 
1751  const labelList& fEdges = faceEdges[facei];
1752 
1753  // Unmark this face for collapse
1754  forAll(fEdges, fEdgeI)
1755  {
1756  label edgeI = fEdges[fEdgeI];
1757 
1758  if (collapseEdge[edgeI])
1759  {
1760  collapseEdge[edgeI] = false;
1761  nUncollapsed++;
1762  }
1763  }
1764 
1765  nFaces += isCollapsedFace[facei] ? 1 : 0;
1766 
1767  // Uncollapsed this face.
1768  isCollapsedFace[facei] = false;
1769  }
1770  }
1771 
1772  if (nFaces < 4)
1773  {
1775  << "Cell " << celli << " " << cFaces << nl
1776  << "is " << nFaces << ", "
1777  << "but cell collapse has been disabled."
1778  << abort(FatalError);
1779  }
1780  }
1781  }
1782 
1783  nUncollapsed += breakStringsAtEdges(collapseEdge, allPointInfo);
1784 
1785  reduce(nUncollapsed, sumOp<label>());
1786 
1787  Info<< " Uncollapsed edges = " << nUncollapsed << " / "
1788  << returnReduce(mesh_.nEdges(), sumOp<label>()) << endl;
1789 
1790  if (nUncollapsed == 0)
1791  {
1792  break;
1793  }
1794  }
1795 }
1796 
1797 
1800  const scalarField& minEdgeLen,
1801  const labelList& pointPriority,
1802  PackedBoolList& collapseEdge,
1803  Map<point>& collapsePointToLocation
1804 ) const
1805 {
1806  // Work out which edges to collapse
1807  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1808 
1809  const pointField& points = mesh_.points();
1810  const edgeList& edges = mesh_.edges();
1811 
1812  label nCollapsed = 0;
1813 
1814  forAll(edges, edgeI)
1815  {
1816  const edge& e = edges[edgeI];
1817 
1818  if (!collapseEdge[edgeI])
1819  {
1820  if (e.mag(points) < minEdgeLen[edgeI])
1821  {
1822  collapseEdge[edgeI] = true;
1823 
1824  label masterPointi = edgeMaster(pointPriority, e);
1825 
1826  if (masterPointi == -1)
1827  {
1828  const point average
1829  = 0.5*(points[e.start()] + points[e.end()]);
1830 
1831  collapsePointToLocation.set(e.start(), average);
1832  }
1833  else
1834  {
1835  const point& collapsePt = points[masterPointi];
1836 
1837  collapsePointToLocation.set(masterPointi, collapsePt);
1838  }
1839 
1840 
1841  nCollapsed++;
1842  }
1843  }
1844  }
1845 
1846  return nCollapsed;
1847 }
1848 
1849 
1852  const scalar maxCos,
1853  const labelList& pointPriority,
1854  PackedBoolList& collapseEdge,
1855  Map<point>& collapsePointToLocation
1856 ) const
1857 {
1858  const edgeList& edges = mesh_.edges();
1859  const pointField& points = mesh_.points();
1860  const labelListList& pointEdges = mesh_.pointEdges();
1861 
1862  // Point removal engine
1863  removePoints pointRemover(mesh_, false);
1864 
1865  // Find out points that can be deleted
1866  boolList pointCanBeDeleted;
1867  label nTotRemove = pointRemover.countPointUsage(maxCos, pointCanBeDeleted);
1868 
1869  // Rework point-to-remove into edge-to-collapse.
1870 
1871  label nCollapsed = 0;
1872 
1873  if (nTotRemove > 0)
1874  {
1875  forAll(pointEdges, pointi)
1876  {
1877  if (pointCanBeDeleted[pointi])
1878  {
1879  const labelList& pEdges = pointEdges[pointi];
1880 
1881  if (pEdges.size() == 2)
1882  {
1883  // Always the case?
1884 
1885  label e0 = pEdges[0];
1886  label e1 = pEdges[1];
1887 
1888  if (!collapseEdge[e0] && !collapseEdge[e1])
1889  {
1890  // Get lengths of both edges and choose the smallest
1891  scalar e0length = mag
1892  (
1893  points[edges[e0][0]] - points[edges[e0][1]]
1894  );
1895 
1896  scalar e1length = mag
1897  (
1898  points[edges[e1][0]] - points[edges[e1][1]]
1899  );
1900 
1901  if (e0length <= e1length)
1902  {
1903  collapseEdge[e0] = true;
1904 
1905  checkBoundaryPointMergeEdges
1906  (
1907  pointi,
1908  edges[e0].otherVertex(pointi),
1909  pointPriority,
1910  collapsePointToLocation
1911  );
1912  }
1913  else
1914  {
1915  collapseEdge[e1] = true;
1916 
1917  checkBoundaryPointMergeEdges
1918  (
1919  pointi,
1920  edges[e1].otherVertex(pointi),
1921  pointPriority,
1922  collapsePointToLocation
1923  );
1924  }
1925 
1926  nCollapsed++;
1927  }
1928  }
1929  }
1930  }
1931  }
1932 
1933  return nCollapsed;
1934 }
1935 
1936 
1939  const scalarField& faceFilterFactor,
1940  const labelList& pointPriority,
1941  PackedBoolList& collapseEdge,
1942  Map<point>& collapsePointToLocation
1943 ) const
1944 {
1945  const faceList& faces = mesh_.faces();
1946 
1947  const scalarField targetFaceSizes = calcTargetFaceSizes();
1948 
1949  // Calculate number of faces that will be collapsed to a point or an edge
1950  label nCollapseToPoint = 0;
1951  label nCollapseToEdge = 0;
1952 
1953  forAll(faces, fI)
1954  {
1955  const face& f = faces[fI];
1956 
1957  if (faceFilterFactor[fI] <= 0)
1958  {
1959  continue;
1960  }
1961 
1962  collapseType flagCollapseFace = collapseFace
1963  (
1964  pointPriority,
1965  f,
1966  fI,
1967  targetFaceSizes[fI],
1968  collapseEdge,
1969  collapsePointToLocation,
1970  faceFilterFactor
1971  );
1972 
1973  if (flagCollapseFace == noCollapse)
1974  {
1975  continue;
1976  }
1977  else if (flagCollapseFace == toPoint)
1978  {
1979  nCollapseToPoint++;
1980  }
1981  else if (flagCollapseFace == toEdge)
1982  {
1983  nCollapseToEdge++;
1984  }
1985  else
1986  {
1988  << "Face is marked to be collapsed to " << flagCollapseFace
1989  << ". Currently can only collapse to point/edge."
1990  << abort(FatalError);
1991  }
1992  }
1993 
1994  return labelPair(nCollapseToPoint, nCollapseToEdge);
1995 }
1996 
1997 
2000  const faceZone& fZone,
2001  const scalarField& faceFilterFactor,
2002  const labelList& pointPriority,
2003  PackedBoolList& collapseEdge,
2004  Map<point>& collapsePointToLocation
2005 ) const
2006 {
2007  const faceList& faces = mesh_.faces();
2008 
2009  const scalarField targetFaceSizes = calcTargetFaceSizes();
2010 
2011  // Calculate number of faces that will be collapsed to a point or an edge
2012  label nCollapseToPoint = 0;
2013  label nCollapseToEdge = 0;
2014 
2015  forAll(faces, fI)
2016  {
2017  if (fZone.whichFace(fI) == -1)
2018  {
2019  continue;
2020  }
2021 
2022  const face& f = faces[fI];
2023 
2024  if (faceFilterFactor[fI] <= 0)
2025  {
2026  continue;
2027  }
2028 
2029  collapseType flagCollapseFace = collapseFace
2030  (
2031  pointPriority,
2032  f,
2033  fI,
2034  targetFaceSizes[fI],
2035  collapseEdge,
2036  collapsePointToLocation,
2037  faceFilterFactor
2038  );
2039 
2040  if (flagCollapseFace == noCollapse)
2041  {
2042  continue;
2043  }
2044  else if (flagCollapseFace == toPoint)
2045  {
2046  nCollapseToPoint++;
2047  }
2048  else if (flagCollapseFace == toEdge)
2049  {
2050  nCollapseToEdge++;
2051  }
2052  else
2053  {
2055  << "Face is marked to be collapsed to " << flagCollapseFace
2056  << ". Currently can only collapse to point/edge."
2057  << abort(FatalError);
2058  }
2059  }
2060 
2061  return labelPair(nCollapseToPoint, nCollapseToEdge);
2062 
2063 // const edgeList& edges = mesh_.edges();
2064 // const pointField& points = mesh_.points();
2065 // const labelListList& edgeFaces = mesh_.edgeFaces();
2066 // const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
2067 //
2068 // forAll(edges, eI)
2069 // {
2070 // const edge& e = edges[eI];
2071 //
2072 // const labelList& eFaces = edgeFaces[eI];
2073 //
2074 // bool keepEdge = false;
2075 //
2076 // label nInternalFaces = 0;
2077 // label nPatchFaces = 0;
2078 // label nIndirectFaces = 0;
2079 //
2080 // bool coupled = false;
2081 //
2082 // forAll(eFaces, eFacei)
2083 // {
2084 // const label eFaceIndex = eFaces[eFacei];
2085 //
2086 // if (mesh_.isInternalFace(eFaceIndex))
2087 // {
2088 // nInternalFaces++;
2089 // }
2090 // else
2091 // {
2092 // const label patchIndex = bMesh.whichPatch(eFaceIndex);
2093 // const polyPatch& pPatch = bMesh[patchIndex];
2094 //
2095 // if (pPatch.coupled())
2096 // {
2097 // coupled = true;
2098 // nInternalFaces++;
2099 // }
2100 // else
2101 // {
2102 // // Keep the edge if an attached face is not in the zone
2103 // if (fZone.whichFace(eFaceIndex) == -1)
2104 // {
2105 // nPatchFaces++;
2106 // }
2107 // else
2108 // {
2109 // nIndirectFaces++;
2110 // }
2111 // }
2112 // }
2113 // }
2114 //
2115 // if (eFaces.size() != nInternalFaces + nPatchFaces + nIndirectFaces)
2116 // {
2117 // Pout<< eFaces.size() << " ("
2118 // << nInternalFaces << "/" << nPatchFaces << "/"
2119 // << nIndirectFaces << ")" << endl;
2120 // }
2121 //
2122 // if
2123 // (
2124 // eFaces.size() == nInternalFaces
2125 // || nIndirectFaces < (coupled ? 1 : 2)
2126 // )
2127 // {
2128 // keepEdge = true;
2129 // }
2130 //
2131 // if (!keepEdge)
2132 // {
2133 // collapseEdge[eI] = true;
2134 //
2135 // const Foam::point collapsePoint =
2136 // 0.5*(points[e.end()] + points[e.start()]);
2137 //
2138 // collapsePointToLocation.insert(e.start(), collapsePoint);
2139 // collapsePointToLocation.insert(e.end(), collapsePoint);
2140 // }
2141 // }
2142 
2143 // OFstream str
2144 // (
2145 // mesh_.time().path()
2146 // /"markedEdges_" + name(collapseEdge.count()) + ".obj"
2147 // );
2148 // label count = 0;
2149 //
2150 // forAll(collapseEdge, eI)
2151 // {
2152 // if (collapseEdge[eI])
2153 // {
2154 // const edge& e = edges[eI];
2155 //
2156 // meshTools::writeOBJ
2157 // (
2158 // str,
2159 // points[e.start()],
2160 // points[e.end()],
2161 // count
2162 // );
2163 // }
2164 // }
2165 }
2166 
2167 
2168 // ************************************************************************* //
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 fvPatchList & patches
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
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:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
edgeCollapser(const polyMesh &mesh)
Construct from mesh.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
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:406
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
const labelListList & faceEdges() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
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:164
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label nCells() const
dimensionedVector eigenValues(const dimensionedTensor &dt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
tensor inertia(const pointField &, const point &refPt=vector::zero, scalar density=1.0) const
Return the inertia tensor, with optional reference.
Definition: face.C:371
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/any.
Definition: Switch.H:60
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
const cellList & cells() const
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:307
void removePoint(const label, const label)
Remove/merge point.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label otherVertex(const label a) const
Given one vertex, return the other.
Definition: edgeI.H:140
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
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
const meshPointZones & pointZones() const
Return point zones.
Definition: polyMesh.H:489
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
label prevLabel(const label i) const
Previous vertex on face.
Definition: faceI.H:106
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:316
const cellShapeList & cells
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
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:340
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)
Synchronise values on all mesh edges.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
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:178
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: faceI.H:81
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
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:97
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:181
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
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
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:83
static const char nl
Definition: Ostream.H:260
static label checkMeshQuality(const polyMesh &mesh, const dictionary &meshQualityDict, PackedBoolList &isErrorPoint)
Check mesh and mark points on faces in error.
Definition: edgeCollapser.C:82
static vector centre(const PointField &ps)
Return centre point given face points.
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:235
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence 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:467
label nEdges() const
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
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: MeshZones.C:221
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:57
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 meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
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:306
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise 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:76
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:436
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
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:311
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
Wave propagation of information through grid. Every iteration information goes through one layer of e...
Definition: PointEdgeWave.H:86
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
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