shellSurfaces.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2015 OpenFOAM Foundation
9 Copyright (C) 2015-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "shellSurfaces.H"
30#include "searchableSurface.H"
31#include "boundBox.H"
32#include "triSurfaceMesh.H"
33#include "refinementSurfaces.H"
34#include "searchableSurfaces.H"
35#include "orientedSurface.H"
36#include "pointIndexHit.H"
37#include "volumeType.H"
39
40// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41
42const Foam::Enum
43<
45>
46Foam::shellSurfaces::refineModeNames_
47({
48 { refineMode::INSIDE, "inside" },
49 { refineMode::OUTSIDE, "outside" },
50 { refineMode::DISTANCE, "distance" },
51});
52
53
54// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
55
56void Foam::shellSurfaces::setAndCheckLevels
57(
58 const label shellI,
59 const List<Tuple2<scalar, label>>& distLevels
60)
61{
62 const searchableSurface& shell = allGeometry_[shells_[shellI]];
63
64 if (modes_[shellI] != DISTANCE && distLevels.size() != 1)
65 {
67 << "For refinement mode "
68 << refineModeNames_[modes_[shellI]]
69 << " specify only one distance+level."
70 << " (its distance gets discarded)"
71 << exit(FatalError);
72 }
73 // Extract information into separate distance and level
74 distances_[shellI].setSize(distLevels.size());
75 levels_[shellI].setSize(distLevels.size());
76
77 forAll(distLevels, j)
78 {
79 distances_[shellI][j] = distLevels[j].first();
80 levels_[shellI][j] = distLevels[j].second();
81
82 if (levels_[shellI][j] < -1)
83 {
85 << "Shell " << shell.name()
86 << " has illegal refinement level "
87 << levels_[shellI][j]
88 << exit(FatalError);
89 }
90
91
92 // Check in incremental order
93 if (j > 0)
94 {
95 if
96 (
97 (distances_[shellI][j] <= distances_[shellI][j-1])
98 || (levels_[shellI][j] > levels_[shellI][j-1])
99 )
100 {
102 << "For refinement mode "
103 << refineModeNames_[modes_[shellI]]
104 << " : Refinement should be specified in order"
105 << " of increasing distance"
106 << " (and decreasing refinement level)." << endl
107 << "Distance:" << distances_[shellI][j]
108 << " refinementLevel:" << levels_[shellI][j]
109 << exit(FatalError);
110 }
111 }
112 }
113
114 if (modes_[shellI] == DISTANCE)
115 {
116 if (!dryRun_)
117 {
118 Info<< "Refinement level according to distance to "
119 << shell.name() << endl;
120 forAll(levels_[shellI], j)
121 {
122 Info<< " level " << levels_[shellI][j]
123 << " for all cells within " << distances_[shellI][j]
124 << " metre." << endl;
125 }
126 }
127 }
128 else
129 {
130 if (!shell.hasVolumeType())
131 {
133 << "Shell " << shell.name()
134 << " does not support testing for "
135 << refineModeNames_[modes_[shellI]] << endl
136 << "Probably it is not closed."
137 << exit(FatalError);
138 }
139
140 if (!dryRun_)
141 {
142 if (modes_[shellI] == INSIDE)
143 {
144 Info<< "Refinement level " << levels_[shellI][0]
145 << " for all cells inside " << shell.name() << endl;
146 }
147 else
148 {
149 Info<< "Refinement level " << levels_[shellI][0]
150 << " for all cells outside " << shell.name() << endl;
151 }
152 }
153 }
154}
155
156
157// Specifically orient triSurfaces using a calculated point outside.
158// Done since quite often triSurfaces not of consistent orientation which
159// is (currently) necessary for sideness calculation
160void Foam::shellSurfaces::orient()
161{
162 // Determine outside point.
163 boundBox overallBb = boundBox::invertedBox;
164
165 bool hasSurface = false;
166
167 forAll(shells_, shellI)
168 {
169 const searchableSurface& s = allGeometry_[shells_[shellI]];
170
171 if
172 (
173 modes_[shellI] != DISTANCE
174 && isA<triSurfaceMesh>(s)
175 && !isA<distributedTriSurfaceMesh>(s)
176 )
177 {
178 hasSurface = true;
179
180 const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
181 if (shell.triSurface::size())
182 {
183 // Assume surface is compact!
184 overallBb.add(s.bounds());
185 }
186 }
187 }
188
189 if (returnReduce(hasSurface, orOp<bool>()))
190 {
191 const point outsidePt = overallBb.max() + overallBb.span();
192
193 //Info<< "Using point " << outsidePt << " to orient shells" << endl;
194
195 forAll(shells_, shellI)
196 {
197 const searchableSurface& s = allGeometry_[shells_[shellI]];
198
199 if
200 (
201 modes_[shellI] != DISTANCE
202 && isA<triSurfaceMesh>(s)
203 && !isA<distributedTriSurfaceMesh>(s)
204 )
205 {
206 triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
207 (
208 refCast<const triSurfaceMesh>(s)
209 );
210
211 // Flip surface so outsidePt is outside.
212 bool anyFlipped = orientedSurface::orient
213 (
214 shell,
215 outsidePt,
216 true
217 );
218
219 if (anyFlipped && !dryRun_)
220 {
221 // orientedSurface will have done a clearOut of the surface.
222 // we could do a clearout of the triSurfaceMeshes::trees()
223 // but these aren't affected by orientation
224 // (except for cached
225 // sideness which should not be set at this point.
226 // !!Should check!)
227
228 Info<< "shellSurfaces : Flipped orientation of surface "
229 << s.name()
230 << " so point " << outsidePt << " is outside." << endl;
231 }
232 }
233 }
234 }
235}
236
237
238// Find maximum level of a shell.
239void Foam::shellSurfaces::findHigherLevel
240(
241 const pointField& pt,
242 const label shellI,
243 labelList& maxLevel
244) const
245{
246 const labelList& levels = levels_[shellI];
247
248 if (modes_[shellI] == DISTANCE)
249 {
250 // Distance mode.
251
252 const scalarField& distances = distances_[shellI];
253
254 // Collect all those points that have a current maxLevel less than
255 // (any of) the shell. Also collect the furthest distance allowable
256 // to any shell with a higher level.
257
258 labelList candidateMap(pt.size());
259 scalarField candidateDistSqr(pt.size());
260 label candidateI = 0;
261
262 forAll(maxLevel, pointi)
263 {
264 forAllReverse(levels, levelI)
265 {
266 if (levels[levelI] > maxLevel[pointi])
267 {
268 candidateMap[candidateI] = pointi;
269 candidateDistSqr[candidateI] = sqr(distances[levelI]);
270 candidateI++;
271 break;
272 }
273 }
274 }
275 candidateMap.setSize(candidateI);
276 candidateDistSqr.setSize(candidateI);
277
278 // Do the expensive nearest test only for the candidate points.
279 List<pointIndexHit> nearInfo;
280 allGeometry_[shells_[shellI]].findNearest
281 (
282 pointField(pt, candidateMap),
283 candidateDistSqr,
284 nearInfo
285 );
286
287 // Update maxLevel
288 forAll(nearInfo, i)
289 {
290 if (nearInfo[i].hit())
291 {
292 const label pointi = candidateMap[i];
293
294 // Check which level it actually is in.
295 const label minDistI = findLower
296 (
297 distances,
298 mag(nearInfo[i].hitPoint()-pt[pointi])
299 );
300
301 // pt is inbetween shell[minDistI] and shell[minDistI+1]
302 maxLevel[pointi] = levels[minDistI+1];
303 }
304 }
305 }
306 else
307 {
308 // Inside/outside mode
309
310 // Collect all those points that have a current maxLevel less than the
311 // shell.
312
313 pointField candidates(pt.size());
314 labelList candidateMap(pt.size());
315 label candidateI = 0;
316
317 forAll(maxLevel, pointi)
318 {
319 if (levels[0] > maxLevel[pointi])
320 {
321 candidates[candidateI] = pt[pointi];
322 candidateMap[candidateI] = pointi;
323 candidateI++;
324 }
325 }
326 candidates.setSize(candidateI);
327 candidateMap.setSize(candidateI);
328
329 // Do the expensive nearest test only for the candidate points.
330 List<volumeType> volType;
331 allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
332
333 forAll(volType, i)
334 {
335 label pointi = candidateMap[i];
336
337 if
338 (
339 (
340 modes_[shellI] == INSIDE
341 && volType[i] == volumeType::INSIDE
342 )
343 || (
344 modes_[shellI] == OUTSIDE
345 && volType[i] == volumeType::OUTSIDE
346 )
347 )
348 {
349 maxLevel[pointi] = levels[0];
350 }
351 }
352 }
353}
354
355
356void Foam::shellSurfaces::findHigherGapLevel
357(
358 const pointField& pt,
359 const labelList& ptLevel,
360 const label shellI,
361 labelList& gapShell,
362 List<FixedList<label, 3>>& gapInfo,
363 List<volumeType>& gapMode
364) const
365{
366 //TBD: hardcoded for region 0 information
367 const FixedList<label, 3>& info = extendedGapLevel_[shellI][0];
368 const volumeType mode = extendedGapMode_[shellI][0];
369
370 if (info[2] == 0)
371 {
372 return;
373 }
374
375 if (modes_[shellI] == DISTANCE)
376 {
377 // Distance mode.
378 const scalar distance = distances_[shellI][0];
379
380 labelList candidateMap(pt.size());
381 scalarField candidateDistSqr(pt.size());
382 label candidateI = 0;
383
384 forAll(ptLevel, pointI)
385 {
386 if (ptLevel[pointI] >= info[1] && ptLevel[pointI] < info[2])
387 {
388 candidateMap[candidateI] = pointI;
389 candidateDistSqr[candidateI] = sqr(distance);
390 candidateI++;
391 }
392 }
393 candidateMap.setSize(candidateI);
394 candidateDistSqr.setSize(candidateI);
395
396 // Do the expensive nearest test only for the candidate points.
397 List<pointIndexHit> nearInfo;
398 allGeometry_[shells_[shellI]].findNearest
399 (
400 pointField(pt, candidateMap),
401 candidateDistSqr,
402 nearInfo
403 );
404
405 // Update info
406 forAll(nearInfo, i)
407 {
408 if (nearInfo[i].hit())
409 {
410 const label pointI = candidateMap[i];
411 gapShell[pointI] = shellI;
412 gapInfo[pointI] = info;
413 gapMode[pointI] = mode;
414 }
415 }
416 }
417 else
418 {
419 // Collect all those points that have a current maxLevel less than the
420 // shell.
421
422 labelList candidateMap(pt.size());
423 label candidateI = 0;
424
425 forAll(ptLevel, pointI)
426 {
427 if (ptLevel[pointI] >= info[1] && ptLevel[pointI] < info[2])
428 {
429 candidateMap[candidateI++] = pointI;
430 }
431 }
432 candidateMap.setSize(candidateI);
433
434 // Do the expensive nearest test only for the candidate points.
435 List<volumeType> volType;
436 allGeometry_[shells_[shellI]].getVolumeType
437 (
438 pointField(pt, candidateMap),
439 volType
440 );
441
442 forAll(volType, i)
443 {
444 const label pointI = candidateMap[i];
445 const bool isInside = (volType[i] == volumeType::INSIDE);
446
447 if
448 (
449 (
450 (modes_[shellI] == INSIDE && isInside)
451 || (modes_[shellI] == OUTSIDE && !isInside)
452 )
453 && info[2] > gapInfo[pointI][2]
454 )
455 {
456 gapShell[pointI] = shellI;
457 gapInfo[pointI] = info;
458 gapMode[pointI] = mode;
459 }
460 }
461 }
462}
463
464
465void Foam::shellSurfaces::findLevel
466(
467 const pointField& pt,
468 const label shellI,
469 labelList& minLevel,
470 labelList& shell
471) const
472{
473 const labelList& levels = levels_[shellI];
474
475 if (modes_[shellI] == DISTANCE)
476 {
477 // Distance mode.
478
479 const scalarField& distances = distances_[shellI];
480
481 // Collect all those points that have a current level equal/greater
482 // (any of) the shell. Also collect the furthest distance allowable
483 // to any shell with a higher level.
484
485 pointField candidates(pt.size());
486 labelList candidateMap(pt.size());
487 scalarField candidateDistSqr(pt.size());
488 label candidateI = 0;
489
490 forAll(shell, pointI)
491 {
492 if (shell[pointI] == -1)
493 {
494 forAllReverse(levels, levelI)
495 {
496 if (levels[levelI] <= minLevel[pointI])
497 {
498 candidates[candidateI] = pt[pointI];
499 candidateMap[candidateI] = pointI;
500 candidateDistSqr[candidateI] = sqr(distances[levelI]);
501 candidateI++;
502 break;
503 }
504 }
505 }
506 }
507 candidates.setSize(candidateI);
508 candidateMap.setSize(candidateI);
509 candidateDistSqr.setSize(candidateI);
510
511 // Do the expensive nearest test only for the candidate points.
512 List<pointIndexHit> nearInfo;
513 allGeometry_[shells_[shellI]].findNearest
514 (
515 candidates,
516 candidateDistSqr,
517 nearInfo
518 );
519
520 // Update maxLevel
521 forAll(nearInfo, i)
522 {
523 if (nearInfo[i].hit())
524 {
525 // Check which level it actually is in.
526 label minDistI = findLower
527 (
528 distances,
529 mag(nearInfo[i].hitPoint()-candidates[i])
530 );
531
532 label pointI = candidateMap[i];
533
534 // pt is inbetween shell[minDistI] and shell[minDistI+1]
535 shell[pointI] = shellI;
536 minLevel[pointI] = levels[minDistI+1];
537 }
538 }
539 }
540 else
541 {
542 // Inside/outside mode
543
544 // Collect all those points that have a current maxLevel less than the
545 // shell.
546
547 pointField candidates(pt.size());
548 labelList candidateMap(pt.size());
549 label candidateI = 0;
550
551 forAll(shell, pointI)
552 {
553 if (shell[pointI] == -1 && levels[0] <= minLevel[pointI])
554 {
555 candidates[candidateI] = pt[pointI];
556 candidateMap[candidateI] = pointI;
557 candidateI++;
558 }
559 }
560 candidates.setSize(candidateI);
561 candidateMap.setSize(candidateI);
562
563 // Do the expensive nearest test only for the candidate points.
564 List<volumeType> volType;
565 allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
566
567 forAll(volType, i)
568 {
569 if
570 (
571 (
572 modes_[shellI] == INSIDE
573 && volType[i] == volumeType::INSIDE
574 )
575 || (
576 modes_[shellI] == OUTSIDE
577 && volType[i] == volumeType::OUTSIDE
578 )
579 )
580 {
581 label pointI = candidateMap[i];
582 shell[pointI] = shellI;
583 minLevel[pointI] = levels[0];
584 }
585 }
586 }
587}
588
589
590// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
591
593(
594 const searchableSurfaces& allGeometry,
595 const dictionary& shellsDict,
596 const bool dryRun
597)
598:
599 allGeometry_(allGeometry),
600 dryRun_(dryRun)
601{
602 // Wildcard specification : loop over all surfaces and try to find a match.
603
604 // Count number of shells.
605 label shellI = 0;
606 for (const word& geomName : allGeometry_.names())
607 {
608 if (shellsDict.found(geomName))
609 {
610 ++shellI;
611 }
612 }
613
614
615 // Size lists
616 shells_.setSize(shellI);
617 modes_.setSize(shellI);
618 distances_.setSize(shellI);
619 levels_.setSize(shellI);
620 dirLevels_.setSize(shellI);
621 smoothDirection_.setSize(shellI);
622 nSmoothExpansion_.setSize(shellI);
623 nSmoothPosition_.setSize(shellI);
624
625 extendedGapLevel_.setSize(shellI);
626 extendedGapMode_.setSize(shellI);
627 selfProximity_.setSize(shellI);
628
629 const FixedList<label, 3> nullGapLevel({0, 0, 0});
630
631 wordHashSet unmatchedKeys(shellsDict.toc());
632 shellI = 0;
633
634 forAll(allGeometry_.names(), geomI)
635 {
636 const word& geomName = allGeometry_.names()[geomI];
637
638 const entry* eptr = shellsDict.findEntry(geomName, keyType::REGEX);
639
640 if (eptr)
641 {
642 const dictionary& dict = eptr->dict();
643 unmatchedKeys.erase(eptr->keyword());
644
645 shells_[shellI] = geomI;
646 modes_[shellI] = refineModeNames_.get("mode", dict);
647
648 // Read pairs of distance+level
649 setAndCheckLevels(shellI, dict.lookup("levels"));
650
651
652 // Directional refinement
653 // ~~~~~~~~~~~~~~~~~~~~~~
654
655 dirLevels_[shellI] = Tuple2<labelPair,labelVector>
656 (
659 );
660 const entry* levelPtr =
661 dict.findEntry("levelIncrement", keyType::REGEX);
662
663 if (levelPtr)
664 {
665 // Do reading ourselves since using labelPair would require
666 // additional bracket pair
667 ITstream& is = levelPtr->stream();
668
669 is.readBegin("levelIncrement");
670 is >> dirLevels_[shellI].first().first()
671 >> dirLevels_[shellI].first().second()
672 >> dirLevels_[shellI].second();
673 is.readEnd("levelIncrement");
674
675 if (modes_[shellI] == INSIDE)
676 {
677 if (!dryRun_)
678 {
679 Info<< "Additional directional refinement level"
680 << " for all cells inside " << geomName << endl;
681 }
682 }
683 else if (modes_[shellI] == OUTSIDE)
684 {
685 if (!dryRun_)
686 {
687 Info<< "Additional directional refinement level"
688 << " for all cells outside " << geomName << endl;
689 }
690 }
691 else
692 {
693 FatalIOErrorInFunction(shellsDict)
694 << "Unsupported mode "
695 << refineModeNames_[modes_[shellI]]
696 << exit(FatalIOError);
697 }
698 }
699
700 // Directional smoothing
701 // ~~~~~~~~~~~~~~~~~~~~~
702
703 nSmoothExpansion_[shellI] = 0;
704 nSmoothPosition_[shellI] = 0;
705 smoothDirection_[shellI] =
706 dict.getOrDefault("smoothDirection", vector::zero);
707
708 if (smoothDirection_[shellI] != vector::zero)
709 {
710 dict.readEntry("nSmoothExpansion", nSmoothExpansion_[shellI]);
711 dict.readEntry("nSmoothPosition", nSmoothPosition_[shellI]);
712 }
713
714
715 // Gap specification
716 // ~~~~~~~~~~~~~~~~~
717
718
719 // Shell-wide gap level specification
720 const searchableSurface& surface = allGeometry_[geomI];
721 const wordList& regionNames = surface.regions();
722
723 FixedList<label, 3> gapSpec
724 (
726 (
727 "gapLevel",
728 nullGapLevel
729 )
730 );
731 extendedGapLevel_[shellI].setSize(regionNames.size());
732 extendedGapLevel_[shellI] = gapSpec;
733
734 extendedGapMode_[shellI].setSize(regionNames.size());
735 extendedGapMode_[shellI] =
736 volumeType("gapMode", dict, volumeType::MIXED);
737
738 // Detect self-intersections
739 selfProximity_[shellI].setSize
740 (
742 dict.getOrDefault<bool>("gapSelf", true)
743 );
744
745
746 // Override on a per-region basis?
747
748 if (dict.found("regions"))
749 {
750 const dictionary& regionsDict = dict.subDict("regions");
751 forAll(regionNames, regionI)
752 {
753 if (regionsDict.found(regionNames[regionI]))
754 {
755 // Get the dictionary for region
756 const dictionary& regionDict = regionsDict.subDict
757 (
758 regionNames[regionI]
759 );
760 FixedList<label, 3> gapSpec
761 (
762 regionDict.getOrDefault
763 (
764 "gapLevel",
765 nullGapLevel
766 )
767 );
768 extendedGapLevel_[shellI][regionI] = gapSpec;
769
770 extendedGapMode_[shellI][regionI] =
772 (
773 "gapMode",
774 regionDict,
776 );
777
778 selfProximity_[shellI][regionI] =
779 regionDict.getOrDefault<bool>
780 (
781 "gapSelf",
782 true
783 );
784 }
785 }
786 }
787
788 if (extendedGapLevel_[shellI][0][0] > 0)
789 {
790 Info<< "Refinement level up to "
791 << extendedGapLevel_[shellI][0][2]
792 << " for all cells in gaps for shell "
793 << geomName << endl;
794
795 if (distances_[shellI].size() > 1)
796 {
797 WarningInFunction << "Using first distance only out of "
798 << distances_[shellI] << " to detect gaps for shell "
799 << geomName << endl;
800 }
801 }
802
803 shellI++;
804 }
805 }
806
807 if (unmatchedKeys.size() > 0)
808 {
809 IOWarningInFunction(shellsDict)
810 << "Not all entries in refinementRegions dictionary were used."
811 << " The following entries were not used : "
812 << unmatchedKeys.sortedToc()
813 << endl;
814 }
815
816 // Orient shell surfaces before any searching is done. Note that this
817 // only needs to be done for inside or outside. Orienting surfaces
818 // constructs lots of addressing which we want to avoid.
819 orient();
820}
821
822
823// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
824
825// Highest shell level
827{
828 label overallMax = 0;
829 forAll(levels_, shellI)
830 {
831 overallMax = max(overallMax, max(levels_[shellI]));
832 }
833 return overallMax;
834}
835
836
838{
839 labelList surfaceMax(extendedGapLevel_.size(), Zero);
840
841 forAll(extendedGapLevel_, shelli)
842 {
843 const List<FixedList<label, 3>>& levels = extendedGapLevel_[shelli];
844 forAll(levels, i)
845 {
846 surfaceMax[shelli] = max(surfaceMax[shelli], levels[i][2]);
847 }
848 }
849 return surfaceMax;
850}
851
852
854{
855 labelPairList levels(dirLevels_.size());
856 forAll(dirLevels_, shelli)
857 {
858 levels[shelli] = dirLevels_[shelli].first();
859 }
860 return levels;
861}
862
863
865{
866 return nSmoothExpansion_;
867}
868
869
871{
872 return smoothDirection_;
873}
874
875
877{
878 return nSmoothPosition_;
879}
880
881
882void Foam::shellSurfaces::findHigherLevel
883(
884 const pointField& pt,
885 const labelList& ptLevel,
886 labelList& maxLevel
887) const
888{
889 // Maximum level of any shell. Start off with level of point.
890 maxLevel = ptLevel;
891
892 forAll(shells_, shelli)
893 {
894 findHigherLevel(pt, shelli, maxLevel);
895 }
896}
897
898
899void Foam::shellSurfaces::findHigherGapLevel
900(
901 const pointField& pt,
902 const labelList& ptLevel,
903 labelList& gapShell,
904 List<FixedList<label, 3>>& gapInfo,
905 List<volumeType>& gapMode
906) const
907{
908 gapShell.setSize(pt.size());
909 gapShell = -1;
910
911 gapInfo.setSize(pt.size());
912 gapInfo = FixedList<label, 3>({0, 0, 0});
913
914 gapMode.setSize(pt.size());
915 gapMode = volumeType::MIXED;
916
917 forAll(shells_, shelli)
918 {
919 findHigherGapLevel(pt, ptLevel, shelli, gapShell, gapInfo, gapMode);
920 }
921}
922
923
924void Foam::shellSurfaces::findHigherGapLevel
925(
926 const pointField& pt,
927 const labelList& ptLevel,
928 List<FixedList<label, 3>>& gapInfo,
929 List<volumeType>& gapMode
930) const
931{
932 labelList gapShell;
933 findHigherGapLevel(pt, ptLevel, gapShell, gapInfo, gapMode);
934}
935
936
937void Foam::shellSurfaces::findLevel
938(
939 const pointField& pt,
940 const labelList& ptLevel,
941 labelList& shell
942) const
943{
944 shell.setSize(pt.size());
945 shell = -1;
946
947 labelList minLevel(ptLevel);
948
949 forAll(shells_, shelli)
950 {
951 findLevel(pt, shelli, minLevel, shell);
952 }
953}
954
955
957(
958 const pointField& pt,
959 const labelList& ptLevel,
960 const labelList& dirLevel, // directional level
961 const direction dir,
962 labelList& shell
963) const
964{
965 shell.setSize(pt.size());
966 shell = -1;
967
968 List<volumeType> volType;
969
970 // Current back to original
971 DynamicList<label> candidateMap(pt.size());
972
973 forAll(shells_, shelli)
974 {
975 if (modes_[shelli] == INSIDE || modes_[shelli] == OUTSIDE)
976 {
977 const labelPair& selectLevels = dirLevels_[shelli].first();
978 const label addLevel = dirLevels_[shelli].second()[dir];
979
980 // Collect the cells that are of the right original level
981 candidateMap.clear();
982 forAll(ptLevel, celli)
983 {
984 label level = ptLevel[celli];
985
986 if
987 (
988 level >= selectLevels.first()
989 && level <= selectLevels.second()
990 && dirLevel[celli] < level+addLevel
991 )
992 {
993 candidateMap.append(celli);
994 }
995 }
996
997 // Do geometric test
998 pointField candidatePt(pt, candidateMap);
999 allGeometry_[shells_[shelli]].getVolumeType(candidatePt, volType);
1000
1001 // Extract selected cells
1002 forAll(candidateMap, i)
1003 {
1004 if
1005 (
1006 (
1007 modes_[shelli] == INSIDE
1008 && volType[i] == volumeType::INSIDE
1009 )
1010 || (
1011 modes_[shelli] == OUTSIDE
1012 && volType[i] == volumeType::OUTSIDE
1013 )
1014 )
1015 {
1016 shell[candidateMap[i]] = shelli;
1017 }
1018 }
1019 }
1020 }
1021}
1022
1023
1024// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
EnumType get(const word &enumName) const
The enumeration corresponding to the given name.
Definition: Enum.C:75
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
static constexpr label size() noexcept
Return the number of elements in the FixedList.
Definition: FixedList.H:416
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
An input stream of tokens.
Definition: ITstream.H:56
bool readEnd(const char *funcName)
End read of data chunk, ends with ')'.
Definition: Istream.C:129
bool readBegin(const char *funcName)
Begin read of data chunk, starts with '('.
Definition: Istream.C:111
void setSize(const label n)
Alias for resize()
Definition: List.H:218
direction orient
The x/y/z orientation (0,1,2)
Definition: PDRobstacle.H:116
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:120
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
Definition: Tuple2.H:58
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static const Form max
Definition: VectorSpace.H:117
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
entry * findEntry(const word &keyword, enum keyType::option matchOpt=keyType::REGEX)
Find for an entry (non-const access) with the given keyword.
Definition: dictionaryI.H:97
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionary.C:386
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
wordList toc() const
Return the table of contents.
Definition: dictionary.C:602
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:70
virtual ITstream & stream() const =0
Return token stream, if entry is a primitive entry.
const keyType & keyword() const noexcept
Return keyword.
Definition: entry.H:195
virtual const dictionary & dict() const =0
Return dictionary, if entry is a dictionary.
@ REGEX
Regular expression.
Definition: keyType.H:82
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
const wordList & names() const
Surface names, not region names.
Encapsulates queries for volume refinement ('refine all cells within shell').
Definition: shellSurfaces.H:58
const vectorField & smoothDirection() const
Per shell the smoothing direction.
labelList maxGapLevel() const
Highest shell gap level.
void findDirectionalLevel(const pointField &pt, const labelList &ptLevel, const labelList &dirLevel, const direction dir, labelList &shell) const
Find any shell (or -1) with higher wanted directional level.
const labelList & nSmoothPosition() const
Per shell the positional smoothing iterations.
label maxLevel() const
Highest shell level.
const labelList & nSmoothExpansion() const
Per shell the directional smoothing iterations.
refineMode
Volume refinement controls.
Definition: shellSurfaces.H:65
labelPairList directionalSelectLevel() const
Min and max cell level for directional refinement.
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
wordList regionNames
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))
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
#define WarningInFunction
Report a warning using Foam::Warning.
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
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)
constexpr label labelMin
Definition: label.H:60
messageStream Info
Information stream (stdout output on master, null elsewhere)
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Definition: MSwindows.C:572
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
label findLower(const ListType &input, const T &val, const label start, const ComparePredicate &comp)
vector point
Point is a vector.
Definition: point.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
constexpr label labelMax
Definition: label.H:61
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
uint8_t direction
Definition: direction.H:56
IOerror FatalIOError
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:346
A non-counting (dummy) refCount.
Definition: refCount.H:59