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

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17