FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

triSurfaceTools.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025 
00026 \*---------------------------------------------------------------------------*/
00027 
00028 #include <meshTools/triSurfaceTools.H>
00029 
00030 #include <triSurface/triSurface.H>
00031 #include <OpenFOAM/OFstream.H>
00032 #include <OpenFOAM/mergePoints.H>
00033 #include <OpenFOAM/polyMesh.H>
00034 #include <OpenFOAM/plane.H>
00035 #include <meshTools/geompack.H>
00036 
00037 
00038 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00039 
00040 const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
00041 const Foam::label Foam::triSurfaceTools::NOEDGE = -2;
00042 const Foam::label Foam::triSurfaceTools::COLLAPSED = -3;
00043 
00044 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00045 
00046 /*
00047     Refine by splitting all three edges of triangle ('red' refinement).
00048     Neighbouring triangles (which are not marked for refinement get split
00049     in half ('green') refinement. (R. Verfuerth, "A review of a posteriori
00050     error estimation and adaptive mesh refinement techniques",
00051     Wiley-Teubner, 1996)
00052 */
00053 
00054 // FaceI gets refined ('red'). Propagate edge refinement.
00055 void Foam::triSurfaceTools::calcRefineStatus
00056 (
00057     const triSurface& surf,
00058     const label faceI,
00059     List<refineType>& refine
00060 )
00061 {
00062     if (refine[faceI] == RED)
00063     {
00064         // Already marked for refinement. Do nothing.
00065     }
00066     else
00067     {
00068         // Not marked or marked for 'green' refinement. Refine.
00069         refine[faceI] = RED;
00070 
00071         const labelList& myNeighbours = surf.faceFaces()[faceI];
00072 
00073         forAll(myNeighbours, myNeighbourI)
00074         {
00075             label neighbourFaceI = myNeighbours[myNeighbourI];
00076 
00077             if (refine[neighbourFaceI] == GREEN)
00078             {
00079                 // Change to red refinement and propagate
00080                 calcRefineStatus(surf, neighbourFaceI, refine);
00081             }
00082             else if (refine[neighbourFaceI] == NONE)
00083             {
00084                 refine[neighbourFaceI] = GREEN;
00085             }
00086         }
00087     }
00088 }
00089 
00090 
00091 // Split faceI along edgeI at position newPointI
00092 void Foam::triSurfaceTools::greenRefine
00093 (
00094     const triSurface& surf,
00095     const label faceI,
00096     const label edgeI,
00097     const label newPointI,
00098     DynamicList<labelledTri>& newFaces
00099 )
00100 {
00101     const labelledTri& f = surf.localFaces()[faceI];
00102     const edge& e = surf.edges()[edgeI];
00103 
00104     // Find index of edge in face.
00105 
00106     label fp0 = findIndex(f, e[0]);
00107     label fp1 = f.fcIndex(fp0);
00108     label fp2 = f.fcIndex(fp1);
00109 
00110     if (f[fp1] == e[1])
00111     {
00112         // Edge oriented like face
00113         newFaces.append
00114         (
00115             labelledTri
00116             (
00117                 f[fp0],
00118                 newPointI,
00119                 f[fp2],
00120                 f.region()
00121             )
00122         );
00123         newFaces.append
00124         (
00125             labelledTri
00126             (
00127                 newPointI,
00128                 f[fp1],
00129                 f[fp2],
00130                 f.region()
00131             )
00132         );
00133     }
00134     else
00135     {
00136         newFaces.append
00137         (
00138             labelledTri
00139             (
00140                 f[fp2],
00141                 newPointI,
00142                 f[fp1],
00143                 f.region()
00144             )
00145         );
00146         newFaces.append
00147         (
00148             labelledTri
00149             (
00150                 newPointI,
00151                 f[fp0],
00152                 f[fp1],
00153                 f.region()
00154             )
00155         );
00156     }
00157 }
00158 
00159 
00160 // Refine all triangles marked for refinement. 
00161 Foam::triSurface Foam::triSurfaceTools::doRefine
00162 (
00163     const triSurface& surf,
00164     const List<refineType>& refineStatus
00165 )
00166 {
00167     // Storage for new points. (start after old points)
00168     DynamicList<point> newPoints(surf.nPoints());
00169     forAll(surf.localPoints(), pointI)
00170     {
00171         newPoints.append(surf.localPoints()[pointI]);
00172     }
00173     label newVertI = surf.nPoints();
00174 
00175     // Storage for new faces
00176     DynamicList<labelledTri> newFaces(surf.size());
00177 
00178 
00179     // Point index for midpoint on edge
00180     labelList edgeMid(surf.nEdges(), -1);
00181 
00182     forAll(refineStatus, faceI)
00183     {
00184         if (refineStatus[faceI] == RED)
00185         {
00186             // Create new vertices on all edges to be refined.
00187             const labelList& fEdges = surf.faceEdges()[faceI];
00188 
00189             forAll(fEdges, i)
00190             {
00191                 label edgeI = fEdges[i];
00192 
00193                 if (edgeMid[edgeI] == -1)
00194                 {
00195                     const edge& e = surf.edges()[edgeI];
00196 
00197                     // Create new point on mid of edge
00198                     newPoints.append
00199                     (
00200                         0.5
00201                       * (
00202                             surf.localPoints()[e.start()]
00203                           + surf.localPoints()[e.end()]
00204                         )
00205                     );
00206                     edgeMid[edgeI] = newVertI++;
00207                 }
00208             }
00209 
00210             // Now we have new mid edge vertices for all edges on face
00211             // so create triangles for RED rerfinement.
00212 
00213             const edgeList& edges = surf.edges();
00214 
00215             // Corner triangles
00216             newFaces.append
00217             (
00218                 labelledTri
00219                 (
00220                     edgeMid[fEdges[0]],
00221                     edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
00222                     edgeMid[fEdges[1]],
00223                     surf[faceI].region()
00224                 )
00225             );
00226 
00227             newFaces.append
00228             (
00229                 labelledTri
00230                 (
00231                     edgeMid[fEdges[1]],
00232                     edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
00233                     edgeMid[fEdges[2]],
00234                     surf[faceI].region()
00235                 )
00236             );
00237 
00238             newFaces.append
00239             (
00240                 labelledTri
00241                 (
00242                     edgeMid[fEdges[2]],
00243                     edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
00244                     edgeMid[fEdges[0]],
00245                     surf[faceI].region()
00246                 )
00247             );
00248 
00249             // Inner triangle
00250             newFaces.append
00251             (
00252                 labelledTri
00253                 (
00254                     edgeMid[fEdges[0]],
00255                     edgeMid[fEdges[1]],
00256                     edgeMid[fEdges[2]],
00257                     surf[faceI].region()
00258                 )
00259             );
00260 
00261 
00262             // Create triangles for GREEN refinement.
00263             forAll(fEdges, i)
00264             {
00265                 const label edgeI = fEdges[i];
00266 
00267                 label otherFaceI = otherFace(surf, faceI, edgeI);
00268 
00269                 if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN))
00270                 {
00271                     greenRefine
00272                     (
00273                         surf,
00274                         otherFaceI,
00275                         edgeI,
00276                         edgeMid[edgeI],
00277                         newFaces
00278                     );
00279                 }
00280             }
00281         }
00282     }
00283 
00284     // Copy unmarked triangles since keep original vertex numbering.
00285     forAll(refineStatus, faceI)
00286     {
00287         if (refineStatus[faceI] == NONE)
00288         {
00289             newFaces.append(surf.localFaces()[faceI]);
00290         }
00291     }
00292 
00293     newFaces.shrink();
00294     newPoints.shrink();
00295 
00296 
00297     // Transfer DynamicLists to straight ones.
00298     pointField allPoints;
00299     allPoints.transfer(newPoints);
00300     newPoints.clear();
00301 
00302     return triSurface(newFaces, surf.patches(), allPoints, true);
00303 }
00304 
00305 
00306 // Angle between two neighbouring triangles,
00307 // angle between shared-edge end points and left and right face end points
00308 Foam::scalar Foam::triSurfaceTools::faceCosAngle
00309 (
00310     const point& pStart,
00311     const point& pEnd,
00312     const point& pLeft,
00313     const point& pRight
00314 )
00315 {
00316     const vector common(pEnd - pStart);
00317     const vector base0(pLeft - pStart);
00318     const vector base1(pRight - pStart);
00319 
00320     vector n0(common ^ base0);
00321     n0 /= Foam::mag(n0);
00322 
00323     vector n1(base1 ^ common);
00324     n1 /= Foam::mag(n1);
00325 
00326     return n0 & n1;
00327 }
00328 
00329 
00330 // Protect edges around vertex from collapsing.
00331 // Moving vertI to a different position
00332 // - affects obviously the cost of the faces using it
00333 // - but also their neighbours since the edge between the faces changes
00334 void Foam::triSurfaceTools::protectNeighbours
00335 (
00336     const triSurface& surf,
00337     const label vertI,
00338     labelList& faceStatus
00339 )
00340 {
00341 //    const labelList& myFaces = surf.pointFaces()[vertI];
00342 //    forAll(myFaces, i)
00343 //    {
00344 //        label faceI = myFaces[i];
00345 //
00346 //        if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
00347 //        {
00348 //            faceStatus[faceI] = NOEDGE;
00349 //        }
00350 //    }
00351 
00352     const labelList& myEdges = surf.pointEdges()[vertI];
00353     forAll(myEdges, i)
00354     {
00355         const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
00356 
00357         forAll(myFaces, myFaceI)
00358         {
00359             label faceI = myFaces[myFaceI];
00360  
00361             if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
00362             {
00363                 faceStatus[faceI] = NOEDGE;
00364             }
00365        }
00366     }
00367 }
00368 
00369 
00370 //
00371 // Edge collapse helper functions
00372 //
00373 
00374 // Get all faces that will get collapsed if edgeI collapses.
00375 Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
00376 (
00377     const triSurface& surf,
00378     label edgeI
00379 )
00380 {
00381     const edge& e = surf.edges()[edgeI];
00382     label v1 = e.start();
00383     label v2 = e.end();
00384 
00385     // Faces using edge will certainly get collapsed.
00386     const labelList& myFaces = surf.edgeFaces()[edgeI];
00387 
00388     labelHashSet facesToBeCollapsed(2*myFaces.size());
00389 
00390     forAll(myFaces, myFaceI)
00391     {
00392         facesToBeCollapsed.insert(myFaces[myFaceI]);
00393     }
00394     
00395     // From faces using v1 check if they share an edge with faces
00396     // using v2.
00397     //  - share edge: are part of 'splay' tree and will collapse if edge    
00398     //    collapses
00399     const labelList& v1Faces = surf.pointFaces()[v1];
00400 
00401     forAll(v1Faces, v1FaceI)
00402     {
00403         label face1I = v1Faces[v1FaceI];
00404 
00405         label otherEdgeI = oppositeEdge(surf, face1I, v1);
00406 
00407         // Step across edge to other face
00408         label face2I = otherFace(surf, face1I, otherEdgeI);
00409 
00410         if (face2I != -1)
00411         {
00412             // found face on other side of edge. Now check if uses v2.
00413             if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
00414             {
00415                 // triangles face1I and face2I form splay tree and will
00416                 // collapse.
00417                 facesToBeCollapsed.insert(face1I);
00418                 facesToBeCollapsed.insert(face2I);
00419             }
00420         }
00421     }
00422 
00423     return facesToBeCollapsed;
00424 }
00425 
00426 
00427 // Return value of faceUsed for faces using vertI
00428 Foam::label Foam::triSurfaceTools::vertexUsesFace
00429 (
00430     const triSurface& surf,
00431     const labelHashSet& faceUsed,
00432     const label vertI
00433 )
00434 {
00435     const labelList& myFaces = surf.pointFaces()[vertI];
00436 
00437     forAll(myFaces, myFaceI)
00438     {
00439         label face1I = myFaces[myFaceI];
00440 
00441         if (faceUsed.found(face1I))
00442         {
00443             return face1I;
00444         }
00445     }
00446     return -1;
00447 }
00448 
00449 
00450 // Calculate new edge-face addressing resulting from edge collapse
00451 void Foam::triSurfaceTools::getMergedEdges
00452 (
00453     const triSurface& surf,
00454     const label edgeI,
00455     const labelHashSet& collapsedFaces,
00456     HashTable<label, label, Hash<label> >& edgeToEdge,
00457     HashTable<label, label, Hash<label> >& edgeToFace
00458 )
00459 {
00460     const edge& e = surf.edges()[edgeI];
00461     label v1 = e.start();
00462     label v2 = e.end();
00463 
00464     const labelList& v1Faces = surf.pointFaces()[v1];
00465     const labelList& v2Faces = surf.pointFaces()[v2];
00466 
00467     // Mark all (non collapsed) faces using v2
00468     labelHashSet v2FacesHash(v2Faces.size());
00469 
00470     forAll(v2Faces, v2FaceI)
00471     {
00472         if (!collapsedFaces.found(v2Faces[v2FaceI]))
00473         {
00474             v2FacesHash.insert(v2Faces[v2FaceI]);
00475         }
00476     }
00477 
00478 
00479     forAll(v1Faces, v1FaceI)
00480     {
00481         label face1I = v1Faces[v1FaceI];
00482 
00483         if (collapsedFaces.found(face1I))
00484         {
00485             continue;
00486         }
00487 
00488         //
00489         // Check if face1I uses a vertex connected to a face using v2
00490         //
00491 
00492         label vert1I = -1;
00493         label vert2I = -1;
00494         otherVertices
00495         (
00496             surf,
00497             face1I,
00498             v1,
00499             vert1I,
00500             vert2I
00501         );
00502         //Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:"
00503         //    << vert1I << ' ' << vert2I << endl;
00504 
00505         // Check vert1, vert2 for usage by v2Face.
00506 
00507         label commonVert = vert1I;
00508         label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
00509         if (face2I == -1)
00510         {
00511             commonVert = vert2I;
00512             face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
00513         }
00514 
00515         if (face2I != -1)
00516         {
00517             // Found one: commonVert is used by both face1I and face2I
00518             label edge1I = getEdge(surf, v1, commonVert);
00519             label edge2I = getEdge(surf, v2, commonVert);
00520 
00521             edgeToEdge.insert(edge1I, edge2I);
00522             edgeToEdge.insert(edge2I, edge1I);
00523 
00524             edgeToFace.insert(edge1I, face2I);
00525             edgeToFace.insert(edge2I, face1I);
00526         }
00527     }
00528 }
00529 
00530 
00531 // Calculates (cos of) angle across edgeI of faceI,
00532 // taking into account updated addressing (resulting from edge collapse)
00533 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
00534 (
00535     const triSurface& surf,
00536     const label v1,
00537     const point& pt,
00538     const labelHashSet& collapsedFaces,
00539     const HashTable<label, label, Hash<label> >& edgeToEdge,
00540     const HashTable<label, label, Hash<label> >& edgeToFace,
00541     const label faceI,
00542     const label edgeI
00543 )
00544 {
00545     const pointField& localPoints = surf.localPoints();
00546 
00547     label A = surf.edges()[edgeI].start();
00548     label B = surf.edges()[edgeI].end();
00549     label C = oppositeVertex(surf, faceI, edgeI);
00550 
00551     label D = -1;
00552 
00553     label face2I = -1;
00554 
00555     if (edgeToEdge.found(edgeI))
00556     {
00557         // Use merged addressing
00558         label edge2I = edgeToEdge[edgeI];
00559         face2I = edgeToFace[edgeI];
00560 
00561         D = oppositeVertex(surf, face2I, edge2I);
00562     }
00563     else
00564     {
00565         // Use normal edge-face addressing
00566         face2I = otherFace(surf, faceI, edgeI);
00567 
00568         if ((face2I != -1) && !collapsedFaces.found(face2I))
00569         {
00570             D = oppositeVertex(surf, face2I, edgeI);
00571         }
00572     }
00573 
00574     scalar cosAngle = 1;
00575 
00576     if (D != -1)
00577     {
00578         if (A == v1)
00579         {
00580             cosAngle = faceCosAngle
00581             (
00582                 pt,
00583                 localPoints[B],
00584                 localPoints[C],
00585                 localPoints[D]
00586             );
00587         }
00588         else if (B == v1)
00589         {
00590             cosAngle = faceCosAngle
00591             (
00592                 localPoints[A],
00593                 pt,
00594                 localPoints[C],
00595                 localPoints[D]
00596             );
00597         }
00598         else if (C == v1)
00599         {
00600             cosAngle = faceCosAngle
00601             (
00602                 localPoints[A],
00603                 localPoints[B],
00604                 pt,
00605                 localPoints[D]
00606             );
00607         }
00608         else if (D == v1)
00609         {
00610             cosAngle = faceCosAngle
00611             (
00612                 localPoints[A],
00613                 localPoints[B],
00614                 localPoints[C],
00615                 pt
00616             );
00617         }
00618         else
00619         {
00620             FatalErrorIn("edgeCosAngle")
00621                 << "face " << faceI << " does not use vertex "
00622                 << v1 << " of collapsed edge" << abort(FatalError);
00623         }
00624     }
00625     return cosAngle;
00626 }
00627 
00628 
00629 //- Calculate minimum (cos of) edge angle using addressing from collapsing
00630 //  edge to v1 at pt.
00631 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
00632 (
00633     const triSurface& surf,
00634     const label v1,
00635     const point& pt,
00636     const labelHashSet& collapsedFaces,
00637     const HashTable<label, label, Hash<label> >& edgeToEdge,
00638     const HashTable<label, label, Hash<label> >& edgeToFace
00639 )
00640 {
00641     const labelList& v1Faces = surf.pointFaces()[v1];
00642 
00643     scalar minCos = 1;
00644 
00645     forAll(v1Faces, v1FaceI)
00646     {
00647         label faceI = v1Faces[v1FaceI];
00648 
00649         if (collapsedFaces.found(faceI))
00650         {
00651             continue;
00652         }
00653 
00654         const labelList& myEdges = surf.faceEdges()[faceI];
00655 
00656         forAll(myEdges, myEdgeI)
00657         {
00658             label edgeI = myEdges[myEdgeI];
00659 
00660             minCos =
00661                 min
00662                 (
00663                     minCos,
00664                     edgeCosAngle
00665                     (
00666                         surf,
00667                         v1,
00668                         pt,
00669                         collapsedFaces,
00670                         edgeToEdge,
00671                         edgeToFace,
00672                         faceI,
00673                         edgeI
00674                     )
00675                 );
00676         }
00677     }
00678 
00679     return minCos;
00680 }
00681 
00682 
00683 // Calculate max dihedral angle after collapsing edge to v1 (at pt).
00684 // Note that all edges of all faces using v1 are affected.
00685 bool Foam::triSurfaceTools::collapseCreatesFold
00686 (
00687     const triSurface& surf,
00688     const label v1,
00689     const point& pt,
00690     const labelHashSet& collapsedFaces,
00691     const HashTable<label, label, Hash<label> >& edgeToEdge,
00692     const HashTable<label, label, Hash<label> >& edgeToFace,
00693     const scalar minCos
00694 )
00695 {
00696     const labelList& v1Faces = surf.pointFaces()[v1];
00697 
00698     forAll(v1Faces, v1FaceI)
00699     {
00700         label faceI = v1Faces[v1FaceI];
00701 
00702         if (collapsedFaces.found(faceI))
00703         {
00704             continue;
00705         }
00706 
00707         const labelList& myEdges = surf.faceEdges()[faceI];
00708 
00709         forAll(myEdges, myEdgeI)
00710         {
00711             label edgeI = myEdges[myEdgeI];
00712 
00713             if
00714             (
00715                 edgeCosAngle
00716                 (
00717                     surf,
00718                     v1,
00719                     pt,
00720                     collapsedFaces,
00721                     edgeToEdge,
00722                     edgeToFace,
00723                     faceI,
00724                     edgeI
00725                 )
00726               < minCos
00727             )
00728             {
00729                 return true;
00730             }
00731         }
00732     }
00733 
00734     return false;
00735 }
00736 
00737 
00740 // a vertex.
00741 //bool Foam::triSurfaceTools::collapseCreatesDuplicates
00742 //(
00743 //    const triSurface& surf,
00744 //    const label edgeI,
00745 //    const labelHashSet& collapsedFaces
00746 //)
00747 //{
00748 //Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI]
00749 //    << " collapsedFaces:" << collapsedFaces.toc() << endl;
00750 //
00751 //    // Mark neighbours of faces to be collapsed, i.e. get the first layer
00752 //    // of triangles outside collapsedFaces.
00753 //    // neighbours actually contains the
00754 //    // edge with which triangle connects to collapsedFaces.
00755 //
00756 //    HashTable<label, label, Hash<label> > neighbours;
00757 //
00758 //    labelList collapsed = collapsedFaces.toc();
00759 //
00760 //    forAll(collapsed, collapseI)
00761 //    {
00762 //        const label faceI = collapsed[collapseI];
00763 //
00764 //        const labelList& myEdges = surf.faceEdges()[faceI];
00765 //
00766 //        Pout<< "collapsing faceI:" << faceI << " uses edges:" << myEdges
00767 //            << endl;
00768 //
00769 //        forAll(myEdges, myEdgeI)
00770 //        {
00771 //            const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]];
00772 //
00773 //            Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces "
00774 //                << myFaces << endl;
00775 //
00776 //            if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2))
00777 //            {
00778 //                // Get the other face
00779 //                label neighbourFaceI = myFaces[0];
00780 //
00781 //                if (neighbourFaceI == faceI)
00782 //                {
00783 //                    neighbourFaceI = myFaces[1];
00784 //                }
00785 //
00786 //                // Is 'outside' face if not in collapsedFaces itself
00787 //                if (!collapsedFaces.found(neighbourFaceI))
00788 //                {
00789 //                    neighbours.insert(neighbourFaceI, myEdges[myEdgeI]);
00790 //                }
00791 //            }
00792 //        }
00793 //    }
00794 //
00795 //    // Now neighbours contains first layer of triangles outside of
00796 //    // collapseFaces
00797 //    // There should be 
00798 //    // -two if edgeI is a boundary edge
00799 //    // since the outside 'edge' of collapseFaces should
00800 //    // form a triangle and the face connected to edgeI is not inserted.
00801 //    // -four if edgeI is not a boundary edge since then the outside edge forms
00802 //    // a diamond.
00803 //
00804 //    // Check if any of the faces in neighbours share an edge. (n^2)
00805 //
00806 //    labelList neighbourList = neighbours.toc();
00807 //
00808 //Pout<< "edgeI:" << edgeI << "  neighbourList:" << neighbourList << endl;
00809 //
00810 //
00811 //    forAll(neighbourList, i)
00812 //    {
00813 //        const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]];
00814 //
00815 //        for (label j = i+1; j < neighbourList.size(); i++)
00816 //        {
00817 //            const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
00818 //            
00819 //            // Check if faceI and faceJ share an edge
00820 //            forAll(faceIEdges, fI)
00821 //            {
00822 //                forAll(faceJEdges, fJ)
00823 //                {
00824 //                    Pout<< " comparing " << faceIEdges[fI] << " to "
00825 //                        << faceJEdges[fJ] << endl;
00826 //
00827 //                    if (faceIEdges[fI] == faceJEdges[fJ])
00828 //                    {
00829 //                        return true;
00830 //                    }
00831 //                }
00832 //            }
00833 //        }
00834 //    }
00835 //    Pout<< "Found no match. Returning false" << endl;
00836 //    return false;
00837 //}
00838 
00839 
00840 // Finds the triangle edge cut by the plane between a point inside / on edge
00841 // of a triangle and a point outside. Returns:
00842 //  - cut through edge/point
00843 //  - miss()
00844 Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
00845 (
00846     const triSurface& s,
00847     const label triI,
00848     const label excludeEdgeI,
00849     const label excludePointI,
00850 
00851     const point& triPoint,
00852     const plane& cutPlane,
00853     const point& toPoint
00854 )
00855 {
00856     const pointField& points = s.points();
00857     const labelledTri& f = s[triI];
00858     const labelList& fEdges = s.faceEdges()[triI];
00859 
00860     // Get normal distance to planeN
00861     FixedList<scalar, 3> d;
00862 
00863     scalar norm = 0;
00864     forAll(d, fp)
00865     {
00866         d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
00867         norm += mag(d[fp]);
00868     }
00869 
00870     // Normalise and truncate
00871     forAll(d, i)
00872     {
00873         d[i] /= norm;
00874 
00875         if (mag(d[i]) < 1E-6)
00876         {
00877             d[i] = 0.0;
00878         }
00879     }
00880 
00881     // Return information
00882     surfaceLocation cut;
00883 
00884     if (excludePointI != -1)
00885     {
00886         // Excluded point. Test only opposite edge.
00887 
00888         label fp0 = findIndex(s.localFaces()[triI], excludePointI);
00889 
00890         if (fp0 == -1)
00891         {
00892             FatalErrorIn("cutEdge(..)") << "excludePointI:" << excludePointI
00893                 << " localF:" << s.localFaces()[triI] << abort(FatalError);
00894         }
00895 
00896         label fp1 = f.fcIndex(fp0);
00897         label fp2 = f.fcIndex(fp1);
00898 
00899 
00900         if (d[fp1] == 0.0)
00901         {
00902             cut.setHit();
00903             cut.setPoint(points[f[fp1]]);
00904             cut.elementType() = triPointRef::POINT;
00905             cut.setIndex(s.localFaces()[triI][fp1]);
00906         }
00907         else if (d[fp2] == 0.0)
00908         {
00909             cut.setHit();
00910             cut.setPoint(points[f[fp2]]);
00911             cut.elementType() = triPointRef::POINT;
00912             cut.setIndex(s.localFaces()[triI][fp2]);
00913         }
00914         else if
00915         (
00916             (d[fp1] < 0 && d[fp2] < 0)
00917          || (d[fp1] > 0 && d[fp2] > 0)
00918         )
00919         {
00920             // Both same sign. Not crossing edge at all.
00921             // cut already set to miss().
00922         }
00923         else
00924         {
00925             cut.setHit();
00926             cut.setPoint
00927             (
00928                 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
00929               / (d[fp2] - d[fp1])
00930             );
00931             cut.elementType() = triPointRef::EDGE;
00932             cut.setIndex(fEdges[fp1]);
00933         }
00934     }
00935     else
00936     {
00937         // Find the two intersections
00938         FixedList<surfaceLocation, 2> inters;
00939         label interI = 0;
00940 
00941         forAll(f, fp0)
00942         {
00943             label fp1 = f.fcIndex(fp0);
00944 
00945             if (d[fp0] == 0)
00946             {
00947                 if (interI >= 2)
00948                 {
00949                     FatalErrorIn("cutEdge(..)")
00950                         << "problem : triangle has three intersections." << nl
00951                         << "triangle:" << f.tri(points)
00952                         << " d:" << d << abort(FatalError);
00953                 }
00954                 inters[interI].setHit();
00955                 inters[interI].setPoint(points[f[fp0]]);
00956                 inters[interI].elementType() = triPointRef::POINT;
00957                 inters[interI].setIndex(s.localFaces()[triI][fp0]);
00958                 interI++;
00959             }
00960             else if
00961             (
00962                 (d[fp0] < 0 && d[fp1] > 0)
00963              || (d[fp0] > 0 && d[fp1] < 0)
00964             )
00965             {
00966                 if (interI >= 2)
00967                 {
00968                     FatalErrorIn("cutEdge(..)")
00969                         << "problem : triangle has three intersections." << nl
00970                         << "triangle:" << f.tri(points)
00971                         << " d:" << d << abort(FatalError);
00972                 }
00973                 inters[interI].setHit();
00974                 inters[interI].setPoint
00975                 (
00976                     (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
00977                   / (d[fp0] - d[fp1])
00978                 );
00979                 inters[interI].elementType() = triPointRef::EDGE;
00980                 inters[interI].setIndex(fEdges[fp0]);
00981                 interI++;
00982             }
00983         }
00984 
00985 
00986         if (interI == 0)
00987         {
00988             // Return miss
00989         }
00990         else if (interI == 1)
00991         {
00992             // Only one intersection. Should not happen!
00993             cut = inters[0];
00994         }
00995         else if (interI == 2)
00996         {
00997             // Handle excludeEdgeI
00998             if
00999             (
01000                 inters[0].elementType() == triPointRef::EDGE
01001              && inters[0].index() == excludeEdgeI
01002             )
01003             {
01004                 cut = inters[1];
01005             }
01006             else if
01007             (
01008                 inters[1].elementType() == triPointRef::EDGE
01009              && inters[1].index() == excludeEdgeI
01010             )
01011             {
01012                 cut = inters[0];
01013             }
01014             else
01015             {
01016                 // Two cuts. Find nearest.
01017                 if
01018                 (
01019                     magSqr(inters[0].rawPoint() - toPoint)
01020                   < magSqr(inters[1].rawPoint() - toPoint)
01021                 )
01022                 {
01023                     cut = inters[0];
01024                 }
01025                 else
01026                 {
01027                     cut = inters[1];
01028                 }
01029             }
01030         }
01031     }
01032     return cut;
01033 }
01034 
01035 
01036 // 'Snap' point on to endPoint.
01037 void Foam::triSurfaceTools::snapToEnd
01038 (
01039     const triSurface& s,
01040     const surfaceLocation& end,
01041     surfaceLocation& current
01042 )
01043 {
01044     if (end.elementType() == triPointRef::NONE)
01045     {
01046         if (current.elementType() == triPointRef::NONE)
01047         {
01048             // endpoint on triangle; current on triangle
01049             if (current.index() == end.index())
01050             {
01051                 //if (debug)
01052                 //{
01053                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01054                 //        << end << endl;
01055                 //}
01056                 current = end;
01057                 current.setHit();
01058             }
01059         }
01060         // No need to handle current on edge/point since tracking handles this.
01061     }
01062     else if (end.elementType() == triPointRef::EDGE)
01063     {
01064         if (current.elementType() == triPointRef::NONE)
01065         {
01066             // endpoint on edge; current on triangle
01067             const labelList& fEdges = s.faceEdges()[current.index()];
01068 
01069             if (findIndex(fEdges, end.index()) != -1)
01070             {
01071                 //if (debug)
01072                 //{
01073                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01074                 //        << end << endl;
01075                 //}
01076                 current = end;
01077                 current.setHit();
01078             }
01079         }
01080         else if (current.elementType() == triPointRef::EDGE)
01081         {
01082             // endpoint on edge; current on edge
01083             if (current.index() == end.index())
01084             {
01085                 //if (debug)
01086                 //{
01087                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01088                 //        << end << endl;
01089                 //}
01090                 current = end;
01091                 current.setHit();
01092             }
01093         }
01094         else
01095         {
01096             // endpoint on edge; current on point
01097             const edge& e = s.edges()[end.index()];
01098 
01099             if (current.index() == e[0] || current.index() == e[1])
01100             {
01101                 //if (debug)
01102                 //{
01103                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01104                 //        << end << endl;
01105                 //}
01106                 current = end;
01107                 current.setHit();
01108             }
01109         }
01110     }
01111     else    // end.elementType() == POINT
01112     {
01113         if (current.elementType() == triPointRef::NONE)
01114         {
01115             // endpoint on point; current on triangle
01116             const labelledTri& f = s.localFaces()[current.index()];
01117 
01118             if (findIndex(f, end.index()) != -1)
01119             {
01120                 //if (debug)
01121                 //{
01122                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01123                 //        << end << endl;
01124                 //}
01125                 current = end;
01126                 current.setHit();
01127             }
01128         }
01129         else if (current.elementType() == triPointRef::EDGE)
01130         {
01131             // endpoint on point; current on edge
01132             const edge& e = s.edges()[current.index()];
01133 
01134             if (end.index() == e[0] || end.index() == e[1])
01135             {
01136                 //if (debug)
01137                 //{
01138                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01139                 //        << end << endl;
01140                 //}
01141                 current = end;
01142                 current.setHit();
01143             }
01144         }
01145         else
01146         {
01147             // endpoint on point; current on point
01148             if (current.index() == end.index())
01149             {
01150                 //if (debug)
01151                 //{
01152                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
01153                 //        << end << endl;
01154                 //}
01155                 current = end;
01156                 current.setHit();
01157             }
01158         }
01159     }
01160 }
01161 
01162 
01163 // Start:
01164 //  - location
01165 //  - element type (triangle/edge/point) and index
01166 //  - triangle to exclude
01167 Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
01168 (
01169     const triSurface& s,
01170     const labelList& eFaces,
01171     const surfaceLocation& start,
01172     const label excludeEdgeI,
01173     const label excludePointI,
01174     const surfaceLocation& end,
01175     const plane& cutPlane
01176 )
01177 {
01178     surfaceLocation nearest;
01179 
01180     scalar minDistSqr = Foam::sqr(GREAT);
01181 
01182     forAll(eFaces, i)
01183     {
01184         label triI = eFaces[i];
01185 
01186         // Make sure we don't revisit previous face
01187         if (triI != start.triangle())
01188         {
01189             if (end.elementType() == triPointRef::NONE && end.index() == triI)
01190             {
01191                 // Endpoint is in this triangle. Jump there.
01192                 nearest = end;
01193                 nearest.setHit();
01194                 nearest.triangle() = triI;
01195                 break;
01196             }            
01197             else
01198             {
01199                // Which edge is cut.
01200 
01201                 surfaceLocation cutInfo = cutEdge
01202                 (
01203                     s,
01204                     triI,
01205                     excludeEdgeI,       // excludeEdgeI
01206                     excludePointI,      // excludePointI
01207                     start.rawPoint(),
01208                     cutPlane,
01209                     end.rawPoint()
01210                 );
01211 
01212                 // If crossing an edge we expect next edge to be cut.
01213                 if (excludeEdgeI != -1 && !cutInfo.hit())
01214                 {
01215                     FatalErrorIn("triSurfaceTools::visitFaces(..)")
01216                         << "Triangle:" << triI
01217                         << " excludeEdge:" << excludeEdgeI
01218                         << " point:" << start.rawPoint()
01219                         << " plane:" << cutPlane
01220                         << " . No intersection!" << abort(FatalError);
01221                 }
01222 
01223                 if (cutInfo.hit())
01224                 {
01225                     scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
01226 
01227                     if (distSqr < minDistSqr)
01228                     {
01229                         minDistSqr = distSqr;
01230                         nearest = cutInfo;
01231                         nearest.triangle() = triI;
01232                         nearest.setMiss();
01233                     }
01234                 }
01235             }
01236         }
01237     }
01238 
01239     if (nearest.triangle() == -1)
01240     {
01241         // Did not move from edge. Give warning? Return something special?
01242         // For now responsability of caller to make sure that nothing has
01243         // moved.
01244     }
01245 
01246     return nearest;
01247 }
01248 
01249 
01250 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
01251 
01252 
01253 // Write pointField to file
01254 void Foam::triSurfaceTools::writeOBJ
01255 (
01256     const fileName& fName,
01257     const pointField& pts
01258 )
01259 {
01260     OFstream outFile(fName);
01261 
01262     forAll(pts, pointI)
01263     {
01264         const point& pt = pts[pointI];
01265 
01266         outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
01267     }
01268     Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
01269 }
01270 
01271 
01272 // Write vertex subset to OBJ format file
01273 void Foam::triSurfaceTools::writeOBJ
01274 (
01275     const triSurface& surf,
01276     const fileName& fName,
01277     const boolList& markedVerts
01278 )
01279 {
01280     OFstream outFile(fName);
01281 
01282     label nVerts = 0;
01283     forAll(markedVerts, vertI)
01284     {
01285         if (markedVerts[vertI])
01286         {
01287             const point& pt = surf.localPoints()[vertI];
01288 
01289             outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
01290 
01291             nVerts++;
01292         }
01293     }
01294     Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
01295 }
01296 
01297 
01298 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
01299 // Addressing helper functions:
01300 
01301 
01302 // Get all triangles using vertices of edge
01303 void Foam::triSurfaceTools::getVertexTriangles
01304 (
01305     const triSurface& surf,
01306     const label edgeI,
01307     labelList& edgeTris
01308 )
01309 {
01310     const edge& e = surf.edges()[edgeI];
01311     const labelList& myFaces = surf.edgeFaces()[edgeI];
01312 
01313     label face1I = myFaces[0];
01314     label face2I = -1;
01315     if (myFaces.size() == 2)
01316     {
01317         face2I = myFaces[1];
01318     }
01319 
01320     const labelList& startFaces = surf.pointFaces()[e.start()];
01321     const labelList& endFaces = surf.pointFaces()[e.end()];
01322 
01323     // Number of triangles is sum of pointfaces - common faces
01324     // (= faces using edge)
01325     edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
01326 
01327     label nTris = 0;
01328     forAll(startFaces, startFaceI)
01329     {
01330         edgeTris[nTris++] = startFaces[startFaceI];
01331     }
01332 
01333     forAll(endFaces, endFaceI)
01334     {
01335         label faceI = endFaces[endFaceI];
01336 
01337         if ((faceI != face1I) && (faceI != face2I))
01338         {
01339             edgeTris[nTris++] = faceI;
01340         }
01341     }
01342 }
01343 
01344 
01345 // Get all vertices connected to vertices of edge
01346 Foam::labelList Foam::triSurfaceTools::getVertexVertices
01347 (
01348     const triSurface& surf,
01349     const edge& e
01350 )
01351 {
01352     const edgeList& edges = surf.edges();
01353 
01354     label v1 = e.start();
01355     label v2 = e.end();
01356 
01357     // Get all vertices connected to v1 or v2 through an edge
01358     labelHashSet vertexNeighbours;
01359 
01360     const labelList& v1Edges = surf.pointEdges()[v1];
01361 
01362     forAll(v1Edges, v1EdgeI)
01363     {
01364         const edge& e = edges[v1Edges[v1EdgeI]];
01365         vertexNeighbours.insert(e.otherVertex(v1));
01366     }
01367 
01368     const labelList& v2Edges = surf.pointEdges()[v2];
01369 
01370     forAll(v2Edges, v2EdgeI)
01371     {
01372         const edge& e = edges[v2Edges[v2EdgeI]];
01373 
01374         label vertI = e.otherVertex(v2);
01375 
01376         vertexNeighbours.insert(vertI);
01377     }
01378     return vertexNeighbours.toc();
01379 }
01380 
01381 
01383 //void Foam::triSurfaceTools::orderVertices
01384 //(
01385 //    const labelledTri& f,
01386 //    const label v1,
01387 //    const label v2,
01388 //    label& vA,
01389 //    label& vB
01390 //)
01391 //{
01392 //    // Order v1, v2 in anticlockwise order.
01393 //    bool reverse = false;
01394 //
01395 //    if (f[0] == v1)
01396 //    {
01397 //        if (f[1] != v2)
01398 //        {
01399 //            reverse = true;
01400 //        }
01401 //    }
01402 //    else if (f[1] == v1)
01403 //    {
01404 //        if (f[2] != v2)
01405 //        {
01406 //            reverse = true;
01407 //        }
01408 //    }
01409 //    else
01410 //    {
01411 //        if (f[0] != v2)
01412 //        {
01413 //            reverse = true;
01414 //        }
01415 //    }
01416 //
01417 //    if (reverse)
01418 //    {
01419 //        vA = v2;
01420 //        vB = v1;
01421 //    }
01422 //    else
01423 //    {
01424 //        vA = v1;
01425 //        vB = v2;
01426 //    }
01427 //}
01428 
01429 
01430 // Get the other face using edgeI
01431 Foam::label Foam::triSurfaceTools::otherFace
01432 (
01433     const triSurface& surf,
01434     const label faceI,
01435     const label edgeI
01436 )
01437 {
01438     const labelList& myFaces = surf.edgeFaces()[edgeI];
01439 
01440     if (myFaces.size() != 2)
01441     {
01442         return -1;
01443     }
01444     else
01445     {
01446         if (faceI == myFaces[0])
01447         {
01448             return myFaces[1];
01449         }
01450         else
01451         {
01452             return myFaces[0];
01453         }
01454     }
01455 }
01456 
01457 
01458 // Get the two edges on faceI counterclockwise after edgeI
01459 void Foam::triSurfaceTools::otherEdges
01460 (
01461     const triSurface& surf,
01462     const label faceI,
01463     const label edgeI,
01464     label& e1,
01465     label& e2
01466 )
01467 {
01468     const labelList& eFaces = surf.faceEdges()[faceI];
01469 
01470     label i0 = findIndex(eFaces, edgeI);
01471 
01472     if (i0 == -1)
01473     {
01474         FatalErrorIn
01475         (
01476             "otherEdges"
01477             "(const triSurface&, const label, const label,"
01478             " label&, label&)"
01479         )   << "Edge " << surf.edges()[edgeI] << " not in face "
01480             << surf.localFaces()[faceI] << abort(FatalError);
01481     }
01482 
01483     label i1 = eFaces.fcIndex(i0);
01484     label i2 = eFaces.fcIndex(i1);
01485 
01486     e1 = eFaces[i1];
01487     e2 = eFaces[i2];
01488 }
01489 
01490 
01491 // Get the two vertices on faceI counterclockwise vertI
01492 void Foam::triSurfaceTools::otherVertices
01493 (
01494     const triSurface& surf,
01495     const label faceI,
01496     const label vertI,
01497     label& vert1I,
01498     label& vert2I
01499 )
01500 {
01501     const labelledTri& f = surf.localFaces()[faceI];
01502 
01503     if (vertI == f[0])
01504     {
01505         vert1I = f[1];
01506         vert2I = f[2];
01507     }
01508     else if (vertI == f[1])
01509     {
01510         vert1I = f[2];
01511         vert2I = f[0];
01512     }
01513     else if (vertI == f[2])
01514     {
01515         vert1I = f[0];
01516         vert2I = f[1];
01517     }
01518     else
01519     {
01520         FatalErrorIn
01521         (
01522             "otherVertices"
01523             "(const triSurface&, const label, const label,"
01524             " label&, label&)"
01525         )   << "Vertex " << vertI << " not in face " << f << abort(FatalError);
01526     }
01527 }
01528 
01529 
01530 // Get edge opposite vertex
01531 Foam::label Foam::triSurfaceTools::oppositeEdge
01532 (
01533     const triSurface& surf,
01534     const label faceI,
01535     const label vertI
01536 )
01537 {
01538     const labelList& myEdges = surf.faceEdges()[faceI];
01539 
01540     forAll(myEdges, myEdgeI)
01541     {
01542         label edgeI = myEdges[myEdgeI];
01543 
01544         const edge& e = surf.edges()[edgeI];
01545 
01546         if ((e.start() != vertI) && (e.end() != vertI))
01547         {
01548             return edgeI;
01549         }
01550     }
01551 
01552     FatalErrorIn
01553     (
01554         "oppositeEdge"
01555         "(const triSurface&, const label, const label)"
01556     )   << "Cannot find vertex " << vertI << " in edges of face " << faceI
01557         << abort(FatalError);
01558 
01559     return -1;
01560 }
01561 
01562 
01563 // Get vertex opposite edge
01564 Foam::label Foam::triSurfaceTools::oppositeVertex
01565 (
01566     const triSurface& surf,
01567     const label faceI,
01568     const label edgeI
01569 )
01570 {
01571     const labelledTri& f = surf.localFaces()[faceI];
01572 
01573     const edge& e = surf.edges()[edgeI];
01574 
01575     forAll(f, fp)
01576     {
01577         label vertI = f[fp];
01578 
01579         if (vertI != e.start() && vertI != e.end())
01580         {
01581             return vertI;
01582         }
01583     }
01584 
01585     FatalErrorIn("triSurfaceTools::oppositeVertex")
01586         << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
01587         << " in face " << faceI << " vertices " << f << abort(FatalError);
01588 
01589     return -1;
01590 }
01591 
01592 
01593 // Returns edge label connecting v1, v2
01594 Foam::label Foam::triSurfaceTools::getEdge
01595 (
01596     const triSurface& surf,
01597     const label v1,
01598     const label v2
01599 )
01600 {
01601     const labelList& v1Edges = surf.pointEdges()[v1];
01602 
01603     forAll(v1Edges, v1EdgeI)
01604     {
01605         label edgeI = v1Edges[v1EdgeI];
01606 
01607         const edge& e = surf.edges()[edgeI];
01608 
01609         if ((e.start() == v2) || (e.end() == v2))
01610         {
01611             return edgeI;
01612         }
01613     }
01614     return -1;
01615 }
01616 
01617 
01618 // Return index of triangle (or -1) using all three edges
01619 Foam::label Foam::triSurfaceTools::getTriangle
01620 (
01621     const triSurface& surf,
01622     const label e0I,
01623     const label e1I,
01624     const label e2I
01625 )
01626 {
01627     if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
01628     {
01629         FatalErrorIn
01630         (
01631             "getTriangle"
01632             "(const triSurface&, const label, const label,"
01633             " const label)"
01634         )   << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
01635             << " e2:" << e2I
01636             << abort(FatalError);
01637     }
01638 
01639     const labelList& eFaces = surf.edgeFaces()[e0I];
01640 
01641     forAll(eFaces, eFaceI)
01642     {
01643         label faceI = eFaces[eFaceI];
01644 
01645         const labelList& myEdges = surf.faceEdges()[faceI];
01646 
01647         if
01648         (
01649             (myEdges[0] == e1I)
01650          || (myEdges[1] == e1I)
01651          || (myEdges[2] == e1I)
01652         )
01653         {
01654             if
01655             (
01656                 (myEdges[0] == e2I)
01657              || (myEdges[1] == e2I)
01658              || (myEdges[2] == e2I)
01659             )
01660             {
01661                 return faceI;
01662             }
01663         }
01664     }
01665     return -1;
01666 }
01667 
01668 
01669 // Collapse indicated edges. Return new tri surface.
01670 Foam::triSurface Foam::triSurfaceTools::collapseEdges
01671 (
01672     const triSurface& surf,
01673     const labelList& collapsableEdges
01674 )
01675 {
01676     pointField edgeMids(surf.nEdges());
01677 
01678     forAll(edgeMids, edgeI)
01679     {
01680         const edge& e = surf.edges()[edgeI];
01681 
01682         edgeMids[edgeI] =
01683             0.5
01684           * (
01685                 surf.localPoints()[e.start()]
01686               + surf.localPoints()[e.end()]
01687             );
01688     }
01689 
01690 
01691     labelList faceStatus(surf.size(), ANYEDGE);
01692 
01694     //forAll(edges, edgeI)
01695     //{
01696     //    const labelList& neighbours = edgeFaces[edgeI];
01697     //
01698     //    if ((neighbours.size() != 2) && (neighbours.size() != 1))
01699     //    {
01700     //        FatalErrorIn("collapseEdges")
01701     //            << abort(FatalError);
01702     //    }
01703     //
01704     //    if (neighbours.size() == 2)
01705     //    {
01706     //        if (surf[neighbours[0]].region() != surf[neighbours[1]].region())
01707     //        {
01708     //            // Neighbours on different regions. For now dont allow
01709     //            // any collapse.
01710     //            //Pout<< "protecting face " << neighbours[0]
01711     //            //    << ' ' << neighbours[1] << endl;
01712     //            faceStatus[neighbours[0]] = NOEDGE;
01713     //            faceStatus[neighbours[1]] = NOEDGE;
01714     //        }
01715     //    }
01716     //}
01717 
01718     return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
01719 }
01720 
01721 
01722 // Collapse indicated edges. Return new tri surface.
01723 Foam::triSurface Foam::triSurfaceTools::collapseEdges
01724 (
01725     const triSurface& surf,
01726     const labelList& collapseEdgeLabels,
01727     const pointField& edgeMids,
01728     labelList& faceStatus
01729 )
01730 {
01731     const labelListList& edgeFaces = surf.edgeFaces();
01732     const pointField& localPoints = surf.localPoints();
01733     const edgeList& edges = surf.edges();
01734 
01735     // Storage for new points
01736     pointField newPoints(localPoints);
01737 
01738     // Map for old to new points
01739     labelList pointMap(localPoints.size());
01740     forAll(localPoints, pointI)
01741     {
01742         pointMap[pointI] = pointI;
01743     }
01744 
01745 
01746     // Do actual 'collapsing' of edges
01747 
01748     forAll(collapseEdgeLabels, collapseEdgeI)
01749     {
01750         const label edgeI = collapseEdgeLabels[collapseEdgeI];
01751 
01752         if ((edgeI < 0) || (edgeI >= surf.nEdges()))
01753         {
01754             FatalErrorIn("collapseEdges")
01755                 << "Edge label outside valid range." << endl
01756                 << "edge label:" << edgeI << endl
01757                 << "total number of edges:" << surf.nEdges() << endl
01758                 << abort(FatalError);
01759         }
01760 
01761         const labelList& neighbours = edgeFaces[edgeI];
01762 
01763         if (neighbours.size() == 2)
01764         {
01765             const label stat0 = faceStatus[neighbours[0]];
01766             const label stat1 = faceStatus[neighbours[1]];
01767 
01768             // Check faceStatus to make sure this one can be collapsed
01769             if
01770             (
01771                 ((stat0 == ANYEDGE) || (stat0 == edgeI))
01772              && ((stat1 == ANYEDGE) || (stat1 == edgeI))
01773             )
01774             {
01775                 const edge& e = edges[edgeI];
01776 
01777                 // Set up mapping to 'collapse' points of edge
01778                 if
01779                 (
01780                     (pointMap[e.start()] != e.start())
01781                  || (pointMap[e.end()] != e.end())
01782                 )
01783                 {
01784                     FatalErrorIn("collapseEdges")
01785                         << "points already mapped. Double collapse." << endl
01786                         << "edgeI:" << edgeI
01787                         << "  start:" << e.start()
01788                         << "  end:" << e.end()
01789                         << "  pointMap[start]:" << pointMap[e.start()]
01790                         << "  pointMap[end]:" << pointMap[e.end()]
01791                         << abort(FatalError);
01792                 }
01793 
01794                 const label minVert = min(e.start(), e.end());
01795                 pointMap[e.start()] = minVert;
01796                 pointMap[e.end()] = minVert;
01797 
01798                 // Move shared vertex to mid of edge
01799                 newPoints[minVert] = edgeMids[edgeI];
01800 
01801                 // Protect neighbouring faces
01802                 protectNeighbours(surf, e.start(), faceStatus);
01803                 protectNeighbours(surf, e.end(), faceStatus);
01804                 protectNeighbours
01805                 (
01806                     surf,
01807                     oppositeVertex(surf, neighbours[0], edgeI),
01808                     faceStatus
01809                 );
01810                 protectNeighbours
01811                 (
01812                     surf,
01813                     oppositeVertex(surf, neighbours[1], edgeI),
01814                     faceStatus
01815                 );
01816 
01817                 // Mark all collapsing faces
01818                 labelList collapseFaces =
01819                     getCollapsedFaces
01820                     (
01821                         surf,
01822                         edgeI
01823                     ).toc();
01824 
01825                 forAll(collapseFaces, collapseI)
01826                 {
01827                      faceStatus[collapseFaces[collapseI]] = COLLAPSED;
01828                 }
01829             }
01830         }
01831     }
01832 
01833 
01834     // Storage for new triangles
01835     List<labelledTri> newTris(surf.size());
01836     label newTriI = 0;
01837 
01838     const List<labelledTri>& localFaces = surf.localFaces();
01839 
01840 
01841     // Get only non-collapsed triangles and renumber vertex labels.
01842     forAll(localFaces, faceI)
01843     {
01844         const labelledTri& f = localFaces[faceI];
01845 
01846         const label a = pointMap[f[0]];
01847         const label b = pointMap[f[1]];
01848         const label c = pointMap[f[2]];
01849 
01850         if
01851         (
01852             (a != b) && (a != c) && (b != c)
01853          && (faceStatus[faceI] != COLLAPSED)
01854         )
01855         {
01856             // uncollapsed triangle
01857             newTris[newTriI++] = labelledTri(a, b, c, f.region());
01858         }
01859         else
01860         {
01861             //Pout<< "Collapsed triangle " << faceI
01862             //    << " vertices:" << f << endl;
01863         }
01864     }
01865     newTris.setSize(newTriI);
01866 
01867 
01868 
01869     // Pack faces
01870 
01871     triSurface tempSurf(newTris, surf.patches(), newPoints);
01872 
01873     return
01874         triSurface
01875         (
01876             tempSurf.localFaces(),
01877             tempSurf.patches(),
01878             tempSurf.localPoints()
01879         );
01880 }
01881 
01882 
01883 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
01884 
01885 Foam::triSurface Foam::triSurfaceTools::redGreenRefine
01886 (
01887     const triSurface& surf,
01888     const labelList& refineFaces
01889 )
01890 {
01891     List<refineType> refineStatus(surf.size(), NONE);
01892 
01893     // Mark & propagate refinement
01894     forAll(refineFaces, refineFaceI)
01895     {
01896         calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus);
01897     }
01898 
01899     // Do actual refinement
01900     return doRefine(surf, refineStatus);
01901 }
01902 
01903 
01904 Foam::triSurface Foam::triSurfaceTools::greenRefine
01905 (
01906     const triSurface& surf,
01907     const labelList& refineEdges
01908 )
01909 {
01910     // Storage for marking faces
01911     List<refineType> refineStatus(surf.size(), NONE);
01912 
01913     // Storage for new faces
01914     DynamicList<labelledTri> newFaces(0);
01915     
01916     pointField newPoints(surf.localPoints());
01917     newPoints.setSize(surf.nPoints() + surf.nEdges());
01918     label newPointI = surf.nPoints();
01919 
01920 
01921     // Refine edges
01922     forAll(refineEdges, refineEdgeI)
01923     {
01924         label edgeI = refineEdges[refineEdgeI];
01925 
01926         const labelList& myFaces = surf.edgeFaces()[edgeI];
01927 
01928         bool neighbourIsRefined= false;
01929 
01930         forAll(myFaces, myFaceI)
01931         {
01932             if (refineStatus[myFaces[myFaceI]] != NONE)
01933             {
01934                 neighbourIsRefined =  true;
01935             }
01936         }
01937 
01938         // Only refine if none of the faces is refined
01939         if (!neighbourIsRefined)
01940         {
01941             // Refine edge
01942             const edge& e = surf.edges()[edgeI];
01943 
01944             point mid =
01945                 0.5
01946               * (
01947                     surf.localPoints()[e.start()]
01948                   + surf.localPoints()[e.end()]
01949                 );
01950 
01951             newPoints[newPointI] = mid;
01952 
01953             // Refine faces using edge
01954             forAll(myFaces, myFaceI)
01955             {
01956                 // Add faces to newFaces
01957                 greenRefine
01958                 (
01959                     surf,
01960                     myFaces[myFaceI],
01961                     edgeI,
01962                     newPointI,
01963                     newFaces
01964                 );
01965 
01966                 // Mark as refined
01967                 refineStatus[myFaces[myFaceI]] = GREEN;
01968             }
01969 
01970             newPointI++;
01971         }
01972     }
01973 
01974     // Add unrefined faces
01975     forAll(surf.localFaces(), faceI)
01976     {
01977         if (refineStatus[faceI] == NONE)
01978         {
01979             newFaces.append(surf.localFaces()[faceI]);
01980         }
01981     }
01982 
01983     newFaces.shrink();
01984     newPoints.setSize(newPointI);
01985 
01986     return triSurface(newFaces, surf.patches(), newPoints, true);
01987 }
01988 
01989 
01990 
01991 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
01992 // Geometric helper functions:
01993 
01994 
01995 // Returns element in edgeIndices with minimum length
01996 Foam::label Foam::triSurfaceTools::minEdge
01997 (
01998     const triSurface& surf,
01999     const labelList& edgeIndices
02000 )
02001 {
02002     scalar minLength = GREAT;
02003     label minIndex = -1;
02004     forAll(edgeIndices, i)
02005     {
02006         const edge& e = surf.edges()[edgeIndices[i]];
02007 
02008         scalar length =
02009             mag
02010             (
02011                 surf.localPoints()[e.end()]
02012               - surf.localPoints()[e.start()]
02013             );
02014 
02015         if (length < minLength)
02016         {
02017             minLength = length;
02018             minIndex = i;
02019         }
02020     }
02021     return edgeIndices[minIndex];
02022 }
02023 
02024 
02025 // Returns element in edgeIndices with maximum length
02026 Foam::label Foam::triSurfaceTools::maxEdge
02027 (
02028     const triSurface& surf,
02029     const labelList& edgeIndices
02030 )
02031 {
02032     scalar maxLength = -GREAT;
02033     label maxIndex = -1;
02034     forAll(edgeIndices, i)
02035     {
02036         const edge& e = surf.edges()[edgeIndices[i]];
02037 
02038         scalar length =
02039             mag
02040             (
02041                 surf.localPoints()[e.end()]
02042               - surf.localPoints()[e.start()]
02043             );
02044 
02045         if (length > maxLength)
02046         {
02047             maxLength = length;
02048             maxIndex = i;
02049         }
02050     }
02051     return edgeIndices[maxIndex];
02052 }
02053 
02054 
02055 // Merge points and reconstruct surface
02056 Foam::triSurface Foam::triSurfaceTools::mergePoints
02057 (
02058     const triSurface& surf,
02059     const scalar mergeTol
02060 )
02061 {
02062     pointField newPoints(surf.nPoints());
02063 
02064     labelList pointMap(surf.nPoints());
02065 
02066     bool hasMerged = Foam::mergePoints
02067     (
02068         surf.localPoints(),
02069         mergeTol,
02070         false,
02071         pointMap,
02072         newPoints
02073     );
02074 
02075     if (hasMerged)
02076     {
02077         // Pack the triangles
02078 
02079         // Storage for new triangles
02080         List<labelledTri> newTriangles(surf.size());
02081         label newTriangleI = 0;
02082 
02083         forAll(surf, faceI)
02084         {
02085             const labelledTri& f = surf.localFaces()[faceI];
02086 
02087             label newA = pointMap[f[0]];
02088             label newB = pointMap[f[1]];
02089             label newC = pointMap[f[2]];
02090 
02091             if ((newA != newB) && (newA != newC) && (newB != newC))
02092             {
02093                 newTriangles[newTriangleI++] =
02094                     labelledTri(newA, newB, newC, f.region());
02095             }
02096         }
02097         newTriangles.setSize(newTriangleI);
02098 
02099         return triSurface
02100         (
02101             newTriangles,
02102             surf.patches(),
02103             newPoints
02104         );
02105     }
02106     else
02107     {
02108         return surf;
02109     }
02110 }
02111 
02112 
02113 // Calculate normal on triangle
02114 Foam::vector Foam::triSurfaceTools::surfaceNormal
02115 (
02116     const triSurface& surf,
02117     const label nearestFaceI,
02118     const point& nearestPt
02119 )
02120 {
02121     const labelledTri& f = surf[nearestFaceI];
02122     const pointField& points = surf.points();
02123 
02124     label nearType;
02125     label nearLabel;
02126     triPointRef
02127     (
02128         points[f[0]],
02129         points[f[1]],
02130         points[f[2]]
02131     ).classify(nearestPt, 1E-6, nearType, nearLabel);
02132 
02133     if (nearType == triPointRef::NONE)
02134     {
02135         // Nearest to face
02136         return surf.faceNormals()[nearestFaceI];
02137     }
02138     else if (nearType == triPointRef::EDGE)
02139     {
02140         // Nearest to edge. Assume order of faceEdges same as face vertices.
02141         label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
02142 
02143         // Calculate edge normal by averaging face normals
02144         const labelList& eFaces = surf.edgeFaces()[edgeI];
02145 
02146         vector edgeNormal(vector::zero);
02147 
02148         forAll(eFaces, i)
02149         {
02150             edgeNormal += surf.faceNormals()[eFaces[i]];
02151         }
02152         return edgeNormal/(mag(edgeNormal) + VSMALL);
02153     }
02154     else
02155     {
02156         // Nearest to point
02157         const labelledTri& localF = surf.localFaces()[nearestFaceI];
02158         return surf.pointNormals()[localF[nearLabel]];
02159     }
02160 }
02161 
02162 
02163 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::edgeSide
02164 (
02165     const triSurface& surf,
02166     const point& sample,
02167     const point& nearestPoint,
02168     const label edgeI
02169 )
02170 {
02171     const labelList& eFaces = surf.edgeFaces()[edgeI];
02172 
02173     if (eFaces.size() != 2)
02174     {
02175         // Surface not closed.
02176         return UNKNOWN;
02177     }
02178     else
02179     {
02180         const vectorField& faceNormals = surf.faceNormals();
02181 
02182         // Compare to bisector. This is actually correct since edge is
02183         // nearest so there is a knife-edge.
02184 
02185         vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
02186 
02187         if (((sample - nearestPoint) & n) > 0)
02188         {
02189             return OUTSIDE;
02190         }
02191         else
02192         {
02193             return INSIDE;
02194         }
02195     }
02196 }
02197 
02198 
02199 // Calculate normal on triangle
02200 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
02201 (
02202     const triSurface& surf,
02203     const point& sample,
02204     const label nearestFaceI,   // nearest face
02205     const point& nearestPoint,  // nearest point on nearest face
02206     const scalar tol
02207 )
02208 {
02209     const labelledTri& f = surf[nearestFaceI];
02210     const pointField& points = surf.points();
02211 
02212     // Find where point is on triangle. Note tolerance needed. Is relative
02213     // one (used in comparing normalized [0..1] triangle coordinates).
02214     label nearType, nearLabel;
02215     triPointRef
02216     (
02217         points[f[0]],
02218         points[f[1]],
02219         points[f[2]]
02220     ).classify(nearestPoint, tol, nearType, nearLabel);
02221 
02222     if (nearType == triPointRef::NONE)
02223     {
02224         // Nearest to face interior. Use faceNormal to determine side
02225         scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
02226 
02227         if (c > 0)
02228         {
02229             return OUTSIDE;
02230         }
02231         else
02232         {
02233             return INSIDE;
02234         }
02235     }
02236     else if (nearType == triPointRef::EDGE)
02237     {
02238         // Nearest to edge nearLabel. Note that this can only be a knife-edge
02239         // situation since otherwise the nearest point could never be the edge.
02240 
02241         // Get the edge. Assume order of faceEdges same as face vertices.
02242         label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
02243 
02244         //if (debug)
02245         //{
02246         //    // Check order of faceEdges same as face vertices.
02247         //    const edge& e = surf.edges()[edgeI];
02248         //    const labelList& meshPoints = surf.meshPoints();
02249         //    const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
02250         //
02251         //    if
02252         //    (
02253         //        meshEdge
02254         //     != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
02255         //    )
02256         //    {
02257         //        FatalErrorIn("triSurfaceTools::surfaceSide")
02258         //            << "Edge:" << edgeI << " local vertices:" << e
02259         //            << " mesh vertices:" << meshEdge
02260         //            << " not at position " << nearLabel
02261         //            << " in face " << f
02262         //            << abort(FatalError);
02263         //    }
02264         //}
02265 
02266         return edgeSide(surf, sample, nearestPoint, edgeI);
02267     }
02268     else
02269     {
02270         // Nearest to point. Could use pointNormal here but is not correct.
02271         // Instead determine which edge using point is nearest and use test
02272         // above (nearType == triPointRef::EDGE).
02273 
02274 
02275         const labelledTri& localF = surf.localFaces()[nearestFaceI];
02276         label nearPointI = localF[nearLabel];
02277 
02278         const edgeList& edges = surf.edges();
02279         const pointField& localPoints = surf.localPoints();
02280         const point& base = localPoints[nearPointI];
02281 
02282         const labelList& pEdges = surf.pointEdges()[nearPointI];
02283 
02284         scalar minDistSqr = Foam::sqr(GREAT);
02285         label minEdgeI = -1;
02286 
02287         forAll(pEdges, i)
02288         {
02289             label edgeI = pEdges[i];
02290 
02291             const edge& e = edges[edgeI];
02292 
02293             label otherPointI = e.otherVertex(nearPointI);
02294 
02295             // Get edge normal.
02296             vector eVec(localPoints[otherPointI] - base);
02297             scalar magEVec = mag(eVec);
02298 
02299             if (magEVec > VSMALL)
02300             {
02301                 eVec /= magEVec;
02302 
02303                 // Get point along vector and determine closest.
02304                 const point perturbPoint = base + eVec;
02305 
02306                 scalar distSqr = Foam::magSqr(sample - perturbPoint);
02307 
02308                 if (distSqr < minDistSqr)
02309                 {
02310                     minDistSqr = distSqr;
02311                     minEdgeI = edgeI;
02312                 }
02313             }
02314         }
02315 
02316         if (minEdgeI == -1)
02317         {
02318             FatalErrorIn("treeDataTriSurface::getSide")
02319                 << "Problem: did not find edge closer than " << minDistSqr
02320                 << abort(FatalError);
02321         }
02322 
02323         return edgeSide(surf, sample, nearestPoint, minEdgeI);
02324     }
02325 }
02326 
02327 
02328 // triangulation of boundaryMesh
02329 Foam::triSurface Foam::triSurfaceTools::triangulate
02330 (
02331     const polyBoundaryMesh& bMesh,
02332     const labelHashSet& includePatches,
02333     const bool verbose
02334 )
02335 {
02336     const polyMesh& mesh = bMesh.mesh();
02337 
02338     // Storage for surfaceMesh. Size estimate.
02339     DynamicList<labelledTri> triangles
02340     (
02341         mesh.nFaces() - mesh.nInternalFaces()
02342     );
02343 
02344     label newPatchI = 0;
02345 
02346     for
02347     (
02348         labelHashSet::const_iterator iter = includePatches.begin();
02349         iter != includePatches.end();
02350         ++iter
02351     )
02352     {
02353         label patchI = iter.key();
02354 
02355         const polyPatch& patch = bMesh[patchI];
02356 
02357         const pointField& points = patch.points();
02358 
02359         label nTriTotal = 0;
02360 
02361         forAll(patch, patchFaceI)
02362         {
02363             const face& f = patch[patchFaceI];
02364 
02365             faceList triFaces(f.nTriangles(points));
02366 
02367             label nTri = 0;
02368 
02369             f.triangles(points, nTri, triFaces);
02370 
02371             forAll(triFaces, triFaceI)
02372             {
02373                 const face& f = triFaces[triFaceI];
02374 
02375                 triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
02376 
02377                 nTriTotal++;
02378             }
02379         }
02380 
02381         if (verbose)
02382         {
02383             Pout<< patch.name() << " : generated " << nTriTotal
02384                 << " triangles from " << patch.size() << " faces with"
02385                 << " new patchid " << newPatchI << endl;
02386         }
02387 
02388         newPatchI++;
02389     }
02390     triangles.shrink();
02391 
02392     // Create globally numbered tri surface
02393     triSurface rawSurface(triangles, mesh.points());
02394 
02395     // Create locally numbered tri surface
02396     triSurface surface
02397     (
02398         rawSurface.localFaces(),
02399         rawSurface.localPoints()
02400     );
02401 
02402     // Add patch names to surface
02403     surface.patches().setSize(newPatchI);
02404 
02405     newPatchI = 0;
02406 
02407     for
02408     (
02409         labelHashSet::const_iterator iter = includePatches.begin();
02410         iter != includePatches.end();
02411         ++iter
02412     )
02413     {
02414         label patchI = iter.key();
02415 
02416         const polyPatch& patch = bMesh[patchI];
02417 
02418         surface.patches()[newPatchI].name() = patch.name();
02419 
02420         surface.patches()[newPatchI].geometricType() = patch.type();
02421 
02422         newPatchI++;
02423     }
02424 
02425     return surface;
02426 }
02427 
02428 
02429 // triangulation of boundaryMesh
02430 Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
02431 (
02432     const polyBoundaryMesh& bMesh,
02433     const labelHashSet& includePatches,
02434     const bool verbose
02435 )
02436 {
02437     const polyMesh& mesh = bMesh.mesh();
02438 
02439     // Storage for new points = meshpoints + face centres.
02440     const pointField& points = mesh.points();
02441     const pointField& faceCentres = mesh.faceCentres();
02442 
02443     pointField newPoints(points.size() + faceCentres.size());
02444 
02445     label newPointI = 0;
02446 
02447     forAll(points, pointI)
02448     {
02449         newPoints[newPointI++] = points[pointI];
02450     }
02451     forAll(faceCentres, faceI)
02452     {
02453         newPoints[newPointI++] = faceCentres[faceI];
02454     }
02455 
02456 
02457     // Count number of faces.
02458     DynamicList<labelledTri> triangles
02459     (
02460         mesh.nFaces() - mesh.nInternalFaces()
02461     );
02462 
02463     label newPatchI = 0;
02464 
02465     for
02466     (
02467         labelHashSet::const_iterator iter = includePatches.begin();
02468         iter != includePatches.end();
02469         ++iter
02470     )
02471     {
02472         label patchI = iter.key();
02473 
02474         const polyPatch& patch = bMesh[patchI];
02475 
02476         label nTriTotal = 0;
02477 
02478         forAll(patch, patchFaceI)
02479         {
02480             // Face in global coords.
02481             const face& f = patch[patchFaceI];
02482 
02483             // Index in newPointI of face centre.
02484             label fc = points.size() + patchFaceI + patch.start();
02485 
02486             forAll(f, fp)
02487             {
02488                 label fp1 = (fp + 1) % f.size();
02489 
02490                 triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchI));
02491 
02492                 nTriTotal++;
02493             }
02494         }
02495 
02496         if (verbose)
02497         {
02498             Pout<< patch.name() << " : generated " << nTriTotal
02499                 << " triangles from " << patch.size() << " faces with"
02500                 << " new patchid " << newPatchI << endl;
02501         }
02502 
02503         newPatchI++;
02504     }
02505     triangles.shrink();
02506 
02507 
02508     // Create globally numbered tri surface
02509     triSurface rawSurface(triangles, newPoints);
02510 
02511     // Create locally numbered tri surface
02512     triSurface surface
02513     (
02514         rawSurface.localFaces(),
02515         rawSurface.localPoints()
02516     );
02517 
02518     // Add patch names to surface
02519     surface.patches().setSize(newPatchI);
02520 
02521     newPatchI = 0;
02522 
02523     for
02524     (
02525         labelHashSet::const_iterator iter = includePatches.begin();
02526         iter != includePatches.end();
02527         ++iter
02528     )
02529     {
02530         label patchI = iter.key();
02531 
02532         const polyPatch& patch = bMesh[patchI];
02533 
02534         surface.patches()[newPatchI].name() = patch.name();
02535 
02536         surface.patches()[newPatchI].geometricType() = patch.type();
02537 
02538         newPatchI++;
02539     }
02540 
02541     return surface;
02542 }
02543 
02544 
02545 Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
02546 {
02547     // Vertices in geompack notation. Note that could probably just use
02548     // pts.begin() if double precision.
02549     List<doubleScalar> geompackVertices(2*pts.size());
02550     label doubleI = 0;
02551     forAll(pts, i)
02552     {
02553         geompackVertices[doubleI++] = pts[i][0];
02554         geompackVertices[doubleI++] = pts[i][1];
02555     }
02556 
02557     // Storage for triangles
02558     int m2 = 3;
02559     List<int> triangle_node(m2*3*pts.size());
02560     List<int> triangle_neighbor(m2*3*pts.size());
02561 
02562     // Triangulate
02563     int nTris = 0;
02564     dtris2
02565     (
02566         pts.size(),
02567         geompackVertices.begin(),
02568         &nTris,
02569         triangle_node.begin(),
02570         triangle_neighbor.begin()
02571     );
02572 
02573     // Trim
02574     triangle_node.setSize(3*nTris);
02575     triangle_neighbor.setSize(3*nTris);
02576 
02577     // Convert to triSurface.
02578     List<labelledTri> faces(nTris);
02579 
02580     forAll(faces, i)
02581     {
02582         faces[i] = labelledTri
02583         (
02584             triangle_node[3*i]-1,
02585             triangle_node[3*i+1]-1,
02586             triangle_node[3*i+2]-1,
02587             0
02588         );
02589     }
02590 
02591     pointField points(pts.size());
02592     forAll(pts, i)
02593     {
02594         points[i][0] = pts[i][0];
02595         points[i][1] = pts[i][1];
02596         points[i][2] = 0.0;
02597     }
02598 
02599     return triSurface(faces, points);
02600 }
02601 
02602 
02604 //#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
02605 //#include <CGAL/Delaunay_triangulation_2.h>
02606 //#include <CGAL/Triangulation_vertex_base_with_info_2.h>
02607 //#include <cassert>
02608 //
02609 //struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
02610 //
02611 //typedef CGAL::Triangulation_vertex_base_with_info_2<Foam::label, K> Vb;
02612 //typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
02613 //typedef CGAL::Delaunay_triangulation_2<K,Tds> Triangulation;
02614 //
02615 //typedef Triangulation::Vertex_handle Vertex_handle;
02616 //typedef Triangulation::Vertex_iterator Vertex_iterator;
02617 //typedef Triangulation::Face_handle Face_handle;
02618 //typedef Triangulation::Finite_faces_iterator Finite_faces_iterator;
02619 //typedef Triangulation::Point Point;
02620 //Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
02621 //{
02622 //    Triangulation T;
02623 //
02624 //    // Insert 2D vertices; building triangulation
02625 //    forAll(pts, i)
02626 //    {
02627 //        const point& pt = pts[i];
02628 //
02629 //        T.insert(Point(pt[0], pt[1]));
02630 //    }
02631 //
02632 //
02633 //    // Number vertices
02634 //    // ~~~~~~~~~~~~~~~
02635 //
02636 //    label vertI = 0;
02637 //    for 
02638 //    (
02639 //        Vertex_iterator it = T.vertices_begin();
02640 //        it != T.vertices_end();
02641 //        ++it
02642 //    )
02643 //    {
02644 //        it->info() = vertI++;
02645 //    }
02646 //
02647 //    // Extract faces
02648 //    // ~~~~~~~~~~~~~
02649 //
02650 //    List<labelledTri> faces(T.number_of_faces());
02651 //
02652 //    label faceI = 0;
02653 //
02654 //    for
02655 //    (
02656 //        Finite_faces_iterator fc = T.finite_faces_begin();
02657 //        fc != T.finite_faces_end();
02658 //        ++fc
02659 //    )
02660 //    {
02661 //        faces[faceI++] = labelledTri
02662 //        (
02663 //            fc->vertex(0)->info(),
02664 //            f[1] = fc->vertex(1)->info(),
02665 //            f[2] = fc->vertex(2)->info()
02666 //        );
02667 //    }
02668 //
02669 //    pointField points(pts.size());
02670 //    forAll(pts, i)
02671 //    {
02672 //        points[i][0] = pts[i][0];
02673 //        points[i][1] = pts[i][1];
02674 //        points[i][2] = 0.0;
02675 //    }
02676 //
02677 //    return triSurface(faces, points);
02678 //}
02679 
02680 
02681 void Foam::triSurfaceTools::calcInterpolationWeights
02682 (
02683     const triPointRef& tri,
02684     const point& p,
02685     FixedList<scalar, 3>& weights
02686 )
02687 {
02688     // calculate triangle edge vectors and triangle face normal
02689     // the 'i':th edge is opposite node i
02690     FixedList<vector, 3> edge;
02691     edge[0] = tri.c()-tri.b();
02692     edge[1] = tri.a()-tri.c();
02693     edge[2] = tri.b()-tri.a();
02694 
02695     vector triangleFaceNormal = edge[1] ^ edge[2];
02696 
02697     // calculate edge normal (pointing inwards)
02698     FixedList<vector, 3> normal;
02699     for(label i=0; i<3; i++)
02700     {
02701         normal[i] = triangleFaceNormal ^ edge[i];
02702         normal[i] /= mag(normal[i]) + VSMALL;
02703     }
02704 
02705     weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
02706     weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
02707     weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
02708 }
02709 
02710 
02711 // Calculate weighting factors from samplePts to triangle it is in.
02712 // Uses linear search.
02713 void Foam::triSurfaceTools::calcInterpolationWeights
02714 (
02715     const triSurface& s,
02716     const pointField& samplePts,
02717     List<FixedList<label, 3> >& allVerts,
02718     List<FixedList<scalar, 3> >& allWeights
02719 )
02720 {
02721     allVerts.setSize(samplePts.size());
02722     allWeights.setSize(samplePts.size());
02723 
02724     const pointField& points = s.points();
02725 
02726     forAll(samplePts, i)
02727     {
02728         const point& samplePt = samplePts[i];
02729 
02730 
02731         FixedList<label, 3>& verts = allVerts[i];
02732         FixedList<scalar, 3>& weights = allWeights[i];
02733 
02734         scalar minDistance = GREAT;
02735 
02736         forAll(s, faceI)
02737         {
02738             const labelledTri& f = s[faceI];
02739 
02740             triPointRef tri(f.tri(points));
02741 
02742             pointHit nearest = tri.nearestPoint(samplePt);
02743 
02744             if (nearest.hit())
02745             {
02746                 // samplePt inside triangle
02747                 verts[0] = f[0];
02748                 verts[1] = f[1];
02749                 verts[2] = f[2];
02750 
02751                 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
02752 
02753                 //Pout<< "calcScalingFactors : samplePt:" << samplePt
02754                 //    << " inside triangle:" << faceI
02755                 //    << " verts:" << verts
02756                 //    << " weights:" << weights
02757                 //    << endl;
02758 
02759                 break;
02760             }
02761             else if (nearest.distance() < minDistance)
02762             {
02763                 minDistance = nearest.distance();
02764 
02765                 // Outside triangle. Store nearest.
02766                 label nearType, nearLabel;
02767                 tri.classify
02768                 (
02769                     nearest.rawPoint(),
02770                     1E-6,               // relative tolerance
02771                     nearType,
02772                     nearLabel
02773                 );
02774 
02775                 if (nearType == triPointRef::POINT)
02776                 {
02777                     verts[0] = f[nearLabel];
02778                     weights[0] = 1;
02779                     verts[1] = -1;
02780                     weights[1] = -GREAT;
02781                     verts[2] = -1;
02782                     weights[2] = -GREAT;
02783 
02784                     //Pout<< "calcScalingFactors : samplePt:" << samplePt
02785                     //    << " distance:" << nearest.distance()
02786                     //    << " from point:" << points[f[nearLabel]]
02787                     //    << endl;
02788                 }
02789                 else if (nearType == triPointRef::EDGE)
02790                 {
02791                     verts[0] = f[nearLabel];
02792                     verts[1] = f[f.fcIndex(nearLabel)];
02793                     verts[2] = -1;
02794 
02795                     const point& p0 = points[verts[0]];
02796                     const point& p1 = points[verts[1]];
02797 
02798                     scalar s = min
02799                     (
02800                         1,
02801                         max
02802                         (
02803                             0,
02804                             mag(nearest.rawPoint()-p0)/mag(p1-p0)
02805                         )
02806                     );
02807 
02808                     // Interpolate
02809                     weights[0] = 1-s;
02810                     weights[1] = s;
02811                     weights[2] = -GREAT;
02812 
02813                     //Pout<< "calcScalingFactors : samplePt:" << samplePt
02814                     //    << " distance:" << nearest.distance()
02815                     //    << " from edge:" << p0 << p1 << " s:" << s
02816                     //    << endl;
02817                 }
02818                 else
02819                 {
02820                     // triangle. Can only happen because of truncation errors.
02821                     verts[0] = f[0];
02822                     verts[1] = f[1];
02823                     verts[2] = f[2];
02824 
02825                     calcInterpolationWeights(tri, nearest.rawPoint(), weights);
02826 
02827                     break;
02828                 }
02829             }
02830         }
02831     }
02832 }
02833 
02834 
02835 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
02836 // Tracking:
02837 
02838 
02839 // Test point on surface to see if is on face,edge or point.
02840 Foam::surfaceLocation Foam::triSurfaceTools::classify
02841 (
02842     const triSurface& s,
02843     const label triI,
02844     const point& trianglePoint
02845 )
02846 {
02847     surfaceLocation nearest;
02848 
02849     // Nearest point could be on point or edge. Retest.
02850     label index, elemType;
02851     //bool inside = 
02852     triPointRef(s[triI].tri(s.points())).classify
02853     (
02854         trianglePoint,
02855         1E-6,
02856         elemType,
02857         index
02858     );
02859 
02860     nearest.setPoint(trianglePoint);
02861 
02862     if (elemType == triPointRef::NONE)
02863     {
02864         nearest.setHit();
02865         nearest.setIndex(triI);
02866         nearest.elementType() = triPointRef::NONE;
02867     }    
02868     else if (elemType == triPointRef::EDGE)
02869     {
02870         nearest.setMiss();
02871         nearest.setIndex(s.faceEdges()[triI][index]);
02872         nearest.elementType() = triPointRef::EDGE;
02873     }
02874     else //if (elemType == triPointRef::POINT)
02875     {
02876         nearest.setMiss();
02877         nearest.setIndex(s.localFaces()[triI][index]);
02878         nearest.elementType() = triPointRef::POINT;
02879     }
02880 
02881     return nearest;
02882 }
02883 
02884 
02885 Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge
02886 (
02887     const triSurface& s,
02888     const surfaceLocation& start,
02889     const surfaceLocation& end,
02890     const plane& cutPlane
02891 )
02892 {
02893     // Start off from starting point
02894     surfaceLocation nearest = start;
02895     nearest.setMiss();
02896 
02897     // See if in same triangle as endpoint. If so snap.
02898     snapToEnd(s, end, nearest);
02899 
02900     if (!nearest.hit())
02901     {
02902         // Not yet at end point
02903 
02904         if (start.elementType() == triPointRef::NONE)
02905         {
02906             // Start point is inside triangle. Trivial cases already handled
02907             // above.
02908 
02909             // end point is on edge or point so cross currrent triangle to
02910             // see which edge is cut.
02911 
02912             nearest = cutEdge
02913             (
02914                 s,
02915                 start.index(),          // triangle
02916                 -1,                     // excludeEdge
02917                 -1,                     // excludePoint
02918                 start.rawPoint(),
02919                 cutPlane,
02920                 end.rawPoint()
02921             );
02922             nearest.elementType() = triPointRef::EDGE;
02923             nearest.triangle() = start.index();
02924             nearest.setMiss();
02925         }
02926         else if (start.elementType() == triPointRef::EDGE)
02927         {
02928             // Pick connected triangle that is most in direction.
02929             const labelList& eFaces = s.edgeFaces()[start.index()];
02930 
02931             nearest = visitFaces
02932             (
02933                 s,
02934                 eFaces,
02935                 start,
02936                 start.index(),      // excludeEdgeI
02937                 -1,                 // excludePointI
02938                 end,
02939                 cutPlane
02940             );
02941         }
02942         else    // start.elementType() == triPointRef::POINT
02943         {
02944             const labelList& pFaces = s.pointFaces()[start.index()];
02945 
02946             nearest = visitFaces
02947             (
02948                 s,
02949                 pFaces,
02950                 start,
02951                 -1,                 // excludeEdgeI
02952                 start.index(),      // excludePointI
02953                 end,
02954                 cutPlane
02955             );
02956         }
02957         snapToEnd(s, end, nearest);
02958     }
02959     return nearest;
02960 }
02961 
02962 
02963 void Foam::triSurfaceTools::track
02964 (
02965     const triSurface& s,
02966     const surfaceLocation& endInfo,
02967     const plane& cutPlane,
02968     surfaceLocation& hitInfo
02969 )
02970 {
02971     //OFstream str("track.obj");
02972     //label vertI = 0;
02973     //meshTools::writeOBJ(str, hitInfo.rawPoint());
02974     //vertI++;
02975 
02976     // Track across surface.
02977     while (true)
02978     {
02979         //Pout<< "Tracking from:" << nl
02980         //    << "    " << hitInfo.info()
02981         //    << endl;
02982 
02983         hitInfo = trackToEdge
02984         (
02985             s,
02986             hitInfo,
02987             endInfo,
02988             cutPlane
02989         );
02990 
02991         //meshTools::writeOBJ(str, hitInfo.rawPoint());
02992         //vertI++;
02993         //str<< "l " << vertI-1 << ' ' << vertI << nl;
02994 
02995         //Pout<< "Tracked to:" << nl
02996         //    << "    " << hitInfo.info() << endl;
02997 
02998         if (hitInfo.hit() || hitInfo.triangle() == -1)
02999         {
03000             break;
03001         }
03002     }
03003 }
03004 
03005 
03006 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines