1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 \*---------------------------------------------------------------------------*/
29 #include "triSurfaceTools.H"
31 #include "triSurface.H"
33 #include "mergePoints.H"
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
42 const Foam::label Foam::triSurfaceTools::NOEDGE = -2;
43 const Foam::label Foam::triSurfaceTools::COLLAPSED = -3;
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 Refine by splitting all three edges of triangle ('red' refinement).
49 Neighbouring triangles (which are not marked for refinement get split
50 in half ('green') refinement. (R. Verfuerth, "A review of a posteriori
51 error estimation and adaptive mesh refinement techniques",
55 // FaceI gets refined ('red'). Propagate edge refinement.
56 void Foam::triSurfaceTools::calcRefineStatus
58 const triSurface& surf,
60 List<refineType>& refine
63 if (refine[faceI] == RED)
65 // Already marked for refinement. Do nothing.
69 // Not marked or marked for 'green' refinement. Refine.
72 const labelList& myNeighbours = surf.faceFaces()[faceI];
74 forAll(myNeighbours, myNeighbourI)
76 label neighbourFaceI = myNeighbours[myNeighbourI];
78 if (refine[neighbourFaceI] == GREEN)
80 // Change to red refinement and propagate
81 calcRefineStatus(surf, neighbourFaceI, refine);
83 else if (refine[neighbourFaceI] == NONE)
85 refine[neighbourFaceI] = GREEN;
92 // Split faceI along edgeI at position newPointI
93 void Foam::triSurfaceTools::greenRefine
95 const triSurface& surf,
98 const label newPointI,
99 DynamicList<labelledTri>& newFaces
102 const labelledTri& f = surf.localFaces()[faceI];
103 const edge& e = surf.edges()[edgeI];
105 // Find index of edge in face.
107 label fp0 = findIndex(f, e[0]);
108 label fp1 = f.fcIndex(fp0);
109 label fp2 = f.fcIndex(fp1);
113 // Edge oriented like face
161 // Refine all triangles marked for refinement.
162 Foam::triSurface Foam::triSurfaceTools::doRefine
164 const triSurface& surf,
165 const List<refineType>& refineStatus
168 // Storage for new points. (start after old points)
169 DynamicList<point> newPoints(surf.nPoints());
170 forAll(surf.localPoints(), pointI)
172 newPoints.append(surf.localPoints()[pointI]);
174 label newVertI = surf.nPoints();
176 // Storage for new faces
177 DynamicList<labelledTri> newFaces(surf.size());
180 // Point index for midpoint on edge
181 labelList edgeMid(surf.nEdges(), -1);
183 forAll(refineStatus, faceI)
185 if (refineStatus[faceI] == RED)
187 // Create new vertices on all edges to be refined.
188 const labelList& fEdges = surf.faceEdges()[faceI];
192 label edgeI = fEdges[i];
194 if (edgeMid[edgeI] == -1)
196 const edge& e = surf.edges()[edgeI];
198 // Create new point on mid of edge
203 surf.localPoints()[e.start()]
204 + surf.localPoints()[e.end()]
207 edgeMid[edgeI] = newVertI++;
211 // Now we have new mid edge vertices for all edges on face
212 // so create triangles for RED rerfinement.
214 const edgeList& edges = surf.edges();
222 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
233 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
244 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
263 // Create triangles for GREEN refinement.
266 const label edgeI = fEdges[i];
268 label otherFaceI = otherFace(surf, faceI, edgeI);
270 if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN))
285 // Copy unmarked triangles since keep original vertex numbering.
286 forAll(refineStatus, faceI)
288 if (refineStatus[faceI] == NONE)
290 newFaces.append(surf.localFaces()[faceI]);
298 // Transfer DynamicLists to straight ones.
299 pointField allPoints;
300 allPoints.transfer(newPoints);
303 return triSurface(newFaces, surf.patches(), allPoints, true);
307 // Angle between two neighbouring triangles,
308 // angle between shared-edge end points and left and right face end points
309 Foam::scalar Foam::triSurfaceTools::faceCosAngle
317 const vector common(pEnd - pStart);
318 const vector base0(pLeft - pStart);
319 const vector base1(pRight - pStart);
321 vector n0(common ^ base0);
324 vector n1(base1 ^ common);
331 // Protect edges around vertex from collapsing.
332 // Moving vertI to a different position
333 // - affects obviously the cost of the faces using it
334 // - but also their neighbours since the edge between the faces changes
335 void Foam::triSurfaceTools::protectNeighbours
337 const triSurface& surf,
339 labelList& faceStatus
342 // const labelList& myFaces = surf.pointFaces()[vertI];
343 // forAll(myFaces, i)
345 // label faceI = myFaces[i];
347 // if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
349 // faceStatus[faceI] = NOEDGE;
353 const labelList& myEdges = surf.pointEdges()[vertI];
356 const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
358 forAll(myFaces, myFaceI)
360 label faceI = myFaces[myFaceI];
362 if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
364 faceStatus[faceI] = NOEDGE;
372 // Edge collapse helper functions
375 // Get all faces that will get collapsed if edgeI collapses.
376 Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
378 const triSurface& surf,
382 const edge& e = surf.edges()[edgeI];
383 label v1 = e.start();
386 // Faces using edge will certainly get collapsed.
387 const labelList& myFaces = surf.edgeFaces()[edgeI];
389 labelHashSet facesToBeCollapsed(2*myFaces.size());
391 forAll(myFaces, myFaceI)
393 facesToBeCollapsed.insert(myFaces[myFaceI]);
396 // From faces using v1 check if they share an edge with faces
398 // - share edge: are part of 'splay' tree and will collapse if edge
400 const labelList& v1Faces = surf.pointFaces()[v1];
402 forAll(v1Faces, v1FaceI)
404 label face1I = v1Faces[v1FaceI];
406 label otherEdgeI = oppositeEdge(surf, face1I, v1);
408 // Step across edge to other face
409 label face2I = otherFace(surf, face1I, otherEdgeI);
413 // found face on other side of edge. Now check if uses v2.
414 if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
416 // triangles face1I and face2I form splay tree and will
418 facesToBeCollapsed.insert(face1I);
419 facesToBeCollapsed.insert(face2I);
424 return facesToBeCollapsed;
428 // Return value of faceUsed for faces using vertI
429 Foam::label Foam::triSurfaceTools::vertexUsesFace
431 const triSurface& surf,
432 const labelHashSet& faceUsed,
436 const labelList& myFaces = surf.pointFaces()[vertI];
438 forAll(myFaces, myFaceI)
440 label face1I = myFaces[myFaceI];
442 if (faceUsed.found(face1I))
451 // Calculate new edge-face addressing resulting from edge collapse
452 void Foam::triSurfaceTools::getMergedEdges
454 const triSurface& surf,
456 const labelHashSet& collapsedFaces,
457 HashTable<label, label, Hash<label> >& edgeToEdge,
458 HashTable<label, label, Hash<label> >& edgeToFace
461 const edge& e = surf.edges()[edgeI];
462 label v1 = e.start();
465 const labelList& v1Faces = surf.pointFaces()[v1];
466 const labelList& v2Faces = surf.pointFaces()[v2];
468 // Mark all (non collapsed) faces using v2
469 labelHashSet v2FacesHash(v2Faces.size());
471 forAll(v2Faces, v2FaceI)
473 if (!collapsedFaces.found(v2Faces[v2FaceI]))
475 v2FacesHash.insert(v2Faces[v2FaceI]);
480 forAll(v1Faces, v1FaceI)
482 label face1I = v1Faces[v1FaceI];
484 if (collapsedFaces.found(face1I))
490 // Check if face1I uses a vertex connected to a face using v2
503 //Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:"
504 // << vert1I << ' ' << vert2I << endl;
506 // Check vert1, vert2 for usage by v2Face.
508 label commonVert = vert1I;
509 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
513 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
518 // Found one: commonVert is used by both face1I and face2I
519 label edge1I = getEdge(surf, v1, commonVert);
520 label edge2I = getEdge(surf, v2, commonVert);
522 edgeToEdge.insert(edge1I, edge2I);
523 edgeToEdge.insert(edge2I, edge1I);
525 edgeToFace.insert(edge1I, face2I);
526 edgeToFace.insert(edge2I, face1I);
532 // Calculates (cos of) angle across edgeI of faceI,
533 // taking into account updated addressing (resulting from edge collapse)
534 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
536 const triSurface& surf,
539 const labelHashSet& collapsedFaces,
540 const HashTable<label, label, Hash<label> >& edgeToEdge,
541 const HashTable<label, label, Hash<label> >& edgeToFace,
546 const pointField& localPoints = surf.localPoints();
548 label A = surf.edges()[edgeI].start();
549 label B = surf.edges()[edgeI].end();
550 label C = oppositeVertex(surf, faceI, edgeI);
556 if (edgeToEdge.found(edgeI))
558 // Use merged addressing
559 label edge2I = edgeToEdge[edgeI];
560 face2I = edgeToFace[edgeI];
562 D = oppositeVertex(surf, face2I, edge2I);
566 // Use normal edge-face addressing
567 face2I = otherFace(surf, faceI, edgeI);
569 if ((face2I != -1) && !collapsedFaces.found(face2I))
571 D = oppositeVertex(surf, face2I, edgeI);
581 cosAngle = faceCosAngle
591 cosAngle = faceCosAngle
601 cosAngle = faceCosAngle
611 cosAngle = faceCosAngle
621 FatalErrorIn("edgeCosAngle")
622 << "face " << faceI << " does not use vertex "
623 << v1 << " of collapsed edge" << abort(FatalError);
630 //- Calculate minimum (cos of) edge angle using addressing from collapsing
632 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
634 const triSurface& surf,
637 const labelHashSet& collapsedFaces,
638 const HashTable<label, label, Hash<label> >& edgeToEdge,
639 const HashTable<label, label, Hash<label> >& edgeToFace
642 const labelList& v1Faces = surf.pointFaces()[v1];
646 forAll(v1Faces, v1FaceI)
648 label faceI = v1Faces[v1FaceI];
650 if (collapsedFaces.found(faceI))
655 const labelList& myEdges = surf.faceEdges()[faceI];
657 forAll(myEdges, myEdgeI)
659 label edgeI = myEdges[myEdgeI];
684 // Calculate max dihedral angle after collapsing edge to v1 (at pt).
685 // Note that all edges of all faces using v1 are affected.
686 bool Foam::triSurfaceTools::collapseCreatesFold
688 const triSurface& surf,
691 const labelHashSet& collapsedFaces,
692 const HashTable<label, label, Hash<label> >& edgeToEdge,
693 const HashTable<label, label, Hash<label> >& edgeToFace,
697 const labelList& v1Faces = surf.pointFaces()[v1];
699 forAll(v1Faces, v1FaceI)
701 label faceI = v1Faces[v1FaceI];
703 if (collapsedFaces.found(faceI))
708 const labelList& myEdges = surf.faceEdges()[faceI];
710 forAll(myEdges, myEdgeI)
712 label edgeI = myEdges[myEdgeI];
739 //// Return true if collapsing edgeI creates triangles on top of each other.
740 //// Is when the triangles neighbouring the collapsed one already share
742 //bool Foam::triSurfaceTools::collapseCreatesDuplicates
744 // const triSurface& surf,
745 // const label edgeI,
746 // const labelHashSet& collapsedFaces
749 //Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI]
750 // << " collapsedFaces:" << collapsedFaces.toc() << endl;
752 // // Mark neighbours of faces to be collapsed, i.e. get the first layer
753 // // of triangles outside collapsedFaces.
754 // // neighbours actually contains the
755 // // edge with which triangle connects to collapsedFaces.
757 // HashTable<label, label, Hash<label> > neighbours;
759 // labelList collapsed = collapsedFaces.toc();
761 // forAll(collapsed, collapseI)
763 // const label faceI = collapsed[collapseI];
765 // const labelList& myEdges = surf.faceEdges()[faceI];
767 // Pout<< "collapsing faceI:" << faceI << " uses edges:" << myEdges
770 // forAll(myEdges, myEdgeI)
772 // const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]];
774 // Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces "
775 // << myFaces << endl;
777 // if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2))
779 // // Get the other face
780 // label neighbourFaceI = myFaces[0];
782 // if (neighbourFaceI == faceI)
784 // neighbourFaceI = myFaces[1];
787 // // Is 'outside' face if not in collapsedFaces itself
788 // if (!collapsedFaces.found(neighbourFaceI))
790 // neighbours.insert(neighbourFaceI, myEdges[myEdgeI]);
796 // // Now neighbours contains first layer of triangles outside of
798 // // There should be
799 // // -two if edgeI is a boundary edge
800 // // since the outside 'edge' of collapseFaces should
801 // // form a triangle and the face connected to edgeI is not inserted.
802 // // -four if edgeI is not a boundary edge since then the outside edge forms
805 // // Check if any of the faces in neighbours share an edge. (n^2)
807 // labelList neighbourList = neighbours.toc();
809 //Pout<< "edgeI:" << edgeI << " neighbourList:" << neighbourList << endl;
812 // forAll(neighbourList, i)
814 // const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]];
816 // for (label j = i+1; j < neighbourList.size(); i++)
818 // const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
820 // // Check if faceI and faceJ share an edge
821 // forAll(faceIEdges, fI)
823 // forAll(faceJEdges, fJ)
825 // Pout<< " comparing " << faceIEdges[fI] << " to "
826 // << faceJEdges[fJ] << endl;
828 // if (faceIEdges[fI] == faceJEdges[fJ])
836 // Pout<< "Found no match. Returning false" << endl;
841 // Finds the triangle edge cut by the plane between a point inside / on edge
842 // of a triangle and a point outside. Returns:
843 // - cut through edge/point
845 Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
849 const label excludeEdgeI,
850 const label excludePointI,
852 const point& triPoint,
853 const plane& cutPlane,
857 const pointField& points = s.points();
858 const labelledTri& f = s[triI];
859 const labelList& fEdges = s.faceEdges()[triI];
861 // Get normal distance to planeN
862 FixedList<scalar, 3> d;
867 d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
871 // Normalise and truncate
876 if (mag(d[i]) < 1E-6)
882 // Return information
885 if (excludePointI != -1)
887 // Excluded point. Test only opposite edge.
889 label fp0 = findIndex(s.localFaces()[triI], excludePointI);
893 FatalErrorIn("cutEdge(..)") << "excludePointI:" << excludePointI
894 << " localF:" << s.localFaces()[triI] << abort(FatalError);
897 label fp1 = f.fcIndex(fp0);
898 label fp2 = f.fcIndex(fp1);
904 cut.setPoint(points[f[fp1]]);
905 cut.elementType() = triPointRef::POINT;
906 cut.setIndex(s.localFaces()[triI][fp1]);
908 else if (d[fp2] == 0.0)
911 cut.setPoint(points[f[fp2]]);
912 cut.elementType() = triPointRef::POINT;
913 cut.setIndex(s.localFaces()[triI][fp2]);
917 (d[fp1] < 0 && d[fp2] < 0)
918 || (d[fp1] > 0 && d[fp2] > 0)
921 // Both same sign. Not crossing edge at all.
922 // cut already set to miss().
929 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
932 cut.elementType() = triPointRef::EDGE;
933 cut.setIndex(fEdges[fp1]);
938 // Find the two intersections
939 FixedList<surfaceLocation, 2> inters;
944 label fp1 = f.fcIndex(fp0);
950 FatalErrorIn("cutEdge(..)")
951 << "problem : triangle has three intersections." << nl
952 << "triangle:" << f.tri(points)
953 << " d:" << d << abort(FatalError);
955 inters[interI].setHit();
956 inters[interI].setPoint(points[f[fp0]]);
957 inters[interI].elementType() = triPointRef::POINT;
958 inters[interI].setIndex(s.localFaces()[triI][fp0]);
963 (d[fp0] < 0 && d[fp1] > 0)
964 || (d[fp0] > 0 && d[fp1] < 0)
969 FatalErrorIn("cutEdge(..)")
970 << "problem : triangle has three intersections." << nl
971 << "triangle:" << f.tri(points)
972 << " d:" << d << abort(FatalError);
974 inters[interI].setHit();
975 inters[interI].setPoint
977 (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
980 inters[interI].elementType() = triPointRef::EDGE;
981 inters[interI].setIndex(fEdges[fp0]);
991 else if (interI == 1)
993 // Only one intersection. Should not happen!
996 else if (interI == 2)
998 // Handle excludeEdgeI
1001 inters[0].elementType() == triPointRef::EDGE
1002 && inters[0].index() == excludeEdgeI
1009 inters[1].elementType() == triPointRef::EDGE
1010 && inters[1].index() == excludeEdgeI
1017 // Two cuts. Find nearest.
1020 magSqr(inters[0].rawPoint() - toPoint)
1021 < magSqr(inters[1].rawPoint() - toPoint)
1037 // 'Snap' point on to endPoint.
1038 void Foam::triSurfaceTools::snapToEnd
1040 const triSurface& s,
1041 const surfaceLocation& end,
1042 surfaceLocation& current
1045 if (end.elementType() == triPointRef::NONE)
1047 if (current.elementType() == triPointRef::NONE)
1049 // endpoint on triangle; current on triangle
1050 if (current.index() == end.index())
1054 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1061 // No need to handle current on edge/point since tracking handles this.
1063 else if (end.elementType() == triPointRef::EDGE)
1065 if (current.elementType() == triPointRef::NONE)
1067 // endpoint on edge; current on triangle
1068 const labelList& fEdges = s.faceEdges()[current.index()];
1070 if (findIndex(fEdges, end.index()) != -1)
1074 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1081 else if (current.elementType() == triPointRef::EDGE)
1083 // endpoint on edge; current on edge
1084 if (current.index() == end.index())
1088 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1097 // endpoint on edge; current on point
1098 const edge& e = s.edges()[end.index()];
1100 if (current.index() == e[0] || current.index() == e[1])
1104 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1112 else // end.elementType() == POINT
1114 if (current.elementType() == triPointRef::NONE)
1116 // endpoint on point; current on triangle
1117 const labelledTri& f = s.localFaces()[current.index()];
1119 if (findIndex(f, end.index()) != -1)
1123 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1130 else if (current.elementType() == triPointRef::EDGE)
1132 // endpoint on point; current on edge
1133 const edge& e = s.edges()[current.index()];
1135 if (end.index() == e[0] || end.index() == e[1])
1139 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1148 // endpoint on point; current on point
1149 if (current.index() == end.index())
1153 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1166 // - element type (triangle/edge/point) and index
1167 // - triangle to exclude
1168 Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
1170 const triSurface& s,
1171 const labelList& eFaces,
1172 const surfaceLocation& start,
1173 const label excludeEdgeI,
1174 const label excludePointI,
1175 const surfaceLocation& end,
1176 const plane& cutPlane
1179 surfaceLocation nearest;
1181 scalar minDistSqr = Foam::sqr(GREAT);
1185 label triI = eFaces[i];
1187 // Make sure we don't revisit previous face
1188 if (triI != start.triangle())
1190 if (end.elementType() == triPointRef::NONE && end.index() == triI)
1192 // Endpoint is in this triangle. Jump there.
1195 nearest.triangle() = triI;
1200 // Which edge is cut.
1202 surfaceLocation cutInfo = cutEdge
1206 excludeEdgeI, // excludeEdgeI
1207 excludePointI, // excludePointI
1213 // If crossing an edge we expect next edge to be cut.
1214 if (excludeEdgeI != -1 && !cutInfo.hit())
1216 FatalErrorIn("triSurfaceTools::visitFaces(..)")
1217 << "Triangle:" << triI
1218 << " excludeEdge:" << excludeEdgeI
1219 << " point:" << start.rawPoint()
1220 << " plane:" << cutPlane
1221 << " . No intersection!" << abort(FatalError);
1226 scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
1228 if (distSqr < minDistSqr)
1230 minDistSqr = distSqr;
1232 nearest.triangle() = triI;
1240 if (nearest.triangle() == -1)
1242 // Did not move from edge. Give warning? Return something special?
1243 // For now responsability of caller to make sure that nothing has
1251 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1254 // Write pointField to file
1255 void Foam::triSurfaceTools::writeOBJ
1257 const fileName& fName,
1258 const pointField& pts
1261 OFstream outFile(fName);
1265 const point& pt = pts[pointI];
1267 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1269 Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
1273 // Write vertex subset to OBJ format file
1274 void Foam::triSurfaceTools::writeOBJ
1276 const triSurface& surf,
1277 const fileName& fName,
1278 const boolList& markedVerts
1281 OFstream outFile(fName);
1284 forAll(markedVerts, vertI)
1286 if (markedVerts[vertI])
1288 const point& pt = surf.localPoints()[vertI];
1290 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1295 Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
1299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1300 // Addressing helper functions:
1303 // Get all triangles using vertices of edge
1304 void Foam::triSurfaceTools::getVertexTriangles
1306 const triSurface& surf,
1311 const edge& e = surf.edges()[edgeI];
1312 const labelList& myFaces = surf.edgeFaces()[edgeI];
1314 label face1I = myFaces[0];
1316 if (myFaces.size() == 2)
1318 face2I = myFaces[1];
1321 const labelList& startFaces = surf.pointFaces()[e.start()];
1322 const labelList& endFaces = surf.pointFaces()[e.end()];
1324 // Number of triangles is sum of pointfaces - common faces
1325 // (= faces using edge)
1326 edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
1329 forAll(startFaces, startFaceI)
1331 edgeTris[nTris++] = startFaces[startFaceI];
1334 forAll(endFaces, endFaceI)
1336 label faceI = endFaces[endFaceI];
1338 if ((faceI != face1I) && (faceI != face2I))
1340 edgeTris[nTris++] = faceI;
1346 // Get all vertices connected to vertices of edge
1347 Foam::labelList Foam::triSurfaceTools::getVertexVertices
1349 const triSurface& surf,
1353 const edgeList& edges = surf.edges();
1355 label v1 = e.start();
1358 // Get all vertices connected to v1 or v2 through an edge
1359 labelHashSet vertexNeighbours;
1361 const labelList& v1Edges = surf.pointEdges()[v1];
1363 forAll(v1Edges, v1EdgeI)
1365 const edge& e = edges[v1Edges[v1EdgeI]];
1366 vertexNeighbours.insert(e.otherVertex(v1));
1369 const labelList& v2Edges = surf.pointEdges()[v2];
1371 forAll(v2Edges, v2EdgeI)
1373 const edge& e = edges[v2Edges[v2EdgeI]];
1375 label vertI = e.otherVertex(v2);
1377 vertexNeighbours.insert(vertI);
1379 return vertexNeighbours.toc();
1383 //// Order vertices consistent with face
1384 //void Foam::triSurfaceTools::orderVertices
1386 // const labelledTri& f,
1393 // // Order v1, v2 in anticlockwise order.
1394 // bool reverse = false;
1403 // else if (f[1] == v1)
1431 // Get the other face using edgeI
1432 Foam::label Foam::triSurfaceTools::otherFace
1434 const triSurface& surf,
1439 const labelList& myFaces = surf.edgeFaces()[edgeI];
1441 if (myFaces.size() != 2)
1447 if (faceI == myFaces[0])
1459 // Get the two edges on faceI counterclockwise after edgeI
1460 void Foam::triSurfaceTools::otherEdges
1462 const triSurface& surf,
1469 const labelList& eFaces = surf.faceEdges()[faceI];
1471 label i0 = findIndex(eFaces, edgeI);
1478 "(const triSurface&, const label, const label,"
1480 ) << "Edge " << surf.edges()[edgeI] << " not in face "
1481 << surf.localFaces()[faceI] << abort(FatalError);
1484 label i1 = eFaces.fcIndex(i0);
1485 label i2 = eFaces.fcIndex(i1);
1492 // Get the two vertices on faceI counterclockwise vertI
1493 void Foam::triSurfaceTools::otherVertices
1495 const triSurface& surf,
1502 const labelledTri& f = surf.localFaces()[faceI];
1509 else if (vertI == f[1])
1514 else if (vertI == f[2])
1524 "(const triSurface&, const label, const label,"
1526 ) << "Vertex " << vertI << " not in face " << f << abort(FatalError);
1531 // Get edge opposite vertex
1532 Foam::label Foam::triSurfaceTools::oppositeEdge
1534 const triSurface& surf,
1539 const labelList& myEdges = surf.faceEdges()[faceI];
1541 forAll(myEdges, myEdgeI)
1543 label edgeI = myEdges[myEdgeI];
1545 const edge& e = surf.edges()[edgeI];
1547 if ((e.start() != vertI) && (e.end() != vertI))
1556 "(const triSurface&, const label, const label)"
1557 ) << "Cannot find vertex " << vertI << " in edges of face " << faceI
1558 << abort(FatalError);
1564 // Get vertex opposite edge
1565 Foam::label Foam::triSurfaceTools::oppositeVertex
1567 const triSurface& surf,
1572 const labelledTri& f = surf.localFaces()[faceI];
1574 const edge& e = surf.edges()[edgeI];
1578 label vertI = f[fp];
1580 if (vertI != e.start() && vertI != e.end())
1586 FatalErrorIn("triSurfaceTools::oppositeVertex")
1587 << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
1588 << " in face " << faceI << " vertices " << f << abort(FatalError);
1594 // Returns edge label connecting v1, v2
1595 Foam::label Foam::triSurfaceTools::getEdge
1597 const triSurface& surf,
1602 const labelList& v1Edges = surf.pointEdges()[v1];
1604 forAll(v1Edges, v1EdgeI)
1606 label edgeI = v1Edges[v1EdgeI];
1608 const edge& e = surf.edges()[edgeI];
1610 if ((e.start() == v2) || (e.end() == v2))
1619 // Return index of triangle (or -1) using all three edges
1620 Foam::label Foam::triSurfaceTools::getTriangle
1622 const triSurface& surf,
1628 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1633 "(const triSurface&, const label, const label,"
1635 ) << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
1637 << abort(FatalError);
1640 const labelList& eFaces = surf.edgeFaces()[e0I];
1642 forAll(eFaces, eFaceI)
1644 label faceI = eFaces[eFaceI];
1646 const labelList& myEdges = surf.faceEdges()[faceI];
1651 || (myEdges[1] == e1I)
1652 || (myEdges[2] == e1I)
1658 || (myEdges[1] == e2I)
1659 || (myEdges[2] == e2I)
1670 // Collapse indicated edges. Return new tri surface.
1671 Foam::triSurface Foam::triSurfaceTools::collapseEdges
1673 const triSurface& surf,
1674 const labelList& collapsableEdges
1677 pointField edgeMids(surf.nEdges());
1679 forAll(edgeMids, edgeI)
1681 const edge& e = surf.edges()[edgeI];
1686 surf.localPoints()[e.start()]
1687 + surf.localPoints()[e.end()]
1692 labelList faceStatus(surf.size(), ANYEDGE);
1694 //// Protect triangles which are on the border of different regions
1695 //forAll(edges, edgeI)
1697 // const labelList& neighbours = edgeFaces[edgeI];
1699 // if ((neighbours.size() != 2) && (neighbours.size() != 1))
1701 // FatalErrorIn("collapseEdges")
1702 // << abort(FatalError);
1705 // if (neighbours.size() == 2)
1707 // if (surf[neighbours[0]].region() != surf[neighbours[1]].region())
1709 // // Neighbours on different regions. For now dont allow
1711 // //Pout<< "protecting face " << neighbours[0]
1712 // // << ' ' << neighbours[1] << endl;
1713 // faceStatus[neighbours[0]] = NOEDGE;
1714 // faceStatus[neighbours[1]] = NOEDGE;
1719 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1723 // Collapse indicated edges. Return new tri surface.
1724 Foam::triSurface Foam::triSurfaceTools::collapseEdges
1726 const triSurface& surf,
1727 const labelList& collapseEdgeLabels,
1728 const pointField& edgeMids,
1729 labelList& faceStatus
1732 const labelListList& edgeFaces = surf.edgeFaces();
1733 const pointField& localPoints = surf.localPoints();
1734 const edgeList& edges = surf.edges();
1736 // Storage for new points
1737 pointField newPoints(localPoints);
1739 // Map for old to new points
1740 labelList pointMap(localPoints.size());
1741 forAll(localPoints, pointI)
1743 pointMap[pointI] = pointI;
1747 // Do actual 'collapsing' of edges
1749 forAll(collapseEdgeLabels, collapseEdgeI)
1751 const label edgeI = collapseEdgeLabels[collapseEdgeI];
1753 if ((edgeI < 0) || (edgeI >= surf.nEdges()))
1755 FatalErrorIn("collapseEdges")
1756 << "Edge label outside valid range." << endl
1757 << "edge label:" << edgeI << endl
1758 << "total number of edges:" << surf.nEdges() << endl
1759 << abort(FatalError);
1762 const labelList& neighbours = edgeFaces[edgeI];
1764 if (neighbours.size() == 2)
1766 const label stat0 = faceStatus[neighbours[0]];
1767 const label stat1 = faceStatus[neighbours[1]];
1769 // Check faceStatus to make sure this one can be collapsed
1772 ((stat0 == ANYEDGE) || (stat0 == edgeI))
1773 && ((stat1 == ANYEDGE) || (stat1 == edgeI))
1776 const edge& e = edges[edgeI];
1778 // Set up mapping to 'collapse' points of edge
1781 (pointMap[e.start()] != e.start())
1782 || (pointMap[e.end()] != e.end())
1785 FatalErrorIn("collapseEdges")
1786 << "points already mapped. Double collapse." << endl
1787 << "edgeI:" << edgeI
1788 << " start:" << e.start()
1789 << " end:" << e.end()
1790 << " pointMap[start]:" << pointMap[e.start()]
1791 << " pointMap[end]:" << pointMap[e.end()]
1792 << abort(FatalError);
1795 const label minVert = min(e.start(), e.end());
1796 pointMap[e.start()] = minVert;
1797 pointMap[e.end()] = minVert;
1799 // Move shared vertex to mid of edge
1800 newPoints[minVert] = edgeMids[edgeI];
1802 // Protect neighbouring faces
1803 protectNeighbours(surf, e.start(), faceStatus);
1804 protectNeighbours(surf, e.end(), faceStatus);
1808 oppositeVertex(surf, neighbours[0], edgeI),
1814 oppositeVertex(surf, neighbours[1], edgeI),
1818 // Mark all collapsing faces
1819 labelList collapseFaces =
1826 forAll(collapseFaces, collapseI)
1828 faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1835 // Storage for new triangles
1836 List<labelledTri> newTris(surf.size());
1839 const List<labelledTri>& localFaces = surf.localFaces();
1842 // Get only non-collapsed triangles and renumber vertex labels.
1843 forAll(localFaces, faceI)
1845 const labelledTri& f = localFaces[faceI];
1847 const label a = pointMap[f[0]];
1848 const label b = pointMap[f[1]];
1849 const label c = pointMap[f[2]];
1853 (a != b) && (a != c) && (b != c)
1854 && (faceStatus[faceI] != COLLAPSED)
1857 // uncollapsed triangle
1858 newTris[newTriI++] = labelledTri(a, b, c, f.region());
1862 //Pout<< "Collapsed triangle " << faceI
1863 // << " vertices:" << f << endl;
1866 newTris.setSize(newTriI);
1872 triSurface tempSurf(newTris, surf.patches(), newPoints);
1877 tempSurf.localFaces(),
1879 tempSurf.localPoints()
1884 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1886 Foam::triSurface Foam::triSurfaceTools::redGreenRefine
1888 const triSurface& surf,
1889 const labelList& refineFaces
1892 List<refineType> refineStatus(surf.size(), NONE);
1894 // Mark & propagate refinement
1895 forAll(refineFaces, refineFaceI)
1897 calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus);
1900 // Do actual refinement
1901 return doRefine(surf, refineStatus);
1905 Foam::triSurface Foam::triSurfaceTools::greenRefine
1907 const triSurface& surf,
1908 const labelList& refineEdges
1911 // Storage for marking faces
1912 List<refineType> refineStatus(surf.size(), NONE);
1914 // Storage for new faces
1915 DynamicList<labelledTri> newFaces(0);
1917 pointField newPoints(surf.localPoints());
1918 newPoints.setSize(surf.nPoints() + surf.nEdges());
1919 label newPointI = surf.nPoints();
1923 forAll(refineEdges, refineEdgeI)
1925 label edgeI = refineEdges[refineEdgeI];
1927 const labelList& myFaces = surf.edgeFaces()[edgeI];
1929 bool neighbourIsRefined= false;
1931 forAll(myFaces, myFaceI)
1933 if (refineStatus[myFaces[myFaceI]] != NONE)
1935 neighbourIsRefined = true;
1939 // Only refine if none of the faces is refined
1940 if (!neighbourIsRefined)
1943 const edge& e = surf.edges()[edgeI];
1948 surf.localPoints()[e.start()]
1949 + surf.localPoints()[e.end()]
1952 newPoints[newPointI] = mid;
1954 // Refine faces using edge
1955 forAll(myFaces, myFaceI)
1957 // Add faces to newFaces
1968 refineStatus[myFaces[myFaceI]] = GREEN;
1975 // Add unrefined faces
1976 forAll(surf.localFaces(), faceI)
1978 if (refineStatus[faceI] == NONE)
1980 newFaces.append(surf.localFaces()[faceI]);
1985 newPoints.setSize(newPointI);
1987 return triSurface(newFaces, surf.patches(), newPoints, true);
1992 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1993 // Geometric helper functions:
1996 // Returns element in edgeIndices with minimum length
1997 Foam::label Foam::triSurfaceTools::minEdge
1999 const triSurface& surf,
2000 const labelList& edgeIndices
2003 scalar minLength = GREAT;
2004 label minIndex = -1;
2005 forAll(edgeIndices, i)
2007 const edge& e = surf.edges()[edgeIndices[i]];
2012 surf.localPoints()[e.end()]
2013 - surf.localPoints()[e.start()]
2016 if (length < minLength)
2022 return edgeIndices[minIndex];
2026 // Returns element in edgeIndices with maximum length
2027 Foam::label Foam::triSurfaceTools::maxEdge
2029 const triSurface& surf,
2030 const labelList& edgeIndices
2033 scalar maxLength = -GREAT;
2034 label maxIndex = -1;
2035 forAll(edgeIndices, i)
2037 const edge& e = surf.edges()[edgeIndices[i]];
2042 surf.localPoints()[e.end()]
2043 - surf.localPoints()[e.start()]
2046 if (length > maxLength)
2052 return edgeIndices[maxIndex];
2056 // Merge points and reconstruct surface
2057 Foam::triSurface Foam::triSurfaceTools::mergePoints
2059 const triSurface& surf,
2060 const scalar mergeTol
2063 pointField newPoints(surf.nPoints());
2065 labelList pointMap(surf.nPoints());
2067 bool hasMerged = Foam::mergePoints
2078 // Pack the triangles
2080 // Storage for new triangles
2081 List<labelledTri> newTriangles(surf.size());
2082 label newTriangleI = 0;
2086 const labelledTri& f = surf.localFaces()[faceI];
2088 label newA = pointMap[f[0]];
2089 label newB = pointMap[f[1]];
2090 label newC = pointMap[f[2]];
2092 if ((newA != newB) && (newA != newC) && (newB != newC))
2094 newTriangles[newTriangleI++] =
2095 labelledTri(newA, newB, newC, f.region());
2098 newTriangles.setSize(newTriangleI);
2114 // Calculate normal on triangle
2115 Foam::vector Foam::triSurfaceTools::surfaceNormal
2117 const triSurface& surf,
2118 const label nearestFaceI,
2119 const point& nearestPt
2122 const labelledTri& f = surf[nearestFaceI];
2123 const pointField& points = surf.points();
2132 ).classify(nearestPt, 1E-6, nearType, nearLabel);
2134 if (nearType == triPointRef::NONE)
2137 return surf.faceNormals()[nearestFaceI];
2139 else if (nearType == triPointRef::EDGE)
2141 // Nearest to edge. Assume order of faceEdges same as face vertices.
2142 label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
2144 // Calculate edge normal by averaging face normals
2145 const labelList& eFaces = surf.edgeFaces()[edgeI];
2147 vector edgeNormal(vector::zero);
2151 edgeNormal += surf.faceNormals()[eFaces[i]];
2153 return edgeNormal/(mag(edgeNormal) + VSMALL);
2158 const labelledTri& localF = surf.localFaces()[nearestFaceI];
2159 return surf.pointNormals()[localF[nearLabel]];
2164 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::edgeSide
2166 const triSurface& surf,
2167 const point& sample,
2168 const point& nearestPoint,
2172 const labelList& eFaces = surf.edgeFaces()[edgeI];
2174 if (eFaces.size() != 2)
2176 // Surface not closed.
2181 const vectorField& faceNormals = surf.faceNormals();
2183 // Compare to bisector. This is actually correct since edge is
2184 // nearest so there is a knife-edge.
2186 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2188 if (((sample - nearestPoint) & n) > 0)
2200 // Calculate normal on triangle
2201 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
2203 const triSurface& surf,
2204 const point& sample,
2205 const label nearestFaceI, // nearest face
2206 const point& nearestPoint, // nearest point on nearest face
2210 const labelledTri& f = surf[nearestFaceI];
2211 const pointField& points = surf.points();
2213 // Find where point is on triangle. Note tolerance needed. Is relative
2214 // one (used in comparing normalized [0..1] triangle coordinates).
2215 label nearType, nearLabel;
2221 ).classify(nearestPoint, tol, nearType, nearLabel);
2223 if (nearType == triPointRef::NONE)
2225 // Nearest to face interior. Use faceNormal to determine side
2226 scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
2237 else if (nearType == triPointRef::EDGE)
2239 // Nearest to edge nearLabel. Note that this can only be a knife-edge
2240 // situation since otherwise the nearest point could never be the edge.
2242 // Get the edge. Assume order of faceEdges same as face vertices.
2243 label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
2247 // // Check order of faceEdges same as face vertices.
2248 // const edge& e = surf.edges()[edgeI];
2249 // const labelList& meshPoints = surf.meshPoints();
2250 // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
2255 // != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
2258 // FatalErrorIn("triSurfaceTools::surfaceSide")
2259 // << "Edge:" << edgeI << " local vertices:" << e
2260 // << " mesh vertices:" << meshEdge
2261 // << " not at position " << nearLabel
2262 // << " in face " << f
2263 // << abort(FatalError);
2267 return edgeSide(surf, sample, nearestPoint, edgeI);
2271 // Nearest to point. Could use pointNormal here but is not correct.
2272 // Instead determine which edge using point is nearest and use test
2273 // above (nearType == triPointRef::EDGE).
2276 const labelledTri& localF = surf.localFaces()[nearestFaceI];
2277 label nearPointI = localF[nearLabel];
2279 const edgeList& edges = surf.edges();
2280 const pointField& localPoints = surf.localPoints();
2281 const point& base = localPoints[nearPointI];
2283 const labelList& pEdges = surf.pointEdges()[nearPointI];
2285 scalar minDistSqr = Foam::sqr(GREAT);
2286 label minEdgeI = -1;
2290 label edgeI = pEdges[i];
2292 const edge& e = edges[edgeI];
2294 label otherPointI = e.otherVertex(nearPointI);
2297 vector eVec(localPoints[otherPointI] - base);
2298 scalar magEVec = mag(eVec);
2300 if (magEVec > VSMALL)
2304 // Get point along vector and determine closest.
2305 const point perturbPoint = base + eVec;
2307 scalar distSqr = Foam::magSqr(sample - perturbPoint);
2309 if (distSqr < minDistSqr)
2311 minDistSqr = distSqr;
2319 FatalErrorIn("treeDataTriSurface::getSide")
2320 << "Problem: did not find edge closer than " << minDistSqr
2321 << abort(FatalError);
2324 return edgeSide(surf, sample, nearestPoint, minEdgeI);
2329 // triangulation of boundaryMesh
2330 Foam::triSurface Foam::triSurfaceTools::triangulate
2332 const polyBoundaryMesh& bMesh,
2333 const labelHashSet& includePatches,
2337 const polyMesh& mesh = bMesh.mesh();
2339 // Storage for surfaceMesh. Size estimate.
2340 DynamicList<labelledTri> triangles
2342 mesh.nFaces() - mesh.nInternalFaces()
2345 label newPatchI = 0;
2349 labelHashSet::const_iterator iter = includePatches.begin();
2350 iter != includePatches.end();
2354 label patchI = iter.key();
2356 const polyPatch& patch = bMesh[patchI];
2358 const pointField& points = patch.points();
2360 label nTriTotal = 0;
2362 forAll(patch, patchFaceI)
2364 const face& f = patch[patchFaceI];
2366 faceList triFaces(f.nTriangles(points));
2370 f.triangles(points, nTri, triFaces);
2372 forAll(triFaces, triFaceI)
2374 const face& f = triFaces[triFaceI];
2376 triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
2384 Pout<< patch.name() << " : generated " << nTriTotal
2385 << " triangles from " << patch.size() << " faces with"
2386 << " new patchid " << newPatchI << endl;
2393 // Create globally numbered tri surface
2394 triSurface rawSurface(triangles, mesh.points());
2396 // Create locally numbered tri surface
2399 rawSurface.localFaces(),
2400 rawSurface.localPoints()
2403 // Add patch names to surface
2404 surface.patches().setSize(newPatchI);
2410 labelHashSet::const_iterator iter = includePatches.begin();
2411 iter != includePatches.end();
2415 label patchI = iter.key();
2417 const polyPatch& patch = bMesh[patchI];
2419 surface.patches()[newPatchI].name() = patch.name();
2421 surface.patches()[newPatchI].geometricType() = patch.type();
2430 // triangulation of boundaryMesh
2431 Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
2433 const polyBoundaryMesh& bMesh,
2434 const labelHashSet& includePatches,
2438 const polyMesh& mesh = bMesh.mesh();
2440 // Storage for new points = meshpoints + face centres.
2441 const pointField& points = mesh.points();
2442 const pointField& faceCentres = mesh.faceCentres();
2444 pointField newPoints(points.size() + faceCentres.size());
2446 label newPointI = 0;
2448 forAll(points, pointI)
2450 newPoints[newPointI++] = points[pointI];
2452 forAll(faceCentres, faceI)
2454 newPoints[newPointI++] = faceCentres[faceI];
2458 // Count number of faces.
2459 DynamicList<labelledTri> triangles
2461 mesh.nFaces() - mesh.nInternalFaces()
2464 label newPatchI = 0;
2468 labelHashSet::const_iterator iter = includePatches.begin();
2469 iter != includePatches.end();
2473 label patchI = iter.key();
2475 const polyPatch& patch = bMesh[patchI];
2477 label nTriTotal = 0;
2479 forAll(patch, patchFaceI)
2481 // Face in global coords.
2482 const face& f = patch[patchFaceI];
2484 // Index in newPointI of face centre.
2485 label fc = points.size() + patchFaceI + patch.start();
2489 label fp1 = (fp + 1) % f.size();
2491 triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchI));
2499 Pout<< patch.name() << " : generated " << nTriTotal
2500 << " triangles from " << patch.size() << " faces with"
2501 << " new patchid " << newPatchI << endl;
2509 // Create globally numbered tri surface
2510 triSurface rawSurface(triangles, newPoints);
2512 // Create locally numbered tri surface
2515 rawSurface.localFaces(),
2516 rawSurface.localPoints()
2519 // Add patch names to surface
2520 surface.patches().setSize(newPatchI);
2526 labelHashSet::const_iterator iter = includePatches.begin();
2527 iter != includePatches.end();
2531 label patchI = iter.key();
2533 const polyPatch& patch = bMesh[patchI];
2535 surface.patches()[newPatchI].name() = patch.name();
2537 surface.patches()[newPatchI].geometricType() = patch.type();
2546 Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
2548 // Vertices in geompack notation. Note that could probably just use
2549 // pts.begin() if double precision.
2550 List<doubleScalar> geompackVertices(2*pts.size());
2554 geompackVertices[doubleI++] = pts[i][0];
2555 geompackVertices[doubleI++] = pts[i][1];
2558 // Storage for triangles
2560 List<int> triangle_node(m2*3*pts.size());
2561 List<int> triangle_neighbor(m2*3*pts.size());
2568 geompackVertices.begin(),
2570 triangle_node.begin(),
2571 triangle_neighbor.begin()
2575 triangle_node.setSize(3*nTris);
2576 triangle_neighbor.setSize(3*nTris);
2578 // Convert to triSurface.
2579 List<labelledTri> faces(nTris);
2583 faces[i] = labelledTri
2585 triangle_node[3*i]-1,
2586 triangle_node[3*i+1]-1,
2587 triangle_node[3*i+2]-1,
2592 pointField points(pts.size());
2595 points[i][0] = pts[i][0];
2596 points[i][1] = pts[i][1];
2600 return triSurface(faces, points);
2604 //// Use CGAL to do Delaunay
2605 //#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
2606 //#include <CGAL/Delaunay_triangulation_2.h>
2607 //#include <CGAL/Triangulation_vertex_base_with_info_2.h>
2608 //#include <cassert>
2610 //struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
2612 //typedef CGAL::Triangulation_vertex_base_with_info_2<Foam::label, K> Vb;
2613 //typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
2614 //typedef CGAL::Delaunay_triangulation_2<K,Tds> Triangulation;
2616 //typedef Triangulation::Vertex_handle Vertex_handle;
2617 //typedef Triangulation::Vertex_iterator Vertex_iterator;
2618 //typedef Triangulation::Face_handle Face_handle;
2619 //typedef Triangulation::Finite_faces_iterator Finite_faces_iterator;
2620 //typedef Triangulation::Point Point;
2621 //Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
2625 // // Insert 2D vertices; building triangulation
2628 // const point& pt = pts[i];
2630 // T.insert(Point(pt[0], pt[1]));
2634 // // Number vertices
2635 // // ~~~~~~~~~~~~~~~
2640 // Vertex_iterator it = T.vertices_begin();
2641 // it != T.vertices_end();
2645 // it->info() = vertI++;
2651 // List<labelledTri> faces(T.number_of_faces());
2657 // Finite_faces_iterator fc = T.finite_faces_begin();
2658 // fc != T.finite_faces_end();
2662 // faces[faceI++] = labelledTri
2664 // fc->vertex(0)->info(),
2665 // f[1] = fc->vertex(1)->info(),
2666 // f[2] = fc->vertex(2)->info()
2670 // pointField points(pts.size());
2673 // points[i][0] = pts[i][0];
2674 // points[i][1] = pts[i][1];
2675 // points[i][2] = 0.0;
2678 // return triSurface(faces, points);
2682 void Foam::triSurfaceTools::calcInterpolationWeights
2684 const triPointRef& tri,
2686 FixedList<scalar, 3>& weights
2689 // calculate triangle edge vectors and triangle face normal
2690 // the 'i':th edge is opposite node i
2691 FixedList<vector, 3> edge;
2692 edge[0] = tri.c()-tri.b();
2693 edge[1] = tri.a()-tri.c();
2694 edge[2] = tri.b()-tri.a();
2696 vector triangleFaceNormal = edge[1] ^ edge[2];
2698 // calculate edge normal (pointing inwards)
2699 FixedList<vector, 3> normal;
2700 for(label i=0; i<3; i++)
2702 normal[i] = triangleFaceNormal ^ edge[i];
2703 normal[i] /= mag(normal[i]) + VSMALL;
2706 weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
2707 weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
2708 weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
2712 // Calculate weighting factors from samplePts to triangle it is in.
2713 // Uses linear search.
2714 void Foam::triSurfaceTools::calcInterpolationWeights
2716 const triSurface& s,
2717 const pointField& samplePts,
2718 List<FixedList<label, 3> >& allVerts,
2719 List<FixedList<scalar, 3> >& allWeights
2722 allVerts.setSize(samplePts.size());
2723 allWeights.setSize(samplePts.size());
2725 const pointField& points = s.points();
2727 forAll(samplePts, i)
2729 const point& samplePt = samplePts[i];
2732 FixedList<label, 3>& verts = allVerts[i];
2733 FixedList<scalar, 3>& weights = allWeights[i];
2735 scalar minDistance = GREAT;
2739 const labelledTri& f = s[faceI];
2741 triPointRef tri(f.tri(points));
2743 pointHit nearest = tri.nearestPoint(samplePt);
2747 // samplePt inside triangle
2752 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2754 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2755 // << " inside triangle:" << faceI
2756 // << " verts:" << verts
2757 // << " weights:" << weights
2762 else if (nearest.distance() < minDistance)
2764 minDistance = nearest.distance();
2766 // Outside triangle. Store nearest.
2767 label nearType, nearLabel;
2771 1E-6, // relative tolerance
2776 if (nearType == triPointRef::POINT)
2778 verts[0] = f[nearLabel];
2781 weights[1] = -GREAT;
2783 weights[2] = -GREAT;
2785 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2786 // << " distance:" << nearest.distance()
2787 // << " from point:" << points[f[nearLabel]]
2790 else if (nearType == triPointRef::EDGE)
2792 verts[0] = f[nearLabel];
2793 verts[1] = f[f.fcIndex(nearLabel)];
2796 const point& p0 = points[verts[0]];
2797 const point& p1 = points[verts[1]];
2805 mag(nearest.rawPoint()-p0)/mag(p1-p0)
2812 weights[2] = -GREAT;
2814 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2815 // << " distance:" << nearest.distance()
2816 // << " from edge:" << p0 << p1 << " s:" << s
2821 // triangle. Can only happen because of truncation errors.
2826 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2836 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2840 // Test point on surface to see if is on face,edge or point.
2841 Foam::surfaceLocation Foam::triSurfaceTools::classify
2843 const triSurface& s,
2845 const point& trianglePoint
2848 surfaceLocation nearest;
2850 // Nearest point could be on point or edge. Retest.
2851 label index, elemType;
2853 triPointRef(s[triI].tri(s.points())).classify
2861 nearest.setPoint(trianglePoint);
2863 if (elemType == triPointRef::NONE)
2866 nearest.setIndex(triI);
2867 nearest.elementType() = triPointRef::NONE;
2869 else if (elemType == triPointRef::EDGE)
2872 nearest.setIndex(s.faceEdges()[triI][index]);
2873 nearest.elementType() = triPointRef::EDGE;
2875 else //if (elemType == triPointRef::POINT)
2878 nearest.setIndex(s.localFaces()[triI][index]);
2879 nearest.elementType() = triPointRef::POINT;
2886 Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge
2888 const triSurface& s,
2889 const surfaceLocation& start,
2890 const surfaceLocation& end,
2891 const plane& cutPlane
2894 // Start off from starting point
2895 surfaceLocation nearest = start;
2898 // See if in same triangle as endpoint. If so snap.
2899 snapToEnd(s, end, nearest);
2903 // Not yet at end point
2905 if (start.elementType() == triPointRef::NONE)
2907 // Start point is inside triangle. Trivial cases already handled
2910 // end point is on edge or point so cross currrent triangle to
2911 // see which edge is cut.
2916 start.index(), // triangle
2923 nearest.elementType() = triPointRef::EDGE;
2924 nearest.triangle() = start.index();
2927 else if (start.elementType() == triPointRef::EDGE)
2929 // Pick connected triangle that is most in direction.
2930 const labelList& eFaces = s.edgeFaces()[start.index()];
2932 nearest = visitFaces
2937 start.index(), // excludeEdgeI
2938 -1, // excludePointI
2943 else // start.elementType() == triPointRef::POINT
2945 const labelList& pFaces = s.pointFaces()[start.index()];
2947 nearest = visitFaces
2953 start.index(), // excludePointI
2958 snapToEnd(s, end, nearest);
2964 void Foam::triSurfaceTools::track
2966 const triSurface& s,
2967 const surfaceLocation& endInfo,
2968 const plane& cutPlane,
2969 surfaceLocation& hitInfo
2972 //OFstream str("track.obj");
2974 //meshTools::writeOBJ(str, hitInfo.rawPoint());
2977 // Track across surface.
2980 //Pout<< "Tracking from:" << nl
2981 // << " " << hitInfo.info()
2984 hitInfo = trackToEdge
2992 //meshTools::writeOBJ(str, hitInfo.rawPoint());
2994 //str<< "l " << vertI-1 << ' ' << vertI << nl;
2996 //Pout<< "Tracked to:" << nl
2997 // << " " << hitInfo.info() << endl;
2999 if (hitInfo.hit() || hitInfo.triangle() == -1)
3007 // ************************************************************************* //