Loading...
Searching...
No Matches
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-2025 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"
48#include "distributedTriSurfaceMesh.H"
49
50#include "holeToFace.H"
51#include "refinementParameters.H"
52#include "indirectPrimitivePatch.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
177void Foam::meshRefinement::markOutsideFaces
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
235Foam::label Foam::meshRefinement::countFaceDirs
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
286void Foam::meshRefinement::growSet
287(
288 const labelList& neiLevel,
289 const bitSet& isOutsideFace,
290 labelList& refineCell,
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.good())
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// IOobject::NO_REGISTER
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
458void Foam::meshRefinement::markProximityCandidateFaces
459(
460 const labelList& blockedSurfaces,
461 const List<scalarList>& regionToBlockSize,
462 const labelList& neiLevel,
463 const pointField& neiCc,
464
465 labelList& testFaces
466) const
467{
468 labelListList faceZones(surfaces_.surfaces().size());
469 {
470 // If triSurface do additional zoning based on connectivity
471 for (const label surfi : blockedSurfaces)
472 {
473 const label geomi = surfaces_.surfaces()[surfi];
474 const searchableSurface& s = surfaces_.geometry()[geomi];
475 if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
476 {
477 const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
478 const labelListList& edFaces = surf.edgeFaces();
479 boolList isOpenEdge(edFaces.size(), false);
480 forAll(edFaces, i)
481 {
482 if (edFaces[i].size() == 1)
483 {
484 isOpenEdge[i] = true;
485 }
486 }
487
488 labelList faceZone;
489 const label nZones = surf.markZones(isOpenEdge, faceZone);
490 if (nZones > 1)
491 {
492 faceZones[surfi].transfer(faceZone);
493 }
494 }
495 }
496 }
497
498
499 // Clever limiting of the number of iterations (= max cells in the channel)
500 // since it has too many problematic issues, e.g. with volume refinement
501 // and the real check uses the proper distance anyway just disable.
502 const label nIters = mesh_.globalData().nTotalCells();
503
504
505 // Collect candidate faces (i.e. intersecting any surface and
506 // owner/neighbour not yet refined)
507 //const labelList testFaces(getRefineCandidateFaces(refineCell));
508
509 // Collect segments
510 pointField start(testFaces.size());
511 pointField end(testFaces.size());
512 labelList minLevel(testFaces.size());
513
514 calcCellCellRays
515 (
516 neiCc,
517 neiLevel,
518 testFaces,
519 start,
520 end,
521 minLevel
522 );
523 // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
524 // might add to cell level). Unfortunately we only have minLevel,
525 // not maxLevel. Problem: what if volume refinement only in middle of
526 // channel? Say channel 1m wide with a 0.1m sphere of refinement
527 // Workaround: have dummy surface with e.g. maxLevel 100 to
528 // force nIters to be high enough.
529
530
531 // Test for all intersections (with surfaces of higher gap level than
532 // minLevel) and cache per cell the max surface level and the local normal
533 // on that surface.
534
535 labelList surface1;
536 List<pointIndexHit> hit1;
537 labelList region1;
538 vectorField normal1;
539
540 labelList surface2;
541 List<pointIndexHit> hit2;
542 labelList region2;
543 vectorField normal2;
544
545 surfaces_.findNearestIntersection
546 (
547 blockedSurfaces,
548 start,
549 end,
550
551 surface1,
552 hit1,
553 region1, // local region
554 normal1,
555
556 surface2,
557 hit2,
558 region2, // local region
559 normal2
560 );
561
562
563 // Detect cells that are using multiple surface regions
564 //bitSet isMultiRegion;
565 //detectMultiRegionCells
566 //(
567 // faceZones,
568 // testFaces,
569 //
570 // surface1,
571 // hit1,
572 // region1,
573 //
574 // surface2,
575 // hit2,
576 // region2,
577 //
578 // isMultiRegion
579 //);
580
581
582 label n = 0;
583 forAll(testFaces, i)
584 {
585 if (hit1[i].hit())
586 {
587 n++;
588 }
589 }
590
591 List<wallPoints> faceDist(n);
592 labelList changedFaces(n);
593 n = 0;
594
595 DynamicList<point> originLocation(2);
596 DynamicList<scalar> originDistSqr(2);
597 DynamicList<FixedList<label, 3>> originSurface(2);
598 //DynamicList<point> originNormal(2);
599
600
601 //- To avoid walking through surfaces we mark all faces that have been
602 // intersected. We can either mark only those faces intersecting
603 // blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
604 // any (refinement) surface (this includes e.g. faceZones). This is
605 // much easier since that information is already cached
606 // (meshRefinement::intersectedFaces())
607
608 //bitSet isBlockedFace(mesh_.nFaces());
609 forAll(testFaces, i)
610 {
611 if (hit1[i].hit())
612 {
613 const label facei = testFaces[i];
614 //isBlockedFace.set(facei);
615 const point& fc = mesh_.faceCentres()[facei];
616 const labelList& fz1 = faceZones[surface1[i]];
617
618 originLocation.clear();
619 originDistSqr.clear();
620 //originNormal.clear();
621 originSurface.clear();
622
623 originLocation.append(hit1[i].point());
624 originDistSqr.append(hit1[i].point().distSqr(fc));
625 //originNormal.append(normal1[i]);
626 originSurface.append
627 (
628 FixedList<label, 3>
629 ({
630 surface1[i],
631 region1[i],
632 (fz1.size() ? fz1[hit1[i].index()] : region1[i])
633 })
634 );
635
636 if (hit2[i].hit() && hit1[i] != hit2[i])
637 {
638 const labelList& fz2 = faceZones[surface2[i]];
639 originLocation.append(hit2[i].point());
640 originDistSqr.append(hit2[i].point().distSqr(fc));
641 //originNormal.append(normal2[i]);
642 originSurface.append
643 (
644 FixedList<label, 3>
645 ({
646 surface2[i],
647 region2[i],
648 (fz2.size() ? fz2[hit2[i].index()] : region2[i])
649 })
650 );
651 }
652
653 // Collect all seed data. Currently walking does not look at
654 // surface direction - if so pass in surface normal as well
655 faceDist[n] = wallPoints
656 (
657 originLocation, // origin
658 originDistSqr, // distance to origin
659 originSurface // surface+region+zone
660 //originNormal // normal at origin
661 );
662 changedFaces[n] = facei;
663 n++;
664 }
665 }
666
667
668 // Clear intersection info
669 surface1.clear();
670 hit1.clear();
671 region1.clear();
672 normal1.clear();
673 surface2.clear();
674 hit2.clear();
675 region2.clear();
676 normal2.clear();
677
678
679 List<wallPoints> allFaceInfo(mesh_.nFaces());
680 List<wallPoints> allCellInfo(mesh_.nCells());
681
682 // Any refinement surface (even a faceZone) should stop the gap walking.
683 // This is exactly the information which is cached in the surfaceIndex_
684 // field.
685 const bitSet isBlockedFace(intersectedFaces());
686
687 wallPoints::trackData td(isBlockedFace, regionToBlockSize);
688 FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
689 (
690 mesh_,
691 changedFaces,
692 faceDist,
693 allFaceInfo,
694 allCellInfo,
695 0, // max iterations
696 td
697 );
698 wallDistCalc.iterate(nIters);
699
700
701 // Extract faces of visited cells
702
703 bitSet isProximityFace(mesh_.nFaces(), false);
704
705 // Make sure to include initial intersected ones
706 isProximityFace.set(testFaces);
707
708 forAll(allCellInfo, celli)
709 {
710 if (allCellInfo[celli].valid(wallDistCalc.data()))
711 {
712 isProximityFace.set(mesh_.cells()[celli]);
713 }
714 }
715
716 syncTools::syncFaceList
717 (
718 mesh_,
719 isProximityFace,
720 orEqOp<unsigned int>(),
721 0u
722 );
723
724 testFaces = isProximityFace.toc();
725}
726
727
728Foam::label Foam::meshRefinement::markFakeGapRefinement
729(
730 const scalar planarCos,
731 const labelList& blockedSurfaces,
732 const label nAllowRefine,
733 const labelList& neiLevel,
734 const pointField& neiCc,
735
736 labelList& refineCell,
737 label& nRefine
738) const
739{
740 // Re-work the blockLevel of the blockedSurfaces into a length scale
741 // and a number of cells to cross
742 List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
743 scalar maxBlockSize(-1);
744 for (const label surfi : blockedSurfaces)
745 {
746 const label geomi = surfaces_.surfaces()[surfi];
747 const searchableSurface& s = surfaces_.geometry()[geomi];
748 const label nRegions = s.regions().size();
749
750 regionToBlockSize[surfi].setSize(nRegions, -1);
751 for (label regioni = 0; regioni < nRegions; regioni++)
752 {
753 const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
754 const label bLevel = surfaces_.blockLevel()[globalRegioni];
755
756 // If valid blockLevel specified
757 if (bLevel != labelMax)
758 {
759 regionToBlockSize[surfi][regioni] =
760 meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
761 }
762
763 maxBlockSize = max(maxBlockSize, regionToBlockSize[surfi][regioni]);
764
765 //const label mLevel = surfaces_.maxLevel()[globalRegioni];
769 //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
770 //{
771 // const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
772 //}
773 }
774 }
775
776 // Collect candidate faces (i.e. intersecting any surface and
777 // owner/neighbour not yet refined) and their cells
778 labelList cellMap;
779 {
780 labelList testFaces(getRefineCandidateFaces(refineCell));
781
782 // Extend the set by faceCellFace wave upto given 3*blockLevel size
783 markProximityCandidateFaces
784 (
785 blockedSurfaces,
786 regionToBlockSize,
787 neiLevel,
788 neiCc,
789 testFaces
790 );
791
792
793 bitSet isTestCell(mesh_.nCells());
794 for (const label facei : testFaces)
795 {
796 const label own = mesh_.faceOwner()[facei];
797 if (refineCell[own] == -1)
798 {
799 isTestCell.set(own);
800 }
801 if (mesh_.isInternalFace(facei))
802 {
803 const label nei = mesh_.faceNeighbour()[facei];
804 if (refineCell[nei] == -1)
805 {
806 isTestCell.set(nei);
807 }
808 }
809 }
810 cellMap = isTestCell.sortedToc();
811 }
812
813
814 // Shoot rays in all 6 directions
815 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
816 //
817 // Almost exact copy of meshRefinement::markInternalGapRefinement.
818 // Difference is only in the length scale (now based on blockLevel)
819
820
821 // Find per cell centre the nearest point and normal on the surfaces
822 List<pointIndexHit> nearInfo;
823 vectorField nearNormal;
824 labelList nearSurface;
825 labelList nearRegion;
826 {
827 surfaces_.findNearestRegion
828 (
829 blockedSurfaces,
830 pointField(mesh_.cellCentres(), cellMap),
831 scalarField(cellMap.size(), maxBlockSize), // use uniform for now
832 nearSurface,
833 nearInfo,
834 nearRegion,
835 nearNormal
836 );
837 }
838
839
840 DynamicList<label> map(nearInfo.size());
841 DynamicField<point> rayStart(nearInfo.size());
842 DynamicField<point> rayEnd(nearInfo.size());
843 DynamicField<scalar> gapSize(nearInfo.size());
844
845 DynamicField<point> rayStart2(nearInfo.size());
846 DynamicField<point> rayEnd2(nearInfo.size());
847 DynamicField<scalar> gapSize2(nearInfo.size());
848
849 label nTestCells = 0;
850
851 forAll(nearInfo, i)
852 {
853 if (nearInfo[i].hit())
854 {
855 const label globalRegioni = surfaces_.globalRegion
856 (
857 nearSurface[i],
858 nearRegion[i]
859 );
860
861 // Combine info from shell and surface
862 const label bLevel = surfaces_.blockLevel()[globalRegioni];
863 const FixedList<label, 3> gapInfo
864 ({
865 100,
866 bLevel,
867 100
868 });
869 const volumeType gapMode(volumeType::MIXED);
870
871 // Store wanted number of cells in gap
872 const label celli = cellMap[i];
873 const label cLevel = meshCutter_.cellLevel()[celli];
874
875 // Construct one or more rays to test for oppositeness
876 const label nRays = generateRays
877 (
878 false, // ignore actual surface normal
879 nearInfo[i].hitPoint(),
880 nearNormal[i],
881 gapInfo,
882 gapMode,
883
884 mesh_.cellCentres()[celli],
885 cLevel,
886
887 rayStart,
888 rayEnd,
889 gapSize,
890
891 rayStart2,
892 rayEnd2,
893 gapSize2
894 );
895 if (nRays > 0)
896 {
897 nTestCells++;
898 for (label j = 0; j < nRays; j++)
899 {
900 map.append(i);
901 }
902 }
903 }
904 }
905
906 Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
907 << " cells for testing out of "
908 << mesh_.globalData().nTotalCells() << endl;
909 map.shrink();
910 rayStart.shrink();
911 rayEnd.shrink();
912 gapSize.shrink();
913
914 rayStart2.shrink();
915 rayEnd2.shrink();
916 gapSize2.shrink();
917
918 cellMap = labelUIndList(cellMap, map)();
919 nearNormal = UIndirectList<vector>(nearNormal, map)();
920 nearInfo.clear();
921 nearSurface.clear();
922 nearRegion.clear();
923
924
925 // Do intersections in pairs
926 labelList surf1;
927 labelList region1;
928 List<pointIndexHit> hit1;
929 vectorField normal1;
930 surfaces_.findNearestIntersection
931 (
932 rayStart,
933 rayEnd,
934 surf1,
935 region1,
936 hit1,
937 normal1
938 );
939
940 labelList surf2;
941 labelList region2;
942 List<pointIndexHit> hit2;
943 vectorField normal2;
944 surfaces_.findNearestIntersection
945 (
946 rayStart2,
947 rayEnd2,
948 surf2,
949 region2,
950 hit2,
951 normal2
952 );
953
954 // Extract cell based gap size
955 const label oldNRefine = nRefine;
956 forAll(surf1, i)
957 {
958 if (surf1[i] != -1 && surf2[i] != -1)
959 {
960 // Found intersections with surface. Check for
961 // - small gap
962 // - coplanar normals
963
964 const label celli = cellMap[i];
965 if (celli != -1 && refineCell[celli] == -1)
966 {
967 const scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
968
969 const scalar maxGapSize
970 (
971 max
972 (
973 regionToBlockSize[surf1[i]][region1[i]],
974 regionToBlockSize[surf2[i]][region2[i]]
975 )
976 );
977
978 if
979 (
980 (mag(normal1[i]&normal2[i]) > planarCos)
981 && (d2 < Foam::sqr(maxGapSize))
982 )
983 {
984 if
985 (
986 !markForRefine
987 (
988 0, // mark level
989 nAllowRefine,
990 refineCell[celli],
991 nRefine
992 )
993 )
994 {
995 if (debug)
996 {
997 Pout<< "Stopped refining since reaching my cell"
998 << " limit of " << mesh_.nCells()+7*nRefine
999 << endl;
1000 }
1001 break;
1002 }
1003 }
1004 }
1005 }
1006 }
1007
1008 if
1009 (
1010 returnReduce(nRefine, sumOp<label>())
1011 > returnReduce(nAllowRefine, sumOp<label>())
1012 )
1013 {
1014 Info<< "Reached refinement limit." << endl;
1015 }
1016
1017 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1018}
1019
1020
1021Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
1022(
1023 const scalar planarAngle,
1024 const labelList& minSurfaceLevel,
1025 const labelList& globalToMasterPatch,
1026 const label growIter
1027)
1028{
1029 // Swap neighbouring cell centres and cell level
1030 labelList neiLevel(mesh_.nBoundaryFaces());
1031 pointField neiCc(mesh_.nBoundaryFaces());
1032 calcNeighbourData(neiLevel, neiCc);
1033
1034 labelList refineCell(mesh_.nCells(), -1);
1035 label nRefine = 0;
1036 //markProximityRefinement
1037 //(
1038 // Foam::cos(degToRad(planarAngle)),
1039 //
1040 // minSurfaceLevel, // surface min level
1041 // labelList(minSurfaceLevel.size(), labelMax), // surfaceGapLevel
1042 //
1043 // labelMax/Pstream::nProcs(), //nAllowRefine,
1044 // neiLevel,
1045 // neiCc,
1046 //
1047 // refineCell,
1048 // nRefine
1049 //);
1050
1051
1052 // Determine minimum blockLevel per surface
1053 Map<label> surfToBlockLevel;
1054
1055 forAll(surfaces_.surfaces(), surfi)
1056 {
1057 const label geomi = surfaces_.surfaces()[surfi];
1058 const searchableSurface& s = surfaces_.geometry()[geomi];
1059 const label nRegions = s.regions().size();
1060
1061 label minBlockLevel = labelMax;
1062 for (label regioni = 0; regioni < nRegions; regioni++)
1063 {
1064 const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
1065 minBlockLevel = min
1066 (
1067 minBlockLevel,
1068 surfaces_.blockLevel()[globalRegioni]
1069 );
1070 }
1071
1072 if (minBlockLevel < labelMax)
1073 {
1074 surfToBlockLevel.insert(surfi, minBlockLevel);
1075 }
1076 }
1077
1078
1079 //markProximityRefinementWave
1080 //(
1081 // Foam::cos(degToRad(planarAngle)),
1082 // surfToBlockLevel.sortedToc(),
1083 //
1084 // labelMax/Pstream::nProcs(), //nAllowRefine,
1085 // neiLevel,
1086 // neiCc,
1087 //
1088 // refineCell,
1089 // nRefine
1090 //);
1091
1092 markFakeGapRefinement
1093 (
1094 Foam::cos(degToRad(planarAngle)),
1095 surfToBlockLevel.sortedToc(),
1096
1097 labelMax/Pstream::nProcs(), //nAllowRefine,
1098 neiLevel,
1099 neiCc,
1100
1101 refineCell,
1102 nRefine
1103 );
1104
1105
1106 Info<< "Marked for blocking due to close opposite surfaces : "
1107 << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1108
1109 // Remove outliers, i.e. cells with all points exposed
1110 if (growIter)
1111 {
1112 labelList oldRefineCell(refineCell);
1113
1114 // Pass1: extend the set to fill in gaps
1115 bitSet isOutsideFace;
1116 for (label iter = 0; iter < growIter; iter++)
1117 {
1118 // Get outside faces
1119 markOutsideFaces
1120 (
1121 meshCutter_.cellLevel(),
1122 neiLevel,
1123 refineCell,
1124 isOutsideFace
1125 );
1126 // Extend with cells with three outside faces
1127 growSet(neiLevel, isOutsideFace, refineCell, nRefine);
1128 }
1129
1130
1131 // Pass2: erode back to original set if pass1 didn't help
1132 for (label iter = 0; iter < growIter; iter++)
1133 {
1134 // Get outside faces. Ignore cell level.
1135 markOutsideFaces
1136 (
1137 labelList(mesh_.nCells(), 0),
1138 labelList(neiLevel.size(), 0),
1139 refineCell,
1140 isOutsideFace
1141 );
1142
1143 // Unmark cells with three or more outside faces
1144 for (label celli = 0; celli < mesh_.nCells(); celli++)
1145 {
1146 if (refineCell[celli] != -1 && oldRefineCell[celli] == -1)
1147 {
1148 if (countFaceDirs(isOutsideFace, celli) >= 3)
1149 {
1150 refineCell[celli] = -1;
1151 --nRefine;
1152 }
1153 }
1154 }
1155 }
1156
1157 Info<< "Marked for blocking after filtering : "
1158 << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1159 }
1160
1161
1162 // Determine patch for every mesh face
1163 const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
1164 labelList unnamedSurfaces(surfaceZonesInfo::getUnnamedSurfaces(surfZones));
1165 const label defaultRegion(surfaces_.globalRegion(unnamedSurfaces[0], 0));
1166
1167 const labelList nearestRegion
1168 (
1169 nearestIntersection
1170 (
1171 unnamedSurfaces,
1172 defaultRegion
1173 )
1174 );
1175
1176 // Pack
1177 labelList cellsToRemove(nRefine);
1178 nRefine = 0;
1179
1180 forAll(refineCell, cellI)
1181 {
1182 if (refineCell[cellI] != -1)
1183 {
1184 cellsToRemove[nRefine++] = cellI;
1185 }
1186 }
1187
1188 // Remove cells
1189 removeCells cellRemover(mesh_);
1190 labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1191
1192 labelList exposedPatches(exposedFaces.size());
1193 forAll(exposedFaces, i)
1194 {
1195 label facei = exposedFaces[i];
1196 exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1197 }
1198
1199 return doRemoveCells
1200 (
1201 cellsToRemove,
1202 exposedFaces,
1203 exposedPatches,
1204 cellRemover
1205 );
1206}
1207
1208
1209void Foam::meshRefinement::selectIntersectedFaces
1210(
1211 const labelList& selectedSurfaces,
1212 boolList& isBlockedFace
1213) const
1214{
1215 // Like meshRefinement::selectSeparatedCoupledFaces. tbd: convert to bitSet
1216
1217 // Check that no connection between inside and outside points
1218 isBlockedFace.setSize(mesh_.nFaces(), false);
1219
1220 // Block off separated couples.
1221 selectSeparatedCoupledFaces(isBlockedFace);
1222
1223 // Block off intersections with selected surfaces
1224
1225 // Mark per face (for efficiency)
1226 boolList isSelectedSurf(surfaces_.surfaces().size(), false);
1227 UIndirectList<bool>(isSelectedSurf, selectedSurfaces) = true;
1228
1229 forAll(surfaceIndex_, facei)
1230 {
1231 const label surfi = surfaceIndex_[facei];
1232 if (surfi != -1 && isSelectedSurf[surfi])
1233 {
1234 isBlockedFace[facei] = true;
1235 }
1236 }
1237}
1238
1239
1240//Foam::labelList Foam::meshRefinement::detectLeakCells
1241//(
1242// const boolList& isBlockedFace,
1243// const labelList& leakFaces,
1244// const labelList& seedCells
1245//) const
1246//{
1247// int dummyTrackData = 0;
1248// List<topoDistanceData<label>> allFaceInfo(mesh_.nFaces());
1249// List<topoDistanceData<label>> allCellInfo(mesh_.nCells());
1250//
1251// // Block faces
1252// forAll(isBlockedFace, facei)
1253// {
1254// if (isBlockedFace[facei])
1255// {
1256// allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1257// }
1258// }
1259// for (const label facei : leakFaces)
1260// {
1261// allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1262// }
1263//
1264// // Walk from inside cell
1265// DynamicList<topoDistanceData<label>> faceDist;
1266// DynamicList<label> cFaces1;
1267// for (const label celli : seedCells)
1268// {
1269// if (celli != -1)
1270// {
1271// const labelList& cFaces = mesh_.cells()[celli];
1272// faceDist.reserve(cFaces.size());
1273// cFaces1.reserve(cFaces.size());
1274//
1275// for (const label facei : cFaces)
1276// {
1277// if (!allFaceInfo[facei].valid(dummyTrackData))
1278// {
1279// cFaces1.append(facei);
1280// faceDist.append(topoDistanceData<label>(0, 123));
1281// }
1282// }
1283// }
1284// }
1285//
1286// // Walk through face-cell wave till all cells are reached
1287// FaceCellWave<topoDistanceData<label>> wallDistCalc
1288// (
1289// mesh_,
1290// cFaces1,
1291// faceDist,
1292// allFaceInfo,
1293// allCellInfo,
1294// mesh_.globalData().nTotalCells()+1 // max iterations
1295// );
1296//
1297// label nRemove = 0;
1298// for (const label facei : leakFaces)
1299// {
1300// const label own = mesh_.faceOwner()[facei];
1301// if (!allCellInfo[own].valid(dummyTrackData))
1302// {
1303// nRemove++;
1304// }
1305// if (mesh_.isInternalFace(facei))
1306// {
1307// const label nei = mesh_.faceNeighbour()[facei];
1308// if (!allCellInfo[nei].valid(dummyTrackData))
1309// {
1310// nRemove++;
1311// }
1312// }
1313// }
1314//
1315// labelList cellsToRemove(nRemove);
1316// nRemove = 0;
1317// for (const label facei : leakFaces)
1318// {
1319// const label own = mesh_.faceOwner()[facei];
1320// if (!allCellInfo[own].valid(dummyTrackData))
1321// {
1322// cellsToRemove[nRemove++] = own;
1323// }
1324// if (mesh_.isInternalFace(facei))
1325// {
1326// const label nei = mesh_.faceNeighbour()[facei];
1327// if (!allCellInfo[nei].valid(dummyTrackData))
1328// {
1329// cellsToRemove[nRemove++] = nei;
1330// }
1331// }
1332// }
1333//
1334// if (debug)
1335// {
1336// volScalarField fld
1337// (
1338// IOobject
1339// (
1340// "cellsToKeep",
1341// mesh_.time().timeName(),
1342// mesh_,
1343// IOobject::NO_READ,
1344// IOobject::NO_WRITE
1345// ),
1346// mesh_,
1347// dimensionedScalar(dimless, Zero)
1348// );
1349// forAll(allCellInfo, celli)
1350// {
1351// if (allCellInfo[celli].valid(dummyTrackData))
1352// {
1353// fld[celli] = allCellInfo[celli].distance();
1354// }
1355// }
1356// forAll(fld.boundaryField(), patchi)
1357// {
1358// const polyPatch& pp = mesh_.boundaryMesh()[patchi];
1359// SubList<topoDistanceData<label>> p(pp.patchSlice(allFaceInfo));
1360// scalarField pfld
1361// (
1362// fld.boundaryField()[patchi].size(),
1363// Zero
1364// );
1365// forAll(pfld, i)
1366// {
1367// if (p[i].valid(dummyTrackData))
1368// {
1369// pfld[i] = 1.0*p[i].distance();
1370// }
1371// }
1372// fld.boundaryFieldRef()[patchi] == pfld;
1373// }
1374// //Note: do not swap cell values so do not do
1375// //fld.correctBoundaryConditions();
1376// Pout<< "Writing distance field for initial cells "
1377// << seedCells << " to " << fld.objectPath() << endl;
1378// fld.write();
1379// }
1380//
1381// return cellsToRemove;
1382//}
1383//
1384//
1385//Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLeakCells
1386//(
1387// const labelList& globalToMasterPatch,
1388// const labelList& globalToSlavePatch,
1389// const pointField& locationsInMesh,
1390// const wordList& zonesInMesh,
1391// const pointField& locationsOutsideMesh,
1392// const labelList& selectedSurfaces
1393//)
1394//{
1395// boolList isBlockedFace;
1396// selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1397//
1398// // Determine cell regions
1399// const regionSplit cellRegion(mesh_, isBlockedFace);
1400//
1401// // Detect locationsInMesh regions
1402// labelList insideCells(locationsInMesh.size(), -1);
1403// labelList insideRegions(locationsInMesh.size(), -1);
1404// forAll(locationsInMesh, i)
1405// {
1406// insideCells[i] = findCell
1407// (
1408// mesh_,
1409// mergeDistance_*vector::one, //perturbVec,
1410// locationsInMesh[i]
1411// );
1412// if (insideCells[i] != -1)
1413// {
1414// insideRegions[i] = cellRegion[insideCells[i]];
1415// }
1416// reduce(insideRegions[i], maxOp<label>());
1417//
1418// if (insideRegions[i] == -1)
1419// {
1420// // See if we can perturb a bit
1421// insideCells[i] = findCell
1422// (
1423// mesh_,
1424// mergeDistance_*vector::one, //perturbVec,
1425// locationsInMesh[i]+mergeDistance_*vector::one
1426// );
1427// if (insideCells[i] != -1)
1428// {
1429// insideRegions[i] = cellRegion[insideCells[i]];
1430// }
1431// reduce(insideRegions[i], maxOp<label>());
1432//
1433// if (insideRegions[i] == -1)
1434// {
1435// FatalErrorInFunction
1436// << "Cannot find locationInMesh " << locationsInMesh[i]
1437// << " on any processor" << exit(FatalError);
1438// }
1439// }
1440// }
1441//
1442//
1443// // Check that all the locations outside the
1444// // mesh do not conflict with those inside
1445//
1446// bool haveLeak = false;
1447// forAll(locationsOutsideMesh, i)
1448// {
1449// // Find the region containing the point
1450// label regioni = findRegion
1451// (
1452// mesh_,
1453// cellRegion,
1454// mergeDistance_*vector::one, //perturbVec,
1455// locationsOutsideMesh[i]
1456// );
1457//
1458// if (regioni != -1)
1459// {
1460// // Check for locationsOutsideMesh overlapping with inside ones
1461// if (insideRegions.find(regioni) != -1)
1462// {
1463// haveLeak = true;
1464// WarningInFunction
1465// << "Outside location " << locationsOutsideMesh[i]
1466// << " in region " << regioni
1467// << " is connected to one of the inside points "
1468// << locationsInMesh << endl;
1469// }
1470// }
1471// }
1472//
1473//
1474// autoPtr<mapPolyMesh> mapPtr;
1475// if (returnReduceOr(haveLeak))
1476// {
1477// // Use shortestPathSet to provide a minimum set of faces needed
1478// // to close hole. Tbd: maybe directly use wave?
1479// meshSearch searchEngine(mesh_);
1480// shortestPathSet leakPath
1481// (
1482// "leakPath",
1483// mesh_,
1484// searchEngine,
1485// coordSet::coordFormatNames[coordSet::coordFormat::DISTANCE],
1486// true,
1487// 50, // tbd. Number of iterations. This is the maximum
1488// // number of faces in the leak hole
1489//
1490// //pbm.groupPatchIDs()["wall"], // patch to grow from
1491// meshedPatches(), // patch to grow from
1492//
1493// locationsInMesh,
1494// locationsOutsideMesh,
1495// isBlockedFace
1496// );
1497//
1498//
1499// // Use leak path to find minimum set of cells to delete
1500// const labelList cellsToRemove
1501// (
1502// detectLeakCells
1503// (
1504// isBlockedFace,
1505// leakPath.leakFaces(),
1506// insideCells
1507// )
1508// );
1509//
1510// // Re-do intersections to find nearest unnamed surface
1511// const label defaultRegion
1512// (
1513// surfaces().globalRegion
1514// (
1515// selectedSurfaces[0],
1516// 0
1517// )
1518// );
1519//
1520// const labelList nearestRegion
1521// (
1522// nearestIntersection
1523// (
1524// selectedSurfaces,
1525// defaultRegion
1526// )
1527// );
1528//
1529//
1530// // Remove cells
1531// removeCells cellRemover(mesh_);
1532// const labelList exposedFaces
1533// (
1534// cellRemover.getExposedFaces(cellsToRemove)
1535// );
1536//
1537// labelList exposedPatches(exposedFaces.size());
1538// forAll(exposedFaces, i)
1539// {
1540// label facei = exposedFaces[i];
1541// exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1542// }
1543//
1544// mapPtr = doRemoveCells
1545// (
1546// cellsToRemove,
1547// exposedFaces,
1548// exposedPatches,
1549// cellRemover
1550// );
1551//
1552//
1553// // Put the exposed faces into a special faceZone
1554// {
1555// // Add to "frozenFaces" zone
1556// faceZoneMesh& faceZones = mesh_.faceZones();
1557//
1558// // Get current frozen faces (if any)
1559// bitSet isFrozenFace(mesh_.nFaces());
1560// label zonei = faceZones.findZoneID("frozenFaces");
1561// if (zonei != -1)
1562// {
1563// const bitSet oldSet(mesh_.nFaces(), faceZones[zonei]);
1564// isFrozenFace.set(oldSet);
1565// }
1566//
1567// // Add newly exposed faces (if not yet in any faceZone!)
1568// const labelList exposed
1569// (
1570// renumber
1571// (
1572// mapPtr().reverseFaceMap(),
1573// exposedFaces
1574// )
1575// );
1576// bitSet isZonedFace(mesh_.nFaces(), faceZones.zoneMap().toc());
1577// for (const label facei : exposed)
1578// {
1579// if (!isZonedFace[facei])
1580// {
1581// isFrozenFace.set(facei);
1582// }
1583// }
1584//
1585// syncTools::syncFaceList
1586// (
1587// mesh_,
1588// isFrozenFace,
1589// orEqOp<unsigned int>(),
1590// 0u
1591// );
1592//
1593// // Add faceZone if non existing
1594// faceZones.clearAddressing();
1595// if (zonei == -1)
1596// {
1597// zonei = faceZones.size();
1598//
1599// faceZones.emplace_back
1600// (
1601// "frozenFaces", // name
1602// zonei, // index
1603// faceZones
1604// );
1605// }
1606//
1607// // Update faceZone with new contents
1608// labelList frozenFaces(isFrozenFace.toc());
1609// boolList frozenFlip(frozenFaces.size(), false);
1610//
1611// faceZones[zonei].resetAddressing
1612// (
1613// std::move(frozenFaces),
1614// std::move(frozenFlip)
1615// );
1616// }
1617//
1618//
1619// //// Put the exposed points into a special pointZone
1620// //if (false)
1621// //{
1622// // const labelList meshFaceIDs
1623// // (
1624// // renumber
1625// // (
1626// // mapPtr().reverseFaceMap(),
1627// // exposedFaces
1628// // )
1629// // );
1630// // const uindirectPrimitivePatch pp
1631// // (
1632// // UIndirectList<face>(mesh_.faces(), meshFaceIDs),
1633// // mesh_.points()
1634// // );
1635// //
1636// // // Count number of faces per edge
1637// // const labelListList& edgeFaces = pp.edgeFaces();
1638// // labelList nEdgeFaces(edgeFaces.size());
1639// // forAll(edgeFaces, edgei)
1640// // {
1641// // nEdgeFaces[edgei] = edgeFaces[edgei].size();
1642// // }
1643// //
1644// // // Sync across processor patches
1645// // if (Pstream::parRun())
1646// // {
1647// // const globalMeshData& globalData = mesh_.globalData();
1648// // const mapDistribute& map = globalData.globalEdgeSlavesMap();
1649// // const indirectPrimitivePatch& cpp =
1650// // globalData.coupledPatch();
1651// //
1652// // // Match pp edges to coupled edges
1653// // labelList patchEdges;
1654// // labelList coupledEdges;
1655// // bitSet sameEdgeOrientation;
1656// // PatchTools::matchEdges
1657// // (
1658// // pp,
1659// // cpp,
1660// // patchEdges,
1661// // coupledEdges,
1662// // sameEdgeOrientation
1663// // );
1664// //
1665// //
1666// // // Convert patch-edge data into cpp-edge data
1667// // labelList coupledNEdgeFaces(map.constructSize(), Zero);
1668// // UIndirectList<label>(coupledNEdgeFaces, coupledEdges) =
1669// // UIndirectList<label>(nEdgeFaces, patchEdges);
1670// //
1671// // // Synchronise
1672// // globalData.syncData
1673// // (
1674// // coupledNEdgeFaces,
1675// // globalData.globalEdgeSlaves(),
1676// // globalData.globalEdgeTransformedSlaves(),
1677// // map,
1678// // plusEqOp<label>()
1679// // );
1680// //
1681// // // Convert back from cpp-edge to patch-edge
1682// // UIndirectList<label>(nEdgeFaces, patchEdges) =
1683// // UIndirectList<label>(coupledNEdgeFaces, coupledEdges);
1684// // }
1685// //
1686// // // Freeze all internal points
1687// // bitSet isFrozenPoint(mesh_.nPoints());
1688// // forAll(nEdgeFaces, edgei)
1689// // {
1690// // if (nEdgeFaces[edgei] != 1)
1691// // {
1692// // const edge& e = pp.edges()[edgei];
1693// // isFrozenPoint.set(pp.meshPoints()[e[0]]);
1694// // isFrozenPoint.set(pp.meshPoints()[e[1]]);
1695// // }
1696// // }
1697// // // Add to "frozenPoints" zone
1698// // pointZoneMesh& pointZones = mesh_.pointZones();
1699// // pointZones.clearAddressing();
1700// // label zonei = pointZones.findZoneID("frozenPoints");
1701// // if (zonei != -1)
1702// // {
1703// // const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
1704// // // Add to isFrozenPoint
1705// // isFrozenPoint.set(oldSet);
1706// // }
1707// //
1708// // syncTools::syncPointList
1709// // (
1710// // mesh_,
1711// // isFrozenPoint,
1712// // orEqOp<unsigned int>(),
1713// // 0u
1714// // );
1715// //
1716// // if (zonei == -1)
1717// // {
1718// // zonei = pointZones.size();
1719// //
1720// // pointZones.emplace_back
1721// // (
1722// // "frozenPoints", // name
1723// // isFrozenPoint.sortedToc(), // addressing
1724// // zonei, // index
1725// // pointZones // pointZoneMesh
1726// // );
1727// // }
1728// //}
1729//
1730//
1731// if (debug&meshRefinement::MESH)
1732// {
1733// const_cast<Time&>(mesh_.time())++;
1734//
1735// Pout<< "Writing current mesh to time "
1736// << timeName() << endl;
1737// write
1738// (
1739// meshRefinement::debugType(debug),
1740// meshRefinement::writeType
1741// (
1742// meshRefinement::writeLevel()
1743// | meshRefinement::WRITEMESH
1744// ),
1745// mesh_.time().path()/timeName()
1746// );
1747// Pout<< "Dumped mesh in = "
1748// << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
1749// }
1750// }
1751// return mapPtr;
1752//}
1753
1754
1755//Foam::autoPtr<Foam::mapDistribute> Foam::meshRefinement::nearestFace
1756//(
1757// const globalIndex& globalSeedFaces,
1758// const labelList& seedFaces,
1759// const labelList& closureFaces
1760//) const
1761//{
1762// // Takes a set of faces for which we have information (seedFaces,
1763// // globalSeedFaces - these are e.g. intersected faces) and walks out
1764// // across nonSeedFace. Returns for
1765// // every nonSeedFace the nearest seed face (in global indexing).
1766// // Used e.g. in hole closing. Assumes that seedFaces, closureFaces
1767// // are a small subset of the master
1768// // so are not large - it uses edge addressing on the closureFaces
1769//
1770//
1771// if (seedFaces.size() != globalSeedFaces.localSize())
1772// {
1773// FatalErrorInFunction << "problem : seedFaces:" << seedFaces.size()
1774// << " globalSeedFaces:" << globalSeedFaces.localSize()
1775// << exit(FatalError);
1776// }
1777//
1778// //// Mark mesh points that are used by any closureFaces. This is for
1779// //// initial filtering
1780// //bitSet isNonSeedPoint(mesh.nPoints());
1781// //for (const label facei : closureFaces)
1782// //{
1783// // isNonSeedPoint.set(mesh_.faces()[facei]);
1784// //}
1785// //syncTools::syncPointList
1786// //(
1787// // mesh_,
1788// // isNonSeedPoint,
1789// // orEqOp<unsigned int>(),
1790// // 0u
1791// //);
1792//
1793// // Make patch
1794// const uindirectPrimitivePatch pp
1795// (
1796// IndirectList<face>(mesh_.faces(), closureFaces),
1797// mesh_.points()
1798// );
1799// const edgeList& edges = pp.edges();
1800// const labelList& mp = pp.meshPoints();
1801// const label nBndEdges = pp.nEdges() - pp.nInternalEdges();
1802//
1803// // For all faces in seedFaces mark the edge with a face. No special
1804// // handling for multiple faces sharing the edge - first one wins
1805// EdgeMap<label> edgeMap(pp.nEdges());
1806// for (const label facei : seedFaces)
1807// {
1808// const label globalFacei = globalSeedFaces.toGlobal(facei);
1809// const face& f = mesh_.faces()[facei];
1810// forAll(f, fp)
1811// {
1812// label nextFp = f.fcIndex(fp);
1813// edgeMap.insert(edge(f[fp], f[nextFp]), globalFacei);
1814// }
1815// }
1816// syncTools::syncEdgeMap(mesh_, edgeMap, maxEqOp<label>());
1817//
1818//
1819//
1820// // Seed
1821// DynamicList<label> initialEdges(2*nBndEdges);
1822// DynamicList<edgeTopoDistanceData<label, uindirectPrimitivePatch>>
1823// initialEdgesInfo(2*nBndEdges);
1824// forAll(edges, edgei)
1825// {
1826// const edge& e = edges[edgei];
1827// const edge meshE = edge(mp[e[0]], mp[e[1]]);
1828//
1829// const auto iter = edgeMap.cfind(meshE);
1830// if (iter.good())
1831// {
1832// initialEdges.append(edgei);
1833// initialEdgesInfo.append
1834// (
1835// edgeTopoDistanceData<label, uindirectPrimitivePatch>
1836// (
1837// 0, // distance
1838// iter.val() // globalFacei
1839// )
1840// );
1841// }
1842// }
1843//
1844// // Data on all edges and faces
1845// List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allEdgeInfo
1846// (
1847// pp.nEdges()
1848// );
1849// List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allFaceInfo
1850// (
1851// pp.size()
1852// );
1853//
1854// // Walk
1855// PatchEdgeFaceWave
1856// <
1857// uindirectPrimitivePatch,
1858// edgeTopoDistanceData<label, uindirectPrimitivePatch>
1859// > calc
1860// (
1861// mesh_,
1862// pp,
1863// initialEdges,
1864// initialEdgesInfo,
1865// allEdgeInfo,
1866// allFaceInfo,
1867// returnReduce(pp.nEdges(), sumOp<label>())
1868// );
1869//
1870//
1871// // Per non-seed face the seed face
1872// labelList closureToSeed(pp.size(), -1);
1873// forAll(allFaceInfo, facei)
1874// {
1875// if (allFaceInfo[facei].valid(calc.data()))
1876// {
1877// closureToSeed[facei] = allFaceInfo[facei].data();
1878// }
1879// }
1880//
1881// List<Map<label>> compactMap;
1882// return autoPtr<mapDistribute>::New
1883// (
1884// globalSeedFaces,
1885// closureToSeed,
1886// compactMap
1887// );
1888//}
1889
1890
1891Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::blockLeakFaces
1892(
1893 const labelList& globalToMasterPatch,
1894 const labelList& globalToSlavePatch,
1895 const pointField& locationsInMesh,
1896 const wordList& zonesInMesh,
1897 const pointField& locationsOutsideMesh,
1898 const labelList& selectedSurfaces
1899)
1900{
1901 // Problem: this is based on cached intersection information. This might
1902 // not include the surface we actually want to use. In which case the
1903 // surface would not be seen as intersected!
1904 boolList isBlockedFace;
1905 selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1906
1907 // Determine cell regions
1908 const regionSplit cellRegion(mesh_, isBlockedFace);
1909
1910 // Detect locationsInMesh regions
1911 labelList insideCells(locationsInMesh.size(), -1);
1912 labelList insideRegions(locationsInMesh.size(), -1);
1913 forAll(locationsInMesh, i)
1914 {
1915 insideCells[i] = findCell
1916 (
1917 mesh_,
1918 mergeDistance_*vector::one, //perturbVec,
1919 locationsInMesh[i]
1920 );
1921 if (insideCells[i] != -1)
1922 {
1923 insideRegions[i] = cellRegion[insideCells[i]];
1924 }
1925 reduce(insideRegions[i], maxOp<label>());
1926
1927 if (insideRegions[i] == -1)
1928 {
1929 // See if we can perturb a bit
1930 insideCells[i] = findCell
1931 (
1932 mesh_,
1933 mergeDistance_*vector::one, //perturbVec,
1934 locationsInMesh[i]+mergeDistance_*vector::one
1935 );
1936 if (insideCells[i] != -1)
1937 {
1938 insideRegions[i] = cellRegion[insideCells[i]];
1939 }
1940 reduce(insideRegions[i], maxOp<label>());
1941
1942 if (insideRegions[i] == -1)
1943 {
1944 FatalErrorInFunction
1945 << "Cannot find locationInMesh " << locationsInMesh[i]
1946 << " on any processor" << exit(FatalError);
1947 }
1948 }
1949 }
1950
1951
1952 // Check that all the locations outside the
1953 // mesh do not conflict with those inside
1954
1955 bool haveLeak = false;
1956 forAll(locationsOutsideMesh, i)
1957 {
1958 // Find the region containing the point
1959 label regioni = findRegion
1960 (
1961 mesh_,
1962 cellRegion,
1963 mergeDistance_*vector::one, //perturbVec,
1964 locationsOutsideMesh[i]
1965 );
1966
1967 if (regioni != -1)
1968 {
1969 // Check for locationsOutsideMesh overlapping with inside ones
1970 if (insideRegions.find(regioni) != -1)
1971 {
1972 haveLeak = true;
1973 WarningInFunction
1974 << "Outside location " << locationsOutsideMesh[i]
1975 << " in region " << regioni
1976 << " is connected to one of the inside points "
1977 << locationsInMesh << endl;
1978 }
1979 }
1980 }
1981
1982
1983 autoPtr<mapPolyMesh> mapPtr;
1984 if (returnReduceOr(haveLeak))
1985 {
1986 // Use holeToFace to provide a minimum set of faces needed
1987 // to close hole.
1988
1989 const List<pointField> allLocations
1990 (
1991 refinementParameters::zonePoints
1992 (
1993 locationsInMesh,
1994 zonesInMesh,
1995 locationsOutsideMesh
1996 )
1997 );
1998
1999 const labelList blockingFaces(findIndices(isBlockedFace, true));
2000
2001 labelList closureFaces;
2002 labelList closureToBlocked;
2003 autoPtr<mapDistribute> closureMapPtr;
2004 {
2005 const globalIndex globalBlockingFaces(blockingFaces.size());
2006
2007 closureMapPtr = holeToFace::calcClosure
2008 (
2009 mesh_,
2010 allLocations,
2011 blockingFaces,
2012 globalBlockingFaces,
2013 true, // allow erosion
2014
2015 closureFaces,
2016 closureToBlocked
2017 );
2018
2019 if (debug)
2020 {
2021 Pout<< "meshRefinement::blockLeakFaces :"
2022 << " found closure faces:" << closureFaces.size()
2023 << " map:" << bool(closureMapPtr) << endl;
2024 }
2025
2026 if (!closureMapPtr)
2027 {
2028 FatalErrorInFunction
2029 << "have leak but did not find any closure faces"
2030 << exit(FatalError);
2031 }
2032 }
2033
2034 // Baffle faces
2035 labelList ownPatch(mesh_.nFaces(), -1);
2036 labelList neiPatch(mesh_.nFaces(), -1);
2037
2038 // Keep (global) boundary faces in their patch
2039 {
2040 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2041 for (label patchi = 0; patchi < patches.nNonProcessor(); ++patchi)
2042 {
2043 const polyPatch& pp = patches[patchi];
2044
2045 forAll(pp, i)
2046 {
2047 ownPatch[pp.start()+i] = patchi;
2048 neiPatch[pp.start()+i] = patchi;
2049 }
2050 }
2051 }
2052
2053 const faceZoneMesh& fzs = mesh_.faceZones();
2054
2055 // Mark zone per face
2056 labelList faceToZone(mesh_.nFaces(), -1);
2057 boolList faceToFlip(mesh_.nFaces(), false);
2058 forAll(fzs, zonei)
2059 {
2060 const labelList& addressing = fzs[zonei];
2061 const boolList& flipMap = fzs[zonei].flipMap();
2062
2063 forAll(addressing, i)
2064 {
2065 faceToZone[addressing[i]] = zonei;
2066 faceToFlip[addressing[i]] = flipMap[i];
2067 }
2068 }
2069
2070
2071 // Fetch patch and zone info from blockingFaces
2072 labelList packedOwnPatch(labelUIndList(ownPatch, blockingFaces));
2073 closureMapPtr->distribute(packedOwnPatch);
2074 labelList packedNeiPatch(labelUIndList(neiPatch, blockingFaces));
2075 closureMapPtr->distribute(packedNeiPatch);
2076 labelList packedZone(labelUIndList(faceToZone, blockingFaces));
2077 closureMapPtr->distribute(packedZone);
2078 boolList packedFlip(UIndirectList<bool>(faceToFlip, blockingFaces));
2079 closureMapPtr->distribute(packedFlip);
2080 forAll(closureFaces, i)
2081 {
2082 const label facei = closureFaces[i];
2083 const label sloti = closureToBlocked[i];
2084
2085 if (sloti != -1)
2086 {
2087 // TBD. how to orient own/nei patch?
2088 ownPatch[facei] = packedOwnPatch[sloti];
2089 neiPatch[facei] = packedNeiPatch[sloti];
2090 faceToZone[facei] = packedZone[sloti];
2091 faceToFlip[facei] = packedFlip[sloti];
2092 }
2093 }
2094
2095
2096 // Add faces to faceZone. For now do this outside of createBaffles
2097 // below
2098 {
2099 List<DynamicList<label>> zoneToFaces(fzs.size());
2100 List<DynamicList<bool>> zoneToFlip(fzs.size());
2101
2102 // Add any to-be-patched face
2103 forAll(faceToZone, facei)
2104 {
2105 const label zonei = faceToZone[facei];
2106 if (zonei != -1)
2107 {
2108 zoneToFaces[zonei].append(facei);
2109 zoneToFlip[zonei].append(faceToFlip[facei]);
2110 }
2111 }
2112
2113 forAll(zoneToFaces, zonei)
2114 {
2115 surfaceZonesInfo::addFaceZone
2116 (
2117 fzs[zonei].name(),
2118 zoneToFaces[zonei],
2119 zoneToFlip[zonei],
2120 mesh_
2121 );
2122 }
2123 }
2124
2125 // Put the points of closureFaces into a special pointZone
2126 {
2127 const uindirectPrimitivePatch pp
2128 (
2129 UIndirectList<face>(mesh_.faces(), closureFaces),
2130 mesh_.points()
2131 );
2132
2133 // Count number of faces per edge
2134 const labelList nEdgeFaces(countEdgeFaces(pp));
2135
2136 // Freeze all internal points
2137 bitSet isFrozenPoint(mesh_.nPoints());
2138 forAll(nEdgeFaces, edgei)
2139 {
2140 if (nEdgeFaces[edgei] != 1)
2141 {
2142 const edge& e = pp.edges()[edgei];
2143 isFrozenPoint.set(pp.meshPoints()[e[0]]);
2144 isFrozenPoint.set(pp.meshPoints()[e[1]]);
2145 }
2146 }
2147
2148 // Lookup/add pointZone and include its points
2149 pointZoneMesh& pointZones =
2150 const_cast<pointZoneMesh&>(mesh_.pointZones());
2151 const label zonei = addPointZone("frozenPoints");
2152 const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
2153 isFrozenPoint.set(oldSet);
2154
2155 syncTools::syncPointList
2156 (
2157 mesh_,
2158 isFrozenPoint,
2159 orEqOp<unsigned int>(),
2160 0u
2161 );
2162
2163 // Override addressing
2164 pointZones.clearAddressing();
2165 pointZones[zonei] = isFrozenPoint.sortedToc();
2166 }
2167
2168
2169
2170 // Create baffles for faces
2171 mapPtr = createBaffles(ownPatch, neiPatch);
2172
2174 //{
2175 // // Add newly exposed faces (if not yet in any faceZone!)
2176 // const labelList exposed
2177 // (
2178 // renumber
2179 // (
2180 // mapPtr().reverseFaceMap(),
2181 // blockingFaces
2182 // )
2183 // );
2184 //
2185 // surfaceZonesInfo::addFaceZone
2186 // (
2187 // "frozenFaces",
2188 // exposed,
2189 // boolList(exposed.size(), false),
2190 // mesh_
2191 // );
2192 //}
2193
2194
2195 if (debug&meshRefinement::MESH)
2196 {
2197 const_cast<Time&>(mesh_.time())++;
2198
2199 Pout<< "Writing current mesh to time "
2200 << timeName() << endl;
2201 write
2202 (
2203 meshRefinement::debugType(debug),
2204 meshRefinement::writeType
2205 (
2206 meshRefinement::writeLevel()
2207 | meshRefinement::WRITEMESH
2208 ),
2209 mesh_.time().path()/timeName()
2210 );
2211 Pout<< "Dumped mesh in = "
2212 << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
2213 }
2214 }
2215 return mapPtr;
2216}
2217
2218
2219
2220/*
2221Foam::label Foam::meshRefinement::markProximityRefinementWave
2222(
2223 const scalar planarCos,
2224 const labelList& blockedSurfaces,
2225 const label nAllowRefine,
2226 const labelList& neiLevel,
2227 const pointField& neiCc,
2228
2229 labelList& refineCell,
2230 label& nRefine
2231) const
2232{
2233 labelListList faceZones(surfaces_.surfaces().size());
2234 {
2235 // If triSurface do additional zoning based on connectivity
2236 for (const label surfi : blockedSurfaces)
2237 {
2238 const label geomi = surfaces_.surfaces()[surfi];
2239 const searchableSurface& s = surfaces_.geometry()[geomi];
2240 if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
2241 {
2242 const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
2243 const labelListList& edFaces = surf.edgeFaces();
2244 boolList isOpenEdge(edFaces.size(), false);
2245 forAll(edFaces, i)
2246 {
2247 if (edFaces[i].size() == 1)
2248 {
2249 isOpenEdge[i] = true;
2250 }
2251 }
2252
2253 labelList faceZone;
2254 const label nZones = surf.markZones(isOpenEdge, faceZone);
2255 if (nZones > 1)
2256 {
2257 faceZones[surfi].transfer(faceZone);
2258 }
2259 }
2260 }
2261 }
2262
2263
2264 // Re-work the blockLevel of the blockedSurfaces into a length scale
2265 // and a number of cells to cross
2266 List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
2267 for (const label surfi : blockedSurfaces)
2268 {
2269 const label geomi = surfaces_.surfaces()[surfi];
2270 const searchableSurface& s = surfaces_.geometry()[geomi];
2271 const label nRegions = s.regions().size();
2272 regionToBlockSize[surfi].setSize(nRegions);
2273 for (label regioni = 0; regioni < nRegions; regioni++)
2274 {
2275 const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
2276 const label bLevel = surfaces_.blockLevel()[globalRegioni];
2277 regionToBlockSize[surfi][regioni] =
2278 meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
2279
2280 //const label mLevel = surfaces_.maxLevel()[globalRegioni];
2284 //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
2285 //{
2286 // const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
2287 //}
2288
2289 //nIters = max(nIters, (2<<(mLevel-bLevel)));
2290 }
2291 }
2292
2293 // Clever limiting of the number of iterations (= max cells in the channel)
2294 // since it has too many problematic issues, e.g. with volume refinement
2295 // and the real check uses the proper distance anyway just disable.
2296 const label nIters = mesh_.globalData().nTotalCells();
2297
2298
2299 // Collect candidate faces (i.e. intersecting any surface and
2300 // owner/neighbour not yet refined)
2301 const labelList testFaces(getRefineCandidateFaces(refineCell));
2302
2303 // Collect segments
2304 pointField start(testFaces.size());
2305 pointField end(testFaces.size());
2306 labelList minLevel(testFaces.size());
2307
2308 calcCellCellRays
2309 (
2310 neiCc,
2311 neiLevel,
2312 testFaces,
2313 start,
2314 end,
2315 minLevel
2316 );
2317 // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
2318 // might add to cell level). Unfortunately we only have minLevel,
2319 // not maxLevel. Problem: what if volume refinement only in middle of
2320 // channel? Say channel 1m wide with a 0.1m sphere of refinement
2321 // Workaround: have dummy surface with e.g. maxLevel 100 to
2322 // force nIters to be high enough.
2323
2324
2325 // Test for all intersections (with surfaces of higher gap level than
2326 // minLevel) and cache per cell the max surface level and the local normal
2327 // on that surface.
2328
2329 labelList surface1;
2330 List<pointIndexHit> hit1;
2331 labelList region1;
2332 vectorField normal1;
2333
2334 labelList surface2;
2335 List<pointIndexHit> hit2;
2336 labelList region2;
2337 vectorField normal2;
2338
2339 surfaces_.findNearestIntersection
2340 (
2341 blockedSurfaces,
2342 start,
2343 end,
2344
2345 surface1,
2346 hit1,
2347 region1, // local region
2348 normal1,
2349
2350 surface2,
2351 hit2,
2352 region2, // local region
2353 normal2
2354 );
2355
2356
2357 // Detect cells that are using multiple surface regions
2358 //bitSet isMultiRegion;
2359 //detectMultiRegionCells
2360 //(
2361 // faceZones,
2362 // testFaces,
2363 //
2364 // surface1,
2365 // hit1,
2366 // region1,
2367 //
2368 // surface2,
2369 // hit2,
2370 // region2,
2371 //
2372 // isMultiRegion
2373 //);
2374
2375
2376 label n = 0;
2377 forAll(testFaces, i)
2378 {
2379 if (hit1[i].hit())
2380 {
2381 n++;
2382 }
2383 }
2384
2385 List<wallPoints> faceDist(n);
2386 labelList changedFaces(n);
2387 n = 0;
2388
2389 DynamicList<point> originLocation(2);
2390 DynamicList<scalar> originDistSqr(2);
2391 DynamicList<FixedList<label, 3>> originSurface(2);
2392 //DynamicList<point> originNormal(2);
2393
2394
2395 //- To avoid walking through surfaces we mark all faces that have been
2396 // intersected. We can either mark only those faces intersecting
2397 // blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
2398 // any (refinement) surface (this includes e.g. faceZones). This is
2399 // much easier since that information is already cached
2400 // (meshRefinement::intersectedFaces())
2401
2402 //bitSet isBlockedFace(mesh_.nFaces());
2403 forAll(testFaces, i)
2404 {
2405 if (hit1[i].hit())
2406 {
2407 const label facei = testFaces[i];
2408 //isBlockedFace.set(facei);
2409 const point& fc = mesh_.faceCentres()[facei];
2410 const labelList& fz1 = faceZones[surface1[i]];
2411
2412 originLocation.clear();
2413 originDistSqr.clear();
2414 //originNormal.clear();
2415 originSurface.clear();
2416
2417 originLocation.append(hit1[i].point());
2418 originDistSqr.append(hit1[i].point().distSqr(fc));
2419 //originNormal.append(normal1[i]);
2420 originSurface.append
2421 (
2422 FixedList<label, 3>
2423 ({
2424 surface1[i],
2425 region1[i],
2426 (fz1.size() ? fz1[hit1[i].index()] : region1[i])
2427 })
2428 );
2429
2430 if (hit2[i].hit() && hit1[i] != hit2[i])
2431 {
2432 const labelList& fz2 = faceZones[surface2[i]];
2433 originLocation.append(hit2[i].point());
2434 originDistSqr.append(hit2[i].point().distSqr(fc));
2435 //originNormal.append(normal2[i]);
2436 originSurface.append
2437 (
2438 FixedList<label, 3>
2439 ({
2440 surface2[i],
2441 region2[i],
2442 (fz2.size() ? fz2[hit2[i].index()] : region2[i])
2443 })
2444 );
2445 }
2446
2447 // Collect all seed data. Currently walking does not look at
2448 // surface direction - if so pass in surface normal as well
2449 faceDist[n] = wallPoints
2450 (
2451 originLocation, // origin
2452 originDistSqr, // distance to origin
2453 originSurface // surface+region+zone
2454 //originNormal // normal at origin
2455 );
2456 changedFaces[n] = facei;
2457 n++;
2458 }
2459 }
2460
2461
2462 // Clear intersection info
2463 surface1.clear();
2464 hit1.clear();
2465 region1.clear();
2466 normal1.clear();
2467 surface2.clear();
2468 hit2.clear();
2469 region2.clear();
2470 normal2.clear();
2471
2472
2473 List<wallPoints> allFaceInfo(mesh_.nFaces());
2474 List<wallPoints> allCellInfo(mesh_.nCells());
2475
2476 // Any refinement surface (even a faceZone) should stop the gap walking.
2477 // This is exactly the information which is cached in the surfaceIndex_
2478 // field.
2479 const bitSet isBlockedFace(intersectedFaces());
2480
2481 wallPoints::trackData td(isBlockedFace, regionToBlockSize);
2482 FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
2483 (
2484 mesh_,
2485 changedFaces,
2486 faceDist,
2487 allFaceInfo,
2488 allCellInfo,
2489 0, // max iterations
2490 td
2491 );
2492 wallDistCalc.iterate(nIters);
2493
2494
2495 if (debug&meshRefinement::MESH)
2496 {
2497 // Dump current nearest opposite surfaces
2498 volScalarField distance
2499 (
2500 IOobject
2501 (
2502 "gapSize",
2503 mesh_.time().timeName(),
2504 mesh_,
2505 IOobject::NO_READ,
2506 IOobject::NO_WRITE,
2507 IOobject::NO_REGISTER
2508 ),
2509 mesh_,
2510 dimensionedScalar
2511 (
2512 "zero",
2513 dimLength, //dimensionSet(0, 1, 0, 0, 0),
2514 -1
2515 )
2516 );
2517
2518 forAll(allCellInfo, celli)
2519 {
2520 if (allCellInfo[celli].valid(wallDistCalc.data()))
2521 {
2522 const point& cc = mesh_.cellCentres()[celli];
2523 // Nearest surface points
2524 const List<point>& origin = allCellInfo[celli].origin();
2525
2526 // Find 'opposite' pair with minimum distance
2527 for (label i = 0; i < origin.size(); i++)
2528 {
2529 for (label j = i + 1; j < origin.size(); j++)
2530 {
2531 if (((cc-origin[i]) & (cc-origin[j])) < 0)
2532 {
2533 const scalar d(mag(origin[i]-origin[j]));
2534 if (distance[celli] < 0)
2535 {
2536 distance[celli] = d;
2537 }
2538 else
2539 {
2540 distance[celli] = min(distance[celli], d);
2541 }
2542 }
2543 }
2544 }
2545 }
2546 }
2547 distance.correctBoundaryConditions();
2548
2549 Info<< "Writing measured gap distance to "
2550 << distance.name() << endl;
2551 distance.write();
2552 }
2553
2554
2555
2556 // Detect tight gaps:
2557 // - cell is inbetween the two surfaces
2558 // - two surfaces are planarish
2559 // - two surfaces are not too far apart
2560 // (number of walking iterations is a too-coarse measure)
2561
2562 scalarField smallGapDistance(mesh_.nCells(), 0.0);
2563 label nMulti = 0;
2564 label nSmallGap = 0;
2565
2566 //OBJstream str(mesh_.time().timePath()/"multiRegion.obj");
2567
2568
2569 forAll(allCellInfo, celli)
2570 {
2571 if (allCellInfo[celli].valid(wallDistCalc.data()))
2572 {
2573 const point& cc = mesh_.cellCentres()[celli];
2574
2575 const List<point>& origin = allCellInfo[celli].origin();
2576 const List<FixedList<label, 3>>& surface =
2577 allCellInfo[celli].surface();
2578
2579 // Find pair with minimum distance
2580 for (label i = 0; i < origin.size(); i++)
2581 {
2582 for (label j = i + 1; j < origin.size(); j++)
2583 {
2584 //if (isMultiRegion[celli])
2585 //{
2586 // // The intersection locations are too inaccurate
2587 // // (since not proper nearest, just a cell-cell ray
2588 // // intersection) so include always
2589 //
2590 // smallGapDistance[celli] =
2591 // max(smallGapDistance[celli], maxDist);
2592 //
2593 //
2594 // str.writeLine(cc, origin[i]);
2595 // str.writeLine(cc, origin[j]);
2596 //
2597 // nMulti++;
2598 //}
2599 //else
2600 if (((cc-origin[i]) & (cc-origin[j])) < 0)
2601 {
2602 const label surfi = surface[i][0];
2603 const label regioni = surface[i][1];
2604
2605 const label surfj = surface[j][0];
2606 const label regionj = surface[j][1];
2607
2608 const scalar maxSize = max
2609 (
2610 regionToBlockSize[surfi][regioni],
2611 regionToBlockSize[surfj][regionj]
2612 );
2613
2614 if
2615 (
2616 magSqr(origin[i]-origin[j])
2617 < Foam::sqr(2*maxSize)
2618 )
2619 {
2620 const scalar maxDist
2621 (
2622 max
2623 (
2624 mag(cc-origin[i]),
2625 mag(cc-origin[j])
2626 )
2627 );
2628
2629 smallGapDistance[celli] =
2630 max(smallGapDistance[celli], maxDist);
2631 nSmallGap++;
2632 }
2633 }
2634 }
2635 }
2636 }
2637 }
2638
2639
2640 if (debug)
2641 {
2642 Info<< "Marked for blocking due to intersecting multiple surfaces : "
2643 << returnReduce(nMulti, sumOp<label>()) << " cells." << endl;
2644 Info<< "Marked for blocking due to close opposite surfaces : "
2645 << returnReduce(nSmallGap, sumOp<label>()) << " cells." << endl;
2646 }
2647
2648 if (debug&meshRefinement::MESH)
2649 {
2650 volScalarField distance
2651 (
2652 IOobject
2653 (
2654 "smallGapDistance",
2655 mesh_.time().timeName(),
2656 mesh_,
2657 IOobject::NO_READ,
2658 IOobject::NO_WRITE,
2659 IOobject::NO_REGISTER
2660 ),
2661 mesh_,
2662 dimensionedScalar(dimLength, Zero)
2663 );
2664 distance.field() = smallGapDistance;
2665 distance.correctBoundaryConditions();
2666
2667 Info<< "Writing all small-gap cells to "
2668 << distance.name() << endl;
2669 distance.write();
2670 }
2671
2672
2673 // Mark refinement
2674 const label oldNRefine = nRefine;
2675 forAll(smallGapDistance, celli)
2676 {
2677 if (smallGapDistance[celli] > SMALL)
2678 {
2679 if
2680 (
2681 !markForRefine
2682 (
2683 0, // mark level
2684 nAllowRefine,
2685 refineCell[celli],
2686 nRefine
2687 )
2688 )
2689 {
2690 if (debug)
2691 {
2692 Pout<< "Stopped refining since reaching my cell"
2693 << " limit of " << mesh_.nCells()+7*nRefine
2694 << endl;
2695 }
2696 break;
2697 }
2698 }
2699 }
2700
2701 if
2702 (
2703 returnReduce(nRefine, sumOp<label>())
2704 > returnReduce(nAllowRefine, sumOp<label>())
2705 )
2706 {
2707 Info<< "Reached refinement limit." << endl;
2708 }
2709
2710 return returnReduce(nRefine-oldNRefine, sumOp<label>());
2711}
2712*/
2713
2714// ************************************************************************* //
n
label n
Definition TABSMDCalcMethod2.H:31
pp
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
Foam::DynamicField
Dynamically sized Field. Similar to DynamicList, but inheriting from a Field instead of a List.
Definition DynamicField.H:62
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition DynamicList.H:68
Foam::FaceCellWave
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition FaceCellWave.H:209
Foam::FixedList
A 1D vector of objects of type <T> with a fixed length <N>.
Definition FixedList.H:73
Foam::HashTable::sortedToc
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition HashTable.C:156
Foam::HashTable::insert
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition HashTableI.H:152
Foam::IOobjectOption::NO_REGISTER
@ NO_REGISTER
Do not request registration (bool: false).
Definition IOobjectOption.H:99
Foam::IOobjectOption::NO_READ
@ NO_READ
Nothing to be read.
Definition IOobjectOption.H:60
Foam::IOobjectOption::NO_WRITE
@ NO_WRITE
Ignore writing from objectRegistry::writeObject().
Definition IOobjectOption.H:86
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition IOobject.H:191
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition List.H:72
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition List.H:497
Foam::List::setSize
void setSize(label n)
Alias for resize().
Definition List.H:536
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition ListI.H:133
Foam::Map
A HashTable to objects of type <T> with a label key.
Definition Map.H:54
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam::PackedList::setSize
void setSize(const label n, unsigned int val=0u)
Alias for resize().
Definition PackedList.H:819
Foam::PrimitivePatch::edgeFaces
const labelListList & edgeFaces() const
Return edge-face addressing.
Definition PrimitivePatch.C:183
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition PtrList.H:67
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition Time.H:75
Foam::UIndirectList
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition IndirectList.H:76
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition UList.H:118
Foam::UList::find
label find(const T &val) const
Find index of the first occurrence of the value.
Definition UList.C:160
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run.
Definition UPstream.H:1697
Foam::UPtrList::size
label size() const noexcept
The number of entries in the list.
Definition UPtrListI.H:106
Foam::VectorSpace::one
static const Form one
Definition VectorSpace.H:129
Foam::VectorSpace::uniform
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition VectorSpaceI.H:164
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition Vector.H:61
Foam::ZoneMesh::clearAddressing
void clearAddressing()
Clear addressing.
Definition ZoneMesh.C:944
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition autoPtr.H:65
Foam::autoPtr::clear
void clear() noexcept
Same as reset(nullptr).
Definition autoPtr.H:255
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition bitSet.H:61
Foam::bitSet::sortedToc
labelList sortedToc() const
The indices of the on bits as a sorted labelList.
Definition bitSetI.H:441
Foam::bitSet::set
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition bitSetI.H:502
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition cell.H:56
Foam::edge
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
Definition edge.H:62
Foam::faceZone
A subset of mesh faces organised as a primitive patch.
Definition faceZone.H:63
Foam::globalIndex
Calculates a non-overlapping list of offsets based on an input size (eg, number of cells) from differ...
Definition globalIndex.H:77
Foam::holeToFace::calcClosure
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:1190
Foam::meshRefinement::blockLeakFaces
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.
Definition meshRefinementBlock.C:1888
Foam::meshRefinement::writeType
writeType
Enumeration for what to write. Used as a bit-pattern.
Definition meshRefinement.H:129
Foam::meshRefinement::countFaceDirs
label countFaceDirs(const bitSet &isOutsideFace, const label celli) const
Count number of faces on cell that are in set.
Definition meshRefinementBlock.C:229
Foam::meshRefinement::removeGapCells
autoPtr< mapPolyMesh > removeGapCells(const scalar planarAngle, const labelList &minSurfaceLevel, const labelList &globalToMasterPatch, const label growIter)
Detect gapRefinement cells and remove them.
Definition meshRefinementBlock.C:1018
Foam::meshRefinement::markOutsideFaces
void markOutsideFaces(const labelList &cellLevel, const labelList &neiLevel, const labelList &refineCell, bitSet &isOutsideFace) const
Mark faces on interface between set and rest.
Definition meshRefinementBlock.C:171
Foam::meshRefinement::doRemoveCells
autoPtr< mapPolyMesh > doRemoveCells(const labelList &cellsToRemove, const labelList &exposedFaces, const labelList &exposedPatchIDs, removeCells &cellRemover)
Remove cells. Put exposedFaces into exposedPatchIDs.
Definition meshRefinement.C:1128
Foam::meshRefinement::findCell
static label findCell(const polyMesh &, const vector &perturbVec, const point &p)
Find cell point is in. Uses optional perturbation to re-test.
Definition meshRefinement.C:2815
Foam::meshRefinement::timeName
word timeName() const
Replacement for Time::timeName() that returns oldInstance (if overwrite_).
Definition meshRefinement.C:3631
Foam::meshRefinement::debugType
debugType
Enumeration for what to debug. Used as a bit-pattern.
Definition meshRefinement.H:107
Foam::meshRefinement::countEdgeFaces
labelList countEdgeFaces(const uindirectPrimitivePatch &pp) const
Count number of faces per patch edge. Parallel consistent.
Definition meshRefinement.C:2758
Foam::meshRefinement::selectSeparatedCoupledFaces
void selectSeparatedCoupledFaces(boolList &) const
Select coupled faces that are not collocated.
Definition meshRefinement.C:2742
Foam::meshRefinement::findRegion
static label findRegion(const polyMesh &, const labelList &cellRegion, const vector &perturbVec, const point &p)
Find region point is in. Uses optional perturbation to re-test.
Definition meshRefinement.C:2836
Foam::meshRefinement::addPointZone
label addPointZone(const word &name)
Add pointZone if does not exist. Return index of zone.
Definition meshRefinement.C:2720
Foam::meshRefinement::selectIntersectedFaces
void selectIntersectedFaces(const labelList &surfaces, boolList &isBlockedFace) const
Faces currently on boundary or intersected by surface.
Definition meshRefinementBlock.C:1206
Foam::meshRefinement::growSet
void growSet(const labelList &neiLevel, const bitSet &isOutsideFace, labelList &refineCell, label &nRefine) const
Add one layer of cells to set.
Definition meshRefinementBlock.C:280
Foam::meshRefinement::createBaffles
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
Definition meshRefinementBaffles.C:485
Foam::meshRefinement::writeLevel
static writeType writeLevel()
Get/set write level.
Definition meshRefinement.C:3829
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with registered IO, a reference to the associated polyMesh,...
Definition polyBoundaryMesh.H:65
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition polyPatch.H:73
Foam::refineCell
Container with cells to refine. Refinement given as single direction.
Definition refineCell.H:53
Foam::refinementParameters::zonePoints
static List< pointField > zonePoints(const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh)
Helper: per zone (entry in zonesInMesh) the locations with.
Definition refinementParameters.C:306
Foam::regionSplit
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition regionSplit.H:140
Foam::removeCells
Given list of cells to remove, insert all the topology changes.
Definition removeCells.H:60
Foam::removeCells::getExposedFaces
labelList getExposedFaces(const bitSet &removedCell) const
Get labels of faces exposed after cells removal.
Definition removeCells.C:77
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Definition searchableSurface.H:68
Foam::searchableSurface::regions
virtual const wordList & regions() const =0
Names of regions.
Foam::surfaceZonesInfo::getUnnamedSurfaces
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName).
Definition surfaceZonesInfo.C:257
Foam::surfaceZonesInfo::addFaceZone
static label addFaceZone(const word &name, const labelList &addressing, const boolList &flipMap, polyMesh &mesh)
Definition surfaceZonesInfo.C:521
Foam::syncTools::swapBoundaryCellList
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData, const bool parRun=UPstream::parRun())
Extract and swap to obtain neighbour cell values for all boundary faces.
Definition syncToolsTemplates.C:1425
Foam::syncTools::syncFaceList
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop, const bool parRun=UPstream::parRun())
Synchronize values on all mesh faces.
Definition syncTools.H:465
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition syncToolsTemplates.C:746
Foam::triSurfaceMesh
IOoject and searching on triSurface.
Definition triSurfaceMesh.H:108
Foam::triSurface::markZones
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
Definition triSurface.C:790
Foam::volumeType
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition volumeType.H:56
Foam::volumeType::MIXED
@ MIXED
A location that is partly inside and outside.
Definition volumeType.H:67
Foam::wallPoints::trackData
Class used to pass additional data in.
Definition wallPoints.H:67
Foam::wallPoints
For use with FaceCellWave. Determines topological distance to starting faces.
Definition wallPoints.H:60
patches
const polyBoundaryMesh & patches
Definition convertProcessorPatches.H:51
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
name
auto & name
Definition getAllFaRegionOptions.H:175
cells
const cellShapeList & cells
Definition gmvOutputHeader.H:3
s
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))
allFaceInfo
List< wallPoints > allFaceInfo(mesh_.nFaces())
td
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
smallGapDistance
scalarField smallGapDistance(mesh_.nCells(), 0.0)
isBlockedFace
const bitSet isBlockedFace(intersectedFaces())
oldNRefine
const label oldNRefine
Definition meshRefinementBlock.C:2673
wallDistCalc
FaceCellWave< wallPoints, wallPoints::trackData > wallDistCalc(mesh_, changedFaces, faceDist, allFaceInfo, allCellInfo, 0, td)
nMulti
label nMulti
Definition meshRefinementBlock.C:2562
allCellInfo
List< wallPoints > allCellInfo(mesh_.nCells())
nSmallGap
label nSmallGap
Definition meshRefinementBlock.C:2563
timeName
word timeName
Definition getTimeIndex.H:3
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition messageStream.H:455
Foam::SVG::end
const char * end
Definition SVGTools.H:223
Foam::debug
Namespace for handling debugging switches.
Definition debug.C:45
Foam::refCast
Type & refCast(U &obj)
A dynamic_cast (for references) to Type reference.
Definition typeInfo.H:172
Foam::distance
scalar distance(const vector &p1, const vector &p2)
Definition curveTools.C:12
Foam::wordList
List< word > wordList
List of word.
Definition fileName.H:60
Foam::returnReduceOr
bool returnReduceOr(const bool value, const int communicator=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
Definition PstreamReduceOps.H:369
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition hashSets.C:40
Foam::pointZoneMesh
ZoneMesh< pointZone, polyMesh > pointZoneMesh
A ZoneMesh with pointZone content on a polyMesh.
Definition pointZoneMeshFwd.H:39
Foam::labelListList
List< labelList > labelListList
List of labelList.
Definition labelList.H:38
Foam::labelList
List< label > labelList
A List of labels.
Definition List.H:62
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition dimensionedSymmTensor.C:44
Foam::dimLength
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition dimensionSets.H:50
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition volFieldsFwd.H:72
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition IndirectList.H:65
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere).
Foam::faceZoneMesh
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with faceZone content on a polyMesh.
Definition faceZoneMeshFwd.H:39
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition primitiveFieldsFwd.H:46
Foam::labelMax
constexpr label labelMax
Definition label.H:55
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition dimensionedScalar.C:68
Foam::returnReduce
T returnReduce(const T &value, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Definition PstreamReduceOps.H:333
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition Ostream.H:519
Foam::degToRad
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
Definition unitConversion.H:60
Foam::cellList
List< cell > cellList
List of cell.
Definition cellListFwd.H:41
Foam::isA
const Type * isA(const U &obj)
Attempt dynamic_cast to Type.
Definition typeInfo.H:87
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::reduce
void reduce(T &value, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce).
Definition PstreamReduceOps.H:51
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition hashSets.C:26
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition primitiveFieldsFwd.H:48
Foam::direction
uint8_t direction
Definition direction.H:49
Foam::point
vector point
Point is a vector.
Definition point.H:37
Foam::boolList
List< bool > boolList
A List of bools.
Definition List.H:60
Foam::Zero
static constexpr const zero Zero
Global zero (0).
Definition zero.H:127
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::uindirectPrimitivePatch
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
Definition indirectPrimitivePatch.H:55
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition pointFieldFwd.H:38
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition dimensionedScalarFwd.H:35
Foam::vector
Vector< scalar > vector
Definition vector.H:57
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition errorManip.H:125
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition dimensionedScalar.C:258
Foam::nl
constexpr char nl
The newline '\n' character (0x0a).
Definition Ostream.H:50
write
runTime write()
e
volScalarField & e
Definition createFields.H:11
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition stdFoam.H:299
unitConversion.H
Unit conversion functions.