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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2019-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "edgeCollapser.H"
30#include "polyMesh.H"
31#include "polyTopoChange.H"
32#include "globalMeshData.H"
33#include "syncTools.H"
34#include "PointEdgeWave.H"
35#include "globalIndex.H"
36#include "removePoints.H"
37#include "motionSmoother.H"
38#include "OFstream.H"
39
40// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41
42namespace Foam
43{
45}
46
47
48// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
49
51(
52 const polyMesh& mesh,
53 const dictionary& meshQualityDict
54)
55{
56 labelHashSet badFaces(mesh.nFaces()/100);
57 DynamicList<label> checkFaces(mesh.nFaces());
58
59 const vectorField& fAreas = mesh.faceAreas();
60
61 scalar faceAreaLimit = SMALL;
62
63 forAll(fAreas, facei)
64 {
65 if (mag(fAreas[facei]) > faceAreaLimit)
66 {
67 checkFaces.append(facei);
68 }
69 }
70
71 Info<< endl;
72
74 (
75 false,
76 mesh,
77 meshQualityDict,
78 checkFaces,
79 badFaces
80 );
81
82 return badFaces;
83}
84
85
87(
88 const polyMesh& mesh,
89 const dictionary& meshQualityDict,
90 bitSet& isErrorPoint
91)
92{
94 (
95 mesh,
96 meshQualityDict
97 );
98
99 label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
100
101 for (const label facei : badFaces)
102 {
103 const face& f = mesh.faces()[facei];
104
105 isErrorPoint.set(f);
106 }
107
109 (
110 mesh,
111 isErrorPoint,
113 0
114 );
115
116 return nBadFaces;
117}
118
119
120// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
121
122Foam::labelList Foam::edgeCollapser::edgesFromPoints
123(
124 const label& facei,
126) const
127{
128 labelList edgeLabels(pointLabels.size() - 1, -1);
129
130 const labelList& faceEdges = mesh_.faceEdges()[facei];
131 const edgeList& edges = mesh_.edges();
132
133 label count = 0;
134
135 forAll(faceEdges, eI)
136 {
137 const label edgeI = faceEdges[eI];
138 const edge& e = edges[edgeI];
139
140 label match = 0;
141
143 {
144 if (e[0] == pointLabels[pI])
145 {
146 match++;
147 }
148
149 if (e[1] == pointLabels[pI])
150 {
151 match++;
152 }
153 }
154
155 if (match == 2)
156 {
157 // Edge found
158 edgeLabels[count++] = edgeI;
159 }
160 }
161
162 if (count != edgeLabels.size())
163 {
164 edgeLabels.setSize(count);
165 }
166
167 return edgeLabels;
168}
169
170
171void Foam::edgeCollapser::collapseToEdge
172(
173 const label facei,
174 const pointField& pts,
175 const labelList& pointPriority,
176 const vector& collapseAxis,
177 const point& fC,
178 const labelList& facePtsNeg,
179 const labelList& facePtsPos,
180 const scalarList& dNeg,
181 const scalarList& dPos,
182 const scalar dShift,
183 bitSet& collapseEdge,
184 Map<point>& collapsePointToLocation
185) const
186{
187 // Negative half
188
189 Foam::point collapseToPtA(GREAT, GREAT, GREAT);
190 //collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
191
192 label maxPriority = labelMin;
193 DynamicList<label> maxPriorityPts(max(dNeg.size(), dPos.size()));
194
195 forAll(facePtsNeg, fPtI)
196 {
197 const label facePointi = facePtsNeg[fPtI];
198 const label facePtPriority = pointPriority[facePointi];
199
200 if (facePtPriority > maxPriority)
201 {
202 maxPriority = facePtPriority;
203 maxPriorityPts.clear();
204 maxPriorityPts.append(facePointi);
205 }
206 else if (facePtPriority == maxPriority)
207 {
208 maxPriorityPts.append(facePointi);
209 }
210 }
211
212 if (!maxPriorityPts.empty())
213 {
214 Foam::point averagePt(Zero);
215
216 forAll(maxPriorityPts, ptI)
217 {
218 averagePt += pts[maxPriorityPts[ptI]];
219 }
220
221 collapseToPtA = averagePt/maxPriorityPts.size();
222// collapseToPtA = pts[maxPriorityPts.first()];
223 }
224
225 maxPriority = labelMin;
226 maxPriorityPts.clear();
227
228 labelList faceEdgesNeg = edgesFromPoints(facei, facePtsNeg);
229
230 collapseEdge.set(faceEdgesNeg);
231
232 forAll(facePtsNeg, pI)
233 {
234 collapsePointToLocation.set(facePtsNeg[pI], collapseToPtA);
235 }
236
237
238 // Positive half
239 Foam::point collapseToPtB(GREAT, GREAT, GREAT);
240// = collapseAxis*(sum(dPos)/dPos.size() - dShift) + fC;
241
242 forAll(facePtsPos, fPtI)
243 {
244 const label facePointi = facePtsPos[fPtI];
245 const label facePtPriority = pointPriority[facePointi];
246
247 if (facePtPriority > maxPriority)
248 {
249 maxPriority = facePtPriority;
250 maxPriorityPts.clear();
251 maxPriorityPts.append(facePointi);
252 }
253 else if (facePtPriority == maxPriority)
254 {
255 maxPriorityPts.append(facePointi);
256 }
257 }
258
259 if (!maxPriorityPts.empty())
260 {
261 Foam::point averagePt(Zero);
262
263 forAll(maxPriorityPts, ptI)
264 {
265 averagePt += pts[maxPriorityPts[ptI]];
266 }
267
268 collapseToPtB = averagePt/maxPriorityPts.size();
269// collapseToPtB = pts[maxPriorityPts.first()];
270 }
271
272 labelList faceEdgesPos = edgesFromPoints(facei, facePtsPos);
273
274 collapseEdge.set(faceEdgesPos);
275
276 forAll(facePtsPos, pI)
277 {
278 collapsePointToLocation.set(facePtsPos[pI], collapseToPtB);
279 }
280}
281
282
283void Foam::edgeCollapser::collapseToPoint
284(
285 const label& facei,
286 const pointField& pts,
287 const labelList& pointPriority,
288 const point& fC,
289 const labelList& facePts,
290 bitSet& collapseEdge,
291 Map<point>& collapsePointToLocation
292) const
293{
294 const face& f = mesh_.faces()[facei];
295
296 Foam::point collapseToPt = fC;
297
298 label maxPriority = labelMin;
299 DynamicList<label> maxPriorityPts(f.size());
300
301 forAll(facePts, fPtI)
302 {
303 const label facePointi = facePts[fPtI];
304 const label facePtPriority = pointPriority[facePointi];
305
306 if (facePtPriority > maxPriority)
307 {
308 maxPriority = facePtPriority;
309 maxPriorityPts.clear();
310 maxPriorityPts.append(facePointi);
311 }
312 else if (facePtPriority == maxPriority)
313 {
314 maxPriorityPts.append(facePointi);
315 }
316 }
317
318 if (!maxPriorityPts.empty())
319 {
320 Foam::point averagePt(Zero);
321
322 forAll(maxPriorityPts, ptI)
323 {
324 averagePt += pts[maxPriorityPts[ptI]];
325 }
326
327 collapseToPt = averagePt/maxPriorityPts.size();
328
329// collapseToPt = pts[maxPriorityPts.first()];
330 }
331
332// DynamicList<label> faceBoundaryPts(f.size());
333// DynamicList<label> faceFeaturePts(f.size());
334//
335// forAll(facePts, fPtI)
336// {
337// if (pointPriority[facePts[fPtI]] == 1)
338// {
339// faceFeaturePts.append(facePts[fPtI]);
340// }
341// else if (pointPriority[facePts[fPtI]] == 0)
342// {
343// faceBoundaryPts.append(facePts[fPtI]);
344// }
345// }
346//
347// if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
348// {
349// if (!faceFeaturePts.empty())
350// {
351// collapseToPt = pts[faceFeaturePts.first()];
352// }
353// else if (faceBoundaryPts.size() == 2)
354// {
355// collapseToPt =
356// 0.5
357// *(
358// pts[faceBoundaryPts[0]]
359// + pts[faceBoundaryPts[1]]
360// );
361// }
362// else if (faceBoundaryPts.size() <= f.size())
363// {
364// face bFace(faceBoundaryPts);
365//
366// collapseToPt = bFace.centre(pts);
367// }
368// }
369
370 const labelList& faceEdges = mesh_.faceEdges()[facei];
371
372 collapseEdge.set(faceEdges);
373
374 forAll(f, pI)
375 {
376 collapsePointToLocation.set(f[pI], collapseToPt);
377 }
378}
379
380
381void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
382(
383 const face& f,
384 const point& fC,
385 vector& collapseAxis,
386 scalar& aspectRatio
387) const
388{
389 const pointField& pts = mesh_.points();
390
391 symmTensor J(symm(f.inertia(pts, fC)));
392
393 // Find the dominant collapse direction by finding the eigenvector
394 // that corresponds to the normal direction, discarding it. The
395 // eigenvector corresponding to the smaller of the two remaining
396 // eigenvalues is the dominant axis in a high aspect ratio face.
397
398 scalar magJ = mag(J);
399
400 scalar detJ = SMALL;
401
402 if (magJ > VSMALL)
403 {
404 // Normalise inertia tensor to remove problems with small values
405
406 J /= mag(J);
407 // J /= cmptMax(J);
408 // J /= max(eigenValues(J).x(), SMALL);
409
410 // Calculating determinant, including stabilisation for zero or
411 // small negative values
412
413 detJ = max(det(J), SMALL);
414 }
415
416 if (detJ < 1e-5)
417 {
418 // It is possible that all the points of a face are the same
419
420 collapseAxis = f.edge(f.longestEdge(pts)).unitVec(pts);
421
422 // Empirical correlation for high aspect ratio faces
423
424 aspectRatio = Foam::sqrt(0.35/detJ);
425 }
426 else
427 {
428 vector eVals(eigenValues(J));
429
430 if (mag(eVals.y() - eVals.x()) < 100*SMALL)
431 {
432 // First two eigenvalues are the same: i.e. a square face
433
434 // Cannot necessarily determine linearly independent
435 // eigenvectors, or any at all, use longest edge direction.
436
437 collapseAxis = f.edge(f.longestEdge(pts)).unitVec(pts);
438
439 aspectRatio = 1.0;
440 }
441 else
442 {
443 // The maximum eigenvalue (z()) must be the direction of the
444 // normal, as it has the greatest value. The minimum eigenvalue
445 // is the dominant collapse axis for high aspect ratio faces.
446
447 collapseAxis = eigenVectors(J, eVals).x();
448
449 // The inertia calculation describes the mass distribution as a
450 // function of distance squared to the axis, so the square root of
451 // the ratio of face-plane moments gives a good indication of the
452 // aspect ratio.
453
454 aspectRatio = Foam::sqrt(eVals.y()/max(eVals.x(), SMALL));
455 }
456 }
457}
458
459
460Foam::scalarField Foam::edgeCollapser::calcTargetFaceSizes() const
461{
462 scalarField targetFaceSizes(mesh_.nFaces(), -1);
463
464 const scalarField& V = mesh_.cellVolumes();
465 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
466
467 const labelList& cellOwner = mesh_.faceOwner();
468 const labelList& cellNeighbour = mesh_.faceNeighbour();
469
470 const label nBoundaryFaces = mesh_.nBoundaryFaces();
471
472 // Calculate face size from cell volumes for internal faces
473 for (label intFacei = 0; intFacei < mesh_.nInternalFaces(); ++intFacei)
474 {
475 const scalar cellOwnerVol = max(0.0, V[cellOwner[intFacei]]);
476 const scalar cellNeighbourVol = max(0.0, V[cellNeighbour[intFacei]]);
477
478 scalar targetFaceSizeA = Foam::cbrt(cellOwnerVol);
479 scalar targetFaceSizeB = Foam::cbrt(cellNeighbourVol);
480
481 targetFaceSizes[intFacei] = 0.5*(targetFaceSizeA + targetFaceSizeB);
482 }
483
484 scalarField neiCellVolumes(nBoundaryFaces, -1);
485
486 // Now do boundary faces
487 forAll(patches, patchi)
488 {
489 const polyPatch& patch = patches[patchi];
490
491 label bFacei = patch.start() - mesh_.nInternalFaces();
492
493 if (patch.coupled())
494 {
495 // Processor boundary face: Need to get the cell volume on the other
496 // processor
497 const labelUList& faceCells = patch.faceCells();
498
499 forAll(faceCells, facei)
500 {
501 neiCellVolumes[bFacei++] = max(0.0, V[faceCells[facei]]);
502 }
503 }
504 else
505 {
506 // Normal boundary face: Just use owner cell volume to calculate
507 // the target face size
508 forAll(patch, patchFacei)
509 {
510 const label extFacei = patchFacei + patch.start();
511 const scalar cellOwnerVol = max(0.0, V[cellOwner[extFacei]]);
512
513 targetFaceSizes[extFacei] = Foam::cbrt(cellOwnerVol);
514 }
515 }
516 }
517
518 syncTools::swapBoundaryFaceList(mesh_, neiCellVolumes);
519
520 forAll(patches, patchi)
521 {
522 const polyPatch& patch = patches[patchi];
523
524 label bFacei = patch.start() - mesh_.nInternalFaces();
525
526 if (patch.coupled())
527 {
528 forAll(patch, patchFacei)
529 {
530 const label localFacei = patchFacei + patch.start();
531 const scalar cellOwnerVol = max(0.0, V[cellOwner[localFacei]]);
532 const scalar cellNeighbourVol = neiCellVolumes[bFacei++];
533
534 scalar targetFaceSizeA = Foam::cbrt(cellOwnerVol);
535 scalar targetFaceSizeB = Foam::cbrt(cellNeighbourVol);
536
537 targetFaceSizes[localFacei]
538 = 0.5*(targetFaceSizeA + targetFaceSizeB);
539 }
540 }
541 }
542
543 // Returns a characteristic length, not an area
544 return targetFaceSizes;
545}
546
547
548Foam::edgeCollapser::collapseType Foam::edgeCollapser::collapseFace
549(
550 const labelList& pointPriority,
551 const face& f,
552 const label facei,
553 const scalar targetFaceSize,
554 bitSet& collapseEdge,
555 Map<point>& collapsePointToLocation,
556 const scalarField& faceFilterFactor
557) const
558{
559 const scalar collapseSizeLimitCoeff = faceFilterFactor[facei];
560
561 const pointField& pts = mesh_.points();
562
563 labelList facePts(f);
564
565 const Foam::point fC = f.centre(pts);
566
567 const scalar fA = f.mag(pts);
568
569 vector collapseAxis = Zero;
570 scalar aspectRatio = 1.0;
571
572 faceCollapseAxisAndAspectRatio(f, fC, collapseAxis, aspectRatio);
573
574 // The signed distance along the collapse axis passing through the
575 // face centre that each vertex projects to.
576
577 scalarField d(f.size());
578
579 forAll(f, fPtI)
580 {
581 const Foam::point& pt = pts[f[fPtI]];
582
583 d[fPtI] = (collapseAxis & (pt - fC));
584 }
585
586 // Sort the projected distances and the corresponding vertex
587 // indices along the collapse axis
588
589 labelList oldToNew(sortedOrder(d));
590
591 oldToNew = invert(oldToNew.size(), oldToNew);
592
593 inplaceReorder(oldToNew, d);
594
595 inplaceReorder(oldToNew, facePts);
596
597 // Shift the points so that they are relative to the centre of the
598 // collapse line.
599
600 scalar dShift = -0.5*(d.first() + d.last());
601
602 d += dShift;
603
604 // Form two lists, one for each half of the set of points
605 // projected along the collapse axis.
606
607 // Middle value, index of first entry in the second half
608 label middle = -1;
609
610 forAll(d, dI)
611 {
612 if (d[dI] > 0)
613 {
614 middle = dI;
615
616 break;
617 }
618 }
619
620 if (middle == -1)
621 {
622// SeriousErrorInFunction
623// << "middle == -1, " << f << " " << d
624// << endl;//abort(FatalError);
625
626 return noCollapse;
627 }
628
629 // Negative half
630 SubList<scalar> dNeg(d, middle, 0);
631 SubList<label> facePtsNeg(facePts, middle, 0);
632
633 // Positive half
634 SubList<scalar> dPos(d, d.size() - middle, middle);
635 SubList<label> facePtsPos(facePts, d.size() - middle, middle);
636
637 // Defining how close to the midpoint (M) of the projected
638 // vertices line a projected vertex (X) can be before making this
639 // an invalid edge collapse
640 //
641 // X---X-g----------------M----X-----------g----X--X
642 //
643 // Only allow a collapse if all projected vertices are outwith
644 // guardFraction (g) of the distance form the face centre to the
645 // furthest vertex in the considered direction
646
647 if (dNeg.size() == 0 || dPos.size() == 0)
648 {
650 << "All points on one side of face centre, not collapsing."
651 << endl;
652 }
653
654// Info<< "Face : " << f << nl
655// << " Collapse Axis: " << collapseAxis << nl
656// << " Aspect Ratio : " << aspectRatio << endl;
657
658 collapseType typeOfCollapse = noCollapse;
659
660 if (magSqr(collapseAxis) < VSMALL)
661 {
662 typeOfCollapse = toPoint;
663 }
664 else if (fA < aspectRatio*sqr(targetFaceSize*collapseSizeLimitCoeff))
665 {
666 if
667 (
668 allowEarlyCollapseToPoint_
669 && (d.last() - d.first())
670 < targetFaceSize
671 *allowEarlyCollapseCoeff_*maxCollapseFaceToPointSideLengthCoeff_
672 )
673 {
674 typeOfCollapse = toPoint;
675 }
676 else if
677 (
678 (dNeg.last() < guardFraction_*dNeg.first())
679 && (dPos.first() > guardFraction_*dPos.last())
680 )
681 {
682 typeOfCollapse = toEdge;
683 }
684 else if
685 (
686 (d.last() - d.first())
687 < targetFaceSize
688 *maxCollapseFaceToPointSideLengthCoeff_
689 )
690 {
691 // If the face can't be collapsed to an edge, and it has a
692 // small enough span, collapse it to a point.
693 typeOfCollapse = toPoint;
694 }
695 }
696
697 if (typeOfCollapse == toPoint)
698 {
699 collapseToPoint
700 (
701 facei,
702 pts,
703 pointPriority,
704 fC,
705 facePts,
707 collapsePointToLocation
708 );
709 }
710 else if (typeOfCollapse == toEdge)
711 {
712 collapseToEdge
713 (
714 facei,
715 pts,
716 pointPriority,
717 collapseAxis,
718 fC,
719 facePtsNeg,
720 facePtsPos,
721 dNeg,
722 dPos,
723 dShift,
725 collapsePointToLocation
726 );
727 }
728
729 return typeOfCollapse;
730}
731
732
733Foam::label Foam::edgeCollapser::edgeMaster
734(
735 const labelList& pointPriority,
736 const edge& e
737) const
738{
739 label masterPoint = -1;
740
741 const label e0 = e.start();
742 const label e1 = e.end();
743
744 const label e0Priority = pointPriority[e0];
745 const label e1Priority = pointPriority[e1];
746
747 if (e0Priority > e1Priority)
748 {
749 masterPoint = e0;
750 }
751 else if (e0Priority < e1Priority)
752 {
753 masterPoint = e1;
754 }
755 else if (e0Priority == e1Priority)
756 {
757 masterPoint = e0;
758 }
759
760// // Collapse edge to point with higher priority.
761// if (pointPriority[e0] >= 0)
762// {
763// if (pointPriority[e1] >= 0)
764// {
765// // Both points have high priority. Choose one to collapse to.
766// // Note: should look at feature edges/points!
767// masterPoint = e0;
768// }
769// else
770// {
771// masterPoint = e0;
772// }
773// }
774// else
775// {
776// if (pointPriority[e1] >= 0)
777// {
778// masterPoint = e1;
779// }
780// else
781// {
782// // None on boundary. Neither is a master.
783// return -1;
784// }
785// }
786
787 return masterPoint;
788}
789
790
791void Foam::edgeCollapser::checkBoundaryPointMergeEdges
792(
793 const label pointi,
794 const label otherPointi,
795 const labelList& pointPriority,
796 Map<point>& collapsePointToLocation
797) const
798{
799 const pointField& points = mesh_.points();
800
801 const label e0Priority = pointPriority[pointi];
802 const label e1Priority = pointPriority[otherPointi];
803
804 if (e0Priority > e1Priority)
805 {
806 collapsePointToLocation.set
807 (
808 otherPointi,
809 points[pointi]
810 );
811 }
812 else if (e0Priority < e1Priority)
813 {
814 collapsePointToLocation.set
815 (
816 pointi,
817 points[otherPointi]
818 );
819 }
820 else // e0Priority == e1Priority
821 {
822 collapsePointToLocation.set
823 (
824 pointi,
825 points[otherPointi]
826 );
827
828// Foam::point averagePt
829// (
830// 0.5*(points[otherPointi] + points[pointi])
831// );
832//
833// collapsePointToLocation.set(pointi, averagePt);
834// collapsePointToLocation.set(otherPointi, averagePt);
835 }
836}
837
838
839Foam::label Foam::edgeCollapser::breakStringsAtEdges
840(
841 const bitSet& markedEdges,
842 bitSet& collapseEdge,
843 List<pointEdgeCollapse>& allPointInfo
844) const
845{
846 const edgeList& edges = mesh_.edges();
847 const labelListList& pointEdges = mesh_.pointEdges();
848
849 label nUncollapsed = 0;
850
851 forAll(edges, eI)
852 {
853 if (markedEdges.test(eI))
854 {
855 const edge& e = edges[eI];
856
857 const label startCollapseIndex
858 = allPointInfo[e.start()].collapseIndex();
859
860 if (startCollapseIndex != -1 && startCollapseIndex != -2)
861 {
862 const label endCollapseIndex
863 = allPointInfo[e.end()].collapseIndex();
864
865 if
866 (
867 !collapseEdge.test(eI)
868 && startCollapseIndex == endCollapseIndex
869 )
870 {
871 const labelList& ptEdgesStart = pointEdges[e.start()];
872
873 forAll(ptEdgesStart, ptEdgeI)
874 {
875 const label edgeI = ptEdgesStart[ptEdgeI];
876
877 const label nbrPointi
878 = edges[edgeI].otherVertex(e.start());
879 const label nbrIndex
880 = allPointInfo[nbrPointi].collapseIndex();
881
882 if
883 (
884 nbrIndex == startCollapseIndex
885 && collapseEdge.test(edgeI)
886 )
887 {
888 collapseEdge.unset(edgeI);
889 nUncollapsed++;
890 }
891 }
892 }
893 }
894 }
895 }
896
897 return nUncollapsed;
898}
899
900
901void Foam::edgeCollapser::determineDuplicatePointsOnFace
902(
903 const face& f,
904 bitSet& markedPoints,
905 labelHashSet& uniqueCollapses,
906 labelHashSet& duplicateCollapses,
907 List<pointEdgeCollapse>& allPointInfo
908) const
909{
910 uniqueCollapses.clear();
911 duplicateCollapses.clear();
912
913 forAll(f, fpI)
914 {
915 label index = allPointInfo[f[fpI]].collapseIndex();
916
917 // Check for consecutive duplicate
918 if (index != allPointInfo[f.prevLabel(fpI)].collapseIndex())
919 {
920 if (!uniqueCollapses.insert(index))
921 {
922 // Failed inserting so duplicate
923 duplicateCollapses.insert(index);
924 }
925 }
926 }
927
928 // Now duplicateCollapses contains duplicate collapse indices.
929 // Convert to points.
930 forAll(f, fpI)
931 {
932 label index = allPointInfo[f[fpI]].collapseIndex();
933 if (duplicateCollapses.found(index))
934 {
935 markedPoints.set(f[fpI]);
936 }
937 }
938}
939
940
941Foam::label Foam::edgeCollapser::countEdgesOnFace
942(
943 const face& f,
944 List<pointEdgeCollapse>& allPointInfo
945) const
946{
947 label nEdges = 0;
948
949 forAll(f, fpI)
950 {
951 const label pointi = f[fpI];
952 const label newPointi = allPointInfo[pointi].collapseIndex();
953
954 if (newPointi == -2)
955 {
956 nEdges++;
957 }
958 else
959 {
960 const label prevPointi = f[f.fcIndex(fpI)];
961 const label prevNewPointi
962 = allPointInfo[prevPointi].collapseIndex();
963
964 if (newPointi != prevNewPointi)
965 {
966 nEdges++;
967 }
968 }
969 }
970
971 return nEdges;
972}
973
974
975bool Foam::edgeCollapser::isFaceCollapsed
976(
977 const face& f,
978 List<pointEdgeCollapse>& allPointInfo
979) const
980{
981 label nEdges = countEdgesOnFace(f, allPointInfo);
982
983 // Polygons must have 3 or more edges to be valid
984 if (nEdges < 3)
985 {
986 return true;
987 }
988
989 return false;
990}
991
992
993// Create consistent set of collapses.
994// collapseEdge : per edge:
995// -1 : do not collapse
996// 0 : collapse to start
997// 1 : collapse to end
998// Note: collapseEdge has to be parallel consistent (in orientation)
999Foam::label Foam::edgeCollapser::syncCollapse
1000(
1001 const globalIndex& globalPoints,
1002 const labelList& pointPriority,
1003 const bitSet& collapseEdge,
1004 const Map<point>& collapsePointToLocation,
1005 List<pointEdgeCollapse>& allPointInfo
1006) const
1007{
1008 const edgeList& edges = mesh_.edges();
1009
1010 label nCollapsed = 0;
1011
1012 DynamicList<label> initPoints(mesh_.nPoints());
1013 DynamicList<pointEdgeCollapse> initPointInfo(mesh_.nPoints());
1014
1015 allPointInfo.clear();
1016 allPointInfo.setSize(mesh_.nPoints());
1017
1018 // Initialise edges to no collapse
1019 List<pointEdgeCollapse> allEdgeInfo
1020 (
1021 mesh_.nEdges(),
1022 pointEdgeCollapse(Zero, -1, -1)
1023 );
1024
1025 // Mark selected edges for collapse
1026 forAll(edges, edgeI)
1027 {
1028 if (collapseEdge[edgeI])
1029 {
1030 const edge& e = edges[edgeI];
1031
1032 label masterPointi = e.start();
1033
1034 // Choose the point on the edge with the highest priority.
1035 if (pointPriority[e.end()] > pointPriority[e.start()])
1036 {
1037 masterPointi = e.end();
1038 }
1039
1040 label masterPointPriority = pointPriority[masterPointi];
1041
1042 label index = globalPoints.toGlobal(masterPointi);
1043
1044 if (!collapsePointToLocation.found(masterPointi))
1045 {
1046 const label otherVertex = e.otherVertex(masterPointi);
1047
1048 if (!collapsePointToLocation.found(otherVertex))
1049 {
1051 << masterPointi << " on edge " << edgeI << " " << e
1052 << " is not marked for collapse."
1053 << abort(FatalError);
1054 }
1055 else
1056 {
1057 masterPointi = otherVertex;
1058 masterPointPriority = pointPriority[masterPointi];
1059 index = globalPoints.toGlobal(masterPointi);
1060 }
1061 }
1062
1063 const point& collapsePoint = collapsePointToLocation[masterPointi];
1064
1065 const pointEdgeCollapse pec
1066 (
1067 collapsePoint,
1068 index,
1069 masterPointPriority
1070 );
1071
1072 // Mark as collapsible but with nonsense master so it gets
1073 // overwritten and starts an update wave
1074 allEdgeInfo[edgeI] = pointEdgeCollapse
1075 (
1076 collapsePoint,
1077 labelMax,
1078 labelMin
1079 );
1080
1081 initPointInfo.append(pec);
1082 initPoints.append(e.start());
1083
1084 initPointInfo.append(pec);
1085 initPoints.append(e.end());
1086
1087 nCollapsed++;
1088 }
1089 }
1090
1091 PointEdgeWave<pointEdgeCollapse> collapsePropagator
1092 (
1093 mesh_,
1094 initPoints,
1095 initPointInfo,
1096 allPointInfo,
1097 allEdgeInfo,
1098 mesh_.globalData().nTotalPoints() // Maximum number of iterations
1099 );
1100
1101 return nCollapsed;
1102}
1103
1104
1105void Foam::edgeCollapser::filterFace
1106(
1107 const Map<DynamicList<label>>& collapseStrings,
1108 const List<pointEdgeCollapse>& allPointInfo,
1109 face& f
1110) const
1111{
1112 label newFp = 0;
1113
1114 face oldFace = f;
1115
1116 forAll(f, fp)
1117 {
1118 label pointi = f[fp];
1119
1120 label collapseIndex = allPointInfo[pointi].collapseIndex();
1121
1122 // Do we have a local point for this index?
1123 if (collapseStrings.found(collapseIndex))
1124 {
1125 const label localPointi = collapseStrings[collapseIndex][0];
1126
1127 if (!SubList<label>(f, newFp).found(localPointi))
1128 {
1129 f[newFp++] = localPointi;
1130 }
1131 }
1132 else if (collapseIndex == -1)
1133 {
1135 << "Point " << pointi << " was not visited by PointEdgeWave"
1136 << endl;
1137 }
1138 else
1139 {
1140 f[newFp++] = pointi;
1141 }
1142 }
1143
1144
1145 // Check for pinched face. Tries to correct
1146 // - consecutive duplicate vertex. Removes duplicate vertex.
1147 // - duplicate vertex with one other vertex in between (spike).
1148 // Both of these should not really occur! and should be checked before
1149 // collapsing edges.
1150
1151 const label size = newFp;
1152
1153 newFp = 2;
1154
1155 for (label fp = 2; fp < size; fp++)
1156 {
1157 label fp1 = fp-1;
1158 label fp2 = fp-2;
1159
1160 label pointi = f[fp];
1161
1162 // Search for previous occurrence.
1163 const label index = SubList<label>(f, fp).find(pointi);
1164
1165 if (index == fp1)
1166 {
1168 << "Removing consecutive duplicate vertex in face "
1169 << f << endl;
1170 // Don't store current pointi
1171 }
1172 else if (index == fp2)
1173 {
1175 << "Removing non-consecutive duplicate vertex in face "
1176 << f << endl;
1177 // Don't store current pointi and remove previous
1178 newFp--;
1179 }
1180 else if (index != -1)
1181 {
1183 << "Pinched face " << f << endl;
1184 f[newFp++] = pointi;
1185 }
1186 else
1187 {
1188 f[newFp++] = pointi;
1189 }
1190 }
1191
1192 f.setSize(newFp);
1193}
1194
1195
1196// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1197
1199:
1200 mesh_(mesh),
1201 guardFraction_(0),
1202 maxCollapseFaceToPointSideLengthCoeff_(0),
1203 allowEarlyCollapseToPoint_(false),
1204 allowEarlyCollapseCoeff_(0)
1205{}
1206
1207
1209(
1210 const polyMesh& mesh,
1211 const dictionary& dict
1212)
1213:
1214 mesh_(mesh),
1215 guardFraction_
1216 (
1217 dict.getOrDefault<scalar>("guardFraction", 0)
1218 ),
1219 maxCollapseFaceToPointSideLengthCoeff_
1220 (
1221 dict.getOrDefault<scalar>("maxCollapseFaceToPointSideLengthCoeff", 0)
1222 ),
1223 allowEarlyCollapseToPoint_
1224 (
1225 dict.getOrDefault("allowEarlyCollapseToPoint", true)
1226 ),
1227 allowEarlyCollapseCoeff_
1228 (
1229 dict.getOrDefault<scalar>("allowEarlyCollapseCoeff", 0)
1230 )
1231{
1232 if (debug)
1233 {
1234 Info<< "Edge Collapser Settings:" << nl
1235 << " Guard Fraction = " << guardFraction_ << nl
1236 << " Max collapse face to point side length = "
1237 << maxCollapseFaceToPointSideLengthCoeff_ << nl
1238 << " " << (allowEarlyCollapseToPoint_ ? "Allow" : "Do not allow")
1239 << " early collapse to point" << nl
1240 << " Early collapse coeff = " << allowEarlyCollapseCoeff_
1241 << endl;
1242 }
1243}
1244
1245
1246// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1247
1249(
1250 const List<pointEdgeCollapse>& allPointInfo,
1251 polyTopoChange& meshMod
1252) const
1253{
1254 const cellList& cells = mesh_.cells();
1255 const labelList& faceOwner = mesh_.faceOwner();
1256 const labelList& faceNeighbour = mesh_.faceNeighbour();
1257 const labelListList& pointFaces = mesh_.pointFaces();
1258 const pointZoneMesh& pointZones = mesh_.pointZones();
1259
1260
1261
1262
1263// // Dump point collapses
1264// label count = 0;
1265// forAll(allPointInfo, ptI)
1266// {
1267// const pointEdgeCollapse& pec = allPointInfo[ptI];
1268//
1269// if (mesh_.points()[ptI] != pec.collapsePoint())
1270// {
1271// count++;
1272// }
1273// }
1274//
1275// OFstream str("collapses_" + name(count) + ".obj");
1276// // Dump point collapses
1277// forAll(allPointInfo, ptI)
1278// {
1279// const pointEdgeCollapse& pec = allPointInfo[ptI];
1280//
1281// if
1282// (
1283// mesh_.points()[ptI] != pec.collapsePoint()
1284// && pec.collapsePoint() != vector(GREAT, GREAT, GREAT)
1285// )
1286// {
1287// meshTools::writeOBJ
1288// (
1289// str,
1290// mesh_.points()[ptI],
1291// pec.collapsePoint()
1292// );
1293// }
1294// }
1295
1296
1297
1298 bool meshChanged = false;
1299
1300 bitSet removedPoints(mesh_.nPoints());
1301
1302 // Create strings of edges.
1303 // Map from collapseIndex(=global master point) to set of points
1304 Map<DynamicList<label>> collapseStrings;
1305 {
1306 // 1. Count elements per collapseIndex
1307 Map<label> nPerIndex(mesh_.nPoints()/10);
1308 forAll(allPointInfo, pointi)
1309 {
1310 label collapseIndex = allPointInfo[pointi].collapseIndex();
1311
1312 if (collapseIndex != -1 && collapseIndex != -2)
1313 {
1314 ++(nPerIndex(collapseIndex, 0));
1315 }
1316 }
1317
1318 // 2. Size
1319 collapseStrings.resize(2*nPerIndex.size());
1320 forAllConstIters(nPerIndex, iter)
1321 {
1322 collapseStrings.insert(iter.key(), DynamicList<label>(iter.val()));
1323 }
1324
1325 // 3. Fill
1326 forAll(allPointInfo, pointi)
1327 {
1328 const label collapseIndex = allPointInfo[pointi].collapseIndex();
1329
1330 if (collapseIndex != -1 && collapseIndex != -2)
1331 {
1332 collapseStrings[collapseIndex].append(pointi);
1333 }
1334 }
1335 }
1336
1337
1338
1339
1340// OFstream str2("collapseStrings_" + name(count) + ".obj");
1341// // Dump point collapses
1342// forAllConstIters(collapseStrings, iter)
1343// {
1344// const label masterPoint = iter.key();
1345// const DynamicList<label>& edgeCollapses = iter();
1346//
1347// forAll(edgeCollapses, eI)
1348// {
1349// meshTools::writeOBJ
1350// (
1351// str2,
1352// mesh_.points()[edgeCollapses[eI]],
1353// mesh_.points()[masterPoint]
1354// );
1355// }
1356// }
1357
1358
1359
1360 // Current faces (is also collapseStatus: f.size() < 3)
1361 faceList newFaces(mesh_.faces());
1362
1363 // Current cellCollapse status
1364 bitSet cellRemoved(mesh_.nCells(), false);
1365
1366 label nUnvisited = 0;
1367 label nUncollapsed = 0;
1368 label nCollapsed = 0;
1369
1370 forAll(allPointInfo, pI)
1371 {
1372 const pointEdgeCollapse& pec = allPointInfo[pI];
1373
1374 if (pec.collapseIndex() == -1)
1375 {
1376 nUnvisited++;
1377 }
1378 else if (pec.collapseIndex() == -2)
1379 {
1380 nUncollapsed++;
1381 }
1382 else
1383 {
1384 nCollapsed++;
1385 }
1386 }
1387
1388 label nPoints = allPointInfo.size();
1389
1391 reduce(nUnvisited, sumOp<label>());
1392 reduce(nUncollapsed, sumOp<label>());
1393 reduce(nCollapsed, sumOp<label>());
1394
1395 Info<< incrIndent;
1396 Info<< indent << "Number of points : " << nPoints << nl
1397 << indent << "Not visited : " << nUnvisited << nl
1398 << indent << "Not collapsed : " << nUncollapsed << nl
1399 << indent << "Collapsed : " << nCollapsed << nl
1400 << endl;
1401 Info<< decrIndent;
1402
1403 do
1404 {
1405 forAll(newFaces, facei)
1406 {
1407 filterFace(collapseStrings, allPointInfo, newFaces[facei]);
1408 }
1409
1410 // Check if faces to be collapsed cause cells to become collapsed.
1411 label nCellCollapsed = 0;
1412
1413 forAll(cells, celli)
1414 {
1415 if (!cellRemoved.test(celli))
1416 {
1417 const cell& cFaces = cells[celli];
1418
1419 label nFaces = cFaces.size();
1420
1421 forAll(cFaces, i)
1422 {
1423 label facei = cFaces[i];
1424
1425 if (newFaces[facei].size() < 3)
1426 {
1427 --nFaces;
1428
1429 if (nFaces < 4)
1430 {
1431 Pout<< "Cell:" << celli
1432 << " uses faces:" << cFaces
1433 << " of which too many are marked for removal:"
1434 << endl
1435 << " ";
1436
1437
1438 forAll(cFaces, j)
1439 {
1440 if (newFaces[cFaces[j]].size() < 3)
1441 {
1442 Pout<< ' '<< cFaces[j];
1443 }
1444 }
1445 Pout<< endl;
1446
1447 cellRemoved.set(celli);
1448
1449 // Collapse all edges of cell to nothing
1450// collapseEdges(cellEdges[celli]);
1451
1452 nCellCollapsed++;
1453
1454 break;
1455 }
1456 }
1457 }
1458 }
1459 }
1460
1461 reduce(nCellCollapsed, sumOp<label>());
1462 Info<< indent << "Collapsing " << nCellCollapsed << " cells" << endl;
1463
1464 if (nCellCollapsed == 0)
1465 {
1466 break;
1467 }
1468
1469 } while (true);
1470
1471
1472 // Keep track of faces that have been done already.
1473 bitSet doneFace(mesh_.nFaces(), false);
1474
1475 {
1476 // Mark points used.
1477 bitSet usedPoint(mesh_.nPoints(), false);
1478
1479 forAll(cellRemoved, celli)
1480 {
1481 if (cellRemoved.test(celli))
1482 {
1483 meshMod.removeCell(celli, -1);
1484 }
1485 }
1486
1487 // Remove faces
1488 forAll(newFaces, facei)
1489 {
1490 const face& f = newFaces[facei];
1491
1492 if (f.size() < 3)
1493 {
1494 meshMod.removeFace(facei, -1);
1495 meshChanged = true;
1496
1497 // Mark face as been done - i.e. ignore it later
1498 doneFace.set(facei);
1499 }
1500 else
1501 {
1502 // Kept face. Mark vertices
1503 usedPoint.set(f);
1504 }
1505 }
1506
1507 // Remove unused vertices that have not been marked for removal already
1508 forAll(usedPoint, pointi)
1509 {
1510 if (!usedPoint.test(pointi))
1511 {
1512 removedPoints.set(pointi);
1513 meshMod.removePoint(pointi, -1);
1514 meshChanged = true;
1515 }
1516 }
1517 }
1518
1519 // Modify the point location of the remaining points
1520 forAll(allPointInfo, pointi)
1521 {
1522 const label collapseIndex = allPointInfo[pointi].collapseIndex();
1523 const point& collapsePoint = allPointInfo[pointi].collapsePoint();
1524
1525 if
1526 (
1527 !removedPoints.test(pointi)
1528 && collapseIndex != -1
1529 && collapseIndex != -2
1530 )
1531 {
1532 meshMod.modifyPoint
1533 (
1534 pointi,
1535 collapsePoint,
1536 pointZones.whichZone(pointi),
1537 false
1538 );
1539 }
1540 }
1541
1542
1543 const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
1544 const faceZoneMesh& faceZones = mesh_.faceZones();
1545
1546 // Renumber faces that use points
1547 forAll(allPointInfo, pointi)
1548 {
1549 if (removedPoints.test(pointi))
1550 {
1551 const labelList& changedFaces = pointFaces[pointi];
1552
1553 forAll(changedFaces, changedFacei)
1554 {
1555 label facei = changedFaces[changedFacei];
1556
1557 if (doneFace.set(facei))
1558 {
1559 // On first visit...
1560
1561 // Get current zone info
1562 label zoneID = faceZones.whichZone(facei);
1563
1564 bool zoneFlip = false;
1565
1566 if (zoneID >= 0)
1567 {
1568 const faceZone& fZone = faceZones[zoneID];
1569
1570 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
1571 }
1572
1573 // Get current connectivity
1574 label own = faceOwner[facei];
1575 label nei = -1;
1576 label patchID = -1;
1577
1578 if (mesh_.isInternalFace(facei))
1579 {
1580 nei = faceNeighbour[facei];
1581 }
1582 else
1583 {
1584 patchID = boundaryMesh.whichPatch(facei);
1585 }
1586
1587 meshMod.modifyFace
1588 (
1589 newFaces[facei], // face
1590 facei, // facei to change
1591 own, // owner
1592 nei, // neighbour
1593 false, // flipFaceFlux
1594 patchID, // patch
1595 zoneID,
1596 zoneFlip
1597 );
1598
1599 meshChanged = true;
1600 }
1601 }
1602 }
1603 }
1604
1605 return meshChanged;
1606}
1607
1608
1610(
1612 const labelList& pointPriority,
1613 const Map<point>& collapsePointToLocation,
1615 List<pointEdgeCollapse>& allPointInfo,
1616 const bool allowCellCollapse
1617) const
1618{
1619 // Make sure we don't collapse cells
1620 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1621 const faceList faces = mesh_.faces();
1622 const edgeList& edges = mesh_.edges();
1623 const labelListList& faceEdges = mesh_.faceEdges();
1624 const labelListList& pointEdges = mesh_.pointEdges();
1625 const cellList& cells = mesh_.cells();
1626
1627 labelHashSet uniqueCollapses;
1628 labelHashSet duplicateCollapses;
1629
1630 while (true)
1631 {
1632 label nUncollapsed = 0;
1633
1635 (
1636 mesh_,
1639 0
1640 );
1641
1642 // Create consistent set of collapses
1643 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1644 // Note: requires collapseEdge to be synchronised.
1645 syncCollapse
1646 (
1648 pointPriority,
1650 collapsePointToLocation,
1651 allPointInfo
1652 );
1653
1654 // Get collapsed faces
1655
1656 bitSet isCollapsedFace(mesh_.nFaces());
1657 bitSet markedPoints(mesh_.nPoints());
1658
1659 forAll(faces, facei)
1660 {
1661 const face& f = faces[facei];
1662
1663 isCollapsedFace[facei] = isFaceCollapsed(f, allPointInfo);
1664
1665 if (!isCollapsedFace.test(facei))
1666 {
1667 determineDuplicatePointsOnFace
1668 (
1669 f,
1670 markedPoints,
1671 uniqueCollapses,
1672 duplicateCollapses,
1673 allPointInfo
1674 );
1675 }
1676 }
1677
1678 // Synchronise the marked points
1680 (
1681 mesh_,
1682 markedPoints,
1684 0
1685 );
1686
1687 // Mark all edges attached to the point for collapse
1688 forAll(markedPoints, pointi)
1689 {
1690 if (markedPoints.test(pointi))
1691 {
1692 const label index = allPointInfo[pointi].collapseIndex();
1693
1694 const labelList& ptEdges = pointEdges[pointi];
1695
1696 forAll(ptEdges, ptEdgeI)
1697 {
1698 const label edgeI = ptEdges[ptEdgeI];
1699 const label nbrPointi = edges[edgeI].otherVertex(pointi);
1700 const label nbrIndex
1701 = allPointInfo[nbrPointi].collapseIndex();
1702
1703 if (nbrIndex == index && collapseEdge.test(edgeI))
1704 {
1705 collapseEdge.unset(edgeI);
1706 nUncollapsed++;
1707 }
1708 }
1709 }
1710 }
1711
1712 bitSet markedEdges(mesh_.nEdges());
1713
1714 if (!allowCellCollapse)
1715 {
1716 // Check collapsed cells
1717 forAll(cells, celli)
1718 {
1719 const cell& cFaces = cells[celli];
1720
1721 label nFaces = cFaces.size();
1722
1723 forAll(cFaces, fI)
1724 {
1725 label facei = cFaces[fI];
1726
1727 if (isCollapsedFace.test(facei))
1728 {
1729 --nFaces;
1730 }
1731 }
1732
1733 if (nFaces < 4)
1734 {
1735 forAll(cFaces, fI)
1736 {
1737 label facei = cFaces[fI];
1738
1739 const labelList& fEdges = faceEdges[facei];
1740
1741 // Unmark this face for collapse
1742 forAll(fEdges, fEdgeI)
1743 {
1744 label edgeI = fEdges[fEdgeI];
1745
1746 if (collapseEdge.unset(edgeI))
1747 {
1748 ++nUncollapsed;
1749 }
1750
1751 markedEdges.set(edgeI);
1752 }
1753
1754 // Uncollapsed this face.
1755 isCollapsedFace.unset(facei);
1756 ++nFaces;
1757 }
1758 }
1759
1760 if (nFaces < 4)
1761 {
1763 << "Cell " << celli << " " << cFaces << nl
1764 << "is " << nFaces << ", "
1765 << "but cell collapse has been disabled."
1766 << abort(FatalError);
1767 }
1768 }
1769 }
1770
1772 (
1773 mesh_,
1774 markedEdges,
1776 0
1777 );
1778
1779 nUncollapsed += breakStringsAtEdges
1780 (
1781 markedEdges,
1783 allPointInfo
1784 );
1785
1786 reduce(nUncollapsed, sumOp<label>());
1787
1788 Info<< " Uncollapsed edges = " << nUncollapsed << " / "
1789 << returnReduce(mesh_.nEdges(), sumOp<label>()) << endl;
1790
1791 if (nUncollapsed == 0)
1792 {
1793 break;
1794 }
1795 }
1796}
1797
1798
1800(
1801 const scalarField& minEdgeLen,
1802 const labelList& pointPriority,
1804 Map<point>& collapsePointToLocation
1805) const
1806{
1807 // Work out which edges to collapse
1808 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1809
1810 const pointField& points = mesh_.points();
1811 const edgeList& edges = mesh_.edges();
1812
1813 label nCollapsed = 0;
1814
1815 forAll(edges, edgeI)
1816 {
1817 const edge& e = edges[edgeI];
1818
1819 if (!collapseEdge.test(edgeI))
1820 {
1821 if (e.mag(points) < minEdgeLen[edgeI])
1822 {
1823 collapseEdge.set(edgeI);
1824
1825 label masterPointi = edgeMaster(pointPriority, e);
1826
1827 if (masterPointi == -1)
1828 {
1829 const point average
1830 = 0.5*(points[e.start()] + points[e.end()]);
1831
1832 collapsePointToLocation.set(e.start(), average);
1833 }
1834 else
1835 {
1836 const point& collapsePt = points[masterPointi];
1837
1838 collapsePointToLocation.set(masterPointi, collapsePt);
1839 }
1840
1841 nCollapsed++;
1842 }
1843 }
1844 }
1845
1846 return nCollapsed;
1847}
1848
1849
1851(
1852 const scalar maxCos,
1853 const labelList& pointPriority,
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.set(e0);
1904
1905 checkBoundaryPointMergeEdges
1906 (
1907 pointi,
1908 edges[e0].otherVertex(pointi),
1909 pointPriority,
1910 collapsePointToLocation
1911 );
1912 }
1913 else
1914 {
1915 collapseEdge.set(e1);
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
1938(
1939 const scalarField& faceFilterFactor,
1940 const labelList& pointPriority,
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],
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
1999(
2000 const faceZone& fZone,
2001 const scalarField& faceFilterFactor,
2002 const labelList& pointPriority,
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],
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.set(eI);
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.test(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// ************************************************************************* //
bool found
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
bool set(const Key &key, const T &obj)
Copy assign a new entry, overwriting existing entries.
Definition: HashTableI.H:202
void resize(const label sz)
Resize the hash table for efficiency.
Definition: HashTable.C:601
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpace.H:176
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
bitSet & unset(const bitSet &other)
Definition: bitSetI.H:628
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Does polyTopoChanges to remove edges. Can remove faces due to edge collapse but can not remove cells ...
Definition: edgeCollapser.H:69
labelPair markSmallSliverFaces(const scalarField &faceFilterFactor, const labelList &pointPriority, bitSet &collapseEdge, Map< point > &collapsePointToLocation) const
Find small faces and sliver faces in the mesh and mark the.
void consistentCollapse(const globalIndex &globalPoints, const labelList &pointPriority, const Map< point > &collapsePointToLocation, bitSet &collapseEdge, List< pointEdgeCollapse > &allPointInfo, const bool allowCellCollapse=false) const
Ensure that the collapse is parallel consistent and update.
bool setRefinement(const List< pointEdgeCollapse > &allPointInfo, polyTopoChange &meshMod) const
Play commands into polyTopoChange to create mesh.
static label checkMeshQuality(const polyMesh &mesh, const dictionary &meshQualityDict, bitSet &isErrorPoint)
Check mesh and mark points on faces in error.
Definition: edgeCollapser.C:87
label markSmallEdges(const scalarField &minEdgeLen, const labelList &pointPriority, bitSet &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to collapse.
static labelHashSet checkBadFaces(const polyMesh &mesh, const dictionary &meshQualityDict)
Calls motionSmoother::checkMesh and returns a set of bad faces.
Definition: edgeCollapser.C:51
label markMergeEdges(const scalar maxCos, const labelList &pointPriority, bitSet &collapseEdge, Map< point > &collapsePointToLocation) const
Mark (in collapseEdge) any edges to merge.
labelPair markFaceZoneEdges(const faceZone &fZone, const scalarField &faceFilterFactor, const labelList &pointPriority, bitSet &collapseEdge, Map< point > &collapsePointToLocation) const
Marks edges in the faceZone indirectPatchFaces for collapse.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:372
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:103
void checkMesh() const
Debug: Check coupled mesh for correctness.
Definition: hexRef8.C:4544
Determines length of string of edges walked to point.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
Direct mesh changes based on v1.3 polyTopoChange syntax.
void removeCell(const label celli, const label mergeCelli)
Remove/merge cell.
void modifyPoint(const label pointi, const point &pt, const label zoneID, const bool inCell)
Modify coordinate.
void removePoint(const label pointi, const label mergePointi)
Remove/merge point.
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.
void removeFace(const label facei, const label mergeFacei)
Remove/merge face.
label nFaces() const noexcept
Number of mesh faces.
const vectorField & faceAreas() const
Removes selected points from mesh and updates faces using these points.
Definition: removePoints.H:61
label countPointUsage(const scalar minCos, boolList &pointCanBeDeleted) const
Mark in pointCanBeDeleted the points that can be deleted.
Definition: removePoints.C:148
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
static void syncEdgeList(const polyMesh &mesh, List< T > &edgeValues, const CombineOp &cop, const T &nullValue, const TransformOp &top, const FlipOp &fop)
Synchronize values on all mesh edges.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
label collapseEdge(triSurface &surf, const scalar minLen)
Keep collapsing all edges < minLen.
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const pointField & points
label nPoints
const cellShapeList & cells
const labelIOList & zoneID
#define WarningInFunction
Report a warning using Foam::Warning.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
const std::string patch
OpenFOAM patch number as a std::string.
Namespace for OpenFOAM.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointFrompoint toPoint(const Foam::point &p)
constexpr label labelMin
Definition: label.H:60
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Definition: symmTensor.H:59
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:64
vector point
Point is a vector.
Definition: point.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr label labelMax
Definition: label.H:61
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
error FatalError
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
dimensionedScalar cbrt(const dimensionedScalar &ds)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
labelList pointLabels(nPoints, -1)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278