FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/Time.H>
29 #include <autoMesh/shellSurfaces.H>
31 #include <OpenFOAM/labelPair.H>
33 #include <OpenFOAM/UPtrList.H>
34 
35 
36 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
37 
38 // Construct from components
39 Foam::refinementSurfaces::refinementSurfaces
40 (
41  const searchableSurfaces& allGeometry,
42  const PtrList<dictionary>& surfaceDicts
43 )
44 :
45  allGeometry_(allGeometry),
46  surfaces_(surfaceDicts.size()),
47  names_(surfaceDicts.size()),
48  faceZoneNames_(surfaceDicts.size()),
49  cellZoneNames_(surfaceDicts.size()),
50  zoneInside_(surfaceDicts.size()),
51  regionOffset_(surfaceDicts.size())
52 {
53  labelList globalMinLevel(surfaceDicts.size(), 0);
54  labelList globalMaxLevel(surfaceDicts.size(), 0);
55  scalarField globalAngle(surfaceDicts.size(), -GREAT);
56  List<Map<label> > regionMinLevel(surfaceDicts.size());
57  List<Map<label> > regionMaxLevel(surfaceDicts.size());
58  List<Map<scalar> > regionAngle(surfaceDicts.size());
59 
60  //wordList globalPatchType(surfaceDicts.size());
61  //List<HashTable<word> > regionPatchType(surfaceDicts.size());
62  //List<HashTable<word> > regionPatchName(surfaceDicts.size());
63 
64  forAll(surfaceDicts, surfI)
65  {
66  const dictionary& dict = surfaceDicts[surfI];
67 
68  dict.lookup("name") >> names_[surfI];
69 
70  surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
71 
72  // Global refinement level
73  globalMinLevel[surfI] = readLabel(dict.lookup("minRefinementLevel"));
74  globalMaxLevel[surfI] = readLabel(dict.lookup("maxRefinementLevel"));
75 
76  // Global zone names per surface
77  if (dict.found("faceZone"))
78  {
79  dict.lookup("faceZone") >> faceZoneNames_[surfI];
80  dict.lookup("cellZone") >> cellZoneNames_[surfI];
81  dict.lookup("zoneInside") >> zoneInside_[surfI];
82  }
83 
84  // Global perpendicular angle
85  if (dict.found("perpendicularAngle"))
86  {
87  globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
88  }
89 
91  //if (dict.found("patchType"))
92  //{
93  // dict.lookup("patchType") >> globalPatchType[surfI];
94  //}
95 
96 
97  if (dict.found("regions"))
98  {
99  PtrList<dictionary> regionDicts(dict.lookup("regions"));
100 
101  const wordList& regionNames =
102  allGeometry_[surfaces_[surfI]].regions();
103 
104  forAll(regionDicts, dictI)
105  {
106  const dictionary& regionDict = regionDicts[dictI];
107 
108  const word regionName(regionDict.lookup("name"));
109 
110  label regionI = findIndex(regionNames, regionName);
111 
112  if (regionI == -1)
113  {
115  (
116  "refinementSurfaces::refinementSurfaces"
117  "(const IOobject&, const PtrList<dictionary>&)"
118  ) << "No region called " << regionName << " on surface "
119  << allGeometry_[surfaces_[surfI]].name() << endl
120  << "Valid regions are " << regionNames
121  << exit(FatalError);
122  }
123 
124 
125  label min = readLabel(regionDict.lookup("minRefinementLevel"));
126  label max = readLabel(regionDict.lookup("maxRefinementLevel"));
127 
128  bool hasInserted = regionMinLevel[surfI].insert(regionI, min);
129  if (!hasInserted)
130  {
132  (
133  "refinementSurfaces::refinementSurfaces"
134  "(const IOobject&, const PtrList<dictionary>&)"
135  ) << "Duplicate region name " << regionName
136  << " on surface " << names_[surfI]
137  << exit(FatalError);
138  }
139  regionMaxLevel[surfI].insert(regionI, max);
140 
141  if (regionDict.found("perpendicularAngle"))
142  {
143  regionAngle[surfI].insert
144  (
145  regionI,
146  readScalar(regionDict.lookup("perpendicularAngle"))
147  );
148  }
149  }
150  }
151  }
152 
153 
154  // Check for duplicate surface names
155  {
156  HashTable<label> surfaceNames(names_.size());
157 
158  forAll(names_, surfI)
159  {
160  if (!surfaceNames.insert(names_[surfI], surfI))
161  {
163  (
164  "refinementSurfaces::refinementSurfaces"
165  "(const IOobject&, const PtrList<dictionary>&)"
166  ) << "Duplicate surface name " << names_[surfI] << endl
167  << "Previous occurrence of name at surface "
168  << surfaceNames[names_[surfI]]
169  << exit(FatalError);
170  }
171  }
172  }
173 
174  // Calculate local to global region offset
175  label nRegions = 0;
176 
177  forAll(surfaceDicts, surfI)
178  {
179  regionOffset_[surfI] = nRegions;
180 
181  nRegions += allGeometry_[surfaces_[surfI]].regions().size();
182  }
183 
184  // Rework surface specific information into information per global region
185  minLevel_.setSize(nRegions);
186  minLevel_ = 0;
187  maxLevel_.setSize(nRegions);
188  maxLevel_ = 0;
189  perpendicularAngle_.setSize(nRegions);
190  perpendicularAngle_ = -GREAT;
191  //patchName_.setSize(nRegions);
192  //patchType_.setSize(nRegions);
193 
194  forAll(surfaceDicts, surfI)
195  {
196  label nRegions = allGeometry_[surfaces_[surfI]].regions().size();
197 
198  // Initialise to global (i.e. per surface)
199  for (label i = 0; i < nRegions; i++)
200  {
201  minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
202  maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
203  perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
204  }
205 
206  // Overwrite with region specific information
207  forAllConstIter(Map<label>, regionMinLevel[surfI], iter)
208  {
209  label globalRegionI = regionOffset_[surfI] + iter.key();
210 
211  minLevel_[globalRegionI] = iter();
212  maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
213 
214  // Check validity
215  if
216  (
217  minLevel_[globalRegionI] < 0
218  || maxLevel_[globalRegionI] < minLevel_[globalRegionI]
219  )
220  {
222  (
223  "refinementSurfaces::refinementSurfaces"
224  "(const IOobject&, const PtrList<dictionary>&)"
225  ) << "Illegal level or layer specification for surface "
226  << names_[surfI]
227  << " : minLevel:" << minLevel_[globalRegionI]
228  << " maxLevel:" << maxLevel_[globalRegionI]
229  << exit(FatalError);
230  }
231  }
232  forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
233  {
234  label globalRegionI = regionOffset_[surfI] + iter.key();
235 
236  perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
237  }
238 
239 
241  //forAllConstIter(HashTable<word>, regionPatchName[surfI], iter)
242  //{
243  // label regionI = findIndex(regionNames, iter.key());
244  // label globalRegionI = regionOffset_[surfI] + regionI;
245  //
246  // patchName_[globalRegionI] = iter();
247  // patchType_[globalRegionI] = regionPatchType[surfI][iter.key()];
248  //}
249  }
250 }
251 
252 
253 Foam::refinementSurfaces::refinementSurfaces
254 (
255  const searchableSurfaces& allGeometry,
256  const dictionary& surfacesDict
257 )
258 :
259  allGeometry_(allGeometry),
260  surfaces_(surfacesDict.size()),
261  names_(surfacesDict.size()),
262  faceZoneNames_(surfacesDict.size()),
263  cellZoneNames_(surfacesDict.size()),
264  zoneInside_(surfacesDict.size()),
265  regionOffset_(surfacesDict.size())
266 {
267  // Wilcard specification : loop over all surface, all regions
268  // and try to find a match.
269 
270  // Count number of surfaces.
271  label surfI = 0;
272  forAll(allGeometry.names(), geomI)
273  {
274  const word& geomName = allGeometry_.names()[geomI];
275 
276  if (surfacesDict.found(geomName))
277  {
278  surfI++;
279  }
280  }
281 
282  // Size lists
283  surfaces_.setSize(surfI);
284  names_.setSize(surfI);
285  faceZoneNames_.setSize(surfI);
286  cellZoneNames_.setSize(surfI);
287  zoneInside_.setSize(surfI);
288  regionOffset_.setSize(surfI);
289 
290  labelList globalMinLevel(surfI, 0);
291  labelList globalMaxLevel(surfI, 0);
292  scalarField globalAngle(surfI, -GREAT);
293  List<Map<label> > regionMinLevel(surfI);
294  List<Map<label> > regionMaxLevel(surfI);
295  List<Map<scalar> > regionAngle(surfI);
296 
297  surfI = 0;
298  forAll(allGeometry.names(), geomI)
299  {
300  const word& geomName = allGeometry_.names()[geomI];
301 
302  if (surfacesDict.found(geomName))
303  {
304  const dictionary& dict = surfacesDict.subDict(geomName);
305 
306  names_[surfI] = geomName;
307  surfaces_[surfI] = geomI;
308 
309  const labelPair refLevel(dict.lookup("level"));
310  globalMinLevel[surfI] = refLevel[0];
311  globalMaxLevel[surfI] = refLevel[1];
312 
313  // Global zone names per surface
314  if (dict.found("faceZone"))
315  {
316  dict.lookup("faceZone") >> faceZoneNames_[surfI];
317  dict.lookup("cellZone") >> cellZoneNames_[surfI];
318  dict.lookup("zoneInside") >> zoneInside_[surfI];
319  }
320 
321  // Global perpendicular angle
322  if (dict.found("perpendicularAngle"))
323  {
324  globalAngle[surfI] = readScalar
325  (
326  dict.lookup("perpendicularAngle")
327  );
328  }
329 
330  if (dict.found("regions"))
331  {
332  const dictionary& regionsDict = dict.subDict("regions");
333  const wordList& regionNames =
334  allGeometry_[surfaces_[surfI]].regions();
335 
336  forAll(regionNames, regionI)
337  {
338  if (regionsDict.found(regionNames[regionI]))
339  {
340  // Get the dictionary for region
341  const dictionary& regionDict = regionsDict.subDict
342  (
343  regionNames[regionI]
344  );
345 
346  const labelPair refLevel(regionDict.lookup("level"));
347 
348  regionMinLevel[surfI].insert(regionI, refLevel[0]);
349  regionMaxLevel[surfI].insert(regionI, refLevel[1]);
350 
351  if (regionDict.found("perpendicularAngle"))
352  {
353  regionAngle[surfI].insert
354  (
355  regionI,
356  readScalar
357  (
358  regionDict.lookup("perpendicularAngle")
359  )
360  );
361  }
362  }
363  }
364  }
365  surfI++;
366  }
367  }
368 
369  // Calculate local to global region offset
370  label nRegions = 0;
371 
372  forAll(surfaces_, surfI)
373  {
374  regionOffset_[surfI] = nRegions;
375  nRegions += allGeometry_[surfaces_[surfI]].regions().size();
376  }
377 
378  // Rework surface specific information into information per global region
379  minLevel_.setSize(nRegions);
380  minLevel_ = 0;
381  maxLevel_.setSize(nRegions);
382  maxLevel_ = 0;
383  perpendicularAngle_.setSize(nRegions);
384  perpendicularAngle_ = -GREAT;
385 
386 
387  forAll(globalMinLevel, surfI)
388  {
389  label nRegions = allGeometry_[surfaces_[surfI]].regions().size();
390 
391  // Initialise to global (i.e. per surface)
392  for (label i = 0; i < nRegions; i++)
393  {
394  minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
395  maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
396  perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
397  }
398 
399  // Overwrite with region specific information
400  forAllConstIter(Map<label>, regionMinLevel[surfI], iter)
401  {
402  label globalRegionI = regionOffset_[surfI] + iter.key();
403 
404  minLevel_[globalRegionI] = iter();
405  maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
406 
407  // Check validity
408  if
409  (
410  minLevel_[globalRegionI] < 0
411  || maxLevel_[globalRegionI] < minLevel_[globalRegionI]
412  )
413  {
415  (
416  "refinementSurfaces::refinementSurfaces"
417  "(const searchableSurfaces&, const dictionary>&"
418  ) << "Illegal level or layer specification for surface "
419  << names_[surfI]
420  << " : minLevel:" << minLevel_[globalRegionI]
421  << " maxLevel:" << maxLevel_[globalRegionI]
422  << exit(FatalError);
423  }
424  }
425  forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
426  {
427  label globalRegionI = regionOffset_[surfI] + iter.key();
428 
429  perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
430  }
431  }
432 }
433 
434 
435 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
436 
437 // Get indices of unnamed surfaces (surfaces without faceZoneName)
439 {
440  labelList anonymousSurfaces(faceZoneNames_.size());
441 
442  label i = 0;
443  forAll(faceZoneNames_, surfI)
444  {
445  if (faceZoneNames_[surfI].empty())
446  {
447  anonymousSurfaces[i++] = surfI;
448  }
449  }
450  anonymousSurfaces.setSize(i);
451 
452  return anonymousSurfaces;
453 }
454 
455 
456 // Get indices of named surfaces (surfaces with faceZoneName)
458 {
459  labelList namedSurfaces(faceZoneNames_.size());
460 
461  label namedI = 0;
462  forAll(faceZoneNames_, surfI)
463  {
464  if (faceZoneNames_[surfI].size())
465  {
466  namedSurfaces[namedI++] = surfI;
467  }
468  }
469  namedSurfaces.setSize(namedI);
470 
471  return namedSurfaces;
472 }
473 
474 
475 // Get indices of closed named surfaces
477 {
478  labelList named(getNamedSurfaces());
479 
480  labelList closed(named.size());
481  label closedI = 0;
482 
483  forAll(named, i)
484  {
485  label surfI = named[i];
486 
487  if (allGeometry_[surfaces_[surfI]].hasVolumeType())
488  {
489  closed[closedI++] = surfI;
490  }
491  }
492  closed.setSize(closedI);
493 
494  return closed;
495 }
496 
497 
498 // // Count number of triangles per surface region
499 // Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
500 // {
501 // const geometricSurfacePatchList& regions = s.patches();
502 //
503 // labelList nTris(regions.size(), 0);
504 //
505 // forAll(s, triI)
506 // {
507 // nTris[s[triI].region()]++;
508 // }
509 // return nTris;
510 // }
511 
512 
513 // // Pre-calculate the refinement level for every element
514 // void Foam::refinementSurfaces::wantedRefinementLevel
515 // (
516 // const shellSurfaces& shells,
517 // const label surfI,
518 // const List<pointIndexHit>& info, // Indices
519 // const pointField& ctrs, // Representative coordinate
520 // labelList& minLevelField
521 // ) const
522 // {
523 // const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
524 //
525 // // Get per element the region
526 // labelList region;
527 // geom.getRegion(info, region);
528 //
529 // // Initialise fields to region wise minLevel
530 // minLevelField.setSize(ctrs.size());
531 // minLevelField = -1;
532 //
533 // forAll(minLevelField, i)
534 // {
535 // if (info[i].hit())
536 // {
537 // minLevelField[i] = minLevel(surfI, region[i]);
538 // }
539 // }
540 //
541 // // Find out if triangle inside shell with higher level
542 // // What level does shell want to refine fc to?
543 // labelList shellLevel;
544 // shells.findHigherLevel(ctrs, minLevelField, shellLevel);
545 //
546 // forAll(minLevelField, i)
547 // {
548 // minLevelField[i] = max(minLevelField[i], shellLevel[i]);
549 // }
550 // }
551 
552 
553 // Precalculate the refinement level for every element of the searchable
554 // surface.
556 (
557  const shellSurfaces& shells
558 )
559 {
560  forAll(surfaces_, surfI)
561  {
562  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
563 
564  // Precalculation only makes sense if there are different regions
565  // (so different refinement levels possible) and there are some
566  // elements. Possibly should have 'enough' elements to have fine
567  // enough resolution but for now just make sure we don't catch e.g.
568  // searchableBox (size=6)
569  if (geom.regions().size() > 1 && geom.globalSize() > 10)
570  {
571  // Representative local coordinates
572  const pointField ctrs = geom.coordinates();
573 
574  labelList minLevelField(ctrs.size(), -1);
575  {
576  // Get the element index in a roundabout way. Problem is e.g.
577  // distributed surface where local indices differ from global
578  // ones (needed for getRegion call)
579  List<pointIndexHit> info;
580  geom.findNearest
581  (
582  ctrs,
583  scalarField(ctrs.size(), sqr(GREAT)),
584  info
585  );
586 
587  // Get per element the region
588  labelList region;
589  geom.getRegion(info, region);
590 
591  // From the region get the surface-wise refinement level
592  forAll(minLevelField, i)
593  {
594  if (info[i].hit())
595  {
596  minLevelField[i] = minLevel(surfI, region[i]);
597  }
598  }
599  }
600 
601  // Find out if triangle inside shell with higher level
602  // What level does shell want to refine fc to?
603  labelList shellLevel;
604  shells.findHigherLevel(ctrs, minLevelField, shellLevel);
605 
606  forAll(minLevelField, i)
607  {
608  minLevelField[i] = max(minLevelField[i], shellLevel[i]);
609  }
610 
611  // Store minLevelField on surface
612  const_cast<searchableSurface&>(geom).setField(minLevelField);
613  }
614  }
615 }
616 
617 
618 // Find intersections of edge. Return -1 or first surface with higher minLevel
619 // number.
621 (
622  const pointField& start,
623  const pointField& end,
624  const labelList& currentLevel, // current cell refinement level
625 
626  labelList& surfaces,
627  labelList& surfaceLevel
628 ) const
629 {
630  surfaces.setSize(start.size());
631  surfaces = -1;
632  surfaceLevel.setSize(start.size());
633  surfaceLevel = -1;
634 
635  if (surfaces_.empty())
636  {
637  return;
638  }
639 
640  if (surfaces_.size() == 1)
641  {
642  // Optimisation: single segmented surface. No need to duplicate
643  // point storage.
644 
645  label surfI = 0;
646 
647  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
648 
649  // Do intersection test
650  List<pointIndexHit> intersectionInfo(start.size());
651  geom.findLineAny(start, end, intersectionInfo);
652 
653  // See if a cached level field available
654  labelList minLevelField;
655  geom.getField(intersectionInfo, minLevelField);
656  bool haveLevelField =
657  (
658  returnReduce(minLevelField.size(), sumOp<label>())
659  > 0
660  );
661 
662  if (!haveLevelField && geom.regions().size() == 1)
663  {
664  minLevelField = labelList
665  (
666  intersectionInfo.size(),
667  minLevel(surfI, 0)
668  );
669  haveLevelField = true;
670  }
671 
672  if (haveLevelField)
673  {
674  forAll(intersectionInfo, i)
675  {
676  if
677  (
678  intersectionInfo[i].hit()
679  && minLevelField[i] > currentLevel[i]
680  )
681  {
682  surfaces[i] = surfI; // index of surface
683  surfaceLevel[i] = minLevelField[i];
684  }
685  }
686  return;
687  }
688  }
689 
690 
691 
692  // Work arrays
693  pointField p0(start);
694  pointField p1(end);
695  labelList intersectionToPoint(identity(start.size()));
696  List<pointIndexHit> intersectionInfo(start.size());
697 
698  forAll(surfaces_, surfI)
699  {
700  const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
701 
702  // Do intersection test
703  geom.findLineAny(p0, p1, intersectionInfo);
704 
705  // See if a cached level field available
706  labelList minLevelField;
707  geom.getField(intersectionInfo, minLevelField);
708 
709  // Copy all hits into arguments, In-place compact misses.
710  label missI = 0;
711  forAll(intersectionInfo, i)
712  {
713  // Get the minLevel for the point
714  label minLocalLevel = -1;
715 
716  if (intersectionInfo[i].hit())
717  {
718  // Check if minLevelField for this surface.
719  if (minLevelField.size())
720  {
721  minLocalLevel = minLevelField[i];
722  }
723  else
724  {
725  // Use the min level for the surface instead. Assume
726  // single region 0.
727  minLocalLevel = minLevel(surfI, 0);
728  }
729  }
730 
731 
732  label pointI = intersectionToPoint[i];
733 
734  if (minLocalLevel > currentLevel[pointI])
735  {
736  // Mark point for refinement
737  surfaces[pointI] = surfI;
738  surfaceLevel[pointI] = minLocalLevel;
739  }
740  else
741  {
742  p0[missI] = start[pointI];
743  p1[missI] = end[pointI];
744  intersectionToPoint[missI] = pointI;
745  missI++;
746  }
747  }
748 
749  // All done? Note that this decision should be synchronised
750  if (returnReduce(missI, sumOp<label>()) == 0)
751  {
752  break;
753  }
754 
755  // Trim misses
756  p0.setSize(missI);
757  p1.setSize(missI);
758  intersectionToPoint.setSize(missI);
759  intersectionInfo.setSize(missI);
760  }
761 }
762 
763 
765 (
766  const pointField& start,
767  const pointField& end,
768  const labelList& currentLevel, // current cell refinement level
769 
770  List<vectorList>& surfaceNormal,
771  labelListList& surfaceLevel
772 ) const
773 {
774  surfaceLevel.setSize(start.size());
775  surfaceNormal.setSize(start.size());
776 
777  if (surfaces_.empty())
778  {
779  return;
780  }
781 
782  // Work arrays
783  List<List<pointIndexHit> > hitInfo;
784  labelList pRegions;
785  vectorField pNormals;
786 
787  forAll(surfaces_, surfI)
788  {
789  allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
790 
791  // Repack hits for surface into flat list
792  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793  // To avoid overhead of calling getRegion for every point
794 
795  label n = 0;
796  forAll(hitInfo, pointI)
797  {
798  n += hitInfo[pointI].size();
799  }
800 
801  List<pointIndexHit> surfInfo(n);
802  labelList pointMap(n);
803  n = 0;
804 
805  forAll(hitInfo, pointI)
806  {
807  const List<pointIndexHit>& pHits = hitInfo[pointI];
808 
809  forAll(pHits, i)
810  {
811  surfInfo[n] = pHits[i];
812  pointMap[n] = pointI;
813  n++;
814  }
815  }
816 
817  labelList surfRegion(n);
818  vectorField surfNormal(n);
819  allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
820  allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
821 
822  surfInfo.clear();
823 
824 
825  // Extract back into pointwise
826  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
827 
828  forAll(surfRegion, i)
829  {
830  label region = globalRegion(surfI, surfRegion[i]);
831  label pointI = pointMap[i];
832 
833  if (maxLevel_[region] > currentLevel[pointI])
834  {
835  // Append to pointI info
836  label sz = surfaceNormal[pointI].size();
837  surfaceNormal[pointI].setSize(sz+1);
838  surfaceNormal[pointI][sz] = surfNormal[i];
839 
840  surfaceLevel[pointI].setSize(sz+1);
841  surfaceLevel[pointI][sz] = maxLevel_[region];
842  }
843  }
844  }
845 }
846 
847 
849 (
850  const labelList& surfacesToTest,
851  const pointField& start,
852  const pointField& end,
853 
854  labelList& surface1,
855  List<pointIndexHit>& hit1,
856  labelList& region1,
857  labelList& surface2,
858  List<pointIndexHit>& hit2,
859  labelList& region2
860 ) const
861 {
862  // 1. intersection from start to end
863  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
864 
865  // Initialize arguments
866  surface1.setSize(start.size());
867  surface1 = -1;
868  hit1.setSize(start.size());
869  region1.setSize(start.size());
870 
871  // Current end of segment to test.
872  pointField nearest(end);
873  // Work array
874  List<pointIndexHit> nearestInfo(start.size());
875  labelList region;
876 
877  forAll(surfacesToTest, testI)
878  {
879  label surfI = surfacesToTest[testI];
880 
881  // See if any intersection between start and current nearest
882  allGeometry_[surfaces_[surfI]].findLine
883  (
884  start,
885  nearest,
886  nearestInfo
887  );
888  allGeometry_[surfaces_[surfI]].getRegion
889  (
890  nearestInfo,
891  region
892  );
893 
894  forAll(nearestInfo, pointI)
895  {
896  if (nearestInfo[pointI].hit())
897  {
898  hit1[pointI] = nearestInfo[pointI];
899  surface1[pointI] = surfI;
900  region1[pointI] = region[pointI];
901  nearest[pointI] = hit1[pointI].hitPoint();
902  }
903  }
904  }
905 
906 
907  // 2. intersection from end to last intersection
908  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
909 
910  // Find the nearest intersection from end to start. Note that we initialize
911  // to the first intersection (if any).
912  surface2 = surface1;
913  hit2 = hit1;
914  region2 = region1;
915 
916  // Set current end of segment to test.
917  forAll(nearest, pointI)
918  {
919  if (hit1[pointI].hit())
920  {
921  nearest[pointI] = hit1[pointI].hitPoint();
922  }
923  else
924  {
925  // Disable testing by setting to end.
926  nearest[pointI] = end[pointI];
927  }
928  }
929 
930  forAll(surfacesToTest, testI)
931  {
932  label surfI = surfacesToTest[testI];
933 
934  // See if any intersection between end and current nearest
935  allGeometry_[surfaces_[surfI]].findLine
936  (
937  end,
938  nearest,
939  nearestInfo
940  );
941  allGeometry_[surfaces_[surfI]].getRegion
942  (
943  nearestInfo,
944  region
945  );
946 
947  forAll(nearestInfo, pointI)
948  {
949  if (nearestInfo[pointI].hit())
950  {
951  hit2[pointI] = nearestInfo[pointI];
952  surface2[pointI] = surfI;
953  region2[pointI] = region[pointI];
954  nearest[pointI] = hit2[pointI].hitPoint();
955  }
956  }
957  }
958 
959 
960  // Make sure that if hit1 has hit something, hit2 will have at least the
961  // same point (due to tolerances it might miss its end point)
962  forAll(hit1, pointI)
963  {
964  if (hit1[pointI].hit() && !hit2[pointI].hit())
965  {
966  hit2[pointI] = hit1[pointI];
967  surface2[pointI] = surface1[pointI];
968  region2[pointI] = region1[pointI];
969  }
970  }
971 }
972 
973 
975 (
976  const pointField& start,
977  const pointField& end,
978 
979  labelList& hitSurface,
980  List<pointIndexHit>& hitInfo
981 ) const
982 {
984  (
985  allGeometry_,
986  surfaces_,
987  start,
988  end,
989  hitSurface,
990  hitInfo
991  );
992 }
993 
994 
996 (
997  const labelList& surfacesToTest,
998  const pointField& samples,
999  const scalarField& nearestDistSqr,
1000  labelList& hitSurface,
1001  List<pointIndexHit>& hitInfo
1002 ) const
1003 {
1004  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1005 
1006  // Do the tests. Note that findNearest returns index in geometries.
1008  (
1009  allGeometry_,
1010  geometries,
1011  samples,
1012  nearestDistSqr,
1013  hitSurface,
1014  hitInfo
1015  );
1016 
1017  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1018  forAll(hitSurface, pointI)
1019  {
1020  if (hitSurface[pointI] != -1)
1021  {
1022  hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1023  }
1024  }
1025 }
1026 
1027 
1030  const labelList& surfacesToTest,
1031  const pointField& samples,
1032  const scalarField& nearestDistSqr,
1033  labelList& hitSurface,
1034  labelList& hitRegion
1035 ) const
1036 {
1037  labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
1038 
1039  // Do the tests. Note that findNearest returns index in geometries.
1040  List<pointIndexHit> hitInfo;
1042  (
1043  allGeometry_,
1044  geometries,
1045  samples,
1046  nearestDistSqr,
1047  hitSurface,
1048  hitInfo
1049  );
1050 
1051  // Rework the hitSurface to be surface (i.e. index into surfaces_)
1052  forAll(hitSurface, pointI)
1053  {
1054  if (hitSurface[pointI] != -1)
1055  {
1056  hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
1057  }
1058  }
1059 
1060  // Collect the region
1061  hitRegion.setSize(hitSurface.size());
1062  hitRegion = -1;
1063 
1064  forAll(surfacesToTest, i)
1065  {
1066  label surfI = surfacesToTest[i];
1067 
1068  // Collect hits for surfI
1069  const labelList localIndices(findIndices(hitSurface, surfI));
1070 
1071  List<pointIndexHit> localHits
1072  (
1074  (
1075  hitInfo,
1076  localIndices
1077  )
1078  );
1079 
1080  labelList localRegion;
1081  allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
1082 
1083  forAll(localIndices, i)
1084  {
1085  hitRegion[localIndices[i]] = localRegion[i];
1086  }
1087  }
1088 }
1089 
1090 
1093 //Foam::label Foam::refinementSurfaces::findHighestIntersection
1094 //(
1095 // const point& start,
1096 // const point& end,
1097 // const label currentLevel, // current cell refinement level
1098 //
1099 // pointIndexHit& maxHit
1100 //) const
1101 //{
1102 // // surface with highest maxlevel
1103 // label maxSurface = -1;
1104 // // maxLevel of maxSurface
1105 // label maxLevel = currentLevel;
1106 //
1107 // forAll(*this, surfI)
1108 // {
1109 // pointIndexHit hit = operator[](surfI).findLineAny(start, end);
1110 //
1111 // if (hit.hit())
1112 // {
1113 // const triSurface& s = operator[](surfI);
1114 //
1115 // label region = globalRegion(surfI, s[hit.index()].region());
1116 //
1117 // if (maxLevel_[region] > maxLevel)
1118 // {
1119 // maxSurface = surfI;
1120 // maxLevel = maxLevel_[region];
1121 // maxHit = hit;
1122 // }
1123 // }
1124 // }
1125 //
1126 // if (maxSurface == -1)
1127 // {
1128 // // maxLevel unchanged. No interesting surface hit.
1129 // maxHit.setMiss();
1130 // }
1131 //
1132 // return maxSurface;
1133 //}
1134 
1135 
1138  const labelList& testSurfaces,
1139  const pointField& pt,
1140  labelList& insideSurfaces
1141 ) const
1142 {
1143  insideSurfaces.setSize(pt.size());
1144  insideSurfaces = -1;
1145 
1146  forAll(testSurfaces, i)
1147  {
1148  label surfI = testSurfaces[i];
1149 
1150  if (allGeometry_[surfaces_[surfI]].hasVolumeType())
1151  {
1153  allGeometry_[surfaces_[surfI]].getVolumeType(pt, volType);
1154 
1155  forAll(volType, pointI)
1156  {
1157  if (insideSurfaces[pointI] == -1)
1158  {
1159  if
1160  (
1161  (
1162  volType[pointI] == triSurfaceMesh::INSIDE
1163  && zoneInside_[surfI]
1164  )
1165  || (
1166  volType[pointI] == triSurfaceMesh::OUTSIDE
1167  && !zoneInside_[surfI]
1168  )
1169  )
1170  {
1171  insideSurfaces[pointI] = surfI;
1172  }
1173  }
1174  }
1175  }
1176  }
1177 }
1178 
1179 
1180 // ************************ vim: set sw=4 sts=4 et: ************************ //