meshRefinementBlock.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) 2018-2022 OpenCFD Ltd.
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 "meshRefinement.H"
29#include "fvMesh.H"
30#include "Time.H"
31#include "refinementSurfaces.H"
32#include "removeCells.H"
33#include "unitConversion.H"
34#include "bitSet.H"
35#include "volFields.H"
36
37// Leak path
38#include "shortestPathSet.H"
39#include "meshSearch.H"
40#include "topoDistanceData.H"
41#include "FaceCellWave.H"
42#include "removeCells.H"
43#include "regionSplit.H"
44
45#include "volFields.H"
46#include "wallPoints.H"
47#include "searchableSurfaces.H"
49
50#include "holeToFace.H"
53#include "OBJstream.H"
54#include "PatchTools.H"
55
56// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
57
58//Foam::label Foam::meshRefinement::markFakeGapRefinement
59//(
60// const scalar planarCos,
61//
62// const label nAllowRefine,
63// const labelList& neiLevel,
64// const pointField& neiCc,
65//
66// labelList& refineCell,
67// label& nRefine
68//) const
69//{
70// label oldNRefine = nRefine;
71//
72//
73// // Collect candidate faces (i.e. intersecting any surface and
74// // owner/neighbour not yet refined.
75// const labelList testFaces(getRefineCandidateFaces(refineCell));
76//
77// // Collect segments
78// pointField start(testFaces.size());
79// pointField end(testFaces.size());
80// labelList minLevel(testFaces.size());
81//
82// calcCellCellRays
83// (
84// neiCc,
85// neiLevel,
86// testFaces,
87// start,
88// end,
89// minLevel
90// );
91//
92//
93// // Re-use the gap shooting methods. This needs:
94// // - shell gapLevel : faked. Set to 0,labelMax
95// // - surface gapLevel : faked by overwriting
96//
97//
98// List<FixedList<label, 3>>& surfGapLevel = const_cast
99// <
100// List<FixedList<label, 3>>&
101// >(surfaces_.extendedGapLevel());
102//
103// List<volumeType>& surfGapMode = const_cast
104// <
105// List<volumeType>&
106// >(surfaces_.extendedGapMode());
107//
108// const List<FixedList<label, 3>> surfOldLevel(surfGapLevel);
109// const List<volumeType> surfOldMode(surfGapMode);
110//
111// // Set the extended gap levels
112// forAll(surfaces_.gapLevel(), regioni)
113// {
114// surfGapLevel[regioni] = FixedList<label, 3>
115// ({
116// 3,
117// -1,
118// surfaces_.gapLevel()[regioni]+1
119// });
120// }
121// surfGapMode = volumeType::MIXED;
122//
123//Pout<< "gapLevel was:" << surfOldLevel << endl;
124//Pout<< "gapLevel now:" << surfGapLevel << endl;
125//Pout<< "gapMode was:" << surfOldMode << endl;
126//Pout<< "gapMode now:" << surfGapMode << endl;
127//Pout<< "nRefine was:" << oldNRefine << endl;
128//
129//
130//
131// List<List<FixedList<label, 3>>>& shellGapLevel = const_cast
132// <
133// List<List<FixedList<label, 3>>>&
134// >(shells_.extendedGapLevel());
135//
136// List<List<volumeType>>& shellGapMode = const_cast
137// <
138// List<List<volumeType>>&
139// >(shells_.extendedGapMode());
140//
141// const List<List<FixedList<label, 3>>> shellOldLevel(shellGapLevel);
142// const List<List<volumeType>> shellOldMode(shellGapMode);
143//
144// // Set the extended gap levels
145// forAll(shellGapLevel, shelli)
146// {
147// shellGapLevel[shelli] = FixedList<label, 3>({3, -1, labelMax});
148// shellGapMode[shelli] = volumeType::MIXED;
149// }
150//Pout<< "shellLevel was:" << shellOldLevel << endl;
151//Pout<< "shellLevel now:" << shellGapLevel << endl;
152//
153// const label nAdditionalRefined = markSurfaceGapRefinement
154// (
155// planarCos,
156//
157// nAllowRefine,
158// neiLevel,
159// neiCc,
160//
161// refineCell,
162// nRefine
163// );
164//
165//Pout<< "nRefine now:" << nRefine << endl;
166//
167// // Restore
168// surfGapLevel = surfOldLevel;
169// surfGapMode = surfOldMode;
170// shellGapLevel = shellOldLevel;
171// shellGapMode = shellOldMode;
172//
173// return nAdditionalRefined;
174//}
175
176
178(
179 const labelList& cellLevel,
180 const labelList& neiLevel,
181 const labelList& refineCell,
182 bitSet& isOutsideFace
183) const
184{
185 // Get faces:
186 // - on outside of cell set
187 // - inbetween same cell level (i.e. quads)
188
189 isOutsideFace.setSize(mesh_.nFaces());
190 isOutsideFace = Zero;
191
192 forAll(mesh_.faceNeighbour(), facei)
193 {
194 label own = mesh_.faceOwner()[facei];
195 label nei = mesh_.faceNeighbour()[facei];
196 if
197 (
198 (cellLevel[own] == cellLevel[nei])
199 && (
200 (refineCell[own] != -1)
201 != (refineCell[nei] != -1)
202 )
203 )
204 {
205 isOutsideFace.set(facei);
206 }
207 }
208 {
209
210 const label nBnd = mesh_.nBoundaryFaces();
211
212 labelList neiRefineCell(nBnd);
213 syncTools::swapBoundaryCellList(mesh_, refineCell, neiRefineCell);
214 for (label bFacei = 0; bFacei < nBnd; ++bFacei)
215 {
216 label facei = mesh_.nInternalFaces()+bFacei;
217 label own = mesh_.faceOwner()[facei];
218
219 if
220 (
221 (cellLevel[own] == neiLevel[bFacei])
222 && (
223 (refineCell[own] != -1)
224 != (neiRefineCell[bFacei] != -1)
225 )
226 )
227 {
228 isOutsideFace.set(facei);
229 }
230 }
231 }
232}
233
234
236(
237 const bitSet& isOutsideFace,
238 const label celli
239) const
240{
241 const cell& cFaces = mesh_.cells()[celli];
242 const vectorField& faceAreas = mesh_.faceAreas();
243
244 Vector<bool> haveDirs(vector::uniform(false));
245
246 forAll(cFaces, cFacei)
247 {
248 const label facei = cFaces[cFacei];
249
250 if (isOutsideFace[facei])
251 {
252 const vector& n = faceAreas[facei];
253 scalar magSqrN = magSqr(n);
254
255 if (magSqrN > ROOTVSMALL)
256 {
257 for
258 (
259 direction dir = 0;
260 dir < pTraits<vector>::nComponents;
261 dir++
262 )
263 {
264 if (Foam::sqr(n[dir]) > 0.99*magSqrN)
265 {
266 haveDirs[dir] = true;
267 break;
268 }
269 }
270 }
271 }
272 }
273
274 label nDirs = 0;
275 forAll(haveDirs, dir)
276 {
277 if (haveDirs[dir])
278 {
279 nDirs++;
280 }
281 }
282 return nDirs;
283}
284
285
287(
288 const labelList& neiLevel,
289 const bitSet& isOutsideFace,
291 label& nRefine
292) const
293{
294 // Get cells with three or more outside faces
295 const cellList& cells = mesh_.cells();
296 forAll(cells, celli)
297 {
298 if (refineCell[celli] == -1)
299 {
300 if (countFaceDirs(isOutsideFace, celli) == 3)
301 {
302 // Mark cell with any value
303 refineCell[celli] = 0;
304 nRefine++;
305 }
306 }
307 }
308}
309
310
311//void Foam::meshRefinement::markMultiRegionCell
312//(
313// const label celli,
314// const FixedList<label, 3>& surface,
315//
316// Map<FixedList<label, 3>>& cellToRegions,
317// bitSet& isMultiRegion
318//) const
319//{
320// if (!isMultiRegion[celli])
321// {
322// Map<FixedList<label, 3>>::iterator fnd = cellToRegions.find(celli);
323// if (!fnd.found())
324// {
325// cellToRegions.insert(celli, surface);
326// }
327// else if (fnd() != surface)
328// {
329// // Found multiple intersections on cell
330// isMultiRegion.set(celli);
331// }
332// }
333//}
334
335
336//void Foam::meshRefinement::detectMultiRegionCells
337//(
338// const labelListList& faceZones,
339// const labelList& testFaces,
340//
341// const labelList& surface1,
342// const List<pointIndexHit>& hit1,
343// const labelList& region1,
344//
345// const labelList& surface2,
346// const List<pointIndexHit>& hit2,
347// const labelList& region2,
348//
349// bitSet& isMultiRegion
350//) const
351//{
352// isMultiRegion.clear();
353// isMultiRegion.setSize(mesh_.nCells());
354//
355// Map<FixedList<label, 3>> cellToRegions(testFaces.size());
356//
357// forAll(testFaces, i)
358// {
359// const pointIndexHit& info1 = hit1[i];
360// if (info1.hit())
361// {
362// const label facei = testFaces[i];
363// const labelList& fz1 = faceZones[surface1[i]];
364// const FixedList<label, 3> surfaceInfo1
365// ({
366// surface1[i],
367// region1[i],
368// (fz1.size() ? fz1[info1.index()] : region1[i])
369// });
370//
371// markMultiRegionCell
372// (
373// mesh_.faceOwner()[facei],
374// surfaceInfo1,
375// cellToRegions,
376// isMultiRegion
377// );
378// if (mesh_.isInternalFace(facei))
379// {
380// markMultiRegionCell
381// (
382// mesh_.faceNeighbour()[facei],
383// surfaceInfo1,
384// cellToRegions,
385// isMultiRegion
386// );
387// }
388//
389// const pointIndexHit& info2 = hit2[i];
390//
391// if (info2.hit() && info1 != info2)
392// {
393// const labelList& fz2 = faceZones[surface2[i]];
394// const FixedList<label, 3> surfaceInfo2
395// ({
396// surface2[i],
397// region2[i],
398// (fz2.size() ? fz2[info2.index()] : region2[i])
399// });
400//
401// markMultiRegionCell
402// (
403// mesh_.faceOwner()[facei],
404// surfaceInfo2,
405// cellToRegions,
406// isMultiRegion
407// );
408// if (mesh_.isInternalFace(facei))
409// {
410// markMultiRegionCell
411// (
412// mesh_.faceNeighbour()[facei],
413// surfaceInfo2,
414// cellToRegions,
415// isMultiRegion
416// );
417// }
418// }
419// }
420// }
421//
422//
423// if (debug&meshRefinement::MESH)
424// {
425// volScalarField multiCell
426// (
427// IOobject
428// (
429// "multiCell",
430// mesh_.time().timeName(),
431// mesh_,
432// IOobject::NO_READ,
433// IOobject::NO_WRITE,
434// false
435// ),
436// mesh_,
437// dimensionedScalar
438// (
439// "zero",
440// dimensionSet(0, 1, 0, 0, 0),
441// 0.0
442// )
443// );
444// forAll(isMultiRegion, celli)
445// {
446// if (isMultiRegion[celli])
447// {
448// multiCell[celli] = 1.0;
449// }
450// }
451//
452// Info<< "Writing all multi cells to " << multiCell.name() << endl;
453// multiCell.write();
454// }
455//}
456
457
458Foam::label Foam::meshRefinement::markProximityRefinementWave
459(
460 const scalar planarCos,
461 const labelList& blockedSurfaces,
462 const label nAllowRefine,
463 const labelList& neiLevel,
464 const pointField& neiCc,
465
467 label& nRefine
468) const
469{
470 labelListList faceZones(surfaces_.surfaces().size());
471 {
472 // If triSurface do additional zoning based on connectivity
473 for (const label surfi : blockedSurfaces)
474 {
475 const label geomi = surfaces_.surfaces()[surfi];
476 const searchableSurface& s = surfaces_.geometry()[geomi];
477 if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
478 {
479 const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
480 const labelListList& edFaces = surf.edgeFaces();
481 boolList isOpenEdge(edFaces.size(), false);
482 forAll(edFaces, i)
483 {
484 if (edFaces[i].size() == 1)
485 {
486 isOpenEdge[i] = true;
487 }
488 }
489
490 labelList faceZone;
491 const label nZones = surf.markZones(isOpenEdge, faceZone);
492 if (nZones > 1)
493 {
494 faceZones[surfi].transfer(faceZone);
495 }
496 }
497 }
498 }
499
500
501 // Re-work the blockLevel of the blockedSurfaces into a length scale
502 // and a number of cells to cross
503 List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
504 for (const label surfi : blockedSurfaces)
505 {
506 const label geomi = surfaces_.surfaces()[surfi];
507 const searchableSurface& s = surfaces_.geometry()[geomi];
508 const label nRegions = s.regions().size();
509 regionToBlockSize[surfi].setSize(nRegions);
510 for (label regioni = 0; regioni < nRegions; regioni++)
511 {
512 const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
513 const label bLevel = surfaces_.blockLevel()[globalRegioni];
514 regionToBlockSize[surfi][regioni] =
515 meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
516
517 //const label mLevel = surfaces_.maxLevel()[globalRegioni];
521 //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
522 //{
523 // const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
524 //}
525
526 //nIters = max(nIters, (2<<(mLevel-bLevel)));
527 }
528 }
529
530 // Clever limiting of the number of iterations (= max cells in the channel)
531 // since it has too many problematic issues, e.g. with volume refinement
532 // and the real check uses the proper distance anyway just disable.
533 const label nIters = mesh_.globalData().nTotalCells();
534
535
536 // Collect candidate faces (i.e. intersecting any surface and
537 // owner/neighbour not yet refined)
538 const labelList testFaces(getRefineCandidateFaces(refineCell));
539
540 // Collect segments
541 pointField start(testFaces.size());
542 pointField end(testFaces.size());
543 labelList minLevel(testFaces.size());
544
545 calcCellCellRays
546 (
547 neiCc,
548 neiLevel,
549 testFaces,
550 start,
551 end,
552 minLevel
553 );
554 // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
555 // might add to cell level). Unfortunately we only have minLevel,
556 // not maxLevel. Problem: what if volume refinement only in middle of
557 // channel? Say channel 1m wide with a 0.1m sphere of refinement
558 // Workaround: have dummy surface with e.g. maxLevel 100 to
559 // force nIters to be high enough.
560
561
562 // Test for all intersections (with surfaces of higher gap level than
563 // minLevel) and cache per cell the max surface level and the local normal
564 // on that surface.
565
566 labelList surface1;
567 List<pointIndexHit> hit1;
568 labelList region1;
569 vectorField normal1;
570
571 labelList surface2;
572 List<pointIndexHit> hit2;
573 labelList region2;
574 vectorField normal2;
575
576 surfaces_.findNearestIntersection
577 (
578 blockedSurfaces,
579 start,
580 end,
581
582 surface1,
583 hit1,
584 region1, // local region
585 normal1,
586
587 surface2,
588 hit2,
589 region2, // local region
590 normal2
591 );
592
593
594 // Detect cells that are using multiple surface regions
595 //bitSet isMultiRegion;
596 //detectMultiRegionCells
597 //(
598 // faceZones,
599 // testFaces,
600 //
601 // surface1,
602 // hit1,
603 // region1,
604 //
605 // surface2,
606 // hit2,
607 // region2,
608 //
609 // isMultiRegion
610 //);
611
612
613 label n = 0;
614 forAll(testFaces, i)
615 {
616 if (hit1[i].hit())
617 {
618 n++;
619 }
620 }
621
622 List<wallPoints> faceDist(n);
623 labelList changedFaces(n);
624 n = 0;
625
626 DynamicList<point> originLocation(2);
627 DynamicList<scalar> originDistSqr(2);
628 DynamicList<FixedList<label, 3>> originSurface(2);
629 //DynamicList<point> originNormal(2);
630
631
632 //- To avoid walking through surfaces we mark all faces that have been
633 // intersected. We can either mark only those faces intersecting
634 // blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
635 // any (refinement) surface (this includes e.g. faceZones). This is
636 // much easier since that information is already cached
637 // (meshRefinement::intersectedFaces())
638
639 //bitSet isBlockedFace(mesh_.nFaces());
640 forAll(testFaces, i)
641 {
642 if (hit1[i].hit())
643 {
644 const label facei = testFaces[i];
645 //isBlockedFace.set(facei);
646 const point& fc = mesh_.faceCentres()[facei];
647 const labelList& fz1 = faceZones[surface1[i]];
648
649 originLocation.clear();
650 originDistSqr.clear();
651 //originNormal.clear();
652 originSurface.clear();
653
654 originLocation.append(hit1[i].hitPoint());
655 originDistSqr.append(magSqr(fc-hit1[i].hitPoint()));
656 //originNormal.append(normal1[i]);
657 originSurface.append
658 (
659 FixedList<label, 3>
660 ({
661 surface1[i],
662 region1[i],
663 (fz1.size() ? fz1[hit1[i].index()] : region1[i])
664 })
665 );
666
667 if (hit2[i].hit() && hit1[i] != hit2[i])
668 {
669 const labelList& fz2 = faceZones[surface2[i]];
670 originLocation.append(hit2[i].hitPoint());
671 originDistSqr.append(magSqr(fc-hit2[i].hitPoint()));
672 //originNormal.append(normal2[i]);
673 originSurface.append
674 (
675 FixedList<label, 3>
676 ({
677 surface2[i],
678 region2[i],
679 (fz2.size() ? fz2[hit2[i].index()] : region2[i])
680 })
681 );
682 }
683
684 // Collect all seed data. Currently walking does not look at
685 // surface direction - if so pass in surface normal as well
686 faceDist[n] = wallPoints
687 (
688 originLocation, // origin
689 originDistSqr, // distance to origin
690 originSurface // surface+region+zone
691 //originNormal // normal at origin
692 );
693 changedFaces[n] = facei;
694 n++;
695 }
696 }
697
698
699 // Clear intersection info
700 surface1.clear();
701 hit1.clear();
702 region1.clear();
703 normal1.clear();
704 surface2.clear();
705 hit2.clear();
706 region2.clear();
707 normal2.clear();
708
709
710 List<wallPoints> allFaceInfo(mesh_.nFaces());
711 List<wallPoints> allCellInfo(mesh_.nCells());
712
713 // Any refinement surface (even a faceZone) should stop the gap walking.
714 // This is exactly the information which is cached in the surfaceIndex_
715 // field.
716 const bitSet isBlockedFace(intersectedFaces());
717
718 wallPoints::trackData td(isBlockedFace, regionToBlockSize);
719 FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
720 (
721 mesh_,
722 changedFaces,
723 faceDist,
724 allFaceInfo,
725 allCellInfo,
726 0, // max iterations
727 td
728 );
729 wallDistCalc.iterate(nIters);
730
731
732 if (debug&meshRefinement::MESH)
733 {
734 // Dump current nearest opposite surfaces
736 (
737 IOobject
738 (
739 "gapSize",
740 mesh_.time().timeName(),
741 mesh_,
744 false
745 ),
746 mesh_,
748 (
749 "zero",
750 dimLength, //dimensionSet(0, 1, 0, 0, 0),
751 -1
752 )
753 );
754
755 forAll(allCellInfo, celli)
756 {
757 if (allCellInfo[celli].valid(wallDistCalc.data()))
758 {
759 const point& cc = mesh_.cellCentres()[celli];
760 // Nearest surface points
761 const List<point>& origin = allCellInfo[celli].origin();
762
763 // Find 'opposite' pair with minimum distance
764 for (label i = 0; i < origin.size(); i++)
765 {
766 for (label j = i + 1; j < origin.size(); j++)
767 {
768 if (((cc-origin[i]) & (cc-origin[j])) < 0)
769 {
770 const scalar d(mag(origin[i]-origin[j]));
771 if (distance[celli] < 0)
772 {
773 distance[celli] = d;
774 }
775 else
776 {
777 distance[celli] = min(distance[celli], d);
778 }
779 }
780 }
781 }
782 }
783 }
784 distance.correctBoundaryConditions();
785
786 Info<< "Writing measured gap distance to "
787 << distance.name() << endl;
788 distance.write();
789 }
790
791
792
793 // Detect tight gaps:
794 // - cell is inbetween the two surfaces
795 // - two surfaces are planarish
796 // - two surfaces are not too far apart
797 // (number of walking iterations is a too-coarse measure)
798
799 scalarField smallGapDistance(mesh_.nCells(), 0.0);
800 label nMulti = 0;
801 label nSmallGap = 0;
802
803 //OBJstream str(mesh_.time().timePath()/"multiRegion.obj");
804
805
806 forAll(allCellInfo, celli)
807 {
808 if (allCellInfo[celli].valid(wallDistCalc.data()))
809 {
810 const point& cc = mesh_.cellCentres()[celli];
811
812 const List<point>& origin = allCellInfo[celli].origin();
813 const List<FixedList<label, 3>>& surface =
814 allCellInfo[celli].surface();
815
816 // Find pair with minimum distance
817 for (label i = 0; i < origin.size(); i++)
818 {
819 for (label j = i + 1; j < origin.size(); j++)
820 {
821 //if (isMultiRegion[celli])
822 //{
823 // // The intersection locations are too inaccurate
824 // // (since not proper nearest, just a cell-cell ray
825 // // intersection) so include always
826 //
827 // smallGapDistance[celli] =
828 // max(smallGapDistance[celli], maxDist);
829 //
830 //
831 // str.write(linePointRef(cc, origin[i]));
832 // str.write(linePointRef(cc, origin[j]));
833 //
834 // nMulti++;
835 //}
836 //else
837 if (((cc-origin[i]) & (cc-origin[j])) < 0)
838 {
839 const label surfi = surface[i][0];
840 const label regioni = surface[i][1];
841
842 const label surfj = surface[j][0];
843 const label regionj = surface[j][1];
844
845 const scalar maxSize = max
846 (
847 regionToBlockSize[surfi][regioni],
848 regionToBlockSize[surfj][regionj]
849 );
850
851 if
852 (
853 magSqr(origin[i]-origin[j])
854 < Foam::sqr(2*maxSize)
855 )
856 {
857 const scalar maxDist
858 (
859 max
860 (
861 mag(cc-origin[i]),
862 mag(cc-origin[j])
863 )
864 );
865
866 smallGapDistance[celli] =
867 max(smallGapDistance[celli], maxDist);
868 nSmallGap++;
869 }
870 }
871 }
872 }
873 }
874 }
875
876
877 if (debug)
878 {
879 Info<< "Marked for blocking due to intersecting multiple surfaces : "
880 << returnReduce(nMulti, sumOp<label>()) << " cells." << endl;
881 Info<< "Marked for blocking due to close opposite surfaces : "
882 << returnReduce(nSmallGap, sumOp<label>()) << " cells." << endl;
883 }
884
885 if (debug&meshRefinement::MESH)
886 {
888 (
889 IOobject
890 (
891 "smallGapDistance",
892 mesh_.time().timeName(),
893 mesh_,
896 false
897 ),
898 mesh_,
900 (
901 "zero",
902 dimensionSet(0, 1, 0, 0, 0),
903 0.0
904 )
905 );
906 distance.field() = smallGapDistance;
907 distance.correctBoundaryConditions();
908
909 Info<< "Writing all small-gap cells to "
910 << distance.name() << endl;
911 distance.write();
912 }
913
914
915 // Mark refinement
916 const label oldNRefine = nRefine;
917 forAll(smallGapDistance, celli)
918 {
919 if (smallGapDistance[celli] > SMALL)
920 {
921 if
922 (
923 !markForRefine
924 (
925 0, // mark level
926 nAllowRefine,
927 refineCell[celli],
928 nRefine
929 )
930 )
931 {
932 if (debug)
933 {
934 Pout<< "Stopped refining since reaching my cell"
935 << " limit of " << mesh_.nCells()+7*nRefine
936 << endl;
937 }
938 break;
939 }
940 }
941 }
942
943 if
944 (
945 returnReduce(nRefine, sumOp<label>())
946 > returnReduce(nAllowRefine, sumOp<label>())
947 )
948 {
949 Info<< "Reached refinement limit." << endl;
950 }
951
952 return returnReduce(nRefine-oldNRefine, sumOp<label>());
953}
954
955
957(
958 const scalar planarAngle,
959 const labelList& minSurfaceLevel,
960 const labelList& globalToMasterPatch,
961 const label growIter
962)
963{
964 // Swap neighbouring cell centres and cell level
965 labelList neiLevel(mesh_.nBoundaryFaces());
966 pointField neiCc(mesh_.nBoundaryFaces());
967 calcNeighbourData(neiLevel, neiCc);
968
969 labelList refineCell(mesh_.nCells(), -1);
970 label nRefine = 0;
971 //markProximityRefinement
972 //(
973 // Foam::cos(degToRad(planarAngle)),
974 //
975 // minSurfaceLevel, // surface min level
976 // labelList(minSurfaceLevel.size(), labelMax), // surfaceGapLevel
977 //
978 // labelMax/Pstream::nProcs(), //nAllowRefine,
979 // neiLevel,
980 // neiCc,
981 //
982 // refineCell,
983 // nRefine
984 //);
985
986
987 // Determine minimum blockLevel per surface
988 Map<label> surfToBlockLevel;
989
990 forAll(surfaces_.surfaces(), surfi)
991 {
992 const label geomi = surfaces_.surfaces()[surfi];
993 const searchableSurface& s = surfaces_.geometry()[geomi];
994 const label nRegions = s.regions().size();
995
996 label minBlockLevel = labelMax;
997 for (label regioni = 0; regioni < nRegions; regioni++)
998 {
999 const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
1000 minBlockLevel = min
1001 (
1002 minBlockLevel,
1003 surfaces_.blockLevel()[globalRegioni]
1004 );
1005 }
1006
1007 if (minBlockLevel < labelMax)
1008 {
1009 surfToBlockLevel.insert(surfi, minBlockLevel);
1010 }
1011 }
1012
1013
1014 markProximityRefinementWave
1015 (
1016 Foam::cos(degToRad(planarAngle)),
1017 surfToBlockLevel.sortedToc(),
1018
1019 labelMax/Pstream::nProcs(), //nAllowRefine,
1020 neiLevel,
1021 neiCc,
1022
1023 refineCell,
1024 nRefine
1025 );
1026
1027
1029 //markFakeGapRefinement
1030 //(
1031 // Foam::cos(degToRad(planarAngle)),
1032 //
1033 // labelMax/Pstream::nProcs(), //nAllowRefine,
1034 // neiLevel,
1035 // neiCc,
1036 //
1037 // refineCell,
1038 // nRefine
1039 //);
1040
1041
1042 Info<< "Marked for blocking due to close opposite surfaces : "
1043 << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1044
1045 // Remove outliers, i.e. cells with all points exposed
1046 if (growIter)
1047 {
1048 labelList oldRefineCell(refineCell);
1049
1050 // Pass1: extend the set to fill in gaps
1051 bitSet isOutsideFace;
1052 for (label iter = 0; iter < growIter; iter++)
1053 {
1054 // Get outside faces
1055 markOutsideFaces
1056 (
1057 meshCutter_.cellLevel(),
1058 neiLevel,
1059 refineCell,
1060 isOutsideFace
1061 );
1062 // Extend with cells with three outside faces
1063 growSet(neiLevel, isOutsideFace, refineCell, nRefine);
1064 }
1065
1066
1067 // Pass2: erode back to original set if pass1 didn't help
1068 for (label iter = 0; iter < growIter; iter++)
1069 {
1070 // Get outside faces. Ignore cell level.
1071 markOutsideFaces
1072 (
1073 labelList(mesh_.nCells(), 0),
1074 labelList(neiLevel.size(), 0),
1075 refineCell,
1076 isOutsideFace
1077 );
1078
1079 // Unmark cells with three or more outside faces
1080 for (label celli = 0; celli < mesh_.nCells(); celli++)
1081 {
1082 if (refineCell[celli] != -1 && oldRefineCell[celli] == -1)
1083 {
1084 if (countFaceDirs(isOutsideFace, celli) >= 3)
1085 {
1086 refineCell[celli] = -1;
1087 --nRefine;
1088 }
1089 }
1090 }
1091 }
1092
1093 Info<< "Marked for blocking after filtering : "
1094 << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1095 }
1096
1097
1098 // Determine patch for every mesh face
1099 const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
1100 labelList unnamedSurfaces(surfaceZonesInfo::getUnnamedSurfaces(surfZones));
1101 const label defaultRegion(surfaces_.globalRegion(unnamedSurfaces[0], 0));
1102
1103 const labelList nearestRegion
1104 (
1105 nearestIntersection
1106 (
1107 unnamedSurfaces,
1108 defaultRegion
1109 )
1110 );
1111
1112 // Pack
1113 labelList cellsToRemove(nRefine);
1114 nRefine = 0;
1115
1116 forAll(refineCell, cellI)
1117 {
1118 if (refineCell[cellI] != -1)
1119 {
1120 cellsToRemove[nRefine++] = cellI;
1121 }
1122 }
1123
1124 // Remove cells
1125 removeCells cellRemover(mesh_);
1126 labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1127
1128 labelList exposedPatches(exposedFaces.size());
1129 forAll(exposedFaces, i)
1130 {
1131 label facei = exposedFaces[i];
1132 exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1133 }
1134
1135 return doRemoveCells
1136 (
1137 cellsToRemove,
1138 exposedFaces,
1139 exposedPatches,
1140 cellRemover
1141 );
1142}
1143
1144
1146(
1147 const labelList& selectedSurfaces,
1148 boolList& isBlockedFace
1149) const
1150{
1151 // Like meshRefinement::selectSeparatedCoupledFaces. tbd: convert to bitSet
1152
1153 // Check that no connection between inside and outside points
1154 isBlockedFace.setSize(mesh_.nFaces(), false);
1155
1156 // Block off separated couples.
1157 selectSeparatedCoupledFaces(isBlockedFace);
1158
1159 // Block off intersections with selected surfaces
1160
1161 // Mark per face (for efficiency)
1162 boolList isSelectedSurf(surfaces_.surfaces().size(), false);
1163 UIndirectList<bool>(isSelectedSurf, selectedSurfaces) = true;
1164
1165 forAll(surfaceIndex_, facei)
1166 {
1167 const label surfi = surfaceIndex_[facei];
1168 if (surfi != -1 && isSelectedSurf[surfi])
1169 {
1170 isBlockedFace[facei] = true;
1171 }
1172 }
1173}
1174
1175
1176//Foam::labelList Foam::meshRefinement::detectLeakCells
1177//(
1178// const boolList& isBlockedFace,
1179// const labelList& leakFaces,
1180// const labelList& seedCells
1181//) const
1182//{
1183// int dummyTrackData = 0;
1184// List<topoDistanceData<label>> allFaceInfo(mesh_.nFaces());
1185// List<topoDistanceData<label>> allCellInfo(mesh_.nCells());
1186//
1187// // Block faces
1188// forAll(isBlockedFace, facei)
1189// {
1190// if (isBlockedFace[facei])
1191// {
1192// allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1193// }
1194// }
1195// for (const label facei : leakFaces)
1196// {
1197// allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1198// }
1199//
1200// // Walk from inside cell
1201// DynamicList<topoDistanceData<label>> faceDist;
1202// DynamicList<label> cFaces1;
1203// for (const label celli : seedCells)
1204// {
1205// if (celli != -1)
1206// {
1207// const labelList& cFaces = mesh_.cells()[celli];
1208// faceDist.reserve(cFaces.size());
1209// cFaces1.reserve(cFaces.size());
1210//
1211// for (const label facei : cFaces)
1212// {
1213// if (!allFaceInfo[facei].valid(dummyTrackData))
1214// {
1215// cFaces1.append(facei);
1216// faceDist.append(topoDistanceData<label>(0, 123));
1217// }
1218// }
1219// }
1220// }
1221//
1222// // Walk through face-cell wave till all cells are reached
1223// FaceCellWave<topoDistanceData<label>> wallDistCalc
1224// (
1225// mesh_,
1226// cFaces1,
1227// faceDist,
1228// allFaceInfo,
1229// allCellInfo,
1230// mesh_.globalData().nTotalCells()+1 // max iterations
1231// );
1232//
1233// label nRemove = 0;
1234// for (const label facei : leakFaces)
1235// {
1236// const label own = mesh_.faceOwner()[facei];
1237// if (!allCellInfo[own].valid(dummyTrackData))
1238// {
1239// nRemove++;
1240// }
1241// if (mesh_.isInternalFace(facei))
1242// {
1243// const label nei = mesh_.faceNeighbour()[facei];
1244// if (!allCellInfo[nei].valid(dummyTrackData))
1245// {
1246// nRemove++;
1247// }
1248// }
1249// }
1250//
1251// labelList cellsToRemove(nRemove);
1252// nRemove = 0;
1253// for (const label facei : leakFaces)
1254// {
1255// const label own = mesh_.faceOwner()[facei];
1256// if (!allCellInfo[own].valid(dummyTrackData))
1257// {
1258// cellsToRemove[nRemove++] = own;
1259// }
1260// if (mesh_.isInternalFace(facei))
1261// {
1262// const label nei = mesh_.faceNeighbour()[facei];
1263// if (!allCellInfo[nei].valid(dummyTrackData))
1264// {
1265// cellsToRemove[nRemove++] = nei;
1266// }
1267// }
1268// }
1269//
1270// if (debug)
1271// {
1272// volScalarField fld
1273// (
1274// IOobject
1275// (
1276// "cellsToKeep",
1277// mesh_.time().timeName(),
1278// mesh_,
1279// IOobject::NO_READ,
1280// IOobject::NO_WRITE
1281// ),
1282// mesh_,
1283// dimensionedScalar(dimless, Zero)
1284// );
1285// forAll(allCellInfo, celli)
1286// {
1287// if (allCellInfo[celli].valid(dummyTrackData))
1288// {
1289// fld[celli] = allCellInfo[celli].distance();
1290// }
1291// }
1292// forAll(fld.boundaryField(), patchi)
1293// {
1294// const polyPatch& pp = mesh_.boundaryMesh()[patchi];
1295// SubList<topoDistanceData<label>> p(pp.patchSlice(allFaceInfo));
1296// scalarField pfld
1297// (
1298// fld.boundaryField()[patchi].size(),
1299// Zero
1300// );
1301// forAll(pfld, i)
1302// {
1303// if (p[i].valid(dummyTrackData))
1304// {
1305// pfld[i] = 1.0*p[i].distance();
1306// }
1307// }
1308// fld.boundaryFieldRef()[patchi] == pfld;
1309// }
1310// //Note: do not swap cell values so do not do
1311// //fld.correctBoundaryConditions();
1312// Pout<< "Writing distance field for initial cells "
1313// << seedCells << " to " << fld.objectPath() << endl;
1314// fld.write();
1315// }
1316//
1317// return cellsToRemove;
1318//}
1319//
1320//
1321//Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLeakCells
1322//(
1323// const labelList& globalToMasterPatch,
1324// const labelList& globalToSlavePatch,
1325// const pointField& locationsInMesh,
1326// const wordList& zonesInMesh,
1327// const pointField& locationsOutsideMesh,
1328// const labelList& selectedSurfaces
1329//)
1330//{
1331// boolList isBlockedFace;
1332// selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1333//
1334// // Determine cell regions
1335// const regionSplit cellRegion(mesh_, isBlockedFace);
1336//
1337// // Detect locationsInMesh regions
1338// labelList insideCells(locationsInMesh.size(), -1);
1339// labelList insideRegions(locationsInMesh.size(), -1);
1340// forAll(locationsInMesh, i)
1341// {
1342// insideCells[i] = findCell
1343// (
1344// mesh_,
1345// mergeDistance_*vector::one, //perturbVec,
1346// locationsInMesh[i]
1347// );
1348// if (insideCells[i] != -1)
1349// {
1350// insideRegions[i] = cellRegion[insideCells[i]];
1351// }
1352// reduce(insideRegions[i], maxOp<label>());
1353//
1354// if (insideRegions[i] == -1)
1355// {
1356// // See if we can perturb a bit
1357// insideCells[i] = findCell
1358// (
1359// mesh_,
1360// mergeDistance_*vector::one, //perturbVec,
1361// locationsInMesh[i]+mergeDistance_*vector::one
1362// );
1363// if (insideCells[i] != -1)
1364// {
1365// insideRegions[i] = cellRegion[insideCells[i]];
1366// }
1367// reduce(insideRegions[i], maxOp<label>());
1368//
1369// if (insideRegions[i] == -1)
1370// {
1371// FatalErrorInFunction
1372// << "Cannot find locationInMesh " << locationsInMesh[i]
1373// << " on any processor" << exit(FatalError);
1374// }
1375// }
1376// }
1377//
1378//
1379// // Check that all the locations outside the
1380// // mesh do not conflict with those inside
1381//
1382// bool haveLeak = false;
1383// forAll(locationsOutsideMesh, i)
1384// {
1385// // Find the region containing the point
1386// label regioni = findRegion
1387// (
1388// mesh_,
1389// cellRegion,
1390// mergeDistance_*vector::one, //perturbVec,
1391// locationsOutsideMesh[i]
1392// );
1393//
1394// if (regioni != -1)
1395// {
1396// // Check for locationsOutsideMesh overlapping with inside ones
1397// if (insideRegions.find(regioni) != -1)
1398// {
1399// haveLeak = true;
1400// WarningInFunction
1401// << "Outside location " << locationsOutsideMesh[i]
1402// << " in region " << regioni
1403// << " is connected to one of the inside points "
1404// << locationsInMesh << endl;
1405// }
1406// }
1407// }
1408//
1409//
1410// autoPtr<mapPolyMesh> mapPtr;
1411// if (returnReduce(haveLeak, orOp<bool>()))
1412// {
1413// // Use shortestPathSet to provide a minimum set of faces needed
1414// // to close hole. Tbd: maybe directly use wave?
1415// meshSearch searchEngine(mesh_);
1416// shortestPathSet leakPath
1417// (
1418// "leakPath",
1419// mesh_,
1420// searchEngine,
1421// coordSet::coordFormatNames[coordSet::coordFormat::DISTANCE],
1422// true,
1423// 50, // tbd. Number of iterations. This is the maximum
1424// // number of faces in the leak hole
1425//
1426// //pbm.groupPatchIDs()["wall"], // patch to grow from
1427// meshedPatches(), // patch to grow from
1428//
1429// locationsInMesh,
1430// locationsOutsideMesh,
1431// isBlockedFace
1432// );
1433//
1434//
1435// // Use leak path to find minimum set of cells to delete
1436// const labelList cellsToRemove
1437// (
1438// detectLeakCells
1439// (
1440// isBlockedFace,
1441// leakPath.leakFaces(),
1442// insideCells
1443// )
1444// );
1445//
1446// // Re-do intersections to find nearest unnamed surface
1447// const label defaultRegion
1448// (
1449// surfaces().globalRegion
1450// (
1451// selectedSurfaces[0],
1452// 0
1453// )
1454// );
1455//
1456// const labelList nearestRegion
1457// (
1458// nearestIntersection
1459// (
1460// selectedSurfaces,
1461// defaultRegion
1462// )
1463// );
1464//
1465//
1466// // Remove cells
1467// removeCells cellRemover(mesh_);
1468// const labelList exposedFaces
1469// (
1470// cellRemover.getExposedFaces(cellsToRemove)
1471// );
1472//
1473// labelList exposedPatches(exposedFaces.size());
1474// forAll(exposedFaces, i)
1475// {
1476// label facei = exposedFaces[i];
1477// exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1478// }
1479//
1480// mapPtr = doRemoveCells
1481// (
1482// cellsToRemove,
1483// exposedFaces,
1484// exposedPatches,
1485// cellRemover
1486// );
1487//
1488//
1489// // Put the exposed faces into a special faceZone
1490// {
1491// // Add to "frozenFaces" zone
1492// faceZoneMesh& faceZones = mesh_.faceZones();
1493//
1494// // Get current frozen faces (if any)
1495// bitSet isFrozenFace(mesh_.nFaces());
1496// label zonei = faceZones.findZoneID("frozenFaces");
1497// if (zonei != -1)
1498// {
1499// const bitSet oldSet(mesh_.nFaces(), faceZones[zonei]);
1500// isFrozenFace.set(oldSet);
1501// }
1502//
1503// // Add newly exposed faces (if not yet in any faceZone!)
1504// const labelList exposed
1505// (
1506// renumber
1507// (
1508// mapPtr().reverseFaceMap(),
1509// exposedFaces
1510// )
1511// );
1512// bitSet isZonedFace(mesh_.nFaces(), faceZones.zoneMap().toc());
1513// for (const label facei : exposed)
1514// {
1515// if (!isZonedFace[facei])
1516// {
1517// isFrozenFace.set(facei);
1518// }
1519// }
1520//
1521// syncTools::syncFaceList
1522// (
1523// mesh_,
1524// isFrozenFace,
1525// orEqOp<unsigned int>(),
1526// 0u
1527// );
1528//
1529// // Add faceZone if non existing
1530// faceZones.clearAddressing();
1531// if (zonei == -1)
1532// {
1533// zonei = faceZones.size();
1534// faceZones.setSize(zonei+1);
1535// faceZones.set
1536// (
1537// zonei,
1538// new faceZone
1539// (
1540// "frozenFaces", // name
1541// labelList(0), // addressing
1542// boolList(0), // flip
1543// zonei, // index
1544// faceZones // faceZoneMesh
1545// )
1546// );
1547// }
1548//
1549// // Update faceZone with new contents
1550// labelList frozenFaces(isFrozenFace.toc());
1551// boolList frozenFlip(frozenFaces.size(), false);
1552//
1553// faceZones[zonei].resetAddressing
1554// (
1555// std::move(frozenFaces),
1556// std::move(frozenFlip)
1557// );
1558// }
1559//
1560//
1561// //// Put the exposed points into a special pointZone
1562// //if (false)
1563// //{
1564// // const labelList meshFaceIDs
1565// // (
1566// // renumber
1567// // (
1568// // mapPtr().reverseFaceMap(),
1569// // exposedFaces
1570// // )
1571// // );
1572// // const uindirectPrimitivePatch pp
1573// // (
1574// // UIndirectList<face>(mesh_.faces(), meshFaceIDs),
1575// // mesh_.points()
1576// // );
1577// //
1578// // // Count number of faces per edge
1579// // const labelListList& edgeFaces = pp.edgeFaces();
1580// // labelList nEdgeFaces(edgeFaces.size());
1581// // forAll(edgeFaces, edgei)
1582// // {
1583// // nEdgeFaces[edgei] = edgeFaces[edgei].size();
1584// // }
1585// //
1586// // // Sync across processor patches
1587// // if (Pstream::parRun())
1588// // {
1589// // const globalMeshData& globalData = mesh_.globalData();
1590// // const mapDistribute& map = globalData.globalEdgeSlavesMap();
1591// // const indirectPrimitivePatch& cpp =
1592// // globalData.coupledPatch();
1593// //
1594// // // Match pp edges to coupled edges
1595// // labelList patchEdges;
1596// // labelList coupledEdges;
1597// // PackedBoolList sameEdgeOrientation;
1598// // PatchTools::matchEdges
1599// // (
1600// // pp,
1601// // cpp,
1602// // patchEdges,
1603// // coupledEdges,
1604// // sameEdgeOrientation
1605// // );
1606// //
1607// //
1608// // // Convert patch-edge data into cpp-edge data
1609// // labelList coupledNEdgeFaces(map.constructSize(), Zero);
1610// // UIndirectList<label>(coupledNEdgeFaces, coupledEdges) =
1611// // UIndirectList<label>(nEdgeFaces, patchEdges);
1612// //
1613// // // Synchronise
1614// // globalData.syncData
1615// // (
1616// // coupledNEdgeFaces,
1617// // globalData.globalEdgeSlaves(),
1618// // globalData.globalEdgeTransformedSlaves(),
1619// // map,
1620// // plusEqOp<label>()
1621// // );
1622// //
1623// // // Convert back from cpp-edge to patch-edge
1624// // UIndirectList<label>(nEdgeFaces, patchEdges) =
1625// // UIndirectList<label>(coupledNEdgeFaces, coupledEdges);
1626// // }
1627// //
1628// // // Freeze all internal points
1629// // bitSet isFrozenPoint(mesh_.nPoints());
1630// // forAll(nEdgeFaces, edgei)
1631// // {
1632// // if (nEdgeFaces[edgei] != 1)
1633// // {
1634// // const edge& e = pp.edges()[edgei];
1635// // isFrozenPoint.set(pp.meshPoints()[e[0]]);
1636// // isFrozenPoint.set(pp.meshPoints()[e[1]]);
1637// // }
1638// // }
1639// // // Add to "frozenPoints" zone
1640// // pointZoneMesh& pointZones = mesh_.pointZones();
1641// // pointZones.clearAddressing();
1642// // label zonei = pointZones.findZoneID("frozenPoints");
1643// // if (zonei != -1)
1644// // {
1645// // const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
1646// // // Add to isFrozenPoint
1647// // isFrozenPoint.set(oldSet);
1648// // }
1649// //
1650// // syncTools::syncPointList
1651// // (
1652// // mesh_,
1653// // isFrozenPoint,
1654// // orEqOp<unsigned int>(),
1655// // 0u
1656// // );
1657// //
1658// // if (zonei == -1)
1659// // {
1660// // zonei = pointZones.size();
1661// // pointZones.setSize(zonei+1);
1662// // pointZones.set
1663// // (
1664// // zonei,
1665// // new pointZone
1666// // (
1667// // "frozenPoints", // name
1668// // isFrozenPoint.sortedToc(), // addressing
1669// // zonei, // index
1670// // pointZones // pointZoneMesh
1671// // )
1672// // );
1673// // }
1674// //}
1675//
1676//
1677// if (debug&meshRefinement::MESH)
1678// {
1679// const_cast<Time&>(mesh_.time())++;
1680//
1681// Pout<< "Writing current mesh to time "
1682// << timeName() << endl;
1683// write
1684// (
1685// meshRefinement::debugType(debug),
1686// meshRefinement::writeType
1687// (
1688// meshRefinement::writeLevel()
1689// | meshRefinement::WRITEMESH
1690// ),
1691// mesh_.time().path()/timeName()
1692// );
1693// Pout<< "Dumped mesh in = "
1694// << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
1695// }
1696// }
1697// return mapPtr;
1698//}
1699
1700
1701//Foam::autoPtr<Foam::mapDistribute> Foam::meshRefinement::nearestFace
1702//(
1703// const globalIndex& globalSeedFaces,
1704// const labelList& seedFaces,
1705// const labelList& closureFaces
1706//) const
1707//{
1708// // Takes a set of faces for which we have information (seedFaces,
1709// // globalSeedFaces - these are e.g. intersected faces) and walks out
1710// // across nonSeedFace. Returns for
1711// // every nonSeedFace the nearest seed face (in global indexing).
1712// // Used e.g. in hole closing. Assumes that seedFaces, closureFaces
1713// // are a small subset of the master
1714// // so are not large - it uses edge addressing on the closureFaces
1715//
1716//
1717// if (seedFaces.size() != globalSeedFaces.localSize())
1718// {
1719// FatalErrorInFunction << "problem : seedFaces:" << seedFaces.size()
1720// << " globalSeedFaces:" << globalSeedFaces.localSize()
1721// << exit(FatalError);
1722// }
1723//
1724// //// Mark mesh points that are used by any closureFaces. This is for
1725// //// initial filtering
1726// //bitSet isNonSeedPoint(mesh.nPoints());
1727// //for (const label facei : closureFaces)
1728// //{
1729// // isNonSeedPoint.set(mesh_.faces()[facei]);
1730// //}
1731// //syncTools::syncPointList
1732// //(
1733// // mesh_,
1734// // isNonSeedPoint,
1735// // orEqOp<unsigned int>(),
1736// // 0u
1737// //);
1738//
1739// // Make patch
1740// const uindirectPrimitivePatch pp
1741// (
1742// IndirectList<face>(mesh_.faces(), closureFaces),
1743// mesh_.points()
1744// );
1745// const edgeList& edges = pp.edges();
1746// const labelList& mp = pp.meshPoints();
1747// const label nBndEdges = pp.nEdges() - pp.nInternalEdges();
1748//
1749// // For all faces in seedFaces mark the edge with a face. No special
1750// // handling for multiple faces sharing the edge - first one wins
1751// EdgeMap<label> edgeMap(pp.nEdges());
1752// for (const label facei : seedFaces)
1753// {
1754// const label globalFacei = globalSeedFaces.toGlobal(facei);
1755// const face& f = mesh_.faces()[facei];
1756// forAll(f, fp)
1757// {
1758// label nextFp = f.fcIndex(fp);
1759// edgeMap.insert(edge(f[fp], f[nextFp]), globalFacei);
1760// }
1761// }
1762// syncTools::syncEdgeMap(mesh_, edgeMap, maxEqOp<label>());
1763//
1764//
1765//
1766// // Seed
1767// DynamicList<label> initialEdges(2*nBndEdges);
1768// DynamicList<edgeTopoDistanceData<label, uindirectPrimitivePatch>>
1769// initialEdgesInfo(2*nBndEdges);
1770// forAll(edges, edgei)
1771// {
1772// const edge& e = edges[edgei];
1773// const edge meshE = edge(mp[e[0]], mp[e[1]]);
1774//
1775// EdgeMap<label>::const_iterator iter = edgeMap.find(meshE);
1776// if (iter.found())
1777// {
1778// initialEdges.append(edgei);
1779// initialEdgesInfo.append
1780// (
1781// edgeTopoDistanceData<label, uindirectPrimitivePatch>
1782// (
1783// 0, // distance
1784// iter() // globalFacei
1785// )
1786// );
1787// }
1788// }
1789//
1790// // Data on all edges and faces
1791// List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allEdgeInfo
1792// (
1793// pp.nEdges()
1794// );
1795// List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allFaceInfo
1796// (
1797// pp.size()
1798// );
1799//
1800// // Walk
1801// PatchEdgeFaceWave
1802// <
1803// uindirectPrimitivePatch,
1804// edgeTopoDistanceData<label, uindirectPrimitivePatch>
1805// > calc
1806// (
1807// mesh_,
1808// pp,
1809// initialEdges,
1810// initialEdgesInfo,
1811// allEdgeInfo,
1812// allFaceInfo,
1813// returnReduce(pp.nEdges(), sumOp<label>())
1814// );
1815//
1816//
1817// // Per non-seed face the seed face
1818// labelList closureToSeed(pp.size(), -1);
1819// forAll(allFaceInfo, facei)
1820// {
1821// if (allFaceInfo[facei].valid(calc.data()))
1822// {
1823// closureToSeed[facei] = allFaceInfo[facei].data();
1824// }
1825// }
1826//
1827// List<Map<label>> compactMap;
1828// return autoPtr<mapDistribute>::New
1829// (
1830// globalSeedFaces,
1831// closureToSeed,
1832// compactMap
1833// );
1834//}
1835
1836
1838(
1839 const labelList& globalToMasterPatch,
1840 const labelList& globalToSlavePatch,
1841 const pointField& locationsInMesh,
1842 const wordList& zonesInMesh,
1843 const pointField& locationsOutsideMesh,
1844 const labelList& selectedSurfaces
1845)
1846{
1847 // Problem: this is based on cached intersection information. This might
1848 // not include the surface we actually want to use. In which case the
1849 // surface would not be seen as intersected!
1850 boolList isBlockedFace;
1851 selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1852
1853 // Determine cell regions
1854 const regionSplit cellRegion(mesh_, isBlockedFace);
1855
1856 // Detect locationsInMesh regions
1857 labelList insideCells(locationsInMesh.size(), -1);
1858 labelList insideRegions(locationsInMesh.size(), -1);
1859 forAll(locationsInMesh, i)
1860 {
1861 insideCells[i] = findCell
1862 (
1863 mesh_,
1864 mergeDistance_*vector::one, //perturbVec,
1865 locationsInMesh[i]
1866 );
1867 if (insideCells[i] != -1)
1868 {
1869 insideRegions[i] = cellRegion[insideCells[i]];
1870 }
1871 reduce(insideRegions[i], maxOp<label>());
1872
1873 if (insideRegions[i] == -1)
1874 {
1875 // See if we can perturb a bit
1876 insideCells[i] = findCell
1877 (
1878 mesh_,
1879 mergeDistance_*vector::one, //perturbVec,
1880 locationsInMesh[i]+mergeDistance_*vector::one
1881 );
1882 if (insideCells[i] != -1)
1883 {
1884 insideRegions[i] = cellRegion[insideCells[i]];
1885 }
1886 reduce(insideRegions[i], maxOp<label>());
1887
1888 if (insideRegions[i] == -1)
1889 {
1891 << "Cannot find locationInMesh " << locationsInMesh[i]
1892 << " on any processor" << exit(FatalError);
1893 }
1894 }
1895 }
1896
1897
1898 // Check that all the locations outside the
1899 // mesh do not conflict with those inside
1900
1901 bool haveLeak = false;
1902 forAll(locationsOutsideMesh, i)
1903 {
1904 // Find the region containing the point
1905 label regioni = findRegion
1906 (
1907 mesh_,
1908 cellRegion,
1909 mergeDistance_*vector::one, //perturbVec,
1910 locationsOutsideMesh[i]
1911 );
1912
1913 if (regioni != -1)
1914 {
1915 // Check for locationsOutsideMesh overlapping with inside ones
1916 if (insideRegions.find(regioni) != -1)
1917 {
1918 haveLeak = true;
1920 << "Outside location " << locationsOutsideMesh[i]
1921 << " in region " << regioni
1922 << " is connected to one of the inside points "
1923 << locationsInMesh << endl;
1924 }
1925 }
1926 }
1927
1928
1929 autoPtr<mapPolyMesh> mapPtr;
1930 if (returnReduce(haveLeak, orOp<bool>()))
1931 {
1932 // Use holeToFace to provide a minimum set of faces needed
1933 // to close hole.
1934
1935 const List<pointField> allLocations
1936 (
1938 (
1939 locationsInMesh,
1940 zonesInMesh,
1941 locationsOutsideMesh
1942 )
1943 );
1944
1945 const labelList blockingFaces(findIndices(isBlockedFace, true));
1946
1947 labelList closureFaces;
1948 labelList closureToBlocked;
1949 autoPtr<mapDistribute> closureMapPtr;
1950 {
1951 const globalIndex globalBlockingFaces(blockingFaces.size());
1952
1953 closureMapPtr = holeToFace::calcClosure
1954 (
1955 mesh_,
1956 allLocations,
1957 blockingFaces,
1958 globalBlockingFaces,
1959 true, // allow erosion
1960
1961 closureFaces,
1962 closureToBlocked
1963 );
1964
1965 if (debug)
1966 {
1967 Pout<< "meshRefinement::blockLeakFaces :"
1968 << " found closure faces:" << closureFaces.size()
1969 << " map:" << closureMapPtr.valid() << endl;
1970 }
1971
1972 if (!closureMapPtr.valid())
1973 {
1975 << "have leak but did not find any closure faces"
1976 << exit(FatalError);
1977 }
1978 }
1979
1980 // Baffle faces
1981 labelList ownPatch(mesh_.nFaces(), -1);
1982 labelList neiPatch(mesh_.nFaces(), -1);
1983
1984 // Keep (global) boundary faces in their patch
1985 {
1986 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1987 for (label patchi = 0; patchi < patches.nNonProcessor(); ++patchi)
1988 {
1989 const polyPatch& pp = patches[patchi];
1990
1991 forAll(pp, i)
1992 {
1993 ownPatch[pp.start()+i] = patchi;
1994 neiPatch[pp.start()+i] = patchi;
1995 }
1996 }
1997 }
1998
1999 const faceZoneMesh& fzs = mesh_.faceZones();
2000
2001 // Mark zone per face
2002 labelList faceToZone(mesh_.nFaces(), -1);
2003 boolList faceToFlip(mesh_.nFaces(), false);
2004 forAll(fzs, zonei)
2005 {
2006 const labelList& addressing = fzs[zonei];
2007 const boolList& flipMap = fzs[zonei].flipMap();
2008
2009 forAll(addressing, i)
2010 {
2011 faceToZone[addressing[i]] = zonei;
2012 faceToFlip[addressing[i]] = flipMap[i];
2013 }
2014 }
2015
2016
2017 // Fetch patch and zone info from blockingFaces
2018 labelList packedOwnPatch(labelUIndList(ownPatch, blockingFaces));
2019 closureMapPtr->distribute(packedOwnPatch);
2020 labelList packedNeiPatch(labelUIndList(neiPatch, blockingFaces));
2021 closureMapPtr->distribute(packedNeiPatch);
2022 labelList packedZone(labelUIndList(faceToZone, blockingFaces));
2023 closureMapPtr->distribute(packedZone);
2024 boolList packedFlip(UIndirectList<bool>(faceToFlip, blockingFaces));
2025 closureMapPtr->distribute(packedFlip);
2026 forAll(closureFaces, i)
2027 {
2028 const label facei = closureFaces[i];
2029 const label sloti = closureToBlocked[i];
2030
2031 if (sloti != -1)
2032 {
2033 // TBD. how to orient own/nei patch?
2034 ownPatch[facei] = packedOwnPatch[sloti];
2035 neiPatch[facei] = packedNeiPatch[sloti];
2036 faceToZone[facei] = packedZone[sloti];
2037 faceToFlip[facei] = packedFlip[sloti];
2038 }
2039 }
2040
2041
2042 // Add faces to faceZone. For now do this outside of createBaffles
2043 // below
2044 {
2045 List<DynamicList<label>> zoneToFaces(fzs.size());
2046 List<DynamicList<bool>> zoneToFlip(fzs.size());
2047
2048 // Add any to-be-patched face
2049 forAll(faceToZone, facei)
2050 {
2051 const label zonei = faceToZone[facei];
2052 if (zonei != -1)
2053 {
2054 zoneToFaces[zonei].append(facei);
2055 zoneToFlip[zonei].append(faceToFlip[facei]);
2056 }
2057 }
2058
2059 forAll(zoneToFaces, zonei)
2060 {
2062 (
2063 fzs[zonei].name(),
2064 zoneToFaces[zonei],
2065 zoneToFlip[zonei],
2066 mesh_
2067 );
2068 }
2069 }
2070
2071 // Put the points of closureFaces into a special pointZone
2072 {
2074 (
2075 UIndirectList<face>(mesh_.faces(), closureFaces),
2076 mesh_.points()
2077 );
2078
2079 // Count number of faces per edge
2080 const labelList nEdgeFaces(countEdgeFaces(pp));
2081
2082 // Freeze all internal points
2083 bitSet isFrozenPoint(mesh_.nPoints());
2084 forAll(nEdgeFaces, edgei)
2085 {
2086 if (nEdgeFaces[edgei] != 1)
2087 {
2088 const edge& e = pp.edges()[edgei];
2089 isFrozenPoint.set(pp.meshPoints()[e[0]]);
2090 isFrozenPoint.set(pp.meshPoints()[e[1]]);
2091 }
2092 }
2093
2094 // Lookup/add pointZone and include its points
2095 pointZoneMesh& pointZones =
2096 const_cast<pointZoneMesh&>(mesh_.pointZones());
2097 const label zonei = addPointZone("frozenPoints");
2098 const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
2099 isFrozenPoint.set(oldSet);
2100
2102 (
2103 mesh_,
2104 isFrozenPoint,
2106 0u
2107 );
2108
2109 // Override addressing
2110 pointZones.clearAddressing();
2111 pointZones[zonei] = isFrozenPoint.sortedToc();
2112 }
2113
2114
2115
2116 // Create baffles for faces
2117 mapPtr = createBaffles(ownPatch, neiPatch);
2118
2120 //{
2121 // // Add newly exposed faces (if not yet in any faceZone!)
2122 // const labelList exposed
2123 // (
2124 // renumber
2125 // (
2126 // mapPtr().reverseFaceMap(),
2127 // blockingFaces
2128 // )
2129 // );
2130 //
2131 // surfaceZonesInfo::addFaceZone
2132 // (
2133 // "frozenFaces",
2134 // exposed,
2135 // boolList(exposed.size(), false),
2136 // mesh_
2137 // );
2138 //}
2139
2140
2141 if (debug&meshRefinement::MESH)
2142 {
2143 const_cast<Time&>(mesh_.time())++;
2144
2145 Pout<< "Writing current mesh to time "
2146 << timeName() << endl;
2147 write
2148 (
2151 (
2154 ),
2155 mesh_.time().path()/timeName()
2156 );
2157 Pout<< "Dumped mesh in = "
2158 << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
2159 }
2160 }
2161 return mapPtr;
2162}
2163
2164
2165// ************************************************************************* //
label n
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition: HashTable.C:137
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
void setSize(const label n, unsigned int val=0u)
Alias for resize()
Definition: PackedList.H:557
A list of faces which address into the list of points.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const labelListList & edgeFaces() const
Return edge-face addressing.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
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 size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:65
void clearAddressing()
Clear addressing.
Definition: ZoneMesh.C:715
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
bool valid() const noexcept
Identical to good(), or bool operator.
Definition: autoPtr.H:272
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
labelList sortedToc() const
The indices of the on bits as a sorted labelList.
Definition: bitSetI.H:533
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
static autoPtr< mapDistribute > calcClosure(const polyMesh &mesh, const List< pointField > &zonePoints, const labelList &blockedFaces, const globalIndex &globalBlockedFaces, const bool erode, labelList &closureFaces, labelList &closureToBlocked)
Optional direct use to generate the set of faces and the method to.
Definition: holeToFace.C:1194
autoPtr< mapPolyMesh > blockLeakFaces(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh, const labelList &selectedSurfaces)
Baffle faces to break any leak from inside to outside.
writeType
Enumeration for what to write. Used as a bit-pattern.
label countFaceDirs(const bitSet &isOutsideFace, const label celli) const
Count number of faces on cell that are in set.
autoPtr< mapPolyMesh > removeGapCells(const scalar planarAngle, const labelList &minSurfaceLevel, const labelList &globalToMasterPatch, const label growIter)
Detect gapRefinement cells and remove them.
void markOutsideFaces(const labelList &cellLevel, const labelList &neiLevel, const labelList &refineCell, bitSet &isOutsideFace) const
Mark faces on interface between set and rest.
debugType
Enumeration for what to debug. Used as a bit-pattern.
void selectIntersectedFaces(const labelList &surfaces, boolList &isBlockedFace) const
Faces currently on boundary or intersected by surface.
void growSet(const labelList &neiLevel, const bitSet &isOutsideFace, labelList &refineCell, label &nRefine) const
Add one layer of cells to set.
static writeType writeLevel()
Get/set write level.
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label nNonProcessor() const
The number of patches before the first processor patch.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
label nInternalFaces() const noexcept
Number of internal faces.
label nFaces() const noexcept
Number of mesh faces.
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:57
static List< pointField > zonePoints(const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh)
Helper: per zone (entry in zonesInMesh) the locations with.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:144
Given list of cells to remove, insert all the topology changes.
Definition: removeCells.H:64
labelList getExposedFaces(const bitSet &removedCell) const
Get labels of faces exposed after cells removal.
Definition: removeCells.C:84
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
static label addFaceZone(const word &name, const labelList &addressing, const boolList &flipMap, polyMesh &mesh)
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
IOoject and searching on triSurface.
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
Definition: triSurface.C:796
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const cellShapeList & cells
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
label nZones
word timeName
Definition: getTimeIndex.H:3
#define WarningInFunction
Report a warning using Foam::Warning.
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
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
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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
Field< vector > vectorField
Specialisation of Field<T> for vector.
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
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.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
runTime write()
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.