53 { intersectionType::SELF,
"self" },
54 { intersectionType::SELF_REGION,
"region" },
55 { intersectionType::NONE,
"none" },
61 void Foam::surfaceIntersection::setOptions(
const dictionary&
dict)
63 dict.readIfPresent(
"tolerance", tolerance_);
64 dict.readIfPresent(
"allowEdgeHits", allowEdgeHits_);
65 dict.readIfPresent(
"snap", snapToEnd_);
66 dict.readIfPresent(
"warnDegenerate", warnDegenerate_);
70 void Foam::surfaceIntersection::storeIntersection
72 const enum intersectionType cutFrom,
75 const UList<point>& allCutPoints,
76 const label cutPointId,
77 DynamicList<edge>& allCutEdges
86 const label faceA = facesA[facesAI];
93 twoFaces.first() = faceA;
99 twoFaces.second() = faceA;
108 twoFaces.first() = faceA;
109 twoFaces.second() = faceB;
113 twoFaces.first() = faceB;
114 twoFaces.second() = faceA;
126 edge& thisEdge = facePairToEdge_(twoFaces);
127 const label pointCount = thisEdge.count();
132 thisEdge.insert(cutPointId);
136 Pout<<
"intersect faces " << twoFaces
137 <<
" point-1: " << cutPointId <<
" = "
138 << allCutPoints[cutPointId] <<
endl;
142 else if (pointCount == 2)
148 Pout<<
"suppressed double intersection " << twoFaces
154 if (thisEdge.insert(cutPointId))
162 if (edgeToId_.found(thisEdge))
166 if (facePairToEdge_.insert(twoFaces, thisEdge))
170 Pout<<
"reuse edge - faces " << twoFaces <<
" edge#"
171 << edgeToId_[thisEdge] <<
" edge " << thisEdge
172 <<
" = " << thisEdge.line(allCutPoints)
177 else if (thisEdge.mag(allCutPoints) < SMALL)
194 if (warnDegenerate_ > 0)
198 <<
"Degenerate edge between faces " << twoFaces
199 <<
" on 1st/2nd surface with points "
200 << thisEdge.line(allCutPoints)
205 Pout<<
"degenerate edge face-pair " << twoFaces <<
" "
206 << thisEdge[0] <<
" point " << allCutPoints[thisEdge[0]]
211 thisEdge.erase(cutPointId);
216 const label edgeId = allCutEdges.size();
218 if (facePairToEdgeId_.insert(twoFaces, edgeId))
222 edgeToId_.insert(thisEdge, edgeId);
223 allCutEdges.append(thisEdge);
227 Pout<<
"create edge - faces " << twoFaces <<
" edge#"
228 << edgeId <<
" edge " << thisEdge
229 <<
" = " << thisEdge.line(allCutPoints)
237 Info<<
"WARN " << twoFaces
238 <<
" already intersected= " << thisEdge <<
endl;
239 thisEdge.erase(cutPointId);
247 facePairToEdge_.erase(twoFaces);
263 void Foam::surfaceIntersection::classifyHit
265 const triSurface& surf1,
267 const triSurface& surf2,
268 const enum intersectionType cutFrom,
272 DynamicList<point>& allCutPoints,
273 DynamicList<edge>& allCutEdges,
274 List<DynamicList<label>>& surfEdgeCuts
277 const edge& e1 = surf1.edges()[edgeI];
279 const labelList& facesA = surf1.edgeFaces()[edgeI];
282 label surf2Facei = pHit.index();
287 const pointField& surf1Pts = surf1.localPoints();
288 const pointField& surf2Pts = surf2.localPoints();
290 label nearType, nearLabel;
292 f2.nearestPointClassify(pHit.hitPoint(), surf2Pts, nearType, nearLabel);
295 const label edgeEnd =
298 surf1PointTol[e1.start()],
299 surf1PointTol[e1.end()],
312 Pout<<
"hit-type[1] " << pHit.hitPoint() <<
" is surf1:"
313 <<
" end point of edge[" << edgeI <<
"] " << e1
314 <<
"==" << e1.line(surf1Pts)
315 <<
" surf2: vertex " << f2[nearLabel]
316 <<
" coord:" << surf2Pts[f2[nearLabel]]
317 <<
" - suppressed" <<
endl;
323 label cutPointId = -1;
324 const label nearVert = f2[nearLabel];
335 const point& nearPt = surf1Pts[nearVert];
339 mag(pHit.hitPoint() - nearPt)
340 < surf1PointTol[nearVert]
343 cutPointId = allCutPoints.
size();
347 if (snappedEnds_.insert(nearVert, cutPointId))
350 allCutPoints.append(nearPt);
355 cutPointId = snappedEnds_[nearVert];
360 allCutPoints.append(nearPt);
367 Pout<<
"hit-type[2] " << pHit.hitPoint() <<
" is surf1:"
368 <<
" from edge[" << edgeI <<
"] " << e1
369 <<
" surf2: vertex " << f2[nearLabel]
370 <<
" coord:" << surf2Pts[f2[nearLabel]]
372 << (cutPointId >= 0 ?
"snapped" :
"stored") <<
endl;
375 if (cutPointId == -1)
377 cutPointId = allCutPoints.size();
378 allCutPoints.append(pHit.hitPoint());
380 surfEdgeCuts[edgeI].append(cutPointId);
382 const labelList& facesB = surf2.pointFaces()[f2[nearLabel]];
408 const label edge2I = getEdge(surf2, surf2Facei, nearLabel);
409 const edge& e2 = surf2.edges()[edge2I];
410 const label nearVert = e1[edgeEnd];
412 label cutPointId = -1;
418 int handling = (allowEdgeHits_ ? 1 : 0);
432 const edge intersect(edgeI, edge2I);
434 if (e2.found(nearVert))
440 else if (edgeEdgeIntersection_.insert(intersect))
442 const point& nearPt = surf1Pts[nearVert];
446 mag(pHit.hitPoint() - nearPt)
447 < surf1PointTol[nearVert]
450 cutPointId = allCutPoints.
size();
454 if (snappedEnds_.insert(nearVert, cutPointId))
457 allCutPoints.append(nearPt);
462 cutPointId = snappedEnds_[nearVert];
468 allCutPoints.append(nearPt);
480 Pout<<
"hit-type[3] " << pHit.hitPoint() <<
" is surf1:"
481 <<
" end point of edge[" << edgeI <<
"] " << e1
482 <<
"==" << e1.line(surf1Pts)
483 <<
" surf2: edge[" << edge2I <<
"] " << e2
484 <<
" coords:" << e2.line(surf2Pts)
488 ?
"cached" : handling
489 ?
"stored" :
"suppressed"
495 if (cutPointId == -1)
497 cutPointId = allCutPoints.size();
498 allCutPoints.append(pHit.hitPoint());
500 surfEdgeCuts[edgeI].append(cutPointId);
505 const labelList& facesB = surf2.edgeFaces()[edge2I];
536 const label edge2I = getEdge(surf2, surf2Facei, nearLabel);
537 const edge& e2 = surf2.edges()[edge2I];
538 label cutPointId = -1;
563 const edge intersect(edgeI, edge2I);
565 if (edgeEdgeIntersection_.insert(intersect))
570 const label endId = e1[edgepti];
571 const point& nearPt = surf1Pts[endId];
575 mag(pHit.hitPoint() - nearPt)
576 < surf1PointTol[endId]
579 cutPointId = allCutPoints.
size();
583 if (snappedEnds_.insert(endId, cutPointId))
586 allCutPoints.append(nearPt);
591 cutPointId = snappedEnds_[endId];
597 allCutPoints.append(nearPt);
615 Pout<<
"hit-type[4] " << pHit.hitPoint() <<
" is surf1:"
616 <<
" from edge[" << edgeI <<
"] " << e1
617 <<
"==" << e1.line(surf1Pts)
618 <<
" surf2: edge[" << edge2I <<
"] " << e2
619 <<
" coords:" << e2.line(surf2Pts)
623 ?
"cut-point" : handling
624 ?
"stored" :
"suppressed"
631 if (cutPointId == -1)
633 cutPointId = allCutPoints.size();
634 allCutPoints.append(pHit.hitPoint());
636 surfEdgeCuts[edgeI].append(cutPointId);
641 const vector eVec = e1.unitVec(surf1Pts);
643 const labelList& facesB = surf2.edgeFaces()[edge2I];
649 mag((surf2.faceNormals()[facesB[faceBI]] & eVec))
679 const label nearVert = (edgeEnd == 0 ? e1.start() : e1.end());
680 const label otherVert = (edgeEnd == 0 ? e1.end() : e1.start());
682 const point& nearPt = surf1Pts[nearVert];
683 const point& otherPt = surf1Pts[otherVert];
685 const vector eVec = otherPt - nearPt;
687 if ((surf2.faceNormals()[surf2Facei] & eVec) > 0)
693 label cutPointId = allCutPoints.
size();
696 if (snappedEnds_.insert(nearVert, cutPointId))
699 allCutPoints.append(nearPt);
704 cutPointId = snappedEnds_[nearVert];
710 allCutPoints.append(nearPt);
713 surfEdgeCuts[edgeI].append(cutPointId);
717 Pout<<
"hit-type[5] " << pHit.hitPoint()
718 <<
" shifted to " << nearPt
719 <<
" from edge[" << edgeI <<
"] " << e1
720 <<
"==" << e1.line(surf1Pts)
721 <<
" hits surf2 face[" << surf2Facei <<
"]"
723 << (cached ?
"cached" :
"stored") <<
endl;
741 Pout<<
"hit-type[5] " << pHit.hitPoint()
742 <<
" from edge[" << edgeI <<
"] " << e1
743 <<
" hits inside of surf2 face[" << surf2Facei <<
"]"
744 <<
" - discarded" <<
endl;
753 Pout<<
"hit-type[6] " << pHit.hitPoint()
754 <<
" from edge[" << edgeI <<
"] " << e1
755 <<
"==" << e1.line(surf1Pts)
756 <<
" hits surf2 face[" << surf2Facei <<
"]"
757 <<
" - stored" <<
endl;
760 const label cutPointId = allCutPoints.size();
761 allCutPoints.append(pHit.hitPoint());
762 surfEdgeCuts[edgeI].append(cutPointId);
787 void Foam::surfaceIntersection::doCutEdges
789 const triSurface& surf1,
790 const triSurfaceSearch& querySurf2,
791 const enum intersectionType cutFrom,
793 DynamicList<point>& allCutPoints,
794 DynamicList<edge>& allCutEdges,
795 List<DynamicList<label>>& surfEdgeCuts
800 const pointField& surf1Pts = surf1.localPoints();
805 forAll(surf1PointTol, pointi)
807 surf1PointTol[pointi] = tolerance_ * minEdgeLen(surf1, pointi);
810 const indexedOctree<treeDataPrimitivePatch<triSurface>>& searchTree
822 DynamicList<label> maskFaces(32);
826 bitSet maskRegions(32);
828 treeDataTriSurface::findAllIntersectOp
829 allIntersectOp(searchTree, maskFaces);
831 forAll(surf1.edges(), edgeI)
833 const edge&
e = surf1.edges()[edgeI];
834 const vector edgeVec =
e.vec(surf1Pts);
837 const point ptStart =
838 surf1Pts[
e.start()] - 0.5*surf1PointTol[
e.start()]*edgeVec;
840 surf1Pts[
e.end()] + 0.5*surf1PointTol[
e.end()]*edgeVec;
845 for (
auto& facei : surf1.edgeFaces()[edgeI])
847 maskRegions.set(surf1[facei].region());
854 maskFaces = surf1.pointFaces()[
e.start()];
855 maskFaces.append(surf1.pointFaces()[
e.end()]);
871 maskFaces.append(pHit.index());
873 if (maskRegions.test(surf1[pHit.index()].region()))
896 const triSurface& surf2 = querySurf2.surface();
898 forAll(surf1.edges(), edgeI)
900 const edge&
e = surf1.edges()[edgeI];
901 const vector edgeVec =
e.vec(surf1Pts);
904 const scalar tolDim =
mag(tolVec);
908 surf1Pts[
e.start()] - 0.5*surf1PointTol[
e.start()]*edgeVec;
910 surf1Pts[
e.end()] + 0.5*surf1PointTol[
e.end()]*edgeVec;
912 bool doTrack =
false;
938 if (
mag(pHit.hitPoint() - ptEnd) < tolDim)
946 ptStart = pHit.hitPoint() + tolVec;
960 edgeEdgeIntersection_.clear();
961 snappedEnds_.clear();
967 void Foam::surfaceIntersection::joinDisconnected
969 DynamicList<edge>& allCutEdges
980 Pair<Map<labelPairHashSet>> missedFacePoint;
988 const edge&
e = iter.val();
993 const label pointId =
e.maxVertex();
995 missedFacePoint[0](twoFaces[0]).
insert
1000 missedFacePoint[1](twoFaces[1]).
insert
1012 forAll(missedFacePoint, sidei)
1014 const auto& mapping = missedFacePoint[sidei];
1018 const auto& connect = iter.val();
1020 if (connect.size() == 2)
1033 if (
e.count() == 2 && !edgeToId_.found(
e))
1041 label edgeId = allCutEdges.size();
1042 edgeList newEdgesLst = newEdges.sortedToc();
1043 for (
const auto&
e : newEdgesLst)
1047 edgeToId_.insert(
e, edgeId);
1058 allowEdgeHits_(true),
1064 facePairToEdgeId_(0),
1078 allowEdgeHits_(
true),
1098 Pout<<
"Cutting surf1 edges" <<
endl;
1120 transfer(edgeCuts1, surf1EdgeCuts_);
1128 Pout<<
"Cutting surf2 edges" <<
endl;
1146 joinDisconnected(allCutEdges);
1149 transfer(edgeCuts2, surf2EdgeCuts_);
1150 cutEdges_.transfer(allCutEdges);
1151 cutPoints_.transfer(allCutPoints);
1155 Pout<<
"surfaceIntersection : Intersection generated:"
1157 <<
" points:" << cutPoints_.size() <<
endl
1158 <<
" edges :" << cutEdges_.size() <<
endl;
1160 Pout<<
"surfaceIntersection : Writing intersection to intEdges.obj"
1166 Pout<<
"Dumping cut edges of surface1 to surf1EdgeCuts.obj" <<
endl;
1167 OFstream edge1Stream(
"surf1EdgeCuts.obj");
1168 writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
1170 Pout<<
"Dumping cut edges of surface2 to surf2EdgeCuts.obj" <<
endl;
1171 OFstream edge2Stream(
"surf2EdgeCuts.obj");
1172 writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
1176 facePairToEdge_.clear();
1190 allowEdgeHits_(
true),
1195 facePairToEdge_(2*query1.
surface().size()),
1196 facePairToEdgeId_(2*query1.
surface().size()),
1204 "intersectionMethod",
1206 intersectionType::SELF
1209 if (cutFrom == intersectionType::NONE)
1213 Pout<<
"Skipping self-intersection (selected: none)" <<
endl;
1217 facePairToEdge_.clear();
1218 facePairToEdgeId_.clear();
1230 Pout<<
"Cutting surf1 edges" <<
endl;
1251 joinDisconnected(allCutEdges);
1254 transfer(edgeCuts1, surf1EdgeCuts_);
1255 cutEdges_.transfer(allCutEdges);
1256 cutPoints_.transfer(allCutPoints);
1259 if (cutPoints_.empty() && cutEdges_.empty())
1263 Pout<<
"Empty intersection" <<
endl;
1270 Pout<<
"surfaceIntersection : Intersection generated and compressed:"
1272 <<
" points:" << cutPoints_.size() <<
endl
1273 <<
" edges :" << cutEdges_.size() <<
endl;
1276 Pout<<
"surfaceIntersection : Writing intersection to intEdges.obj"
1282 Pout<<
"Dumping cut edges of surface1 to surf1EdgeCuts.obj" <<
endl;
1283 OFstream edge1Stream(
"surf1EdgeCuts.obj");
1284 writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
1288 facePairToEdge_.clear();
1304 allowEdgeHits_(
true),
1309 facePairToEdge_(2*
max(surf1.size(), surf2.size())),
1310 facePairToEdgeId_(2*
max(surf1.size(), surf2.size())),
1325 Pout<<
"Storing surf1 intersections" <<
endl;
1332 forAll(intersections1, edgeI)
1340 const label cutPointId = allCutPoints.size();
1343 edgeCuts1[edgeI].append(cutPointId);
1358 transfer(edgeCuts1, surf1EdgeCuts_);
1367 Pout<<
"Storing surf2 intersections" <<
endl;
1374 forAll(intersections2, edgeI)
1382 const label cutPointId = allCutPoints.size();
1385 edgeCuts2[edgeI].append(cutPointId);
1400 transfer(edgeCuts2, surf2EdgeCuts_);
1405 cutEdges_.transfer(allCutEdges);
1406 cutPoints_.transfer(allCutPoints);
1411 Pout<<
"surfaceIntersection : Intersection generated:"
1413 <<
" points:" << cutPoints_.size() <<
endl
1414 <<
" edges :" << cutEdges_.size() <<
endl;
1416 Pout<<
"surfaceIntersection : Writing intersection to intEdges.obj"
1422 Pout<<
"Dumping cut edges of surface1 to surf1EdgeCuts.obj" <<
endl;
1423 OFstream edge1Stream(
"surf1EdgeCuts.obj");
1424 writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
1426 Pout<<
"Dumping cut edges of surface2 to surf2EdgeCuts.obj" <<
endl;
1427 OFstream edge2Stream(
"surf2EdgeCuts.obj");
1428 writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
1432 facePairToEdge_.clear();
1455 return facePairToEdgeId_;
1461 const bool isFirstSurf
1466 return surf1EdgeCuts_;
1470 return surf2EdgeCuts_;
1477 return surf1EdgeCuts_;
1483 return surf2EdgeCuts_;
1503 cutPoints_.transfer(newPoints);
1507 edge&
e = cutEdges_[edgei];
1509 e[0] = pointMap[
e[0]];
1510 e[1] = pointMap[
e[1]];
1513 forAll(surf1EdgeCuts_, edgei)
1518 forAll(surf2EdgeCuts_, edgei)
1533 label nUniqEdges = 0;
1534 labelList edgeNumbering(cutEdges_.size(), -1);
1538 const edge&
e = cutEdges_[edgeI];
1542 if (
e[0] !=
e[1] && uniqEdges.insert(
e))
1544 edgeNumbering[edgeI] = nUniqEdges;
1545 if (nUniqEdges != edgeI)
1547 cutEdges_[nUniqEdges] =
e;
1549 cutEdges_[nUniqEdges].sort();
1563 cutEdges_.setSize(nUniqEdges);