conformationSurfaces.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) 2012-2018 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 "conformationSurfaces.H"
27 #include "conformalVoronoiMesh.H"
28 #include "triSurface.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 defineTypeNameAndDebug(conformationSurfaces, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
40 
41 void Foam::conformationSurfaces::hasBoundedVolume
42 (
43  List<volumeType>& referenceVolumeTypes
44 ) const
45 {
46  vector sum(Zero);
47  label totalTriangles = 0;
48 
49  forAll(surfaces_, s)
50  {
51  const searchableSurface& surface(allGeometry_[surfaces_[s]]);
52 
53  if
54  (
55  surface.hasVolumeType()
56  && (
57  normalVolumeTypes_[regionOffset_[s]]
59  )
60  )
61  {
62  pointField pts(1, locationInMesh_);
63 
64  List<volumeType> vTypes
65  (
66  pts.size(),
68  );
69 
70  surface.getVolumeType(pts, vTypes);
71 
72  referenceVolumeTypes[s] = vTypes[0];
73 
74  Info<< " is "
75  << volumeType::names[referenceVolumeTypes[s]]
76  << " surface " << surface.name()
77  << endl;
78  }
79 
80  if (isA<triSurface>(surface))
81  {
82  const triSurface& triSurf = refCast<const triSurface>(surface);
83 
84  const pointField& surfPts = triSurf.points();
85 
86  Info<< " Checking " << surface.name() << endl;
87 
88  label nBaffles = 0;
89 
90  Info<< " Index = " << surfaces_[s] << endl;
91  Info<< " Offset = " << regionOffset_[s] << endl;
92 
93  forAll(triSurf, sI)
94  {
95  const label patchID =
96  triSurf[sI].region()
97  + regionOffset_[s];
98 
99  // Don't include baffle surfaces in the calculation
100  if
101  (
102  normalVolumeTypes_[patchID]
104  )
105  {
106  sum += triSurf[sI].area(surfPts);
107  }
108  else
109  {
110  nBaffles++;
111  }
112  }
113  Info<< " has " << nBaffles << " baffles out of "
114  << triSurf.size() << " triangles" << endl;
115 
116  totalTriangles += triSurf.size();
117  }
118  }
119 
120  Info<< " Sum of all the surface normals (if near zero, surface is"
121  << " probably closed):" << nl
122  << " Note: Does not include baffle surfaces in calculation" << nl
123  << " Sum = " << sum/(totalTriangles + small) << nl
124  << " mag(Sum) = " << mag(sum)/(totalTriangles + small)
125  << endl;
126 }
127 
128 
129 void Foam::conformationSurfaces::readFeatures
130 (
131  const label surfI,
132  const dictionary& featureDict,
133  const word& surfaceName,
134  label& featureIndex
135 )
136 {
137  word featureMethod =
138  featureDict.lookupOrDefault<word>("featureMethod", "none");
139 
140  if (featureMethod == "extendedFeatureEdgeMesh")
141  {
142  fileName feMeshName(featureDict.lookup("extendedFeatureEdgeMesh"));
143 
144  Info<< " features: " << feMeshName << endl;
145 
146  features_.set
147  (
148  featureIndex,
149  new extendedFeatureEdgeMesh
150  (
151  IOobject
152  (
153  feMeshName,
154  runTime_.time().constant(),
155  "extendedFeatureEdgeMesh",
156  runTime_.time(),
159  )
160  )
161  );
162 
163  featureIndex++;
164  }
165  else if (featureMethod == "extractFeatures")
166  {
167  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
168 
169  Info<< " features: " << surface.name()
170  << " of type " << surface.type()
171  << ", id: " << featureIndex << endl;
172 
173  autoPtr<searchableSurfaceFeatures> ssFeatures
174  (
175  searchableSurfaceFeatures::New(surface, featureDict)
176  );
177 
178  if (ssFeatures().hasFeatures())
179  {
180  features_.set
181  (
182  featureIndex,
183  ssFeatures().features()
184  );
185 
186  featureIndex++;
187  }
188  else
189  {
191  << surface.name() << " of type "
192  << surface.type() << " does not have features"
193  << endl;
194  }
195  }
196  else if (featureMethod == "none")
197  {
198  // Currently nothing to do
199  }
200  else
201  {
203  << "No valid featureMethod found for surface " << surfaceName
204  << nl << "Use \"extendedFeatureEdgeMesh\" "
205  << "or \"extractFeatures\"."
206  << exit(FatalError);
207  }
208 }
209 
210 void Foam::conformationSurfaces::readFeatures
211 (
212  const dictionary& featureDict,
213  const word& surfaceName,
214  label& featureIndex
215 )
216 {
217  word featureMethod =
218  featureDict.lookupOrDefault<word>("featureMethod", "none");
219 
220  if (featureMethod == "extendedFeatureEdgeMesh")
221  {
222  fileName feMeshName(featureDict.lookup("extendedFeatureEdgeMesh"));
223 
224  Info<< " features: " << feMeshName << ", id: " << featureIndex
225  << endl;
226 
227  features_.set
228  (
229  featureIndex,
230  new extendedFeatureEdgeMesh
231  (
232  IOobject
233  (
234  feMeshName,
235  runTime_.time().constant(),
236  "extendedFeatureEdgeMesh",
237  runTime_.time(),
240  )
241  )
242  );
243 
244  featureIndex++;
245  }
246  else if (featureMethod == "none")
247  {
248  // Currently nothing to do
249  }
250  else
251  {
253  << "No valid featureMethod found for surface " << surfaceName
254  << nl << "Use \"extendedFeatureEdgeMesh\" "
255  << "or \"extractFeatures\"."
256  << exit(FatalError);
257  }
258 }
259 
260 
261 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
262 
264 (
265  const Time& runTime,
266  Random& rndGen,
267  const searchableSurfaces& allGeometry,
268  const dictionary& surfaceConformationDict
269 )
270 :
271  runTime_(runTime),
272  rndGen_(rndGen),
273  allGeometry_(allGeometry),
274  features_(),
275  locationInMesh_(surfaceConformationDict.lookup("locationInMesh")),
276  surfaces_(),
277  allGeometryToSurfaces_(),
278  normalVolumeTypes_(),
279  patchNames_(),
280  surfZones_(),
281  regionOffset_(),
282  patchInfo_(),
283  globalBounds_(),
284  referenceVolumeTypes_(0)
285 {
286  const dictionary& surfacesDict
287  (
288  surfaceConformationDict.subDict("geometryToConformTo")
289  );
290 
291  const dictionary& additionalFeaturesDict
292  (
293  surfaceConformationDict.subDict("additionalFeatures")
294  );
295 
296 
297  // Wildcard specification : loop over all surface, all regions
298  // and try to find a match.
299 
300  // Count number of surfaces.
301  label surfI = 0;
302  forAll(allGeometry.names(), geomI)
303  {
304  const word& geomName = allGeometry_.names()[geomI];
305 
306  if (surfacesDict.found(geomName))
307  {
308  surfI++;
309  }
310  }
311 
312  const label nAddFeat = additionalFeaturesDict.size();
313 
314  Info<< nl << "Reading geometryToConformTo" << endl;
315 
316  allGeometryToSurfaces_.setSize(allGeometry_.size(), -1);
317 
318  normalVolumeTypes_.setSize(surfI);
319  surfaces_.setSize(surfI);
320  surfZones_.setSize(surfI);
321 
322  // Features may be attached to host surfaces or independent
323  features_.setSize(surfI + nAddFeat);
324 
325  label featureI = 0;
326 
327  regionOffset_.setSize(surfI, 0);
328 
329  PtrList<dictionary> globalPatchInfo(surfI);
330  List<Map<autoPtr<dictionary>>> regionPatchInfo(surfI);
331  List<sideVolumeType> globalVolumeTypes(surfI);
332  List<Map<sideVolumeType>> regionVolumeTypes(surfI);
333 
334  HashSet<word> unmatchedKeys(surfacesDict.toc());
335 
336  surfI = 0;
337  forAll(allGeometry_.names(), geomI)
338  {
339  const word& geomName = allGeometry_.names()[geomI];
340 
341  const entry* ePtr = surfacesDict.lookupEntryPtr(geomName, false, true);
342 
343  if (ePtr)
344  {
345  const dictionary& dict = ePtr->dict();
346  unmatchedKeys.erase(ePtr->keyword());
347 
348  surfaces_[surfI] = geomI;
349 
350  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
351 
352  // Surface zones
353  if (dict.found("faceZone"))
354  {
355  surfZones_.set(surfI, new surfaceZonesInfo(surface, dict));
356  }
357 
358  allGeometryToSurfaces_[surfaces_[surfI]] = surfI;
359 
360  Info<< nl << " " << geomName << endl;
361 
362  const wordList& regionNames =
363  allGeometry_.regionNames()[surfaces_[surfI]];
364 
365  patchNames_.append(regionNames);
366 
367  globalVolumeTypes[surfI] =
368  (
370  [
371  dict.lookupOrDefault<word>
372  (
373  "meshableSide",
374  "inside"
375  )
376  ]
377  );
378 
379  if (!globalVolumeTypes[surfI])
380  {
381  if (!surface.hasVolumeType())
382  {
384  << "Non-baffle surface "
385  << surface.name()
386  << " does not allow inside/outside queries."
387  << " This usually is an error." << endl;
388  }
389  }
390 
391  // Load patch info
392  if (dict.found("patchInfo"))
393  {
394  globalPatchInfo.set
395  (
396  surfI,
397  dict.subDict("patchInfo").clone()
398  );
399  }
400 
401  readFeatures
402  (
403  surfI,
404  dict,
405  geomName,
406  featureI
407  );
408 
409  const wordList& rNames = surface.regions();
410 
411  if (dict.found("regions"))
412  {
413  const dictionary& regionsDict = dict.subDict("regions");
414 
415  forAll(rNames, regionI)
416  {
417  const word& regionName = rNames[regionI];
418 
419  if (regionsDict.found(regionName))
420  {
421  Info<< " region " << regionName << endl;
422 
423  // Get the dictionary for region
424  const dictionary& regionDict = regionsDict.subDict
425  (
426  regionName
427  );
428 
429  if (regionDict.found("patchInfo"))
430  {
431  regionPatchInfo[surfI].insert
432  (
433  regionI,
434  regionDict.subDict("patchInfo").clone()
435  );
436  }
437 
438  regionVolumeTypes[surfI].insert
439  (
440  regionI,
442  [
443  regionDict.lookupOrDefault<word>
444  (
445  "meshableSide",
446  extendedFeatureEdgeMesh::
447  sideVolumeTypeNames_
448  [
449  globalVolumeTypes[surfI]
450  ]
451  )
452  ]
453  );
454 
455  readFeatures(regionDict, regionName, featureI);
456  }
457  }
458  }
459 
460  surfI++;
461  }
462  }
463 
464 
465  if (unmatchedKeys.size() > 0)
466  {
468  (
469  surfacesDict
470  ) << "Not all entries in conformationSurfaces dictionary were used."
471  << " The following entries were not used : "
472  << unmatchedKeys.sortedToc()
473  << endl;
474  }
475 
476 
477  // Calculate local to global region offset
478  label nRegions = 0;
479 
480  forAll(surfaces_, surfI)
481  {
482  regionOffset_[surfI] = nRegions;
483 
484  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
485  nRegions += surface.regions().size();
486  }
487 
488  // Rework surface specific information into information per global region
489  patchInfo_.setSize(nRegions);
490  normalVolumeTypes_.setSize(nRegions);
491 
492  forAll(surfaces_, surfI)
493  {
494  const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
495 
496  label nRegions = surface.regions().size();
497 
498  // Initialise to global (i.e. per surface)
499  for (label i = 0; i < nRegions; i++)
500  {
501  label globalRegionI = regionOffset_[surfI] + i;
502  normalVolumeTypes_[globalRegionI] = globalVolumeTypes[surfI];
503  if (globalPatchInfo.set(surfI))
504  {
505  patchInfo_.set
506  (
507  globalRegionI,
508  globalPatchInfo[surfI].clone()
509  );
510  }
511  }
512 
513  forAllConstIter(Map<sideVolumeType>, regionVolumeTypes[surfI], iter)
514  {
515  label globalRegionI = regionOffset_[surfI] + iter.key();
516 
517  normalVolumeTypes_[globalRegionI] =
518  regionVolumeTypes[surfI][iter.key()];
519  }
520 
521  const Map<autoPtr<dictionary>>& localInfo = regionPatchInfo[surfI];
522  forAllConstIter(Map<autoPtr<dictionary>>, localInfo, iter)
523  {
524  label globalRegionI = regionOffset_[surfI] + iter.key();
525 
526  patchInfo_.set(globalRegionI, iter()().clone());
527  }
528  }
529 
530 
531 
532  if (!additionalFeaturesDict.empty())
533  {
534  Info<< nl << "Reading additionalFeatures" << endl;
535  }
536 
537  forAllConstIter(dictionary, additionalFeaturesDict, iter)
538  {
539  word featureName = iter().keyword();
540 
541  Info<< nl << " " << iter().keyword() << endl;
542 
543  const dictionary& featureSubDict
544  (
545  additionalFeaturesDict.subDict(featureName)
546  );
547 
548  readFeatures(featureSubDict, featureName, featureI);
549  }
550 
551  // Remove unnecessary space from the features list
552  features_.setSize(featureI);
553 
554  globalBounds_ = treeBoundBox
555  (
556  searchableSurfacesQueries::bounds(allGeometry_, surfaces_)
557  );
558 
559  // Extend the global bounds to stop the bound box sitting on the surfaces
560  // to be conformed to
561  vector newSpan = 1e-4*globalBounds_.span();
562  globalBounds_.min() -= newSpan;
563  globalBounds_.max() += newSpan;
564 
565  // Look at all surfaces at determine whether the locationInMesh point is
566  // inside or outside each, to establish a signature for the domain to be
567  // meshed.
568 
569  referenceVolumeTypes_.setSize
570  (
571  surfaces_.size(),
573  );
574 
575  Info<< endl
576  << "Testing for locationInMesh " << locationInMesh_ << endl;
577 
578  hasBoundedVolume(referenceVolumeTypes_);
579 
580  if (debug)
581  {
582  Info<< "Names = " << allGeometry_.names() << endl;
583  Info<< "Surfaces = " << surfaces_ << endl;
584  Info<< "AllGeom to Surfaces = " << allGeometryToSurfaces_ << endl;
585  Info<< "Volume types = " << normalVolumeTypes_ << endl;
586  Info<< "Patch names = " << patchNames_ << endl;
587  Info<< "Region Offset = " << regionOffset_ << endl;
588 
589  forAll(features_, fI)
590  {
591  Info<< features_[fI].name() << endl;
592  }
593  }
594 }
595 
596 
597 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
598 
600 {}
601 
602 
603 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
604 
605 bool Foam::conformationSurfaces::overlaps(const treeBoundBox& bb) const
606 {
607  forAll(surfaces_, s)
608  {
609  if (allGeometry_[surfaces_[s]].overlaps(bb))
610  {
611  return true;
612  }
613  }
614 
615  return false;
616 }
617 
618 
620 (
621  const pointField& samplePts
622 ) const
623 {
624  return wellInside(samplePts, scalarField(samplePts.size(), 0.0));
625 }
626 
627 
629 (
630  const point& samplePt
631 ) const
632 {
633  return wellInside(pointField(1, samplePt), scalarField(1, 0))[0];
634 }
635 
636 
638 (
639  const pointField& samplePts
640 ) const
641 {
642  return wellOutside(samplePts, scalarField(samplePts.size(), 0.0));
643 }
644 
645 
647 (
648  const point& samplePt
649 ) const
650 {
651  return wellOutside(pointField(1, samplePt), scalarField(1, 0))[0];
652  // return !inside(samplePt);
653 }
654 
655 
657 (
658  const pointField& samplePts,
659  const scalarField& testDistSqr,
660  const bool testForInside
661 ) const
662 {
663  List<List<volumeType>> surfaceVolumeTests
664  (
665  surfaces_.size(),
666  List<volumeType>
667  (
668  samplePts.size(),
670  )
671  );
672 
673  // Get lists for the volumeTypes for each sample wrt each surface
674  forAll(surfaces_, s)
675  {
676  const searchableSurface& surface(allGeometry_[surfaces_[s]]);
677 
678  const label regionI = regionOffset_[s];
679 
680  if (normalVolumeTypes_[regionI] != extendedFeatureEdgeMesh::BOTH)
681  {
682  surface.getVolumeType(samplePts, surfaceVolumeTests[s]);
683  }
684  }
685 
686  // Compare the volumeType result for each point wrt to each surface with the
687  // reference value and if the points are inside the surface by a given
688  // distanceSquared
689 
690  // Assume that the point is wellInside until demonstrated otherwise.
691  Field<bool> insideOutsidePoint(samplePts.size(), testForInside);
692 
693  // Check if the points are inside the surface by the given distance squared
694 
695  labelList hitSurfaces;
696  List<pointIndexHit> hitInfo;
698  (
699  allGeometry_,
700  surfaces_,
701  samplePts,
702  testDistSqr,
703  hitSurfaces,
704  hitInfo
705  );
706 
707  forAll(samplePts, i)
708  {
709  const pointIndexHit& pHit = hitInfo[i];
710 
711  if (pHit.hit())
712  {
713  // If the point is within range of the surface, then it can't be
714  // well (in|out)side
715  insideOutsidePoint[i] = false;
716 
717  continue;
718  }
719 
720  forAll(surfaces_, s)
721  {
722  const label regionI = regionOffset_[s];
723 
724  if (normalVolumeTypes_[regionI] == extendedFeatureEdgeMesh::BOTH)
725  {
726  continue;
727  }
728 
729  const searchableSurface& surface(allGeometry_[surfaces_[s]]);
730 
731  if
732  (
733  !surface.hasVolumeType()
734  //&& surfaceVolumeTests[s][i] == volumeType::unknown
735  )
736  {
737  pointField sample(1, samplePts[i]);
738  scalarField nearestDistSqr(1, great);
739  List<pointIndexHit> info;
740 
741  surface.findNearest(sample, nearestDistSqr, info);
742 
743  vector hitDir = info[0].rawPoint() - samplePts[i];
744  hitDir /= mag(hitDir) + small;
745 
746  pointIndexHit surfHit;
747  label hitSurface;
748 
749  findSurfaceNearestIntersection
750  (
751  samplePts[i],
752  info[0].rawPoint() - 1e-3*mag(hitDir)*hitDir,
753  surfHit,
754  hitSurface
755  );
756 
757  if (surfHit.hit() && hitSurface != surfaces_[s])
758  {
759  continue;
760  }
761  }
762 
763  if (surfaceVolumeTests[s][i] == volumeType::outside)
764  {
765  if
766  (
767  normalVolumeTypes_[regionI]
769  )
770  {
771  insideOutsidePoint[i] = !testForInside;
772  break;
773  }
774  }
775  else if (surfaceVolumeTests[s][i] == volumeType::inside)
776  {
777  if
778  (
779  normalVolumeTypes_[regionI]
781  )
782  {
783  insideOutsidePoint[i] = !testForInside;
784  break;
785  }
786  }
787  }
788  }
789 
790  return insideOutsidePoint;
791 }
792 
793 
795 (
796  const pointField& samplePts,
797  const scalarField& testDistSqr
798 ) const
799 {
800  return wellInOutSide(samplePts, testDistSqr, true);
801 }
802 
803 
805 (
806  const point& samplePt,
807  scalar testDistSqr
808 ) const
809 {
810  return wellInside(pointField(1, samplePt), scalarField(1, testDistSqr))[0];
811 }
812 
813 
815 (
816  const pointField& samplePts,
817  const scalarField& testDistSqr
818 ) const
819 {
820  return wellInOutSide(samplePts, testDistSqr, false);
821 }
822 
823 
825 (
826  const point& samplePt,
827  scalar testDistSqr
828 ) const
829 {
830  return wellOutside(pointField(1, samplePt), scalarField(1, testDistSqr))[0];
831 }
832 
833 
835 (
836  const point& start,
837  const point& end
838 ) const
839 {
840  labelList hitSurfaces;
841  List<pointIndexHit> hitInfo;
842 
844  (
845  allGeometry_,
846  surfaces_,
847  pointField(1, start),
848  pointField(1, end),
849  hitSurfaces,
850  hitInfo
851  );
852 
853  return hitInfo[0].hit();
854 }
855 
856 
858 (
859  const point& start,
860  const point& end,
861  pointIndexHit& surfHit,
862  label& hitSurface
863 ) const
864 {
865  labelList hitSurfaces;
866  List<pointIndexHit> hitInfo;
867 
869  (
870  allGeometry_,
871  surfaces_,
872  pointField(1, start),
873  pointField(1, end),
874  hitSurfaces,
875  hitInfo
876  );
877 
878  surfHit = hitInfo[0];
879 
880  if (surfHit.hit())
881  {
882  // hitSurfaces has returned the index of the entry in surfaces_ that was
883  // found, not the index of the surface in allGeometry_, translating this
884  // to allGeometry_
885 
886  hitSurface = surfaces_[hitSurfaces[0]];
887  }
888 }
889 
890 
892 (
893  const point& start,
894  const point& end,
895  List<pointIndexHit>& surfHit,
896  labelList& hitSurface
897 ) const
898 {
899  labelListList hitSurfaces;
900  List<List<pointIndexHit>> hitInfo;
901 
903  (
904  allGeometry_,
905  surfaces_,
906  pointField(1, start),
907  pointField(1, end),
908  hitSurfaces,
909  hitInfo
910  );
911 
912  surfHit = hitInfo[0];
913 
914  hitSurface.setSize(hitSurfaces[0].size());
915 
916  forAll(hitSurfaces[0], surfI)
917  {
918  // hitSurfaces has returned the index of the entry in surfaces_ that was
919  // found, not the index of the surface in allGeometry_, translating this
920  // to allGeometry_
921 
922  hitSurface[surfI] = surfaces_[hitSurfaces[0][surfI]];
923  }
924 }
925 
926 
928 (
929  const point& start,
930  const point& end,
931  pointIndexHit& surfHit,
932  label& hitSurface
933 ) const
934 {
935  labelList hitSurfacesStart;
936  List<pointIndexHit> hitInfoStart;
937  labelList hitSurfacesEnd;
938  List<pointIndexHit> hitInfoEnd;
939 
941  (
942  allGeometry_,
943  surfaces_,
944  pointField(1, start),
945  pointField(1, end),
946  hitSurfacesStart,
947  hitInfoStart,
948  hitSurfacesEnd,
949  hitInfoEnd
950  );
951 
952  surfHit = hitInfoStart[0];
953 
954  if (surfHit.hit())
955  {
956  // hitSurfaces has returned the index of the entry in surfaces_ that was
957  // found, not the index of the surface in allGeometry_, translating this
958  // to allGeometry_
959 
960  hitSurface = surfaces_[hitSurfacesStart[0]];
961  }
962 }
963 
964 
966 (
967  const point& sample,
968  scalar nearestDistSqr,
969  pointIndexHit& surfHit,
970  label& hitSurface
971 ) const
972 {
973  labelList hitSurfaces;
974  List<pointIndexHit> surfaceHits;
975 
977  (
978  allGeometry_,
979  surfaces_,
980  pointField(1, sample),
981  scalarField(1, nearestDistSqr),
982  hitSurfaces,
983  surfaceHits
984  );
985 
986  surfHit = surfaceHits[0];
987 
988  if (surfHit.hit())
989  {
990  // hitSurfaces has returned the index of the entry in surfaces_ that was
991  // found, not the index of the surface in allGeometry_, translating this
992  // to allGeometry_
993 
994  hitSurface = surfaces_[hitSurfaces[0]];
995  }
996 }
997 
998 
1000 (
1001  const pointField& samples,
1002  const scalarField& nearestDistSqr,
1003  List<pointIndexHit>& surfaceHits,
1004  labelList& hitSurfaces
1005 ) const
1006 {
1008  (
1009  allGeometry_,
1010  surfaces_,
1011  samples,
1012  nearestDistSqr,
1013  hitSurfaces,
1014  surfaceHits
1015  );
1016 
1017  forAll(surfaceHits, i)
1018  {
1019  if (surfaceHits[i].hit())
1020  {
1021  // hitSurfaces has returned the index of the entry in surfaces_ that
1022  // was found, not the index of the surface in allGeometry_,
1023  // translating this to the surface in allGeometry_.
1024 
1025  hitSurfaces[i] = surfaces_[hitSurfaces[i]];
1026  }
1027  }
1028 }
1029 
1030 
1032 (
1033  const point& sample,
1034  scalar nearestDistSqr,
1035  pointIndexHit& fpHit,
1036  label& featureHit
1037 ) const
1038 {
1039  // Work arrays
1040  scalar minDistSqr = nearestDistSqr;
1041  pointIndexHit hitInfo;
1042 
1043  forAll(features_, testI)
1044  {
1045  features_[testI].nearestFeaturePoint
1046  (
1047  sample,
1048  minDistSqr,
1049  hitInfo
1050  );
1051 
1052  if (hitInfo.hit())
1053  {
1054  minDistSqr = magSqr(hitInfo.hitPoint()- sample);
1055  fpHit = hitInfo;
1056  featureHit = testI;
1057  }
1058  }
1059 }
1060 
1061 
1063 (
1064  const point& sample,
1065  scalar nearestDistSqr,
1066  pointIndexHit& edgeHit,
1067  label& featureHit
1068 ) const
1069 {
1070  pointField samples(1, sample);
1071  scalarField nearestDistsSqr(1, nearestDistSqr);
1072 
1073  List<pointIndexHit> edgeHits;
1074  labelList featuresHit;
1075 
1076  findEdgeNearest
1077  (
1078  samples,
1079  nearestDistsSqr,
1080  edgeHits,
1081  featuresHit
1082  );
1083 
1084  edgeHit = edgeHits[0];
1085  featureHit = featuresHit[0];
1086 }
1087 
1088 
1090 (
1091  const pointField& samples,
1092  const scalarField& nearestDistsSqr,
1093  List<pointIndexHit>& edgeHits,
1094  labelList& featuresHit
1095 ) const
1096 {
1097  // Initialise
1098  featuresHit.setSize(samples.size());
1099  featuresHit = -1;
1100  edgeHits.setSize(samples.size());
1101 
1102  // Work arrays
1103  scalarField minDistSqr(nearestDistsSqr);
1104  List<pointIndexHit> hitInfo(samples.size());
1105 
1106  forAll(features_, testI)
1107  {
1108  features_[testI].nearestFeatureEdge
1109  (
1110  samples,
1111  minDistSqr,
1112  hitInfo
1113  );
1114 
1115  // Update minDistSqr and arguments
1116  forAll(hitInfo, pointi)
1117  {
1118  if (hitInfo[pointi].hit())
1119  {
1120  minDistSqr[pointi] = magSqr
1121  (
1122  hitInfo[pointi].hitPoint()
1123  - samples[pointi]
1124  );
1125  edgeHits[pointi] = hitInfo[pointi];
1126  featuresHit[pointi] = testI;
1127  }
1128  }
1129  }
1130 }
1131 
1132 
1134 (
1135  const point& sample,
1136  scalar nearestDistSqr,
1137  List<pointIndexHit>& edgeHits,
1138  List<label>& featuresHit
1139 ) const
1140 {
1141  // Initialise
1142  featuresHit.setSize(extendedFeatureEdgeMesh::nEdgeTypes);
1143  featuresHit = -1;
1144  edgeHits.setSize(extendedFeatureEdgeMesh::nEdgeTypes);
1145 
1146  // Work arrays
1147  scalarField minDistSqr(extendedFeatureEdgeMesh::nEdgeTypes, nearestDistSqr);
1148  List<pointIndexHit> hitInfo(extendedFeatureEdgeMesh::nEdgeTypes);
1149 
1150  forAll(features_, testI)
1151  {
1152  features_[testI].nearestFeatureEdgeByType
1153  (
1154  sample,
1155  minDistSqr,
1156  hitInfo
1157  );
1158 
1159  // Update minDistSqr and arguments
1160  forAll(hitInfo, typeI)
1161  {
1162  if (hitInfo[typeI].hit())
1163  {
1164  minDistSqr[typeI] = magSqr(hitInfo[typeI].hitPoint() - sample);
1165  edgeHits[typeI] = hitInfo[typeI];
1166  featuresHit[typeI] = testI;
1167  }
1168  }
1169  }
1170 }
1171 
1172 
1174 (
1175  const point& sample,
1176  const scalar searchRadiusSqr,
1177  List<List<pointIndexHit>>& edgeHitsByFeature,
1178  List<label>& featuresHit
1179 ) const
1180 {
1181  // Initialise
1182  // featuresHit.setSize(features_.size());
1183  // featuresHit = -1;
1184  // edgeHitsByFeature.setSize(features_.size());
1185 
1186  // Work arrays
1187  List<pointIndexHit> hitInfo(extendedFeatureEdgeMesh::nEdgeTypes);
1188 
1189  forAll(features_, testI)
1190  {
1191  features_[testI].allNearestFeatureEdges
1192  (
1193  sample,
1194  searchRadiusSqr,
1195  hitInfo
1196  );
1197 
1198  bool anyHit = false;
1199  forAll(hitInfo, hitI)
1200  {
1201  if (hitInfo[hitI].hit())
1202  {
1203  anyHit = true;
1204  }
1205  }
1206 
1207  if (anyHit)
1208  {
1209  edgeHitsByFeature.append(hitInfo);
1210  featuresHit.append(testI);
1211  }
1212  }
1213 }
1214 
1215 
1216 void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const
1217 {
1218  OFstream ftStr(runTime_.time().path()/prefix + "_allFeatures.obj");
1219 
1220  Pout<< nl << "Writing all features to " << ftStr.name() << endl;
1221 
1222  label verti = 0;
1223 
1224  forAll(features_, i)
1225  {
1226  const extendedFeatureEdgeMesh& fEM(features_[i]);
1227  const pointField pts(fEM.points());
1228  const edgeList eds(fEM.edges());
1229 
1230  ftStr << "g " << fEM.name() << endl;
1231 
1232  forAll(eds, j)
1233  {
1234  const edge& e = eds[j];
1235 
1236  meshTools::writeOBJ(ftStr, pts[e[0]]); verti++;
1237  meshTools::writeOBJ(ftStr, pts[e[1]]); verti++;
1238  ftStr << "l " << verti-1 << ' ' << verti << endl;
1239  }
1240  }
1241 }
1242 
1243 
1245 (
1246  const point& ptA,
1247  const point& ptB
1248 ) const
1249 {
1250  pointIndexHit surfHit;
1251  label hitSurface;
1252 
1253  findSurfaceAnyIntersection(ptA, ptB, surfHit, hitSurface);
1254 
1255  return getPatchID(hitSurface, surfHit);
1256 }
1257 
1258 
1260 {
1261  pointIndexHit surfHit;
1262  label hitSurface;
1263 
1264  findSurfaceNearest(pt, sqr(great), surfHit, hitSurface);
1265 
1266  return getPatchID(hitSurface, surfHit);
1267 }
1268 
1269 
1271 (
1272  const label hitSurface,
1273  const pointIndexHit& surfHit
1274 ) const
1275 {
1276  if (!surfHit.hit())
1277  {
1278  return -1;
1279  }
1280 
1281  labelList surfLocalRegion;
1282 
1283  allGeometry_[hitSurface].getRegion
1284  (
1285  List<pointIndexHit>(1, surfHit),
1286  surfLocalRegion
1287  );
1288 
1289  const label patchID =
1290  surfLocalRegion[0]
1291  + regionOffset_[allGeometryToSurfaces_[hitSurface]];
1292 
1293  return patchID;
1294 }
1295 
1296 
1299 (
1300  const label hitSurface,
1301  const pointIndexHit& surfHit
1302 ) const
1303 {
1304  const label patchID = getPatchID(hitSurface, surfHit);
1305 
1306  if (patchID == -1)
1307  {
1309  }
1310 
1311  return normalVolumeTypes_[patchID];
1312 }
1313 
1314 
1316 (
1317  const label hitSurface,
1318  const List<pointIndexHit>& surfHit,
1319  vectorField& normal
1320 ) const
1321 {
1322  allGeometry_[hitSurface].getNormal(surfHit, normal);
1323 
1324  const label patchID = regionOffset_[allGeometryToSurfaces_[hitSurface]];
1325 
1326  // Now flip sign of normal depending on mesh side
1327  if (normalVolumeTypes_[patchID] == extendedFeatureEdgeMesh::OUTSIDE)
1328  {
1329  normal *= -1;
1330  }
1331 }
1332 
1333 
1334 // ************************************************************************* //
Field< bool > inside(const pointField &samplePts) const
Check if points are inside surfaces to conform to.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void findEdgeNearestByType(const point &sample, scalar nearestDistSqr, List< pointIndexHit > &edgeHit, List< label > &featuresHit) const
Find the nearest point on each type of feature edge.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool overlaps(const treeBoundBox &bb) const
Check if the supplied bound box overlaps any part of any of.
Field< bool > outside(const pointField &samplePts) const
Check if points are outside surfaces to conform to.
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
conformationSurfaces(const Time &runTime, Random &rndGen, const searchableSurfaces &allGeometry, const dictionary &surfaceConformationDict)
Construct from dictionary and references to conformalVoronoiMesh and.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Field< bool > wellInside(const pointField &samplePts, const scalarField &testDistSqr) const
Check if point is inside surfaces to conform to by at least.
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.
static autoPtr< searchableSurfaceFeatures > New(const searchableSurface &surface, const dictionary &dict)
Return a reference to the selected searchableSurfaceFeatures.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static const Foam::NamedEnum< sideVolumeType, 4 > sideVolumeTypeNames_
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
void findFeaturePointNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &fpHit, label &featureHit) const
Find the nearest point on any feature edge.
extendedFeatureEdgeMesh::sideVolumeType meshableSide(const label hitSurface, const pointIndexHit &surfHit) const
Is the surface a baffle.
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static void findAllIntersections(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelListList &surfaces, List< List< pointIndexHit >> &surfaceHits)
Find all intersections in order from start to end. Returns for.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
scalarField samples(nIntervals, 0)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
T clone(const T &t)
Definition: List.H:54
static void findNearestIntersection(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &surface2, List< pointIndexHit > &hit2)
Find intersections of edge nearest to both endpoints.
List< edge > edgeList
Definition: edgeList.H:38
static label nEdgeTypes
Number of possible feature edge types (i.e. number of slices)
bool findSurfaceAnyIntersection(const point &start, const point &end) const
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:177
void findAllNearestEdges(const point &sample, const scalar searchRadiusSqr, List< List< pointIndexHit >> &edgeHitsByFeature, List< label > &featuresHit) const
Find the nearest points on each feature edge that is within.
void getNormal(const label hitSurface, const List< pointIndexHit > &surfHit, vectorField &normal) const
sideVolumeType
Normals point to the outside.
List< label > labelList
A List of labels.
Definition: labelList.H:56
void writeFeatureObj(const fileName &prefix) const
Write all components of all the extendedFeatureEdgeMeshes as.
static const zero Zero
Definition: zero.H:97
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
void findEdgeNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &edgeHit, label &featureHit) const
Find the nearest point on any feature edge.
Field< bool > wellInOutSide(const pointField &samplePts, const scalarField &testDistSqr, bool testForInside) const
Check if point is closer to the surfaces to conform to than.
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.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:265
defineTypeNameAndDebug(combustionModel, 0)
Field< bool > wellOutside(const pointField &samplePts, const scalarField &testDistSqr) const
Check if point is outside surfaces to conform to by at least.
~conformationSurfaces()
Destructor.
void findSurfaceAllIntersections(const point &start, const point &end, List< pointIndexHit > &surfHit, labelList &hitSurface) const
List< word > wordList
A List of words.
Definition: fileName.H:54
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
void findSurfaceNearestIntersection(const point &start, const point &end, pointIndexHit &surfHit, label &hitSurface) const
Finding the nearestIntersection of the surface to start.
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
static const NamedEnum< volumeType, 4 > names
Definition: volumeType.H:80
label getPatchID(const label hitSurface, const pointIndexHit &surfHit) const
Get the region number of a hit surface.
static boundBox bounds(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest)
Find the boundBox of the selected surfaces.
Namespace for OpenFOAM.
label findPatch(const point &ptA, const point &ptB) const
Find which patch is intersected by the line from one point to.