/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . \*---------------------------------------------------------------------------*/ #include "triSurfaceTools.H" #include "triSurface.H" #include "OFstream.H" #include "mergePoints.H" #include "polyMesh.H" #include "plane.H" #include "geompack.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // const Foam::label Foam::triSurfaceTools::ANYEDGE = -1; const Foam::label Foam::triSurfaceTools::NOEDGE = -2; const Foam::label Foam::triSurfaceTools::COLLAPSED = -3; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // /* Refine by splitting all three edges of triangle ('red' refinement). Neighbouring triangles (which are not marked for refinement get split in half ('green') refinement. (R. Verfuerth, "A review of a posteriori error estimation and adaptive mesh refinement techniques", Wiley-Teubner, 1996) */ // FaceI gets refined ('red'). Propagate edge refinement. void Foam::triSurfaceTools::calcRefineStatus ( const triSurface& surf, const label faceI, List& refine ) { if (refine[faceI] == RED) { // Already marked for refinement. Do nothing. } else { // Not marked or marked for 'green' refinement. Refine. refine[faceI] = RED; const labelList& myNeighbours = surf.faceFaces()[faceI]; forAll(myNeighbours, myNeighbourI) { label neighbourFaceI = myNeighbours[myNeighbourI]; if (refine[neighbourFaceI] == GREEN) { // Change to red refinement and propagate calcRefineStatus(surf, neighbourFaceI, refine); } else if (refine[neighbourFaceI] == NONE) { refine[neighbourFaceI] = GREEN; } } } } // Split faceI along edgeI at position newPointI void Foam::triSurfaceTools::greenRefine ( const triSurface& surf, const label faceI, const label edgeI, const label newPointI, DynamicList& newFaces ) { const labelledTri& f = surf.localFaces()[faceI]; const edge& e = surf.edges()[edgeI]; // Find index of edge in face. label fp0 = findIndex(f, e[0]); label fp1 = f.fcIndex(fp0); label fp2 = f.fcIndex(fp1); if (f[fp1] == e[1]) { // Edge oriented like face newFaces.append ( labelledTri ( f[fp0], newPointI, f[fp2], f.region() ) ); newFaces.append ( labelledTri ( newPointI, f[fp1], f[fp2], f.region() ) ); } else { newFaces.append ( labelledTri ( f[fp2], newPointI, f[fp1], f.region() ) ); newFaces.append ( labelledTri ( newPointI, f[fp0], f[fp1], f.region() ) ); } } // Refine all triangles marked for refinement. Foam::triSurface Foam::triSurfaceTools::doRefine ( const triSurface& surf, const List& refineStatus ) { // Storage for new points. (start after old points) DynamicList newPoints(surf.nPoints()); forAll(surf.localPoints(), pointI) { newPoints.append(surf.localPoints()[pointI]); } label newVertI = surf.nPoints(); // Storage for new faces DynamicList newFaces(surf.size()); // Point index for midpoint on edge labelList edgeMid(surf.nEdges(), -1); forAll(refineStatus, faceI) { if (refineStatus[faceI] == RED) { // Create new vertices on all edges to be refined. const labelList& fEdges = surf.faceEdges()[faceI]; forAll(fEdges, i) { label edgeI = fEdges[i]; if (edgeMid[edgeI] == -1) { const edge& e = surf.edges()[edgeI]; // Create new point on mid of edge newPoints.append ( 0.5 * ( surf.localPoints()[e.start()] + surf.localPoints()[e.end()] ) ); edgeMid[edgeI] = newVertI++; } } // Now we have new mid edge vertices for all edges on face // so create triangles for RED rerfinement. const edgeList& edges = surf.edges(); // Corner triangles newFaces.append ( labelledTri ( edgeMid[fEdges[0]], edges[fEdges[0]].commonVertex(edges[fEdges[1]]), edgeMid[fEdges[1]], surf[faceI].region() ) ); newFaces.append ( labelledTri ( edgeMid[fEdges[1]], edges[fEdges[1]].commonVertex(edges[fEdges[2]]), edgeMid[fEdges[2]], surf[faceI].region() ) ); newFaces.append ( labelledTri ( edgeMid[fEdges[2]], edges[fEdges[2]].commonVertex(edges[fEdges[0]]), edgeMid[fEdges[0]], surf[faceI].region() ) ); // Inner triangle newFaces.append ( labelledTri ( edgeMid[fEdges[0]], edgeMid[fEdges[1]], edgeMid[fEdges[2]], surf[faceI].region() ) ); // Create triangles for GREEN refinement. forAll(fEdges, i) { const label edgeI = fEdges[i]; label otherFaceI = otherFace(surf, faceI, edgeI); if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN)) { greenRefine ( surf, otherFaceI, edgeI, edgeMid[edgeI], newFaces ); } } } } // Copy unmarked triangles since keep original vertex numbering. forAll(refineStatus, faceI) { if (refineStatus[faceI] == NONE) { newFaces.append(surf.localFaces()[faceI]); } } newFaces.shrink(); newPoints.shrink(); // Transfer DynamicLists to straight ones. pointField allPoints; allPoints.transfer(newPoints); newPoints.clear(); return triSurface(newFaces, surf.patches(), allPoints, true); } // Angle between two neighbouring triangles, // angle between shared-edge end points and left and right face end points Foam::scalar Foam::triSurfaceTools::faceCosAngle ( const point& pStart, const point& pEnd, const point& pLeft, const point& pRight ) { const vector common(pEnd - pStart); const vector base0(pLeft - pStart); const vector base1(pRight - pStart); vector n0(common ^ base0); n0 /= Foam::mag(n0); vector n1(base1 ^ common); n1 /= Foam::mag(n1); return n0 & n1; } // Protect edges around vertex from collapsing. // Moving vertI to a different position // - affects obviously the cost of the faces using it // - but also their neighbours since the edge between the faces changes void Foam::triSurfaceTools::protectNeighbours ( const triSurface& surf, const label vertI, labelList& faceStatus ) { // const labelList& myFaces = surf.pointFaces()[vertI]; // forAll(myFaces, i) // { // label faceI = myFaces[i]; // // if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0)) // { // faceStatus[faceI] = NOEDGE; // } // } const labelList& myEdges = surf.pointEdges()[vertI]; forAll(myEdges, i) { const labelList& myFaces = surf.edgeFaces()[myEdges[i]]; forAll(myFaces, myFaceI) { label faceI = myFaces[myFaceI]; if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0)) { faceStatus[faceI] = NOEDGE; } } } } // // Edge collapse helper functions // // Get all faces that will get collapsed if edgeI collapses. Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces ( const triSurface& surf, label edgeI ) { const edge& e = surf.edges()[edgeI]; label v1 = e.start(); label v2 = e.end(); // Faces using edge will certainly get collapsed. const labelList& myFaces = surf.edgeFaces()[edgeI]; labelHashSet facesToBeCollapsed(2*myFaces.size()); forAll(myFaces, myFaceI) { facesToBeCollapsed.insert(myFaces[myFaceI]); } // From faces using v1 check if they share an edge with faces // using v2. // - share edge: are part of 'splay' tree and will collapse if edge // collapses const labelList& v1Faces = surf.pointFaces()[v1]; forAll(v1Faces, v1FaceI) { label face1I = v1Faces[v1FaceI]; label otherEdgeI = oppositeEdge(surf, face1I, v1); // Step across edge to other face label face2I = otherFace(surf, face1I, otherEdgeI); if (face2I != -1) { // found face on other side of edge. Now check if uses v2. if (oppositeVertex(surf, face2I, otherEdgeI) == v2) { // triangles face1I and face2I form splay tree and will // collapse. facesToBeCollapsed.insert(face1I); facesToBeCollapsed.insert(face2I); } } } return facesToBeCollapsed; } // Return value of faceUsed for faces using vertI Foam::label Foam::triSurfaceTools::vertexUsesFace ( const triSurface& surf, const labelHashSet& faceUsed, const label vertI ) { const labelList& myFaces = surf.pointFaces()[vertI]; forAll(myFaces, myFaceI) { label face1I = myFaces[myFaceI]; if (faceUsed.found(face1I)) { return face1I; } } return -1; } // Calculate new edge-face addressing resulting from edge collapse void Foam::triSurfaceTools::getMergedEdges ( const triSurface& surf, const label edgeI, const labelHashSet& collapsedFaces, HashTable >& edgeToEdge, HashTable >& edgeToFace ) { const edge& e = surf.edges()[edgeI]; label v1 = e.start(); label v2 = e.end(); const labelList& v1Faces = surf.pointFaces()[v1]; const labelList& v2Faces = surf.pointFaces()[v2]; // Mark all (non collapsed) faces using v2 labelHashSet v2FacesHash(v2Faces.size()); forAll(v2Faces, v2FaceI) { if (!collapsedFaces.found(v2Faces[v2FaceI])) { v2FacesHash.insert(v2Faces[v2FaceI]); } } forAll(v1Faces, v1FaceI) { label face1I = v1Faces[v1FaceI]; if (collapsedFaces.found(face1I)) { continue; } // // Check if face1I uses a vertex connected to a face using v2 // label vert1I = -1; label vert2I = -1; otherVertices ( surf, face1I, v1, vert1I, vert2I ); //Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:" // << vert1I << ' ' << vert2I << endl; // Check vert1, vert2 for usage by v2Face. label commonVert = vert1I; label face2I = vertexUsesFace(surf, v2FacesHash, commonVert); if (face2I == -1) { commonVert = vert2I; face2I = vertexUsesFace(surf, v2FacesHash, commonVert); } if (face2I != -1) { // Found one: commonVert is used by both face1I and face2I label edge1I = getEdge(surf, v1, commonVert); label edge2I = getEdge(surf, v2, commonVert); edgeToEdge.insert(edge1I, edge2I); edgeToEdge.insert(edge2I, edge1I); edgeToFace.insert(edge1I, face2I); edgeToFace.insert(edge2I, face1I); } } } // Calculates (cos of) angle across edgeI of faceI, // taking into account updated addressing (resulting from edge collapse) Foam::scalar Foam::triSurfaceTools::edgeCosAngle ( const triSurface& surf, const label v1, const point& pt, const labelHashSet& collapsedFaces, const HashTable >& edgeToEdge, const HashTable >& edgeToFace, const label faceI, const label edgeI ) { const pointField& localPoints = surf.localPoints(); label A = surf.edges()[edgeI].start(); label B = surf.edges()[edgeI].end(); label C = oppositeVertex(surf, faceI, edgeI); label D = -1; label face2I = -1; if (edgeToEdge.found(edgeI)) { // Use merged addressing label edge2I = edgeToEdge[edgeI]; face2I = edgeToFace[edgeI]; D = oppositeVertex(surf, face2I, edge2I); } else { // Use normal edge-face addressing face2I = otherFace(surf, faceI, edgeI); if ((face2I != -1) && !collapsedFaces.found(face2I)) { D = oppositeVertex(surf, face2I, edgeI); } } scalar cosAngle = 1; if (D != -1) { if (A == v1) { cosAngle = faceCosAngle ( pt, localPoints[B], localPoints[C], localPoints[D] ); } else if (B == v1) { cosAngle = faceCosAngle ( localPoints[A], pt, localPoints[C], localPoints[D] ); } else if (C == v1) { cosAngle = faceCosAngle ( localPoints[A], localPoints[B], pt, localPoints[D] ); } else if (D == v1) { cosAngle = faceCosAngle ( localPoints[A], localPoints[B], localPoints[C], pt ); } else { FatalErrorIn("edgeCosAngle") << "face " << faceI << " does not use vertex " << v1 << " of collapsed edge" << abort(FatalError); } } return cosAngle; } //- Calculate minimum (cos of) edge angle using addressing from collapsing // edge to v1 at pt. Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle ( const triSurface& surf, const label v1, const point& pt, const labelHashSet& collapsedFaces, const HashTable >& edgeToEdge, const HashTable >& edgeToFace ) { const labelList& v1Faces = surf.pointFaces()[v1]; scalar minCos = 1; forAll(v1Faces, v1FaceI) { label faceI = v1Faces[v1FaceI]; if (collapsedFaces.found(faceI)) { continue; } const labelList& myEdges = surf.faceEdges()[faceI]; forAll(myEdges, myEdgeI) { label edgeI = myEdges[myEdgeI]; minCos = min ( minCos, edgeCosAngle ( surf, v1, pt, collapsedFaces, edgeToEdge, edgeToFace, faceI, edgeI ) ); } } return minCos; } // Calculate max dihedral angle after collapsing edge to v1 (at pt). // Note that all edges of all faces using v1 are affected. bool Foam::triSurfaceTools::collapseCreatesFold ( const triSurface& surf, const label v1, const point& pt, const labelHashSet& collapsedFaces, const HashTable >& edgeToEdge, const HashTable >& edgeToFace, const scalar minCos ) { const labelList& v1Faces = surf.pointFaces()[v1]; forAll(v1Faces, v1FaceI) { label faceI = v1Faces[v1FaceI]; if (collapsedFaces.found(faceI)) { continue; } const labelList& myEdges = surf.faceEdges()[faceI]; forAll(myEdges, myEdgeI) { label edgeI = myEdges[myEdgeI]; if ( edgeCosAngle ( surf, v1, pt, collapsedFaces, edgeToEdge, edgeToFace, faceI, edgeI ) < minCos ) { return true; } } } return false; } //// Return true if collapsing edgeI creates triangles on top of each other. //// Is when the triangles neighbouring the collapsed one already share // a vertex. //bool Foam::triSurfaceTools::collapseCreatesDuplicates //( // const triSurface& surf, // const label edgeI, // const labelHashSet& collapsedFaces //) //{ //Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI] // << " collapsedFaces:" << collapsedFaces.toc() << endl; // // // Mark neighbours of faces to be collapsed, i.e. get the first layer // // of triangles outside collapsedFaces. // // neighbours actually contains the // // edge with which triangle connects to collapsedFaces. // // HashTable > neighbours; // // labelList collapsed = collapsedFaces.toc(); // // forAll(collapsed, collapseI) // { // const label faceI = collapsed[collapseI]; // // const labelList& myEdges = surf.faceEdges()[faceI]; // // Pout<< "collapsing faceI:" << faceI << " uses edges:" << myEdges // << endl; // // forAll(myEdges, myEdgeI) // { // const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]]; // // Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces " // << myFaces << endl; // // if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2)) // { // // Get the other face // label neighbourFaceI = myFaces[0]; // // if (neighbourFaceI == faceI) // { // neighbourFaceI = myFaces[1]; // } // // // Is 'outside' face if not in collapsedFaces itself // if (!collapsedFaces.found(neighbourFaceI)) // { // neighbours.insert(neighbourFaceI, myEdges[myEdgeI]); // } // } // } // } // // // Now neighbours contains first layer of triangles outside of // // collapseFaces // // There should be // // -two if edgeI is a boundary edge // // since the outside 'edge' of collapseFaces should // // form a triangle and the face connected to edgeI is not inserted. // // -four if edgeI is not a boundary edge since then the outside edge forms // // a diamond. // // // Check if any of the faces in neighbours share an edge. (n^2) // // labelList neighbourList = neighbours.toc(); // //Pout<< "edgeI:" << edgeI << " neighbourList:" << neighbourList << endl; // // // forAll(neighbourList, i) // { // const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]]; // // for (label j = i+1; j < neighbourList.size(); i++) // { // const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]]; // // // Check if faceI and faceJ share an edge // forAll(faceIEdges, fI) // { // forAll(faceJEdges, fJ) // { // Pout<< " comparing " << faceIEdges[fI] << " to " // << faceJEdges[fJ] << endl; // // if (faceIEdges[fI] == faceJEdges[fJ]) // { // return true; // } // } // } // } // } // Pout<< "Found no match. Returning false" << endl; // return false; //} // Finds the triangle edge cut by the plane between a point inside / on edge // of a triangle and a point outside. Returns: // - cut through edge/point // - miss() Foam::surfaceLocation Foam::triSurfaceTools::cutEdge ( const triSurface& s, const label triI, const label excludeEdgeI, const label excludePointI, const point& triPoint, const plane& cutPlane, const point& toPoint ) { const pointField& points = s.points(); const labelledTri& f = s[triI]; const labelList& fEdges = s.faceEdges()[triI]; // Get normal distance to planeN FixedList d; scalar norm = 0; forAll(d, fp) { d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal(); norm += mag(d[fp]); } // Normalise and truncate forAll(d, i) { d[i] /= norm; if (mag(d[i]) < 1E-6) { d[i] = 0.0; } } // Return information surfaceLocation cut; if (excludePointI != -1) { // Excluded point. Test only opposite edge. label fp0 = findIndex(s.localFaces()[triI], excludePointI); if (fp0 == -1) { FatalErrorIn("cutEdge(..)") << "excludePointI:" << excludePointI << " localF:" << s.localFaces()[triI] << abort(FatalError); } label fp1 = f.fcIndex(fp0); label fp2 = f.fcIndex(fp1); if (d[fp1] == 0.0) { cut.setHit(); cut.setPoint(points[f[fp1]]); cut.elementType() = triPointRef::POINT; cut.setIndex(s.localFaces()[triI][fp1]); } else if (d[fp2] == 0.0) { cut.setHit(); cut.setPoint(points[f[fp2]]); cut.elementType() = triPointRef::POINT; cut.setIndex(s.localFaces()[triI][fp2]); } else if ( (d[fp1] < 0 && d[fp2] < 0) || (d[fp1] > 0 && d[fp2] > 0) ) { // Both same sign. Not crossing edge at all. // cut already set to miss(). } else { cut.setHit(); cut.setPoint ( (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]]) / (d[fp2] - d[fp1]) ); cut.elementType() = triPointRef::EDGE; cut.setIndex(fEdges[fp1]); } } else { // Find the two intersections FixedList inters; label interI = 0; forAll(f, fp0) { label fp1 = f.fcIndex(fp0); if (d[fp0] == 0) { if (interI >= 2) { FatalErrorIn("cutEdge(..)") << "problem : triangle has three intersections." << nl << "triangle:" << f.tri(points) << " d:" << d << abort(FatalError); } inters[interI].setHit(); inters[interI].setPoint(points[f[fp0]]); inters[interI].elementType() = triPointRef::POINT; inters[interI].setIndex(s.localFaces()[triI][fp0]); interI++; } else if ( (d[fp0] < 0 && d[fp1] > 0) || (d[fp0] > 0 && d[fp1] < 0) ) { if (interI >= 2) { FatalErrorIn("cutEdge(..)") << "problem : triangle has three intersections." << nl << "triangle:" << f.tri(points) << " d:" << d << abort(FatalError); } inters[interI].setHit(); inters[interI].setPoint ( (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]]) / (d[fp0] - d[fp1]) ); inters[interI].elementType() = triPointRef::EDGE; inters[interI].setIndex(fEdges[fp0]); interI++; } } if (interI == 0) { // Return miss } else if (interI == 1) { // Only one intersection. Should not happen! cut = inters[0]; } else if (interI == 2) { // Handle excludeEdgeI if ( inters[0].elementType() == triPointRef::EDGE && inters[0].index() == excludeEdgeI ) { cut = inters[1]; } else if ( inters[1].elementType() == triPointRef::EDGE && inters[1].index() == excludeEdgeI ) { cut = inters[0]; } else { // Two cuts. Find nearest. if ( magSqr(inters[0].rawPoint() - toPoint) < magSqr(inters[1].rawPoint() - toPoint) ) { cut = inters[0]; } else { cut = inters[1]; } } } } return cut; } // 'Snap' point on to endPoint. void Foam::triSurfaceTools::snapToEnd ( const triSurface& s, const surfaceLocation& end, surfaceLocation& current ) { if (end.elementType() == triPointRef::NONE) { if (current.elementType() == triPointRef::NONE) { // endpoint on triangle; current on triangle if (current.index() == end.index()) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } // No need to handle current on edge/point since tracking handles this. } else if (end.elementType() == triPointRef::EDGE) { if (current.elementType() == triPointRef::NONE) { // endpoint on edge; current on triangle const labelList& fEdges = s.faceEdges()[current.index()]; if (findIndex(fEdges, end.index()) != -1) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } else if (current.elementType() == triPointRef::EDGE) { // endpoint on edge; current on edge if (current.index() == end.index()) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } else { // endpoint on edge; current on point const edge& e = s.edges()[end.index()]; if (current.index() == e[0] || current.index() == e[1]) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } } else // end.elementType() == POINT { if (current.elementType() == triPointRef::NONE) { // endpoint on point; current on triangle const triSurface::FaceType& f = s.localFaces()[current.index()]; if (findIndex(f, end.index()) != -1) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } else if (current.elementType() == triPointRef::EDGE) { // endpoint on point; current on edge const edge& e = s.edges()[current.index()]; if (end.index() == e[0] || end.index() == e[1]) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } else { // endpoint on point; current on point if (current.index() == end.index()) { //if (debug) //{ // Pout<< "snapToEnd : snapping:" << current << " onto:" // << end << endl; //} current = end; current.setHit(); } } } } // Start: // - location // - element type (triangle/edge/point) and index // - triangle to exclude Foam::surfaceLocation Foam::triSurfaceTools::visitFaces ( const triSurface& s, const labelList& eFaces, const surfaceLocation& start, const label excludeEdgeI, const label excludePointI, const surfaceLocation& end, const plane& cutPlane ) { surfaceLocation nearest; scalar minDistSqr = Foam::sqr(GREAT); forAll(eFaces, i) { label triI = eFaces[i]; // Make sure we don't revisit previous face if (triI != start.triangle()) { if (end.elementType() == triPointRef::NONE && end.index() == triI) { // Endpoint is in this triangle. Jump there. nearest = end; nearest.setHit(); nearest.triangle() = triI; break; } else { // Which edge is cut. surfaceLocation cutInfo = cutEdge ( s, triI, excludeEdgeI, // excludeEdgeI excludePointI, // excludePointI start.rawPoint(), cutPlane, end.rawPoint() ); // If crossing an edge we expect next edge to be cut. if (excludeEdgeI != -1 && !cutInfo.hit()) { FatalErrorIn("triSurfaceTools::visitFaces(..)") << "Triangle:" << triI << " excludeEdge:" << excludeEdgeI << " point:" << start.rawPoint() << " plane:" << cutPlane << " . No intersection!" << abort(FatalError); } if (cutInfo.hit()) { scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint()); if (distSqr < minDistSqr) { minDistSqr = distSqr; nearest = cutInfo; nearest.triangle() = triI; nearest.setMiss(); } } } } } if (nearest.triangle() == -1) { // Did not move from edge. Give warning? Return something special? // For now responsability of caller to make sure that nothing has // moved. } return nearest; } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // Write pointField to file void Foam::triSurfaceTools::writeOBJ ( const fileName& fName, const pointField& pts ) { OFstream outFile(fName); forAll(pts, pointI) { const point& pt = pts[pointI]; outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl; } Pout<< "Written " << pts.size() << " vertices to file " << fName << endl; } // Write vertex subset to OBJ format file void Foam::triSurfaceTools::writeOBJ ( const triSurface& surf, const fileName& fName, const boolList& markedVerts ) { OFstream outFile(fName); label nVerts = 0; forAll(markedVerts, vertI) { if (markedVerts[vertI]) { const point& pt = surf.localPoints()[vertI]; outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl; nVerts++; } } Pout<< "Written " << nVerts << " vertices to file " << fName << endl; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Addressing helper functions: // Get all triangles using vertices of edge void Foam::triSurfaceTools::getVertexTriangles ( const triSurface& surf, const label edgeI, labelList& edgeTris ) { const edge& e = surf.edges()[edgeI]; const labelList& myFaces = surf.edgeFaces()[edgeI]; label face1I = myFaces[0]; label face2I = -1; if (myFaces.size() == 2) { face2I = myFaces[1]; } const labelList& startFaces = surf.pointFaces()[e.start()]; const labelList& endFaces = surf.pointFaces()[e.end()]; // Number of triangles is sum of pointfaces - common faces // (= faces using edge) edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size()); label nTris = 0; forAll(startFaces, startFaceI) { edgeTris[nTris++] = startFaces[startFaceI]; } forAll(endFaces, endFaceI) { label faceI = endFaces[endFaceI]; if ((faceI != face1I) && (faceI != face2I)) { edgeTris[nTris++] = faceI; } } } // Get all vertices connected to vertices of edge Foam::labelList Foam::triSurfaceTools::getVertexVertices ( const triSurface& surf, const edge& e ) { const edgeList& edges = surf.edges(); label v1 = e.start(); label v2 = e.end(); // Get all vertices connected to v1 or v2 through an edge labelHashSet vertexNeighbours; const labelList& v1Edges = surf.pointEdges()[v1]; forAll(v1Edges, v1EdgeI) { const edge& e = edges[v1Edges[v1EdgeI]]; vertexNeighbours.insert(e.otherVertex(v1)); } const labelList& v2Edges = surf.pointEdges()[v2]; forAll(v2Edges, v2EdgeI) { const edge& e = edges[v2Edges[v2EdgeI]]; label vertI = e.otherVertex(v2); vertexNeighbours.insert(vertI); } return vertexNeighbours.toc(); } //// Order vertices consistent with face //void Foam::triSurfaceTools::orderVertices //( // const labelledTri& f, // const label v1, // const label v2, // label& vA, // label& vB //) //{ // // Order v1, v2 in anticlockwise order. // bool reverse = false; // // if (f[0] == v1) // { // if (f[1] != v2) // { // reverse = true; // } // } // else if (f[1] == v1) // { // if (f[2] != v2) // { // reverse = true; // } // } // else // { // if (f[0] != v2) // { // reverse = true; // } // } // // if (reverse) // { // vA = v2; // vB = v1; // } // else // { // vA = v1; // vB = v2; // } //} // Get the other face using edgeI Foam::label Foam::triSurfaceTools::otherFace ( const triSurface& surf, const label faceI, const label edgeI ) { const labelList& myFaces = surf.edgeFaces()[edgeI]; if (myFaces.size() != 2) { return -1; } else { if (faceI == myFaces[0]) { return myFaces[1]; } else { return myFaces[0]; } } } // Get the two edges on faceI counterclockwise after edgeI void Foam::triSurfaceTools::otherEdges ( const triSurface& surf, const label faceI, const label edgeI, label& e1, label& e2 ) { const labelList& eFaces = surf.faceEdges()[faceI]; label i0 = findIndex(eFaces, edgeI); if (i0 == -1) { FatalErrorIn ( "otherEdges" "(const triSurface&, const label, const label," " label&, label&)" ) << "Edge " << surf.edges()[edgeI] << " not in face " << surf.localFaces()[faceI] << abort(FatalError); } label i1 = eFaces.fcIndex(i0); label i2 = eFaces.fcIndex(i1); e1 = eFaces[i1]; e2 = eFaces[i2]; } // Get the two vertices on faceI counterclockwise vertI void Foam::triSurfaceTools::otherVertices ( const triSurface& surf, const label faceI, const label vertI, label& vert1I, label& vert2I ) { const labelledTri& f = surf.localFaces()[faceI]; if (vertI == f[0]) { vert1I = f[1]; vert2I = f[2]; } else if (vertI == f[1]) { vert1I = f[2]; vert2I = f[0]; } else if (vertI == f[2]) { vert1I = f[0]; vert2I = f[1]; } else { FatalErrorIn ( "otherVertices" "(const triSurface&, const label, const label," " label&, label&)" ) << "Vertex " << vertI << " not in face " << f << abort(FatalError); } } // Get edge opposite vertex Foam::label Foam::triSurfaceTools::oppositeEdge ( const triSurface& surf, const label faceI, const label vertI ) { const labelList& myEdges = surf.faceEdges()[faceI]; forAll(myEdges, myEdgeI) { label edgeI = myEdges[myEdgeI]; const edge& e = surf.edges()[edgeI]; if ((e.start() != vertI) && (e.end() != vertI)) { return edgeI; } } FatalErrorIn ( "oppositeEdge" "(const triSurface&, const label, const label)" ) << "Cannot find vertex " << vertI << " in edges of face " << faceI << abort(FatalError); return -1; } // Get vertex opposite edge Foam::label Foam::triSurfaceTools::oppositeVertex ( const triSurface& surf, const label faceI, const label edgeI ) { const triSurface::FaceType& f = surf.localFaces()[faceI]; const edge& e = surf.edges()[edgeI]; forAll(f, fp) { label vertI = f[fp]; if (vertI != e.start() && vertI != e.end()) { return vertI; } } FatalErrorIn("triSurfaceTools::oppositeVertex") << "Cannot find vertex opposite edge " << edgeI << " vertices " << e << " in face " << faceI << " vertices " << f << abort(FatalError); return -1; } // Returns edge label connecting v1, v2 Foam::label Foam::triSurfaceTools::getEdge ( const triSurface& surf, const label v1, const label v2 ) { const labelList& v1Edges = surf.pointEdges()[v1]; forAll(v1Edges, v1EdgeI) { label edgeI = v1Edges[v1EdgeI]; const edge& e = surf.edges()[edgeI]; if ((e.start() == v2) || (e.end() == v2)) { return edgeI; } } return -1; } // Return index of triangle (or -1) using all three edges Foam::label Foam::triSurfaceTools::getTriangle ( const triSurface& surf, const label e0I, const label e1I, const label e2I ) { if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I)) { FatalErrorIn ( "getTriangle" "(const triSurface&, const label, const label," " const label)" ) << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I << " e2:" << e2I << abort(FatalError); } const labelList& eFaces = surf.edgeFaces()[e0I]; forAll(eFaces, eFaceI) { label faceI = eFaces[eFaceI]; const labelList& myEdges = surf.faceEdges()[faceI]; if ( (myEdges[0] == e1I) || (myEdges[1] == e1I) || (myEdges[2] == e1I) ) { if ( (myEdges[0] == e2I) || (myEdges[1] == e2I) || (myEdges[2] == e2I) ) { return faceI; } } } return -1; } // Collapse indicated edges. Return new tri surface. Foam::triSurface Foam::triSurfaceTools::collapseEdges ( const triSurface& surf, const labelList& collapsableEdges ) { pointField edgeMids(surf.nEdges()); forAll(edgeMids, edgeI) { const edge& e = surf.edges()[edgeI]; edgeMids[edgeI] = 0.5 * ( surf.localPoints()[e.start()] + surf.localPoints()[e.end()] ); } labelList faceStatus(surf.size(), ANYEDGE); //// Protect triangles which are on the border of different regions //forAll(edges, edgeI) //{ // const labelList& neighbours = edgeFaces[edgeI]; // // if ((neighbours.size() != 2) && (neighbours.size() != 1)) // { // FatalErrorIn("collapseEdges") // << abort(FatalError); // } // // if (neighbours.size() == 2) // { // if (surf[neighbours[0]].region() != surf[neighbours[1]].region()) // { // // Neighbours on different regions. For now dont allow // // any collapse. // //Pout<< "protecting face " << neighbours[0] // // << ' ' << neighbours[1] << endl; // faceStatus[neighbours[0]] = NOEDGE; // faceStatus[neighbours[1]] = NOEDGE; // } // } //} return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus); } // Collapse indicated edges. Return new tri surface. Foam::triSurface Foam::triSurfaceTools::collapseEdges ( const triSurface& surf, const labelList& collapseEdgeLabels, const pointField& edgeMids, labelList& faceStatus ) { const labelListList& edgeFaces = surf.edgeFaces(); const pointField& localPoints = surf.localPoints(); const edgeList& edges = surf.edges(); // Storage for new points pointField newPoints(localPoints); // Map for old to new points labelList pointMap(localPoints.size()); forAll(localPoints, pointI) { pointMap[pointI] = pointI; } // Do actual 'collapsing' of edges forAll(collapseEdgeLabels, collapseEdgeI) { const label edgeI = collapseEdgeLabels[collapseEdgeI]; if ((edgeI < 0) || (edgeI >= surf.nEdges())) { FatalErrorIn("collapseEdges") << "Edge label outside valid range." << endl << "edge label:" << edgeI << endl << "total number of edges:" << surf.nEdges() << endl << abort(FatalError); } const labelList& neighbours = edgeFaces[edgeI]; if (neighbours.size() == 2) { const label stat0 = faceStatus[neighbours[0]]; const label stat1 = faceStatus[neighbours[1]]; // Check faceStatus to make sure this one can be collapsed if ( ((stat0 == ANYEDGE) || (stat0 == edgeI)) && ((stat1 == ANYEDGE) || (stat1 == edgeI)) ) { const edge& e = edges[edgeI]; // Set up mapping to 'collapse' points of edge if ( (pointMap[e.start()] != e.start()) || (pointMap[e.end()] != e.end()) ) { FatalErrorIn("collapseEdges") << "points already mapped. Double collapse." << endl << "edgeI:" << edgeI << " start:" << e.start() << " end:" << e.end() << " pointMap[start]:" << pointMap[e.start()] << " pointMap[end]:" << pointMap[e.end()] << abort(FatalError); } const label minVert = min(e.start(), e.end()); pointMap[e.start()] = minVert; pointMap[e.end()] = minVert; // Move shared vertex to mid of edge newPoints[minVert] = edgeMids[edgeI]; // Protect neighbouring faces protectNeighbours(surf, e.start(), faceStatus); protectNeighbours(surf, e.end(), faceStatus); protectNeighbours ( surf, oppositeVertex(surf, neighbours[0], edgeI), faceStatus ); protectNeighbours ( surf, oppositeVertex(surf, neighbours[1], edgeI), faceStatus ); // Mark all collapsing faces labelList collapseFaces = getCollapsedFaces ( surf, edgeI ).toc(); forAll(collapseFaces, collapseI) { faceStatus[collapseFaces[collapseI]] = COLLAPSED; } } } } // Storage for new triangles List newTris(surf.size()); label newTriI = 0; const List& localFaces = surf.localFaces(); // Get only non-collapsed triangles and renumber vertex labels. forAll(localFaces, faceI) { const labelledTri& f = localFaces[faceI]; const label a = pointMap[f[0]]; const label b = pointMap[f[1]]; const label c = pointMap[f[2]]; if ( (a != b) && (a != c) && (b != c) && (faceStatus[faceI] != COLLAPSED) ) { // uncollapsed triangle newTris[newTriI++] = labelledTri(a, b, c, f.region()); } else { //Pout<< "Collapsed triangle " << faceI // << " vertices:" << f << endl; } } newTris.setSize(newTriI); // Pack faces triSurface tempSurf(newTris, surf.patches(), newPoints); return triSurface ( tempSurf.localFaces(), tempSurf.patches(), tempSurf.localPoints() ); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Foam::triSurface Foam::triSurfaceTools::redGreenRefine ( const triSurface& surf, const labelList& refineFaces ) { List refineStatus(surf.size(), NONE); // Mark & propagate refinement forAll(refineFaces, refineFaceI) { calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus); } // Do actual refinement return doRefine(surf, refineStatus); } Foam::triSurface Foam::triSurfaceTools::greenRefine ( const triSurface& surf, const labelList& refineEdges ) { // Storage for marking faces List refineStatus(surf.size(), NONE); // Storage for new faces DynamicList newFaces(0); pointField newPoints(surf.localPoints()); newPoints.setSize(surf.nPoints() + surf.nEdges()); label newPointI = surf.nPoints(); // Refine edges forAll(refineEdges, refineEdgeI) { label edgeI = refineEdges[refineEdgeI]; const labelList& myFaces = surf.edgeFaces()[edgeI]; bool neighbourIsRefined= false; forAll(myFaces, myFaceI) { if (refineStatus[myFaces[myFaceI]] != NONE) { neighbourIsRefined = true; } } // Only refine if none of the faces is refined if (!neighbourIsRefined) { // Refine edge const edge& e = surf.edges()[edgeI]; point mid = 0.5 * ( surf.localPoints()[e.start()] + surf.localPoints()[e.end()] ); newPoints[newPointI] = mid; // Refine faces using edge forAll(myFaces, myFaceI) { // Add faces to newFaces greenRefine ( surf, myFaces[myFaceI], edgeI, newPointI, newFaces ); // Mark as refined refineStatus[myFaces[myFaceI]] = GREEN; } newPointI++; } } // Add unrefined faces forAll(surf.localFaces(), faceI) { if (refineStatus[faceI] == NONE) { newFaces.append(surf.localFaces()[faceI]); } } newFaces.shrink(); newPoints.setSize(newPointI); return triSurface(newFaces, surf.patches(), newPoints, true); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Geometric helper functions: // Returns element in edgeIndices with minimum length Foam::label Foam::triSurfaceTools::minEdge ( const triSurface& surf, const labelList& edgeIndices ) { scalar minLength = GREAT; label minIndex = -1; forAll(edgeIndices, i) { const edge& e = surf.edges()[edgeIndices[i]]; scalar length = mag ( surf.localPoints()[e.end()] - surf.localPoints()[e.start()] ); if (length < minLength) { minLength = length; minIndex = i; } } return edgeIndices[minIndex]; } // Returns element in edgeIndices with maximum length Foam::label Foam::triSurfaceTools::maxEdge ( const triSurface& surf, const labelList& edgeIndices ) { scalar maxLength = -GREAT; label maxIndex = -1; forAll(edgeIndices, i) { const edge& e = surf.edges()[edgeIndices[i]]; scalar length = mag ( surf.localPoints()[e.end()] - surf.localPoints()[e.start()] ); if (length > maxLength) { maxLength = length; maxIndex = i; } } return edgeIndices[maxIndex]; } // Merge points and reconstruct surface Foam::triSurface Foam::triSurfaceTools::mergePoints ( const triSurface& surf, const scalar mergeTol ) { pointField newPoints(surf.nPoints()); labelList pointMap(surf.nPoints()); bool hasMerged = Foam::mergePoints ( surf.localPoints(), mergeTol, false, pointMap, newPoints ); if (hasMerged) { // Pack the triangles // Storage for new triangles List newTriangles(surf.size()); label newTriangleI = 0; forAll(surf, faceI) { const labelledTri& f = surf.localFaces()[faceI]; label newA = pointMap[f[0]]; label newB = pointMap[f[1]]; label newC = pointMap[f[2]]; if ((newA != newB) && (newA != newC) && (newB != newC)) { newTriangles[newTriangleI++] = labelledTri(newA, newB, newC, f.region()); } } newTriangles.setSize(newTriangleI); return triSurface ( newTriangles, surf.patches(), newPoints, true //reuse storage ); } else { return surf; } } // Calculate normal on triangle Foam::vector Foam::triSurfaceTools::surfaceNormal ( const triSurface& surf, const label nearestFaceI, const point& nearestPt ) { const triSurface::FaceType& f = surf[nearestFaceI]; const pointField& points = surf.points(); label nearType, nearLabel; f.nearestPointClassify(nearestPt, points, nearType, nearLabel); if (nearType == triPointRef::NONE) { // Nearest to face return surf.faceNormals()[nearestFaceI]; } else if (nearType == triPointRef::EDGE) { // Nearest to edge. Assume order of faceEdges same as face vertices. label edgeI = surf.faceEdges()[nearestFaceI][nearLabel]; // Calculate edge normal by averaging face normals const labelList& eFaces = surf.edgeFaces()[edgeI]; vector edgeNormal(vector::zero); forAll(eFaces, i) { edgeNormal += surf.faceNormals()[eFaces[i]]; } return edgeNormal/(mag(edgeNormal) + VSMALL); } else { // Nearest to point const triSurface::FaceType& localF = surf.localFaces()[nearestFaceI]; return surf.pointNormals()[localF[nearLabel]]; } } Foam::triSurfaceTools::sideType Foam::triSurfaceTools::edgeSide ( const triSurface& surf, const point& sample, const point& nearestPoint, const label edgeI ) { const labelList& eFaces = surf.edgeFaces()[edgeI]; if (eFaces.size() != 2) { // Surface not closed. return UNKNOWN; } else { const vectorField& faceNormals = surf.faceNormals(); // Compare to bisector. This is actually correct since edge is // nearest so there is a knife-edge. vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]); if (((sample - nearestPoint) & n) > 0) { return OUTSIDE; } else { return INSIDE; } } } // Calculate normal on triangle Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide ( const triSurface& surf, const point& sample, const label nearestFaceI ) { const triSurface::FaceType& f = surf[nearestFaceI]; const pointField& points = surf.points(); // Find where point is on face label nearType, nearLabel; pointHit pHit = f.nearestPointClassify(sample, points, nearType, nearLabel); const point& nearestPoint(pHit.rawPoint()); if (nearType == triPointRef::NONE) { vector sampleNearestVec = (sample - nearestPoint); // Nearest to face interior. Use faceNormal to determine side scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI]; // // If the sample is essentially on the face, do not check for // // it being perpendicular. // scalar magSampleNearestVec = mag(sampleNearestVec); // if (magSampleNearestVec > SMALL) // { // c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]); // if (mag(c) < 0.99) // { // FatalErrorIn("triSurfaceTools::surfaceSide") // << "nearestPoint identified as being on triangle face " // << "but vector from nearestPoint to sample is not " // << "perpendicular to the normal." << nl // << "sample: " << sample << nl // << "nearestPoint: " << nearestPoint << nl // << "sample - nearestPoint: " // << sample - nearestPoint << nl // << "normal: " << surf.faceNormals()[nearestFaceI] << nl // << "mag(sample - nearestPoint): " // << mag(sample - nearestPoint) << nl // << "normalised dot product: " << c << nl // << "triangle vertices: " << nl // << " " << points[f[0]] << nl // << " " << points[f[1]] << nl // << " " << points[f[2]] << nl // << abort(FatalError); // } // } if (c > 0) { return OUTSIDE; } else { return INSIDE; } } else if (nearType == triPointRef::EDGE) { // Nearest to edge nearLabel. Note that this can only be a knife-edge // situation since otherwise the nearest point could never be the edge. // Get the edge. Assume order of faceEdges same as face vertices. label edgeI = surf.faceEdges()[nearestFaceI][nearLabel]; // if (debug) // { // // Check order of faceEdges same as face vertices. // const edge& e = surf.edges()[edgeI]; // const labelList& meshPoints = surf.meshPoints(); // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]); // if // ( // meshEdge // != edge(f[nearLabel], f[f.fcIndex(nearLabel)]) // ) // { // FatalErrorIn("triSurfaceTools::surfaceSide") // << "Edge:" << edgeI << " local vertices:" << e // << " mesh vertices:" << meshEdge // << " not at position " << nearLabel // << " in face " << f // << abort(FatalError); // } // } return edgeSide(surf, sample, nearestPoint, edgeI); } else { // Nearest to point. Could use pointNormal here but is not correct. // Instead determine which edge using point is nearest and use test // above (nearType == triPointRef::EDGE). const triSurface::FaceType& localF = surf.localFaces()[nearestFaceI]; label nearPointI = localF[nearLabel]; const edgeList& edges = surf.edges(); const pointField& localPoints = surf.localPoints(); const point& base = localPoints[nearPointI]; const labelList& pEdges = surf.pointEdges()[nearPointI]; scalar minDistSqr = Foam::sqr(GREAT); label minEdgeI = -1; forAll(pEdges, i) { label edgeI = pEdges[i]; const edge& e = edges[edgeI]; label otherPointI = e.otherVertex(nearPointI); // Get edge normal. vector eVec(localPoints[otherPointI] - base); scalar magEVec = mag(eVec); if (magEVec > VSMALL) { eVec /= magEVec; // Get point along vector and determine closest. const point perturbPoint = base + eVec; scalar distSqr = Foam::magSqr(sample - perturbPoint); if (distSqr < minDistSqr) { minDistSqr = distSqr; minEdgeI = edgeI; } } } if (minEdgeI == -1) { FatalErrorIn("treeDataTriSurface::getSide") << "Problem: did not find edge closer than " << minDistSqr << abort(FatalError); } return edgeSide(surf, sample, nearestPoint, minEdgeI); } } // triangulation of boundaryMesh Foam::triSurface Foam::triSurfaceTools::triangulate ( const polyBoundaryMesh& bMesh, const labelHashSet& includePatches, const bool verbose ) { const polyMesh& mesh = bMesh.mesh(); // Storage for surfaceMesh. Size estimate. DynamicList triangles ( mesh.nFaces() - mesh.nInternalFaces() ); label newPatchI = 0; forAllConstIter(labelHashSet, includePatches, iter) { const label patchI = iter.key(); const polyPatch& patch = bMesh[patchI]; const pointField& points = patch.points(); label nTriTotal = 0; forAll(patch, patchFaceI) { const face& f = patch[patchFaceI]; faceList triFaces(f.nTriangles(points)); label nTri = 0; f.triangles(points, nTri, triFaces); forAll(triFaces, triFaceI) { const face& f = triFaces[triFaceI]; triangles.append(labelledTri(f[0], f[1], f[2], newPatchI)); nTriTotal++; } } if (verbose) { Pout<< patch.name() << " : generated " << nTriTotal << " triangles from " << patch.size() << " faces with" << " new patchid " << newPatchI << endl; } newPatchI++; } triangles.shrink(); // Create globally numbered tri surface triSurface rawSurface(triangles, mesh.points()); // Create locally numbered tri surface triSurface surface ( rawSurface.localFaces(), rawSurface.localPoints() ); // Add patch names to surface surface.patches().setSize(newPatchI); newPatchI = 0; forAllConstIter(labelHashSet, includePatches, iter) { const label patchI = iter.key(); const polyPatch& patch = bMesh[patchI]; surface.patches()[newPatchI].name() = patch.name(); surface.patches()[newPatchI].geometricType() = patch.type(); newPatchI++; } return surface; } // triangulation of boundaryMesh Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre ( const polyBoundaryMesh& bMesh, const labelHashSet& includePatches, const bool verbose ) { const polyMesh& mesh = bMesh.mesh(); // Storage for new points = meshpoints + face centres. const pointField& points = mesh.points(); const pointField& faceCentres = mesh.faceCentres(); pointField newPoints(points.size() + faceCentres.size()); label newPointI = 0; forAll(points, pointI) { newPoints[newPointI++] = points[pointI]; } forAll(faceCentres, faceI) { newPoints[newPointI++] = faceCentres[faceI]; } // Count number of faces. DynamicList triangles ( mesh.nFaces() - mesh.nInternalFaces() ); label newPatchI = 0; forAllConstIter(labelHashSet, includePatches, iter) { const label patchI = iter.key(); const polyPatch& patch = bMesh[patchI]; label nTriTotal = 0; forAll(patch, patchFaceI) { // Face in global coords. const face& f = patch[patchFaceI]; // Index in newPointI of face centre. label fc = points.size() + patchFaceI + patch.start(); forAll(f, fp) { label fp1 = f.fcIndex(fp); triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchI)); nTriTotal++; } } if (verbose) { Pout<< patch.name() << " : generated " << nTriTotal << " triangles from " << patch.size() << " faces with" << " new patchid " << newPatchI << endl; } newPatchI++; } triangles.shrink(); // Create globally numbered tri surface triSurface rawSurface(triangles, newPoints); // Create locally numbered tri surface triSurface surface ( rawSurface.localFaces(), rawSurface.localPoints() ); // Add patch names to surface surface.patches().setSize(newPatchI); newPatchI = 0; forAllConstIter(labelHashSet, includePatches, iter) { const label patchI = iter.key(); const polyPatch& patch = bMesh[patchI]; surface.patches()[newPatchI].name() = patch.name(); surface.patches()[newPatchI].geometricType() = patch.type(); newPatchI++; } return surface; } Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List& pts) { // Vertices in geompack notation. Note that could probably just use // pts.begin() if double precision. List geompackVertices(2*pts.size()); label doubleI = 0; forAll(pts, i) { geompackVertices[doubleI++] = pts[i][0]; geompackVertices[doubleI++] = pts[i][1]; } // Storage for triangles int m2 = 3; List triangle_node(m2*3*pts.size()); List triangle_neighbor(m2*3*pts.size()); // Triangulate int nTris = 0; int err = dtris2 ( pts.size(), geompackVertices.begin(), &nTris, triangle_node.begin(), triangle_neighbor.begin() ); if (err != 0) { FatalErrorIn("triSurfaceTools::delaunay2D(const List&)") << "Failed dtris2 with vertices:" << pts.size() << abort(FatalError); } // Trim triangle_node.setSize(3*nTris); triangle_neighbor.setSize(3*nTris); // Convert to triSurface. List faces(nTris); forAll(faces, i) { faces[i] = labelledTri ( triangle_node[3*i]-1, triangle_node[3*i+1]-1, triangle_node[3*i+2]-1, 0 ); } pointField points(pts.size()); forAll(pts, i) { points[i][0] = pts[i][0]; points[i][1] = pts[i][1]; points[i][2] = 0.0; } return triSurface(faces, points); } void Foam::triSurfaceTools::calcInterpolationWeights ( const triPointRef& tri, const point& p, FixedList& weights ) { // calculate triangle edge vectors and triangle face normal // the 'i':th edge is opposite node i FixedList edge; edge[0] = tri.c()-tri.b(); edge[1] = tri.a()-tri.c(); edge[2] = tri.b()-tri.a(); vector triangleFaceNormal = edge[1] ^ edge[2]; // calculate edge normal (pointing inwards) FixedList normal; for (label i=0; i<3; i++) { normal[i] = triangleFaceNormal ^ edge[i]; normal[i] /= mag(normal[i]) + VSMALL; } weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]); weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]); weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]); } // Calculate weighting factors from samplePts to triangle it is in. // Uses linear search. void Foam::triSurfaceTools::calcInterpolationWeights ( const triSurface& s, const pointField& samplePts, List >& allVerts, List >& allWeights ) { allVerts.setSize(samplePts.size()); allWeights.setSize(samplePts.size()); const pointField& points = s.points(); forAll(samplePts, i) { const point& samplePt = samplePts[i]; FixedList& verts = allVerts[i]; FixedList& weights = allWeights[i]; scalar minDistance = GREAT; forAll(s, faceI) { const labelledTri& f = s[faceI]; triPointRef tri(f.tri(points)); label nearType, nearLabel; pointHit nearest = tri.nearestPointClassify ( samplePt, nearType, nearLabel ); if (nearest.hit()) { // samplePt inside triangle verts[0] = f[0]; verts[1] = f[1]; verts[2] = f[2]; calcInterpolationWeights(tri, nearest.rawPoint(), weights); //Pout<< "calcScalingFactors : samplePt:" << samplePt // << " inside triangle:" << faceI // << " verts:" << verts // << " weights:" << weights // << endl; break; } else if (nearest.distance() < minDistance) { minDistance = nearest.distance(); // Outside triangle. Store nearest. if (nearType == triPointRef::POINT) { verts[0] = f[nearLabel]; weights[0] = 1; verts[1] = -1; weights[1] = -GREAT; verts[2] = -1; weights[2] = -GREAT; //Pout<< "calcScalingFactors : samplePt:" << samplePt // << " distance:" << nearest.distance() // << " from point:" << points[f[nearLabel]] // << endl; } else if (nearType == triPointRef::EDGE) { verts[0] = f[nearLabel]; verts[1] = f[f.fcIndex(nearLabel)]; verts[2] = -1; const point& p0 = points[verts[0]]; const point& p1 = points[verts[1]]; scalar s = min ( 1, max ( 0, mag(nearest.rawPoint() - p0)/mag(p1 - p0) ) ); // Interpolate weights[0] = 1 - s; weights[1] = s; weights[2] = -GREAT; //Pout<< "calcScalingFactors : samplePt:" << samplePt // << " distance:" << nearest.distance() // << " from edge:" << p0 << p1 << " s:" << s // << endl; } else { // triangle. Can only happen because of truncation errors. verts[0] = f[0]; verts[1] = f[1]; verts[2] = f[2]; calcInterpolationWeights(tri, nearest.rawPoint(), weights); break; } } } } } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Tracking: // Test point on surface to see if is on face,edge or point. Foam::surfaceLocation Foam::triSurfaceTools::classify ( const triSurface& s, const label triI, const point& trianglePoint ) { surfaceLocation nearest; // Nearest point could be on point or edge. Retest. label index, elemType; //bool inside = triPointRef(s[triI].tri(s.points())).classify ( trianglePoint, elemType, index ); nearest.setPoint(trianglePoint); if (elemType == triPointRef::NONE) { nearest.setHit(); nearest.setIndex(triI); nearest.elementType() = triPointRef::NONE; } else if (elemType == triPointRef::EDGE) { nearest.setMiss(); nearest.setIndex(s.faceEdges()[triI][index]); nearest.elementType() = triPointRef::EDGE; } else //if (elemType == triPointRef::POINT) { nearest.setMiss(); nearest.setIndex(s.localFaces()[triI][index]); nearest.elementType() = triPointRef::POINT; } return nearest; } Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge ( const triSurface& s, const surfaceLocation& start, const surfaceLocation& end, const plane& cutPlane ) { // Start off from starting point surfaceLocation nearest = start; nearest.setMiss(); // See if in same triangle as endpoint. If so snap. snapToEnd(s, end, nearest); if (!nearest.hit()) { // Not yet at end point if (start.elementType() == triPointRef::NONE) { // Start point is inside triangle. Trivial cases already handled // above. // end point is on edge or point so cross currrent triangle to // see which edge is cut. nearest = cutEdge ( s, start.index(), // triangle -1, // excludeEdge -1, // excludePoint start.rawPoint(), cutPlane, end.rawPoint() ); nearest.elementType() = triPointRef::EDGE; nearest.triangle() = start.index(); nearest.setMiss(); } else if (start.elementType() == triPointRef::EDGE) { // Pick connected triangle that is most in direction. const labelList& eFaces = s.edgeFaces()[start.index()]; nearest = visitFaces ( s, eFaces, start, start.index(), // excludeEdgeI -1, // excludePointI end, cutPlane ); } else // start.elementType() == triPointRef::POINT { const labelList& pFaces = s.pointFaces()[start.index()]; nearest = visitFaces ( s, pFaces, start, -1, // excludeEdgeI start.index(), // excludePointI end, cutPlane ); } snapToEnd(s, end, nearest); } return nearest; } void Foam::triSurfaceTools::track ( const triSurface& s, const surfaceLocation& endInfo, const plane& cutPlane, surfaceLocation& hitInfo ) { //OFstream str("track.obj"); //label vertI = 0; //meshTools::writeOBJ(str, hitInfo.rawPoint()); //vertI++; // Track across surface. while (true) { //Pout<< "Tracking from:" << nl // << " " << hitInfo.info() // << endl; hitInfo = trackToEdge ( s, hitInfo, endInfo, cutPlane ); //meshTools::writeOBJ(str, hitInfo.rawPoint()); //vertI++; //str<< "l " << vertI-1 << ' ' << vertI << nl; //Pout<< "Tracked to:" << nl // << " " << hitInfo.info() << endl; if (hitInfo.hit() || hitInfo.triangle() == -1) { break; } } } // ************************************************************************* //