polyMeshFilter.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) 2012-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "polyMeshFilter.H"
29#include "polyMesh.H"
30#include "fvMesh.H"
31#include "unitConversion.H"
32#include "edgeCollapser.H"
33#include "syncTools.H"
34#include "polyTopoChange.H"
35#include "globalIndex.H"
36#include "bitSet.H"
37#include "pointSet.H"
38#include "faceSet.H"
39#include "cellSet.H"
40
41// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
42
43namespace Foam
44{
46}
47
48
49void Foam::polyMeshFilter::updateSets(const mapPolyMesh& map)
50{
51 updateSets<pointSet>(map);
52 updateSets<faceSet>(map);
53 updateSets<cellSet>(map);
54}
55
56
57void Foam::polyMeshFilter::copySets
58(
59 const polyMesh& oldMesh,
60 const polyMesh& newMesh
61)
62{
63 copySets<pointSet>(oldMesh, newMesh);
64 copySets<faceSet>(oldMesh, newMesh);
65 copySets<cellSet>(oldMesh, newMesh);
66}
67
68
70{
71 polyTopoChange originalMeshToNewMesh(mesh);
72
73 autoPtr<fvMesh> meshCopy;
74 autoPtr<mapPolyMesh> mapPtr = originalMeshToNewMesh.makeMesh
75 (
76 meshCopy,
78 (
79 mesh.name(),
80 mesh.polyMesh::instance(),
81 mesh.time(),
82 IOobject::READ_IF_PRESENT, // read fv* if present
84 false
85 ),
86 mesh,
87 true // parallel sync
88 );
89
90 const mapPolyMesh& map = mapPtr();
91
92 // Update fields
93 meshCopy().updateMesh(map);
94 if (map.hasMotionPoints())
95 {
96 meshCopy().movePoints(map.preMotionPoints());
97 }
98
99 copySets(mesh, meshCopy());
100
101 return meshCopy;
102}
103
104
105// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
106
107Foam::label Foam::polyMeshFilter::filterFacesLoop(const label nOriginalBadFaces)
108{
109 label nBadFaces = labelMax;
110 label nOuterIterations = 0;
111
112 // Maintain the number of times a point has been part of a bad face
113 labelList pointErrorCount(mesh_.nPoints(), Zero);
114
115 bitSet newErrorPoint(mesh_.nPoints());
117 (
118 mesh_,
119 meshQualityCoeffDict(),
120 newErrorPoint
121 );
122
123 bool newBadFaces = true;
124
125 // Main loop
126 // ~~~~~~~~~
127 // It tries and do some collapses, checks the resulting mesh and
128 // 'freezes' some edges (by marking in minEdgeLen) and tries again.
129 // This will iterate ultimately to the situation where every edge is
130 // frozen and nothing gets collapsed.
131 while
132 (
133 nOuterIterations < maxIterations()
134 //&& nBadFaces > nOriginalBadFaces
135 && newBadFaces
136 )
137 {
138 Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
139 << endl;
140
141 printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
142 printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
143
144 // Reset the new mesh to the old mesh
145 newMeshPtr_ = copyMesh(mesh_);
146 fvMesh& newMesh = newMeshPtr_();
147
148 scalarField newMeshFaceFilterFactor = faceFilterFactor_;
149 pointPriority_.reset(new labelList(originalPointPriority_));
150
151 labelList origToCurrentPointMap(identity(newMesh.nPoints()));
152
153 {
154 label nInnerIterations = 0;
155 label nPrevLocalCollapse = labelMax;
156
158
159 while (true)
160 {
161 Info<< nl << indent << "Inner iteration = "
162 << nInnerIterations++ << nl << incrIndent << endl;
163
164 label nLocalCollapse = filterFaces
165 (
166 newMesh,
167 newMeshFaceFilterFactor,
168 origToCurrentPointMap
169 );
171
172 if
173 (
174 nLocalCollapse >= nPrevLocalCollapse
175 || nLocalCollapse == 0
176 )
177 {
179 break;
180 }
181 else
182 {
183 nPrevLocalCollapse = nLocalCollapse;
184 }
185 }
186 }
187
188 scalarField newMeshMinEdgeLen = minEdgeLen_;
189
190 {
191 label nInnerIterations = 0;
192 label nPrevLocalCollapse = labelMax;
193
195
196 while (true)
197 {
198 Info<< nl << indent << "Inner iteration = "
199 << nInnerIterations++ << nl << incrIndent << endl;
200
201 label nLocalCollapse = filterEdges
202 (
203 newMesh,
204 newMeshMinEdgeLen,
205 origToCurrentPointMap
206 );
208
209 if
210 (
211 nLocalCollapse >= nPrevLocalCollapse
212 || nLocalCollapse == 0
213 )
214 {
216 break;
217 }
218 else
219 {
220 nPrevLocalCollapse = nLocalCollapse;
221 }
222 }
223 }
224
225
226 // Mesh check
227 // ~~~~~~~~~~~~~~~~~~
228 // Do not allow collapses in regions of error.
229 // Updates minEdgeLen, nRelaxedEdges
230
231 if (controlMeshQuality())
232 {
233 bitSet isErrorPoint(newMesh.nPoints());
235 (
236 newMesh,
237 meshQualityCoeffDict(),
238 isErrorPoint
239 );
240
241 Info<< nl << " Number of bad faces : " << nBadFaces << nl
242 << " Number of marked points : "
243 << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
244 << endl;
245
246 updatePointErrorCount
247 (
248 isErrorPoint,
249 origToCurrentPointMap,
250 pointErrorCount
251 );
252
253 checkMeshEdgesAndRelaxEdges
254 (
255 newMesh,
256 origToCurrentPointMap,
257 isErrorPoint,
258 pointErrorCount
259 );
260
261 checkMeshFacesAndRelaxEdges
262 (
263 newMesh,
264 origToCurrentPointMap,
265 isErrorPoint,
266 pointErrorCount
267 );
268
269 newBadFaces = false;
270 forAll(mesh_.points(), pI)
271 {
272 if
273 (
274 origToCurrentPointMap[pI] >= 0
275 && isErrorPoint[origToCurrentPointMap[pI]]
276 )
277 {
278 if (!newErrorPoint[pI])
279 {
280 newBadFaces = true;
281 break;
282 }
283 }
284 }
285
286 reduce(newBadFaces, orOp<bool>());
287 }
288 else
289 {
290 return -1;
291 }
292 }
293
294 return nBadFaces;
295}
296
297
299(
300 polyMesh& newMesh,
301 scalarField& newMeshFaceFilterFactor,
302 labelList& origToCurrentPointMap
303)
304{
305 // Per edge collapse status
306 bitSet collapseEdge(newMesh.nEdges());
307
308 Map<point> collapsePointToLocation(newMesh.nPoints());
309
310 edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
311
312 {
313 // Collapse faces
314 labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
315 (
316 newMeshFaceFilterFactor,
317 pointPriority_(),
319 collapsePointToLocation
320 );
321
322 label nCollapsed = 0;
323 forAll(nCollapsedPtEdge, collapseTypeI)
324 {
325 nCollapsed += nCollapsedPtEdge[collapseTypeI];
326 }
327
328 reduce(nCollapsed, sumOp<label>());
329
330 label nToPoint = returnReduce(nCollapsedPtEdge.first(), sumOp<label>());
331 label nToEdge = returnReduce(nCollapsedPtEdge.second(), sumOp<label>());
332
333 Info<< indent
334 << "Collapsing " << nCollapsed << " faces "
335 << "(to point = " << nToPoint << ", to edge = " << nToEdge << ")"
336 << endl;
337
338 if (nCollapsed == 0)
339 {
340 return 0;
341 }
342 }
343
344 // Merge edge collapses into consistent collapse-network.
345 // Make sure no cells get collapsed.
346 List<pointEdgeCollapse> allPointInfo;
347 const globalIndex globalPoints(newMesh.nPoints());
348
349 collapser.consistentCollapse
350 (
351 globalPoints,
352 pointPriority_(),
353 collapsePointToLocation,
355 allPointInfo
356 );
357
358 label nLocalCollapse = collapseEdge.count();
359
360 reduce(nLocalCollapse, sumOp<label>());
361 Info<< nl << indent << "Collapsing " << nLocalCollapse
362 << " edges after synchronisation and PointEdgeWave" << endl;
363
364 if (nLocalCollapse == 0)
365 {
366 return 0;
367 }
368
369 {
370 // Apply collapses to current mesh
371 polyTopoChange newMeshMod(newMesh);
372
373 // Insert mesh refinement into polyTopoChange.
374 collapser.setRefinement(allPointInfo, newMeshMod);
375
376 Info<< indent << "Apply changes to the current mesh" << endl;
377
378 // Apply changes to current mesh
379 autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
380 (
381 newMesh,
382 false
383 );
384 const mapPolyMesh& newMap = newMapPtr();
385
386 // Update fields
387 newMesh.updateMesh(newMap);
388 if (newMap.hasMotionPoints())
389 {
390 newMesh.movePoints(newMap.preMotionPoints());
391 }
392 updateSets(newMap);
393
394 updatePointPriorities(newMesh, newMap.pointMap());
395
396 mapOldMeshFaceFieldToNewMesh
397 (
398 newMesh,
399 newMap.faceMap(),
400 newMeshFaceFilterFactor
401 );
402
403 updateOldToNewPointMap
404 (
405 newMap.reversePointMap(),
406 origToCurrentPointMap
407 );
408 }
409
410 return nLocalCollapse;
411}
412
413
415(
416 polyMesh& newMesh,
417 scalarField& newMeshMinEdgeLen,
418 labelList& origToCurrentPointMap
419)
420{
421 // Per edge collapse status
422 bitSet collapseEdge(newMesh.nEdges());
423
424 Map<point> collapsePointToLocation(newMesh.nPoints());
425
426 edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
427
428 // Work out which edges to collapse
429 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
430 // This is by looking at minEdgeLen (to avoid frozen edges)
431 // and marking in collapseEdge.
432 label nSmallCollapsed = collapser.markSmallEdges
433 (
434 newMeshMinEdgeLen,
435 pointPriority_(),
437 collapsePointToLocation
438 );
439
440 reduce(nSmallCollapsed, sumOp<label>());
441 Info<< indent << "Collapsing " << nSmallCollapsed
442 << " small edges" << endl;
443
444 // Merge inline edges
445 label nMerged = collapser.markMergeEdges
446 (
447 maxCos(),
448 pointPriority_(),
450 collapsePointToLocation
451 );
452
453 reduce(nMerged, sumOp<label>());
454 Info<< indent << "Collapsing " << nMerged << " in line edges"
455 << endl;
456
457 if (nMerged + nSmallCollapsed == 0)
458 {
459 return 0;
460 }
461
462 // Merge edge collapses into consistent collapse-network.
463 // Make sure no cells get collapsed.
464 List<pointEdgeCollapse> allPointInfo;
465 const globalIndex globalPoints(newMesh.nPoints());
466
467 collapser.consistentCollapse
468 (
469 globalPoints,
470 pointPriority_(),
471 collapsePointToLocation,
473 allPointInfo
474 );
475
476 label nLocalCollapse = collapseEdge.count();
477
478 reduce(nLocalCollapse, sumOp<label>());
479 Info<< nl << indent << "Collapsing " << nLocalCollapse
480 << " edges after synchronisation and PointEdgeWave" << endl;
481
482 if (nLocalCollapse == 0)
483 {
484 return 0;
485 }
486
487 // Apply collapses to current mesh
488 polyTopoChange newMeshMod(newMesh);
489
490 // Insert mesh refinement into polyTopoChange.
491 collapser.setRefinement(allPointInfo, newMeshMod);
492
493 Info<< indent << "Apply changes to the current mesh" << endl;
494
495 // Apply changes to current mesh
496 autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
497 (
498 newMesh,
499 false
500 );
501 const mapPolyMesh& newMap = newMapPtr();
502
503 // Update fields
504 newMesh.updateMesh(newMap);
505 if (newMap.hasMotionPoints())
506 {
507 newMesh.movePoints(newMap.preMotionPoints());
508 }
509 updateSets(newMap);
510
511 // Synchronise the factors
512 mapOldMeshEdgeFieldToNewMesh
513 (
514 newMesh,
515 newMap.pointMap(),
516 newMeshMinEdgeLen
517 );
518
519 updateOldToNewPointMap
520 (
521 newMap.reversePointMap(),
522 origToCurrentPointMap
523 );
524
525 updatePointPriorities(newMesh, newMap.pointMap());
526
527 return nLocalCollapse;
528}
529
530
531void Foam::polyMeshFilter::updatePointErrorCount
532(
533 const bitSet& isErrorPoint,
534 const labelList& oldToNewMesh,
535 labelList& pointErrorCount
536) const
537{
538 forAll(mesh_.points(), pI)
539 {
540 if (isErrorPoint[oldToNewMesh[pI]])
541 {
542 pointErrorCount[pI]++;
543 }
544 }
545}
546
547
548void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
549(
550 const polyMesh& newMesh,
551 const labelList& oldToNewMesh,
552 const bitSet& isErrorPoint,
553 const labelList& pointErrorCount
554)
555{
556 const edgeList& edges = mesh_.edges();
557
558 forAll(edges, edgeI)
559 {
560 const edge& e = edges[edgeI];
561 label newStart = oldToNewMesh[e[0]];
562 label newEnd = oldToNewMesh[e[1]];
563
564 if
565 (
566 pointErrorCount[e[0]] >= maxPointErrorCount()
567 || pointErrorCount[e[1]] >= maxPointErrorCount()
568 )
569 {
570 minEdgeLen_[edgeI] = -1;
571 }
572
573 if
574 (
575 (newStart >= 0 && isErrorPoint[newStart])
576 || (newEnd >= 0 && isErrorPoint[newEnd])
577 )
578 {
579 minEdgeLen_[edgeI] *= edgeReductionFactor();
580 }
581 }
582
583 syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), scalar(0));
584
585 for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
586 {
587 // Smooth minEdgeLen
588 forAll(mesh_.edges(), edgeI)
589 {
590 const edge& e = mesh_.edges()[edgeI];
591
592 scalar sumMinEdgeLen = 0;
593 label nEdges = 0;
594
595 forAll(e, pointi)
596 {
597 const labelList& pEdges = mesh_.pointEdges()[e[pointi]];
598
599 forAll(pEdges, pEdgeI)
600 {
601 const label pEdge = pEdges[pEdgeI];
602 sumMinEdgeLen += minEdgeLen_[pEdge];
603 nEdges++;
604 }
605 }
606
607 minEdgeLen_[edgeI] = min
608 (
609 minEdgeLen_[edgeI],
610 sumMinEdgeLen/nEdges
611 );
612 }
613
615 (
616 mesh_,
617 minEdgeLen_,
618 minEqOp<scalar>(),
619 scalar(0)
620 );
621 }
622}
623
624
625void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
626(
627 const polyMesh& newMesh,
628 const labelList& oldToNewMesh,
629 const bitSet& isErrorPoint,
630 const labelList& pointErrorCount
631)
632{
633 const faceList& faces = mesh_.faces();
634
635 forAll(faces, facei)
636 {
637 const face& f = faces[facei];
638
639 forAll(f, fpI)
640 {
641 const label ptIndex = oldToNewMesh[f[fpI]];
642
643 if (pointErrorCount[f[fpI]] >= maxPointErrorCount())
644 {
645 faceFilterFactor_[facei] = -1;
646 }
647
648 if (isErrorPoint[ptIndex])
649 {
650 faceFilterFactor_[facei] *= faceReductionFactor();
651
652 break;
653 }
654 }
655 }
656
657 syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
658
659 for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
660 {
661 // Smooth faceFilterFactor
662 forAll(faces, facei)
663 {
664 const labelList& fEdges = mesh_.faceEdges()[facei];
665
666 scalar sumFaceFilterFactors = 0;
667 label nFaces = 0;
668
669 // This is important: Only smooth around faces that share an
670 // edge with a bad face
671 bool skipFace = true;
672
673 forAll(fEdges, fEdgeI)
674 {
675 const labelList& eFaces = mesh_.edgeFaces()[fEdges[fEdgeI]];
676
677 forAll(eFaces, eFacei)
678 {
679 const label eFace = eFaces[eFacei];
680
681 const face& f = faces[eFace];
682
683 forAll(f, fpI)
684 {
685 const label ptIndex = oldToNewMesh[f[fpI]];
686
687 if (isErrorPoint[ptIndex])
688 {
689 skipFace = false;
690 break;
691 }
692 }
693
694 if (eFace != facei)
695 {
696 sumFaceFilterFactors += faceFilterFactor_[eFace];
697 nFaces++;
698 }
699 }
700 }
701
702 if (skipFace)
703 {
704 continue;
705 }
706
707 faceFilterFactor_[facei] = min
708 (
709 faceFilterFactor_[facei],
710 sumFaceFilterFactors/nFaces
711 );
712 }
713
714 // Face filter factor needs to be synchronised!
715 syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
716 }
717}
718
719
720void Foam::polyMeshFilter::updatePointPriorities
721(
722 const polyMesh& newMesh,
723 const labelList& pointMap
724)
725{
726 labelList newPointPriority(newMesh.nPoints(), labelMin);
727 const labelList& currPointPriority = pointPriority_();
728
729 forAll(newPointPriority, ptI)
730 {
731 const label newPointToOldPoint = pointMap[ptI];
732 const label origPointPriority = currPointPriority[newPointToOldPoint];
733
734 newPointPriority[ptI] = max(origPointPriority, newPointPriority[ptI]);
735 }
736
738 (
739 newMesh,
740 newPointPriority,
741 maxEqOp<label>(),
743 );
744
745 pointPriority_.reset(new labelList(newPointPriority));
746}
747
748
749void Foam::polyMeshFilter::printScalarFieldStats
750(
751 const string desc,
752 const scalarField& fld
753) const
754{
755 scalar sum = 0;
756 scalar validElements = 0;
757 scalar min = GREAT;
758 scalar max = -GREAT;
759
760 forAll(fld, i)
761 {
762 const scalar fldElement = fld[i];
763
764 if (fldElement >= 0)
765 {
766 sum += fldElement;
767
768 if (fldElement < min)
769 {
770 min = fldElement;
771 }
772
773 if (fldElement > max)
774 {
775 max = fldElement;
776 }
777
778 validElements++;
779 }
780 }
781
782 reduce(sum, sumOp<scalar>());
783 reduce(min, minOp<scalar>());
784 reduce(max, maxOp<scalar>());
785 reduce(validElements, sumOp<label>());
786 const label totFieldSize = returnReduce(fld.size(), sumOp<label>());
787
788 Info<< incrIndent << indent << desc
789 << ": min = " << min
790 << " av = " << sum/(validElements + SMALL)
791 << " max = " << max << nl
792 << indent
793 << " " << validElements << " / " << totFieldSize << " elements used"
794 << decrIndent << endl;
795}
796
797
798void Foam::polyMeshFilter::mapOldMeshEdgeFieldToNewMesh
799(
800 const polyMesh& newMesh,
801 const labelList& pointMap,
802 scalarField& newMeshMinEdgeLen
803) const
804{
805 scalarField tmp(newMesh.nEdges());
806
807 const edgeList& newEdges = newMesh.edges();
808
809 forAll(newEdges, newEdgeI)
810 {
811 const edge& newEdge = newEdges[newEdgeI];
812 const label pStart = newEdge.start();
813 const label pEnd = newEdge.end();
814
815 tmp[newEdgeI] = min
816 (
817 newMeshMinEdgeLen[pointMap[pStart]],
818 newMeshMinEdgeLen[pointMap[pEnd]]
819 );
820 }
821
822 newMeshMinEdgeLen.transfer(tmp);
823
825 (
826 newMesh,
827 newMeshMinEdgeLen,
828 maxEqOp<scalar>(),
829 scalar(0)
830 );
831}
832
833
834void Foam::polyMeshFilter::mapOldMeshFaceFieldToNewMesh
835(
836 const polyMesh& newMesh,
837 const labelList& faceMap,
838 scalarField& newMeshFaceFilterFactor
839) const
840{
841 scalarField tmp(newMesh.nFaces());
842
843 forAll(faceMap, newFacei)
844 {
845 const label oldFacei = faceMap[newFacei];
846
847 tmp[newFacei] = newMeshFaceFilterFactor[oldFacei];
848 }
849
850 newMeshFaceFilterFactor.transfer(tmp);
851
853 (
854 newMesh,
855 newMeshFaceFilterFactor,
856 maxEqOp<scalar>()
857 );
858}
859
860
861void Foam::polyMeshFilter::updateOldToNewPointMap
862(
863 const labelList& currToNew,
864 labelList& origToCurrentPointMap
865) const
866{
867 forAll(origToCurrentPointMap, origPointi)
868 {
869 label oldPointi = origToCurrentPointMap[origPointi];
870
871 if (oldPointi != -1)
872 {
873 label newPointi = currToNew[oldPointi];
874
875 if (newPointi >= 0)
876 {
877 origToCurrentPointMap[origPointi] = newPointi;
878 }
879 else if (newPointi == -1)
880 {
881 origToCurrentPointMap[origPointi] = -1;
882 }
883 else
884 {
885 origToCurrentPointMap[origPointi] = -newPointi-2;
886 }
887 }
888 }
889}
890
891
892// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
893
895:
897 (
899 (
901 (
902 "collapseDict",
903 mesh.time().system(),
904 mesh.time(),
905 IOobject::MUST_READ,
906 IOobject::NO_WRITE
907 )
908 )
909 ),
910 mesh_(mesh),
911 newMeshPtr_(),
912 originalPointPriority_(mesh.nPoints(), labelMin),
913 pointPriority_(),
914 minEdgeLen_(),
915 faceFilterFactor_()
916{
918}
919
920
922(
923 const fvMesh& mesh,
924 const labelList& pointPriority
925)
926:
928 (
930 (
932 (
933 "collapseDict",
934 mesh.time().system(),
935 mesh.time(),
936 IOobject::MUST_READ,
937 IOobject::NO_WRITE
938 )
939 )
940 ),
941 mesh_(mesh),
942 newMeshPtr_(),
943 originalPointPriority_(pointPriority),
944 pointPriority_(),
945 minEdgeLen_(),
946 faceFilterFactor_()
947{
949}
950
952(
953 const fvMesh& mesh,
954 const labelList& pointPriority,
955 const dictionary& dict
956)
957:
959 mesh_(mesh),
960 newMeshPtr_(),
961 originalPointPriority_(pointPriority),
962 pointPriority_(),
963 minEdgeLen_(),
964 faceFilterFactor_()
965{
967}
968
969
970// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
971
972Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
973{
974 minEdgeLen_.resize(mesh_.nEdges(), minLen());
975 faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
976
977 return filterFacesLoop(nOriginalBadFaces);
978}
979
980
981Foam::label Foam::polyMeshFilter::filter(const faceSet& fSet)
982{
983 minEdgeLen_.resize(mesh_.nEdges(), minLen());
984 faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
985
986 forAll(faceFilterFactor_, fI)
987 {
988 if (!fSet.found(fI))
989 {
990 faceFilterFactor_[fI] = -1;
991 }
992 }
993
994 return filterFacesLoop(0);
995}
996
997
999(
1000 const label nOriginalBadFaces
1001)
1002{
1003 label nBadFaces = labelMax/2;
1004 label nPreviousBadFaces = labelMax;
1005 label nOuterIterations = 0;
1006
1007 minEdgeLen_.resize(mesh_.nEdges(), minLen());
1008 faceFilterFactor_.resize(0);
1009
1010 labelList pointErrorCount(mesh_.nPoints(), Zero);
1011
1012 // Main loop
1013 // ~~~~~~~~~
1014 // It tries and do some collapses, checks the resulting mesh and
1015 // 'freezes' some edges (by marking in minEdgeLen) and tries again.
1016 // This will iterate ultimately to the situation where every edge is
1017 // frozen and nothing gets collapsed.
1018 while
1019 (
1020 nOuterIterations < maxIterations()
1021 && nBadFaces > nOriginalBadFaces
1022 && nBadFaces < nPreviousBadFaces
1023 )
1024 {
1025 Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
1026 << endl;
1027
1028 printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
1029
1030 nPreviousBadFaces = nBadFaces;
1031
1032 // Reset the new mesh to the old mesh
1033 newMeshPtr_ = copyMesh(mesh_);
1034 fvMesh& newMesh = newMeshPtr_();
1035
1036 scalarField newMeshMinEdgeLen = minEdgeLen_;
1037 pointPriority_.reset(new labelList(originalPointPriority_));
1038
1039 labelList origToCurrentPointMap(identity(newMesh.nPoints()));
1040
1041 Info<< incrIndent;
1042
1043 label nInnerIterations = 0;
1044 label nPrevLocalCollapse = labelMax;
1045
1046 while (true)
1047 {
1048 Info<< nl
1049 << indent << "Inner iteration = " << nInnerIterations++ << nl
1050 << incrIndent << endl;
1051
1052 label nLocalCollapse = filterEdges
1053 (
1054 newMesh,
1055 newMeshMinEdgeLen,
1056 origToCurrentPointMap
1057 );
1058
1059 Info<< decrIndent;
1060
1061 if
1062 (
1063 nLocalCollapse >= nPrevLocalCollapse
1064 || nLocalCollapse == 0
1065 )
1066 {
1067 Info<< decrIndent;
1068 break;
1069 }
1070 else
1071 {
1072 nPrevLocalCollapse = nLocalCollapse;
1073 }
1074 }
1075
1076 // Mesh check
1077 // ~~~~~~~~~~~~~~~~~~
1078 // Do not allow collapses in regions of error.
1079 // Updates minEdgeLen, nRelaxedEdges
1080
1081 if (controlMeshQuality())
1082 {
1083 bitSet isErrorPoint(newMesh.nPoints());
1085 (
1086 newMesh,
1087 meshQualityCoeffDict(),
1088 isErrorPoint
1089 );
1090
1091 Info<< nl << " Number of bad faces : " << nBadFaces << nl
1092 << " Number of marked points : "
1093 << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
1094 << endl;
1095
1096 updatePointErrorCount
1097 (
1098 isErrorPoint,
1099 origToCurrentPointMap,
1100 pointErrorCount
1101 );
1102
1103 checkMeshEdgesAndRelaxEdges
1104 (
1105 newMesh,
1106 origToCurrentPointMap,
1107 isErrorPoint,
1108 pointErrorCount
1109 );
1110 }
1111 else
1112 {
1113 return -1;
1114 }
1115 }
1116
1117 return nBadFaces;
1118}
1119
1120
1122{
1123 return newMeshPtr_;
1124}
1125
1126
1129{
1130 return pointPriority_;
1131}
1132
1133
1134// ************************************************************************* //
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:500
Switch filterFaces() const
Filter faces at output time.
Definition: cvControlsI.H:210
Switch filterEdges() const
Filter edges at output time.
Definition: cvControlsI.H:205
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
static label checkMeshQuality(const polyMesh &mesh, const dictionary &meshQualityDict, bitSet &isErrorPoint)
Check mesh and mark points on faces in error.
Definition: edgeCollapser.C:87
A list of face labels.
Definition: faceSet.H:54
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
const word & name() const
Return reference to name.
Definition: fvMesh.H:310
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:613
bool hasMotionPoints() const
Has valid preMotionPoints?
Definition: mapPolyMesh.H:619
Class to store the settings for the polyMeshFilter class.
void writeSettings(Ostream &os) const
Write the settings to a stream.
Remove the edges and faces of a polyMesh whilst satisfying the given mesh quality criteria.
const autoPtr< labelList > & pointPriority() const
Return the new pointPriority list.
static autoPtr< fvMesh > copyMesh(const fvMesh &mesh)
Return a copy of an fvMesh.
const autoPtr< fvMesh > & filteredMesh() const
Return reference to the filtered mesh. Does not check if the.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Direct mesh changes based on v1.3 polyTopoChange syntax.
autoPtr< mapPolyMesh > makeMesh(autoPtr< Type > &newMesh, const IOobject &io, const polyMesh &mesh, const labelUList &patchMap, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Create new mesh with old mesh patches. Additional dictionaries.
label nPoints() const noexcept
Number of mesh points.
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
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.
virtual bool found(const label id) const
Has the given index?
Definition: topoSet.C:508
#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.
dynamicFvMesh & mesh
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
label nPoints
Namespace for OpenFOAM.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
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
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
int system(const std::string &command, const bool bg=false)
Execute the specified command via the shell.
Definition: MSwindows.C:1158
List< label > labelList
A List of labels.
Definition: List.H:66
constexpr label labelMin
Definition: label.H:60
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
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 reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Functor wrapper of allow/deny lists of wordRe for filtering.
Definition: wordRes.H:174
Unit conversion functions.