55 void Foam::triSurfaceTools::calcRefineStatus
57 const triSurface& surf,
59 List<refineType>& refine
62 if (refine[faceI] == RED)
71 const labelList& myNeighbours = surf.faceFaces()[faceI];
73 forAll(myNeighbours, myNeighbourI)
75 label neighbourFaceI = myNeighbours[myNeighbourI];
77 if (refine[neighbourFaceI] == GREEN)
80 calcRefineStatus(surf, neighbourFaceI, refine);
82 else if (refine[neighbourFaceI] == NONE)
84 refine[neighbourFaceI] = GREEN;
92 void Foam::triSurfaceTools::greenRefine
94 const triSurface& surf,
97 const label newPointI,
98 DynamicList<labelledTri>& newFaces
101 const labelledTri&
f = surf.localFaces()[faceI];
102 const edge&
e = surf.edges()[edgeI];
107 label fp1 = f.fcIndex(fp0);
108 label fp2 = f.fcIndex(fp1);
163 const triSurface& surf,
164 const List<refineType>& refineStatus
168 DynamicList<point> newPoints(surf.nPoints());
169 forAll(surf.localPoints(), pointI)
171 newPoints.append(surf.localPoints()[pointI]);
173 label newVertI = surf.nPoints();
176 DynamicList<labelledTri> newFaces(surf.size());
182 forAll(refineStatus, faceI)
184 if (refineStatus[faceI] == RED)
187 const labelList& fEdges = surf.faceEdges()[faceI];
191 label edgeI = fEdges[i];
193 if (edgeMid[edgeI] == -1)
195 const edge& e = surf.edges()[edgeI];
202 surf.localPoints()[e.start()]
203 + surf.localPoints()[e.end()]
206 edgeMid[edgeI] = newVertI++;
213 const edgeList& edges = surf.edges();
221 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
232 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
243 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
265 const label edgeI = fEdges[i];
267 label otherFaceI =
otherFace(surf, faceI, edgeI);
269 if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN))
285 forAll(refineStatus, faceI)
287 if (refineStatus[faceI] == NONE)
289 newFaces.append(surf.localFaces()[faceI]);
302 return triSurface(newFaces, surf.patches(), allPoints,
true);
308 Foam::scalar Foam::triSurfaceTools::faceCosAngle
316 const vector common(pEnd - pStart);
317 const vector base0(pLeft - pStart);
318 const vector base1(pRight - pStart);
320 vector n0(common ^ base0);
323 vector n1(base1 ^ common);
334 void Foam::triSurfaceTools::protectNeighbours
336 const triSurface& surf,
352 const labelList& myEdges = surf.pointEdges()[vertI];
355 const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
359 label faceI = myFaces[myFaceI];
361 if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
363 faceStatus[faceI] = NOEDGE;
377 const triSurface& surf,
381 const edge& e = surf.edges()[edgeI];
382 label
v1 = e.start();
386 const labelList& myFaces = surf.edgeFaces()[edgeI];
392 facesToBeCollapsed.insert(myFaces[myFaceI]);
403 label face1I = v1Faces[v1FaceI];
405 label otherEdgeI = oppositeEdge(surf, face1I, v1);
408 label face2I =
otherFace(surf, face1I, otherEdgeI);
413 if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
417 facesToBeCollapsed.insert(face1I);
418 facesToBeCollapsed.insert(face2I);
423 return facesToBeCollapsed;
428 Foam::label Foam::triSurfaceTools::vertexUsesFace
430 const triSurface& surf,
435 const labelList& myFaces = surf.pointFaces()[vertI];
439 label face1I = myFaces[myFaceI];
441 if (faceUsed.found(face1I))
451 void Foam::triSurfaceTools::getMergedEdges
453 const triSurface& surf,
456 HashTable<label, label, Hash<label> >& edgeToEdge,
457 HashTable<label, label, Hash<label> >& edgeToFace
460 const edge& e = surf.edges()[edgeI];
461 label v1 = e.start();
472 if (!collapsedFaces.found(v2Faces[v2FaceI]))
474 v2FacesHash.insert(v2Faces[v2FaceI]);
481 label face1I = v1Faces[v1FaceI];
483 if (collapsedFaces.found(face1I))
507 label commonVert = vert1I;
508 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
512 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
518 label edge1I = getEdge(surf, v1, commonVert);
519 label edge2I = getEdge(surf, v2, commonVert);
521 edgeToEdge.insert(edge1I, edge2I);
522 edgeToEdge.insert(edge2I, edge1I);
524 edgeToFace.insert(edge1I, face2I);
525 edgeToFace.insert(edge2I, face1I);
533 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
535 const triSurface& surf,
539 const HashTable<label, label, Hash<label> >& edgeToEdge,
540 const HashTable<label, label, Hash<label> >& edgeToFace,
545 const pointField& localPoints = surf.localPoints();
547 label
A = surf.edges()[edgeI].start();
548 label B = surf.edges()[edgeI].
end();
549 label C = oppositeVertex(surf, faceI, edgeI);
555 if (edgeToEdge.found(edgeI))
558 label edge2I = edgeToEdge[edgeI];
559 face2I = edgeToFace[edgeI];
561 D = oppositeVertex(surf, face2I, edge2I);
568 if ((face2I != -1) && !collapsedFaces.found(face2I))
570 D = oppositeVertex(surf, face2I, edgeI);
580 cosAngle = faceCosAngle
590 cosAngle = faceCosAngle
600 cosAngle = faceCosAngle
610 cosAngle = faceCosAngle
621 <<
"face " << faceI <<
" does not use vertex "
631 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
633 const triSurface& surf,
637 const HashTable<label, label, Hash<label> >& edgeToEdge,
638 const HashTable<label, label, Hash<label> >& edgeToFace
647 label faceI = v1Faces[v1FaceI];
649 if (collapsedFaces.found(faceI))
654 const labelList& myEdges = surf.faceEdges()[faceI];
658 label edgeI = myEdges[myEdgeI];
685 bool Foam::triSurfaceTools::collapseCreatesFold
687 const triSurface& surf,
691 const HashTable<label, label, Hash<label> >& edgeToEdge,
692 const HashTable<label, label, Hash<label> >& edgeToFace,
700 label faceI = v1Faces[v1FaceI];
702 if (collapsedFaces.found(faceI))
707 const labelList& myEdges = surf.faceEdges()[faceI];
711 label edgeI = myEdges[myEdgeI];
848 const label excludeEdgeI,
849 const label excludePointI,
851 const point& triPoint,
852 const plane& cutPlane,
857 const labelledTri& f = s[triI];
858 const labelList& fEdges = s.faceEdges()[triI];
861 FixedList<scalar, 3>
d;
866 d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
875 if (
mag(d[i]) < 1
E-6)
884 if (excludePointI != -1)
888 label fp0 =
findIndex(s.localFaces()[triI], excludePointI);
892 FatalErrorIn(
"cutEdge(..)") <<
"excludePointI:" << excludePointI
896 label fp1 = f.fcIndex(fp0);
897 label fp2 = f.fcIndex(fp1);
903 cut.setPoint(points[f[fp1]]);
905 cut.setIndex(s.localFaces()[triI][fp1]);
907 else if (d[fp2] == 0.0)
910 cut.setPoint(points[f[fp2]]);
912 cut.setIndex(s.localFaces()[triI][fp2]);
916 (d[fp1] < 0 && d[fp2] < 0)
917 || (d[fp1] > 0 && d[fp2] > 0)
928 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
932 cut.setIndex(fEdges[fp1]);
938 FixedList<surfaceLocation, 2> inters;
943 label fp1 = f.fcIndex(fp0);
950 <<
"problem : triangle has three intersections." <<
nl
951 <<
"triangle:" << f.tri(points)
954 inters[interI].setHit();
955 inters[interI].setPoint(points[f[fp0]]);
957 inters[interI].setIndex(s.localFaces()[triI][fp0]);
962 (d[fp0] < 0 && d[fp1] > 0)
963 || (d[fp0] > 0 && d[fp1] < 0)
969 <<
"problem : triangle has three intersections." <<
nl
970 <<
"triangle:" << f.tri(points)
973 inters[interI].setHit();
974 inters[interI].setPoint
976 (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
980 inters[interI].setIndex(fEdges[fp0]);
990 else if (interI == 1)
995 else if (interI == 2)
1001 && inters[0].index() == excludeEdgeI
1009 && inters[1].index() == excludeEdgeI
1019 magSqr(inters[0].rawPoint() - toPoint)
1020 <
magSqr(inters[1].rawPoint() - toPoint)
1037 void Foam::triSurfaceTools::snapToEnd
1039 const triSurface& s,
1040 const surfaceLocation& end,
1041 surfaceLocation& current
1049 if (current.index() == end.index())
1067 const labelList& fEdges = s.faceEdges()[current.index()];
1069 if (
findIndex(fEdges, end.index()) != -1)
1083 if (current.index() == end.index())
1097 const edge& e = s.edges()[end.index()];
1099 if (current.index() == e[0] || current.index() == e[1])
1116 const labelledTri& f = s.localFaces()[current.index()];
1132 const edge& e = s.edges()[current.index()];
1134 if (end.index() == e[0] || end.index() == e[1])
1148 if (current.index() == end.index())
1169 const triSurface& s,
1171 const surfaceLocation& start,
1172 const label excludeEdgeI,
1173 const label excludePointI,
1174 const surfaceLocation& end,
1175 const plane& cutPlane
1178 surfaceLocation nearest;
1184 label triI = eFaces[i];
1187 if (triI != start.triangle())
1194 nearest.triangle() = triI;
1201 surfaceLocation cutInfo = cutEdge
1213 if (excludeEdgeI != -1 && !cutInfo.hit())
1216 <<
"Triangle:" << triI
1217 <<
" excludeEdge:" << excludeEdgeI
1218 <<
" point:" << start.rawPoint()
1219 <<
" plane:" << cutPlane
1225 scalar distSqr =
magSqr(cutInfo.rawPoint()-end.rawPoint());
1227 if (distSqr < minDistSqr)
1229 minDistSqr = distSqr;
1231 nearest.triangle() = triI;
1239 if (nearest.triangle() == -1)
1264 const point& pt = pts[pointI];
1266 outFile<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
1268 Pout<<
"Written " << pts.
size() <<
" vertices to file " << fName <<
endl;
1283 forAll(markedVerts, vertI)
1285 if (markedVerts[vertI])
1289 outFile<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
1294 Pout<<
"Written " << nVerts <<
" vertices to file " << fName <<
endl;
1313 label face1I = myFaces[0];
1315 if (myFaces.
size() == 2)
1317 face2I = myFaces[1];
1328 forAll(startFaces, startFaceI)
1330 edgeTris[nTris++] = startFaces[startFaceI];
1333 forAll(endFaces, endFaceI)
1335 label faceI = endFaces[endFaceI];
1337 if ((faceI != face1I) && (faceI != face2I))
1339 edgeTris[nTris++] = faceI;
1354 label v1 = e.
start();
1364 const edge& e = edges[v1Edges[v1EdgeI]];
1372 const edge& e = edges[v2Edges[v2EdgeI]];
1376 vertexNeighbours.
insert(vertI);
1378 return vertexNeighbours.
toc();
1440 if (myFaces.
size() != 2)
1446 if (faceI == myFaces[0])
1477 "(const triSurface&, const label, const label,"
1479 ) <<
"Edge " << surf.
edges()[edgeI] <<
" not in face "
1483 label i1 = eFaces.
fcIndex(i0);
1484 label i2 = eFaces.
fcIndex(i1);
1508 else if (vertI == f[1])
1513 else if (vertI == f[2])
1523 "(const triSurface&, const label, const label,"
1542 label edgeI = myEdges[myEdgeI];
1546 if ((e.
start() != vertI) && (e.
end() != vertI))
1555 "(const triSurface&, const label, const label)"
1556 ) <<
"Cannot find vertex " << vertI <<
" in edges of face " << faceI
1577 label vertI = f[fp];
1579 if (vertI != e.
start() && vertI != e.
end())
1586 <<
"Cannot find vertex opposite edge " << edgeI <<
" vertices " << e
1605 label edgeI = v1Edges[v1EdgeI];
1627 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1632 "(const triSurface&, const label, const label,"
1634 ) <<
"Duplicate edge labels : e0:" << e0I <<
" e1:" << e1I
1643 label faceI = eFaces[eFaceI];
1650 || (myEdges[1] == e1I)
1651 || (myEdges[2] == e1I)
1657 || (myEdges[1] == e2I)
1658 || (myEdges[2] == e2I)
1718 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1740 forAll(localPoints, pointI)
1742 pointMap[pointI] = pointI;
1748 forAll(collapseEdgeLabels, collapseEdgeI)
1750 const label edgeI = collapseEdgeLabels[collapseEdgeI];
1752 if ((edgeI < 0) || (edgeI >= surf.
nEdges()))
1755 <<
"Edge label outside valid range." <<
endl
1756 <<
"edge label:" << edgeI <<
endl
1757 <<
"total number of edges:" << surf.
nEdges() <<
endl
1761 const labelList& neighbours = edgeFaces[edgeI];
1763 if (neighbours.
size() == 2)
1765 const label stat0 = faceStatus[neighbours[0]];
1766 const label stat1 = faceStatus[neighbours[1]];
1771 ((stat0 == ANYEDGE) || (stat0 == edgeI))
1772 && ((stat1 == ANYEDGE) || (stat1 == edgeI))
1775 const edge& e = edges[edgeI];
1781 || (pointMap[e.
end()] != e.
end())
1785 <<
"points already mapped. Double collapse." <<
endl
1786 <<
"edgeI:" << edgeI
1787 <<
" start:" << e.
start()
1788 <<
" end:" << e.
end()
1789 <<
" pointMap[start]:" << pointMap[e.
start()]
1790 <<
" pointMap[end]:" << pointMap[e.
end()]
1795 pointMap[e.
start()] = minVert;
1796 pointMap[e.
end()] = minVert;
1799 newPoints[minVert] = edgeMids[edgeI];
1802 protectNeighbours(surf, e.
start(), faceStatus);
1803 protectNeighbours(surf, e.
end(), faceStatus);
1807 oppositeVertex(surf, neighbours[0], edgeI),
1813 oppositeVertex(surf, neighbours[1], edgeI),
1825 forAll(collapseFaces, collapseI)
1827 faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1842 forAll(localFaces, faceI)
1846 const label a = pointMap[f[0]];
1847 const label
b = pointMap[f[1]];
1848 const label c = pointMap[f[2]];
1852 (a != b) && (a != c) && (b != c)
1853 && (faceStatus[faceI] != COLLAPSED)
1857 newTris[newTriI++] =
labelledTri(a, b, c, f.region());
1876 tempSurf.localFaces(),
1878 tempSurf.localPoints()
1894 forAll(refineFaces, refineFaceI)
1896 calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus);
1900 return doRefine(surf, refineStatus);
1918 label newPointI = surf.
nPoints();
1922 forAll(refineEdges, refineEdgeI)
1924 label edgeI = refineEdges[refineEdgeI];
1928 bool neighbourIsRefined=
false;
1932 if (refineStatus[myFaces[myFaceI]] != NONE)
1934 neighbourIsRefined =
true;
1939 if (!neighbourIsRefined)
1951 newPoints[newPointI] = mid;
1967 refineStatus[myFaces[myFaceI]] = GREEN;
1977 if (refineStatus[faceI] == NONE)
2002 scalar minLength = GREAT;
2003 label minIndex = -1;
2006 const edge& e = surf.
edges()[edgeIndices[i]];
2015 if (length < minLength)
2021 return edgeIndices[minIndex];
2032 scalar maxLength = -GREAT;
2033 label maxIndex = -1;
2036 const edge& e = surf.
edges()[edgeIndices[i]];
2045 if (length > maxLength)
2051 return edgeIndices[maxIndex];
2059 const scalar mergeTol
2081 label newTriangleI = 0;
2087 label newA = pointMap[f[0]];
2088 label newB = pointMap[f[1]];
2089 label newC = pointMap[f[2]];
2091 if ((newA != newB) && (newA != newC) && (newB != newC))
2093 newTriangles[newTriangleI++] =
2097 newTriangles.
setSize(newTriangleI);
2117 const label nearestFaceI,
2118 const point& nearestPt
2131 ).classify(nearestPt, 1
E-6, nearType, nearLabel);
2141 label edgeI = surf.
faceEdges()[nearestFaceI][nearLabel];
2152 return edgeNormal/(
mag(edgeNormal) + VSMALL);
2166 const point& sample,
2167 const point& nearestPoint,
2173 if (eFaces.
size() != 2)
2185 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2187 if (((sample - nearestPoint) & n) > 0)
2203 const point& sample,
2204 const label nearestFaceI,
2205 const point& nearestPoint,
2214 label nearType, nearLabel;
2220 ).classify(nearestPoint, tol, nearType, nearLabel);
2225 scalar c = (sample - nearestPoint) & surf.
faceNormals()[nearestFaceI];
2242 label edgeI = surf.
faceEdges()[nearestFaceI][nearLabel];
2266 return edgeSide(surf, sample, nearestPoint, edgeI);
2276 label nearPointI = localF[nearLabel];
2280 const point& base = localPoints[nearPointI];
2285 label minEdgeI = -1;
2289 label edgeI = pEdges[i];
2291 const edge& e = edges[edgeI];
2296 vector eVec(localPoints[otherPointI] - base);
2297 scalar magEVec =
mag(eVec);
2299 if (magEVec > VSMALL)
2304 const point perturbPoint = base + eVec;
2308 if (distSqr < minDistSqr)
2310 minDistSqr = distSqr;
2319 <<
"Problem: did not find edge closer than " << minDistSqr
2323 return edgeSide(surf, sample, nearestPoint, minEdgeI);
2344 label newPatchI = 0;
2349 iter != includePatches.
end();
2353 label patchI = iter.key();
2359 label nTriTotal = 0;
2361 forAll(patch, patchFaceI)
2363 const face& f = patch[patchFaceI];
2371 forAll(triFaces, triFaceI)
2373 const face& f = triFaces[triFaceI];
2383 Pout<< patch.
name() <<
" : generated " << nTriTotal
2384 <<
" triangles from " << patch.size() <<
" faces with"
2385 <<
" new patchid " << newPatchI <<
endl;
2398 rawSurface.localFaces(),
2399 rawSurface.localPoints()
2410 iter != includePatches.
end();
2414 label patchI = iter.key();
2418 surface.patches()[newPatchI].
name() = patch.
name();
2420 surface.patches()[newPatchI].geometricType() = patch.type();
2445 label newPointI = 0;
2449 newPoints[newPointI++] = points[pointI];
2451 forAll(faceCentres, faceI)
2453 newPoints[newPointI++] = faceCentres[faceI];
2463 label newPatchI = 0;
2468 iter != includePatches.
end();
2472 label patchI = iter.key();
2476 label nTriTotal = 0;
2478 forAll(patch, patchFaceI)
2481 const face& f = patch[patchFaceI];
2484 label fc = points.
size() + patchFaceI + patch.
start();
2488 label fp1 = (fp + 1) % f.
size();
2490 triangles.append(
labelledTri(f[fp], f[fp1], fc, newPatchI));
2498 Pout<< patch.
name() <<
" : generated " << nTriTotal
2499 <<
" triangles from " << patch.size() <<
" faces with"
2500 <<
" new patchid " << newPatchI <<
endl;
2526 iter != includePatches.
end();
2530 label patchI = iter.key();
2534 surface.patches()[newPatchI].
name() = patch.
name();
2536 surface.patches()[newPatchI].geometricType() = patch.type();
2553 geompackVertices[doubleI++] = pts[i][0];
2554 geompackVertices[doubleI++] = pts[i][1];
2567 geompackVertices.begin(),
2569 triangle_node.
begin(),
2570 triangle_neighbor.begin()
2574 triangle_node.setSize(3*nTris);
2575 triangle_neighbor.setSize(3*nTris);
2584 triangle_node[3*i]-1,
2585 triangle_node[3*i+1]-1,
2586 triangle_node[3*i+2]-1,
2594 points[i][0] = pts[i][0];
2595 points[i][1] = pts[i][1];
2691 edge[0] = tri.
c()-tri.
b();
2692 edge[1] = tri.
a()-tri.
c();
2693 edge[2] = tri.
b()-tri.
a();
2695 vector triangleFaceNormal = edge[1] ^ edge[2];
2699 for(label i=0; i<3; i++)
2701 normal[i] = triangleFaceNormal ^ edge[i];
2702 normal[i] /=
mag(normal[i]) + VSMALL;
2705 weights[0] = ((p-tri.
b()) & normal[0]) /
max(VSMALL, normal[0] & edge[1]);
2706 weights[1] = ((p-tri.
c()) & normal[1]) /
max(VSMALL, normal[1] & edge[2]);
2707 weights[2] = ((p-tri.
a()) & normal[2]) /
max(VSMALL, normal[2] & edge[0]);
2721 allVerts.setSize(samplePts.
size());
2722 allWeights.setSize(samplePts.
size());
2728 const point& samplePt = samplePts[i];
2734 scalar minDistance = GREAT;
2742 pointHit nearest = tri.nearestPoint(samplePt);
2751 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2761 else if (nearest.distance() < minDistance)
2763 minDistance = nearest.distance();
2766 label nearType, nearLabel;
2777 verts[0] = f[nearLabel];
2780 weights[1] = -GREAT;
2782 weights[2] = -GREAT;
2791 verts[0] = f[nearLabel];
2792 verts[1] = f[f.
fcIndex(nearLabel)];
2795 const point& p0 = points[verts[0]];
2796 const point& p1 = points[verts[1]];
2804 mag(nearest.rawPoint()-p0)/
mag(p1-p0)
2811 weights[2] = -GREAT;
2825 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2844 const point& trianglePoint
2850 label index, elemType;
2890 const plane& cutPlane
2898 snapToEnd(s, end, nearest);
2931 nearest = visitFaces
2946 nearest = visitFaces
2957 snapToEnd(s, end, nearest);
2967 const plane& cutPlane,
2983 hitInfo = trackToEdge