57 void Foam::intersectedSurface::writeOBJ
66 os <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
69 for (
const edge&
e :
edges)
71 os <<
"l " <<
e.start()+1 <<
' ' <<
e.end()+1 <<
nl;
76 void Foam::intersectedSurface::writeOBJ
86 os <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
88 for (
const label edgei : faceEdges)
90 const edge&
e = edges[edgei];
92 os <<
"l " <<
e.start()+1 <<
' ' <<
e.end()+1 <<
nl;
97 void Foam::intersectedSurface::writeLocalOBJ
102 const fileName& fName
111 for (
const label edgei : faceEdges)
113 const edge&
e = edges[edgei];
117 const label pointi =
e[i];
119 if (pointMap[pointi] == -1)
123 os <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
125 pointMap[pointi] = maxVerti++;
130 for (
const label edgei : faceEdges)
132 const edge&
e = edges[edgei];
134 os <<
"l " << pointMap[
e.start()]+1 <<
' ' << pointMap[
e.end()]+1
140 void Foam::intersectedSurface::writeOBJ
149 os <<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
nl;
154 for (
const label pointi :
f)
156 os <<
' ' << pointi+1;
162 void Foam::intersectedSurface::printVisit
166 const Map<label>& visited
170 for (
const label edgei : edgeLabels)
172 const edge&
e = edges[edgei];
174 label stat = visited[edgei];
176 if (stat == UNVISITED)
178 Pout<<
" edge:" << edgei <<
" verts:" <<
e 179 <<
" unvisited" <<
nl;
181 else if (stat == STARTTOEND)
183 Pout<<
" edge:" << edgei <<
" from " <<
e[0]
184 <<
" to " <<
e[1] <<
nl;
186 else if (stat == ENDTOSTART)
188 Pout<<
" edge:" << edgei <<
" from " <<
e[1]
189 <<
" to " <<
e[0] <<
nl;
193 Pout<<
" edge:" << edgei <<
" both " <<
e 201 bool Foam::intersectedSurface::sameEdgeOrder
203 const labelledTri& fA,
204 const labelledTri& fB
209 label fpB = fB.find(fA[fpA]);
214 label vA1 = fA[fA.fcIndex(fpA)];
215 label vAMin1 = fA[fA.rcIndex(fpA)];
218 label vB1 = fB[fB.fcIndex(fpB)];
219 label vBMin1 = fB[fB.rcIndex(fpB)];
221 if (vA1 == vB1 || vAMin1 == vBMin1)
225 else if (vA1 == vBMin1 || vAMin1 == vB1)
233 <<
"Triangle:" << fA <<
" and triangle:" << fB
234 <<
" share a point but not an edge" 241 <<
"Triangle:" << fA <<
" and triangle:" << fB
242 <<
" do not share an edge" 249 void Foam::intersectedSurface::incCount
256 visited(
key, 0) += offset;
261 Foam::intersectedSurface::calcPointEdgeAddressing
263 const edgeSurface& eSurf,
268 const edgeList& edges = eSurf.edges();
270 const labelList& fEdges = eSurf.faceEdges()[facei];
272 Map<DynamicList<label>> facePointEdges(4*fEdges.size());
274 for (
const label edgei : fEdges)
276 const edge&
e = edges[edgei];
279 facePointEdges(
e.start()).
append(edgei);
282 facePointEdges(
e.end()).
append(edgei);
291 if (iter.val().empty())
294 <<
"Point:" << iter.key() <<
" used by too few edges:" 302 Pout<<
"calcPointEdgeAddressing: face consisting of edges:" <<
nl;
303 for (
const label edgei : fEdges)
305 const edge&
e = edges[edgei];
306 Pout<<
" " << edgei <<
' ' <<
e 311 Pout<<
" Constructed point-edge addressing:" <<
nl;
314 Pout<<
" vertex " << iter.key() <<
" is connected to edges " 320 return facePointEdges;
324 Foam::label Foam::intersectedSurface::nextEdge
326 const edgeSurface& eSurf,
327 const Map<label>& visited,
330 const Map<DynamicList<label>>& facePointEdges,
331 const label prevEdgei,
332 const label prevVerti
336 const edgeList& edges = eSurf.edges();
337 const labelList& fEdges = eSurf.faceEdges()[facei];
341 const DynamicList<label>& connectedEdges = facePointEdges[prevVerti];
343 if (connectedEdges.size() <= 1)
347 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
348 OFstream str(
"face.obj");
351 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
352 writeLocalOBJ(
points, edges, connectedEdges,
"faceEdges.obj");
356 <<
"Problem: prevVertI:" << prevVerti <<
" on edge " << prevEdgei
357 <<
" has less than 2 connected edges." 363 if (connectedEdges.size() == 2)
366 if (connectedEdges[0] == prevEdgei)
370 Pout<<
"Stepped from edge:" << edges[prevEdgei]
371 <<
" to edge:" << edges[connectedEdges[1]]
372 <<
" chosen from candidates:" << connectedEdges <<
endl;
374 return connectedEdges[1];
380 Pout<<
"Stepped from edge:" << edges[prevEdgei]
381 <<
" to edge:" << edges[connectedEdges[0]]
382 <<
" chosen from candidates:" << connectedEdges <<
endl;
384 return connectedEdges[0];
391 const edge& prevE = edges[prevEdgei];
397 n ^ (
points[prevE.otherVertex(prevVerti)] -
points[prevVerti])
403 if (
mag(
mag(e1) - 1) > SMALL)
406 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
407 OFstream str(
"face.obj");
410 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
411 writeLocalOBJ(
points, edges, connectedEdges,
"faceEdges.obj");
415 <<
"Unnormalized normal e1:" << e1
416 <<
" formed from cross product of e0:" << e0 <<
" n:" <<
n 425 scalar maxAngle = -GREAT;
428 for (
const label edgei : connectedEdges)
430 if (edgei != prevEdgei)
432 label stat = visited[edgei];
434 const edge&
e = edges[edgei];
460 if (angle > maxAngle)
474 Pout<<
"Writing face:" << facei <<
" to face.obj" <<
endl;
475 OFstream str(
"face.obj");
478 Pout<<
"Writing connectedEdges edges to faceEdges.obj" <<
endl;
479 writeLocalOBJ(
points, edges, connectedEdges,
"faceEdges.obj");
483 <<
"Trying to step from edge " << edges[prevEdgei]
484 <<
", vertex " << prevVerti
485 <<
" but cannot find 'unvisited' edges among candidates:" 492 Pout<<
"Stepped from edge:" << edges[prevEdgei]
493 <<
" to edge:" << maxEdgei <<
" angle:" << edges[maxEdgei]
494 <<
" with angle:" << maxAngle
504 const edgeSurface& eSurf,
507 const Map<DynamicList<label>>& facePointEdges,
509 const label startEdgei,
510 const label startVerti,
516 const edgeList& edges = eSurf.edges();
519 face
f(eSurf.faceEdges()[facei].size());
523 label verti = startVerti;
524 label edgei = startEdgei;
528 const edge&
e = edges[edgei];
533 <<
" edge:" << edgei <<
" vertices:" <<
e 535 <<
" vertex:" << verti <<
endl;
541 visited[edgei] |= STARTTOEND;
545 visited[edgei] |= ENDTOSTART;
553 verti =
e.otherVertex(verti);
555 if (verti == startVerti)
579 void Foam::intersectedSurface::findNearestVisited
581 const edgeSurface& eSurf,
583 const Map<DynamicList<label>>& facePointEdges,
584 const Map<label>& pointVisited,
586 const label excludePointi,
597 const label pointi = iter.key();
598 const label nVisits = iter.val();
600 if (pointi != excludePointi)
602 if (nVisits == 2*facePointEdges[pointi].size())
605 scalar dist =
mag(eSurf.points()[pointi] - pt);
618 const labelList& fEdges = eSurf.faceEdges()[facei];
621 <<
"Dumping face edges to faceEdges.obj" <<
endl;
623 writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges,
"faceEdges.obj");
626 <<
"No fully visited edge found for pt " << pt
634 const triSurface& surf,
636 const Map<DynamicList<label>>& facePointEdges,
637 const Map<label>& visited,
643 Map<label> pointVisited(2*facePointEdges.size());
647 const label edgei = iter.key();
648 const label stat = iter.val();
650 const edge&
e = eSurf.edges()[edgei];
652 if (stat == STARTTOEND || stat == ENDTOSTART)
654 incCount(pointVisited,
e[0], 1);
655 incCount(pointVisited,
e[1], 1);
657 else if (stat == BOTH)
659 incCount(pointVisited,
e[0], 2);
660 incCount(pointVisited,
e[1], 2);
662 else if (stat == UNVISITED)
664 incCount(pointVisited,
e[0], 0);
665 incCount(pointVisited,
e[1], 0);
674 const label pointi = iter.key();
675 const label nVisits = iter.val();
677 Pout<<
"point:" << pointi
678 <<
" nVisited:" << nVisits
679 <<
" pointEdges:" << facePointEdges[pointi].size() <<
endl;
687 label visitedVert0 = -1;
688 label unvisitedVert0 = -1;
691 scalar minDist = GREAT;
695 const label pointi = iter.key();
697 const DynamicList<label>& pEdges = iter.val();
699 const label nVisits = pointVisited[pointi];
701 if (nVisits < 2*pEdges.size())
714 eSurf.points()[pointi],
721 if (nearDist < minDist)
724 visitedVert0 = nearVerti;
725 unvisitedVert0 = pointi;
734 scalar minDist = GREAT;
738 const label pointi = iter.key();
740 const DynamicList<label>& pEdges = iter.val();
742 if (pointi != unvisitedVert0)
744 const label nVisits = pointVisited[pointi];
746 if (nVisits < 2*pEdges.size())
759 eSurf.points()[pointi],
766 if (nearDist < minDist)
778 additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
780 eSurf.addIntersectionEdges(facei, additionalEdges);
784 fileName newFName(
"face_" +
Foam::name(facei) +
"_newEdges.obj");
785 Pout<<
"Dumping face:" << facei <<
" to " << newFName <<
endl;
790 eSurf.faceEdges()[facei],
796 return splitFace(surf, facei, eSurf);
802 const triSurface& surf,
809 const edgeList& edges = eSurf.edges();
810 const labelList& fEdges = eSurf.faceEdges()[facei];
813 Map<DynamicList<label>> facePointEdges
815 calcPointEdgeAddressing
823 Map<label> visited(fEdges.size()*2);
827 label edgei = fEdges[i];
829 if (eSurf.isSurfaceEdge(edgei))
833 label surfEdgei = eSurf.parentEdge(edgei);
835 label owner = surf.edgeOwner()[surfEdgei];
842 surf.localFaces()[owner],
843 surf.localFaces()[facei]
849 visited.insert(edgei, ENDTOSTART);
855 visited.insert(edgei, STARTTOEND);
860 visited.insert(edgei, UNVISITED);
867 DynamicList<face> faces(fEdges.size());
875 label startEdgei = -1;
876 label startVerti = -1;
880 label edgei = fEdges[i];
882 const edge&
e = edges[edgei];
884 label stat = visited[edgei];
886 if (stat == STARTTOEND)
891 if (eSurf.isSurfaceEdge(edgei))
896 else if (stat == ENDTOSTART)
901 if (eSurf.isSurfaceEdge(edgei))
908 if (startEdgei == -1)
931 surf.faceNormals()[facei],
946 label edgei = fEdges[i];
948 label stat = visited[edgei];
950 if (eSurf.isSurfaceEdge(edgei) && stat != BOTH)
953 <<
"Dumping face edges to faceEdges.obj" <<
endl;
955 writeLocalOBJ(
points, edges, fEdges,
"faceEdges.obj");
958 <<
"Problem: edge " << edgei <<
" vertices " 959 << edges[edgei] <<
" on face " << facei
960 <<
" has visited status " << stat <<
" from a " 961 <<
"right-handed walk along all" 962 <<
" of the triangle edges. Are the original surfaces" 963 <<
" closed and non-intersecting?" 966 else if (stat != BOTH)
985 const vector n = faces[0].areaNormal(eSurf.points());
987 if ((
n & surf.faceNormals()[facei]) < 0)
989 for (face&
f : faces)
1004 intersectionEdges_(0),
1013 intersectionEdges_(0),
1022 const bool isFirstSurface,
1027 intersectionEdges_(0),
1029 nSurfacePoints_(surf.
nPoints())
1041 faceMap_[facei] = facei;
1048 edgeSurface eSurf(surf, isFirstSurface, inter);
1052 nSurfacePoints_ = eSurf.nSurfacePoints();
1058 DynamicList<labelledTri> newTris(eSurf.edges().size()/2);
1065 startTrii[facei] = newTris.size();
1067 if (eSurf.faceEdges()[facei].size() != surf.
faceEdges()[facei].
size())
1083 forAll(newFaces, newFacei)
1085 const face& newF = newFaces[newFacei];
1087 const label region = surf[facei].region();
1089 faceTriangulation tris(eSurf.points(), newF,
n);
1093 const triFace& t = tris[trii];
1097 if (t[i] < 0 || t[i] >= eSurf.points().size())
1100 <<
"Face triangulation of face " << facei
1101 <<
" uses points outside range 0.." 1102 << eSurf.points().size()-1 <<
endl 1108 newTris.append(labelledTri(t[0], t[1], t[2], region));
1125 triSurface::operator=
1138 for (label facei = 0; facei < surf.
size()-1; facei++)
1140 for (label trii = startTrii[facei]; trii < startTrii[facei+1]; trii++)
1142 faceMap_[trii] = facei;
1145 for (label trii = startTrii[surf.
size()-1]; trii <
size(); trii++)
1147 faceMap_[trii] = surf.
size()-1;
1154 intersectionEdges_.
setSize(eSurf.edges().size() - eSurf.nSurfaceEdges());
1156 label intersectionEdgei = 0;
1160 label edgei = eSurf.nSurfaceEdges();
1161 edgei < eSurf.edges().size();
1166 const edge&
e = eSurf.edges()[edgei];
1171 label surfEdgei = -1;
1177 const edge& surfE =
edges()[pEdges[i]];
1181 if (surfE.found(surfEndi))
1183 surfEdgei = pEdges[i];
1189 if (surfEdgei != -1)
1191 intersectionEdges_[intersectionEdgei++] = surfEdgei;
1196 <<
"Cannot find edge among candidates " << pEdges
1197 <<
" which uses points " << surfStarti
1198 <<
" and " << surfEndi
const labelListList & pointEdges() const
Return point-edge addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
void size(const label n)
Older name for setAddressableSize.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
constexpr char nl
The newline '\n' character (0x0a)
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
List< face > faceList
A List of faces.
intersectedSurface()
Construct null.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static const label UNVISITED
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
#define forAll(list, i)
Loop across all elements in list.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const edgeList & cutEdges() const
The list of created edges.
const Map< label > & meshPointMap() const
Mesh point map.
const geometricSurfacePatchList & patches() const noexcept
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
List< edge > edgeList
A List of edges.
#define forAllIters(container, iter)
Iterate across all elements in the container object.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
const Field< point_type > & points() const noexcept
Return reference to global points.
errorManip< error > abort(error &err)
int debug
Static debugging option.
static const label STARTTOEND
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
const pointField & cutPoints() const
The list of cut points.
rAUs append(new volScalarField(IOobject::groupName("rAU", phase1.name()), 1.0/(U1Eqn.A()+byDt(max(phase1.residualAlpha() - alpha1, scalar(0)) *rho1))))
triSurface()
Default construct.
vector point
Point is a vector.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
label size() const noexcept
The number of elements in the UList.
const labelListList & faceEdges() const
Return face-edge addressing.
List< label > labelList
A List of labels.
Triangulated surface description with patch information.
void operator=(const triSurface &surf)
Copy assignment.
static const label ENDTOSTART
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
forAllConstIters(mixture.phases(), phase)
A HashTable to objects of type <T> with a label key.
static constexpr const zero Zero
Global zero (0)