conformalVoronoiMeshConformToSurface.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-2016 OpenFOAM Foundation
9 Copyright (C) 2020-2022 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 "vectorTools.H"
32#include "indexedCellChecks.H"
33#include "IOmanip.H"
34#include "OBJstream.H"
35
36using namespace Foam::vectorTools;
37
38const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle
39 = Foam::cos(degToRad(30.0));
40
41const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface
42 = Foam::cos(degToRad(150.0));
43
44
45// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
46
47void Foam::conformalVoronoiMesh::conformToSurface()
48{
49 this->resetCellCount();
50 // Index the cells
51 for
52 (
53 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
54 cit != finite_cells_end();
55 ++cit
56 )
57 {
58 cit->cellIndex() = Cb::ctUnassigned;
59 }
60
61 if (!reconformToSurface())
62 {
63 // Reinsert stored surface conformation
64 reinsertSurfaceConformation();
65
66 if (Pstream::parRun())
67 {
68 sync(decomposition().procBounds());
69 }
70 }
71 else
72 {
73 ptPairs_.clear();
74
75 // Rebuild, insert and store new surface conformation
76 buildSurfaceConformation();
77
78 if (distributeBackground(*this))
79 {
80 if (Pstream::parRun())
81 {
82 sync(decomposition().procBounds());
83 }
84 }
85
86 // Do not store the surface conformation until after it has been
87 // (potentially) redistributed.
88 storeSurfaceConformation();
89 }
90
91 // reportSurfaceConformationQuality();
92}
93
94
95bool Foam::conformalVoronoiMesh::reconformToSurface() const
96{
97 if
98 (
99 runTime_.timeIndex()
100 % foamyHexMeshControls().surfaceConformationRebuildFrequency() == 0
101 )
102 {
103 return true;
104 }
105
106 return false;
107}
108
109
110// TODO: Investigate topological tests
111Foam::label Foam::conformalVoronoiMesh::findVerticesNearBoundaries()
112{
113 label countNearBoundaryVertices = 0;
114
115 for
116 (
117 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
118 fit != finite_facets_end();
119 ++fit
120 )
121 {
122 Cell_handle c1 = fit->first;
123 Cell_handle c2 = fit->first->neighbor(fit->second);
124
125 if (is_infinite(c1) || is_infinite(c2))
126 {
127 continue;
128 }
129
130 pointFromPoint dE0 = c1->dual();
131 pointFromPoint dE1 = c2->dual();
132
133 if (!geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
134 {
135 continue;
136 }
137
138 for (label celli = 0; celli < 4; ++celli)
139 {
140 Vertex_handle v = c1->vertex(celli);
141
142 if
143 (
144 !is_infinite(v)
145 && v->internalPoint()
146 && fit->second != celli
147 )
148 {
149 v->setNearBoundary();
150 }
151 }
152
153 for (label celli = 0; celli < 4; ++celli)
154 {
155 Vertex_handle v = c2->vertex(celli);
156
157 if
158 (
159 !is_infinite(v)
160 && v->internalPoint()
161 && fit->second != celli
162 )
163 {
164 v->setNearBoundary();
165 }
166 }
167 }
168
169 for
170 (
171 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
172 vit != finite_vertices_end();
173 ++vit
174 )
175 {
176 if (vit->nearBoundary())
177 {
178 countNearBoundaryVertices++;
179 }
180 }
181
182 // Geometric test.
183// for
184// (
185// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
186// vit != finite_vertices_end();
187// ++vit
188// )
189// {
190// if (vit->internalPoint() && !vit->nearBoundary())
191// {
192// pointFromPoint pt = topoint(vit->point());
193//
194// const scalar range = sqr
195// (
196// foamyHexMeshControls().nearBoundaryDistanceCoeff()
197// *targetCellSize(pt)
198// );
199//
200// pointIndexHit pHit;
201// label hitSurface;
202//
203// geometryToConformTo_.findSurfaceNearest
204// (
205// pt,
206// range,
207// pHit,
208// hitSurface
209// );
210//
211// if (pHit.hit())
212// {
213// vit->setNearBoundary();
214// countNearBoundaryVertices++;
215// }
216// }
217// }
218
219 return countNearBoundaryVertices;
220}
221
222
223void Foam::conformalVoronoiMesh::buildSurfaceConformation()
224{
225 timeCheck("Start buildSurfaceConformation");
226
227 Info<< nl
228 << "Rebuilding surface conformation for more iterations"
229 << endl;
230
231 existingEdgeLocations_.clearStorage();
232 existingSurfacePtLocations_.clearStorage();
233
234 buildEdgeLocationTree(existingEdgeLocations_);
235 buildSurfacePtLocationTree(existingSurfacePtLocations_);
236
237 label initialTotalHits = 0;
238
239 // Surface protrusion conformation is done in two steps.
240 // 1. the dual edges (of all internal vertices) can stretch to
241 // 'infinity' so any intersection would be badly behaved. So
242 // just find the nearest point on the geometry and insert point
243 // pairs.
244 // Now most of the surface conformation will be done with some
245 // residual protrusions / incursions.
246 // 2. find any segments of dual edges outside the geometry. Shoot
247 // ray from Delaunay vertex to middle of this segment and introduce
248 // point pairs. This will handle e.g.
249
250 // protruding section of face:
251 //
252 // internal
253 // \ /
254 // -+-----------+-- boundary
255 // \ /
256 // --------
257 //
258 // Shoot ray and find intersection with outside segment (x) and
259 // introduce point pair (..)
260 //
261 // |
262 // \ . /
263 // -+-----|-----+-- boundary
264 // \ . /
265 // ---x----
266
267 // Find vertices near boundaries to speed up subsequent checks.
268 label countNearBoundaryVertices = findVerticesNearBoundaries();
269
270 Info<< " Vertices marked as being near a boundary: "
271 << returnReduce(countNearBoundaryVertices, sumOp<label>())
272 << " (estimated)" << endl;
273
274 timeCheck("After set near boundary");
275
276 const scalar edgeSearchDistCoeffSqr =
277 foamyHexMeshControls().edgeSearchDistCoeffSqr();
278
279 const scalar surfacePtReplaceDistCoeffSqr =
280 foamyHexMeshControls().surfacePtReplaceDistCoeffSqr();
281
282 const label AtoV = label(6/Foam::pow(scalar(number_of_vertices()), 3));
283
284 // Initial surface protrusion conformation - nearest surface point
285 {
286 pointIndexHitAndFeatureDynList featureEdgeHits(AtoV/4);
287 pointIndexHitAndFeatureDynList surfaceHits(AtoV);
288 DynamicList<label> edgeToTreeShape(AtoV/4);
289 DynamicList<label> surfaceToTreeShape(AtoV);
290
291 Map<scalar> surfacePtToEdgePtDist(AtoV/4);
292
293 for
294 (
295 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
296 vit != finite_vertices_end();
297 vit++
298 )
299 {
300 if (vit->nearBoundary())
301 {
302 pointIndexHitAndFeatureDynList surfaceIntersections(AtoV);
303
304 if
305 (
306 dualCellSurfaceAllIntersections
307 (
308 vit,
309 surfaceIntersections
310 )
311 )
312 {
313 // meshTools::writeOBJ(Pout, vert);
314 // meshTools::writeOBJ(Pout, surfHit.hitPoint());
315 // Pout<< "l cr0 cr1" << endl;
316
317 addSurfaceAndEdgeHits
318 (
319 topoint(vit->point()),
320 surfaceIntersections,
321 surfacePtReplaceDistCoeffSqr,
322 edgeSearchDistCoeffSqr,
323 surfaceHits,
324 featureEdgeHits,
325 surfaceToTreeShape,
326 edgeToTreeShape,
327 surfacePtToEdgePtDist,
328 true
329 );
330 }
331 else
332 {
333 vit->setInternal();
334 countNearBoundaryVertices--;
335 }
336 }
337 }
338
339 Info<< " Vertices marked as being near a boundary: "
340 << returnReduce(countNearBoundaryVertices, sumOp<label>())
341 << " (after dual surface intersection)" << endl;
342
343 label nVerts = number_of_vertices();
344 label nSurfHits = surfaceHits.size();
345 label nFeatEdHits = featureEdgeHits.size();
346
347 if (Pstream::parRun())
348 {
349 reduce(nVerts, sumOp<label>());
350 reduce(nSurfHits, sumOp<label>());
351 reduce(nFeatEdHits, sumOp<label>());
352 }
353
354 Info<< nl << "Initial conformation" << nl
355 << " Number of vertices " << nVerts << nl
356 << " Number of surface hits " << nSurfHits << nl
357 << " Number of edge hits " << nFeatEdHits
358 << endl;
359
360 // In parallel, synchronise the surface trees
361 if (Pstream::parRun())
362 {
363 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
364 }
365
366 DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
367
368 insertSurfacePointPairs
369 (
370 surfaceHits,
371 "surfaceConformationLocations_initial.obj",
372 pts
373 );
374
375 // In parallel, synchronise the edge trees
376 if (Pstream::parRun())
377 {
378 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
379 }
380
381 insertEdgePointGroups
382 (
383 featureEdgeHits,
384 "edgeConformationLocations_initial.obj",
385 pts
386 );
387
388 pts.shrink();
389
390 Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
391
392 // Re-index the point pairs
393 ptPairs_.reIndex(oldToNewIndices);
394
395 //writePointPairs("pointPairs_initial.obj");
396
397 // Remove location from surface/edge tree
398
399 timeCheck("After initial conformation");
400
401 initialTotalHits = nSurfHits + nFeatEdHits;
402 }
403
404 // Remember which vertices were referred to each processor so only updates
405 // are sent.
406 PtrList<labelPairHashSet> referralVertices(Pstream::nProcs());
407
408 // Store the vertices that have been received and added from each processor
409 // already so that there is no attempt to add them more than once.
410 autoPtr<labelPairHashSet> receivedVertices;
411
412 if (Pstream::parRun())
413 {
414 forAll(referralVertices, proci)
415 {
416 if (proci != Pstream::myProcNo())
417 {
418 referralVertices.set
419 (
420 proci,
421 new labelPairHashSet(number_of_vertices()/Pstream::nProcs())
422 );
423 }
424 }
425
426 receivedVertices.reset
427 (
428 new labelPairHashSet(number_of_vertices()/Pstream::nProcs())
429 );
430
431 // Build the parallel interface the initial surface conformation
432 sync
433 (
434 decomposition_().procBounds(),
435 referralVertices,
436 receivedVertices()
437 );
438 }
439
440 label iterationNo = 0;
441
442 label maxIterations = foamyHexMeshControls().maxConformationIterations();
443
444 scalar iterationToInitialHitRatioLimit =
445 foamyHexMeshControls().iterationToInitialHitRatioLimit();
446
447 label hitLimit = label(iterationToInitialHitRatioLimit*initialTotalHits);
448
449 Info<< nl << "Stopping iterations when: " << nl
450 << " total number of hits drops below "
451 << iterationToInitialHitRatioLimit
452 << " of initial hits (" << hitLimit << ")" << nl
453 << " or " << nl
454 << " maximum number of iterations (" << maxIterations
455 << ") is reached"
456 << endl;
457
458 // Set totalHits to a large enough positive value to enter the while loop on
459 // the first iteration
460 label totalHits = initialTotalHits;
461
462 while
463 (
464 totalHits > 0
465 && totalHits >= hitLimit
466 && iterationNo < maxIterations
467 )
468 {
469 pointIndexHitAndFeatureDynList surfaceHits(0.5*AtoV);
470 pointIndexHitAndFeatureDynList featureEdgeHits(0.25*AtoV);
471 DynamicList<label> surfaceToTreeShape(AtoV/2);
472 DynamicList<label> edgeToTreeShape(AtoV/4);
473
474 Map<scalar> surfacePtToEdgePtDist;
475
476 for
477 (
478 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
479 vit != finite_vertices_end();
480 ++vit
481 )
482 {
483 // The initial surface conformation has already identified the
484 // nearBoundary set of vertices. Previously inserted boundary
485 // points and referred internal vertices from other processors can
486 // also generate protrusions and must be assessed too.
487 if
488 (
489 vit->nearBoundary()
490 || vit->internalBoundaryPoint()
491 || (vit->internalOrBoundaryPoint() && vit->referred())
492 )
493 {
494 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
495
496 pointIndexHit surfHit;
497 label hitSurface;
498
499 // Find segments of dual face outside the geometry and find the
500 // the middle of this
501 dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
502
503 if (surfHit.hit())
504 {
505 surfaceIntersections.append
506 (
507 pointIndexHitAndFeature(surfHit, hitSurface)
508 );
509
510 addSurfaceAndEdgeHits
511 (
512 topoint(vit->point()),
513 surfaceIntersections,
514 surfacePtReplaceDistCoeffSqr,
515 edgeSearchDistCoeffSqr,
516 surfaceHits,
517 featureEdgeHits,
518 surfaceToTreeShape,
519 edgeToTreeShape,
520 surfacePtToEdgePtDist,
521 false
522 );
523 }
524 else
525 {
526 // No surface hit detected so if internal then don't check
527 // again
528 if (vit->nearBoundary())
529 {
530 vit->setInternal();
531 }
532 }
533 }
534 else if
535 (
536 vit->externalBoundaryPoint()
537 || (vit->externalBoundaryPoint() && vit->referred())
538 )
539 {
540 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
541
542 pointIndexHit surfHit;
543 label hitSurface;
544
545 // Detect slave (external vertices) whose dual face incurs
546 // into nearby (other than originating) geometry
547 dualCellLargestSurfaceIncursion(vit, surfHit, hitSurface);
548
549 if (surfHit.hit())
550 {
551 surfaceIntersections.append
552 (
553 pointIndexHitAndFeature(surfHit, hitSurface)
554 );
555
556 addSurfaceAndEdgeHits
557 (
558 topoint(vit->point()),
559 surfaceIntersections,
560 surfacePtReplaceDistCoeffSqr,
561 edgeSearchDistCoeffSqr,
562 surfaceHits,
563 featureEdgeHits,
564 surfaceToTreeShape,
565 edgeToTreeShape,
566 surfacePtToEdgePtDist,
567 false
568 );
569 }
570 }
571 }
572
573 label nVerts = number_of_vertices();
574 label nSurfHits = surfaceHits.size();
575 label nFeatEdHits = featureEdgeHits.size();
576
577 if (Pstream::parRun())
578 {
579 reduce(nVerts, sumOp<label>());
580 reduce(nSurfHits, sumOp<label>());
581 reduce(nFeatEdHits, sumOp<label>());
582 }
583
584 Info<< nl << "Conformation iteration " << iterationNo << nl
585 << " Number of vertices " << nVerts << nl
586 << " Number of surface hits " << nSurfHits << nl
587 << " Number of edge hits " << nFeatEdHits
588 << endl;
589
590 totalHits = nSurfHits + nFeatEdHits;
591
592 label nNotInserted = 0;
593
594 if (totalHits > 0)
595 {
596 // In parallel, synchronise the surface trees
597 if (Pstream::parRun())
598 {
599 nNotInserted +=
600 synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
601 }
602
604 (
605 2*surfaceHits.size() + 3*featureEdgeHits.size()
606 );
607
608 insertSurfacePointPairs
609 (
610 surfaceHits,
611 "surfaceConformationLocations_" + name(iterationNo) + ".obj",
612 pts
613 );
614
615 // In parallel, synchronise the edge trees
616 if (Pstream::parRun())
617 {
618 nNotInserted +=
619 synchroniseEdgeTrees(edgeToTreeShape, featureEdgeHits);
620 }
621
622 insertEdgePointGroups
623 (
624 featureEdgeHits,
625 "edgeConformationLocations_" + name(iterationNo) + ".obj",
626 pts
627 );
628
629 pts.shrink();
630
631 Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
632
633 // Reindex the point pairs
634 ptPairs_.reIndex(oldToNewIndices);
635
636 //writePointPairs("pointPairs_" + name(iterationNo) + ".obj");
637
638 if (Pstream::parRun())
639 {
640 sync
641 (
642 decomposition_().procBounds(),
643 referralVertices,
644 receivedVertices()
645 );
646 }
647 }
648
649 timeCheck("Conformation iteration " + name(iterationNo));
650
651 iterationNo++;
652
653 if (iterationNo == maxIterations)
654 {
656 << "Maximum surface conformation iterations ("
657 << maxIterations << ") reached." << endl;
658 }
659
660 if (totalHits <= nNotInserted)
661 {
662 Info<< nl << "Total hits (" << totalHits
663 << ") less than number of failed insertions (" << nNotInserted
664 << "), stopping iterations" << endl;
665 break;
666 }
667
668 if (totalHits < hitLimit)
669 {
670 Info<< nl << "Total hits (" << totalHits
671 << ") less than limit (" << hitLimit
672 << "), stopping iterations" << endl;
673 }
674 }
675
676 edgeLocationTreePtr_.clear();
677 surfacePtLocationTreePtr_.clear();
678}
679
680
681Foam::label Foam::conformalVoronoiMesh::synchroniseSurfaceTrees
682(
683 const DynamicList<label>& surfaceToTreeShape,
684 pointIndexHitAndFeatureList& surfaceHits
685)
686{
687 Info<< " Surface tree synchronisation" << endl;
688
689 pointIndexHitAndFeatureDynList synchronisedSurfLocations
690 (
691 surfaceHits.size()
692 );
693
695 procSurfLocations[Pstream::myProcNo()] = surfaceHits;
696 Pstream::allGatherList(procSurfLocations);
697
699
700 label nStoppedInsertion = 0;
701
702 // Do the nearness tests here
703 for (const int proci : Pstream::allProcs())
704 {
705 // Skip own points
706 if (proci >= Pstream::myProcNo())
707 {
708 continue;
709 }
710
711 const pointIndexHitAndFeatureList& otherSurfEdges =
712 procSurfLocations[proci];
713
714 forAll(otherSurfEdges, peI)
715 {
716 const Foam::point& pt = otherSurfEdges[peI].first().hitPoint();
717
718 pointIndexHit nearest;
719 pointIsNearSurfaceLocation(pt, nearest);
720
721 pointIndexHit nearestEdge;
722 pointIsNearFeatureEdgeLocation(pt, nearestEdge);
723
724 if (nearest.hit() || nearestEdge.hit())
725 {
726 nStoppedInsertion++;
727 hits[proci].insert(peI);
728 }
729 }
730 }
731
733
734 forAll(surfaceHits, eI)
735 {
736 if (!hits[Pstream::myProcNo()].found(eI))
737 {
738 synchronisedSurfLocations.append(surfaceHits[eI]);
739 }
740 else
741 {
742 surfacePtLocationTreePtr_().remove(surfaceToTreeShape[eI]);
743 }
744 }
745
746// forAll(synchronisedSurfLocations, pI)
747// {
748// appendToSurfacePtTree
749// (
750// synchronisedSurfLocations[pI].first().hitPoint()
751// );
752// }
753
754 const label nNotInserted = returnReduce(nStoppedInsertion, sumOp<label>());
755
756 Info<< " Not inserting total of " << nNotInserted << " locations"
757 << endl;
758
759 surfaceHits = synchronisedSurfLocations;
760
761 return nNotInserted;
762}
763
764
765Foam::label Foam::conformalVoronoiMesh::synchroniseEdgeTrees
766(
767 const DynamicList<label>& edgeToTreeShape,
768 pointIndexHitAndFeatureList& featureEdgeHits
769)
770{
771 Info<< " Edge tree synchronisation" << endl;
772
773 pointIndexHitAndFeatureDynList synchronisedEdgeLocations
774 (
775 featureEdgeHits.size()
776 );
777
779 procEdgeLocations[Pstream::myProcNo()] = featureEdgeHits;
780 Pstream::allGatherList(procEdgeLocations);
781
783
784 label nStoppedInsertion = 0;
785
786 // Do the nearness tests here
787 for (const int proci : Pstream::allProcs())
788 {
789 // Skip own points
790 if (proci >= Pstream::myProcNo())
791 {
792 continue;
793 }
794
795 pointIndexHitAndFeatureList& otherProcEdges = procEdgeLocations[proci];
796
797 forAll(otherProcEdges, peI)
798 {
799 const Foam::point& pt = otherProcEdges[peI].first().hitPoint();
800
801 pointIndexHit nearest;
802 pointIsNearFeatureEdgeLocation(pt, nearest);
803
804 if (nearest.hit())
805 {
806// Pout<< "Not inserting " << peI << " " << pt << " "
807// << nearest.rawPoint() << " on proc " << proci
808// << ", near edge = " << nearest
809// << " near ftPt = "<< info
810// << " " << featureEdgeExclusionDistanceSqr(pt)
811// << endl;
812
813 nStoppedInsertion++;
814 hits[proci].insert(peI);
815 }
816 }
817 }
818
820
821 forAll(featureEdgeHits, eI)
822 {
823 if (!hits[Pstream::myProcNo()].found(eI))
824 {
825 synchronisedEdgeLocations.append(featureEdgeHits[eI]);
826 }
827 else
828 {
829 edgeLocationTreePtr_().remove(edgeToTreeShape[eI]);
830 }
831 }
832
833// forAll(synchronisedEdgeLocations, pI)
834// {
835// appendToEdgeLocationTree
836// (
837// synchronisedEdgeLocations[pI].first().hitPoint()
838// );
839// }
840
841 const label nNotInserted = returnReduce(nStoppedInsertion, sumOp<label>());
842
843 Info<< " Not inserting total of " << nNotInserted << " locations"
844 << endl;
845
846 featureEdgeHits = synchronisedEdgeLocations;
847
848 return nNotInserted;
849}
850
851
852bool Foam::conformalVoronoiMesh::surfaceLocationConformsToInside
853(
854 const pointIndexHitAndFeature& info
855) const
856{
857 if (info.first().hit())
858 {
859 vectorField norm(1);
860
861 geometryToConformTo_.getNormal
862 (
863 info.second(),
864 List<pointIndexHit>(1, info.first()),
865 norm
866 );
867
868 const vector& n = norm[0];
869
870 const scalar ppDist = pointPairDistance(info.first().hitPoint());
871
872 const Foam::point innerPoint = info.first().hitPoint() - ppDist*n;
873
874 if (!geometryToConformTo_.inside(innerPoint))
875 {
876 return false;
877 }
878
879 return true;
880 }
881
882 return false;
883}
884
885
886bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
887(
888 const Delaunay::Finite_vertices_iterator& vit
889) const
890{
891 std::list<Facet> facets;
892 incident_facets(vit, std::back_inserter(facets));
893
894 for
895 (
896 std::list<Facet>::iterator fit=facets.begin();
897 fit != facets.end();
898 ++fit
899 )
900 {
901 if
902 (
903 is_infinite(fit->first)
904 || is_infinite(fit->first->neighbor(fit->second))
905 || !fit->first->hasInternalPoint()
906 || !fit->first->neighbor(fit->second)->hasInternalPoint()
907 )
908 {
909 continue;
910 }
911
912 Foam::point dE0 = fit->first->dual();
913 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
914
915 if (Pstream::parRun())
916 {
917 Foam::point& a = dE0;
918 Foam::point& b = dE1;
919
920 bool inProc = clipLineToProc(topoint(vit->point()), a, b);
921
922 // Check for the edge passing through a surface
923 if
924 (
925 inProc
926 && geometryToConformTo_.findSurfaceAnyIntersection(a, b)
927 )
928 {
929 return true;
930 }
931 }
932 else
933 {
934 if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
935 {
936 return true;
937 }
938 }
939 }
940
941 return false;
942}
943
944
945bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
946(
947 const Delaunay::Finite_vertices_iterator& vit,
948 pointIndexHitAndFeatureDynList& infoList
949) const
950{
951 bool flagIntersection = false;
952
953 std::list<Facet> facets;
954 incident_facets(vit, std::back_inserter(facets));
955
956 for
957 (
958 std::list<Facet>::iterator fit = facets.begin();
959 fit != facets.end();
960 ++fit
961 )
962 {
963 if
964 (
965 is_infinite(fit->first)
966 || is_infinite(fit->first->neighbor(fit->second))
967 || !fit->first->hasInternalPoint()
968 || !fit->first->neighbor(fit->second)->hasInternalPoint()
969 )
970 {
971 continue;
972 }
973
974 // Construct the dual edge and search for intersections of the edge
975 // with the surface
976 Foam::point dE0 = fit->first->dual();
977 Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
978
979 pointIndexHit infoIntersection;
980 label hitSurfaceIntersection = -1;
981
982 if (Pstream::parRun())
983 {
984 bool inProc = clipLineToProc(topoint(vit->point()), dE0, dE1);
985
986 if (!inProc)
987 {
988 continue;
989 }
990 }
991
992 geometryToConformTo_.findSurfaceNearestIntersection
993 (
994 dE0,
995 dE1,
996 infoIntersection,
997 hitSurfaceIntersection
998 );
999
1000 if (infoIntersection.hit())
1001 {
1002 vectorField norm(1);
1003
1004 geometryToConformTo_.getNormal
1005 (
1006 hitSurfaceIntersection,
1007 List<pointIndexHit>(1, infoIntersection),
1008 norm
1009 );
1010
1011 const vector& n = norm[0];
1012
1013 pointFromPoint vertex = topoint(vit->point());
1014
1015 const plane p(infoIntersection.hitPoint(), n);
1016
1017 const plane::ray r(vertex, n);
1018
1019 const scalar d = p.normalIntersect(r);
1020
1021 Foam::point newPoint = vertex + d*n;
1022
1023 pointIndexHitAndFeature info;
1024 geometryToConformTo_.findSurfaceNearest
1025 (
1026 newPoint,
1027 4.0*magSqr(newPoint - vertex),
1028 info.first(),
1029 info.second()
1030 );
1031
1032 bool rejectPoint = false;
1033
1034 if (!surfaceLocationConformsToInside(info))
1035 {
1036 rejectPoint = true;
1037 }
1038
1039 if (!rejectPoint && info.first().hit())
1040 {
1041 if (!infoList.empty())
1042 {
1043 forAll(infoList, hitI)
1044 {
1045 // Reject point if the point is already added
1046 if
1047 (
1048 infoList[hitI].first().index()
1049 == info.first().index()
1050 )
1051 {
1052 rejectPoint = true;
1053 break;
1054 }
1055
1056 const Foam::point& p
1057 = infoList[hitI].first().hitPoint();
1058
1059 const scalar separationDistance =
1060 mag(p - info.first().hitPoint());
1061
1062 const scalar minSepDist =
1063 sqr
1064 (
1065 foamyHexMeshControls().removalDistCoeff()
1066 *targetCellSize(p)
1067 );
1068
1069 // Reject the point if it is too close to another
1070 // surface point.
1071 // Could merge the points?
1072 if (separationDistance < minSepDist)
1073 {
1074 rejectPoint = true;
1075 break;
1076 }
1077 }
1078 }
1079 }
1080
1081 // The normal ray from the vertex will not always result in a hit
1082 // because another surface may be in the way.
1083 if (!rejectPoint && info.first().hit())
1084 {
1085 flagIntersection = true;
1086 infoList.append(info);
1087 }
1088 }
1089 }
1090
1091 return flagIntersection;
1092}
1093
1094
1095bool Foam::conformalVoronoiMesh::clipLineToProc
1096(
1097 const Foam::point& pt,
1098 Foam::point& a,
1099 Foam::point& b
1100) const
1101{
1102 bool inProc = false;
1103
1104 pointIndexHit findAnyIntersection = decomposition_().findLine(a, b);
1105
1106 if (!findAnyIntersection.hit())
1107 {
1108 pointIndexHit info = decomposition_().findLine(a, pt);
1109
1110 if (!info.hit())
1111 {
1112 inProc = true;
1113 }
1114 else
1115 {
1116 inProc = false;
1117 }
1118 }
1119 else
1120 {
1121 pointIndexHit info = decomposition_().findLine(a, pt);
1122
1123 if (!info.hit())
1124 {
1125 inProc = true;
1126 b = findAnyIntersection.hitPoint();
1127 }
1128 else
1129 {
1130 inProc = true;
1131 a = findAnyIntersection.hitPoint();
1132 }
1133 }
1134
1135 return inProc;
1136}
1137
1138
1139void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
1140(
1141 const Delaunay::Finite_vertices_iterator& vit,
1142 pointIndexHit& surfHitLargest,
1143 label& hitSurfaceLargest
1144) const
1145{
1146 // Set no-hit data
1147 surfHitLargest = pointIndexHit();
1148 hitSurfaceLargest = -1;
1149
1150 std::list<Facet> facets;
1151 finite_incident_facets(vit, std::back_inserter(facets));
1152
1153 pointFromPoint vert = topoint(vit->point());
1154
1155 scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
1156
1157 for
1158 (
1159 std::list<Facet>::iterator fit = facets.begin();
1160 fit != facets.end();
1161 ++fit
1162 )
1163 {
1164 Cell_handle c1 = fit->first;
1165 Cell_handle c2 = fit->first->neighbor(fit->second);
1166
1167 if
1168 (
1169 is_infinite(c1) || is_infinite(c2)
1170 || (
1171 !c1->internalOrBoundaryDualVertex()
1172 || !c2->internalOrBoundaryDualVertex()
1173 )
1174 || !c1->real() || !c2->real()
1175 )
1176 {
1177 continue;
1178 }
1179
1180// Foam::point endPt = 0.5*(c1->dual() + c2->dual());
1181 Foam::point endPt = c1->dual();
1182
1183 if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
1184 {
1185 endPt = c2->dual();
1186 }
1187
1188 if
1189 (
1190 magSqr(vert - endPt)
1191 > magSqr(geometryToConformTo().globalBounds().mag())
1192 )
1193 {
1194 continue;
1195 }
1196
1197 pointIndexHit surfHit;
1198 label hitSurface;
1199
1200 geometryToConformTo_.findSurfaceNearestIntersection
1201 (
1202 vert,
1203 endPt,
1204 surfHit,
1205 hitSurface
1206 );
1207
1208 if (surfHit.hit())
1209 {
1210 vectorField norm(1);
1211
1212 allGeometry_[hitSurface].getNormal
1213 (
1214 List<pointIndexHit>(1, surfHit),
1215 norm
1216 );
1217
1218 const vector& n = norm[0];
1219
1220 const scalar normalProtrusionDistance
1221 (
1222 (endPt - surfHit.hitPoint()) & n
1223 );
1224
1225 if (normalProtrusionDistance > maxProtrusionDistance)
1226 {
1227 const plane p(surfHit.hitPoint(), n);
1228
1229 const plane::ray r(endPt, -n);
1230
1231 const scalar d = p.normalIntersect(r);
1232
1233 Foam::point newPoint = endPt - d*n;
1234
1235 pointIndexHitAndFeature info;
1236 geometryToConformTo_.findSurfaceNearest
1237 (
1238 newPoint,
1239 4.0*magSqr(newPoint - endPt),
1240 info.first(),
1241 info.second()
1242 );
1243
1244 if (info.first().hit())
1245 {
1246 if
1247 (
1248 surfaceLocationConformsToInside
1249 (
1250 pointIndexHitAndFeature(info.first(), info.second())
1251 )
1252 )
1253 {
1254 surfHitLargest = info.first();
1255 hitSurfaceLargest = info.second();
1256
1257 maxProtrusionDistance = normalProtrusionDistance;
1258 }
1259 }
1260 }
1261 }
1262 }
1263
1264 // Relying on short-circuit evaluation to not call for hitPoint when this
1265 // is a miss
1266 if
1267 (
1268 surfHitLargest.hit()
1269 && (
1271 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1272 )
1273 )
1274 {
1275 // A protrusion was identified, but not penetrating on this processor,
1276 // so set no-hit data and allow the other that should have this point
1277 // referred to generate it.
1278 surfHitLargest = pointIndexHit();
1279 hitSurfaceLargest = -1;
1280 }
1281}
1282
1283
1284void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
1285(
1286 const Delaunay::Finite_vertices_iterator& vit,
1287 pointIndexHit& surfHitLargest,
1288 label& hitSurfaceLargest
1289) const
1290{
1291 // Set no-hit data
1292 surfHitLargest = pointIndexHit();
1293 hitSurfaceLargest = -1;
1294
1295 std::list<Facet> facets;
1296 finite_incident_facets(vit, std::back_inserter(facets));
1297
1298 pointFromPoint vert = topoint(vit->point());
1299
1300 scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
1301
1302 for
1303 (
1304 std::list<Facet>::iterator fit = facets.begin();
1305 fit != facets.end();
1306 ++fit
1307 )
1308 {
1309 Cell_handle c1 = fit->first;
1310 Cell_handle c2 = fit->first->neighbor(fit->second);
1311
1312 if
1313 (
1314 is_infinite(c1) || is_infinite(c2)
1315 || (
1316 !c1->internalOrBoundaryDualVertex()
1317 || !c2->internalOrBoundaryDualVertex()
1318 )
1319 || !c1->real() || !c2->real()
1320 )
1321 {
1322 continue;
1323 }
1324
1325// Foam::point endPt = 0.5*(c1->dual() + c2->dual());
1326 Foam::point endPt = c1->dual();
1327
1328 if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
1329 {
1330 endPt = c2->dual();
1331 }
1332
1333 if
1334 (
1335 magSqr(vert - endPt)
1336 > magSqr(geometryToConformTo().globalBounds().mag())
1337 )
1338 {
1339 continue;
1340 }
1341
1342 pointIndexHit surfHit;
1343 label hitSurface;
1344
1345 geometryToConformTo_.findSurfaceNearestIntersection
1346 (
1347 vert,
1348 endPt,
1349 surfHit,
1350 hitSurface
1351 );
1352
1353 if (surfHit.hit())
1354 {
1355 vectorField norm(1);
1356
1357 allGeometry_[hitSurface].getNormal
1358 (
1359 List<pointIndexHit>(1, surfHit),
1360 norm
1361 );
1362
1363 const vector& n = norm[0];
1364
1365 scalar normalIncursionDistance
1366 (
1367 (endPt - surfHit.hitPoint()) & n
1368 );
1369
1370 if (normalIncursionDistance < minIncursionDistance)
1371 {
1372 const plane p(surfHit.hitPoint(), n);
1373
1374 const plane::ray r(endPt, n);
1375
1376 const scalar d = p.normalIntersect(r);
1377
1378 Foam::point newPoint = endPt + d*n;
1379
1380 pointIndexHitAndFeature info;
1381 geometryToConformTo_.findSurfaceNearest
1382 (
1383 newPoint,
1384 4.0*magSqr(newPoint - endPt),
1385 info.first(),
1386 info.second()
1387 );
1388
1389 if (info.first().hit())
1390 {
1391 if
1392 (
1393 surfaceLocationConformsToInside
1394 (
1395 pointIndexHitAndFeature(info.first(), info.second())
1396 )
1397 )
1398 {
1399 surfHitLargest = info.first();
1400 hitSurfaceLargest = info.second();
1401
1402 minIncursionDistance = normalIncursionDistance;
1403 }
1404 }
1405 }
1406 }
1407 }
1408
1409 // Relying on short-circuit evaluation to not call for hitPoint when this
1410 // is a miss
1411 if
1412 (
1413 surfHitLargest.hit()
1414 && (
1416 && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
1417 )
1418 )
1419 {
1420 // A protrusion was identified, but not penetrating on this processor,
1421 // so set no-hit data and allow the other that should have this point
1422 // referred to generate it.
1423 surfHitLargest = pointIndexHit();
1424 hitSurfaceLargest = -1;
1425 }
1426}
1427
1428
1429void Foam::conformalVoronoiMesh::reportProcessorOccupancy()
1430{
1431 for
1432 (
1433 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1434 vit != finite_vertices_end();
1435 vit++
1436 )
1437 {
1438 if (vit->real())
1439 {
1440 if
1441 (
1443 && !decomposition().positionOnThisProcessor(topoint(vit->point()))
1444 )
1445 {
1446 Pout<< topoint(vit->point()) << " is not on this processor "
1447 << endl;
1448 }
1449 }
1450 }
1451}
1452
1453
1454//void Foam::conformalVoronoiMesh::reportSurfaceConformationQuality()
1455//{
1456// Info<< nl << "Check surface conformation quality" << endl;
1457//
1458// for
1459// (
1460// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1461// vit != finite_vertices_end();
1462// vit++
1463// )
1464// {
1465// if (vit->internalOrBoundaryPoint())
1466// {
1467// Foam::point vert(topoint(vit->point()));
1468// pointIndexHit surfHit;
1469// label hitSurface;
1470//
1471// dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
1472//
1473// if (surfHit.hit())
1474// {
1475// Pout<< nl << "Residual penetration: " << nl
1476// << vit->index() << nl
1477// << vit->type() << nl
1478// << vit->ppMaster() << nl
1479// << "nearFeaturePt "
1480// << nearFeaturePt(surfHit.hitPoint()) << nl
1481// << vert << nl
1482// << surfHit.hitPoint()
1483// << endl;
1484// }
1485// }
1486// }
1487//
1488// {
1489// // Assess close surface points
1490//
1491// setVertexSizeAndAlignment();
1492//
1493// for
1494// (
1495// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1496// vit != finite_vertices_end();
1497// vit++
1498// )
1499// {
1500// if (vit->ppMaster())
1501// {
1502// std::list<Vertex_handle> adjacentVertices;
1503//
1504// adjacent_vertices(vit, std::back_inserter(adjacentVertices));
1505//
1506// Foam::point pt = topoint(vit->point());
1507//
1508// // Pout<< nl << "vit: " << vit->index() << " "
1509// // << topoint(vit->point())
1510// // << endl;
1511//
1512// // Pout<< adjacentVertices.size() << endl;
1513//
1514// for
1515// (
1516// std::list<Vertex_handle>::iterator
1517// avit = adjacentVertices.begin();
1518// avit != adjacentVertices.end();
1519// ++avit
1520// )
1521// {
1522// Vertex_handle avh = *avit;
1523//
1524// // The lower indexed vertex will perform the assessment
1525// if
1526// (
1527// avh->ppMaster()
1528// && vit->index() < avh->index()
1529// && vit->type() != avh->type()
1530// )
1531// {
1532// scalar targetSize = 0.2*averageAnyCellSize(vit, avh);
1533//
1534// // Pout<< "diff " << mag(pt - topoint(avh->point()))
1535// // << " " << targetSize << endl;
1536//
1537// if
1538// (
1539// magSqr(pt - topoint(avh->point()))
1540// < sqr(targetSize)
1541// )
1542// {
1543// Pout<< nl << "vit: " << vit->index() << " "
1544// << topoint(vit->point())
1545// << endl;
1546//
1547// Pout<< " adjacent too close: "
1548// << avh->index() << " "
1549// << topoint(avh->point())
1550// << endl;
1551// }
1552// }
1553// }
1554// }
1555// }
1556// }
1557//}
1558
1559void Foam::conformalVoronoiMesh::limitDisplacement
1560(
1561 const Delaunay::Finite_vertices_iterator& vit,
1562 vector& displacement,
1563 label callCount
1564) const
1565{
1566 callCount++;
1567
1568 // Do not allow infinite recursion
1569 if (callCount > 7)
1570 {
1571 displacement = Zero;
1572 return;
1573 }
1574
1575 pointFromPoint pt = topoint(vit->point());
1576 Foam::point dispPt = pt + displacement;
1577
1578 bool limit = false;
1579
1580 pointIndexHit surfHit;
1581 label hitSurface;
1582
1583 if (!geometryToConformTo_.globalBounds().contains(dispPt))
1584 {
1585 // If dispPt is outside bounding box then displacement cuts boundary
1586 limit = true;
1587 }
1588 else if (geometryToConformTo_.findSurfaceAnyIntersection(pt, dispPt))
1589 {
1590 // Full surface penetration test
1591 limit = true;
1592 }
1593 else
1594 {
1595 // Testing if the displaced position is too close to the surface.
1596 // Within twice the local surface point pair insertion distance is
1597 // considered "too close"
1598
1599 scalar searchDistanceSqr = sqr
1600 (
1601 2*vit->targetCellSize()
1602 *foamyHexMeshControls().pointPairDistanceCoeff()
1603 );
1604
1605 geometryToConformTo_.findSurfaceNearest
1606 (
1607 dispPt,
1608 searchDistanceSqr,
1609 surfHit,
1610 hitSurface
1611 );
1612
1613 if (surfHit.hit())
1614 {
1615 limit = true;
1616
1617 if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
1618 {
1619 // Cannot limit displacement, point closer than tolerance
1620 displacement = Zero;
1621 return;
1622 }
1623 }
1624 }
1625
1626 if (limit)
1627 {
1628 // Halve the displacement and call this function again. Will continue
1629 // recursively until the displacement is small enough.
1630
1631 displacement *= 0.5;
1632
1633 limitDisplacement(vit, displacement, callCount);
1634 }
1635}
1636
1637
1638Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
1639(
1640 Foam::point pA,
1641 Foam::point pB
1642) const
1643{
1644 pointIndexHit pAhit;
1645 label pAsurfaceHit = -1;
1646
1647 const scalar searchDist = 5.0*targetCellSize(pA);
1648
1649 geometryToConformTo_.findSurfaceNearest
1650 (
1651 pA,
1652 searchDist,
1653 pAhit,
1654 pAsurfaceHit
1655 );
1656
1657 if (!pAhit.hit())
1658 {
1660 }
1661
1662 vectorField norm(1);
1663
1664 allGeometry_[pAsurfaceHit].getNormal
1665 (
1666 List<pointIndexHit>(1, pAhit),
1667 norm
1668 );
1669
1670 const vector nA = norm[0];
1671
1672 pointIndexHit pBhit;
1673 label pBsurfaceHit = -1;
1674
1675 geometryToConformTo_.findSurfaceNearest
1676 (
1677 pB,
1678 searchDist,
1679 pBhit,
1680 pBsurfaceHit
1681 );
1682
1683 if (!pBhit.hit())
1684 {
1686 }
1687
1688 allGeometry_[pBsurfaceHit].getNormal
1689 (
1690 List<pointIndexHit>(1, pBhit),
1691 norm
1692 );
1693
1694 const vector nB = norm[0];
1695
1696 return vectorTools::cosPhi(nA, nB);
1697}
1698
1699
1700bool Foam::conformalVoronoiMesh::nearSurfacePoint
1701(
1702 pointIndexHitAndFeature& pHit
1703) const
1704{
1705 const Foam::point& pt = pHit.first().hitPoint();
1706
1707 pointIndexHit closePoint;
1708 const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
1709
1710 if
1711 (
1712 closeToSurfacePt
1713 && (
1714 magSqr(pt - closePoint.hitPoint())
1715 > sqr(pointPairDistance(pt))
1716 )
1717 )
1718 {
1719 const scalar cosAngle =
1720 angleBetweenSurfacePoints(pt, closePoint.hitPoint());
1721
1722 // TODO: make this tolerance run-time selectable?
1723 if (cosAngle < searchAngleOppositeSurface)
1724 {
1725 pointIndexHit pCloseHit;
1726 label pCloseSurfaceHit = -1;
1727
1728 const scalar searchDist = targetCellSize(closePoint.hitPoint());
1729
1730 geometryToConformTo_.findSurfaceNearest
1731 (
1732 closePoint.hitPoint(),
1733 searchDist,
1734 pCloseHit,
1735 pCloseSurfaceHit
1736 );
1737
1738 vectorField norm(1);
1739
1740 allGeometry_[pCloseSurfaceHit].getNormal
1741 (
1742 List<pointIndexHit>(1, pCloseHit),
1743 norm
1744 );
1745
1746 const vector& nA = norm[0];
1747
1748 pointIndexHit oppositeHit;
1749 label oppositeSurfaceHit = -1;
1750
1751 geometryToConformTo_.findSurfaceNearestIntersection
1752 (
1753 closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
1754 closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
1755 oppositeHit,
1756 oppositeSurfaceHit
1757 );
1758
1759 if (oppositeHit.hit())
1760 {
1761 // Replace point
1762 pHit.first() = oppositeHit;
1763 pHit.second() = oppositeSurfaceHit;
1764
1765 return !closeToSurfacePt;
1766 }
1767 }
1768 }
1769
1770 return closeToSurfacePt;
1771}
1772
1773
1774bool Foam::conformalVoronoiMesh::appendToSurfacePtTree
1775(
1776 const Foam::point& pt
1777) const
1778{
1779 label startIndex = existingSurfacePtLocations_.size();
1780
1781 existingSurfacePtLocations_.append(pt);
1782
1783 label endIndex = existingSurfacePtLocations_.size();
1784
1785 return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
1786}
1787
1788
1789bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
1790(
1791 const Foam::point& pt
1792) const
1793{
1794 label startIndex = existingEdgeLocations_.size();
1795
1796 existingEdgeLocations_.append(pt);
1797
1798 label endIndex = existingEdgeLocations_.size();
1799
1800 return edgeLocationTreePtr_().insert(startIndex, endIndex);
1801}
1802
1803
1805Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
1806(
1807 const Foam::point& pt
1808) const
1809{
1810 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1811
1812 labelList elems
1813 = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
1814
1815 DynamicList<pointIndexHit> dynPointHit;
1816
1817 forAll(elems, elemI)
1818 {
1819 label index = elems[elemI];
1820
1821 const Foam::point& pointi
1822 = edgeLocationTreePtr_().shapes().shapePoints()[index];
1823
1824 pointIndexHit nearHit(true, pointi, index);
1825
1826 dynPointHit.append(nearHit);
1827 }
1828
1829 return dynPointHit;
1830}
1831
1832
1833bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1834(
1835 const Foam::point& pt
1836) const
1837{
1838 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1839
1840 pointIndexHit info
1841 = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1842
1843 return info.hit();
1844}
1845
1846
1847bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
1848(
1849 const Foam::point& pt,
1850 pointIndexHit& info
1851) const
1852{
1853 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1854
1855 info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1856
1857 return info.hit();
1858}
1859
1860
1861bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1862(
1863 const Foam::point& pt
1864) const
1865{
1866 pointIndexHit info;
1867
1868 pointIsNearSurfaceLocation(pt, info);
1869
1870 return info.hit();
1871}
1872
1873
1874bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
1875(
1876 const Foam::point& pt,
1877 pointIndexHit& info
1878) const
1879{
1880 const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
1881
1882 info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
1883
1884 return info.hit();
1885}
1886
1887
1888bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
1889(
1890 const pointIndexHit& pHit,
1891 pointIndexHit& nearestEdgeHit
1892) const
1893{
1894 const Foam::point& pt = pHit.hitPoint();
1895
1896 const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
1897
1898 bool closeToFeatureEdge =
1899 pointIsNearFeatureEdgeLocation(pt, nearestEdgeHit);
1900
1901 if (closeToFeatureEdge)
1902 {
1903 List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
1904
1905 forAll(nearHits, elemI)
1906 {
1907 pointIndexHit& info = nearHits[elemI];
1908
1909 // Check if the edge location that the new edge location is near to
1910 // "might" be on a different edge. If so, add it anyway.
1911 pointIndexHit edgeHit;
1912 label featureHit = -1;
1913
1914 geometryToConformTo_.findEdgeNearest
1915 (
1916 pt,
1917 exclusionRangeSqr,
1918 edgeHit,
1919 featureHit
1920 );
1921
1922 const extendedFeatureEdgeMesh& eMesh
1923 = geometryToConformTo_.features()[featureHit];
1924
1925 const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
1926
1927 const vector lineBetweenPoints = pt - info.hitPoint();
1928
1929 const scalar cosAngle
1930 = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
1931
1932 // Allow the point to be added if it is almost at right angles to
1933 // the other point. Also check it is not the same point.
1934 // Info<< cosAngle<< " "
1935 // << radToDeg(acos(cosAngle)) << " "
1936 // << searchConeAngle << " "
1937 // << radToDeg(acos(searchConeAngle)) << endl;
1938
1939 if
1940 (
1941 mag(cosAngle) < searchConeAngle
1942 && (mag(lineBetweenPoints) > pointPairDistance(pt))
1943 )
1944 {
1945 //pt = edgeHit.hitPoint();
1946 //pHit.setPoint(pt);
1947 closeToFeatureEdge = false;
1948 }
1949 else
1950 {
1951 closeToFeatureEdge = true;
1952 break;
1953 }
1954 }
1955 }
1956
1957 return closeToFeatureEdge;
1958}
1959
1960
1961void Foam::conformalVoronoiMesh::buildEdgeLocationTree
1962(
1963 const DynamicList<Foam::point>& existingEdgeLocations
1964) const
1965{
1966 treeBoundBox overallBb
1967 (
1968 geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
1969 );
1970
1971 overallBb.min() -= Foam::point::uniform(ROOTVSMALL);
1972 overallBb.max() += Foam::point::uniform(ROOTVSMALL);
1973
1974 edgeLocationTreePtr_.reset
1975 (
1977 (
1978 dynamicTreeDataPoint(existingEdgeLocations),
1979 overallBb, // overall search domain
1980 10, // max levels, n/a
1981 20.0, // maximum ratio of cubes v.s. cells
1982 100.0 // max. duplicity; n/a since no bounding boxes.
1983 )
1984 );
1985}
1986
1987
1988void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
1989(
1990 const DynamicList<Foam::point>& existingSurfacePtLocations
1991) const
1992{
1993 treeBoundBox overallBb
1994 (
1995 geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
1996 );
1997
1998 overallBb.min() -= Foam::point::uniform(ROOTVSMALL);
1999 overallBb.max() += Foam::point::uniform(ROOTVSMALL);
2000
2001 surfacePtLocationTreePtr_.reset
2002 (
2004 (
2005 dynamicTreeDataPoint(existingSurfacePtLocations),
2006 overallBb, // overall search domain
2007 10, // max levels, n/a
2008 20.0, // maximum ratio of cubes v.s. cells
2009 100.0 // max. duplicity; n/a since no bounding boxes.
2010 )
2011 );
2012}
2013
2014
2015void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
2016(
2017 const Foam::point& vit,
2018 const pointIndexHitAndFeatureDynList& surfaceIntersections,
2019 scalar surfacePtReplaceDistCoeffSqr,
2020 scalar edgeSearchDistCoeffSqr,
2021 pointIndexHitAndFeatureDynList& surfaceHits,
2022 pointIndexHitAndFeatureDynList& featureEdgeHits,
2023 DynamicList<label>& surfaceToTreeShape,
2024 DynamicList<label>& edgeToTreeShape,
2025 Map<scalar>& surfacePtToEdgePtDist,
2026 bool firstPass
2027) const
2028{
2029 const scalar cellSize = targetCellSize(vit);
2030 const scalar cellSizeSqr = sqr(cellSize);
2031
2032 forAll(surfaceIntersections, sI)
2033 {
2034 pointIndexHitAndFeature surfHitI = surfaceIntersections[sI];
2035
2036 bool keepSurfacePoint = true;
2037
2038 if (!surfHitI.first().hit())
2039 {
2040 continue;
2041 }
2042
2043 const Foam::point& surfPt = surfHitI.first().hitPoint();
2044
2045 bool isNearFeaturePt = nearFeaturePt(surfPt);
2046
2047 bool isNearFeatureEdge = surfacePtNearFeatureEdge(surfPt);
2048
2049 bool isNearSurfacePt = nearSurfacePoint(surfHitI);
2050
2051 if (isNearFeaturePt || isNearSurfacePt || isNearFeatureEdge)
2052 {
2053 keepSurfacePoint = false;
2054 }
2055
2056 List<List<pointIndexHit>> edHitsByFeature;
2057
2058 labelList featuresHit;
2059
2060 const scalar searchRadiusSqr = edgeSearchDistCoeffSqr*cellSizeSqr;
2061
2062 geometryToConformTo_.findAllNearestEdges
2063 (
2064 surfPt,
2065 searchRadiusSqr,
2066 edHitsByFeature,
2067 featuresHit
2068 );
2069
2070 forAll(edHitsByFeature, i)
2071 {
2072 const label featureHit = featuresHit[i];
2073
2074 List<pointIndexHit>& edHits = edHitsByFeature[i];
2075
2076 forAll(edHits, eHitI)
2077 {
2078 pointIndexHit& edHit = edHits[eHitI];
2079
2080 if (edHit.hit())
2081 {
2082 const Foam::point& edPt = edHit.hitPoint();
2083
2084 if
2085 (
2087 && !decomposition().positionOnThisProcessor(edPt)
2088 )
2089 {
2090 // Do not insert
2091 continue;
2092 }
2093
2094 if (!nearFeaturePt(edPt))
2095 {
2096 if
2097 (
2098 magSqr(edPt - surfPt)
2099 < surfacePtReplaceDistCoeffSqr*cellSizeSqr
2100 )
2101 {
2102 // If the point is within a given distance of a
2103 // feature edge, give control to edge control points
2104 // instead, this will prevent "pits" forming.
2105
2106 // Allow if different surfaces
2107
2108
2109 keepSurfacePoint = false;
2110 }
2111
2112 pointIndexHit nearestEdgeHit;
2113
2114 if
2115 (
2116// !pointIsNearFeatureEdgeLocation
2117// (
2118// edPt,
2119// nearestEdgeHit
2120// )
2121 !nearFeatureEdgeLocation(edHit, nearestEdgeHit)
2122 )
2123 {
2124 appendToEdgeLocationTree(edPt);
2125
2126 edgeToTreeShape.append
2127 (
2128 existingEdgeLocations_.size() - 1
2129 );
2130
2131 // Do not place edge control points too close to a
2132 // feature point or existing edge control points
2133 featureEdgeHits.append
2134 (
2135 pointIndexHitAndFeature(edHit, featureHit)
2136 );
2137
2138// Info<< "Add " << existingEdgeLocations_.size() - 1
2139// << " " << magSqr(edPt - surfPt) << endl;
2140
2141 surfacePtToEdgePtDist.insert
2142 (
2143 existingEdgeLocations_.size() - 1,
2144 magSqr(edPt - surfPt)
2145 );
2146 }
2147 else if (firstPass)
2148 {
2149 label hitIndex = nearestEdgeHit.index();
2150
2151// Info<< "Close to " << nearestEdgeHit << endl;
2152
2153 if
2154 (
2155 magSqr(edPt - surfPt)
2156 < surfacePtToEdgePtDist[hitIndex]
2157 )
2158 {
2159 featureEdgeHits[hitIndex] =
2160 pointIndexHitAndFeature(edHit, featureHit);
2161
2162 existingEdgeLocations_[hitIndex] =
2163 edHit.hitPoint();
2164 surfacePtToEdgePtDist[hitIndex] =
2165 magSqr(edPt - surfPt);
2166
2167 // Change edge location in featureEdgeHits
2168 // remove index from edge tree
2169 // reinsert new point into tree
2170 edgeLocationTreePtr_().remove(hitIndex);
2171 edgeLocationTreePtr_().insert
2172 (
2173 hitIndex,
2174 hitIndex + 1
2175 );
2176 }
2177 }
2178 }
2179 }
2180 }
2181 }
2182
2183 if (keepSurfacePoint)
2184 {
2185 surfaceHits.append(surfHitI);
2186 appendToSurfacePtTree(surfPt);
2187 surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
2188
2189// addedPoints.write(surfPt);
2190 }
2191 else
2192 {
2193// removedPoints.write(surfPt);
2194 }
2195 }
2196}
2197
2198
2199void Foam::conformalVoronoiMesh::storeSurfaceConformation()
2200{
2201 Info<< nl << "Storing surface conformation" << endl;
2202
2203 surfaceConformationVertices_.clear();
2204
2205 // Use a temporary dynamic list to speed up insertion.
2206 DynamicList<Vb> tempSurfaceVertices(number_of_vertices()/10);
2207
2208 for
2209 (
2210 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
2211 vit != finite_vertices_end();
2212 vit++
2213 )
2214 {
2215 // Store points that are not referred, part of a pair, but not feature
2216 // points
2217 if
2218 (
2219 !vit->referred()
2220 && vit->boundaryPoint()
2221 && !vit->featurePoint()
2222 && !vit->constrained()
2223 )
2224 {
2225 tempSurfaceVertices.append
2226 (
2227 Vb
2228 (
2229 vit->point(),
2230 vit->index(),
2231 vit->type(),
2233 )
2234 );
2235 }
2236 }
2237
2238 tempSurfaceVertices.shrink();
2239
2240 surfaceConformationVertices_.transfer(tempSurfaceVertices);
2241
2242 Info<< " Stored "
2243 << returnReduce
2244 (
2245 label(surfaceConformationVertices_.size()),
2246 sumOp<label>()
2247 )
2248 << " vertices" << nl << endl;
2249}
2250
2251
2252void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
2253{
2254 Info<< nl << "Reinserting stored surface conformation" << endl;
2255
2256 Map<label> oldToNewIndices =
2257 insertPointPairs(surfaceConformationVertices_, true, true);
2258
2259 ptPairs_.reIndex(oldToNewIndices);
2260
2261 bitSet selectedElems(surfaceConformationVertices_.size(), true);
2262
2263 forAll(surfaceConformationVertices_, vI)
2264 {
2265 Vb& v = surfaceConformationVertices_[vI];
2266 label& vIndex = v.index();
2267
2268 const auto iter = oldToNewIndices.cfind(vIndex);
2269
2270 if (iter.found())
2271 {
2272 const label newIndex = *iter;
2273
2274 if (newIndex != -1)
2275 {
2276 vIndex = newIndex;
2277 }
2278 else
2279 {
2280 selectedElems.unset(vI);
2281 }
2282 }
2283 }
2284
2285 inplaceSubset<bitSet, List<Vb>>
2286 (
2287 selectedElems,
2288 surfaceConformationVertices_
2289 );
2290}
2291
2292
2293// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
bool found
label n
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
Foam::label & index()
void resetCellCount()
Set the cell count to zero.
void sync(const boundBox &bb)
Refer vertices so that the processor interfaces are consistent.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
Definition: HashTableI.H:141
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
void clear()
Clear all entries from table.
Definition: HashTable.C:678
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label index() const noexcept
Return the hit index.
bool hit() const noexcept
Is there a hit?
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
UPstream::rangeType allProcs() const noexcept
Range of ranks indices associated with PstreamBuffers.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
T & first()
Return the first element of the list.
Definition: UListI.H:202
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
const backgroundMeshDecomposition & decomposition() const
Return the backgroundMeshDecomposition.
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
Non-pointer based hierarchical recursive searching. Storage is dynamic, so elements can be deleted.
Holds (reference to) pointField. Encapsulation of data needed for octree searches....
const vectorField & edgeDirections() const
Return the edgeDirection vectors.
A reference point and direction.
Definition: plane.H:110
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
int myProcNo() const noexcept
Return processor number.
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar pi(M_PI)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
Collection of functions for testing relationships between two vectors.
Definition: vectorTools.H:50
T cosPhi(const Vector< T > &a, const Vector< T > &b, const T &tolerance=SMALL)
Calculate angle between a and b in radians.
Definition: vectorTools.H:107
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
complex limit(const complex &, const complex &)
Definition: complexI.H:263
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
HashSet< labelPair, Foam::Hash< labelPair > > labelPairHashSet
A HashSet for a labelPair. The hashing is based on labelPair (FixedList) and is thus non-commutative.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333