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