43 const extendedFeatureEdgeMesh& feMesh,
48 if (foamyHexMeshControls().circulateEdges())
50 createEdgePointGroupByCirculating(feMesh, edHit, pts);
54 label edgeI = edHit.index();
57 feMesh.getEdgeStatus(edgeI);
63 createExternalEdgePointGroup(feMesh, edHit, pts);
68 createInternalEdgePointGroup(feMesh, edHit, pts);
73 createFlatEdgePointGroup(feMesh, edHit, pts);
78 createOpenEdgePointGroup(feMesh, edHit, pts);
83 createMultipleEdgePointGroup(feMesh, edHit, pts);
95 bool Foam::conformalVoronoiMesh::meshableRegion
127 bool Foam::conformalVoronoiMesh::regionIsInside
146 const bool meshableRegionA = meshableRegion(sideA, volTypeA);
147 const bool meshableRegionB = meshableRegion(sideB, volTypeB);
149 if (meshableRegionA == meshableRegionB)
151 return meshableRegionA;
160 void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
162 const extendedFeatureEdgeMesh& feMesh,
171 const label edgeI = edHit.index();
173 scalar ppDist = pointPairDistance(edgePt);
176 const vector& edDir = feMesh.edgeDirections()[edgeI];
177 const labelList& edNormalIs = feMesh.edgeNormals()[edgeI];
178 const labelList& feNormalDirections = feMesh.normalDirections()[edgeI];
180 const List<sideVolumeType>& normalVolumeTypes = feMesh.normalVolumeTypes();
182 ConstCirculator<labelList> circ(edNormalIs);
183 ConstCirculator<labelList> circNormalDirs(feNormalDirections);
185 Map<Foam::point> masterPoints;
186 Map<vertexType> masterPointsTypes;
187 Map<plane> masterPointReflectionsPrev;
188 Map<plane> masterPointReflectionsNext;
193 bool addedMasterPreviously =
false;
194 label initialRegion = -1;
198 const sideVolumeType volType = normalVolumeTypes[circ()];
199 const sideVolumeType nextVolType = normalVolumeTypes[circ.next()];
201 const vector& normal = feNormals[circ()];
202 const vector& nextNormal = feNormals[circ.next()];
204 vector normalDir = (normal ^ edDir);
205 normalDir *= circNormalDirs()/
mag(normalDir);
207 vector nextNormalDir = (nextNormal ^ edDir);
208 nextNormalDir *= circNormalDirs.next()/
mag(nextNormalDir);
226 ((normalDir ^ nextNormalDir) & edDir) < SMALL
227 ||
mag(masterPtVec) < SMALL
231 addedMasterPreviously =
false;
236 &&
mag((normal & nextNormal) - 1) < SMALL
239 const vector n = 0.5*(normal + nextNormal);
241 const vector s = ppDist*(edDir ^
n);
245 createBafflePointPair(ppDist, edgePt +
s,
n,
true, pts);
246 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
251 <<
"Failed to insert flat/open edge as volType is "
265 const Foam::point masterPt = edgePt + ppDist*masterPtVec;
279 if (
mag((normalDir & nextNormalDir) - 1) < SMALL)
285 vector s = ppDist*(edDir ^ normal);
287 plane facePlane(edgePt, normal);
295 pts.append(
Vb(pt1, Vb::vtInternalFeatureEdge));
296 pts.append(
Vb(pt2, Vb::vtInternalFeatureEdge));
297 pts.append(
Vb(pt3, Vb::vtInternalFeatureEdge));
298 pts.append(
Vb(pt4, Vb::vtInternalFeatureEdge));
305 <<
"Faces are parallel but master point is not inside"
310 if (!addedMasterPreviously)
312 if (initialRegion == -1)
314 initialRegion = circ.nRotations();
317 addedMasterPreviously =
true;
319 masterPoints.insert(circ.nRotations(), masterPt);
320 masterPointsTypes.insert
324 ? Vb::vtInternalFeatureEdge
325 : Vb::vtExternalFeatureEdge
328 masterPointReflectionsPrev.insert
331 plane(edgePt, normal)
334 masterPointReflectionsNext.insert
337 plane(edgePt, nextNormal)
340 else if (addedMasterPreviously)
342 addedMasterPreviously =
true;
344 masterPointReflectionsNext.erase(circ.nRotations() - 1);
351 plane
p(masterPoints[circ.nRotations() - 1], normalDir);
352 plane::ray r(edgePt, masterPt - edgePt);
354 scalar cutPoint =
p.normalIntersect(r);
359 edgePt + cutPoint*(masterPt - edgePt)
362 masterPointsTypes.insert
366 ? Vb::vtInternalFeatureEdge
367 : Vb::vtExternalFeatureEdge
370 masterPointReflectionsNext.insert
373 plane(edgePt, nextNormal)
379 masterPoints.size() > 1
381 && circ.nRotations() == circ.size() - 1
384 if (initialRegion == 0)
386 plane
p(masterPoints[initialRegion], nextNormalDir);
387 plane::ray r(edgePt, masterPt - edgePt);
389 scalar cutPoint =
p.normalIntersect(r);
391 masterPoints[circ.nRotations()] =
392 edgePt + cutPoint*(masterPt - edgePt);
397 masterPointReflectionsPrev.erase(initialRegion);
398 masterPointReflectionsNext.erase(circ.nRotations());
416 const vertexType ptType = masterPointsTypes[iter.key()];
421 pts.append(
Vb(pt, ptType));
423 const vertexType reflectedPtType =
425 ptType == Vb::vtInternalFeatureEdge
426 ? Vb::vtExternalFeatureEdge
427 : Vb::vtInternalFeatureEdge
430 if (masterPointReflectionsPrev.found(iter.key()))
433 masterPointReflectionsPrev[iter.key()].mirror(pt);
439 pts.append(
Vb(reflectedPt, reflectedPtType));
442 if (masterPointReflectionsNext.found(iter.key()))
445 masterPointReflectionsNext[iter.key()].mirror(pt);
451 pts.append(
Vb(reflectedPt, reflectedPtType));
459 void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
461 const extendedFeatureEdgeMesh& feMesh,
468 scalar ppDist = pointPairDistance(edgePt);
471 const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
472 const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
473 feMesh.normalVolumeTypes();
476 const vector& nA = feNormals[edNormalIs[0]];
477 const vector& nB = feNormals[edNormalIs[1]];
480 normalVolumeTypes[edNormalIs[0]];
483 normalVolumeTypes[edNormalIs[1]];
493 vector refVec((nA + nB)/(1 + (nA & nB)));
498 ppDist *= 5.0/
mag(refVec);
516 if (!geometryToConformTo_.inside(refPt))
526 vertexCount() + pts.size(),
527 Vb::vtInternalFeatureEdge,
541 vertexCount() + pts.size(),
544 ? Vb::vtInternalFeatureEdge
545 : Vb::vtExternalFeatureEdge
557 vertexCount() + pts.size(),
560 ? Vb::vtInternalFeatureEdge
561 : Vb::vtExternalFeatureEdge
567 ptPairs_.addPointPair
569 pts[pts.size() - 3].index(),
570 pts[pts.size() - 1].index()
573 ptPairs_.addPointPair
575 pts[pts.size() - 3].index(),
576 pts[pts.size() - 2].index()
581 void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
583 const extendedFeatureEdgeMesh& feMesh,
590 scalar ppDist = pointPairDistance(edgePt);
593 const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
594 const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
595 feMesh.normalVolumeTypes();
598 const vector& nA = feNormals[edNormalIs[0]];
599 const vector& nB = feNormals[edNormalIs[1]];
602 normalVolumeTypes[edNormalIs[0]];
615 vector refVec((nA + nB)/(1 + (nA & nB)));
620 ppDist *= 5.0/
mag(refVec);
637 Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
640 Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
642 Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
648 int nQuads = int(totalAngle/foamyHexMeshControls().maxQuadAngle()) + 1;
652 int nAddPoints =
min(
max(nQuads - 2, 0), 2);
678 !geometryToConformTo_.inside(reflectedA)
679 || !geometryToConformTo_.inside(reflectedB)
691 vertexCount() + pts.size(),
692 Vb::vtInternalFeatureEdge,
703 vertexCount() + pts.size(),
704 Vb::vtInternalFeatureEdge,
715 vertexCount() + pts.size(),
718 ? Vb::vtInternalFeatureEdge
719 : Vb::vtExternalFeatureEdge
725 ptPairs_.addPointPair
727 pts[pts.size() - 2].index(),
728 pts[pts.size() - 1].index()
731 ptPairs_.addPointPair
733 pts[pts.size() - 3].index(),
734 pts[pts.size() - 1].index()
746 vertexCount() + pts.size(),
747 Vb::vtInternalFeatureEdge,
752 else if (nAddPoints == 2)
760 vertexCount() + pts.size(),
761 Vb::vtInternalFeatureEdge,
772 vertexCount() + pts.size(),
773 Vb::vtInternalFeatureEdge,
781 void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
783 const extendedFeatureEdgeMesh& feMesh,
790 const scalar ppDist = pointPairDistance(edgePt);
793 const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
794 const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
795 feMesh.normalVolumeTypes();
798 const vector& nA = feNormals[edNormalIs[0]];
799 const vector& nB = feNormals[edNormalIs[1]];
803 const vector n = 0.5*(nA + nB);
808 vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^
n);
812 createPointPair(ppDist, edgePt +
s, -
n,
true, pts);
813 createPointPair(ppDist, edgePt -
s, -
n,
true, pts);
817 createBafflePointPair(ppDist, edgePt +
s,
n,
true, pts);
818 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
822 createPointPair(ppDist, edgePt +
s,
n,
true, pts);
823 createPointPair(ppDist, edgePt -
s,
n,
true, pts);
828 void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
830 const extendedFeatureEdgeMesh& feMesh,
838 const scalar ppDist = pointPairDistance(edgePt);
841 const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
843 if (edNormalIs.size() == 1)
850 const vector&
n = feNormals[edNormalIs[0]];
852 vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^
n);
854 plane facePlane(edgePt,
n);
856 const label initialPtsSize = pts.size();
860 !geometryToConformTo_.inside(edgePt)
866 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
867 createBafflePointPair(ppDist, edgePt +
s,
n,
false, pts);
869 for (label ptI = initialPtsSize; ptI < pts.size(); ++ptI)
871 pts[ptI].type() = Vb::vtInternalFeatureEdge;
876 Info<<
"NOT INSERTING OPEN EDGE POINT GROUP WITH MORE THAN 1 "
877 <<
"EDGE NORMAL, NOT IMPLEMENTED" <<
endl;
882 void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
884 const extendedFeatureEdgeMesh& feMesh,
893 const scalar ppDist = pointPairDistance(edgePt);
895 const vector edDir = feMesh.edgeDirections()[edHit.index()];
898 const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
899 const labelList& normalDirs = feMesh.normalDirections()[edHit.index()];
901 const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
902 feMesh.normalVolumeTypes();
906 forAll(edNormalIs, edgeNormalI)
909 normalVolumeTypes[edNormalIs[edgeNormalI]];
911 nNormalTypes[sType]++;
916 label masterEdgeNormalIndex = -1;
918 forAll(edNormalIs, edgeNormalI)
921 normalVolumeTypes[edNormalIs[edgeNormalI]];
925 masterEdgeNormalIndex = edgeNormalI;
930 const vector&
n = feNormals[edNormalIs[masterEdgeNormalIndex]];
932 label nDir = normalDirs[masterEdgeNormalIndex];
935 (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
936 normalDir *= nDir/
mag(normalDir);
941 plane plane3(edgePt, normalDir);
951 vertexCount() + pts.size(),
952 Vb::vtInternalSurface,
961 vertexCount() + pts.size(),
962 Vb::vtInternalSurface,
967 ptPairs_.addPointPair
969 pts[pts.size() - 2].index(),
970 pts[pts.size() - 1].index()
978 vertexCount() + pts.size(),
979 Vb::vtInternalSurface,
984 ptPairs_.addPointPair
986 pts[pts.size() - 3].index(),
987 pts[pts.size() - 1].index()
995 vertexCount() + pts.size(),
996 Vb::vtInternalSurface,
1001 ptPairs_.addPointPair
1003 pts[pts.size() - 3].index(),
1004 pts[pts.size() - 1].index()
1007 ptPairs_.addPointPair
1009 pts[pts.size() - 2].index(),
1010 pts[pts.size() - 1].index()
1019 label masterEdgeNormalIndex = -1;
1021 forAll(edNormalIs, edgeNormalI)
1024 normalVolumeTypes[edNormalIs[edgeNormalI]];
1028 masterEdgeNormalIndex = edgeNormalI;
1033 const vector&
n = feNormals[edNormalIs[masterEdgeNormalIndex]];
1035 label nDir = normalDirs[masterEdgeNormalIndex];
1038 (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
1039 normalDir *= nDir/
mag(normalDir);
1041 const label nextNormalI =
1042 (masterEdgeNormalIndex + 1) % edNormalIs.size();
1043 if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
1048 Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*
n;
1049 Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*
n;
1051 plane plane3(edgePt, normalDir);
1061 vertexCount() + pts.size(),
1062 Vb::vtInternalSurface,
1071 vertexCount() + pts.size(),
1072 Vb::vtInternalSurface,
1077 ptPairs_.addPointPair
1079 pts[pts.size() - 2].index(),
1080 pts[pts.size() - 1].index()
1088 vertexCount() + pts.size(),
1089 Vb::vtExternalSurface,
1094 ptPairs_.addPointPair
1096 pts[pts.size() - 3].index(),
1097 pts[pts.size() - 1].index()
1105 vertexCount() + pts.size(),
1106 Vb::vtExternalSurface,
1111 ptPairs_.addPointPair
1113 pts[pts.size() - 3].index(),
1114 pts[pts.size() - 1].index()
1137 void Foam::conformalVoronoiMesh::insertFeaturePoints(
bool distribute)
1140 <<
"Inserting feature points" <<
endl;
1142 const label preFeaturePointSize(number_of_vertices());
1146 ftPtConformer_.distribute(decomposition());
1149 const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
1152 Map<label> oldToNewIndices =
1155 ftPtConformer_.reIndexPointPairs(oldToNewIndices);
1157 label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
1158 reduce(nFeatureVertices, sumOp<label>());
1160 Info<<
" Inserted " << nFeatureVertices <<
" feature vertices" <<
endl;