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