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