refinementSurfaces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "refinementSurfaces.H"
27 #include "Time.H"
28 #include "searchableSurfaces.H"
29 #include "refinementRegions.H"
30 #include "triSurfaceMesh.H"
31 #include "labelPair.H"
33 #include "UPtrList.H"
34 #include "volumeType.H"
35 
36 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
37 
39 (
40  const searchableSurfaces& allGeometry,
41  const dictionary& surfacesDict,
42  const label gapLevelIncrement
43 )
44 :
45  allGeometry_(allGeometry),
46  surfaces_(surfacesDict.size()),
47  names_(surfacesDict.size()),
48  surfZones_(surfacesDict.size()),
49  regionOffset_(surfacesDict.size())
50 {
51  // Wildcard specification : loop over all surface, all regions
52  // and try to find a match.
53 
54  // Count number of surfaces.
55  label surfi = 0;
56  forAll(allGeometry_.names(), geomi)
57  {
58  const word& geomName = allGeometry_.names()[geomi];
59 
60  if (surfacesDict.found(geomName))
61  {
62  surfi++;
63  }
64  }
65 
66  // Size lists
67  surfaces_.setSize(surfi);
68  names_.setSize(surfi);
69  surfZones_.setSize(surfi);
70  regionOffset_.setSize(surfi);
71 
72  labelList globalMinLevel(surfi, 0);
73  labelList globalMaxLevel(surfi, 0);
74  labelList globalLevelIncr(surfi, 0);
75  scalarField globalAngle(surfi, -great);
76  PtrList<dictionary> globalPatchInfo(surfi);
77  List<Map<label>> regionMinLevel(surfi);
78  List<Map<label>> regionMaxLevel(surfi);
79  List<Map<label>> regionLevelIncr(surfi);
80  List<Map<scalar>> regionAngle(surfi);
81  List<Map<autoPtr<dictionary>>> regionPatchInfo(surfi);
82 
83 
84  HashSet<word> unmatchedKeys(surfacesDict.toc());
85 
86  surfi = 0;
87  forAll(allGeometry_.names(), geomi)
88  {
89  const word& geomName = allGeometry_.names()[geomi];
90 
91  const entry* ePtr = surfacesDict.lookupEntryPtr(geomName, false, true);
92 
93  if (ePtr)
94  {
95  const dictionary& dict = ePtr->dict();
96  unmatchedKeys.erase(ePtr->keyword());
97 
98  names_[surfi] = geomName;
99  surfaces_[surfi] = geomi;
100 
101  const labelPair refLevel(dict.lookup("level"));
102  globalMinLevel[surfi] = refLevel[0];
103  globalMaxLevel[surfi] = refLevel[1];
104  globalLevelIncr[surfi] = dict.lookupOrDefault
105  (
106  "gapLevelIncrement",
107  gapLevelIncrement
108  );
109 
110  if
111  (
112  globalMinLevel[surfi] < 0
113  || globalMaxLevel[surfi] < globalMinLevel[surfi]
114  || globalLevelIncr[surfi] < 0
115  )
116  {
118  (
119  dict
120  ) << "Illegal level specification for surface "
121  << names_[surfi]
122  << " : minLevel:" << globalMinLevel[surfi]
123  << " maxLevel:" << globalMaxLevel[surfi]
124  << " levelincrement:" << globalLevelIncr[surfi]
125  << exit(FatalIOError);
126  }
127 
128  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
129 
130  // Surface zones
131  surfZones_.set(surfi, new surfaceZonesInfo(surface, dict));
132 
133  // Global perpendicular angle
134  if (dict.found("patchInfo"))
135  {
136  globalPatchInfo.set
137  (
138  surfi,
139  dict.subDict("patchInfo").clone()
140  );
141  }
142  dict.readIfPresent("perpendicularAngle", globalAngle[surfi]);
143 
144  if (dict.found("regions"))
145  {
146  const dictionary& regionsDict = dict.subDict("regions");
147  const wordList& regionNames = surface.regions();
148 
149  forAll(regionNames, regioni)
150  {
151  if (regionsDict.found(regionNames[regioni]))
152  {
153  // Get the dictionary for region
154  const dictionary& regionDict = regionsDict.subDict
155  (
156  regionNames[regioni]
157  );
158 
159  const labelPair refLevel(regionDict.lookup("level"));
160 
161  regionMinLevel[surfi].insert(regioni, refLevel[0]);
162  regionMaxLevel[surfi].insert(regioni, refLevel[1]);
163  label levelincr = regionDict.lookupOrDefault
164  (
165  "gapLevelIncrement",
166  gapLevelIncrement
167  );
168  regionLevelIncr[surfi].insert(regioni, levelincr);
169 
170  if
171  (
172  refLevel[0] < 0
173  || refLevel[1] < refLevel[0]
174  || levelincr < 0
175  )
176  {
178  (
179  dict
180  ) << "Illegal level specification for surface "
181  << names_[surfi] << " region "
182  << regionNames[regioni]
183  << " : minLevel:" << refLevel[0]
184  << " maxLevel:" << refLevel[1]
185  << " levelincrement:" << levelincr
186  << exit(FatalIOError);
187  }
188 
189  if (regionDict.found("perpendicularAngle"))
190  {
191  regionAngle[surfi].insert
192  (
193  regioni,
194  regionDict.lookup<scalar>("perpendicularAngle")
195  );
196  }
197 
198  if (regionDict.found("patchInfo"))
199  {
200  regionPatchInfo[surfi].insert
201  (
202  regioni,
203  regionDict.subDict("patchInfo").clone()
204  );
205  }
206  }
207  }
208  }
209  surfi++;
210  }
211  }
212 
213  if (unmatchedKeys.size() > 0)
214  {
216  (
217  surfacesDict
218  ) << "Not all entries in refinementSurfaces dictionary were used."
219  << " The following entries were not used : "
220  << unmatchedKeys.sortedToc()
221  << endl;
222  }
223 
224 
225  // Calculate local to global region offset
226  label nRegions = 0;
227 
228  forAll(surfaces_, surfi)
229  {
230  regionOffset_[surfi] = nRegions;
231  nRegions += allGeometry_[surfaces_[surfi]].regions().size();
232  }
233 
234  // Rework surface specific information into information per global region
235  minLevel_.setSize(nRegions);
236  minLevel_ = 0;
237  maxLevel_.setSize(nRegions);
238  maxLevel_ = 0;
239  gapLevel_.setSize(nRegions);
240  gapLevel_ = -1;
241  perpendicularAngle_.setSize(nRegions);
242  perpendicularAngle_ = -great;
243  patchInfo_.setSize(nRegions);
244 
245 
246  forAll(globalMinLevel, surfi)
247  {
248  label nRegions = allGeometry_[surfaces_[surfi]].regions().size();
249 
250  // Initialise to global (i.e. per surface)
251  for (label i = 0; i < nRegions; i++)
252  {
253  label globalRegioni = regionOffset_[surfi] + i;
254  minLevel_[globalRegioni] = globalMinLevel[surfi];
255  maxLevel_[globalRegioni] = globalMaxLevel[surfi];
256  gapLevel_[globalRegioni] =
257  maxLevel_[globalRegioni]
258  + globalLevelIncr[surfi];
259 
260  perpendicularAngle_[globalRegioni] = globalAngle[surfi];
261  if (globalPatchInfo.set(surfi))
262  {
263  patchInfo_.set
264  (
265  globalRegioni,
266  globalPatchInfo[surfi].clone()
267  );
268  }
269  }
270 
271  // Overwrite with region specific information
272  forAllConstIter(Map<label>, regionMinLevel[surfi], iter)
273  {
274  label globalRegioni = regionOffset_[surfi] + iter.key();
275 
276  minLevel_[globalRegioni] = iter();
277  maxLevel_[globalRegioni] = regionMaxLevel[surfi][iter.key()];
278  gapLevel_[globalRegioni] =
279  maxLevel_[globalRegioni]
280  + regionLevelIncr[surfi][iter.key()];
281  }
282  forAllConstIter(Map<scalar>, regionAngle[surfi], iter)
283  {
284  label globalRegioni = regionOffset_[surfi] + iter.key();
285 
286  perpendicularAngle_[globalRegioni] = regionAngle[surfi][iter.key()];
287  }
288 
289  const Map<autoPtr<dictionary>>& localInfo = regionPatchInfo[surfi];
290  forAllConstIter(Map<autoPtr<dictionary>>, localInfo, iter)
291  {
292  label globalRegioni = regionOffset_[surfi] + iter.key();
293 
294  patchInfo_.set(globalRegioni, iter()().clone());
295  }
296  }
297 }
298 
299 
301 (
302  const searchableSurfaces& allGeometry,
303  const labelList& surfaces,
304  const wordList& names,
305  const PtrList<surfaceZonesInfo>& surfZones,
306  const labelList& regionOffset,
307  const labelList& minLevel,
308  const labelList& maxLevel,
309  const labelList& gapLevel,
310  const scalarField& perpendicularAngle,
311  PtrList<dictionary>& patchInfo
312 )
313 :
314  allGeometry_(allGeometry),
315  surfaces_(surfaces),
316  names_(names),
317  surfZones_(surfZones),
318  regionOffset_(regionOffset),
319  minLevel_(minLevel),
320  maxLevel_(maxLevel),
321  gapLevel_(gapLevel),
322  perpendicularAngle_(perpendicularAngle),
323  patchInfo_(patchInfo.size())
324 {
325  forAll(patchInfo_, pi)
326  {
327  if (patchInfo.set(pi))
328  {
329  patchInfo_.set(pi, patchInfo.set(pi, nullptr));
330  }
331  }
332 }
333 
334 
335 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
336 
338 (
339  const refinementRegions& shells,
340  const scalar level0EdgeLength,
341  const bool extendedRefinementSpan
342 )
343 {
344  forAll(surfaces_, surfi)
345  {
346  const searchableSurface& geom = allGeometry_[surfaces_[surfi]];
347 
348  // Precalculation only makes sense if there are different regions
349  // (so different refinement levels possible) and there are some
350  // elements. Possibly should have 'enough' elements to have fine
351  // enough resolution but for now just make sure we don't catch e.g.
352  // searchableBox (size=6)
353  if (geom.regions().size() > 1 && geom.globalSize() > 10)
354  {
355  // Representative local coordinates and bounding sphere
356  pointField ctrs;
357  scalarField radiusSqr;
358  geom.boundingSpheres(ctrs, radiusSqr);
359 
360  labelList minLevelField(ctrs.size(), -1);
361  {
362  // Get the element index in a roundabout way. Problem is e.g.
363  // distributed surface where local indices differ from global
364  // ones (needed for getRegion call)
365  List<pointIndexHit> info;
366  geom.findNearest(ctrs, radiusSqr, info);
367 
368  // Get per element the region
369  labelList region;
370  geom.getRegion(info, region);
371 
372  // From the region get the surface-wise refinement level
373  forAll(minLevelField, i)
374  {
375  if (info[i].hit()) // Note: should not be necessary
376  {
377  minLevelField[i] = minLevel(surfi, region[i]);
378  }
379  }
380  }
381 
382  // Find out if triangle inside shell with higher level
383  // What level does shell want to refine fc to?
384  //
385  // Note: for triangulated surfaces with triangles that
386  // span refinement regions it introduces unnecessary refinement
387  if (extendedRefinementSpan)
388  {
389  labelList shellLevel;
390  shells.findHigherLevel
391  (
392  ctrs,
393  minLevelField,
394  level0EdgeLength,
395  shellLevel
396  );
397 
398  forAll(minLevelField, i)
399  {
400  minLevelField[i] = max(minLevelField[i], shellLevel[i]);
401  }
402  }
403 
404  // Store minLevelField on surface
405  const_cast<searchableSurface&>(geom).setField(minLevelField);
406  }
407  }
408 }
409 
410 
412 (
413  const pointField& start,
414  const pointField& end,
415  const labelList& currentLevel, // current cell refinement level
416 
417  labelList& surfaces,
418  labelList& surfaceLevel
419 ) const
420 {
421  surfaces.setSize(start.size());
422  surfaces = -1;
423  surfaceLevel.setSize(start.size());
424  surfaceLevel = -1;
425 
426  if (surfaces_.empty())
427  {
428  return;
429  }
430 
431  if (surfaces_.size() == 1)
432  {
433  // Optimisation: single segmented surface. No need to duplicate
434  // point storage.
435 
436  label surfi = 0;
437 
438  const searchableSurface& geom = allGeometry_[surfaces_[surfi]];
439 
440  // Do intersection test
441  List<pointIndexHit> intersectionInfo(start.size());
442  geom.findLineAny(start, end, intersectionInfo);
443 
444  // See if a cached level field available
445  labelList minLevelField;
446  geom.getField(intersectionInfo, minLevelField);
447  bool haveLevelField =
448  (
449  returnReduce(minLevelField.size(), sumOp<label>())
450  > 0
451  );
452 
453  if (!haveLevelField && geom.regions().size() == 1)
454  {
455  minLevelField = labelList
456  (
457  intersectionInfo.size(),
458  minLevel(surfi, 0)
459  );
460  haveLevelField = true;
461  }
462 
463  if (haveLevelField)
464  {
465  forAll(intersectionInfo, i)
466  {
467  if
468  (
469  intersectionInfo[i].hit()
470  && minLevelField[i] > currentLevel[i]
471  )
472  {
473  surfaces[i] = surfi; // index of surface
474  surfaceLevel[i] = minLevelField[i];
475  }
476  }
477  return;
478  }
479  }
480 
481 
482  // Work arrays
483  pointField p0(start);
484  pointField p1(end);
485  labelList intersectionToPoint(identity(start.size()));
486  List<pointIndexHit> intersectionInfo(start.size());
487 
488  forAll(surfaces_, surfi)
489  {
490  const searchableSurface& geom = allGeometry_[surfaces_[surfi]];
491 
492  // Do intersection test
493  geom.findLineAny(p0, p1, intersectionInfo);
494 
495  // See if a cached level field available
496  labelList minLevelField;
497  geom.getField(intersectionInfo, minLevelField);
498 
499  // Copy all hits into arguments, In-place compact misses.
500  label missI = 0;
501  forAll(intersectionInfo, i)
502  {
503  // Get the minLevel for the point
504  label minLocalLevel = -1;
505 
506  if (intersectionInfo[i].hit())
507  {
508  // Check if minLevelField for this surface.
509  if (minLevelField.size())
510  {
511  minLocalLevel = minLevelField[i];
512  }
513  else
514  {
515  // Use the min level for the surface instead. Assume
516  // single region 0.
517  minLocalLevel = minLevel(surfi, 0);
518  }
519  }
520 
521 
522  label pointi = intersectionToPoint[i];
523 
524  if (minLocalLevel > currentLevel[pointi])
525  {
526  // Mark point for refinement
527  surfaces[pointi] = surfi;
528  surfaceLevel[pointi] = minLocalLevel;
529  }
530  else
531  {
532  p0[missI] = start[pointi];
533  p1[missI] = end[pointi];
534  intersectionToPoint[missI] = pointi;
535  missI++;
536  }
537  }
538 
539  // All done? Note that this decision should be synchronised
540  if (returnReduce(missI, sumOp<label>()) == 0)
541  {
542  break;
543  }
544 
545  // Trim misses
546  p0.setSize(missI);
547  p1.setSize(missI);
548  intersectionToPoint.setSize(missI);
549  intersectionInfo.setSize(missI);
550  }
551 }
552 
553 
555 (
556  const pointField& start,
557  const pointField& end,
558  const labelList& currentLevel, // current cell refinement level
559 
560  const labelList& globalRegionLevel,
561 
562  List<vectorList>& surfaceNormal,
563  labelListList& surfaceLevel
564 ) const
565 {
566  surfaceLevel.setSize(start.size());
567  surfaceNormal.setSize(start.size());
568 
569  if (surfaces_.empty())
570  {
571  return;
572  }
573 
574  // Work arrays
575  List<List<pointIndexHit>> hitinfo;
576  labelList pRegions;
577  vectorField pNormals;
578 
579  forAll(surfaces_, surfi)
580  {
581  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
582 
583  surface.findLineAll(start, end, hitinfo);
584 
585  // Repack hits for surface into flat list
586  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
587  // To avoid overhead of calling getRegion for every point
588 
589  label n = 0;
590  forAll(hitinfo, pointi)
591  {
592  n += hitinfo[pointi].size();
593  }
594 
595  List<pointIndexHit> surfinfo(n);
596  labelList pointMap(n);
597  n = 0;
598 
599  forAll(hitinfo, pointi)
600  {
601  const List<pointIndexHit>& pHits = hitinfo[pointi];
602 
603  forAll(pHits, i)
604  {
605  surfinfo[n] = pHits[i];
606  pointMap[n] = pointi;
607  n++;
608  }
609  }
610 
611  labelList surfRegion(n);
612  vectorField surfNormal(n);
613  surface.getRegion(surfinfo, surfRegion);
614  surface.getNormal(surfinfo, surfNormal);
615 
616  surfinfo.clear();
617 
618 
619  // Extract back into pointwise
620  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
621 
622  forAll(surfRegion, i)
623  {
624  label region = globalRegion(surfi, surfRegion[i]);
625  label pointi = pointMap[i];
626 
627  if (globalRegionLevel[region] > currentLevel[pointi])
628  {
629  // Append to pointi info
630  label sz = surfaceNormal[pointi].size();
631  surfaceNormal[pointi].setSize(sz+1);
632  surfaceNormal[pointi][sz] = surfNormal[i];
633 
634  surfaceLevel[pointi].setSize(sz+1);
635  surfaceLevel[pointi][sz] = globalRegionLevel[region];
636  }
637  }
638  }
639 }
640 
641 
643 (
644  const pointField& start,
645  const pointField& end,
646  const labelList& currentLevel, // current cell refinement level
647 
648  const labelList& globalRegionLevel,
649 
650  List<pointList>& surfaceLocation,
651  List<vectorList>& surfaceNormal,
652  labelListList& surfaceLevel
653 ) const
654 {
655  surfaceLevel.setSize(start.size());
656  surfaceNormal.setSize(start.size());
657  surfaceLocation.setSize(start.size());
658 
659  if (surfaces_.empty())
660  {
661  return;
662  }
663 
664  // Work arrays
665  List<List<pointIndexHit>> hitinfo;
666  labelList pRegions;
667  vectorField pNormals;
668 
669  forAll(surfaces_, surfi)
670  {
671  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
672 
673  surface.findLineAll(start, end, hitinfo);
674 
675  // Repack hits for surface into flat list
676  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
677  // To avoid overhead of calling getRegion for every point
678 
679  label n = 0;
680  forAll(hitinfo, pointi)
681  {
682  n += hitinfo[pointi].size();
683  }
684 
685  List<pointIndexHit> surfinfo(n);
686  labelList pointMap(n);
687  n = 0;
688 
689  forAll(hitinfo, pointi)
690  {
691  const List<pointIndexHit>& pHits = hitinfo[pointi];
692 
693  forAll(pHits, i)
694  {
695  surfinfo[n] = pHits[i];
696  pointMap[n] = pointi;
697  n++;
698  }
699  }
700 
701  labelList surfRegion(n);
702  vectorField surfNormal(n);
703  surface.getRegion(surfinfo, surfRegion);
704  surface.getNormal(surfinfo, surfNormal);
705 
706  // Extract back into pointwise
707  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
708 
709  forAll(surfRegion, i)
710  {
711  label region = globalRegion(surfi, surfRegion[i]);
712  label pointi = pointMap[i];
713 
714  if (globalRegionLevel[region] > currentLevel[pointi])
715  {
716  // Append to pointi info
717  label sz = surfaceNormal[pointi].size();
718  surfaceLocation[pointi].setSize(sz+1);
719  surfaceLocation[pointi][sz] = surfinfo[i].hitPoint();
720 
721  surfaceNormal[pointi].setSize(sz+1);
722  surfaceNormal[pointi][sz] = surfNormal[i];
723 
724  surfaceLevel[pointi].setSize(sz+1);
725  surfaceLevel[pointi][sz] = globalRegionLevel[region];
726  }
727  }
728  }
729 }
730 
731 
733 (
734  const labelList& surfacesToTest,
735  const pointField& start,
736  const pointField& end,
737 
738  labelList& surface1,
739  List<pointIndexHit>& hit1,
740  labelList& region1,
741  labelList& surface2,
742  List<pointIndexHit>& hit2,
743  labelList& region2
744 ) const
745 {
746  // 1. intersection from start to end
747  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748 
749  // Initialise arguments
750  surface1.setSize(start.size());
751  surface1 = -1;
752  hit1.setSize(start.size());
753  region1.setSize(start.size());
754 
755  // Current end of segment to test.
756  pointField nearest(end);
757  // Work array
758  List<pointIndexHit> nearestInfo(start.size());
759  labelList region;
760 
761  forAll(surfacesToTest, testi)
762  {
763  label surfi = surfacesToTest[testi];
764 
765  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
766 
767  // See if any intersection between start and current nearest
768  surface.findLine
769  (
770  start,
771  nearest,
772  nearestInfo
773  );
774  surface.getRegion
775  (
776  nearestInfo,
777  region
778  );
779 
780  forAll(nearestInfo, pointi)
781  {
782  if (nearestInfo[pointi].hit())
783  {
784  hit1[pointi] = nearestInfo[pointi];
785  surface1[pointi] = surfi;
786  region1[pointi] = region[pointi];
787  nearest[pointi] = hit1[pointi].hitPoint();
788  }
789  }
790  }
791 
792 
793  // 2. intersection from end to last intersection
794  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
795 
796  // Find the nearest intersection from end to start. Note that we initialise
797  // to the first intersection (if any).
798  surface2 = surface1;
799  hit2 = hit1;
800  region2 = region1;
801 
802  // Set current end of segment to test.
803  forAll(nearest, pointi)
804  {
805  if (hit1[pointi].hit())
806  {
807  nearest[pointi] = hit1[pointi].hitPoint();
808  }
809  else
810  {
811  // Disable testing by setting to end.
812  nearest[pointi] = end[pointi];
813  }
814  }
815 
816  forAll(surfacesToTest, testi)
817  {
818  label surfi = surfacesToTest[testi];
819 
820  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
821 
822  // See if any intersection between end and current nearest
823  surface.findLine
824  (
825  end,
826  nearest,
827  nearestInfo
828  );
829  surface.getRegion
830  (
831  nearestInfo,
832  region
833  );
834 
835  forAll(nearestInfo, pointi)
836  {
837  if (nearestInfo[pointi].hit())
838  {
839  hit2[pointi] = nearestInfo[pointi];
840  surface2[pointi] = surfi;
841  region2[pointi] = region[pointi];
842  nearest[pointi] = hit2[pointi].hitPoint();
843  }
844  }
845  }
846 
847 
848  // Make sure that if hit1 has hit something, hit2 will have at least the
849  // same point (due to tolerances it might miss its end point)
850  forAll(hit1, pointi)
851  {
852  if (hit1[pointi].hit() && !hit2[pointi].hit())
853  {
854  hit2[pointi] = hit1[pointi];
855  surface2[pointi] = surface1[pointi];
856  region2[pointi] = region1[pointi];
857  }
858  }
859 }
860 
861 
863 (
864  const labelList& surfacesToTest,
865  const pointField& start,
866  const pointField& end,
867 
868  labelList& surface1,
869  List<pointIndexHit>& hit1,
870  labelList& region1,
871  vectorField& normal1,
872 
873  labelList& surface2,
874  List<pointIndexHit>& hit2,
875  labelList& region2,
876  vectorField& normal2
877 ) const
878 {
879  // 1. intersection from start to end
880  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
881 
882  // Initialise arguments
883  surface1.setSize(start.size());
884  surface1 = -1;
885  hit1.setSize(start.size());
886  region1.setSize(start.size());
887  region1 = -1;
888  normal1.setSize(start.size());
889  normal1 = Zero;
890 
891  // Current end of segment to test.
892  pointField nearest(end);
893  // Work array
894  List<pointIndexHit> nearestInfo(start.size());
895  labelList region;
896  vectorField normal;
897 
898  forAll(surfacesToTest, testi)
899  {
900  label surfi = surfacesToTest[testi];
901  const searchableSurface& geom = allGeometry_[surfaces_[surfi]];
902 
903  // See if any intersection between start and current nearest
904  geom.findLine(start, nearest, nearestInfo);
905  geom.getRegion(nearestInfo, region);
906  geom.getNormal(nearestInfo, normal);
907 
908  forAll(nearestInfo, pointi)
909  {
910  if (nearestInfo[pointi].hit())
911  {
912  hit1[pointi] = nearestInfo[pointi];
913  surface1[pointi] = surfi;
914  region1[pointi] = region[pointi];
915  normal1[pointi] = normal[pointi];
916  nearest[pointi] = hit1[pointi].hitPoint();
917  }
918  }
919  }
920 
921 
922  // 2. intersection from end to last intersection
923  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
924 
925  // Find the nearest intersection from end to start. Note that we initialise
926  // to the first intersection (if any).
927  surface2 = surface1;
928  hit2 = hit1;
929  region2 = region1;
930  normal2 = normal1;
931 
932  // Set current end of segment to test.
933  forAll(nearest, pointi)
934  {
935  if (hit1[pointi].hit())
936  {
937  nearest[pointi] = hit1[pointi].hitPoint();
938  }
939  else
940  {
941  // Disable testing by setting to end.
942  nearest[pointi] = end[pointi];
943  }
944  }
945 
946  forAll(surfacesToTest, testi)
947  {
948  label surfi = surfacesToTest[testi];
949  const searchableSurface& geom = allGeometry_[surfaces_[surfi]];
950 
951  // See if any intersection between end and current nearest
952  geom.findLine(end, nearest, nearestInfo);
953  geom.getRegion(nearestInfo, region);
954  geom.getNormal(nearestInfo, normal);
955 
956  forAll(nearestInfo, pointi)
957  {
958  if (nearestInfo[pointi].hit())
959  {
960  hit2[pointi] = nearestInfo[pointi];
961  surface2[pointi] = surfi;
962  region2[pointi] = region[pointi];
963  normal2[pointi] = normal[pointi];
964  nearest[pointi] = hit2[pointi].hitPoint();
965  }
966  }
967  }
968 
969 
970  // Make sure that if hit1 has hit something, hit2 will have at least the
971  // same point (due to tolerances it might miss its end point)
972  forAll(hit1, pointi)
973  {
974  if (hit1[pointi].hit() && !hit2[pointi].hit())
975  {
976  hit2[pointi] = hit1[pointi];
977  surface2[pointi] = surface1[pointi];
978  region2[pointi] = region1[pointi];
979  normal2[pointi] = normal1[pointi];
980  }
981  }
982 }
983 
984 
986 (
987  const pointField& start,
988  const pointField& end,
989 
990  labelList& hitSurface,
991  List<pointIndexHit>& hitinfo
992 ) const
993 {
994  searchableSurfacesQueries::findAnyIntersection
995  (
996  allGeometry_,
997  surfaces_,
998  start,
999  end,
1000  hitSurface,
1001  hitinfo
1002  );
1003 }
1004 
1005 
1008  const labelList& surfacesToTest,
1009  const pointField& samples,
1010  const scalarField& nearestDistSqr,
1011  labelList& hitSurface,
1012  List<pointIndexHit>& hitinfo
1013 ) const
1014 {
1015  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1016 
1017  // Do the tests. Note that findNearest returns index in geometries.
1018  searchableSurfacesQueries::findNearest
1019  (
1020  allGeometry_,
1021  geometries,
1022  samples,
1023  nearestDistSqr,
1024  hitSurface,
1025  hitinfo
1026  );
1027 
1028  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1029  forAll(hitSurface, pointi)
1030  {
1031  if (hitSurface[pointi] != -1)
1032  {
1033  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1034  }
1035  }
1036 }
1037 
1038 
1041  const labelList& surfacesToTest,
1042  const pointField& samples,
1043  const scalarField& nearestDistSqr,
1044  labelList& hitSurface,
1045  labelList& hitRegion
1046 ) const
1047 {
1048  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1049 
1050  // Do the tests. Note that findNearest returns index in geometries.
1051  List<pointIndexHit> hitinfo;
1052  searchableSurfacesQueries::findNearest
1053  (
1054  allGeometry_,
1055  geometries,
1056  samples,
1057  nearestDistSqr,
1058  hitSurface,
1059  hitinfo
1060  );
1061 
1062  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1063  forAll(hitSurface, pointi)
1064  {
1065  if (hitSurface[pointi] != -1)
1066  {
1067  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1068  }
1069  }
1070 
1071  // Collect the region
1072  hitRegion.setSize(hitSurface.size());
1073  hitRegion = -1;
1074 
1075  forAll(surfacesToTest, i)
1076  {
1077  label surfi = surfacesToTest[i];
1078 
1079  // Collect hits for surfi
1080  const labelList localIndices(findIndices(hitSurface, surfi));
1081 
1082  List<pointIndexHit> localHits
1083  (
1085  (
1086  hitinfo,
1087  localIndices
1088  )
1089  );
1090 
1091  labelList localRegion;
1092  allGeometry_[surfaces_[surfi]].getRegion(localHits, localRegion);
1093 
1094  forAll(localIndices, i)
1095  {
1096  hitRegion[localIndices[i]] = localRegion[i];
1097  }
1098  }
1099 }
1100 
1101 
1104  const labelList& surfacesToTest,
1105  const pointField& samples,
1106  const scalarField& nearestDistSqr,
1107  labelList& hitSurface,
1108  List<pointIndexHit>& hitinfo,
1109  labelList& hitRegion,
1110  vectorField& hitNormal
1111 ) const
1112 {
1113  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1114 
1115  // Do the tests. Note that findNearest returns index in geometries.
1116  searchableSurfacesQueries::findNearest
1117  (
1118  allGeometry_,
1119  geometries,
1120  samples,
1121  nearestDistSqr,
1122  hitSurface,
1123  hitinfo
1124  );
1125 
1126  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1127  forAll(hitSurface, pointi)
1128  {
1129  if (hitSurface[pointi] != -1)
1130  {
1131  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1132  }
1133  }
1134 
1135  // Collect the region
1136  hitRegion.setSize(hitSurface.size());
1137  hitRegion = -1;
1138  hitNormal.setSize(hitSurface.size());
1139  hitNormal = Zero;
1140 
1141  forAll(surfacesToTest, i)
1142  {
1143  label surfi = surfacesToTest[i];
1144 
1145  // Collect hits for surfi
1146  const labelList localIndices(findIndices(hitSurface, surfi));
1147 
1148  List<pointIndexHit> localHits
1149  (
1151  (
1152  hitinfo,
1153  localIndices
1154  )
1155  );
1156 
1157  // Region
1158  labelList localRegion;
1159  allGeometry_[surfaces_[surfi]].getRegion(localHits, localRegion);
1160 
1161  forAll(localIndices, i)
1162  {
1163  hitRegion[localIndices[i]] = localRegion[i];
1164  }
1165 
1166  // Normal
1167  vectorField localNormal;
1168  allGeometry_[surfaces_[surfi]].getNormal(localHits, localNormal);
1169 
1170  forAll(localIndices, i)
1171  {
1172  hitNormal[localIndices[i]] = localNormal[i];
1173  }
1174  }
1175 }
1176 
1177 
1180  const labelList& testSurfaces,
1181  const pointField& pt,
1182  labelList& insideSurfaces
1183 ) const
1184 {
1185  insideSurfaces.setSize(pt.size());
1186  insideSurfaces = -1;
1187 
1188  forAll(testSurfaces, i)
1189  {
1190  label surfi = testSurfaces[i];
1191 
1192  const searchableSurface& surface = allGeometry_[surfaces_[surfi]];
1193 
1194  const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
1195  surfZones_[surfi].zoneInside();
1196 
1197  if
1198  (
1199  selectionMethod != surfaceZonesInfo::INSIDE
1200  && selectionMethod != surfaceZonesInfo::OUTSIDE
1201  )
1202  {
1204  << "Trying to use surface "
1205  << surface.name()
1206  << " which has non-geometric inside selection method "
1207  << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
1208  << exit(FatalError);
1209  }
1210 
1211  if (surface.hasVolumeType())
1212  {
1213  List<volumeType> volType;
1214  surface.getVolumeType(pt, volType);
1215 
1216  forAll(volType, pointi)
1217  {
1218  if (insideSurfaces[pointi] == -1)
1219  {
1220  if
1221  (
1222  (
1223  volType[pointi] == volumeType::inside
1224  && selectionMethod == surfaceZonesInfo::INSIDE
1225  )
1226  || (
1227  volType[pointi] == volumeType::outside
1228  && selectionMethod == surfaceZonesInfo::OUTSIDE
1229  )
1230  )
1231  {
1232  insideSurfaces[pointi] = surfi;
1233  }
1234  }
1235  }
1236  }
1237  }
1238 }
1239 
1240 
1241 // ************************************************************************* //
A HashTable with keys but without contents.
Definition: HashSet.H:59
dictionary dict
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:663
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const word & name() const
Return name.
Definition: IOobject.H:315
const entry * lookupEntryPtr(const word &, bool recursive, bool patternMatch) const
Find and return an entry data stream pointer if present.
Definition: dictionary.C:702
virtual const wordList & regions() const =0
Names of regions.
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
const keyType & keyword() const
Return keyword.
Definition: entry.H:123
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
void findAllHigherIntersections(const pointField &start, const pointField &end, const labelList &currentLevel, const labelList &globalRegionLevel, List< vectorList > &surfaceNormal, labelListList &surfaceLevel) const
Find all intersections of edge. Unsorted order.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const dictionary & dict() const =0
Return dictionary if this entry is a dictionary.
void findAnyIntersection(const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &) const
Used for debugging only: find intersection of edge.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Encapsulates queries for volume refinement (&#39;refine all cells within shell&#39;).
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const =0
Return any intersection on segment from start to end.
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
virtual void getField(const List< pointIndexHit > &, labelList &values) const
WIP. From a set of hits (points and.
void findNearest(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest point on surfaces.
bool erase(T *p)
Remove the specified element from the list and delete it.
Definition: ILList.C:105
autoPtr< dictionary > clone() const
Construct and return clone.
Definition: dictionary.C:552
wordList toc() const
Return the table of contents.
Definition: dictionary.C:1115
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:1002
T clone(const T &t)
Definition: List.H:54
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
surfacesMesh setField(triSurfaceToAgglom)
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.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit >> &) const =0
Get all intersections in order from start to end.
A class for handling words, derived from string.
Definition: word.H:59
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
Container for searchableSurfaces.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual label globalSize() const
Range of global indices that can be returned.
void setMinLevelFields(const refinementRegions &shells, const scalar level0EdgeLength, const bool extendedRefinementSpan)
Calculate the refinement level for every element.
refinementSurfaces(const searchableSurfaces &allGeometry, const dictionary &, const label gapLevelIncrement)
Construct from surfaces and dictionary.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setSize(const label)
Reset size of List.
Definition: List.C:281
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
A List with indirect addressing.
Definition: fvMatrix.H:106
virtual void getVolumeType(const pointField &, List< volumeType > &) const =0
Determine type (inside/outside) for point. unknown if.
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const =0
From a set of points and indices get the region.
label n
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
void findHigherIntersection(const pointField &start, const pointField &end, const labelList &currentLevel, labelList &surfaces, labelList &surfaceLevel) const
Find intersection of edge. Return -1 or first surface.
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.
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const =0
From a set of points and indices get the normal.
areaSelectionAlgo
Types of selection of area.
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const =0
Find first intersection on segment from start to end.
virtual bool hasVolumeType() const =0
Whether supports volume type below.
A keyword and a list of tokens is an &#39;entry&#39;.
Definition: entry.H:65
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
void findInside(const labelList &surfacesToTest, const pointField &pt, labelList &insideSurfaces) const
Detect if a point is &#39;inside&#39; (closed) surfaces.
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const =0
Get bounding spheres (centre and radius squared), one per element.
IOerror FatalIOError