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