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.point() <<
" 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 pHit.hitPoint().dist(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.point() <<
" 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 pHit.hitPoint().dist(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.point() <<
" 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 pHit.hitPoint().dist(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.point() <<
" 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.point()
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.point()
742 <<
" from edge[" << edgeI <<
"] " << e1
743 <<
" hits inside of surf2 face[" << surf2Facei <<
"]" 744 <<
" - discarded" <<
endl;
753 Pout<<
"hit-type[6] " << pHit.point()
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 (pHit.hitPoint().dist(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;
1102 DynamicList<edge> allCutEdges(surf1.
nEdges()/20);
1103 DynamicList<point> allCutPoints(surf1.
nPoints()/20);
1107 List<DynamicList<label>> edgeCuts1(query1.
surface().
nEdges());
1120 transfer(edgeCuts1, surf1EdgeCuts_);
1128 Pout<<
"Cutting surf2 edges" <<
endl;
1132 List<DynamicList<label>> edgeCuts2(query2.
surface().
nEdges());
1146 joinDisconnected(allCutEdges);
1149 transfer(edgeCuts2, surf2EdgeCuts_);
1155 Pout<<
"surfaceIntersection : Intersection generated:" 1157 <<
" points:" << cutPoints_.
size() <<
endl 1158 <<
" edges :" << cutEdges_.
size() <<
endl;
1160 Pout<<
"surfaceIntersection : Writing intersection to intEdges.obj" 1163 OBJstream(
"intEdges.obj").
write(cutEdges_, cutPoints_);
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();
1223 const triSurface& surf1 = query1.
surface();
1230 Pout<<
"Cutting surf1 edges" <<
endl;
1233 DynamicList<edge> allCutEdges;
1234 DynamicList<point> allCutPoints;
1237 List<DynamicList<label>> edgeCuts1(query1.
surface().
nEdges());
1251 joinDisconnected(allCutEdges);
1254 transfer(edgeCuts1, surf1EdgeCuts_);
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" 1279 OBJstream(
"intEdges.obj").
write(cutEdges_, cutPoints_);
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;
1330 List<DynamicList<label>> edgeCuts1(surf1.
nEdges());
1332 forAll(intersections1, edgeI)
1334 const List<pointIndexHit>& intersections = intersections1[edgeI];
1340 const label cutPointId = allCutPoints.
size();
1342 allCutPoints.
append(pHit.hitPoint());
1343 edgeCuts1[edgeI].append(cutPointId);
1358 transfer(edgeCuts1, surf1EdgeCuts_);
1367 Pout<<
"Storing surf2 intersections" <<
endl;
1372 List<DynamicList<label>> edgeCuts2(surf2.
nEdges());
1374 forAll(intersections2, edgeI)
1376 const List<pointIndexHit>& intersections = intersections2[edgeI];
1382 const label cutPointId = allCutPoints.
size();
1384 allCutPoints.
append(pHit.hitPoint());
1385 edgeCuts2[edgeI].append(cutPointId);
1400 transfer(edgeCuts2, surf2EdgeCuts_);
1411 Pout<<
"surfaceIntersection : Intersection generated:" 1413 <<
" points:" << cutPoints_.
size() <<
endl 1414 <<
" edges :" << cutEdges_.
size() <<
endl;
1416 Pout<<
"surfaceIntersection : Writing intersection to intEdges.obj" 1419 OBJstream(
"intEdges.obj").
write(cutEdges_, cutPoints_);
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 edge&
e = cutEdges_[edgei];
1505 e[0] = pointMap[
e[0]];
1506 e[1] = pointMap[
e[1]];
1509 forAll(surf1EdgeCuts_, edgei)
1514 forAll(surf2EdgeCuts_, edgei)
1529 label nUniqEdges = 0;
1530 labelList edgeNumbering(cutEdges_.size(), -1);
1534 const edge&
e = cutEdges_[edgeI];
1538 if (
e[0] !=
e[1] && uniqEdges.insert(
e))
1540 edgeNumbering[edgeI] = nUniqEdges;
1541 if (nUniqEdges != edgeI)
1543 cutEdges_[nUniqEdges] =
e;
1545 cutEdges_[nUniqEdges].sort();
1559 cutEdges_.setSize(nUniqEdges);
surfaceIntersection()
Construct null.
label nPoints() const
Number of points supporting patch faces.
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
A list of keyword definitions, which are a keyword followed by a number of values (eg...
void append(const T &val)
Append an element at the end of the list.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
const labelListList & surf2EdgeCuts() const
List of cut points on edges of surface2.
#define forAll(list, i)
Loop across all elements in list.
static scalar setPlanarTol(const scalar t)
Set the planar tolerance, returning the previous value.
Helper class to search on triSurface.
const edgeList & cutEdges() const
The list of created edges.
static scalar planarTol()
Return planar tolerance.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
List< edge > edgeList
A List of edges.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Self-intersection, region-wise only.
void clear()
Clear all entries from table.
labelledTri face_type
The face type (same as the underlying PrimitivePatch)
void mergeEdges()
Merge duplicate edges.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const labelListList & edgeFaces() const
Return edge-face addressing.
const labelListList & edgeCuts(const bool isFirstSurf) const
Access either surf1EdgeCuts (isFirstSurface = true) or.
void append(const T &val)
Copy append an element to the end of this list.
void inplaceUniqueSort(ListType &input)
Inplace sorting and removal of duplicates.
label inplaceMergePoints(PointList &points, const scalar mergeTol, const bool verbose, labelList &pointToUnique)
Inplace merge points, preserving the original point order. All points closer/equal mergeTol are to be...
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
HashSet< edge, Hash< edge > > edgeHashSet
A HashSet with edge for its key.
label nEdges() const
Number of edges in patch.
int debug
Static debugging option.
label facePoint(const int facei, const block &block, const label i, const label j)
Pair< label > labelPair
A pair of labels.
defineTypeNameAndDebug(combustionModel, 0)
Geometric merging of points. See below.
const pointField & cutPoints() const
The list of cut points.
const triSurface & surface() const
Return reference to the surface.
const labelPairLookup & facePairToEdgeId() const
Lookup of pairs of faces to created edges.
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
vector point
Point is a vector.
static constexpr direction size() noexcept
The number of elements in the VectorSpace = Ncmpts.
const labelListList & surf1EdgeCuts() const
List of cut points on edges of surface1.
#define WarningInFunction
Report a warning using Foam::Warning.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
void mergePoints(const scalar mergeDist)
Geometric merge points (points within mergeDist) prior to.
messageStream Info
Information stream (stdout output on master, null elsewhere)
intersectionType
Surface intersection types for classify, doCutEdges.
None = invalid (for input only)
List< label > labelList
A List of labels.
Triangulated surface description with patch information.
EnumType getOrDefault(const word &key, const dictionary &dict, const EnumType deflt, const bool failsafe=false) const
Find the key in the dictionary and return the corresponding enumeration element based on its name...
Holder of intersections of edges of a surface with another surface. Optionally shuffles around points...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAllConstIters(mixture.phases(), phase)
static const Enum< intersectionType > selfIntersectionNames
The user-selectable self-intersection enumeration names.