refinementSurfaces.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-2020 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 "refinementSurfaces.H"
30#include "Time.H"
31#include "searchableSurfaces.H"
32#include "shellSurfaces.H"
33#include "triSurfaceMesh.H"
34#include "labelPair.H"
35#include "searchableSurfacesQueries.H"
36#include "UPtrList.H"
37#include "volumeType.H"
38// For dictionary::get wrapper
39#include "meshRefinement.H"
40
41// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42
43Foam::labelList Foam::refinementSurfaces::findHigherLevel
44(
45 const searchableSurface& geom,
46 const shellSurfaces& shells,
47 const List<pointIndexHit>& intersectionInfo,
48 const labelList& surfaceLevel // current level
49) const
50{
51 // See if a cached level field available
52 labelList minLevelField;
53 geom.getField(intersectionInfo, minLevelField);
54
55
56 // Detect any uncached values and do proper search
57 labelList localLevel(surfaceLevel);
58 {
59 // Check hits:
60 // 1. cached value == -1 : store for re-testing
61 // 2. cached value != -1 : use
62 // 3. uncached : use region 0 value
63
64 DynamicList<label> retestSet;
65 label nHits = 0;
66
67 forAll(intersectionInfo, i)
68 {
69 if (intersectionInfo[i].hit())
70 {
71 nHits++;
72
73 // Check if minLevelField for this surface.
74 if (minLevelField.size())
75 {
76 if (minLevelField[i] == -1)
77 {
78 retestSet.append(i);
79 }
80 else
81 {
82 localLevel[i] = max(localLevel[i], minLevelField[i]);
83 }
84 }
85 else
86 {
87 retestSet.append(i);
88 }
89 }
90 }
91
92 label nRetest = returnReduce(retestSet.size(), sumOp<label>());
93 if (nRetest > 0)
94 {
95 reduce(nHits, sumOp<label>());
96
97 //Info<< "Retesting " << nRetest
98 // << " out of " << nHits
99 // << " intersections on uncached elements on geometry "
100 // << geom.name() << endl;
101
102 pointField samples(retestSet.size());
103 forAll(retestSet, i)
104 {
105 samples[i] = intersectionInfo[retestSet[i]].hitPoint();
106 }
107 labelList shellLevel;
108 shells.findHigherLevel
109 (
110 samples,
111 labelUIndList(surfaceLevel, retestSet)(),
112 shellLevel
113 );
114 forAll(retestSet, i)
115 {
116 label sampleI = retestSet[i];
117 localLevel[sampleI] = max(localLevel[sampleI], shellLevel[i]);
118 }
119 }
120 }
121
122 return localLevel;
123}
124
125
126Foam::labelList Foam::refinementSurfaces::calcSurfaceIndex
127(
128 const searchableSurfaces& allGeometry,
129 const labelList& surfaces
130)
131{
132 // Determine overall number of global regions
133 label globalI = 0;
134 forAll(surfaces, surfI)
135 {
136 globalI += allGeometry[surfaces[surfI]].regions().size();
137 }
138
139 labelList regionToSurface(globalI);
140 globalI = 0;
141 forAll(surfaces, surfI)
142 {
143 const label nLocal = allGeometry[surfaces[surfI]].regions().size();
144 for (label i = 0; i < nLocal; i++)
145 {
146 regionToSurface[globalI++] = surfI;
147 }
148 }
149
150 return regionToSurface;
151}
152
153
154// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
155
156Foam::refinementSurfaces::refinementSurfaces
157(
158 const searchableSurfaces& allGeometry,
159 const dictionary& surfacesDict,
160 const label gapLevelIncrement,
161 const bool dryRun
162)
163:
164 allGeometry_(allGeometry),
165 surfaces_(surfacesDict.size()),
166 names_(surfacesDict.size()),
167 surfZones_(surfacesDict.size()),
168 regionOffset_(surfacesDict.size()),
169 dryRun_(dryRun)
170{
171 // Wildcard specification : loop over all surface, all regions
172 // and try to find a match.
173
174 // Count number of surfaces.
175 label surfI = 0;
176 forAll(allGeometry_.names(), geomI)
177 {
178 const word& geomName = allGeometry_.names()[geomI];
179
180 if (surfacesDict.found(geomName))
181 {
182 surfI++;
183 }
184 }
185
186 // Size lists
187 surfaces_.setSize(surfI);
188 names_.setSize(surfI);
189 surfZones_.setSize(surfI);
190 regionOffset_.setSize(surfI);
191
192 // Per surface data
193 labelList globalMinLevel(surfI, Zero);
194 labelList globalMaxLevel(surfI, Zero);
195 labelList globalLevelIncr(surfI, Zero);
196
197 FixedList<label, 3> nullGapLevel;
198 nullGapLevel[0] = 0;
199 nullGapLevel[1] = 0;
200 nullGapLevel[2] = 0;
201
202 List<FixedList<label, 3>> globalGapLevel(surfI);
203 List<volumeType> globalGapMode(surfI);
204 boolList globalGapSelf(surfI);
205
206 scalarField globalAngle(surfI, -GREAT);
207 PtrList<dictionary> globalPatchInfo(surfI);
208
209 labelList globalBlockLevel(surfI, labelMax);
210 labelList globalLeakLevel(surfI, labelMax);
211
212 // Per surface, per region data
213 List<Map<label>> regionMinLevel(surfI);
214 List<Map<label>> regionMaxLevel(surfI);
215 List<Map<label>> regionLevelIncr(surfI);
216 List<Map<FixedList<label, 3>>> regionGapLevel(surfI);
217 List<Map<volumeType>> regionGapMode(surfI);
218 List<Map<bool>> regionGapSelf(surfI);
219 List<Map<scalar>> regionAngle(surfI);
220 List<Map<autoPtr<dictionary>>> regionPatchInfo(surfI);
221 List<Map<label>> regionBlockLevel(surfI);
222 List<Map<label>> regionLeakLevel(surfI);
223
224 wordHashSet unmatchedKeys(surfacesDict.toc());
225
226 surfI = 0;
227 forAll(allGeometry_.names(), geomI)
228 {
229 const word& geomName = allGeometry_.names()[geomI];
230
231 const entry* ePtr =
232 surfacesDict.findEntry(geomName, keyType::REGEX);
233
234 if (ePtr)
235 {
236 const dictionary& dict = ePtr->dict();
237 unmatchedKeys.erase(ePtr->keyword());
238
239 names_[surfI] = geomName;
240 surfaces_[surfI] = geomI;
241
242 const labelPair refLevel
243 (
244 meshRefinement::get<labelPair>
245 (
246 dict,
247 "level",
248 dryRun_,
249 keyType::REGEX,
250 labelPair(0, 0)
251 )
252 );
253
254 globalMinLevel[surfI] = refLevel[0];
255 globalMaxLevel[surfI] = refLevel[1];
256 globalLevelIncr[surfI] = dict.getOrDefault
257 (
258 "gapLevelIncrement",
259 gapLevelIncrement
260 );
261
262 if
263 (
264 globalMinLevel[surfI] < 0
265 || globalMaxLevel[surfI] < globalMinLevel[surfI]
266 || globalMaxLevel[surfI] < 0
267 || globalLevelIncr[surfI] < 0
268 )
269 {
270 FatalIOErrorInFunction(dict)
271 << "Illegal level specification for surface "
272 << names_[surfI]
273 << " : minLevel:" << globalMinLevel[surfI]
274 << " maxLevel:" << globalMaxLevel[surfI]
275 << " levelIncrement:" << globalLevelIncr[surfI]
276 << exit(FatalIOError);
277 }
278
279
280 // Optional gapLevel specification
281
282 globalGapLevel[surfI] = dict.getOrDefault
283 (
284 "gapLevel",
285 nullGapLevel
286 );
287 globalGapMode[surfI] =
288 volumeType("gapMode", dict, volumeType::MIXED);
289
290 if
291 (
292 globalGapMode[surfI] == volumeType::UNKNOWN
293 || globalGapLevel[surfI][0] < 0
294 || globalGapLevel[surfI][1] < 0
295 || globalGapLevel[surfI][2] < 0
296 || globalGapLevel[surfI][1] > globalGapLevel[surfI][2]
297 )
298 {
299 FatalIOErrorInFunction(dict)
300 << "Illegal gapLevel specification for surface "
301 << names_[surfI]
302 << " : gapLevel:" << globalGapLevel[surfI]
303 << " gapMode:" << globalGapMode[surfI].str()
304 << exit(FatalIOError);
305 }
306
307 globalGapSelf[surfI] =
308 dict.getOrDefault<bool>("gapSelf", true);
309
310 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
311
312 // Surface zones
313 surfZones_.set
314 (
315 surfI,
316 new surfaceZonesInfo
317 (
318 surface,
319 dict,
320 allGeometry_.regionNames()[surfaces_[surfI]]
321 )
322 );
323
324 // Global perpendicular angle
325 if (dict.found("patchInfo"))
326 {
327 globalPatchInfo.set
328 (
329 surfI,
330 dict.subDict("patchInfo").clone()
331 );
332 }
333 dict.readIfPresent("perpendicularAngle", globalAngle[surfI]);
334 dict.readIfPresent("blockLevel", globalBlockLevel[surfI]);
335 dict.readIfPresent("leakLevel", globalLeakLevel[surfI]);
336
337
338 if (dict.found("regions"))
339 {
340 const dictionary& regionsDict = dict.subDict("regions");
341 const wordList& regionNames = surface.regions();
342
343 forAll(regionNames, regionI)
344 {
345 if (regionsDict.found(regionNames[regionI]))
346 {
347 // Get the dictionary for region
348 const dictionary& regionDict = regionsDict.subDict
349 (
350 regionNames[regionI]
351 );
352
353 const labelPair refLevel
354 (
355 meshRefinement::get<labelPair>
356 (
357 regionDict,
358 "level",
359 dryRun_,
360 keyType::REGEX,
361 labelPair(0, 0)
362 )
363 );
364
365
366 regionMinLevel[surfI].insert(regionI, refLevel[0]);
367 regionMaxLevel[surfI].insert(regionI, refLevel[1]);
368 label levelIncr = regionDict.getOrDefault
369 (
370 "gapLevelIncrement",
371 gapLevelIncrement
372 );
373 regionLevelIncr[surfI].insert(regionI, levelIncr);
374
375 if
376 (
377 refLevel[0] < 0
378 || refLevel[1] < refLevel[0]
379 || levelIncr < 0
380 )
381 {
382 FatalIOErrorInFunction(dict)
383 << "Illegal level specification for surface "
384 << names_[surfI] << " region "
385 << regionNames[regionI]
386 << " : minLevel:" << refLevel[0]
387 << " maxLevel:" << refLevel[1]
388 << " levelIncrement:" << levelIncr
389 << exit(FatalIOError);
390 }
391
392
393
394 // Optional gapLevel specification
395
396 FixedList<label, 3> gapSpec
397 (
398 regionDict.getOrDefault
399 (
400 "gapLevel",
401 nullGapLevel
402 )
403 );
404 regionGapLevel[surfI].insert(regionI, gapSpec);
405 volumeType gapModeSpec
406 (
407 "gapMode",
408 regionDict,
409 volumeType::MIXED
410 );
411 regionGapMode[surfI].insert(regionI, gapModeSpec);
412 if
413 (
414 gapModeSpec == volumeType::UNKNOWN
415 || gapSpec[0] < 0
416 || gapSpec[1] < 0
417 || gapSpec[2] < 0
418 || gapSpec[1] > gapSpec[2]
419 )
420 {
421 FatalIOErrorInFunction(dict)
422 << "Illegal gapLevel specification for surface "
423 << names_[surfI]
424 << " : gapLevel:" << gapSpec
425 << " gapMode:" << gapModeSpec.str()
426 << exit(FatalIOError);
427 }
428 regionGapSelf[surfI].insert
429 (
430 regionI,
431 regionDict.getOrDefault<bool>
432 (
433 "gapSelf",
434 true
435 )
436 );
437
438 if (regionDict.found("perpendicularAngle"))
439 {
440 regionAngle[surfI].insert
441 (
442 regionI,
443 regionDict.get<scalar>("perpendicularAngle")
444 );
445 }
446
447 if (regionDict.found("patchInfo"))
448 {
449 regionPatchInfo[surfI].insert
450 (
451 regionI,
452 regionDict.subDict("patchInfo").clone()
453 );
454 }
455
456 label l;
457 if (regionDict.readIfPresent<label>("blockLevel", l))
458 {
459 regionBlockLevel[surfI].insert(regionI, l);
460 }
461 if (regionDict.readIfPresent<label>("leakLevel", l))
462 {
463 regionLeakLevel[surfI].insert(regionI, l);
464 }
465 }
466 }
467 }
468 surfI++;
469 }
470 }
471
472 if (unmatchedKeys.size() > 0)
473 {
474 IOWarningInFunction(surfacesDict)
475 << "Not all entries in refinementSurfaces dictionary were used."
476 << " The following entries were not used : "
477 << unmatchedKeys.sortedToc()
478 << endl;
479 }
480
481
482 // Calculate local to global region offset
483 label nRegions = 0;
484
485 forAll(surfaces_, surfI)
486 {
487 regionOffset_[surfI] = nRegions;
488 nRegions += allGeometry_[surfaces_[surfI]].regions().size();
489 }
490
491 // Rework surface specific information into information per global region
492
493 regionToSurface_ = calcSurfaceIndex(allGeometry_, surfaces_);
494
495
496 minLevel_.setSize(nRegions);
497 minLevel_ = 0;
498 maxLevel_.setSize(nRegions);
499 maxLevel_ = 0;
500 gapLevel_.setSize(nRegions);
501 gapLevel_ = -1;
502 extendedGapLevel_.setSize(nRegions);
503 extendedGapLevel_ = nullGapLevel;
504 extendedGapMode_.setSize(nRegions);
505 extendedGapMode_ = volumeType::UNKNOWN;
506 selfProximity_.setSize(nRegions);
507 selfProximity_ = true;
508 perpendicularAngle_.setSize(nRegions);
509 perpendicularAngle_ = -GREAT;
510 patchInfo_.setSize(nRegions);
511 blockLevel_.setSize(nRegions);
512 blockLevel_ = labelMax;
513 leakLevel_.setSize(nRegions);
514 leakLevel_ = labelMax;
515
516
517 forAll(globalMinLevel, surfI)
518 {
519 label nRegions = allGeometry_[surfaces_[surfI]].regions().size();
520
521 // Initialise to global (i.e. per surface)
522 for (label i = 0; i < nRegions; i++)
523 {
524 const label globalRegionI = regionOffset_[surfI] + i;
525
526 minLevel_[globalRegionI] = globalMinLevel[surfI];
527 maxLevel_[globalRegionI] = globalMaxLevel[surfI];
528 gapLevel_[globalRegionI] =
529 maxLevel_[globalRegionI]
530 + globalLevelIncr[surfI];
531 extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
532 extendedGapMode_[globalRegionI] = globalGapMode[surfI];
533 selfProximity_[globalRegionI] = globalGapSelf[surfI];
534 perpendicularAngle_[globalRegionI] = globalAngle[surfI];
535 if (globalPatchInfo.set(surfI))
536 {
537 patchInfo_.set
538 (
539 globalRegionI,
540 globalPatchInfo[surfI].clone()
541 );
542 }
543 blockLevel_[globalRegionI] = globalBlockLevel[surfI];
544 leakLevel_[globalRegionI] = globalLeakLevel[surfI];
545 }
546
547 // Overwrite with region specific information
548 forAllConstIters(regionMinLevel[surfI], iter)
549 {
550 const label globalRegionI = regionOffset_[surfI] + iter.key();
551
552 minLevel_[globalRegionI] = iter.val();
553 maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
554 gapLevel_[globalRegionI] =
555 maxLevel_[globalRegionI]
556 + regionLevelIncr[surfI][iter.key()];
557 extendedGapLevel_[globalRegionI] =
558 regionGapLevel[surfI][iter.key()];
559 extendedGapMode_[globalRegionI] =
560 regionGapMode[surfI][iter.key()];
561 selfProximity_[globalRegionI] =
562 regionGapSelf[surfI][iter.key()];
563 }
564 forAllConstIters(regionAngle[surfI], iter)
565 {
566 const label globalRegionI = regionOffset_[surfI] + iter.key();
567
568 perpendicularAngle_[globalRegionI] = iter.val();
569 }
570
571 const Map<autoPtr<dictionary>>& localInfo = regionPatchInfo[surfI];
572 forAllConstIters(localInfo, iter)
573 {
574 const label globalRegionI = regionOffset_[surfI] + iter.key();
575 const dictionary& dict = *(iter.val());
576
577 patchInfo_.set(globalRegionI, dict.clone());
578 }
579
580 forAllConstIters(regionBlockLevel[surfI], iter)
581 {
582 const label globalRegionI = regionOffset_[surfI] + iter.key();
583
584 blockLevel_[globalRegionI] = iter.val();
585 leakLevel_[globalRegionI] = iter.val();
586 }
587 }
588}
589
590
591Foam::refinementSurfaces::refinementSurfaces
592(
593 const searchableSurfaces& allGeometry,
594 const labelList& surfaces,
595 const wordList& names,
596 const PtrList<surfaceZonesInfo>& surfZones,
597 const labelList& regionOffset,
598 const labelList& minLevel,
599 const labelList& maxLevel,
600 const labelList& gapLevel,
601 const scalarField& perpendicularAngle,
602 PtrList<dictionary>& patchInfo,
603 const bool dryRun
604)
605:
606 allGeometry_(allGeometry),
607 surfaces_(surfaces),
608 names_(names),
609 surfZones_(surfZones),
610 regionOffset_(regionOffset),
611 regionToSurface_(calcSurfaceIndex(allGeometry, surfaces)),
612 minLevel_(minLevel),
613 maxLevel_(maxLevel),
614 gapLevel_(gapLevel),
615 perpendicularAngle_(perpendicularAngle),
616 patchInfo_(patchInfo.size()),
617 dryRun_(dryRun)
618{
619 forAll(patchInfo_, pI)
620 {
621 if (patchInfo.set(pI))
622 {
623 patchInfo_.set(pI, patchInfo.set(pI, nullptr));
624 }
625 }
626}
627
628
629// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
630
631Foam::labelPair Foam::refinementSurfaces::whichSurface
632(
633 const label globalRegionI
634) const
635{
636 const label surfI = regionToSurface_[globalRegionI];
637 const label localI = globalRegionI-regionOffset_[surfI];
638 return labelPair(surfI, localI);
639}
640
641
642// // Count number of triangles per surface region
643// Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
644// {
645// const geometricSurfacePatchList& regions = s.patches();
646//
647// labelList nTris(regions.size(), Zero);
648//
649// forAll(s, triI)
650// {
651// nTris[s[triI].region()]++;
652// }
653// return nTris;
654// }
655
656
657// // Pre-calculate the refinement level for every element
658// void Foam::refinementSurfaces::wantedRefinementLevel
659// (
660// const shellSurfaces& shells,
661// const label surfI,
662// const List<pointIndexHit>& info, // Indices
663// const pointField& ctrs, // Representative coordinate
664// labelList& minLevelField
665// ) const
666// {
667// const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
668//
669// // Get per element the region
670// labelList region;
671// geom.getRegion(info, region);
672//
673// // Initialise fields to region wise minLevel
674// minLevelField.setSize(ctrs.size());
675// minLevelField = -1;
676//
677// forAll(minLevelField, i)
678// {
679// if (info[i].hit())
680// {
681// minLevelField[i] = minLevel(surfI, region[i]);
682// }
683// }
684//
685// // Find out if triangle inside shell with higher level
686// // What level does shell want to refine fc to?
687// labelList shellLevel;
688// shells.findHigherLevel(ctrs, minLevelField, shellLevel);
689//
690// forAll(minLevelField, i)
691// {
692// minLevelField[i] = max(minLevelField[i], shellLevel[i]);
693// }
694// }
695
696
697Foam::labelList Foam::refinementSurfaces::maxGapLevel() const
698{
699 labelList surfaceMax(surfaces_.size(), Zero);
700
701 forAll(surfaces_, surfI)
702 {
703 const wordList& regionNames = allGeometry_[surfaces_[surfI]].regions();
704
705 forAll(regionNames, regionI)
706 {
707 label globalI = globalRegion(surfI, regionI);
708 const FixedList<label, 3>& gapInfo = extendedGapLevel_[globalI];
709 surfaceMax[surfI] = max(surfaceMax[surfI], gapInfo[2]);
710 }
711 }
712 return surfaceMax;
713}
714
715
716// Precalculate the refinement level for every element of the searchable
717// surface.
718void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
719{
720 forAll(surfaces_, surfI)
721 {
722 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
723
724 // Cache the refinement level (max of surface level and shell level)
725 // on a per-element basis. Only makes sense if there are lots of
726 // elements. Possibly should have 'enough' elements to have fine
727 // enough resolution but for now just make sure we don't catch e.g.
728 // searchableBox (size=6)
729 if (geom.globalSize() > 10)
730 {
731 // Representative local coordinates and bounding sphere
732 pointField ctrs;
733 scalarField radiusSqr;
734 geom.boundingSpheres(ctrs, radiusSqr);
735
736 labelList minLevelField(ctrs.size(), Zero);
737 {
738 // Get the element index in a roundabout way. Problem is e.g.
739 // distributed surface where local indices differ from global
740 // ones (needed for getRegion call)
741 List<pointIndexHit> info;
742 geom.findNearest(ctrs, radiusSqr, info);
743
744 // Get per element the region
745 labelList region;
746 geom.getRegion(info, region);
747
748 // From the region get the surface-wise refinement level
749 forAll(minLevelField, i)
750 {
751 if (info[i].hit()) //Note: should not be necessary
752 {
753 minLevelField[i] = minLevel(surfI, region[i]);
754 }
755 }
756 }
757
758 // Find out if triangle inside shell with higher level
759 // What level does shell want to refine fc to?
760 labelList shellLevel;
761 shells.findHigherLevel(ctrs, minLevelField, shellLevel);
762
763
764 // In case of triangulated surfaces only cache value if triangle
765 // centre and vertices are in same shell
766 if (isA<triSurface>(geom))
767 {
768 label nUncached = 0;
769
770 // Check if points differing from ctr level
771
772 const triSurface& ts = refCast<const triSurface>(geom);
773 const pointField& points = ts.points();
774
775 // Determine minimum expected level to avoid having to
776 // test lots of points
777 labelList minPointLevel(points.size(), labelMax);
778 forAll(shellLevel, triI)
779 {
780 const labelledTri& t = ts[triI];
781 label level = shellLevel[triI];
782 forAll(t, tI)
783 {
784 minPointLevel[t[tI]] = min(minPointLevel[t[tI]], level);
785 }
786 }
787
788
789 // See if inside any shells with higher refinement level
790 labelList pointLevel;
791 shells.findHigherLevel(points, minPointLevel, pointLevel);
792
793
794 // See if triangle centre values differ from triangle points
795 forAll(shellLevel, triI)
796 {
797 const labelledTri& t = ts[triI];
798 label fLevel = shellLevel[triI];
799 if
800 (
801 (pointLevel[t[0]] != fLevel)
802 || (pointLevel[t[1]] != fLevel)
803 || (pointLevel[t[2]] != fLevel)
804 )
805 {
806 //Pout<< "Detected triangle " << t.tri(ts.points())
807 // << " partially inside/partially outside" << endl;
808
809 // Mark as uncached
810 shellLevel[triI] = -1;
811 nUncached++;
812 }
813 }
814
815 if (!dryRun_)
816 {
817 Info<< "For geometry " << geom.name()
818 << " detected "
819 << returnReduce(nUncached, sumOp<label>())
820 << " uncached triangles out of " << geom.globalSize()
821 << endl;
822 }
823 }
824
825
826 // Combine overall level field with current shell level. Make sure
827 // to preserve -1 (from triSurfaceMeshes with triangles partly
828 // inside/outside
829 forAll(minLevelField, i)
830 {
831 if (min(minLevelField[i], shellLevel[i]) < 0)
832 {
833 minLevelField[i] = -1;
834 }
835 else
836 {
837 minLevelField[i] = max(minLevelField[i], shellLevel[i]);
838 }
839 }
840
841 // Store minLevelField on surface
842 const_cast<searchableSurface&>(geom).setField(minLevelField);
843 }
844 }
845}
846
847
848// Find intersections of edge. Return -1 or first surface with higher minLevel
849// number.
850void Foam::refinementSurfaces::findHigherIntersection
851(
852 const shellSurfaces& shells,
853
854 const pointField& start,
855 const pointField& end,
856 const labelList& currentLevel, // current cell refinement level
857
858 labelList& surfaces,
859 labelList& surfaceLevel
860) const
861{
862 surfaces.setSize(start.size());
863 surfaces = -1;
864 surfaceLevel.setSize(start.size());
865 surfaceLevel = -1;
866
867 if (surfaces_.empty())
868 {
869 return;
870 }
871
872 if (surfaces_.size() == 1)
873 {
874 // Optimisation: single segmented surface. No need to duplicate
875 // point storage.
876
877 label surfI = 0;
878
879 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
880
881 // Do intersection test
882 List<pointIndexHit> intersectionInfo(start.size());
883 geom.findLineAny(start, end, intersectionInfo);
884
885
886 // Surface-based refinement level
887 labelList surfaceOnlyLevel(start.size(), -1);
888 {
889 // Get per intersection the region
890 labelList region;
891 geom.getRegion(intersectionInfo, region);
892
893 forAll(intersectionInfo, i)
894 {
895 if (intersectionInfo[i].hit())
896 {
897 surfaceOnlyLevel[i] = minLevel(surfI, region[i]);
898 }
899 }
900 }
901
902
903 // Get shell refinement level if higher
904 const labelList localLevel
905 (
906 findHigherLevel
907 (
908 geom,
909 shells,
910 intersectionInfo,
911 surfaceOnlyLevel // starting level
912 )
913 );
914
915
916 // Combine localLevel with current level
917 forAll(localLevel, i)
918 {
919 if (localLevel[i] > currentLevel[i])
920 {
921 surfaces[i] = surfI; // index of surface
922 surfaceLevel[i] = localLevel[i];
923 }
924 }
925
926 return;
927 }
928
929
930
931 // Work arrays
932 pointField p0(start);
933 pointField p1(end);
934 labelList intersectionToPoint(identity(start.size()));
935 List<pointIndexHit> intersectionInfo(start.size());
936
937 forAll(surfaces_, surfI)
938 {
939 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
940
941 // Do intersection test
942 geom.findLineAny(p0, p1, intersectionInfo);
943
944
945 // Surface-based refinement level
946 labelList surfaceOnlyLevel(intersectionInfo.size(), -1);
947 {
948 // Get per intersection the region
949 labelList region;
950 geom.getRegion(intersectionInfo, region);
951
952 forAll(intersectionInfo, i)
953 {
954 if (intersectionInfo[i].hit())
955 {
956 surfaceOnlyLevel[i] = minLevel(surfI, region[i]);
957 }
958 }
959 }
960
961
962 // Get shell refinement level if higher
963 const labelList localLevel
964 (
965 findHigherLevel
966 (
967 geom,
968 shells,
969 intersectionInfo,
970 surfaceOnlyLevel
971 )
972 );
973
974
975 // Combine localLevel with current level
976 label missI = 0;
977 forAll(localLevel, i)
978 {
979 label pointI = intersectionToPoint[i];
980
981 if (localLevel[i] > currentLevel[pointI])
982 {
983 // Mark point for refinement
984 surfaces[pointI] = surfI;
985 surfaceLevel[pointI] = localLevel[i];
986 }
987 else
988 {
989 p0[missI] = start[pointI];
990 p1[missI] = end[pointI];
991 intersectionToPoint[missI] = pointI;
992 missI++;
993 }
994 }
995
996
997 // All done? Note that this decision should be synchronised
998 if (returnReduce(missI, sumOp<label>()) == 0)
999 {
1000 break;
1001 }
1002
1003 // Trim misses
1004 p0.setSize(missI);
1005 p1.setSize(missI);
1006 intersectionToPoint.setSize(missI);
1007 intersectionInfo.setSize(missI);
1008 }
1009}
1010
1011
1012void Foam::refinementSurfaces::findAllIntersections
1013(
1014 const pointField& start,
1015 const pointField& end,
1016 const labelList& currentLevel, // current cell refinement level
1017
1018 const labelList& globalMinLevel,
1019 const labelList& globalMaxLevel,
1020
1021 List<vectorList>& surfaceNormal,
1022 labelListList& surfaceLevel
1023) const
1024{
1025 surfaceLevel.setSize(start.size());
1026 surfaceNormal.setSize(start.size());
1027
1028 if (surfaces_.empty())
1029 {
1030 return;
1031 }
1032
1033 // Work arrays
1034 List<List<pointIndexHit>> hitInfo;
1035 labelList pRegions;
1036 vectorField pNormals;
1037
1038 forAll(surfaces_, surfI)
1039 {
1040 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1041
1042 surface.findLineAll(start, end, hitInfo);
1043
1044 // Repack hits for surface into flat list
1045 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046 // To avoid overhead of calling getRegion for every point
1047
1048 label n = 0;
1049 forAll(hitInfo, pointI)
1050 {
1051 n += hitInfo[pointI].size();
1052 }
1053
1054 List<pointIndexHit> surfInfo(n);
1055 labelList pointMap(n);
1056 n = 0;
1057
1058 forAll(hitInfo, pointI)
1059 {
1060 const List<pointIndexHit>& pHits = hitInfo[pointI];
1061
1062 forAll(pHits, i)
1063 {
1064 surfInfo[n] = pHits[i];
1065 pointMap[n] = pointI;
1066 n++;
1067 }
1068 }
1069
1070 labelList surfRegion(n);
1071 vectorField surfNormal(n);
1072 surface.getRegion(surfInfo, surfRegion);
1073 surface.getNormal(surfInfo, surfNormal);
1074
1075 surfInfo.clear();
1076
1077
1078 // Extract back into pointwise
1079 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1080
1081 forAll(surfRegion, i)
1082 {
1083 label region = globalRegion(surfI, surfRegion[i]);
1084 label pointI = pointMap[i];
1085
1086 if
1087 (
1088 currentLevel[pointI] >= globalMinLevel[region]
1089 && currentLevel[pointI] < globalMaxLevel[region]
1090 )
1091 {
1092 // Append to pointI info
1093 label sz = surfaceNormal[pointI].size();
1094 surfaceNormal[pointI].setSize(sz+1);
1095 surfaceNormal[pointI][sz] = surfNormal[i];
1096
1097 surfaceLevel[pointI].setSize(sz+1);
1098 surfaceLevel[pointI][sz] = globalMaxLevel[region];
1099 }
1100 }
1101 }
1102}
1103
1104
1105void Foam::refinementSurfaces::findAllIntersections
1106(
1107 const pointField& start,
1108 const pointField& end,
1109 const labelList& currentLevel, // current cell refinement level
1110
1111 const labelList& globalMinLevel,
1112 const labelList& globalMaxLevel,
1113
1114 List<pointList>& surfaceLocation,
1115 List<vectorList>& surfaceNormal,
1116 labelListList& surfaceLevel
1117) const
1118{
1119 surfaceLevel.setSize(start.size());
1120 surfaceNormal.setSize(start.size());
1121 surfaceLocation.setSize(start.size());
1122
1123 if (surfaces_.empty())
1124 {
1125 return;
1126 }
1127
1128 // Work arrays
1129 List<List<pointIndexHit>> hitInfo;
1130 labelList pRegions;
1131 vectorField pNormals;
1132
1133 forAll(surfaces_, surfI)
1134 {
1135 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1136
1137 surface.findLineAll(start, end, hitInfo);
1138
1139 // Repack hits for surface into flat list
1140 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1141 // To avoid overhead of calling getRegion for every point
1142
1143 label n = 0;
1144 forAll(hitInfo, pointI)
1145 {
1146 n += hitInfo[pointI].size();
1147 }
1148
1149 List<pointIndexHit> surfInfo(n);
1150 labelList pointMap(n);
1151 n = 0;
1152
1153 forAll(hitInfo, pointI)
1154 {
1155 const List<pointIndexHit>& pHits = hitInfo[pointI];
1156
1157 forAll(pHits, i)
1158 {
1159 surfInfo[n] = pHits[i];
1160 pointMap[n] = pointI;
1161 n++;
1162 }
1163 }
1164
1165 labelList surfRegion(n);
1166 vectorField surfNormal(n);
1167 surface.getRegion(surfInfo, surfRegion);
1168 surface.getNormal(surfInfo, surfNormal);
1169
1170 // Extract back into pointwise
1171 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1172
1173 forAll(surfRegion, i)
1174 {
1175 label region = globalRegion(surfI, surfRegion[i]);
1176 label pointI = pointMap[i];
1177
1178 if
1179 (
1180 currentLevel[pointI] >= globalMinLevel[region]
1181 && currentLevel[pointI] < globalMaxLevel[region]
1182 )
1183 {
1184 // Append to pointI info
1185 label sz = surfaceNormal[pointI].size();
1186 surfaceLocation[pointI].setSize(sz+1);
1187 surfaceLocation[pointI][sz] = surfInfo[i].hitPoint();
1188
1189 surfaceNormal[pointI].setSize(sz+1);
1190 surfaceNormal[pointI][sz] = surfNormal[i];
1191
1192 surfaceLevel[pointI].setSize(sz+1);
1193 // Level should just be higher than provided point level.
1194 // Actual value is not important.
1195 surfaceLevel[pointI][sz] = globalMaxLevel[region];
1196 }
1197 }
1198 }
1199}
1200
1201
1202void Foam::refinementSurfaces::findNearestIntersection
1203(
1204 const labelList& surfacesToTest,
1205 const pointField& start,
1206 const pointField& end,
1207
1208 labelList& surface1,
1209 List<pointIndexHit>& hit1,
1210 labelList& region1,
1211 labelList& surface2,
1212 List<pointIndexHit>& hit2,
1213 labelList& region2
1214) const
1215{
1216 // 1. intersection from start to end
1217 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1218
1219 // Initialize arguments
1220 surface1.setSize(start.size());
1221 surface1 = -1;
1222 hit1.setSize(start.size());
1223 region1.setSize(start.size());
1224
1225 // Current end of segment to test.
1226 pointField nearest(end);
1227 // Work array
1228 List<pointIndexHit> nearestInfo(start.size());
1229 labelList region;
1230
1231 forAll(surfacesToTest, testI)
1232 {
1233 label surfI = surfacesToTest[testI];
1234
1235 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1236
1237 // See if any intersection between start and current nearest
1238 surface.findLine
1239 (
1240 start,
1241 nearest,
1242 nearestInfo
1243 );
1244 surface.getRegion
1245 (
1246 nearestInfo,
1247 region
1248 );
1249
1250 forAll(nearestInfo, pointI)
1251 {
1252 if (nearestInfo[pointI].hit())
1253 {
1254 hit1[pointI] = nearestInfo[pointI];
1255 surface1[pointI] = surfI;
1256 region1[pointI] = region[pointI];
1257 nearest[pointI] = hit1[pointI].hitPoint();
1258 }
1259 }
1260 }
1261
1262
1263 // 2. intersection from end to last intersection
1264 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1265
1266 // Find the nearest intersection from end to start. Note that we initialize
1267 // to the first intersection (if any).
1268 surface2 = surface1;
1269 hit2 = hit1;
1270 region2 = region1;
1271
1272 // Set current end of segment to test.
1273 forAll(nearest, pointI)
1274 {
1275 if (hit1[pointI].hit())
1276 {
1277 nearest[pointI] = hit1[pointI].hitPoint();
1278 }
1279 else
1280 {
1281 // Disable testing by setting to end.
1282 nearest[pointI] = end[pointI];
1283 }
1284 }
1285
1286 forAll(surfacesToTest, testI)
1287 {
1288 label surfI = surfacesToTest[testI];
1289
1290 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1291
1292 // See if any intersection between end and current nearest
1293 surface.findLine
1294 (
1295 end,
1296 nearest,
1297 nearestInfo
1298 );
1299 surface.getRegion
1300 (
1301 nearestInfo,
1302 region
1303 );
1304
1305 forAll(nearestInfo, pointI)
1306 {
1307 if (nearestInfo[pointI].hit())
1308 {
1309 hit2[pointI] = nearestInfo[pointI];
1310 surface2[pointI] = surfI;
1311 region2[pointI] = region[pointI];
1312 nearest[pointI] = hit2[pointI].hitPoint();
1313 }
1314 }
1315 }
1316
1317
1318 // Make sure that if hit1 has hit something, hit2 will have at least the
1319 // same point (due to tolerances it might miss its end point)
1320 forAll(hit1, pointI)
1321 {
1322 if (hit1[pointI].hit() && !hit2[pointI].hit())
1323 {
1324 hit2[pointI] = hit1[pointI];
1325 surface2[pointI] = surface1[pointI];
1326 region2[pointI] = region1[pointI];
1327 }
1328 }
1329}
1330
1331
1332void Foam::refinementSurfaces::findNearestIntersection
1333(
1334 const labelList& surfacesToTest,
1335 const pointField& start,
1336 const pointField& end,
1337
1338 labelList& surface1,
1339 List<pointIndexHit>& hit1,
1340 labelList& region1,
1341 vectorField& normal1,
1342
1343 labelList& surface2,
1344 List<pointIndexHit>& hit2,
1345 labelList& region2,
1346 vectorField& normal2
1347) const
1348{
1349 // 1. intersection from start to end
1350 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1351
1352 // Initialize arguments
1353 surface1.setSize(start.size());
1354 surface1 = -1;
1355 hit1.setSize(start.size());
1356 region1.setSize(start.size());
1357 region1 = -1;
1358 normal1.setSize(start.size());
1359 normal1 = Zero;
1360
1361 // Current end of segment to test.
1362 pointField nearest(end);
1363 // Work array
1364 List<pointIndexHit> nearestInfo(start.size());
1365 labelList region;
1366 vectorField normal;
1367
1368 forAll(surfacesToTest, testI)
1369 {
1370 label surfI = surfacesToTest[testI];
1371 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
1372
1373 // See if any intersection between start and current nearest
1374 geom.findLine(start, nearest, nearestInfo);
1375 geom.getRegion(nearestInfo, region);
1376 geom.getNormal(nearestInfo, normal);
1377
1378 forAll(nearestInfo, pointI)
1379 {
1380 if (nearestInfo[pointI].hit())
1381 {
1382 hit1[pointI] = nearestInfo[pointI];
1383 surface1[pointI] = surfI;
1384 region1[pointI] = region[pointI];
1385 normal1[pointI] = normal[pointI];
1386 nearest[pointI] = hit1[pointI].hitPoint();
1387 }
1388 }
1389 }
1390
1391
1392 // 2. intersection from end to last intersection
1393 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1394
1395 // Find the nearest intersection from end to start. Note that we initialize
1396 // to the first intersection (if any).
1397 surface2 = surface1;
1398 hit2 = hit1;
1399 region2 = region1;
1400 normal2 = normal1;
1401
1402 // Set current end of segment to test.
1403 forAll(nearest, pointI)
1404 {
1405 if (hit1[pointI].hit())
1406 {
1407 nearest[pointI] = hit1[pointI].hitPoint();
1408 }
1409 else
1410 {
1411 // Disable testing by setting to end.
1412 nearest[pointI] = end[pointI];
1413 }
1414 }
1415
1416 forAll(surfacesToTest, testI)
1417 {
1418 label surfI = surfacesToTest[testI];
1419 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
1420
1421 // See if any intersection between end and current nearest
1422 geom.findLine(end, nearest, nearestInfo);
1423 geom.getRegion(nearestInfo, region);
1424 geom.getNormal(nearestInfo, normal);
1425
1426 forAll(nearestInfo, pointI)
1427 {
1428 if (nearestInfo[pointI].hit())
1429 {
1430 hit2[pointI] = nearestInfo[pointI];
1431 surface2[pointI] = surfI;
1432 region2[pointI] = region[pointI];
1433 normal2[pointI] = normal[pointI];
1434 nearest[pointI] = hit2[pointI].hitPoint();
1435 }
1436 }
1437 }
1438
1439
1440 // Make sure that if hit1 has hit something, hit2 will have at least the
1441 // same point (due to tolerances it might miss its end point)
1442 forAll(hit1, pointI)
1443 {
1444 if (hit1[pointI].hit() && !hit2[pointI].hit())
1445 {
1446 hit2[pointI] = hit1[pointI];
1447 surface2[pointI] = surface1[pointI];
1448 region2[pointI] = region1[pointI];
1449 normal2[pointI] = normal1[pointI];
1450 }
1451 }
1452}
1453
1454
1455void Foam::refinementSurfaces::findNearestIntersection
1456(
1457 const pointField& start,
1458 const pointField& end,
1459
1460 labelList& surface1,
1461 vectorField& normal1
1462) const
1463{
1464 // Initialize arguments
1465 surface1.setSize(start.size());
1466 surface1 = -1;
1467 normal1.setSize(start.size());
1468 normal1 = Zero;
1469
1470 // Current end of segment to test.
1471 pointField nearest(end);
1472 // Work array
1473 List<pointIndexHit> nearestInfo(start.size());
1474 labelList region;
1475 vectorField normal;
1476
1477 forAll(surfaces_, surfI)
1478 {
1479 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
1480
1481 // See if any intersection between start and current nearest
1482 geom.findLine(start, nearest, nearestInfo);
1483 geom.getNormal(nearestInfo, normal);
1484
1485 forAll(nearestInfo, pointI)
1486 {
1487 if (nearestInfo[pointI].hit())
1488 {
1489 surface1[pointI] = surfI;
1490 normal1[pointI] = normal[pointI];
1491 nearest[pointI] = nearestInfo[pointI].hitPoint();
1492 }
1493 }
1494 }
1495}
1496
1497
1498void Foam::refinementSurfaces::findNearestIntersection
1499(
1500 const pointField& start,
1501 const pointField& end,
1502
1503 labelList& surface1,
1504 List<pointIndexHit>& hitInfo1,
1505 vectorField& normal1
1506) const
1507{
1508 // 1. intersection from start to end
1509 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1510
1511 // Initialize arguments
1512 surface1.setSize(start.size());
1513 surface1 = -1;
1514 hitInfo1.setSize(start.size());
1515 hitInfo1 = pointIndexHit();
1516 normal1.setSize(start.size());
1517 normal1 = Zero;
1518
1519 // Current end of segment to test.
1520 pointField nearest(end);
1521 // Work array
1522 List<pointIndexHit> nearestInfo(start.size());
1523 labelList region;
1524 vectorField normal;
1525
1526 forAll(surfaces_, surfI)
1527 {
1528 const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
1529
1530 // See if any intersection between start and current nearest
1531 geom.findLine(start, nearest, nearestInfo);
1532 geom.getNormal(nearestInfo, normal);
1533
1534 forAll(nearestInfo, pointI)
1535 {
1536 if (nearestInfo[pointI].hit())
1537 {
1538 surface1[pointI] = surfI;
1539 hitInfo1[pointI] = nearestInfo[pointI];
1540 normal1[pointI] = normal[pointI];
1541 nearest[pointI] = nearestInfo[pointI].hitPoint();
1542 }
1543 }
1544 }
1545}
1546
1547
1548void Foam::refinementSurfaces::findAnyIntersection
1549(
1550 const pointField& start,
1551 const pointField& end,
1552
1553 labelList& hitSurface,
1554 List<pointIndexHit>& hitInfo
1555) const
1556{
1557 searchableSurfacesQueries::findAnyIntersection
1558 (
1559 allGeometry_,
1560 surfaces_,
1561 start,
1562 end,
1563 hitSurface,
1564 hitInfo
1565 );
1566}
1567
1568
1569void Foam::refinementSurfaces::findNearest
1570(
1571 const labelList& surfacesToTest,
1572 const pointField& samples,
1573 const scalarField& nearestDistSqr,
1574 labelList& hitSurface,
1575 List<pointIndexHit>& hitInfo
1576) const
1577{
1578 labelList geometries(labelUIndList(surfaces_, surfacesToTest));
1579
1580 // Do the tests. Note that findNearest returns index in geometries.
1581 searchableSurfacesQueries::findNearest
1582 (
1583 allGeometry_,
1584 geometries,
1585 samples,
1586 nearestDistSqr,
1587 hitSurface,
1588 hitInfo
1589 );
1590
1591 // Rework the hitSurface to be surface (i.e. index into surfaces_)
1592 forAll(hitSurface, pointI)
1593 {
1594 if (hitSurface[pointI] != -1)
1595 {
1596 hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1597 }
1598 }
1599}
1600
1601
1602void Foam::refinementSurfaces::findNearestRegion
1603(
1604 const labelList& surfacesToTest,
1605 const pointField& samples,
1606 const scalarField& nearestDistSqr,
1607 labelList& hitSurface,
1608 labelList& hitRegion
1609) const
1610{
1611 labelList geometries(labelUIndList(surfaces_, surfacesToTest));
1612
1613 // Do the tests. Note that findNearest returns index in geometries.
1614 List<pointIndexHit> hitInfo;
1615 searchableSurfacesQueries::findNearest
1616 (
1617 allGeometry_,
1618 geometries,
1619 samples,
1620 nearestDistSqr,
1621 hitSurface,
1622 hitInfo
1623 );
1624
1625 // Rework the hitSurface to be surface (i.e. index into surfaces_)
1626 forAll(hitSurface, pointI)
1627 {
1628 if (hitSurface[pointI] != -1)
1629 {
1630 hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1631 }
1632 }
1633
1634 // Collect the region
1635 hitRegion.setSize(hitSurface.size());
1636 hitRegion = -1;
1637
1638 forAll(surfacesToTest, i)
1639 {
1640 label surfI = surfacesToTest[i];
1641
1642 // Collect hits for surfI
1643 const labelList localIndices(findIndices(hitSurface, surfI));
1644
1645 List<pointIndexHit> localHits
1646 (
1647 UIndirectList<pointIndexHit>
1648 (
1649 hitInfo,
1650 localIndices
1651 )
1652 );
1653
1654 labelList localRegion;
1655 allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1656
1657 forAll(localIndices, i)
1658 {
1659 hitRegion[localIndices[i]] = localRegion[i];
1660 }
1661 }
1662}
1663
1664
1665void Foam::refinementSurfaces::findNearestRegion
1666(
1667 const labelList& surfacesToTest,
1668 const pointField& samples,
1669 const scalarField& nearestDistSqr,
1670 labelList& hitSurface,
1671 List<pointIndexHit>& hitInfo,
1672 labelList& hitRegion,
1673 vectorField& hitNormal
1674) const
1675{
1676 labelList geometries(labelUIndList(surfaces_, surfacesToTest));
1677
1678 // Do the tests. Note that findNearest returns index in geometries.
1679 searchableSurfacesQueries::findNearest
1680 (
1681 allGeometry_,
1682 geometries,
1683 samples,
1684 nearestDistSqr,
1685 hitSurface,
1686 hitInfo
1687 );
1688
1689 // Rework the hitSurface to be surface (i.e. index into surfaces_)
1690 forAll(hitSurface, pointI)
1691 {
1692 if (hitSurface[pointI] != -1)
1693 {
1694 hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1695 }
1696 }
1697
1698 // Collect the region
1699 hitRegion.setSize(hitSurface.size());
1700 hitRegion = -1;
1701 hitNormal.setSize(hitSurface.size());
1702 hitNormal = Zero;
1703
1704 forAll(surfacesToTest, i)
1705 {
1706 label surfI = surfacesToTest[i];
1707
1708 // Collect hits for surfI
1709 const labelList localIndices(findIndices(hitSurface, surfI));
1710
1711 List<pointIndexHit> localHits
1712 (
1713 UIndirectList<pointIndexHit>
1714 (
1715 hitInfo,
1716 localIndices
1717 )
1718 );
1719
1720 // Region
1721 labelList localRegion;
1722 allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1723
1724 forAll(localIndices, i)
1725 {
1726 hitRegion[localIndices[i]] = localRegion[i];
1727 }
1728
1729 // Normal
1730 vectorField localNormal;
1731 allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
1732
1733 forAll(localIndices, i)
1734 {
1735 hitNormal[localIndices[i]] = localNormal[i];
1736 }
1737 }
1738}
1739
1740
1743//Foam::label Foam::refinementSurfaces::findHighestIntersection
1744//(
1745// const point& start,
1746// const point& end,
1747// const label currentLevel, // current cell refinement level
1748//
1749// pointIndexHit& maxHit
1750//) const
1751//{
1752// // surface with highest maxlevel
1753// label maxSurface = -1;
1754// // maxLevel of maxSurface
1755// label maxLevel = currentLevel;
1756//
1757// forAll(*this, surfI)
1758// {
1759// pointIndexHit hit = operator[](surfI).findLineAny(start, end);
1760//
1761// if (hit.hit())
1762// {
1763// const triSurface& s = operator[](surfI);
1764//
1765// label region = globalRegion(surfI, s[hit.index()].region());
1766//
1767// if (maxLevel_[region] > maxLevel)
1768// {
1769// maxSurface = surfI;
1770// maxLevel = maxLevel_[region];
1771// maxHit = hit;
1772// }
1773// }
1774// }
1775//
1776// if (maxSurface == -1)
1777// {
1778// // maxLevel unchanged. No interesting surface hit.
1779// maxHit.setMiss();
1780// }
1781//
1782// return maxSurface;
1783//}
1784
1785
1786void Foam::refinementSurfaces::findInside
1787(
1788 const labelList& testSurfaces,
1789 const pointField& pt,
1790 labelList& insideSurfaces
1791) const
1792{
1793 insideSurfaces.setSize(pt.size());
1794 insideSurfaces = -1;
1795
1796 forAll(testSurfaces, i)
1797 {
1798 label surfI = testSurfaces[i];
1799
1800 const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1801
1802 const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
1803 surfZones_[surfI].zoneInside();
1804
1805 if
1806 (
1807 selectionMethod != surfaceZonesInfo::INSIDE
1808 && selectionMethod != surfaceZonesInfo::OUTSIDE
1809 )
1810 {
1811 FatalErrorInFunction
1812 << "Trying to use surface "
1813 << surface.name()
1814 << " which has non-geometric inside selection method "
1815 << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
1816 << exit(FatalError);
1817 }
1818
1819 if (surface.hasVolumeType())
1820 {
1821 List<volumeType> volType;
1822 surface.getVolumeType(pt, volType);
1823
1824 forAll(volType, pointI)
1825 {
1826 if (insideSurfaces[pointI] == -1)
1827 {
1828 if
1829 (
1830 (
1831 volType[pointI] == volumeType::INSIDE
1832 && selectionMethod == surfaceZonesInfo::INSIDE
1833 )
1834 || (
1835 volType[pointI] == volumeType::OUTSIDE
1836 && selectionMethod == surfaceZonesInfo::OUTSIDE
1837 )
1838 )
1839 {
1840 insideSurfaces[pointI] = surfI;
1841 }
1842 }
1843 }
1844 }
1845 }
1846}
1847
1848
1849void Foam::refinementSurfaces::findNearest
1850(
1851 const labelList& surfacesToTest,
1852 const labelListList& regions,
1853
1854 const pointField& samples,
1855 const scalarField& nearestDistSqr,
1856
1857 labelList& hitSurface,
1858 List<pointIndexHit>& hitInfo
1859) const
1860{
1861 labelList geometries(labelUIndList(surfaces_, surfacesToTest));
1862
1863 // Do the tests. Note that findNearest returns index in geometries.
1864 searchableSurfacesQueries::findNearest
1865 (
1866 allGeometry_,
1867 geometries,
1868 regions,
1869 samples,
1870 nearestDistSqr,
1871 hitSurface,
1872 hitInfo
1873 );
1874
1875 // Rework the hitSurface to be surface (i.e. index into surfaces_)
1876 forAll(hitSurface, pointI)
1877 {
1878 if (hitSurface[pointI] != -1)
1879 {
1880 hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1881 }
1882 }
1883}
1884
1885
1886void Foam::refinementSurfaces::findNearestRegion
1887(
1888 const labelList& surfacesToTest,
1889 const labelListList& regions,
1890
1891 const pointField& samples,
1892 const scalarField& nearestDistSqr,
1893
1894 labelList& hitSurface,
1895 List<pointIndexHit>& hitInfo,
1896 labelList& hitRegion,
1897 vectorField& hitNormal
1898) const
1899{
1900 labelList geometries(labelUIndList(surfaces_, surfacesToTest));
1901
1902 // Do the tests. Note that findNearest returns index in geometries.
1903 searchableSurfacesQueries::findNearest
1904 (
1905 allGeometry_,
1906 geometries,
1907 regions,
1908 samples,
1909 nearestDistSqr,
1910 hitSurface,
1911 hitInfo
1912 );
1913
1914 // Rework the hitSurface to be surface (i.e. index into surfaces_)
1915 forAll(hitSurface, pointI)
1916 {
1917 if (hitSurface[pointI] != -1)
1918 {
1919 hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1920 }
1921 }
1922
1923 // Collect the region
1924 hitRegion.setSize(hitSurface.size());
1925 hitRegion = -1;
1926 hitNormal.setSize(hitSurface.size());
1927 hitNormal = Zero;
1928
1929 forAll(surfacesToTest, i)
1930 {
1931 label surfI = surfacesToTest[i];
1932
1933 // Collect hits for surfI
1934 const labelList localIndices(findIndices(hitSurface, surfI));
1935
1936 List<pointIndexHit> localHits
1937 (
1938 UIndirectList<pointIndexHit>
1939 (
1940 hitInfo,
1941 localIndices
1942 )
1943 );
1944
1945 // Region
1946 labelList localRegion;
1947 allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1948
1949 forAll(localIndices, i)
1950 {
1951 hitRegion[localIndices[i]] = localRegion[i];
1952 }
1953
1954 // Normal
1955 vectorField localNormal;
1956 allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
1957
1958 forAll(localIndices, i)
1959 {
1960 hitNormal[localIndices[i]] = localNormal[i];
1961 }
1962 }
1963}
1964
1965
1966// ************************************************************************* //
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::FixedList
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
Foam::HashTable::sortedToc
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition: HashTable.C:137
Foam::HashTable::size
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Foam::HashTable::erase
bool erase(const iterator &iter)
Erase an entry specified by given iterator.
Definition: HashTable.C:440
Foam::IOobject::name
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
Foam::Map
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
Foam::PrimitivePatch::points
const Field< point_type > & points() const noexcept
Return reference to global points.
Definition: PrimitivePatch.H:299
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
Foam::PtrList::set
const T * set(const label i) const
Definition: PtrList.H:138
Foam::UIndirectList
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Foam::UPtrList::size
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Foam::dictionary::get
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:107
Foam::dictionary::subDict
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
Foam::dictionary::clone
autoPtr< dictionary > clone() const
Construct and return clone.
Definition: dictionary.C:172
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::dictionary::findEntry
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
Foam::dictionary::found
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
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:405
Foam::dictionary::toc
wordList toc() const
Return the table of contents.
Definition: dictionary.C:602
Foam::entry
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:70
Foam::entry::keyword
const keyType & keyword() const noexcept
Return keyword.
Definition: entry.H:195
Foam::entry::dict
virtual const dictionary & dict() const =0
Return dictionary, if entry is a dictionary.
Foam::keyType::REGEX
@ REGEX
Regular expression.
Definition: keyType.H:82
Foam::labelledTri
A triFace with additional (region) index.
Definition: labelledTri.H:60
Foam::refinementSurfaces
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
Definition: refinementSurfaces.H:64
Foam::refinementSurfaces::maxGapLevel
labelList maxGapLevel() const
Per surface the maximum extendedGapLevel over all its regions.
Definition: refinementSurfaces.C:697
Foam::refinementSurfaces::findAllIntersections
void findAllIntersections(const pointField &start, const pointField &end, const labelList &currentLevel, const labelList &globalMinLevel, const labelList &globalMaxLevel, List< vectorList > &surfaceNormal, labelListList &surfaceLevel) const
Find all intersections of edge with any surface with applicable.
Definition: refinementSurfaces.C:1013
Foam::refinementSurfaces::patchInfo
const PtrList< dictionary > & patchInfo() const
From global region number to patch type.
Definition: refinementSurfaces.H:272
Foam::refinementSurfaces::whichSurface
labelPair whichSurface(const label globalRegionI) const
From global region to surface + region.
Definition: refinementSurfaces.C:632
Foam::refinementSurfaces::findNearest
void findNearest(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest point on surfaces.
Definition: refinementSurfaces.C:1570
Foam::refinementSurfaces::findInside
void findInside(const labelList &surfacesToTest, const pointField &pt, labelList &insideSurfaces) const
Detect if a point is 'inside' (closed) surfaces.
Definition: refinementSurfaces.C:1787
Foam::refinementSurfaces::findNearestRegion
void findNearestRegion(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &hitSurface, labelList &hitRegion) const
Find nearest point on surfaces. Return surface and region on.
Definition: refinementSurfaces.C:1603
Foam::refinementSurfaces::findAnyIntersection
void findAnyIntersection(const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &) const
Used for debugging only: find intersection of edge.
Definition: refinementSurfaces.C:1549
Foam::refinementSurfaces::findHigherIntersection
void findHigherIntersection(const shellSurfaces &shells, const pointField &start, const pointField &end, const labelList &currentLevel, labelList &surfaces, labelList &surfaceLevel) const
Find intersection of edge. Return -1 or first surface.
Definition: refinementSurfaces.C:851
Foam::refinementSurfaces::findNearestIntersection
void findNearestIntersection(const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &region1, labelList &surface2, List< pointIndexHit > &hit2, labelList &region2) const
Find intersection nearest to the endpoints. surface1,2 are.
Definition: refinementSurfaces.C:1203
Foam::refinementSurfaces::setMinLevelFields
void setMinLevelFields(const shellSurfaces &shells)
Calculate minLevelFields.
Definition: refinementSurfaces.C:718
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Definition: searchableSurface.H:72
Foam::searchableSurface::getRegion
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const =0
From a set of points and indices get the region.
Foam::searchableSurface::findLine
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const =0
Find first intersection on segment from start to end.
Foam::searchableSurface::findNearest
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
Foam::searchableSurface::findLineAny
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const =0
Return any intersection on segment from start to end.
Foam::searchableSurface::globalSize
virtual label globalSize() const
Range of global indices that can be returned.
Definition: searchableSurface.H:205
Foam::searchableSurface::getNormal
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const =0
From a set of points and indices get the normal.
Foam::searchableSurface::boundingSpheres
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const =0
Get bounding spheres (centre and radius squared), one per element.
Foam::searchableSurfacesQueries::findAnyIntersection
static void findAnyIntersection(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &)
Find any intersection. Return hit point information and.
Definition: searchableSurfacesQueries.C:122
Foam::searchableSurfacesQueries::findNearest
static void findNearest(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &)
Find nearest. Return -1 (and a miss()) or surface and nearest.
Definition: searchableSurfacesQueries.C:349
Foam::searchableSurfaces
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
Definition: searchableSurfaces.H:96
Foam::searchableSurfaces::names
const wordList & names() const
Surface names, not region names.
Definition: searchableSurfaces.H:157
Foam::searchableSurfaces::regionNames
const List< wordList > & regionNames() const
Region names per surface.
Definition: searchableSurfaces.H:169
Foam::shellSurfaces
Encapsulates queries for volume refinement ('refine all cells within shell').
Definition: shellSurfaces.H:58
Foam::surfaceLocation
Contains information about location on a triSurface.
Definition: surfaceLocation.H:74
Foam::surfaceZonesInfo::areaSelectionAlgo
areaSelectionAlgo
Types of selection of area.
Definition: surfaceZonesInfo.H:66
Foam::surfaceZonesInfo::areaSelectionAlgoNames
static const Enum< areaSelectionAlgo > areaSelectionAlgoNames
Definition: surfaceZonesInfo.H:73
Foam::triSurface
Triangulated surface description with patch information.
Definition: triSurface.H:79
Foam::volumeType
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
Foam::volumeType::str
const word & str() const
The string representation of the volume type enumeration.
Definition: volumeType.C:62
Foam::volumeType::OUTSIDE
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
Foam::volumeType::MIXED
@ MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
Foam::volumeType::UNKNOWN
@ UNKNOWN
Unknown state.
Definition: volumeType.H:67
Foam::volumeType::INSIDE
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:68
p0
const volScalarField & p0
Definition: EEqn.H:36
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
regionNames
wordList regionNames
Definition: getAllRegionOptions.H:37
points
const pointField & points
Definition: gmvOutputHeader.H:1
IOWarningInFunction
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
Definition: messageStream.H:340
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:66
Foam::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::labelMax
constexpr label labelMax
Definition: label.H:61
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:52
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::FatalIOError
IOerror FatalIOError
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::FatalError
error FatalError
Foam::returnReduce
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
Definition: PstreamReduceOps.H:91
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:130
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
dict
dictionary dict
Definition: searchingEngine.H:14
setField
surfacesMesh setField(triSurfaceToAgglom)
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
forAllConstIters
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
samples
scalarField samples(nIntervals, Zero)