50 void Foam::intersectedSurface::writeOBJ
59 const point& pt = points[pointI];
61 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
65 const edge&
e = edges[edgeI];
67 os <<
"l " << e.start()+1 <<
' ' << e.end()+1 <<
nl;
83 const point& pt = points[pointI];
85 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
89 const edge& e = edges[faceEdges[i]];
91 os <<
"l " << e.start()+1 <<
' ' << e.end()+1 <<
nl;
97 void Foam::intersectedSurface::writeLocalOBJ
102 const fileName& fName
113 const edge& e = edges[faceEdges[i]];
119 if (pointMap[pointI] == -1)
121 const point& pt = points[pointI];
123 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
125 pointMap[pointI] = maxVertI++;
132 const edge& e = edges[faceEdges[i]];
134 os <<
"l " << pointMap[e.start()]+1 <<
' ' << pointMap[e.end()]+1
150 const point& pt = points[pointI];
152 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
159 os <<
' ' << f[fp]+1;
166 void Foam::intersectedSurface::printVisit
170 const Map<label>& visited
176 label edgeI = edgeLabels[i];
178 const edge& e = edges[edgeI];
180 label stat = visited[edgeI];
182 if (stat == UNVISITED)
184 Pout<<
" edge:" << edgeI <<
" verts:" << e
185 <<
" unvisited" <<
nl;
187 else if (stat == STARTTOEND)
189 Pout<<
" edge:" << edgeI <<
" from " << e[0]
190 <<
" to " << e[1] <<
nl;
192 else if (stat == ENDTOSTART)
194 Pout<<
" edge:" << edgeI <<
" from " << e[1]
195 <<
" to " << e[0] <<
nl;
199 Pout<<
" edge:" << edgeI <<
" both " << e
210 bool Foam::intersectedSurface::sameEdgeOrder
212 const labelledTri& fA,
213 const labelledTri& fB
223 label vA1 = fA[(fpA + 1) % 3];
224 label vAMin1 = fA[fpA ? fpA-1 : 2];
227 label vB1 = fB[(fpB + 1) % 3];
228 label vBMin1 = fB[fpB ? fpB-1 : 2];
230 if (vA1 == vB1 || vAMin1 == vBMin1)
234 else if (vA1 == vBMin1 || vAMin1 == vB1)
242 <<
"Triangle:" << fA <<
" and triangle:" << fB
243 <<
" share a point but not an edge"
250 <<
"Triangle:" << fA <<
" and triangle:" << fB
251 <<
" do not share an edge"
258 void Foam::intersectedSurface::incCount
267 if (iter == visited.end())
269 visited.insert(key, offset);
282 Foam::intersectedSurface::calcPointEdgeAddressing
284 const edgeSurface& eSurf,
289 const edgeList& edges = eSurf.edges();
291 const labelList& fEdges = eSurf.faceEdges()[faceI];
293 Map<DynamicList<label> > facePointEdges(4*fEdges.size());
297 label edgeI = fEdges[i];
299 const edge& e = edges[edgeI];
302 Map<DynamicList<label> >::iterator iter =
303 facePointEdges.find(e.start());
305 if (iter == facePointEdges.end())
307 DynamicList<label> oneEdge;
308 oneEdge.append(edgeI);
309 facePointEdges.insert(e.start(), oneEdge);
313 iter().append(edgeI);
317 Map<DynamicList<label> >::iterator iter2 =
318 facePointEdges.find(e.end());
320 if (iter2 == facePointEdges.end())
322 DynamicList<label> oneEdge;
323 oneEdge.append(edgeI);
324 facePointEdges.insert(e.end(), oneEdge);
328 iter2().append(edgeI);
335 Map<DynamicList<label> >::iterator iter = facePointEdges.begin();
336 iter != facePointEdges.end();
347 "intersectedSurface::calcPointEdgeAddressing"
348 "(const edgeSurface&, const label)"
349 ) <<
"Point:" << iter.key() <<
" used by too few edges:"
357 Pout<<
"calcPointEdgeAddressing: face consisting of edges:" <<
endl;
360 label edgeI = fEdges[i];
361 const edge& e = edges[edgeI];
362 Pout<<
" " << edgeI <<
' ' << e << points[e.start()]
363 << points[e.end()] <<
endl;
366 Pout<<
" Constructed point-edge adressing:" <<
endl;
369 Map<DynamicList<label> >::iterator iter = facePointEdges.begin();
370 iter != facePointEdges.end();
374 Pout<<
" vertex " << iter.key() <<
" is connected to edges "
380 return facePointEdges;
388 Foam::label Foam::intersectedSurface::nextEdge
390 const edgeSurface& eSurf,
391 const Map<label>& visited,
394 const Map<DynamicList<label> >& facePointEdges,
395 const label prevEdgeI,
396 const label prevVertI
400 const edgeList& edges = eSurf.edges();
401 const labelList& fEdges = eSurf.faceEdges()[faceI];
405 const DynamicList<label>& connectedEdges = facePointEdges[prevVertI];
407 if (connectedEdges.size() <= 1)
411 Pout<<
"Writing face:" << faceI <<
" to face.obj" <<
endl;
412 OFstream str(
"face.obj");
413 writeOBJ(points, edges, fEdges, str);
415 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
416 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
421 "intersectedSurface::nextEdge(const pointField&, const edgeList&"
422 ", const vector&, Map<DynamicList<label> >, const label"
424 ) <<
"Problem: prevVertI:" << prevVertI <<
" on edge " << prevEdgeI
425 <<
" has less than 2 connected edges."
431 if (connectedEdges.size() == 2)
434 if (connectedEdges[0] == prevEdgeI)
438 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
439 <<
" to edge:" << edges[connectedEdges[1]]
440 <<
" chosen from candidates:" << connectedEdges <<
endl;
442 return connectedEdges[1];
448 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
449 <<
" to edge:" << edges[connectedEdges[0]]
450 <<
" chosen from candidates:" << connectedEdges <<
endl;
452 return connectedEdges[0];
459 const edge& prevE = edges[prevEdgeI];
462 vector e0 = n ^ (points[prevE.otherVertex(prevVertI)] - points[prevVertI]);
463 e0 /=
mag(e0) + VSMALL;
468 if (
mag(
mag(e1) - 1) > SMALL)
471 Pout<<
"Writing face:" << faceI <<
" to face.obj" <<
endl;
472 OFstream str(
"face.obj");
473 writeOBJ(points, edges, fEdges, str);
475 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
476 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
480 <<
"Unnormalized normal e1:" << e1
481 <<
" formed from cross product of e0:" << e0 <<
" n:" << n
490 scalar maxAngle = -GREAT;
493 forAll(connectedEdges, connI)
495 label edgeI = connectedEdges[connI];
497 if (edgeI != prevEdgeI)
499 label stat = visited[edgeI];
501 const edge& e = edges[edgeI];
519 n ^ (points[e.otherVertex(prevVertI)] - points[prevVertI]);
521 vec /=
mag(vec) + VSMALL;
525 if (angle > maxAngle)
539 Pout<<
"Writing face:" << faceI <<
" to face.obj" <<
endl;
540 OFstream str(
"face.obj");
541 writeOBJ(points, edges, fEdges, str);
543 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
544 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
549 "intersectedSurface::nextEdge(const pointField&, const edgeList&"
550 ", const Map<label>&, const vector&"
551 ", const Map<DynamicList<label> >&"
552 ", const label, const label"
553 ) <<
"Trying to step from edge " << edges[prevEdgeI]
554 <<
", vertex " << prevVertI
555 <<
" but cannot find 'unvisited' edges among candidates:"
562 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
563 <<
" to edge:" << maxEdgeI <<
" angle:" << edges[maxEdgeI]
564 <<
" with angle:" << maxAngle
578 const edgeSurface& eSurf,
581 const Map<DynamicList<label> >& facePointEdges,
583 const label startEdgeI,
584 const label startVertI,
590 const edgeList& edges = eSurf.edges();
593 face
f(eSurf.faceEdges()[faceI].size());
597 label vertI = startVertI;
598 label edgeI = startEdgeI;
602 const edge& e = edges[edgeI];
606 Pout<<
"Now at:" << endl
607 <<
" edge:" << edgeI <<
" vertices:" << e
608 <<
" positions:" << points[e.start()] <<
' ' << points[e.end()]
609 <<
" vertex:" << vertI <<
endl;
615 visited[edgeI] |= STARTTOEND;
619 visited[edgeI] |= ENDTOSTART;
627 vertI = e.otherVertex(vertI);
629 if (vertI == startVertI)
653 void Foam::intersectedSurface::findNearestVisited
655 const edgeSurface& eSurf,
657 const Map<DynamicList<label> >& facePointEdges,
658 const Map<label>& pointVisited,
660 const label excludePointI,
671 label pointI = iter.key();
673 if (pointI != excludePointI)
675 label nVisits = iter();
677 if (nVisits == 2*facePointEdges[pointI].size())
680 scalar
dist =
mag(eSurf.points()[pointI] - pt);
693 const labelList& fEdges = eSurf.faceEdges()[faceI];
696 <<
"Dumping face edges to faceEdges.obj" <<
endl;
698 writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges,
"faceEdges.obj");
701 <<
"No fully visited edge found for pt " << pt
717 const triSurface& surf,
719 const Map<DynamicList<label> >& facePointEdges,
720 const Map<label>& visited,
726 Pout<<
"Writing face:" << faceI <<
" to face.obj" <<
endl;
727 OFstream str(
"face.obj");
728 writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[faceI], str);
734 Map<label> pointVisited(2*facePointEdges.size());
738 OFstream str0(
"visitedNone.obj");
739 OFstream str1(
"visitedOnce.obj");
740 OFstream str2(
"visitedTwice.obj");
741 forAll(eSurf.points(), pointI)
743 const point& pt = eSurf.points()[pointI];
745 str0 <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
746 str1 <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
747 str2 <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
753 label edgeI = iter.key();
755 const edge& e = eSurf.edges()[edgeI];
759 if (stat == STARTTOEND || stat == ENDTOSTART)
761 incCount(pointVisited, e[0], 1);
762 incCount(pointVisited, e[1], 1);
764 str1 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
766 else if (stat == BOTH)
768 incCount(pointVisited, e[0], 2);
769 incCount(pointVisited, e[1], 2);
771 str2 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
773 else if (stat == UNVISITED)
775 incCount(pointVisited, e[0], 0);
776 incCount(pointVisited, e[1], 0);
778 str0 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
787 label pointI = iter.key();
789 label nVisits = iter();
791 Pout<<
"point:" << pointI <<
" nVisited:" << nVisits
792 <<
" pointEdges:" << facePointEdges[pointI].size() <<
endl;
800 label visitedVert0 = -1;
801 label unvisitedVert0 = -1;
804 scalar minDist = GREAT;
808 label pointI = iter.key();
810 label nVisits = pointVisited[pointI];
812 const DynamicList<label>& pEdges = iter();
814 if (nVisits < 2*pEdges.size())
827 eSurf.points()[pointI],
834 if (nearDist < minDist)
837 visitedVert0 = nearVertI;
838 unvisitedVert0 = pointI;
846 label visitedVert1 = -1;
847 label unvisitedVert1 = -1;
849 scalar minDist = GREAT;
853 label pointI = iter.key();
855 if (pointI != unvisitedVert0)
857 label nVisits = pointVisited[pointI];
859 const DynamicList<label>& pEdges = iter();
861 if (nVisits < 2*pEdges.size())
874 eSurf.points()[pointI],
881 if (nearDist < minDist)
884 visitedVert1 = nearVertI;
885 unvisitedVert1 = pointI;
893 Pout<<
"resplitFace : adding intersection from " << visitedVert0
894 <<
" to " << unvisitedVert0 << endl
895 <<
" and from " << visitedVert1
896 <<
" to " << unvisitedVert1 <<
endl;
906 additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
908 eSurf.addIntersectionEdges(faceI, additionalEdges);
910 fileName newFName(
"face_" +
Foam::name(faceI) +
"_newEdges.obj");
911 Pout<<
"Dumping face:" << faceI <<
" to " << newFName <<
endl;
916 eSurf.faceEdges()[faceI],
921 return splitFace(surf, faceI, eSurf);
927 const triSurface& surf,
934 const edgeList& edges = eSurf.edges();
935 const labelList& fEdges = eSurf.faceEdges()[faceI];
938 Map<DynamicList<label> > facePointEdges
940 calcPointEdgeAddressing
948 Map<label> visited(fEdges.size()*2);
952 label edgeI = fEdges[i];
954 if (eSurf.isSurfaceEdge(edgeI))
958 label surfEdgeI = eSurf.parentEdge(edgeI);
960 label owner = surf.edgeOwner()[surfEdgeI];
967 surf.localFaces()[owner],
968 surf.localFaces()[faceI]
974 visited.insert(edgeI, ENDTOSTART);
980 visited.insert(edgeI, STARTTOEND);
985 visited.insert(edgeI, UNVISITED);
992 DynamicList<face> faces(fEdges.size());
1000 label startEdgeI = -1;
1001 label startVertI = -1;
1005 label edgeI = fEdges[i];
1007 const edge& e = edges[edgeI];
1009 label stat = visited[edgeI];
1011 if (stat == STARTTOEND)
1016 if (eSurf.isSurfaceEdge(edgeI))
1021 else if (stat == ENDTOSTART)
1026 if (eSurf.isSurfaceEdge(edgeI))
1033 if (startEdgeI == -1)
1056 surf.faceNormals()[faceI],
1071 label edgeI = fEdges[i];
1073 label stat = visited[edgeI];
1075 if (eSurf.isSurfaceEdge(edgeI) && stat != BOTH)
1078 <<
"Dumping face edges to faceEdges.obj" <<
endl;
1080 writeLocalOBJ(points, edges, fEdges,
"faceEdges.obj");
1083 <<
"Problem: edge " << edgeI <<
" vertices "
1084 << edges[edgeI] <<
" on face " << faceI
1085 <<
" has visited status " << stat <<
" from a "
1086 <<
"righthanded walk along all"
1087 <<
" of the triangle edges. Are the original surfaces"
1088 <<
" closed and non-intersecting?"
1091 else if (stat != BOTH)
1105 Pout<<
"** Resplitting **" <<
endl;
1124 vector n = faces[0].normal(eSurf.points());
1126 if ((n & surf.faceNormals()[faceI]) < 0)
1144 intersectionEdges_(0),
1154 intersectionEdges_(0),
1164 const bool isFirstSurface,
1169 intersectionEdges_(0),
1171 nSurfacePoints_(surf.
nPoints())
1179 faceMap_.setSize(size());
1183 faceMap_[faceI] = faceI;
1207 startTriI[faceI] = newTris.size();
1225 forAll(newFaces, newFaceI)
1227 const face& newF = newFaces[newFaceI];
1257 const label region = surf[faceI].region();
1263 const triFace& t = tris[triI];
1267 if (t[i] < 0 || t[i] >= eSurf.
points().
size())
1271 "intersectedSurface::intersectedSurface"
1272 ) <<
"Face triangulation of face " << faceI
1273 <<
" uses points outside range 0.."
1280 newTris.append(
labelledTri(t[0], t[1], t[2], region));
1308 faceMap_.setSize(size());
1310 for (label faceI = 0; faceI < surf.
size()-1; faceI++)
1312 for (label triI = startTriI[faceI]; triI < startTriI[faceI+1]; triI++)
1314 faceMap_[triI] = faceI;
1317 for (label triI = startTriI[surf.
size()-1]; triI < size(); triI++)
1319 faceMap_[triI] = surf.
size()-1;
1328 label intersectionEdgeI = 0;
1339 label surfStartI = meshPointMap()[e.
start()];
1340 label surfEndI = meshPointMap()[e.
end()];
1343 label surfEdgeI = -1;
1345 const labelList& pEdges = pointEdges()[surfStartI];
1349 const edge& surfE = edges()[pEdges[i]];
1353 if (surfE.
start() == surfEndI || surfE.
end() == surfEndI)
1355 surfEdgeI = pEdges[i];
1361 if (surfEdgeI != -1)
1363 intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1368 <<
"Cannot find edge among candidates " << pEdges
1369 <<
" which uses points " << surfStartI
1370 <<
" and " << surfEndI