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 
337 // // Count number of triangles per surface region
338 // Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
339 // {
340 // const geometricSurfacePatchList& regions = s.patches();
341 //
342 // labelList nTris(regions.size(), 0);
343 //
344 // forAll(s, triI)
345 // {
346 // nTris[s[triI].region()]++;
347 // }
348 // return nTris;
349 // }
350 
351 
352 // // Pre-calculate the refinement level for every element
353 // void Foam::refinementSurfaces::wantedRefinementLevel
354 // (
355 // const refinementRegions& shells,
356 // const label surfI,
357 // const List<pointIndexHit>& info, // Indices
358 // const pointField& ctrs, // Representative coordinate
359 // labelList& minLevelField
360 // ) const
361 // {
362 // const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
363 //
364 // // Get per element the region
365 // labelList region;
366 // geom.getRegion(info, region);
367 //
368 // // Initialise fields to region wise minLevel
369 // minLevelField.setSize(ctrs.size());
370 // minLevelField = -1;
371 //
372 // forAll(minLevelField, i)
373 // {
374 // if (info[i].hit())
375 // {
376 // minLevelField[i] = minLevel(surfI, region[i]);
377 // }
378 // }
379 //
380 // // Find out if triangle inside shell with higher level
381 // // What level does shell want to refine fc to?
382 // labelList shellLevel;
383 // shells.findHigherLevel(ctrs, minLevelField, shellLevel);
384 //
385 // forAll(minLevelField, i)
386 // {
387 // minLevelField[i] = max(minLevelField[i], shellLevel[i]);
388 // }
389 // }
390 
391 
392 // Precalculate the refinement level for every element of the searchable
393 // surface.
395 (
396  const refinementRegions& shells
397 )
398 {
399  forAll(surfaces_, surfI)
400  {
401  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
402 
403  // Precalculation only makes sense if there are different regions
404  // (so different refinement levels possible) and there are some
405  // elements. Possibly should have 'enough' elements to have fine
406  // enough resolution but for now just make sure we don't catch e.g.
407  // searchableBox (size=6)
408  if (geom.regions().size() > 1 && geom.globalSize() > 10)
409  {
410  // Representative local coordinates and bounding sphere
411  pointField ctrs;
412  scalarField radiusSqr;
413  geom.boundingSpheres(ctrs, radiusSqr);
414 
415  labelList minLevelField(ctrs.size(), -1);
416  {
417  // Get the element index in a roundabout way. Problem is e.g.
418  // distributed surface where local indices differ from global
419  // ones (needed for getRegion call)
420  List<pointIndexHit> info;
421  geom.findNearest(ctrs, radiusSqr, info);
422 
423  // Get per element the region
424  labelList region;
425  geom.getRegion(info, region);
426 
427  // From the region get the surface-wise refinement level
428  forAll(minLevelField, i)
429  {
430  if (info[i].hit()) // Note: should not be necessary
431  {
432  minLevelField[i] = minLevel(surfI, region[i]);
433  }
434  }
435  }
436 
437  // Find out if triangle inside shell with higher level
438  // What level does shell want to refine fc to?
439  //
440  // Note: it is not clear for what cases this additional requirement
441  // is beneficial but for triangulated surfaces with triangles that
442  // span refinement regions it introduces unnecessary refinement so
443  // it has been removed.
444  //
445  // This option can be reinstated under a switch if cases are
446  // provided which demonstrate the benefit.
447  /*
448  labelList shellLevel;
449  shells.findHigherLevel(ctrs, minLevelField, shellLevel);
450 
451  forAll(minLevelField, i)
452  {
453  minLevelField[i] = max(minLevelField[i], shellLevel[i]);
454  }
455  */
456 
457  // Store minLevelField on surface
458  const_cast<searchableSurface&>(geom).setField(minLevelField);
459  }
460  }
461 }
462 
463 
464 // Find intersections of edge. Return -1 or first surface with higher minLevel
465 // number.
467 (
468  const pointField& start,
469  const pointField& end,
470  const labelList& currentLevel, // current cell refinement level
471 
472  labelList& surfaces,
473  labelList& surfaceLevel
474 ) const
475 {
476  surfaces.setSize(start.size());
477  surfaces = -1;
478  surfaceLevel.setSize(start.size());
479  surfaceLevel = -1;
480 
481  if (surfaces_.empty())
482  {
483  return;
484  }
485 
486  if (surfaces_.size() == 1)
487  {
488  // Optimisation: single segmented surface. No need to duplicate
489  // point storage.
490 
491  label surfI = 0;
492 
493  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
494 
495  // Do intersection test
496  List<pointIndexHit> intersectionInfo(start.size());
497  geom.findLineAny(start, end, intersectionInfo);
498 
499  // See if a cached level field available
500  labelList minLevelField;
501  geom.getField(intersectionInfo, minLevelField);
502  bool haveLevelField =
503  (
504  returnReduce(minLevelField.size(), sumOp<label>())
505  > 0
506  );
507 
508  if (!haveLevelField && geom.regions().size() == 1)
509  {
510  minLevelField = labelList
511  (
512  intersectionInfo.size(),
513  minLevel(surfI, 0)
514  );
515  haveLevelField = true;
516  }
517 
518  if (haveLevelField)
519  {
520  forAll(intersectionInfo, i)
521  {
522  if
523  (
524  intersectionInfo[i].hit()
525  && minLevelField[i] > currentLevel[i]
526  )
527  {
528  surfaces[i] = surfI; // index of surface
529  surfaceLevel[i] = minLevelField[i];
530  }
531  }
532  return;
533  }
534  }
535 
536 
537 
538  // Work arrays
539  pointField p0(start);
540  pointField p1(end);
541  labelList intersectionToPoint(identity(start.size()));
542  List<pointIndexHit> intersectionInfo(start.size());
543 
544  forAll(surfaces_, surfI)
545  {
546  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
547 
548  // Do intersection test
549  geom.findLineAny(p0, p1, intersectionInfo);
550 
551  // See if a cached level field available
552  labelList minLevelField;
553  geom.getField(intersectionInfo, minLevelField);
554 
555  // Copy all hits into arguments, In-place compact misses.
556  label missI = 0;
557  forAll(intersectionInfo, i)
558  {
559  // Get the minLevel for the point
560  label minLocalLevel = -1;
561 
562  if (intersectionInfo[i].hit())
563  {
564  // Check if minLevelField for this surface.
565  if (minLevelField.size())
566  {
567  minLocalLevel = minLevelField[i];
568  }
569  else
570  {
571  // Use the min level for the surface instead. Assume
572  // single region 0.
573  minLocalLevel = minLevel(surfI, 0);
574  }
575  }
576 
577 
578  label pointi = intersectionToPoint[i];
579 
580  if (minLocalLevel > currentLevel[pointi])
581  {
582  // Mark point for refinement
583  surfaces[pointi] = surfI;
584  surfaceLevel[pointi] = minLocalLevel;
585  }
586  else
587  {
588  p0[missI] = start[pointi];
589  p1[missI] = end[pointi];
590  intersectionToPoint[missI] = pointi;
591  missI++;
592  }
593  }
594 
595  // All done? Note that this decision should be synchronised
596  if (returnReduce(missI, sumOp<label>()) == 0)
597  {
598  break;
599  }
600 
601  // Trim misses
602  p0.setSize(missI);
603  p1.setSize(missI);
604  intersectionToPoint.setSize(missI);
605  intersectionInfo.setSize(missI);
606  }
607 }
608 
609 
611 (
612  const pointField& start,
613  const pointField& end,
614  const labelList& currentLevel, // current cell refinement level
615 
616  const labelList& globalRegionLevel,
617 
618  List<vectorList>& surfaceNormal,
619  labelListList& surfaceLevel
620 ) const
621 {
622  surfaceLevel.setSize(start.size());
623  surfaceNormal.setSize(start.size());
624 
625  if (surfaces_.empty())
626  {
627  return;
628  }
629 
630  // Work arrays
631  List<List<pointIndexHit>> hitInfo;
632  labelList pRegions;
633  vectorField pNormals;
634 
635  forAll(surfaces_, surfI)
636  {
637  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
638 
639  surface.findLineAll(start, end, hitInfo);
640 
641  // Repack hits for surface into flat list
642  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
643  // To avoid overhead of calling getRegion for every point
644 
645  label n = 0;
646  forAll(hitInfo, pointi)
647  {
648  n += hitInfo[pointi].size();
649  }
650 
651  List<pointIndexHit> surfInfo(n);
652  labelList pointMap(n);
653  n = 0;
654 
655  forAll(hitInfo, pointi)
656  {
657  const List<pointIndexHit>& pHits = hitInfo[pointi];
658 
659  forAll(pHits, i)
660  {
661  surfInfo[n] = pHits[i];
662  pointMap[n] = pointi;
663  n++;
664  }
665  }
666 
667  labelList surfRegion(n);
668  vectorField surfNormal(n);
669  surface.getRegion(surfInfo, surfRegion);
670  surface.getNormal(surfInfo, surfNormal);
671 
672  surfInfo.clear();
673 
674 
675  // Extract back into pointwise
676  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
677 
678  forAll(surfRegion, i)
679  {
680  label region = globalRegion(surfI, surfRegion[i]);
681  label pointi = pointMap[i];
682 
683  if (globalRegionLevel[region] > currentLevel[pointi])
684  {
685  // Append to pointi info
686  label sz = surfaceNormal[pointi].size();
687  surfaceNormal[pointi].setSize(sz+1);
688  surfaceNormal[pointi][sz] = surfNormal[i];
689 
690  surfaceLevel[pointi].setSize(sz+1);
691  surfaceLevel[pointi][sz] = globalRegionLevel[region];
692  }
693  }
694  }
695 }
696 
697 
699 (
700  const pointField& start,
701  const pointField& end,
702  const labelList& currentLevel, // current cell refinement level
703 
704  const labelList& globalRegionLevel,
705 
706  List<pointList>& surfaceLocation,
707  List<vectorList>& surfaceNormal,
708  labelListList& surfaceLevel
709 ) const
710 {
711  surfaceLevel.setSize(start.size());
712  surfaceNormal.setSize(start.size());
713  surfaceLocation.setSize(start.size());
714 
715  if (surfaces_.empty())
716  {
717  return;
718  }
719 
720  // Work arrays
721  List<List<pointIndexHit>> hitInfo;
722  labelList pRegions;
723  vectorField pNormals;
724 
725  forAll(surfaces_, surfI)
726  {
727  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
728 
729  surface.findLineAll(start, end, hitInfo);
730 
731  // Repack hits for surface into flat list
732  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
733  // To avoid overhead of calling getRegion for every point
734 
735  label n = 0;
736  forAll(hitInfo, pointi)
737  {
738  n += hitInfo[pointi].size();
739  }
740 
741  List<pointIndexHit> surfInfo(n);
742  labelList pointMap(n);
743  n = 0;
744 
745  forAll(hitInfo, pointi)
746  {
747  const List<pointIndexHit>& pHits = hitInfo[pointi];
748 
749  forAll(pHits, i)
750  {
751  surfInfo[n] = pHits[i];
752  pointMap[n] = pointi;
753  n++;
754  }
755  }
756 
757  labelList surfRegion(n);
758  vectorField surfNormal(n);
759  surface.getRegion(surfInfo, surfRegion);
760  surface.getNormal(surfInfo, surfNormal);
761 
762  // Extract back into pointwise
763  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
764 
765  forAll(surfRegion, i)
766  {
767  label region = globalRegion(surfI, surfRegion[i]);
768  label pointi = pointMap[i];
769 
770  if (globalRegionLevel[region] > currentLevel[pointi])
771  {
772  // Append to pointi info
773  label sz = surfaceNormal[pointi].size();
774  surfaceLocation[pointi].setSize(sz+1);
775  surfaceLocation[pointi][sz] = surfInfo[i].hitPoint();
776 
777  surfaceNormal[pointi].setSize(sz+1);
778  surfaceNormal[pointi][sz] = surfNormal[i];
779 
780  surfaceLevel[pointi].setSize(sz+1);
781  surfaceLevel[pointi][sz] = globalRegionLevel[region];
782  }
783  }
784  }
785 }
786 
787 
789 (
790  const labelList& surfacesToTest,
791  const pointField& start,
792  const pointField& end,
793 
794  labelList& surface1,
795  List<pointIndexHit>& hit1,
796  labelList& region1,
797  labelList& surface2,
798  List<pointIndexHit>& hit2,
799  labelList& region2
800 ) const
801 {
802  // 1. intersection from start to end
803  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804 
805  // Initialise arguments
806  surface1.setSize(start.size());
807  surface1 = -1;
808  hit1.setSize(start.size());
809  region1.setSize(start.size());
810 
811  // Current end of segment to test.
812  pointField nearest(end);
813  // Work array
814  List<pointIndexHit> nearestInfo(start.size());
815  labelList region;
816 
817  forAll(surfacesToTest, testI)
818  {
819  label surfI = surfacesToTest[testI];
820 
821  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
822 
823  // See if any intersection between start and current nearest
824  surface.findLine
825  (
826  start,
827  nearest,
828  nearestInfo
829  );
830  surface.getRegion
831  (
832  nearestInfo,
833  region
834  );
835 
836  forAll(nearestInfo, pointi)
837  {
838  if (nearestInfo[pointi].hit())
839  {
840  hit1[pointi] = nearestInfo[pointi];
841  surface1[pointi] = surfI;
842  region1[pointi] = region[pointi];
843  nearest[pointi] = hit1[pointi].hitPoint();
844  }
845  }
846  }
847 
848 
849  // 2. intersection from end to last intersection
850  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
851 
852  // Find the nearest intersection from end to start. Note that we initialise
853  // to the first intersection (if any).
854  surface2 = surface1;
855  hit2 = hit1;
856  region2 = region1;
857 
858  // Set current end of segment to test.
859  forAll(nearest, pointi)
860  {
861  if (hit1[pointi].hit())
862  {
863  nearest[pointi] = hit1[pointi].hitPoint();
864  }
865  else
866  {
867  // Disable testing by setting to end.
868  nearest[pointi] = end[pointi];
869  }
870  }
871 
872  forAll(surfacesToTest, testI)
873  {
874  label surfI = surfacesToTest[testI];
875 
876  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
877 
878  // See if any intersection between end and current nearest
879  surface.findLine
880  (
881  end,
882  nearest,
883  nearestInfo
884  );
885  surface.getRegion
886  (
887  nearestInfo,
888  region
889  );
890 
891  forAll(nearestInfo, pointi)
892  {
893  if (nearestInfo[pointi].hit())
894  {
895  hit2[pointi] = nearestInfo[pointi];
896  surface2[pointi] = surfI;
897  region2[pointi] = region[pointi];
898  nearest[pointi] = hit2[pointi].hitPoint();
899  }
900  }
901  }
902 
903 
904  // Make sure that if hit1 has hit something, hit2 will have at least the
905  // same point (due to tolerances it might miss its end point)
906  forAll(hit1, pointi)
907  {
908  if (hit1[pointi].hit() && !hit2[pointi].hit())
909  {
910  hit2[pointi] = hit1[pointi];
911  surface2[pointi] = surface1[pointi];
912  region2[pointi] = region1[pointi];
913  }
914  }
915 }
916 
917 
919 (
920  const labelList& surfacesToTest,
921  const pointField& start,
922  const pointField& end,
923 
924  labelList& surface1,
925  List<pointIndexHit>& hit1,
926  labelList& region1,
927  vectorField& normal1,
928 
929  labelList& surface2,
930  List<pointIndexHit>& hit2,
931  labelList& region2,
932  vectorField& normal2
933 ) const
934 {
935  // 1. intersection from start to end
936  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
937 
938  // Initialise arguments
939  surface1.setSize(start.size());
940  surface1 = -1;
941  hit1.setSize(start.size());
942  region1.setSize(start.size());
943  region1 = -1;
944  normal1.setSize(start.size());
945  normal1 = Zero;
946 
947  // Current end of segment to test.
948  pointField nearest(end);
949  // Work array
950  List<pointIndexHit> nearestInfo(start.size());
951  labelList region;
952  vectorField normal;
953 
954  forAll(surfacesToTest, testI)
955  {
956  label surfI = surfacesToTest[testI];
957  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
958 
959  // See if any intersection between start and current nearest
960  geom.findLine(start, nearest, nearestInfo);
961  geom.getRegion(nearestInfo, region);
962  geom.getNormal(nearestInfo, normal);
963 
964  forAll(nearestInfo, pointi)
965  {
966  if (nearestInfo[pointi].hit())
967  {
968  hit1[pointi] = nearestInfo[pointi];
969  surface1[pointi] = surfI;
970  region1[pointi] = region[pointi];
971  normal1[pointi] = normal[pointi];
972  nearest[pointi] = hit1[pointi].hitPoint();
973  }
974  }
975  }
976 
977 
978  // 2. intersection from end to last intersection
979  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
980 
981  // Find the nearest intersection from end to start. Note that we initialise
982  // to the first intersection (if any).
983  surface2 = surface1;
984  hit2 = hit1;
985  region2 = region1;
986  normal2 = normal1;
987 
988  // Set current end of segment to test.
989  forAll(nearest, pointi)
990  {
991  if (hit1[pointi].hit())
992  {
993  nearest[pointi] = hit1[pointi].hitPoint();
994  }
995  else
996  {
997  // Disable testing by setting to end.
998  nearest[pointi] = end[pointi];
999  }
1000  }
1001 
1002  forAll(surfacesToTest, testI)
1003  {
1004  label surfI = surfacesToTest[testI];
1005  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
1006 
1007  // See if any intersection between end and current nearest
1008  geom.findLine(end, nearest, nearestInfo);
1009  geom.getRegion(nearestInfo, region);
1010  geom.getNormal(nearestInfo, normal);
1011 
1012  forAll(nearestInfo, pointi)
1013  {
1014  if (nearestInfo[pointi].hit())
1015  {
1016  hit2[pointi] = nearestInfo[pointi];
1017  surface2[pointi] = surfI;
1018  region2[pointi] = region[pointi];
1019  normal2[pointi] = normal[pointi];
1020  nearest[pointi] = hit2[pointi].hitPoint();
1021  }
1022  }
1023  }
1024 
1025 
1026  // Make sure that if hit1 has hit something, hit2 will have at least the
1027  // same point (due to tolerances it might miss its end point)
1028  forAll(hit1, pointi)
1029  {
1030  if (hit1[pointi].hit() && !hit2[pointi].hit())
1031  {
1032  hit2[pointi] = hit1[pointi];
1033  surface2[pointi] = surface1[pointi];
1034  region2[pointi] = region1[pointi];
1035  normal2[pointi] = normal1[pointi];
1036  }
1037  }
1038 }
1039 
1040 
1043  const pointField& start,
1044  const pointField& end,
1045 
1046  labelList& hitSurface,
1047  List<pointIndexHit>& hitInfo
1048 ) const
1049 {
1050  searchableSurfacesQueries::findAnyIntersection
1051  (
1052  allGeometry_,
1053  surfaces_,
1054  start,
1055  end,
1056  hitSurface,
1057  hitInfo
1058  );
1059 }
1060 
1061 
1064  const labelList& surfacesToTest,
1065  const pointField& samples,
1066  const scalarField& nearestDistSqr,
1067  labelList& hitSurface,
1068  List<pointIndexHit>& hitInfo
1069 ) const
1070 {
1071  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1072 
1073  // Do the tests. Note that findNearest returns index in geometries.
1074  searchableSurfacesQueries::findNearest
1075  (
1076  allGeometry_,
1077  geometries,
1078  samples,
1079  nearestDistSqr,
1080  hitSurface,
1081  hitInfo
1082  );
1083 
1084  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1085  forAll(hitSurface, pointi)
1086  {
1087  if (hitSurface[pointi] != -1)
1088  {
1089  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1090  }
1091  }
1092 }
1093 
1094 
1097  const labelList& surfacesToTest,
1098  const pointField& samples,
1099  const scalarField& nearestDistSqr,
1100  labelList& hitSurface,
1101  labelList& hitRegion
1102 ) const
1103 {
1104  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1105 
1106  // Do the tests. Note that findNearest returns index in geometries.
1107  List<pointIndexHit> hitInfo;
1108  searchableSurfacesQueries::findNearest
1109  (
1110  allGeometry_,
1111  geometries,
1112  samples,
1113  nearestDistSqr,
1114  hitSurface,
1115  hitInfo
1116  );
1117 
1118  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1119  forAll(hitSurface, pointi)
1120  {
1121  if (hitSurface[pointi] != -1)
1122  {
1123  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1124  }
1125  }
1126 
1127  // Collect the region
1128  hitRegion.setSize(hitSurface.size());
1129  hitRegion = -1;
1130 
1131  forAll(surfacesToTest, i)
1132  {
1133  label surfI = surfacesToTest[i];
1134 
1135  // Collect hits for surfI
1136  const labelList localIndices(findIndices(hitSurface, surfI));
1137 
1138  List<pointIndexHit> localHits
1139  (
1141  (
1142  hitInfo,
1143  localIndices
1144  )
1145  );
1146 
1147  labelList localRegion;
1148  allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1149 
1150  forAll(localIndices, i)
1151  {
1152  hitRegion[localIndices[i]] = localRegion[i];
1153  }
1154  }
1155 }
1156 
1157 
1160  const labelList& surfacesToTest,
1161  const pointField& samples,
1162  const scalarField& nearestDistSqr,
1163  labelList& hitSurface,
1164  List<pointIndexHit>& hitInfo,
1165  labelList& hitRegion,
1166  vectorField& hitNormal
1167 ) const
1168 {
1169  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1170 
1171  // Do the tests. Note that findNearest returns index in geometries.
1172  searchableSurfacesQueries::findNearest
1173  (
1174  allGeometry_,
1175  geometries,
1176  samples,
1177  nearestDistSqr,
1178  hitSurface,
1179  hitInfo
1180  );
1181 
1182  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1183  forAll(hitSurface, pointi)
1184  {
1185  if (hitSurface[pointi] != -1)
1186  {
1187  hitSurface[pointi] = surfacesToTest[hitSurface[pointi]];
1188  }
1189  }
1190 
1191  // Collect the region
1192  hitRegion.setSize(hitSurface.size());
1193  hitRegion = -1;
1194  hitNormal.setSize(hitSurface.size());
1195  hitNormal = Zero;
1196 
1197  forAll(surfacesToTest, i)
1198  {
1199  label surfI = surfacesToTest[i];
1200 
1201  // Collect hits for surfI
1202  const labelList localIndices(findIndices(hitSurface, surfI));
1203 
1204  List<pointIndexHit> localHits
1205  (
1207  (
1208  hitInfo,
1209  localIndices
1210  )
1211  );
1212 
1213  // Region
1214  labelList localRegion;
1215  allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1216 
1217  forAll(localIndices, i)
1218  {
1219  hitRegion[localIndices[i]] = localRegion[i];
1220  }
1221 
1222  // Normal
1223  vectorField localNormal;
1224  allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
1225 
1226  forAll(localIndices, i)
1227  {
1228  hitNormal[localIndices[i]] = localNormal[i];
1229  }
1230  }
1231 }
1232 
1233 
1236 //Foam::label Foam::refinementSurfaces::findHighestIntersection
1237 //(
1238 // const point& start,
1239 // const point& end,
1240 // const label currentLevel, // current cell refinement level
1241 //
1242 // pointIndexHit& maxHit
1243 //) const
1244 //{
1245 // // surface with highest maxlevel
1246 // label maxSurface = -1;
1247 // // maxLevel of maxSurface
1248 // label maxLevel = currentLevel;
1249 //
1250 // forAll(*this, surfI)
1251 // {
1252 // pointIndexHit hit = operator[](surfI).findLineAny(start, end);
1253 //
1254 // if (hit.hit())
1255 // {
1256 // const triSurface& s = operator[](surfI);
1257 //
1258 // label region = globalRegion(surfI, s[hit.index()].region());
1259 //
1260 // if (maxLevel_[region] > maxLevel)
1261 // {
1262 // maxSurface = surfI;
1263 // maxLevel = maxLevel_[region];
1264 // maxHit = hit;
1265 // }
1266 // }
1267 // }
1268 //
1269 // if (maxSurface == -1)
1270 // {
1271 // // maxLevel unchanged. No interesting surface hit.
1272 // maxHit.setMiss();
1273 // }
1274 //
1275 // return maxSurface;
1276 //}
1277 
1278 
1281  const labelList& testSurfaces,
1282  const pointField& pt,
1283  labelList& insideSurfaces
1284 ) const
1285 {
1286  insideSurfaces.setSize(pt.size());
1287  insideSurfaces = -1;
1288 
1289  forAll(testSurfaces, i)
1290  {
1291  label surfI = testSurfaces[i];
1292 
1293  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
1294 
1295  const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
1296  surfZones_[surfI].zoneInside();
1297 
1298  if
1299  (
1300  selectionMethod != surfaceZonesInfo::INSIDE
1301  && selectionMethod != surfaceZonesInfo::OUTSIDE
1302  )
1303  {
1305  << "Trying to use surface "
1306  << surface.name()
1307  << " which has non-geometric inside selection method "
1308  << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
1309  << exit(FatalError);
1310  }
1311 
1312  if (surface.hasVolumeType())
1313  {
1314  List<volumeType> volType;
1315  surface.getVolumeType(pt, volType);
1316 
1317  forAll(volType, pointi)
1318  {
1319  if (insideSurfaces[pointi] == -1)
1320  {
1321  if
1322  (
1323  (
1324  volType[pointi] == volumeType::inside
1325  && selectionMethod == surfaceZonesInfo::INSIDE
1326  )
1327  || (
1328  volType[pointi] == volumeType::outside
1329  && selectionMethod == surfaceZonesInfo::OUTSIDE
1330  )
1331  )
1332  {
1333  insideSurfaces[pointi] = surfI;
1334  }
1335  }
1336  }
1337  }
1338  }
1339 }
1340 
1341 
1342 // ************************************************************************* //
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:643
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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:59
const word & name() const
Return name.
Definition: IOobject.H:303
const entry * lookupEntryPtr(const word &, bool recursive, bool patternMatch) const
Find and return an entry data stream pointer if present.
Definition: dictionary.C:682
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:323
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:1095
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:982
T clone(const T &t)
Definition: List.H:54
void setMinLevelFields(const refinementRegions &shells)
Calculate minLevelFields.
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.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
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:335
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:844
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