conformalVoronoiMeshFeaturePoints.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"
27 #include "vectorTools.H"
28 #include "triangle.H"
29 #include "tetrahedron.H"
30 #include "ConstCirculator.H"
31 #include "DelaunayMeshTools.H"
32 #include "OBJstream.H"
33 
34 using namespace Foam::vectorTools;
35 
36 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
37 
39 (
40  const extendedFeatureEdgeMesh& feMesh,
41  const pointIndexHit& edHit,
42  DynamicList<Vb>& pts
43 ) const
44 {
45  if (foamyHexMeshControls().circulateEdges())
46  {
47  createEdgePointGroupByCirculating(feMesh, edHit, pts);
48  }
49  else
50  {
51  label edgeI = edHit.index();
52 
54  feMesh.getEdgeStatus(edgeI);
55 
56  switch (edStatus)
57  {
59  {
60  createExternalEdgePointGroup(feMesh, edHit, pts);
61  break;
62  }
64  {
65  createInternalEdgePointGroup(feMesh, edHit, pts);
66  break;
67  }
69  {
70  createFlatEdgePointGroup(feMesh, edHit, pts);
71  break;
72  }
74  {
75  createOpenEdgePointGroup(feMesh, edHit, pts);
76  break;
77  }
79  {
80  createMultipleEdgePointGroup(feMesh, edHit, pts);
81  break;
82  }
84  {
85  break;
86  }
87  }
88  }
89 }
90 
91 
92 bool Foam::conformalVoronoiMesh::meshableRegion
93 (
94  const plane::side side,
96 ) const
97 {
98  switch (volType)
99  {
101  {
102  return (side == plane::FLIP) ? true : false;
103  }
105  {
106  return (side == plane::NORMAL) ? true : false;
107  }
109  {
110  return true;
111  }
113  {
114  return false;
115  }
116  default:
117  {
118  return false;
119  }
120  }
121 }
122 
123 
124 bool Foam::conformalVoronoiMesh::regionIsInside
125 (
127  const vector& normalA,
129  const vector& normalB,
130  const vector& masterPtVec
131 ) const
132 {
133  plane::side sideA
134  (
135  ((masterPtVec & normalA) <= 0) ? plane::FLIP : plane::NORMAL
136  );
137 
138  plane::side sideB
139  (
140  ((masterPtVec & normalB) <= 0) ? plane::FLIP : plane::NORMAL
141  );
142 
143  const bool meshableRegionA = meshableRegion(sideA, volTypeA);
144  const bool meshableRegionB = meshableRegion(sideB, volTypeB);
145 
146  if (meshableRegionA == meshableRegionB)
147  {
148  return meshableRegionA;
149  }
150  else
151  {
153  << ""
154  << endl;
155 
156  return false;
157  }
158 }
159 
160 
161 void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
162 (
163  const extendedFeatureEdgeMesh& feMesh,
164  const pointIndexHit& edHit,
165  DynamicList<Vb>& pts
166 ) const
167 {
168  typedef Foam::indexedVertexEnum::vertexType vertexType;
169  typedef extendedFeatureEdgeMesh::sideVolumeType sideVolumeType;
170 
171  const Foam::point& edgePt = edHit.hitPoint();
172  const label edgeI = edHit.index();
173 
174  scalar ppDist = pointPairDistance(edgePt);
175 
176  const vectorField& feNormals = feMesh.normals();
177  const vector& edDir = feMesh.edgeDirections()[edgeI];
178  const labelList& edNormalIs = feMesh.edgeNormals()[edgeI];
179  const labelList& feNormalDirections = feMesh.normalDirections()[edgeI];
180 
181  const List<sideVolumeType>& normalVolumeTypes = feMesh.normalVolumeTypes();
182 
183  ConstCirculator<labelList> circ(edNormalIs);
184  ConstCirculator<labelList> circNormalDirs(feNormalDirections);
185 
186  Map<Foam::point> masterPoints;
187  Map<vertexType> masterPointsTypes;
188  Map<plane> masterPointReflectionsPrev;
189  Map<plane> masterPointReflectionsNext;
190 
191 // Info<< "Edge = " << edHit << ", edDir = " << edDir << endl;
192 // Info<< " edNorms = " << edNormalIs.size() << endl;
193 
194  bool addedMasterPreviously = false;
195  label initialRegion = -1;
196 
197  if (circ.size()) do
198  {
199  const sideVolumeType volType = normalVolumeTypes[circ()];
200  const sideVolumeType nextVolType = normalVolumeTypes[circ.next()];
201 
202  const vector& normal = feNormals[circ()];
203  const vector& nextNormal = feNormals[circ.next()];
204 
205  vector normalDir = (normal ^ edDir);
206  normalDir *= circNormalDirs()/mag(normalDir);
207 
208  vector nextNormalDir = (nextNormal ^ edDir);
209  nextNormalDir *= circNormalDirs.next()/mag(nextNormalDir);
210 
211 // Info<< " " << circ() << " " << circ.next() << nl
212 // << " " << circNormalDirs() << " " << circNormalDirs.next()
213 // << nl
214 // << " normals = " << normal << " " << nextNormal << nl
215 // << " normalDirs = " << normalDir << " " << nextNormalDir << nl
216 // << " cross = " << (normalDir ^ nextNormalDir) << nl
217 // << " "
218 // << extendedFeatureEdgeMesh::sideVolumeTypeNames_[volType] << " "
219 // << extendedFeatureEdgeMesh::sideVolumeTypeNames_[nextVolType]
220 // << endl;
221 
222  // Calculate master point
223  vector masterPtVec(normalDir + nextNormalDir);
224  masterPtVec /= mag(masterPtVec) + small;
225 
226  if
227  (
228  ((normalDir ^ nextNormalDir) & edDir) < small
229  || mag(masterPtVec) < small
230  )
231  {
232 // Info<< " IGNORE REGION" << endl;
233  addedMasterPreviously = false;
234 
235  if
236  (
237  circ.size() == 2
238  && mag((normal & nextNormal) - 1) < small
239  )
240  {
241  const vector n = 0.5*(normal + nextNormal);
242 
243  const vector s = ppDist*(edDir ^ n);
244 
245  if (volType == extendedFeatureEdgeMesh::BOTH)
246  {
247  createBafflePointPair(ppDist, edgePt + s, n, true, pts);
248  createBafflePointPair(ppDist, edgePt - s, n, true, pts);
249  }
250  else
251  {
253  << "Failed to insert flat/open edge as volType is "
255  [
256  volType
257  ]
258  << endl;
259  }
260 
261  break;
262  }
263 
264  continue;
265  }
266 
267  const Foam::point masterPt = edgePt + ppDist*masterPtVec;
268 
269  // Check that region is inside or outside
270  const bool inside =
271  regionIsInside
272  (
273  volType,
274  normal,
275  nextVolType,
276  nextNormal,
277  masterPtVec
278  );
279 
280  // Specialize for size = 1 && baffle
281  if (mag((normalDir & nextNormalDir) - 1) < small)
282  {
283  if (inside)
284  {
285 // Info<< "Specialize for size 1 and baffle" << endl;
286 
287  vector s = ppDist*(edDir ^ normal);
288 
289  plane facePlane(edgePt, normal);
290 
291  Foam::point pt1 = edgePt + s + ppDist*normal;
292  Foam::point pt2 = edgePt - s + ppDist*normal;
293 
294  Foam::point pt3 = facePlane.mirror(pt1);
295  Foam::point pt4 = facePlane.mirror(pt2);
296 
297  pts.append(Vb(pt1, Vb::vtInternalFeatureEdge));
298  pts.append(Vb(pt2, Vb::vtInternalFeatureEdge));
299  pts.append(Vb(pt3, Vb::vtInternalFeatureEdge));
300  pts.append(Vb(pt4, Vb::vtInternalFeatureEdge));
301 
302  break;
303  }
304  else
305  {
307  << "Faces are parallel but master point is not inside"
308  << endl;
309  }
310  }
311 
312  if (!addedMasterPreviously)
313  {
314  if (initialRegion == -1)
315  {
316  initialRegion = circ.nRotations();
317  }
318 
319  addedMasterPreviously = true;
320 
321  masterPoints.insert(circ.nRotations(), masterPt);
322  masterPointsTypes.insert
323  (
324  circ.nRotations(),
325  inside
328  );
329 
330  masterPointReflectionsPrev.insert
331  (
332  circ.nRotations(),
333  plane(edgePt, normal)
334  );
335 
336  masterPointReflectionsNext.insert
337  (
338  circ.nRotations(),
339  plane(edgePt, nextNormal)
340  );
341  }
342  else if (addedMasterPreviously)
343  {
344  addedMasterPreviously = true;
345 
346  masterPointReflectionsNext.erase(circ.nRotations() - 1);
347 
348  // Shift the master point to be normal to the plane between it and
349  // the previous master point
350  // Should be the intersection of the normal and the plane with the
351  // new master point in it.
352 
353  plane p(masterPoints[circ.nRotations() - 1], normalDir);
354  plane::ray r(edgePt, masterPt - edgePt);
355 
356  scalar cutPoint = p.normalIntersect(r);
357 
358  masterPoints.insert
359  (
360  circ.nRotations(),
361  edgePt + cutPoint*(masterPt - edgePt)
362  );
363 
364  masterPointsTypes.insert
365  (
366  circ.nRotations(),
367  inside
370  );
371 
372  masterPointReflectionsNext.insert
373  (
374  circ.nRotations(),
375  plane(edgePt, nextNormal)
376  );
377  }
378 
379  if
380  (
381  masterPoints.size() > 1
382  && inside
383  && circ.nRotations() == circ.size() - 1
384  )
385  {
386  if (initialRegion == 0)
387  {
388  plane p(masterPoints[initialRegion], nextNormalDir);
389  plane::ray r(edgePt, masterPt - edgePt);
390 
391  scalar cutPoint = p.normalIntersect(r);
392 
393  masterPoints[circ.nRotations()] =
394  edgePt + cutPoint*(masterPt - edgePt);
395 
396  // Remove the first reflection plane if we are no longer
397  // circulating
398 
399  masterPointReflectionsPrev.erase(initialRegion);
400  masterPointReflectionsNext.erase(circ.nRotations());
401  }
402  else
403  {
404 
405  }
406  }
407  }
408  while
409  (
410  circ.circulate(CirculatorBase::CLOCKWISE),
411  circNormalDirs.circulate(CirculatorBase::CLOCKWISE)
412  );
413 
414 
415  forAllConstIter(Map<Foam::point>, masterPoints, iter)
416  {
417  const Foam::point& pt = masterPoints[iter.key()];
418  const vertexType ptType = masterPointsTypes[iter.key()];
419 
420 // Info<< " Adding Master " << iter.key() << " " << pt << " "
421 // << indexedVertexEnum::vertexTypeNames_[ptType] << endl;
422 
423  pts.append(Vb(pt, ptType));
424 
425  const vertexType reflectedPtType =
426  (
427  ptType == Vb::vtInternalFeatureEdge
430  );
431 
432  if (masterPointReflectionsPrev.found(iter.key()))
433  {
434  const Foam::point reflectedPt =
435  masterPointReflectionsPrev[iter.key()].mirror(pt);
436 
437 // Info<< " Adding Prev " << reflectedPt << " "
438 // << indexedVertexEnum::vertexTypeNames_[reflectedPtType]
439 // << endl;
440 
441  pts.append(Vb(reflectedPt, reflectedPtType));
442  }
443 
444  if (masterPointReflectionsNext.found(iter.key()))
445  {
446  const Foam::point reflectedPt =
447  masterPointReflectionsNext[iter.key()].mirror(pt);
448 
449 // Info<< " Adding Next " << reflectedPt << " "
450 // << indexedVertexEnum::vertexTypeNames_[reflectedPtType]
451 // << endl;
452 
453  pts.append(Vb(reflectedPt, reflectedPtType));
454  }
455  }
456 
457 // pts.append(Vb(edgePt, Vb::vtExternalFeatureEdge));
458 }
459 
460 
461 void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
462 (
463  const extendedFeatureEdgeMesh& feMesh,
464  const pointIndexHit& edHit,
465  DynamicList<Vb>& pts
466 ) const
467 {
468  const Foam::point& edgePt = edHit.hitPoint();
469 
470  scalar ppDist = pointPairDistance(edgePt);
471 
472  const vectorField& feNormals = feMesh.normals();
473  const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
474  const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
475  feMesh.normalVolumeTypes();
476 
477  // As this is an external edge, there are two normals by definition
478  const vector& nA = feNormals[edNormalIs[0]];
479  const vector& nB = feNormals[edNormalIs[1]];
480 
482  normalVolumeTypes[edNormalIs[0]];
483 
485  normalVolumeTypes[edNormalIs[1]];
486 
487  if (areParallel(nA, nB))
488  {
489  // The normals are nearly parallel, so this is too sharp a feature to
490  // conform to.
491  return;
492  }
493 
494  // Normalised distance of reference point from edge point
495  vector refVec((nA + nB)/(1 + (nA & nB)));
496 
497  if (magSqr(refVec) > sqr(5.0))
498  {
499  // Limit the size of the conformation
500  ppDist *= 5.0/mag(refVec);
501 
502  // Pout<< nl << "createExternalEdgePointGroup limit "
503  // << "edgePt " << edgePt << nl
504  // << "refVec " << refVec << nl
505  // << "mag(refVec) " << mag(refVec) << nl
506  // << "ppDist " << ppDist << nl
507  // << "nA " << nA << nl
508  // << "nB " << nB << nl
509  // << "(nA & nB) " << (nA & nB) << nl
510  // << endl;
511  }
512 
513  // Convex. So refPt will be inside domain and hence a master point
514  Foam::point refPt = edgePt - ppDist*refVec;
515 
516  // Insert the master point pairing the first slave
517 
518  if (!geometryToConformTo_.inside(refPt))
519  {
520  return;
521  }
522 
523  pts.append
524  (
525  Vb
526  (
527  refPt,
528  vertexCount() + pts.size(),
531  )
532  );
533 
534  // Insert the slave points by reflecting refPt in both faces.
535  // with each slave referring to the master
536 
537  Foam::point reflectedA = refPt + 2*ppDist*nA;
538  pts.append
539  (
540  Vb
541  (
542  reflectedA,
543  vertexCount() + pts.size(),
544  (
548  ),
550  )
551  );
552 
553  Foam::point reflectedB = refPt + 2*ppDist*nB;
554  pts.append
555  (
556  Vb
557  (
558  reflectedB,
559  vertexCount() + pts.size(),
560  (
564  ),
566  )
567  );
568 
569  ptPairs_.addPointPair
570  (
571  pts[pts.size() - 3].index(),
572  pts[pts.size() - 1].index()
573  );
574 
575  ptPairs_.addPointPair
576  (
577  pts[pts.size() - 3].index(),
578  pts[pts.size() - 2].index()
579  );
580 }
581 
582 
583 void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
584 (
585  const extendedFeatureEdgeMesh& feMesh,
586  const pointIndexHit& edHit,
587  DynamicList<Vb>& pts
588 ) const
589 {
590  const Foam::point& edgePt = edHit.hitPoint();
591 
592  scalar ppDist = pointPairDistance(edgePt);
593 
594  const vectorField& feNormals = feMesh.normals();
595  const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
596  const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
597  feMesh.normalVolumeTypes();
598 
599  // As this is an external edge, there are two normals by definition
600  const vector& nA = feNormals[edNormalIs[0]];
601  const vector& nB = feNormals[edNormalIs[1]];
602 
604  normalVolumeTypes[edNormalIs[0]];
605 
606 // const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
607 // normalVolumeTypes[edNormalIs[1]];
608 
609  if (areParallel(nA, nB))
610  {
611  // The normals are nearly parallel, so this is too sharp a feature to
612  // conform to.
613  return;
614  }
615 
616  // Normalised distance of reference point from edge point
617  vector refVec((nA + nB)/(1 + (nA & nB)));
618 
619  if (magSqr(refVec) > sqr(5.0))
620  {
621  // Limit the size of the conformation
622  ppDist *= 5.0/mag(refVec);
623 
624  // Pout<< nl << "createInternalEdgePointGroup limit "
625  // << "edgePt " << edgePt << nl
626  // << "refVec " << refVec << nl
627  // << "mag(refVec) " << mag(refVec) << nl
628  // << "ppDist " << ppDist << nl
629  // << "nA " << nA << nl
630  // << "nB " << nB << nl
631  // << "(nA & nB) " << (nA & nB) << nl
632  // << endl;
633  }
634 
635  // Concave. master and reflected points inside the domain.
636  Foam::point refPt = edgePt - ppDist*refVec;
637 
638  // Generate reflected master to be outside.
639  Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
640 
641  // Reflect reflMasterPt in both faces.
642  Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
643 
644  Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
645 
646  scalar totalAngle =
648 
649  // Number of quadrants the angle should be split into
650  int nQuads = int(totalAngle/foamyHexMeshControls().maxQuadAngle()) + 1;
651 
652  // The number of additional master points needed to obtain the
653  // required number of quadrants.
654  int nAddPoints = min(max(nQuads - 2, 0), 2);
655 
656  // Add number_of_vertices() at insertion of first vertex to all numbers:
657  // Result for nAddPoints 1 when the points are eventually inserted
658  // pt index type
659  // reflectedA 0 2
660  // reflectedB 1 2
661  // reflMasterPt 2 0
662 
663  // Result for nAddPoints 1 when the points are eventually inserted
664  // pt index type
665  // reflectedA 0 3
666  // reflectedB 1 3
667  // refPt 2 3
668  // reflMasterPt 3 0
669 
670  // Result for nAddPoints 2 when the points are eventually inserted
671  // pt index type
672  // reflectedA 0 4
673  // reflectedB 1 4
674  // reflectedAa 2 4
675  // reflectedBb 3 4
676  // reflMasterPt 4 0
677 
678  if
679  (
680  !geometryToConformTo_.inside(reflectedA)
681  || !geometryToConformTo_.inside(reflectedB)
682  )
683  {
684  return;
685  }
686 
687  // Master A is inside.
688  pts.append
689  (
690  Vb
691  (
692  reflectedA,
693  vertexCount() + pts.size(),
696  )
697  );
698 
699  // Master B is inside.
700  pts.append
701  (
702  Vb
703  (
704  reflectedB,
705  vertexCount() + pts.size(),
708  )
709  );
710 
711  // Slave is outside.
712  pts.append
713  (
714  Vb
715  (
716  reflMasterPt,
717  vertexCount() + pts.size(),
718  (
722  ),
724  )
725  );
726 
727  ptPairs_.addPointPair
728  (
729  pts[pts.size() - 2].index(),
730  pts[pts.size() - 1].index()
731  );
732 
733  ptPairs_.addPointPair
734  (
735  pts[pts.size() - 3].index(),
736  pts[pts.size() - 1].index()
737  );
738 
739  if (nAddPoints == 1)
740  {
741  // One additional point is the reflection of the slave point,
742  // i.e. the original reference point
743  pts.append
744  (
745  Vb
746  (
747  refPt,
748  vertexCount() + pts.size(),
751  )
752  );
753  }
754  else if (nAddPoints == 2)
755  {
756  Foam::point reflectedAa = refPt + ppDist*nB;
757  pts.append
758  (
759  Vb
760  (
761  reflectedAa,
762  vertexCount() + pts.size(),
765  )
766  );
767 
768  Foam::point reflectedBb = refPt + ppDist*nA;
769  pts.append
770  (
771  Vb
772  (
773  reflectedBb,
774  vertexCount() + pts.size(),
777  )
778  );
779  }
780 }
781 
782 
783 void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
784 (
785  const extendedFeatureEdgeMesh& feMesh,
786  const pointIndexHit& edHit,
787  DynamicList<Vb>& pts
788 ) const
789 {
790  const Foam::point& edgePt = edHit.hitPoint();
791 
792  const scalar ppDist = pointPairDistance(edgePt);
793 
794  const vectorField& feNormals = feMesh.normals();
795  const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
796  const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
797  feMesh.normalVolumeTypes();
798 
799  // As this is a flat edge, there are two normals by definition
800  const vector& nA = feNormals[edNormalIs[0]];
801  const vector& nB = feNormals[edNormalIs[1]];
802 
803  // Average normal to remove any bias to one side, although as this
804  // is a flat edge, the normals should be essentially the same
805  const vector n = 0.5*(nA + nB);
806 
807  // Direction along the surface to the control point, sense of edge
808  // direction not important, as +s and -s can be used because this
809  // is a flat edge
810  vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
811 
812  if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::OUTSIDE)
813  {
814  createPointPair(ppDist, edgePt + s, -n, true, pts);
815  createPointPair(ppDist, edgePt - s, -n, true, pts);
816  }
817  else if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::BOTH)
818  {
819  createBafflePointPair(ppDist, edgePt + s, n, true, pts);
820  createBafflePointPair(ppDist, edgePt - s, n, true, pts);
821  }
822  else
823  {
824  createPointPair(ppDist, edgePt + s, n, true, pts);
825  createPointPair(ppDist, edgePt - s, n, true, pts);
826  }
827 }
828 
829 
830 void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
831 (
832  const extendedFeatureEdgeMesh& feMesh,
833  const pointIndexHit& edHit,
834  DynamicList<Vb>& pts
835 ) const
836 {
837  // Assume it is a baffle and insert flat edge point pairs
838  const Foam::point& edgePt = edHit.hitPoint();
839 
840  const scalar ppDist = pointPairDistance(edgePt);
841 
842  const vectorField& feNormals = feMesh.normals();
843  const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
844 
845  if (edNormalIs.size() == 1)
846  {
847 // Info<< "Inserting open edge point group around " << edgePt << endl;
848 // Info<< " ppDist = " << ppDist << nl
849 // << " edNormals = " << edNormalIs
850 // << endl;
851 
852  const vector& n = feNormals[edNormalIs[0]];
853 
854  vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
855 
856  plane facePlane(edgePt, n);
857 
858  const label initialPtsSize = pts.size();
859 
860  if
861  (
862  !geometryToConformTo_.inside(edgePt)
863  )
864  {
865  return;
866  }
867 
868  createBafflePointPair(ppDist, edgePt - s, n, true, pts);
869  createBafflePointPair(ppDist, edgePt + s, n, false, pts);
870 
871  for (label ptI = initialPtsSize; ptI < pts.size(); ++ptI)
872  {
873  pts[ptI].type() = Vb::vtInternalFeatureEdge;
874  }
875  }
876  else
877  {
878  Info<< "NOT INSERTING OPEN EDGE POINT GROUP WITH MORE THAN 1 "
879  << "EDGE NORMAL, NOT IMPLEMENTED" << endl;
880  }
881 }
882 
883 
884 void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
885 (
886  const extendedFeatureEdgeMesh& feMesh,
887  const pointIndexHit& edHit,
888  DynamicList<Vb>& pts
889 ) const
890 {
891 // Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
892 
893  const Foam::point& edgePt = edHit.hitPoint();
894 
895  const scalar ppDist = pointPairDistance(edgePt);
896 
897  const vector edDir = feMesh.edgeDirections()[edHit.index()];
898 
899  const vectorField& feNormals = feMesh.normals();
900  const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
901  const labelList& normalDirs = feMesh.normalDirections()[edHit.index()];
902 
903  const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
904  feMesh.normalVolumeTypes();
905 
906  labelList nNormalTypes(4, label(0));
907 
908  forAll(edNormalIs, edgeNormalI)
909  {
911  normalVolumeTypes[edNormalIs[edgeNormalI]];
912 
913  nNormalTypes[sType]++;
914  }
915 
916  if (nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 4)
917  {
918  label masterEdgeNormalIndex = -1;
919 
920  forAll(edNormalIs, edgeNormalI)
921  {
923  normalVolumeTypes[edNormalIs[edgeNormalI]];
924 
925  if (sType == extendedFeatureEdgeMesh::BOTH)
926  {
927  masterEdgeNormalIndex = edgeNormalI;
928  break;
929  }
930  }
931 
932  const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
933 
934  label nDir = normalDirs[masterEdgeNormalIndex];
935 
936  vector normalDir =
937  (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
938  normalDir *= nDir/mag(normalDir);
939 
940  Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
941  Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
942 
943  plane plane3(edgePt, normalDir);
944 
945  Foam::point pt3 = plane3.mirror(pt1);
946  Foam::point pt4 = plane3.mirror(pt2);
947 
948  pts.append
949  (
950  Vb
951  (
952  pt1,
953  vertexCount() + pts.size(),
956  )
957  );
958  pts.append
959  (
960  Vb
961  (
962  pt2,
963  vertexCount() + pts.size(),
966  )
967  );
968 
969  ptPairs_.addPointPair
970  (
971  pts[pts.size() - 2].index(), // external 0 -> slave
972  pts[pts.size() - 1].index()
973  );
974 
975  pts.append
976  (
977  Vb
978  (
979  pt3,
980  vertexCount() + pts.size(),
983  )
984  );
985 
986  ptPairs_.addPointPair
987  (
988  pts[pts.size() - 3].index(), // external 0 -> slave
989  pts[pts.size() - 1].index()
990  );
991 
992  pts.append
993  (
994  Vb
995  (
996  pt4,
997  vertexCount() + pts.size(),
1000  )
1001  );
1002 
1003  ptPairs_.addPointPair
1004  (
1005  pts[pts.size() - 3].index(), // external 0 -> slave
1006  pts[pts.size() - 1].index()
1007  );
1008 
1009  ptPairs_.addPointPair
1010  (
1011  pts[pts.size() - 2].index(), // external 0 -> slave
1012  pts[pts.size() - 1].index()
1013  );
1014  }
1015  else if
1016  (
1017  nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 1
1018  && nNormalTypes[extendedFeatureEdgeMesh::INSIDE] == 2
1019  )
1020  {
1021  label masterEdgeNormalIndex = -1;
1022 
1023  forAll(edNormalIs, edgeNormalI)
1024  {
1026  normalVolumeTypes[edNormalIs[edgeNormalI]];
1027 
1028  if (sType == extendedFeatureEdgeMesh::BOTH)
1029  {
1030  masterEdgeNormalIndex = edgeNormalI;
1031  break;
1032  }
1033  }
1034 
1035  const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
1036 
1037  label nDir = normalDirs[masterEdgeNormalIndex];
1038 
1039  vector normalDir =
1040  (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
1041  normalDir *= nDir/mag(normalDir);
1042 
1043  const label nextNormalI =
1044  (masterEdgeNormalIndex + 1) % edNormalIs.size();
1045  if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
1046  {
1047  normalDir *= -1;
1048  }
1049 
1050  Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
1051  Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
1052 
1053  plane plane3(edgePt, normalDir);
1054 
1055  Foam::point pt3 = plane3.mirror(pt1);
1056  Foam::point pt4 = plane3.mirror(pt2);
1057 
1058  pts.append
1059  (
1060  Vb
1061  (
1062  pt1,
1063  vertexCount() + pts.size(),
1066  )
1067  );
1068  pts.append
1069  (
1070  Vb
1071  (
1072  pt2,
1073  vertexCount() + pts.size(),
1076  )
1077  );
1078 
1079  ptPairs_.addPointPair
1080  (
1081  pts[pts.size() - 2].index(), // external 0 -> slave
1082  pts[pts.size() - 1].index()
1083  );
1084 
1085  pts.append
1086  (
1087  Vb
1088  (
1089  pt3,
1090  vertexCount() + pts.size(),
1093  )
1094  );
1095 
1096  ptPairs_.addPointPair
1097  (
1098  pts[pts.size() - 3].index(), // external 0 -> slave
1099  pts[pts.size() - 1].index()
1100  );
1101 
1102  pts.append
1103  (
1104  Vb
1105  (
1106  pt4,
1107  vertexCount() + pts.size(),
1110  )
1111  );
1112 
1113  ptPairs_.addPointPair
1114  (
1115  pts[pts.size() - 3].index(), // external 0 -> slave
1116  pts[pts.size() - 1].index()
1117  );
1118  }
1119 
1120 
1121 // // As this is a flat edge, there are two normals by definition
1122 // const vector& nA = feNormals[edNormalIs[0]];
1123 // const vector& nB = feNormals[edNormalIs[1]];
1124 //
1125 // // Average normal to remove any bias to one side, although as this
1126 // // is a flat edge, the normals should be essentially the same
1127 // const vector n = 0.5*(nA + nB);
1128 //
1129 // // Direction along the surface to the control point, sense of edge
1130 // // direction not important, as +s and -s can be used because this
1131 // // is a flat edge
1132 // vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
1133 //
1134 // createBafflePointPair(ppDist, edgePt + s, n, true, pts);
1135 // createBafflePointPair(ppDist, edgePt - s, n, true, pts);
1136 }
1137 
1138 
1139 void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
1140 {
1141  Info<< nl
1142  << "Inserting feature points" << endl;
1143 
1144  const label preFeaturePointSize(number_of_vertices());
1145 
1146  if (Pstream::parRun() && distribute)
1147  {
1148  ftPtConformer_.distribute(decomposition());
1149  }
1150 
1151  const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
1152 
1153  // Insert the created points directly as already distributed.
1154  Map<label> oldToNewIndices =
1155  this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices, true);
1156 
1157  ftPtConformer_.reIndexPointPairs(oldToNewIndices);
1158 
1159  label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
1160  reduce(nFeatureVertices, sumOp<label>());
1161 
1162  Info<< " Inserted " << nFeatureVertices << " feature vertices" << endl;
1163 }
1164 
1165 
1166 //Foam::scalar Foam::conformalVoronoiMesh::pyramidVolume
1167 //(
1168 // const Foam::point& apex,
1169 // const Foam::point& a,
1170 // const Foam::point& b,
1171 // const Foam::point& c,
1172 // const bool printInfo
1173 //) const
1174 //{
1175 // triPointRef tri(a, b, c);
1176 //
1177 // tetPointRef tet(tri.a(), tri.b(), tri.c(), apex);
1178 //
1179 // scalar volume = tet.mag();
1180 //
1188 //
1189 // if (printInfo)
1190 // {
1191 // Info<< "Calculating volume of pyramid..." << nl
1192 // << " Apex : " << apex << nl
1193 // << " Point a : " << a << nl
1194 // << " Point b : " << b << nl
1195 // << " Point c : " << c << nl
1196 // << " Center : " << tri.centre() << nl
1197 // << " Volume : " << volume << endl;
1198 // }
1199 //
1200 // return volume;
1201 //}
1202 
1203 
1204 //void Foam::conformalVoronoiMesh::createPyramidMasterPoint
1205 //(
1206 // const Foam::point& apex,
1207 // const vectorField& edgeDirections,
1208 // Foam::point& masterPoint,
1209 // vectorField& norms
1210 //) const
1211 //{
1212 // pointField basePoints(edgeDirections.size() + 1);
1213 //
1214 // forAll(edgeDirections, eI)
1215 // {
1216 // basePoints[eI] = edgeDirections[eI] + apex;
1217 // }
1218 //
1219 // basePoints[edgeDirections.size() + 1] = apex;
1220 //
1221 // face f(identity(edgeDirections.size()));
1222 //
1223 // pyramidPointFaceRef p(f, apex);
1224 //
1225 // const scalar ppDist = pointPairDistance(apex);
1226 //
1227 //
1228 // vector unitDir = f.centre();
1229 // unitDir /= mag(unitDir);
1230 //
1231 // masterPoint = apex + ppDist*unitDir;
1232 //
1233 // norms.setSize(edgeDirections.size());
1234 //
1235 // forAll(norms, nI)
1236 // {
1237 // norms[nI] =
1238 // }
1239 //}
1240 
1241 
1242 //void Foam::conformalVoronoiMesh::createConvexConcaveFeaturePoints
1243 //(
1244 // DynamicList<Foam::point>& pts,
1245 // DynamicList<label>& indices,
1246 // DynamicList<label>& types
1247 //)
1248 //{
1249 // const PtrList<extendedFeatureEdgeMesh>& feMeshes
1250 // (
1251 // geometryToConformTo_.features()
1252 // );
1253 //
1254 // forAll(feMeshes, i)
1255 // {
1256 // const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
1257 //
1258 // for
1259 // (
1260 // label ptI = feMesh.convexStart();
1261 // ptI < feMesh.mixedStart();
1262 // ptI++
1263 // )
1264 // {
1265 // const Foam::point& apex = feMesh.points()[ptI];
1266 //
1267 // if (!positionOnThisProc(apex))
1268 // {
1269 // continue;
1270 // }
1271 //
1272 // const vectorField& featPtEdgeDirections
1273 // = feMesh.featurePointEdgeDirections(ptI);
1274 //
1275 // Foam::point masterPoint;
1276 // vectorField tetNorms;
1277 //
1278 // createPyramidMasterPoint
1279 // (
1280 // apex,
1281 // featPtEdgeDirections,
1282 // masterPoint,
1283 // tetNorms
1284 // );
1285 //
1286 //
1287 //
1288 // // Result when the points are eventually inserted (example n = 4)
1289 // // Add number_of_vertices() at insertion of first vertex to all
1290 // // numbers:
1291 // // pt index type
1292 // // internalPt 0 1
1293 // // externalPt0 1 0
1294 // // externalPt1 2 0
1295 // // externalPt2 3 0
1296 // // externalPt3 4 0
1297 //
1298 // // Result when the points are eventually inserted (example n = 5)
1299 // // Add number_of_vertices() at insertion of first vertex to all
1300 // // numbers:
1301 // // pt index type
1302 // // internalPt0 0 5
1303 // // internalPt1 1 5
1304 // // internalPt2 2 5
1305 // // internalPt3 3 5
1306 // // internalPt4 4 5
1307 // // externalPt 5 4
1308 //
1309 // if (geometryToConformTo_.inside(masterPoint))
1310 // {
1311 //
1312 // }
1313 // else
1314 // {
1315 //
1316 // }
1317 //
1318 // pts.append(masterPoint);
1319 // indices.append(0);
1320 // types.append(1);
1321 //
1322 // label internalPtIndex = -1;
1323 //
1324 // forAll(tetNorms, nI)
1325 // {
1326 // const vector& n = tetNorms[nI];
1327 //
1328 // Foam::point reflectedPoint
1329 // = reflectPoint(featPt, masterPoint, n);
1330 //
1331 // pts.append(reflectedPoint);
1332 // indices.append(0);
1333 // types.append(internalPtIndex--);
1334 // }
1335 // }
1336 // }
1337 //}
#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
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:427
static const Foam::NamedEnum< sideVolumeType, 4 > sideVolumeTypeNames_
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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
void createEdgePointGroup(const extendedFeatureEdgeMesh &feMesh, const pointIndexHit &edHit, DynamicList< Vb > &pts) const
Call the appropriate function to conform to an edge.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:51
CGAL::indexedVertex< K > Vb
sideVolumeType
Normals point to the outside.
List< label > labelList
A List of labels.
Definition: labelList.H:56
T radAngleBetween(const Vector< T > &a, const Vector< T > &b, const T &tolerance=small)
Calculate angle between a and b in radians.
Definition: vectorTools.H:120
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
Definition: VectorSpaceI.H:83
static const char nl
Definition: Ostream.H:265
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:397
side
Side of the plane.
Definition: plane.H:65
#define WarningInFunction
Report a warning using Foam::Warning.
bool areParallel(const Vector< T > &a, const Vector< T > &b, const T &tolerance=small)
Test if a and b are parallel: a^b = 0.
Definition: vectorTools.H:54
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
Map< label > insertPoints(const List< Vb > &vertices, const bool reIndex)
Insert the list of vertices (calls rangeInsertWithInfo)
volScalarField & p