36 bool Foam::featurePointConformer::createSpecialisedFeaturePoint
38 const extendedFeatureEdgeMesh& feMesh,
40 const pointFeatureEdgesTypes& pFEdgesTypes,
41 const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
62 if (debug)
Info<<
"nExternal == 2 && nInternal == 1" <<
endl;
69 && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
75 label nVert = foamyHexMesh_.number_of_vertices();
77 const label initialNumOfPoints = pts.size();
79 const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
85 label concaveEdgeI = -1;
95 concaveEdgeI = pEds[i];
99 convexEdgesI[nConvex++] = pEds[i];
104 <<
"Edge " << eS <<
" is flat" 110 <<
"Edge " << eS <<
" not concave/convex" 115 const vector& concaveEdgePlaneANormal =
116 normals[edgeNormals[concaveEdgeI][0]];
118 const vector& concaveEdgePlaneBNormal =
119 normals[edgeNormals[concaveEdgeI][1]];
125 featPt + ppDist*concaveEdgePlaneANormal,
126 concaveEdgePlaneANormal
131 featPt + ppDist*concaveEdgePlaneBNormal,
132 concaveEdgePlaneBNormal
135 const vector& concaveEdgeDir = feMesh.edgeDirection
143 featPt + ppDist*concaveEdgeDir;
148 plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
150 const Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
158 planeA = plane(featPt, concaveEdgePlaneANormal);
160 planeB = plane(featPt, concaveEdgePlaneBNormal);
163 concaveEdgeExternalPt
164 - 2.0*planeA.distance(concaveEdgeExternalPt)
165 *concaveEdgePlaneANormal;
172 foamyHexMesh_.vertexCount() + pts.
size(),
178 const label internalPtAIndex(pts.last().index());
181 concaveEdgeExternalPt
182 - 2.0*planeB.distance(concaveEdgeExternalPt)
183 *concaveEdgePlaneBNormal;
190 foamyHexMesh_.vertexCount() + pts.size(),
196 const label internalPtBIndex(pts.last().index());
206 const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
207 const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
208 const labelList& convexEdgeBNormals = edgeNormals[convexEdgesI[1]];
210 forAll(concaveEdgeNormals, edgeNormalI)
212 bool convexEdgeA =
false;
213 bool convexEdgeB =
false;
215 forAll(convexEdgeANormals, edgeAnormalI)
217 const vector& concaveNormal
218 = normals[concaveEdgeNormals[edgeNormalI]];
219 const vector& convexNormal
220 = normals[convexEdgeANormals[edgeAnormalI]];
224 Info<<
"Angle between vectors = " 230 if (
areParallel(concaveNormal, convexNormal, tolParallel))
236 forAll(convexEdgeBNormals, edgeBnormalI)
238 const vector& concaveNormal
239 = normals[concaveEdgeNormals[edgeNormalI]];
240 const vector& convexNormal
241 = normals[convexEdgeBNormals[edgeBnormalI]];
245 Info<<
"Angle between vectors = " 251 if (
areParallel(concaveNormal, convexNormal, tolParallel))
257 if ((convexEdgeA && convexEdgeB) || (!convexEdgeA && !convexEdgeB))
260 <<
"Both or neither of the convex edges share the concave " 262 <<
" convexEdgeA = " << convexEdgeA
263 <<
" convexEdgeB = " << convexEdgeB
267 for (
label i = 0; i < 2; ++i)
278 forAll(convexEdgeANormals, edgeAnormalI)
280 const vector& concaveNormal
281 = normals[concaveEdgeNormals[edgeNormalI]];
282 const vector& convexNormal
283 = normals[convexEdgeANormals[edgeAnormalI]];
287 !
areParallel(concaveNormal, convexNormal, tolParallel)
290 convexEdgePlaneCNormal = convexNormal;
292 plane planeC(featPt, convexEdgePlaneCNormal);
296 + 2.0*planeC.distance(internalPtA)
297 *convexEdgePlaneCNormal;
304 foamyHexMesh_.vertexCount() + pts.size(),
310 ftPtPairs_.addPointPair
321 forAll(convexEdgeBNormals, edgeBnormalI)
323 const vector& concaveNormal
324 = normals[concaveEdgeNormals[edgeNormalI]];
325 const vector& convexNormal
326 = normals[convexEdgeBNormals[edgeBnormalI]];
330 !
areParallel(concaveNormal, convexNormal, tolParallel)
333 convexEdgePlaneDNormal = convexNormal;
335 plane planned(featPt, convexEdgePlaneDNormal);
339 + 2.0*planned.distance(internalPtB)
340 *convexEdgePlaneDNormal;
347 foamyHexMesh_.vertexCount() + pts.
size(),
353 ftPtPairs_.addPointPair
367 concaveEdgeExternalPt,
368 foamyHexMesh_.vertexCount() + pts.
size(),
374 ftPtPairs_.addPointPair
380 ftPtPairs_.addPointPair
386 const label concaveEdgeExternalPtIndex(pts.last().index());
394 if (totalAngle > foamyHexMeshControls_.maxQuadAngle())
400 vector convexEdgesPlaneNormal =
401 0.5*(convexEdgePlaneCNormal + convexEdgePlaneDNormal);
403 plane planeM(featPt, convexEdgesPlaneNormal);
428 concaveEdgeExternalPt
430 + 2.0*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
437 foamyHexMesh_.vertexCount() + pts.
size(),
443 const label internalPtFIndex(pts.last().index());
445 ftPtPairs_.addPointPair
447 concaveEdgeExternalPtIndex,
453 + 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal;
460 foamyHexMesh_.vertexCount() + pts.
size(),
466 ftPtPairs_.addPointPair
475 for (
label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI)
477 Info<<
"Point " << ptI <<
" : ";
493 Info<<
"nExternal == 1 && nInternal == 2" <<
endl;
501 && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
507 label nVert = foamyHexMesh_.number_of_vertices();
509 const label initialNumOfPoints = pts.size();
511 const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
517 label convexEdgeI = -1;
527 convexEdgeI = pEds[i];
531 concaveEdgesI[nConcave++] = pEds[i];
536 <<
"Edge " << eS <<
" is flat" 542 <<
"Edge " << eS <<
" not concave/convex" 547 const vector& convexEdgePlaneANormal =
548 normals[edgeNormals[convexEdgeI][0]];
550 const vector& convexEdgePlaneBNormal =
551 normals[edgeNormals[convexEdgeI][1]];
557 featPt - ppDist*convexEdgePlaneANormal,
558 convexEdgePlaneANormal
563 featPt - ppDist*convexEdgePlaneBNormal,
564 convexEdgePlaneBNormal
567 const vector& convexEdgeDir = feMesh.edgeDirection
575 featPt + ppDist*convexEdgeDir;
580 plane planeF(convexEdgeLocalFeatPt, convexEdgeDir);
582 const Foam::point convexEdgeExternalPt = planeF.planePlaneIntersect
590 planeA = plane(featPt, convexEdgePlaneANormal);
592 planeB = plane(featPt, convexEdgePlaneBNormal);
596 + 2.0*planeA.distance(convexEdgeExternalPt)
597 *convexEdgePlaneANormal;
604 foamyHexMesh_.vertexCount() + pts.
size(),
610 const label internalPtAIndex(pts.last().index());
614 + 2.0*planeB.distance(convexEdgeExternalPt)
615 *convexEdgePlaneBNormal;
622 foamyHexMesh_.vertexCount() + pts.size(),
628 const label internalPtBIndex(pts.last().index());
638 const labelList& convexEdgeNormals = edgeNormals[convexEdgeI];
639 const labelList& concaveEdgeANormals = edgeNormals[concaveEdgesI[0]];
640 const labelList& concaveEdgeBNormals = edgeNormals[concaveEdgesI[1]];
642 forAll(convexEdgeNormals, edgeNormalI)
644 bool concaveEdgeA =
false;
645 bool concaveEdgeB =
false;
647 forAll(concaveEdgeANormals, edgeAnormalI)
649 const vector& convexNormal
650 = normals[convexEdgeNormals[edgeNormalI]];
651 const vector& concaveNormal
652 = normals[concaveEdgeANormals[edgeAnormalI]];
656 Info<<
"Angle between vectors = " 662 if (
areParallel(convexNormal, concaveNormal, tolParallel))
668 forAll(concaveEdgeBNormals, edgeBnormalI)
670 const vector& convexNormal
671 = normals[convexEdgeNormals[edgeNormalI]];
672 const vector& concaveNormal
673 = normals[concaveEdgeBNormals[edgeBnormalI]];
677 Info<<
"Angle between vectors = " 683 if (
areParallel(convexNormal, concaveNormal, tolParallel))
691 (concaveEdgeA && concaveEdgeB)
692 || (!concaveEdgeA && !concaveEdgeB)
696 <<
"Both or neither of the concave edges share the convex " 698 <<
" concaveEdgeA = " << concaveEdgeA
699 <<
" concaveEdgeB = " << concaveEdgeB
703 for (
label i = 0; i < 2; ++i)
714 forAll(concaveEdgeANormals, edgeAnormalI)
716 const vector& convexNormal
717 = normals[convexEdgeNormals[edgeNormalI]];
718 const vector& concaveNormal
719 = normals[concaveEdgeANormals[edgeAnormalI]];
723 !
areParallel(convexNormal, concaveNormal, tolParallel)
726 concaveEdgePlaneCNormal = concaveNormal;
728 plane planeC(featPt, concaveEdgePlaneCNormal);
732 - 2.0*planeC.distance(internalPtA)
733 *concaveEdgePlaneCNormal;
740 foamyHexMesh_.vertexCount() + pts.size(),
746 ftPtPairs_.addPointPair
757 forAll(concaveEdgeBNormals, edgeBnormalI)
759 const vector& convexNormal
760 = normals[convexEdgeNormals[edgeNormalI]];
761 const vector& concaveNormal
762 = normals[concaveEdgeBNormals[edgeBnormalI]];
766 !
areParallel(convexNormal, concaveNormal, tolParallel)
769 concaveEdgePlaneDNormal = concaveNormal;
771 plane planned(featPt, concaveEdgePlaneDNormal);
775 - 2.0*planned.distance(internalPtB)
776 *concaveEdgePlaneDNormal;
783 foamyHexMesh_.vertexCount() + pts.
size(),
789 ftPtPairs_.addPointPair
803 convexEdgeExternalPt,
804 foamyHexMesh_.vertexCount() + pts.
size(),
810 ftPtPairs_.addPointPair
816 ftPtPairs_.addPointPair
828 if (totalAngle > foamyHexMeshControls_.maxQuadAngle())
834 vector convexEdgesPlaneNormal =
835 0.5*(concaveEdgePlaneCNormal + concaveEdgePlaneDNormal);
837 plane planeM(featPt, convexEdgesPlaneNormal);
864 + 2.0*(convexEdgeLocalFeatPt - convexEdgeExternalPt);
871 foamyHexMesh_.vertexCount() + pts.
size(),
877 ftPtPairs_.addPointPair
879 pts[pts.size() - 2].index(),
885 - 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal;
892 foamyHexMesh_.vertexCount() + pts.
size(),
898 ftPtPairs_.addPointPair
900 pts[pts.size() - 2].index(),
907 for (
label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI)
909 Info<<
"Point " << ptI <<
" " List< labelList > labelListList
A List of labelList.
static const Foam::NamedEnum< vertexType, 15 > vertexTypeNames_
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
pointFromPoint topoint(const Point &P)
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
List< label > labelList
A List of labels.
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
static bool & parRun()
Is this a parallel run?
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
Field< vector > vectorField
Specialisation of Field<T> for vector.