55 void Foam::intersectedSurface::writeOBJ
64 const point& pt = points[pointi];
66 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
70 const edge& e = edges[edgeI];
72 os <<
"l " << e.start()+1 <<
' ' << e.end()+1 <<
nl;
78 void Foam::intersectedSurface::writeOBJ
88 const point& pt = points[pointi];
90 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
94 const edge& e = edges[faceEdges[i]];
96 os <<
"l " << e.start()+1 <<
' ' << e.end()+1 <<
nl;
102 void Foam::intersectedSurface::writeLocalOBJ
107 const fileName& fName
118 const edge& e = edges[faceEdges[i]];
124 if (pointMap[pointi] == -1)
126 const point& pt = points[pointi];
128 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
130 pointMap[pointi] = maxVertI++;
137 const edge& e = edges[faceEdges[i]];
139 os <<
"l " << pointMap[e.start()]+1 <<
' ' << pointMap[e.end()]+1
146 void Foam::intersectedSurface::writeOBJ
155 const point& pt = points[pointi];
157 os <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
164 os <<
' ' << f[fp]+1;
171 void Foam::intersectedSurface::printVisit
175 const Map<label>& visited
181 label edgeI = edgeLabels[i];
183 const edge& e = edges[edgeI];
185 label stat = visited[edgeI];
187 if (stat == UNVISITED)
189 Pout<<
" edge:" << edgeI <<
" verts:" << e
190 <<
" unvisited" <<
nl;
192 else if (stat == STARTTOEND)
194 Pout<<
" edge:" << edgeI <<
" from " << e[0]
195 <<
" to " << e[1] <<
nl;
197 else if (stat == ENDTOSTART)
199 Pout<<
" edge:" << edgeI <<
" from " << e[1]
200 <<
" to " << e[0] <<
nl;
204 Pout<<
" edge:" << edgeI <<
" both " << e
214 bool Foam::intersectedSurface::sameEdgeOrder
216 const labelledTri& fA,
217 const labelledTri& fB
227 label vA1 = fA[fA.fcIndex(fpA)];
228 label vAMin1 = fA[fA.rcIndex(fpA)];
231 label vB1 = fB[fB.fcIndex(fpB)];
232 label vBMin1 = fB[fB.rcIndex(fpB)];
234 if (vA1 == vB1 || vAMin1 == vBMin1)
238 else if (vA1 == vBMin1 || vAMin1 == vB1)
246 <<
"Triangle:" << fA <<
" and triangle:" << fB
247 <<
" share a point but not an edge" 254 <<
"Triangle:" << fA <<
" and triangle:" << fB
255 <<
" do not share an edge" 262 void Foam::intersectedSurface::incCount
271 if (iter == visited.end())
273 visited.insert(key, offset);
286 Foam::intersectedSurface::calcPointEdgeAddressing
288 const edgeSurface& eSurf,
293 const edgeList& edges = eSurf.edges();
295 const labelList& fEdges = eSurf.faceEdges()[facei];
297 Map<DynamicList<label>> facePointEdges(4*fEdges.size());
301 label edgeI = fEdges[i];
303 const edge& e = edges[edgeI];
306 Map<DynamicList<label>>::iterator iter =
307 facePointEdges.find(e.start());
309 if (iter == facePointEdges.end())
311 DynamicList<label> oneEdge;
312 oneEdge.append(edgeI);
313 facePointEdges.insert(e.start(), oneEdge);
317 iter().append(edgeI);
321 Map<DynamicList<label>>::iterator iter2 =
322 facePointEdges.find(e.end());
324 if (iter2 == facePointEdges.end())
326 DynamicList<label> oneEdge;
327 oneEdge.append(edgeI);
328 facePointEdges.insert(e.end(), oneEdge);
332 iter2().append(edgeI);
337 forAllIter(Map<DynamicList<label>>, facePointEdges, iter)
345 <<
"Point:" << iter.key() <<
" used by too few edges:" 353 Pout<<
"calcPointEdgeAddressing: face consisting of edges:" <<
endl;
356 label edgeI = fEdges[i];
357 const edge& e = edges[edgeI];
358 Pout<<
" " << edgeI <<
' ' << e
360 << points[e.end()] <<
endl;
363 Pout<<
" Constructed point-edge addressing:" <<
endl;
366 Pout<<
" vertex " << iter.key() <<
" is connected to edges " 372 return facePointEdges;
382 const edgeSurface& eSurf,
383 const Map<label>& visited,
386 const Map<DynamicList<label>>& facePointEdges,
387 const label prevEdgeI,
388 const label prevVertI
392 const edgeList& edges = eSurf.edges();
393 const labelList& fEdges = eSurf.faceEdges()[facei];
397 const DynamicList<label>& connectedEdges = facePointEdges[prevVertI];
399 if (connectedEdges.size() <= 1)
403 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
404 OFstream str(
"face.obj");
405 writeOBJ(points, edges, fEdges, str);
407 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
408 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
412 <<
"Problem: prevVertI:" << prevVertI <<
" on edge " << prevEdgeI
413 <<
" has less than 2 connected edges." 419 if (connectedEdges.size() == 2)
422 if (connectedEdges[0] == prevEdgeI)
426 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
427 <<
" to edge:" << edges[connectedEdges[1]]
428 <<
" chosen from candidates:" << connectedEdges <<
endl;
430 return connectedEdges[1];
436 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
437 <<
" to edge:" << edges[connectedEdges[0]]
438 <<
" chosen from candidates:" << connectedEdges <<
endl;
440 return connectedEdges[0];
447 const edge& prevE = edges[prevEdgeI];
450 vector e0 = n ^ (points[prevE.otherVertex(prevVertI)] - points[prevVertI]);
451 e0 /=
mag(e0) + vSmall;
456 if (
mag(
mag(e1) - 1) > small)
459 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
460 OFstream str(
"face.obj");
461 writeOBJ(points, edges, fEdges, str);
463 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
464 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
468 <<
"Unnormalized normal e1:" << e1
469 <<
" formed from cross product of e0:" << e0 <<
" n:" << n
478 scalar maxAngle = -great;
481 forAll(connectedEdges, connI)
483 label edgeI = connectedEdges[connI];
485 if (edgeI != prevEdgeI)
487 label stat = visited[edgeI];
489 const edge& e = edges[edgeI];
507 n ^ (points[e.otherVertex(prevVertI)] - points[prevVertI]);
509 vec /=
mag(vec) + vSmall;
513 if (angle > maxAngle)
527 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
528 OFstream str(
"face.obj");
529 writeOBJ(points, edges, fEdges, str);
531 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
532 writeLocalOBJ(points, edges, connectedEdges,
"faceEdges.obj");
536 <<
"Trying to step from edge " << edges[prevEdgeI]
537 <<
", vertex " << prevVertI
538 <<
" but cannot find 'unvisited' edges among candidates:" 545 Pout<<
"Stepped from edge:" << edges[prevEdgeI]
546 <<
" to edge:" << maxEdgeI <<
" angle:" << edges[maxEdgeI]
547 <<
" with angle:" << maxAngle
561 const edgeSurface& eSurf,
564 const Map<DynamicList<label>>& facePointEdges,
566 const label startEdgeI,
567 const label startVertI,
573 const edgeList& edges = eSurf.edges();
576 face
f(eSurf.faceEdges()[facei].size());
580 label vertI = startVertI;
581 label edgeI = startEdgeI;
585 const edge& e = edges[edgeI];
589 Pout<<
"Now at:" << endl
590 <<
" edge:" << edgeI <<
" vertices:" << e
591 <<
" positions:" << points[e.start()] <<
' ' << points[e.end()]
592 <<
" vertex:" << vertI <<
endl;
598 visited[edgeI] |= STARTTOEND;
602 visited[edgeI] |= ENDTOSTART;
610 vertI = e.otherVertex(vertI);
612 if (vertI == startVertI)
636 void Foam::intersectedSurface::findNearestVisited
638 const edgeSurface& eSurf,
640 const Map<DynamicList<label>>& facePointEdges,
641 const Map<label>& pointVisited,
643 const label excludePointi,
654 label pointi = iter.key();
656 if (pointi != excludePointi)
658 label nVisits = iter();
660 if (nVisits == 2*facePointEdges[pointi].size())
663 scalar dist =
mag(eSurf.points()[pointi] - pt);
676 const labelList& fEdges = eSurf.faceEdges()[facei];
679 <<
"Dumping face edges to faceEdges.obj" <<
endl;
681 writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges,
"faceEdges.obj");
684 <<
"No fully visited edge found for pt " << pt
700 const triSurface& surf,
702 const Map<DynamicList<label>>& facePointEdges,
703 const Map<label>& visited,
709 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
710 OFstream str(
"face.obj");
711 writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[facei], str);
717 Map<label> pointVisited(2*facePointEdges.size());
721 OFstream str0(
"visitedNone.obj");
722 OFstream str1(
"visitedOnce.obj");
723 OFstream str2(
"visitedTwice.obj");
724 forAll(eSurf.points(), pointi)
726 const point& pt = eSurf.points()[pointi];
728 str0 <<
"v " << pt.
x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
729 str1 <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
730 str2 <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
736 label edgeI = iter.key();
738 const edge& e = eSurf.edges()[edgeI];
742 if (stat == STARTTOEND || stat == ENDTOSTART)
744 incCount(pointVisited, e[0], 1);
745 incCount(pointVisited, e[1], 1);
747 str1 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
749 else if (stat == BOTH)
751 incCount(pointVisited, e[0], 2);
752 incCount(pointVisited, e[1], 2);
754 str2 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
756 else if (stat == UNVISITED)
758 incCount(pointVisited, e[0], 0);
759 incCount(pointVisited, e[1], 0);
761 str0 <<
"l " << e[0]+1 <<
' ' << e[1]+1 <<
nl;
770 label pointi = iter.key();
772 label nVisits = iter();
774 Pout<<
"point:" << pointi <<
" nVisited:" << nVisits
775 <<
" pointEdges:" << facePointEdges[pointi].size() <<
endl;
783 label visitedVert0 = -1;
784 label unvisitedVert0 = -1;
787 scalar minDist = great;
791 label pointi = iter.key();
793 label nVisits = pointVisited[pointi];
795 const DynamicList<label>& pEdges = iter();
797 if (nVisits < 2*pEdges.size())
810 eSurf.points()[pointi],
817 if (nearDist < minDist)
820 visitedVert0 = nearVertI;
821 unvisitedVert0 = pointi;
829 label visitedVert1 = -1;
830 label unvisitedVert1 = -1;
832 scalar minDist = great;
836 label pointi = iter.key();
838 if (pointi != unvisitedVert0)
840 label nVisits = pointVisited[pointi];
842 const DynamicList<label>& pEdges = iter();
844 if (nVisits < 2*pEdges.size())
857 eSurf.points()[pointi],
864 if (nearDist < minDist)
867 visitedVert1 = nearVertI;
868 unvisitedVert1 = pointi;
876 Pout<<
"resplitFace : adding intersection from " << visitedVert0
877 <<
" to " << unvisitedVert0 << endl
878 <<
" and from " << visitedVert1
879 <<
" to " << unvisitedVert1 <<
endl;
889 additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
891 eSurf.addIntersectionEdges(facei, additionalEdges);
893 fileName newFName(
"face_" +
Foam::name(facei) +
"_newEdges.obj");
894 Pout<<
"Dumping face:" << facei <<
" to " << newFName <<
endl;
899 eSurf.faceEdges()[facei],
904 return splitFace(surf, facei, eSurf);
910 const triSurface& surf,
917 const edgeList& edges = eSurf.edges();
918 const labelList& fEdges = eSurf.faceEdges()[facei];
921 Map<DynamicList<label>> facePointEdges
923 calcPointEdgeAddressing
931 Map<label> visited(fEdges.size()*2);
935 label edgeI = fEdges[i];
937 if (eSurf.isSurfaceEdge(edgeI))
941 label surfEdgeI = eSurf.parentEdge(edgeI);
943 label owner = surf.edgeOwner()[surfEdgeI];
950 surf.localFaces()[owner],
951 surf.localFaces()[facei]
957 visited.insert(edgeI, ENDTOSTART);
963 visited.insert(edgeI, STARTTOEND);
968 visited.insert(edgeI, UNVISITED);
975 DynamicList<face> faces(fEdges.size());
983 label startEdgeI = -1;
984 label startVertI = -1;
988 label edgeI = fEdges[i];
990 const edge& e = edges[edgeI];
992 label stat = visited[edgeI];
994 if (stat == STARTTOEND)
999 if (eSurf.isSurfaceEdge(edgeI))
1004 else if (stat == ENDTOSTART)
1009 if (eSurf.isSurfaceEdge(edgeI))
1016 if (startEdgeI == -1)
1039 surf.faceNormals()[facei],
1054 label edgeI = fEdges[i];
1056 label stat = visited[edgeI];
1058 if (eSurf.isSurfaceEdge(edgeI) && stat != BOTH)
1061 <<
"Dumping face edges to faceEdges.obj" <<
endl;
1063 writeLocalOBJ(points, edges, fEdges,
"faceEdges.obj");
1066 <<
"Problem: edge " << edgeI <<
" vertices " 1067 << edges[edgeI] <<
" on face " << facei
1068 <<
" has visited status " << stat <<
" from a " 1069 <<
"righthanded walk along all" 1070 <<
" of the triangle edges. Are the original surfaces" 1071 <<
" closed and non-intersecting?" 1074 else if (stat != BOTH)
1088 Pout<<
"** Resplitting **" <<
endl;
1107 vector a = faces[0].area(eSurf.points());
1109 if ((a & surf.faceNormals()[facei]) < 0)
1127 intersectionEdges_(0),
1137 intersectionEdges_(0),
1147 const bool isFirstSurface,
1152 intersectionEdges_(0),
1154 nSurfacePoints_(surf.
nPoints())
1166 faceMap_[facei] = facei;
1190 startTriI[facei] = newTris.size();
1208 forAll(newFaces, newFacei)
1210 const face& newF = newFaces[newFacei];
1240 const label region = surf[facei].region();
1246 const triFace& t = tris[triI];
1250 if (t[i] < 0 || t[i] >= eSurf.
points().
size())
1253 <<
"Face triangulation of face " << facei
1254 <<
" uses points outside range 0.." 1261 newTris.append(
labelledTri(t[0], t[1], t[2], region));
1291 for (
label facei = 0; facei < surf.
size()-1; facei++)
1293 for (
label triI = startTriI[facei]; triI < startTriI[facei+1]; triI++)
1295 faceMap_[triI] = facei;
1298 for (
label triI = startTriI[surf.
size()-1]; triI <
size(); triI++)
1300 faceMap_[triI] = surf.
size()-1;
1309 label intersectionEdgeI = 0;
1324 label surfEdgeI = -1;
1334 if (surfE.
start() == surfEndI || surfE.
end() == surfEndI)
1336 surfEdgeI = pEdges[i];
1342 if (surfEdgeI != -1)
1344 intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1349 <<
"Cannot find edge among candidates " << pEdges
1350 <<
" which uses points " << surfStartI
1351 <<
" and " << surfEndI
const labelListList & pointEdges() const
Return point-edge addressing.
label nPoints() const
Return number of points supporting patch faces.
Triangulation of faces. Handles concave polygons as well (inefficiently)
#define forAll(list, i)
Loop across all elements in list.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
label nSurfacePoints() const
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
friend Ostream & operator(Ostream &, const UList< T > &)
A face is a list of labels corresponding to mesh vertices.
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
void operator=(const triSurface &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nSurfaceEdges() const
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
HashTable< label, label, Hash< label > >::iterator iterator
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
void size(const label)
Override size to be inconsistent with allocated storage.
intersectedSurface()
Construct null.
const pointField & points() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
static const label UNVISITED
Description of surface in form of 'cloud of edges'.
Vector< scalar > vector
A scalar version of the templated Vector.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
const edgeList & cutEdges() const
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
vectorField pointField
pointField is a vectorField.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A triangular face using a FixedList of labels corresponding to mesh vertices.
const edgeList & edges() const
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Triangle with additional region number.
const Field< PointType > & faceNormals() const
Return face normals for patch.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
void reverse(UList< T > &, const label n)
static const label STARTTOEND
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
const pointField & cutPoints() const
word name(const complex &)
Return a string representation of a complex.
void setSize(const label)
Reset size of List.
const geometricSurfacePatchList & patches() const
const labelListList & faceEdges() const
From face to our edges_.
triSurface()
Construct null.
vector point
Point is a vector.
label end() const
Return end vertex label.
prefixOSstream Pout(cout, "Pout")
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Triangulated surface description with patch information.
label size() const
Return the number of elements in the UList.
static const label ENDTOSTART
label start() const
Return start vertex label.
A HashTable to objects of type <T> with a label key.