55 void Foam::intersectedSurface::writeOBJ
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
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)
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
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
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");
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];
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");
463 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
464 writeLocalOBJ(
points, edges, connectedEdges,
"faceEdges.obj");
468 <<
"Unnormalised 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];
509 vec /=
mag(vec) + vSmall;
513 if (angle > maxAngle)
527 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
528 OFstream str(
"face.obj");
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];
590 <<
" edge:" << edgeI <<
" vertices:" <<
e
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;
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],
820 visitedVert0 = nearVertI;
821 unvisitedVert0 = pointi;
829 label visitedVert1 = -1;
830 label unvisitedVert1 = -1;
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],
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;
1193 startTriI[facei] = newTris.
size();
1211 forAll(newFaces, newFacei)
1213 const face& newF = newFaces[newFacei];
1216 const label region = surf[facei].region();
1218 triEngine.triangulate
1259 for (
label facei = 0; facei < surf.
size()-1; facei++)
1261 for (
label triI = startTriI[facei]; triI < startTriI[facei+1]; triI++)
1263 faceMap_[triI] = facei;
1266 for (
label triI = startTriI[surf.
size()-1]; triI <
size(); triI++)
1268 faceMap_[triI] = surf.
size()-1;
1277 label intersectionEdgeI = 0;
1292 label surfEdgeI = -1;
1302 if (surfE.
start() == surfEndI || surfE.
end() == surfEndI)
1304 surfEdgeI = pEdges[i];
1310 if (surfEdgeI != -1)
1312 intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1317 <<
"Cannot find edge among candidates " << pEdges
1318 <<
" which uses points " << surfStartI
1319 <<
" and " << surfEndI
#define forAll(list, i)
Loop across all elements in list.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
friend class iterator
Declare friendship with the iterator.
label size() const
Return the number of elements in the UList.
void size(const label)
Override size to be inconsistent with allocated storage.
void setSize(const label)
Reset size of List.
A HashTable to objects of type <T> with a label key.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
const Field< PointType > & points() const
Return reference to global points.
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
const labelListList & faceEdges() const
Return face-edge addressing.
const Field< PointType > & faceNormals() const
Return face normals for patch.
A List with indirect addressing.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
friend Ostream & operator(Ostream &, const UList< T > &)
Description of surface in form of 'cloud of edges'.
label nSurfaceEdges() const
const labelListList & faceEdges() const
From face to our edges_.
label nSurfacePoints() const
const edgeList & edges() const
const pointField & points() const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
label end() const
Return end vertex label.
label start() const
Return start vertex label.
A face is a list of labels corresponding to mesh vertices.
Given triSurface and intersection creates the intersected (properly triangulated) surface....
intersectedSurface()
Construct null.
static const label ENDTOSTART
static const label UNVISITED
static const label STARTTOEND
Triangle with additional region number.
Triangulation of three-dimensional polygons.
const UList< triFace > & triPoints() const
Get the triangles' points.
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
const edgeList & cutEdges() const
const pointField & cutPoints() const
Triangulated surface description with patch information.
triSurface()
Construct null.
void operator=(const triSurface &)
const geometricSurfacePatchList & patches() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
scalar minDist(const List< pointIndexHit > &hitList)
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
word name(const bool)
Return a word representation of a bool.
errorManip< error > abort(error &err)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
void reverse(UList< T > &, const label n)
Vector< scalar > vector
A scalar version of the templated Vector.
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
prefixOSstream Pout(cout, "Pout")
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void offset(label &lst, const label o)