src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "triSurfaceTools.H"
  30
  31 #include "triSurface.H"
  32 #include "OFstream.H"
  33 #include "mergePoints.H"
  34 #include "polyMesh.H"
  35 #include "plane.H"
  36 #include "geompack.H"
  37
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
  42 const Foam::label Foam::triSurfaceTools::NOEDGE = -2;
  43 const Foam::label Foam::triSurfaceTools::COLLAPSED = -3;
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 /*
  48     Refine by splitting all three edges of triangle ('red' refinement).
  49     Neighbouring triangles (which are not marked for refinement get split
  50     in half ('green') refinement. (R. Verfuerth, "A review of a posteriori
  51     error estimation and adaptive mesh refinement techniques",
  52     Wiley-Teubner, 1996)
  53 */
  54
  55 // FaceI gets refined ('red'). Propagate edge refinement.
  56 void Foam::triSurfaceTools::calcRefineStatus
  57 (
  58     const triSurface& surf,
  59     const label faceI,
  60     List<refineType>& refine
  61 )
  62 {
  63     if (refine[faceI] == RED)
  64     {
  65         // Already marked for refinement. Do nothing.
  66     }
  67     else
  68     {
  69         // Not marked or marked for 'green' refinement. Refine.
  70         refine[faceI] = RED;
  71
  72         const labelList& myNeighbours = surf.faceFaces()[faceI];
  73
  74         forAll(myNeighbours, myNeighbourI)
  75         {
  76             label neighbourFaceI = myNeighbours[myNeighbourI];
  77
  78             if (refine[neighbourFaceI] == GREEN)
  79             {
  80                 // Change to red refinement and propagate
  81                 calcRefineStatus(surf, neighbourFaceI, refine);
  82             }
  83             else if (refine[neighbourFaceI] == NONE)
  84             {
  85                 refine[neighbourFaceI] = GREEN;
  86             }
  87         }
  88     }
  89 }
  90
  91
  92 // Split faceI along edgeI at position newPointI
  93 void Foam::triSurfaceTools::greenRefine
  94 (
  95     const triSurface& surf,
  96     const label faceI,
  97     const label edgeI,
  98     const label newPointI,
  99     DynamicList<labelledTri>& newFaces
 100 )
 101 {
 102     const labelledTri& f = surf.localFaces()[faceI];
 103     const edge& e = surf.edges()[edgeI];
 104
 105     // Find index of edge in face.
 106
 107     label fp0 = findIndex(f, e[0]);
 108     label fp1 = f.fcIndex(fp0);
 109     label fp2 = f.fcIndex(fp1);
 110
 111     if (f[fp1] == e[1])
 112     {
 113         // Edge oriented like face
 114         newFaces.append
 115         (
 116             labelledTri
 117             (
 118                 f[fp0],
 119                 newPointI,
 120                 f[fp2],
 121                 f.region()
 122             )
 123         );
 124         newFaces.append
 125         (
 126             labelledTri
 127             (
 128                 newPointI,
 129                 f[fp1],
 130                 f[fp2],
 131                 f.region()
 132             )
 133         );
 134     }
 135     else
 136     {
 137         newFaces.append
 138         (
 139             labelledTri
 140             (
 141                 f[fp2],
 142                 newPointI,
 143                 f[fp1],
 144                 f.region()
 145             )
 146         );
 147         newFaces.append
 148         (
 149             labelledTri
 150             (
 151                 newPointI,
 152                 f[fp0],
 153                 f[fp1],
 154                 f.region()
 155             )
 156         );
 157     }
 158 }
 159
 160
 161 // Refine all triangles marked for refinement.
 162 Foam::triSurface Foam::triSurfaceTools::doRefine
 163 (
 164     const triSurface& surf,
 165     const List<refineType>& refineStatus
 166 )
 167 {
 168     // Storage for new points. (start after old points)
 169     DynamicList<point> newPoints(surf.nPoints());
 170     forAll(surf.localPoints(), pointI)
 171     {
 172         newPoints.append(surf.localPoints()[pointI]);
 173     }
 174     label newVertI = surf.nPoints();
 175
 176     // Storage for new faces
 177     DynamicList<labelledTri> newFaces(surf.size());
 178
 179
 180     // Point index for midpoint on edge
 181     labelList edgeMid(surf.nEdges(), -1);
 182
 183     forAll(refineStatus, faceI)
 184     {
 185         if (refineStatus[faceI] == RED)
 186         {
 187             // Create new vertices on all edges to be refined.
 188             const labelList& fEdges = surf.faceEdges()[faceI];
 189
 190             forAll(fEdges, i)
 191             {
 192                 label edgeI = fEdges[i];
 193
 194                 if (edgeMid[edgeI] == -1)
 195                 {
 196                     const edge& e = surf.edges()[edgeI];
 197
 198                     // Create new point on mid of edge
 199                     newPoints.append
 200                     (
 201                         0.5
 202                       * (
 203                             surf.localPoints()[e.start()]
 204                           + surf.localPoints()[e.end()]
 205                         )
 206                     );
 207                     edgeMid[edgeI] = newVertI++;
 208                 }
 209             }
 210
 211             // Now we have new mid edge vertices for all edges on face
 212             // so create triangles for RED rerfinement.
 213
 214             const edgeList& edges = surf.edges();
 215
 216             // Corner triangles
 217             newFaces.append
 218             (
 219                 labelledTri
 220                 (
 221                     edgeMid[fEdges[0]],
 222                     edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
 223                     edgeMid[fEdges[1]],
 224                     surf[faceI].region()
 225                 )
 226             );
 227
 228             newFaces.append
 229             (
 230                 labelledTri
 231                 (
 232                     edgeMid[fEdges[1]],
 233                     edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
 234                     edgeMid[fEdges[2]],
 235                     surf[faceI].region()
 236                 )
 237             );
 238
 239             newFaces.append
 240             (
 241                 labelledTri
 242                 (
 243                     edgeMid[fEdges[2]],
 244                     edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
 245                     edgeMid[fEdges[0]],
 246                     surf[faceI].region()
 247                 )
 248             );
 249
 250             // Inner triangle
 251             newFaces.append
 252             (
 253                 labelledTri
 254                 (
 255                     edgeMid[fEdges[0]],
 256                     edgeMid[fEdges[1]],
 257                     edgeMid[fEdges[2]],
 258                     surf[faceI].region()
 259                 )
 260             );
 261
 262
 263             // Create triangles for GREEN refinement.
 264             forAll(fEdges, i)
 265             {
 266                 const label edgeI = fEdges[i];
 267
 268                 label otherFaceI = otherFace(surf, faceI, edgeI);
 269
 270                 if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN))
 271                 {
 272                     greenRefine
 273                     (
 274                         surf,
 275                         otherFaceI,
 276                         edgeI,
 277                         edgeMid[edgeI],
 278                         newFaces
 279                     );
 280                 }
 281             }
 282         }
 283     }
 284
 285     // Copy unmarked triangles since keep original vertex numbering.
 286     forAll(refineStatus, faceI)
 287     {
 288         if (refineStatus[faceI] == NONE)
 289         {
 290             newFaces.append(surf.localFaces()[faceI]);
 291         }
 292     }
 293
 294     newFaces.shrink();
 295     newPoints.shrink();
 296
 297
 298     // Transfer DynamicLists to straight ones.
 299     pointField allPoints;
 300     allPoints.transfer(newPoints);
 301     newPoints.clear();
 302
 303     List<labelledTri> allFaces;
 304     allFaces.transfer(newFaces);
 305     newFaces.clear();
 306
 307     return triSurface(allFaces, surf.patches(), allPoints);
 308 }
 309
 310
 311 // Angle between two neighbouring triangles,
 312 // angle between shared-edge end points and left and right face end points
 313 Foam::scalar Foam::triSurfaceTools::faceCosAngle
 314 (
 315     const point& pStart,
 316     const point& pEnd,
 317     const point& pLeft,
 318     const point& pRight
 319 )
 320 {
 321     const vector common(pEnd - pStart);
 322     const vector base0(pLeft - pStart);
 323     const vector base1(pRight - pStart);
 324
 325     vector n0(common ^ base0);
 326     n0 /= Foam::mag(n0);
 327
 328     vector n1(base1 ^ common);
 329     n1 /= Foam::mag(n1);
 330
 331     return n0 & n1;
 332 }
 333
 334
 335 // Protect edges around vertex from collapsing.
 336 // Moving vertI to a different position
 337 // - affects obviously the cost of the faces using it
 338 // - but also their neighbours since the edge between the faces changes
 339 void Foam::triSurfaceTools::protectNeighbours
 340 (
 341     const triSurface& surf,
 342     const label vertI,
 343     labelList& faceStatus
 344 )
 345 {
 346 //    const labelList& myFaces = surf.pointFaces()[vertI];
 347 //    forAll(myFaces, i)
 348 //    {
 349 //        label faceI = myFaces[i];
 350 //
 351 //        if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
 352 //        {
 353 //            faceStatus[faceI] = NOEDGE;
 354 //        }
 355 //    }
 356
 357     const labelList& myEdges = surf.pointEdges()[vertI];
 358     forAll(myEdges, i)
 359     {
 360         const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
 361
 362         forAll(myFaces, myFaceI)
 363         {
 364             label faceI = myFaces[myFaceI];
 365
 366             if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
 367             {
 368                 faceStatus[faceI] = NOEDGE;
 369             }
 370        }
 371     }
 372 }
 373
 374
 375 //
 376 // Edge collapse helper functions
 377 //
 378
 379 // Get all faces that will get collapsed if edgeI collapses.
 380 Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
 381 (
 382     const triSurface& surf,
 383     label edgeI
 384 )
 385 {
 386     const edge& e = surf.edges()[edgeI];
 387     label v1 = e.start();
 388     label v2 = e.end();
 389
 390     // Faces using edge will certainly get collapsed.
 391     const labelList& myFaces = surf.edgeFaces()[edgeI];
 392
 393     labelHashSet facesToBeCollapsed(2*myFaces.size());
 394
 395     forAll(myFaces, myFaceI)
 396     {
 397         facesToBeCollapsed.insert(myFaces[myFaceI]);
 398     }
 399
 400     // From faces using v1 check if they share an edge with faces
 401     // using v2.
 402     //  - share edge: are part of 'splay' tree and will collapse if edge
 403     //    collapses
 404     const labelList& v1Faces = surf.pointFaces()[v1];
 405
 406     forAll(v1Faces, v1FaceI)
 407     {
 408         label face1I = v1Faces[v1FaceI];
 409
 410         label otherEdgeI = oppositeEdge(surf, face1I, v1);
 411
 412         // Step across edge to other face
 413         label face2I = otherFace(surf, face1I, otherEdgeI);
 414
 415         if (face2I != -1)
 416         {
 417             // found face on other side of edge. Now check if uses v2.
 418             if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
 419             {
 420                 // triangles face1I and face2I form splay tree and will
 421                 // collapse.
 422                 facesToBeCollapsed.insert(face1I);
 423                 facesToBeCollapsed.insert(face2I);
 424             }
 425         }
 426     }
 427
 428     return facesToBeCollapsed;
 429 }
 430
 431
 432 // Return value of faceUsed for faces using vertI
 433 Foam::label Foam::triSurfaceTools::vertexUsesFace
 434 (
 435     const triSurface& surf,
 436     const labelHashSet& faceUsed,
 437     const label vertI
 438 )
 439 {
 440     const labelList& myFaces = surf.pointFaces()[vertI];
 441
 442     forAll(myFaces, myFaceI)
 443     {
 444         label face1I = myFaces[myFaceI];
 445
 446         if (faceUsed.found(face1I))
 447         {
 448             return face1I;
 449         }
 450     }
 451     return -1;
 452 }
 453
 454
 455 // Calculate new edge-face addressing resulting from edge collapse
 456 void Foam::triSurfaceTools::getMergedEdges
 457 (
 458     const triSurface& surf,
 459     const label edgeI,
 460     const labelHashSet& collapsedFaces,
 461     HashTable<label, label, Hash<label> >& edgeToEdge,
 462     HashTable<label, label, Hash<label> >& edgeToFace
 463 )
 464 {
 465     const edge& e = surf.edges()[edgeI];
 466     label v1 = e.start();
 467     label v2 = e.end();
 468
 469     const labelList& v1Faces = surf.pointFaces()[v1];
 470     const labelList& v2Faces = surf.pointFaces()[v2];
 471
 472     // Mark all (non collapsed) faces using v2
 473     labelHashSet v2FacesHash(v2Faces.size());
 474
 475     forAll(v2Faces, v2FaceI)
 476     {
 477         if (!collapsedFaces.found(v2Faces[v2FaceI]))
 478         {
 479             v2FacesHash.insert(v2Faces[v2FaceI]);
 480         }
 481     }
 482
 483
 484     forAll(v1Faces, v1FaceI)
 485     {
 486         label face1I = v1Faces[v1FaceI];
 487
 488         if (collapsedFaces.found(face1I))
 489         {
 490             continue;
 491         }
 492
 493         //
 494         // Check if face1I uses a vertex connected to a face using v2
 495         //
 496
 497         label vert1I = -1;
 498         label vert2I = -1;
 499         otherVertices
 500         (
 501             surf,
 502             face1I,
 503             v1,
 504             vert1I,
 505             vert2I
 506         );
 507         //Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:"
 508         //    << vert1I << ' ' << vert2I << endl;
 509
 510         // Check vert1, vert2 for usage by v2Face.
 511
 512         label commonVert = vert1I;
 513         label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
 514         if (face2I == -1)
 515         {
 516             commonVert = vert2I;
 517             face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
 518         }
 519
 520         if (face2I != -1)
 521         {
 522             // Found one: commonVert is used by both face1I and face2I
 523             label edge1I = getEdge(surf, v1, commonVert);
 524             label edge2I = getEdge(surf, v2, commonVert);
 525
 526             edgeToEdge.insert(edge1I, edge2I);
 527             edgeToEdge.insert(edge2I, edge1I);
 528
 529             edgeToFace.insert(edge1I, face2I);
 530             edgeToFace.insert(edge2I, face1I);
 531         }
 532     }
 533 }
 534
 535
 536 // Calculates (cos of) angle across edgeI of faceI,
 537 // taking into account updated addressing (resulting from edge collapse)
 538 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
 539 (
 540     const triSurface& surf,
 541     const label v1,
 542     const point& pt,
 543     const labelHashSet& collapsedFaces,
 544     const HashTable<label, label, Hash<label> >& edgeToEdge,
 545     const HashTable<label, label, Hash<label> >& edgeToFace,
 546     const label faceI,
 547     const label edgeI
 548 )
 549 {
 550     const pointField& localPoints = surf.localPoints();
 551
 552     label A = surf.edges()[edgeI].start();
 553     label B = surf.edges()[edgeI].end();
 554     label C = oppositeVertex(surf, faceI, edgeI);
 555
 556     label D = -1;
 557
 558     label face2I = -1;
 559
 560     if (edgeToEdge.found(edgeI))
 561     {
 562         // Use merged addressing
 563         label edge2I = edgeToEdge[edgeI];
 564         face2I = edgeToFace[edgeI];
 565
 566         D = oppositeVertex(surf, face2I, edge2I);
 567     }
 568     else
 569     {
 570         // Use normal edge-face addressing
 571         face2I = otherFace(surf, faceI, edgeI);
 572
 573         if ((face2I != -1) && !collapsedFaces.found(face2I))
 574         {
 575             D = oppositeVertex(surf, face2I, edgeI);
 576         }
 577     }
 578
 579     scalar cosAngle = 1;
 580
 581     if (D != -1)
 582     {
 583         if (A == v1)
 584         {
 585             cosAngle = faceCosAngle
 586             (
 587                 pt,
 588                 localPoints[B],
 589                 localPoints[C],
 590                 localPoints[D]
 591             );
 592         }
 593         else if (B == v1)
 594         {
 595             cosAngle = faceCosAngle
 596             (
 597                 localPoints[A],
 598                 pt,
 599                 localPoints[C],
 600                 localPoints[D]
 601             );
 602         }
 603         else if (C == v1)
 604         {
 605             cosAngle = faceCosAngle
 606             (
 607                 localPoints[A],
 608                 localPoints[B],
 609                 pt,
 610                 localPoints[D]
 611             );
 612         }
 613         else if (D == v1)
 614         {
 615             cosAngle = faceCosAngle
 616             (
 617                 localPoints[A],
 618                 localPoints[B],
 619                 localPoints[C],
 620                 pt
 621             );
 622         }
 623         else
 624         {
 625             FatalErrorIn("edgeCosAngle")
 626                 << "face " << faceI << " does not use vertex "
 627                 << v1 << " of collapsed edge" << abort(FatalError);
 628         }
 629     }
 630     return cosAngle;
 631 }
 632
 633
 634 //- Calculate minimum (cos of) edge angle using addressing from collapsing
 635 //  edge to v1 at pt.
 636 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
 637 (
 638     const triSurface& surf,
 639     const label v1,
 640     const point& pt,
 641     const labelHashSet& collapsedFaces,
 642     const HashTable<label, label, Hash<label> >& edgeToEdge,
 643     const HashTable<label, label, Hash<label> >& edgeToFace
 644 )
 645 {
 646     const labelList& v1Faces = surf.pointFaces()[v1];
 647
 648     scalar minCos = 1;
 649
 650     forAll(v1Faces, v1FaceI)
 651     {
 652         label faceI = v1Faces[v1FaceI];
 653
 654         if (collapsedFaces.found(faceI))
 655         {
 656             continue;
 657         }
 658
 659         const labelList& myEdges = surf.faceEdges()[faceI];
 660
 661         forAll(myEdges, myEdgeI)
 662         {
 663             label edgeI = myEdges[myEdgeI];
 664
 665             minCos =
 666                 min
 667                 (
 668                     minCos,
 669                     edgeCosAngle
 670                     (
 671                         surf,
 672                         v1,
 673                         pt,
 674                         collapsedFaces,
 675                         edgeToEdge,
 676                         edgeToFace,
 677                         faceI,
 678                         edgeI
 679                     )
 680                 );
 681         }
 682     }
 683
 684     return minCos;
 685 }
 686
 687
 688 // Calculate max dihedral angle after collapsing edge to v1 (at pt).
 689 // Note that all edges of all faces using v1 are affected.
 690 bool Foam::triSurfaceTools::collapseCreatesFold
 691 (
 692     const triSurface& surf,
 693     const label v1,
 694     const point& pt,
 695     const labelHashSet& collapsedFaces,
 696     const HashTable<label, label, Hash<label> >& edgeToEdge,
 697     const HashTable<label, label, Hash<label> >& edgeToFace,
 698     const scalar minCos
 699 )
 700 {
 701     const labelList& v1Faces = surf.pointFaces()[v1];
 702
 703     forAll(v1Faces, v1FaceI)
 704     {
 705         label faceI = v1Faces[v1FaceI];
 706
 707         if (collapsedFaces.found(faceI))
 708         {
 709             continue;
 710         }
 711
 712         const labelList& myEdges = surf.faceEdges()[faceI];
 713
 714         forAll(myEdges, myEdgeI)
 715         {
 716             label edgeI = myEdges[myEdgeI];
 717
 718             if
 719             (
 720                 edgeCosAngle
 721                 (
 722                     surf,
 723                     v1,
 724                     pt,
 725                     collapsedFaces,
 726                     edgeToEdge,
 727                     edgeToFace,
 728                     faceI,
 729                     edgeI
 730                 )
 731               < minCos
 732             )
 733             {
 734                 return true;
 735             }
 736         }
 737     }
 738
 739     return false;
 740 }
 741
 742
 743 //// Return true if collapsing edgeI creates triangles on top of each other.
 744 //// Is when the triangles neighbouring the collapsed one already share
 745 // a vertex.
 746 //bool Foam::triSurfaceTools::collapseCreatesDuplicates
 747 //(
 748 //    const triSurface& surf,
 749 //    const label edgeI,
 750 //    const labelHashSet& collapsedFaces
 751 //)
 752 //{
 753 //Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI]
 754 //    << " collapsedFaces:" << collapsedFaces.toc() << endl;
 755 //
 756 //    // Mark neighbours of faces to be collapsed, i.e. get the first layer
 757 //    // of triangles outside collapsedFaces.
 758 //    // neighbours actually contains the
 759 //    // edge with which triangle connects to collapsedFaces.
 760 //
 761 //    HashTable<label, label, Hash<label> > neighbours;
 762 //
 763 //    labelList collapsed = collapsedFaces.toc();
 764 //
 765 //    forAll(collapsed, collapseI)
 766 //    {
 767 //        const label faceI = collapsed[collapseI];
 768 //
 769 //        const labelList& myEdges = surf.faceEdges()[faceI];
 770 //
 771 //        Pout<< "collapsing faceI:" << faceI << " uses edges:" << myEdges
 772 //            << endl;
 773 //
 774 //        forAll(myEdges, myEdgeI)
 775 //        {
 776 //            const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]];
 777 //
 778 //            Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces "
 779 //                << myFaces << endl;
 780 //
 781 //            if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2))
 782 //            {
 783 //                // Get the other face
 784 //                label neighbourFaceI = myFaces[0];
 785 //
 786 //                if (neighbourFaceI == faceI)
 787 //                {
 788 //                    neighbourFaceI = myFaces[1];
 789 //                }
 790 //
 791 //                // Is 'outside' face if not in collapsedFaces itself
 792 //                if (!collapsedFaces.found(neighbourFaceI))
 793 //                {
 794 //                    neighbours.insert(neighbourFaceI, myEdges[myEdgeI]);
 795 //                }
 796 //            }
 797 //        }
 798 //    }
 799 //
 800 //    // Now neighbours contains first layer of triangles outside of
 801 //    // collapseFaces
 802 //    // There should be
 803 //    // -two if edgeI is a boundary edge
 804 //    // since the outside 'edge' of collapseFaces should
 805 //    // form a triangle and the face connected to edgeI is not inserted.
 806 //    // -four if edgeI is not a boundary edge since then the outside edge forms
 807 //    // a diamond.
 808 //
 809 //    // Check if any of the faces in neighbours share an edge. (n^2)
 810 //
 811 //    labelList neighbourList = neighbours.toc();
 812 //
 813 //Pout<< "edgeI:" << edgeI << "  neighbourList:" << neighbourList << endl;
 814 //
 815 //
 816 //    forAll(neighbourList, i)
 817 //    {
 818 //        const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]];
 819 //
 820 //        for (label j = i+1; j < neighbourList.size(); i++)
 821 //        {
 822 //            const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
 823 //
 824 //            // Check if faceI and faceJ share an edge
 825 //            forAll(faceIEdges, fI)
 826 //            {
 827 //                forAll(faceJEdges, fJ)
 828 //                {
 829 //                    Pout<< " comparing " << faceIEdges[fI] << " to "
 830 //                        << faceJEdges[fJ] << endl;
 831 //
 832 //                    if (faceIEdges[fI] == faceJEdges[fJ])
 833 //                    {
 834 //                        return true;
 835 //                    }
 836 //                }
 837 //            }
 838 //        }
 839 //    }
 840 //    Pout<< "Found no match. Returning false" << endl;
 841 //    return false;
 842 //}
 843
 844
 845 // Finds the triangle edge cut by the plane between a point inside / on edge
 846 // of a triangle and a point outside. Returns:
 847 //  - cut through edge/point
 848 //  - miss()
 849 Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
 850 (
 851     const triSurface& s,
 852     const label triI,
 853     const label excludeEdgeI,
 854     const label excludePointI,
 855
 856     const point& triPoint,
 857     const plane& cutPlane,
 858     const point& toPoint
 859 )
 860 {
 861     const pointField& points = s.points();
 862     const labelledTri& f = s[triI];
 863     const labelList& fEdges = s.faceEdges()[triI];
 864
 865     // Get normal distance to planeN
 866     FixedList<scalar, 3> d;
 867
 868     scalar norm = 0;
 869     forAll(d, fp)
 870     {
 871         d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
 872         norm += mag(d[fp]);
 873     }
 874
 875     // Normalise and truncate
 876     forAll(d, i)
 877     {
 878         d[i] /= norm;
 879
 880         if (mag(d[i]) < 1E-6)
 881         {
 882             d[i] = 0.0;
 883         }
 884     }
 885
 886     // Return information
 887     surfaceLocation cut;
 888
 889     if (excludePointI != -1)
 890     {
 891         // Excluded point. Test only opposite edge.
 892
 893         label fp0 = findIndex(s.localFaces()[triI], excludePointI);
 894
 895         if (fp0 == -1)
 896         {
 897             FatalErrorIn("cutEdge(..)") << "excludePointI:" << excludePointI
 898                 << " localF:" << s.localFaces()[triI] << abort(FatalError);
 899         }
 900
 901         label fp1 = f.fcIndex(fp0);
 902         label fp2 = f.fcIndex(fp1);
 903
 904
 905         if (d[fp1] == 0.0)
 906         {
 907             cut.setHit();
 908             cut.setPoint(points[f[fp1]]);
 909             cut.elementType() = triPointRef::POINT;
 910             cut.setIndex(s.localFaces()[triI][fp1]);
 911         }
 912         else if (d[fp2] == 0.0)
 913         {
 914             cut.setHit();
 915             cut.setPoint(points[f[fp2]]);
 916             cut.elementType() = triPointRef::POINT;
 917             cut.setIndex(s.localFaces()[triI][fp2]);
 918         }
 919         else if
 920         (
 921             (d[fp1] < 0 && d[fp2] < 0)
 922          || (d[fp1] > 0 && d[fp2] > 0)
 923         )
 924         {
 925             // Both same sign. Not crossing edge at all.
 926             // cut already set to miss().
 927         }
 928         else
 929         {
 930             cut.setHit();
 931             cut.setPoint
 932             (
 933                 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
 934               / (d[fp2] - d[fp1])
 935             );
 936             cut.elementType() = triPointRef::EDGE;
 937             cut.setIndex(fEdges[fp1]);
 938         }
 939     }
 940     else
 941     {
 942         // Find the two intersections
 943         FixedList<surfaceLocation, 2> inters;
 944         label interI = 0;
 945
 946         forAll(f, fp0)
 947         {
 948             label fp1 = f.fcIndex(fp0);
 949
 950             if (d[fp0] == 0)
 951             {
 952                 if (interI >= 2)
 953                 {
 954                     FatalErrorIn("cutEdge(..)")
 955                         << "problem : triangle has three intersections." << nl
 956                         << "triangle:" << f.tri(points)
 957                         << " d:" << d << abort(FatalError);
 958                 }
 959                 inters[interI].setHit();
 960                 inters[interI].setPoint(points[f[fp0]]);
 961                 inters[interI].elementType() = triPointRef::POINT;
 962                 inters[interI].setIndex(s.localFaces()[triI][fp0]);
 963                 interI++;
 964             }
 965             else if
 966             (
 967                 (d[fp0] < 0 && d[fp1] > 0)
 968              || (d[fp0] > 0 && d[fp1] < 0)
 969             )
 970             {
 971                 if (interI >= 2)
 972                 {
 973                     FatalErrorIn("cutEdge(..)")
 974                         << "problem : triangle has three intersections." << nl
 975                         << "triangle:" << f.tri(points)
 976                         << " d:" << d << abort(FatalError);
 977                 }
 978                 inters[interI].setHit();
 979                 inters[interI].setPoint
 980                 (
 981                     (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
 982                   / (d[fp0] - d[fp1])
 983                 );
 984                 inters[interI].elementType() = triPointRef::EDGE;
 985                 inters[interI].setIndex(fEdges[fp0]);
 986                 interI++;
 987             }
 988         }
 989
 990
 991         if (interI == 0)
 992         {
 993             // Return miss
 994         }
 995         else if (interI == 1)
 996         {
 997             // Only one intersection. Should not happen!
 998             cut = inters[0];
 999         }
1000         else if (interI == 2)
1001         {
1002             // Handle excludeEdgeI
1003             if
1004             (
1005                 inters[0].elementType() == triPointRef::EDGE
1006              && inters[0].index() == excludeEdgeI
1007             )
1008             {
1009                 cut = inters[1];
1010             }
1011             else if
1012             (
1013                 inters[1].elementType() == triPointRef::EDGE
1014              && inters[1].index() == excludeEdgeI
1015             )
1016             {
1017                 cut = inters[0];
1018             }
1019             else
1020             {
1021                 // Two cuts. Find nearest.
1022                 if
1023                 (
1024                     magSqr(inters[0].rawPoint() - toPoint)
1025                   < magSqr(inters[1].rawPoint() - toPoint)
1026                 )
1027                 {
1028                     cut = inters[0];
1029                 }
1030                 else
1031                 {
1032                     cut = inters[1];
1033                 }
1034             }
1035         }
1036     }
1037     return cut;
1038 }
1039
1040
1041 // 'Snap' point on to endPoint.
1042 void Foam::triSurfaceTools::snapToEnd
1043 (
1044     const triSurface& s,
1045     const surfaceLocation& end,
1046     surfaceLocation& current
1047 )
1048 {
1049     if (end.elementType() == triPointRef::NONE)
1050     {
1051         if (current.elementType() == triPointRef::NONE)
1052         {
1053             // endpoint on triangle; current on triangle
1054             if (current.index() == end.index())
1055             {
1056                 //if (debug)
1057                 //{
1058                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1059                 //        << end << endl;
1060                 //}
1061                 current = end;
1062                 current.setHit();
1063             }
1064         }
1065         // No need to handle current on edge/point since tracking handles this.
1066     }
1067     else if (end.elementType() == triPointRef::EDGE)
1068     {
1069         if (current.elementType() == triPointRef::NONE)
1070         {
1071             // endpoint on edge; current on triangle
1072             const labelList& fEdges = s.faceEdges()[current.index()];
1073
1074             if (findIndex(fEdges, end.index()) != -1)
1075             {
1076                 //if (debug)
1077                 //{
1078                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1079                 //        << end << endl;
1080                 //}
1081                 current = end;
1082                 current.setHit();
1083             }
1084         }
1085         else if (current.elementType() == triPointRef::EDGE)
1086         {
1087             // endpoint on edge; current on edge
1088             if (current.index() == end.index())
1089             {
1090                 //if (debug)
1091                 //{
1092                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1093                 //        << end << endl;
1094                 //}
1095                 current = end;
1096                 current.setHit();
1097             }
1098         }
1099         else
1100         {
1101             // endpoint on edge; current on point
1102             const edge& e = s.edges()[end.index()];
1103
1104             if (current.index() == e[0] || current.index() == e[1])
1105             {
1106                 //if (debug)
1107                 //{
1108                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1109                 //        << end << endl;
1110                 //}
1111                 current = end;
1112                 current.setHit();
1113             }
1114         }
1115     }
1116     else    // end.elementType() == POINT
1117     {
1118         if (current.elementType() == triPointRef::NONE)
1119         {
1120             // endpoint on point; current on triangle
1121             const labelledTri& f = s.localFaces()[current.index()];
1122
1123             if (findIndex(f, end.index()) != -1)
1124             {
1125                 //if (debug)
1126                 //{
1127                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1128                 //        << end << endl;
1129                 //}
1130                 current = end;
1131                 current.setHit();
1132             }
1133         }
1134         else if (current.elementType() == triPointRef::EDGE)
1135         {
1136             // endpoint on point; current on edge
1137             const edge& e = s.edges()[current.index()];
1138
1139             if (end.index() == e[0] || end.index() == e[1])
1140             {
1141                 //if (debug)
1142                 //{
1143                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1144                 //        << end << endl;
1145                 //}
1146                 current = end;
1147                 current.setHit();
1148             }
1149         }
1150         else
1151         {
1152             // endpoint on point; current on point
1153             if (current.index() == end.index())
1154             {
1155                 //if (debug)
1156                 //{
1157                 //    Pout<< "snapToEnd : snapping:" << current << " onto:"
1158                 //        << end << endl;
1159                 //}
1160                 current = end;
1161                 current.setHit();
1162             }
1163         }
1164     }
1165 }
1166
1167
1168 // Start:
1169 //  - location
1170 //  - element type (triangle/edge/point) and index
1171 //  - triangle to exclude
1172 Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
1173 (
1174     const triSurface& s,
1175     const labelList& eFaces,
1176     const surfaceLocation& start,
1177     const label excludeEdgeI,
1178     const label excludePointI,
1179     const surfaceLocation& end,
1180     const plane& cutPlane
1181 )
1182 {
1183     surfaceLocation nearest;
1184
1185     scalar minDistSqr = Foam::sqr(GREAT);
1186
1187     forAll(eFaces, i)
1188     {
1189         label triI = eFaces[i];
1190
1191         // Make sure we don't revisit previous face
1192         if (triI != start.triangle())
1193         {
1194             if (end.elementType() == triPointRef::NONE && end.index() == triI)
1195             {
1196                 // Endpoint is in this triangle. Jump there.
1197                 nearest = end;
1198                 nearest.setHit();
1199                 nearest.triangle() = triI;
1200                 break;
1201             }
1202             else
1203             {
1204                // Which edge is cut.
1205
1206                 surfaceLocation cutInfo = cutEdge
1207                 (
1208                     s,
1209                     triI,
1210                     excludeEdgeI,       // excludeEdgeI
1211                     excludePointI,      // excludePointI
1212                     start.rawPoint(),
1213                     cutPlane,
1214                     end.rawPoint()
1215                 );
1216
1217                 // If crossing an edge we expect next edge to be cut.
1218                 if (excludeEdgeI != -1 && !cutInfo.hit())
1219                 {
1220                     FatalErrorIn("triSurfaceTools::visitFaces(..)")
1221                         << "Triangle:" << triI
1222                         << " excludeEdge:" << excludeEdgeI
1223                         << " point:" << start.rawPoint()
1224                         << " plane:" << cutPlane
1225                         << " . No intersection!" << abort(FatalError);
1226                 }
1227
1228                 if (cutInfo.hit())
1229                 {
1230                     scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
1231
1232                     if (distSqr < minDistSqr)
1233                     {
1234                         minDistSqr = distSqr;
1235                         nearest = cutInfo;
1236                         nearest.triangle() = triI;
1237                         nearest.setMiss();
1238                     }
1239                 }
1240             }
1241         }
1242     }
1243
1244     if (nearest.triangle() == -1)
1245     {
1246         // Did not move from edge. Give warning? Return something special?
1247         // For now responsability of caller to make sure that nothing has
1248         // moved.
1249     }
1250
1251     return nearest;
1252 }
1253
1254
1255 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1256
1257
1258 // Write pointField to file
1259 void Foam::triSurfaceTools::writeOBJ
1260 (
1261     const fileName& fName,
1262     const pointField& pts
1263 )
1264 {
1265     OFstream outFile(fName);
1266
1267     forAll(pts, pointI)
1268     {
1269         const point& pt = pts[pointI];
1270
1271         outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1272     }
1273     Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
1274 }
1275
1276
1277 // Write vertex subset to OBJ format file
1278 void Foam::triSurfaceTools::writeOBJ
1279 (
1280     const triSurface& surf,
1281     const fileName& fName,
1282     const boolList& markedVerts
1283 )
1284 {
1285     OFstream outFile(fName);
1286
1287     label nVerts = 0;
1288     forAll(markedVerts, vertI)
1289     {
1290         if (markedVerts[vertI])
1291         {
1292             const point& pt = surf.localPoints()[vertI];
1293
1294             outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1295
1296             nVerts++;
1297         }
1298     }
1299     Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
1300 }
1301
1302
1303 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1304 // Addressing helper functions:
1305
1306
1307 // Get all triangles using vertices of edge
1308 void Foam::triSurfaceTools::getVertexTriangles
1309 (
1310     const triSurface& surf,
1311     const label edgeI,
1312     labelList& edgeTris
1313 )
1314 {
1315     const edge& e = surf.edges()[edgeI];
1316     const labelList& myFaces = surf.edgeFaces()[edgeI];
1317
1318     label face1I = myFaces[0];
1319     label face2I = -1;
1320     if (myFaces.size() == 2)
1321     {
1322         face2I = myFaces[1];
1323     }
1324
1325     const labelList& startFaces = surf.pointFaces()[e.start()];
1326     const labelList& endFaces = surf.pointFaces()[e.end()];
1327
1328     // Number of triangles is sum of pointfaces - common faces
1329     // (= faces using edge)
1330     edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
1331
1332     label nTris = 0;
1333     forAll(startFaces, startFaceI)
1334     {
1335         edgeTris[nTris++] = startFaces[startFaceI];
1336     }
1337
1338     forAll(endFaces, endFaceI)
1339     {
1340         label faceI = endFaces[endFaceI];
1341
1342         if ((faceI != face1I) && (faceI != face2I))
1343         {
1344             edgeTris[nTris++] = faceI;
1345         }
1346     }
1347 }
1348
1349
1350 // Get all vertices connected to vertices of edge
1351 Foam::labelList Foam::triSurfaceTools::getVertexVertices
1352 (
1353     const triSurface& surf,
1354     const edge& e
1355 )
1356 {
1357     const edgeList& edges = surf.edges();
1358
1359     label v1 = e.start();
1360     label v2 = e.end();
1361
1362     // Get all vertices connected to v1 or v2 through an edge
1363     labelHashSet vertexNeighbours;
1364
1365     const labelList& v1Edges = surf.pointEdges()[v1];
1366
1367     forAll(v1Edges, v1EdgeI)
1368     {
1369         const edge& e = edges[v1Edges[v1EdgeI]];
1370         vertexNeighbours.insert(e.otherVertex(v1));
1371     }
1372
1373     const labelList& v2Edges = surf.pointEdges()[v2];
1374
1375     forAll(v2Edges, v2EdgeI)
1376     {
1377         const edge& e = edges[v2Edges[v2EdgeI]];
1378
1379         label vertI = e.otherVertex(v2);
1380
1381         vertexNeighbours.insert(vertI);
1382     }
1383     return vertexNeighbours.toc();
1384 }
1385
1386
1387 //// Order vertices consistent with face
1388 //void Foam::triSurfaceTools::orderVertices
1389 //(
1390 //    const labelledTri& f,
1391 //    const label v1,
1392 //    const label v2,
1393 //    label& vA,
1394 //    label& vB
1395 //)
1396 //{
1397 //    // Order v1, v2 in anticlockwise order.
1398 //    bool reverse = false;
1399 //
1400 //    if (f[0] == v1)
1401 //    {
1402 //        if (f[1] != v2)
1403 //        {
1404 //            reverse = true;
1405 //        }
1406 //    }
1407 //    else if (f[1] == v1)
1408 //    {
1409 //        if (f[2] != v2)
1410 //        {
1411 //            reverse = true;
1412 //        }
1413 //    }
1414 //    else
1415 //    {
1416 //        if (f[0] != v2)
1417 //        {
1418 //            reverse = true;
1419 //        }
1420 //    }
1421 //
1422 //    if (reverse)
1423 //    {
1424 //        vA = v2;
1425 //        vB = v1;
1426 //    }
1427 //    else
1428 //    {
1429 //        vA = v1;
1430 //        vB = v2;
1431 //    }
1432 //}
1433
1434
1435 // Get the other face using edgeI
1436 Foam::label Foam::triSurfaceTools::otherFace
1437 (
1438     const triSurface& surf,
1439     const label faceI,
1440     const label edgeI
1441 )
1442 {
1443     const labelList& myFaces = surf.edgeFaces()[edgeI];
1444
1445     if (myFaces.size() != 2)
1446     {
1447         return -1;
1448     }
1449     else
1450     {
1451         if (faceI == myFaces[0])
1452         {
1453             return myFaces[1];
1454         }
1455         else
1456         {
1457             return myFaces[0];
1458         }
1459     }
1460 }
1461
1462
1463 // Get the two edges on faceI counterclockwise after edgeI
1464 void Foam::triSurfaceTools::otherEdges
1465 (
1466     const triSurface& surf,
1467     const label faceI,
1468     const label edgeI,
1469     label& e1,
1470     label& e2
1471 )
1472 {
1473     const labelList& eFaces = surf.faceEdges()[faceI];
1474
1475     label i0 = findIndex(eFaces, edgeI);
1476
1477     if (i0 == -1)
1478     {
1479         FatalErrorIn
1480         (
1481             "otherEdges"
1482             "(const triSurface&, const label, const label,"
1483             " label&, label&)"
1484         )   << "Edge " << surf.edges()[edgeI] << " not in face "
1485             << surf.localFaces()[faceI] << abort(FatalError);
1486     }
1487
1488     label i1 = eFaces.fcIndex(i0);
1489     label i2 = eFaces.fcIndex(i1);
1490
1491     e1 = eFaces[i1];
1492     e2 = eFaces[i2];
1493 }
1494
1495
1496 // Get the two vertices on faceI counterclockwise vertI
1497 void Foam::triSurfaceTools::otherVertices
1498 (
1499     const triSurface& surf,
1500     const label faceI,
1501     const label vertI,
1502     label& vert1I,
1503     label& vert2I
1504 )
1505 {
1506     const labelledTri& f = surf.localFaces()[faceI];
1507
1508     if (vertI == f[0])
1509     {
1510         vert1I = f[1];
1511         vert2I = f[2];
1512     }
1513     else if (vertI == f[1])
1514     {
1515         vert1I = f[2];
1516         vert2I = f[0];
1517     }
1518     else if (vertI == f[2])
1519     {
1520         vert1I = f[0];
1521         vert2I = f[1];
1522     }
1523     else
1524     {
1525         FatalErrorIn
1526         (
1527             "otherVertices"
1528             "(const triSurface&, const label, const label,"
1529             " label&, label&)"
1530         )   << "Vertex " << vertI << " not in face " << f << abort(FatalError);
1531     }
1532 }
1533
1534
1535 // Get edge opposite vertex
1536 Foam::label Foam::triSurfaceTools::oppositeEdge
1537 (
1538     const triSurface& surf,
1539     const label faceI,
1540     const label vertI
1541 )
1542 {
1543     const labelList& myEdges = surf.faceEdges()[faceI];
1544
1545     forAll(myEdges, myEdgeI)
1546     {
1547         label edgeI = myEdges[myEdgeI];
1548
1549         const edge& e = surf.edges()[edgeI];
1550
1551         if ((e.start() != vertI) && (e.end() != vertI))
1552         {
1553             return edgeI;
1554         }
1555     }
1556
1557     FatalErrorIn
1558     (
1559         "oppositeEdge"
1560         "(const triSurface&, const label, const label)"
1561     )   << "Cannot find vertex " << vertI << " in edges of face " << faceI
1562         << abort(FatalError);
1563
1564     return -1;
1565 }
1566
1567
1568 // Get vertex opposite edge
1569 Foam::label Foam::triSurfaceTools::oppositeVertex
1570 (
1571     const triSurface& surf,
1572     const label faceI,
1573     const label edgeI
1574 )
1575 {
1576     const labelledTri& f = surf.localFaces()[faceI];
1577
1578     const edge& e = surf.edges()[edgeI];
1579
1580     forAll(f, fp)
1581     {
1582         label vertI = f[fp];
1583
1584         if (vertI != e.start() && vertI != e.end())
1585         {
1586             return vertI;
1587         }
1588     }
1589
1590     FatalErrorIn("triSurfaceTools::oppositeVertex")
1591         << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
1592         << " in face " << faceI << " vertices " << f << abort(FatalError);
1593
1594     return -1;
1595 }
1596
1597
1598 // Returns edge label connecting v1, v2
1599 Foam::label Foam::triSurfaceTools::getEdge
1600 (
1601     const triSurface& surf,
1602     const label v1,
1603     const label v2
1604 )
1605 {
1606     const labelList& v1Edges = surf.pointEdges()[v1];
1607
1608     forAll(v1Edges, v1EdgeI)
1609     {
1610         label edgeI = v1Edges[v1EdgeI];
1611
1612         const edge& e = surf.edges()[edgeI];
1613
1614         if ((e.start() == v2) || (e.end() == v2))
1615         {
1616             return edgeI;
1617         }
1618     }
1619     return -1;
1620 }
1621
1622
1623 // Return index of triangle (or -1) using all three edges
1624 Foam::label Foam::triSurfaceTools::getTriangle
1625 (
1626     const triSurface& surf,
1627     const label e0I,
1628     const label e1I,
1629     const label e2I
1630 )
1631 {
1632     if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1633     {
1634         FatalErrorIn
1635         (
1636             "getTriangle"
1637             "(const triSurface&, const label, const label,"
1638             " const label)"
1639         )   << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
1640             << " e2:" << e2I
1641             << abort(FatalError);
1642     }
1643
1644     const labelList& eFaces = surf.edgeFaces()[e0I];
1645
1646     forAll(eFaces, eFaceI)
1647     {
1648         label faceI = eFaces[eFaceI];
1649
1650         const labelList& myEdges = surf.faceEdges()[faceI];
1651
1652         if
1653         (
1654             (myEdges[0] == e1I)
1655          || (myEdges[1] == e1I)
1656          || (myEdges[2] == e1I)
1657         )
1658         {
1659             if
1660             (
1661                 (myEdges[0] == e2I)
1662              || (myEdges[1] == e2I)
1663              || (myEdges[2] == e2I)
1664             )
1665             {
1666                 return faceI;
1667             }
1668         }
1669     }
1670     return -1;
1671 }
1672
1673
1674 // Collapse indicated edges. Return new tri surface.
1675 Foam::triSurface Foam::triSurfaceTools::collapseEdges
1676 (
1677     const triSurface& surf,
1678     const labelList& collapsableEdges
1679 )
1680 {
1681     pointField edgeMids(surf.nEdges());
1682
1683     forAll(edgeMids, edgeI)
1684     {
1685         const edge& e = surf.edges()[edgeI];
1686
1687         edgeMids[edgeI] =
1688             0.5
1689           * (
1690                 surf.localPoints()[e.start()]
1691               + surf.localPoints()[e.end()]
1692             );
1693     }
1694
1695
1696     labelList faceStatus(surf.size(), ANYEDGE);
1697
1698     //// Protect triangles which are on the border of different regions
1699     //forAll(edges, edgeI)
1700     //{
1701     //    const labelList& neighbours = edgeFaces[edgeI];
1702     //
1703     //    if ((neighbours.size() != 2) && (neighbours.size() != 1))
1704     //    {
1705     //        FatalErrorIn("collapseEdges")
1706     //            << abort(FatalError);
1707     //    }
1708     //
1709     //    if (neighbours.size() == 2)
1710     //    {
1711     //        if (surf[neighbours[0]].region() != surf[neighbours[1]].region())
1712     //        {
1713     //            // Neighbours on different regions. For now dont allow
1714     //            // any collapse.
1715     //            //Pout<< "protecting face " << neighbours[0]
1716     //            //    << ' ' << neighbours[1] << endl;
1717     //            faceStatus[neighbours[0]] = NOEDGE;
1718     //            faceStatus[neighbours[1]] = NOEDGE;
1719     //        }
1720     //    }
1721     //}
1722
1723     return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1724 }
1725
1726
1727 // Collapse indicated edges. Return new tri surface.
1728 Foam::triSurface Foam::triSurfaceTools::collapseEdges
1729 (
1730     const triSurface& surf,
1731     const labelList& collapseEdgeLabels,
1732     const pointField& edgeMids,
1733     labelList& faceStatus
1734 )
1735 {
1736     const labelListList& edgeFaces = surf.edgeFaces();
1737     const pointField& localPoints = surf.localPoints();
1738     const edgeList& edges = surf.edges();
1739
1740     // Storage for new points
1741     pointField newPoints(localPoints);
1742
1743     // Map for old to new points
1744     labelList pointMap(localPoints.size());
1745     forAll(localPoints, pointI)
1746     {
1747         pointMap[pointI] = pointI;
1748     }
1749
1750
1751     // Do actual 'collapsing' of edges
1752
1753     forAll(collapseEdgeLabels, collapseEdgeI)
1754     {
1755         const label edgeI = collapseEdgeLabels[collapseEdgeI];
1756
1757         if ((edgeI < 0) || (edgeI >= surf.nEdges()))
1758         {
1759             FatalErrorIn("collapseEdges")
1760                 << "Edge label outside valid range." << endl
1761                 << "edge label:" << edgeI << endl
1762                 << "total number of edges:" << surf.nEdges() << endl
1763                 << abort(FatalError);
1764         }
1765
1766         const labelList& neighbours = edgeFaces[edgeI];
1767
1768         if (neighbours.size() == 2)
1769         {
1770             const label stat0 = faceStatus[neighbours[0]];
1771             const label stat1 = faceStatus[neighbours[1]];
1772
1773             // Check faceStatus to make sure this one can be collapsed
1774             if
1775             (
1776                 ((stat0 == ANYEDGE) || (stat0 == edgeI))
1777              && ((stat1 == ANYEDGE) || (stat1 == edgeI))
1778             )
1779             {
1780                 const edge& e = edges[edgeI];
1781
1782                 // Set up mapping to 'collapse' points of edge
1783                 if
1784                 (
1785                     (pointMap[e.start()] != e.start())
1786                  || (pointMap[e.end()] != e.end())
1787                 )
1788                 {
1789                     FatalErrorIn("collapseEdges")
1790                         << "points already mapped. Double collapse." << endl
1791                         << "edgeI:" << edgeI
1792                         << "  start:" << e.start()
1793                         << "  end:" << e.end()
1794                         << "  pointMap[start]:" << pointMap[e.start()]
1795                         << "  pointMap[end]:" << pointMap[e.end()]
1796                         << abort(FatalError);
1797                 }
1798
1799                 const label minVert = min(e.start(), e.end());
1800                 pointMap[e.start()] = minVert;
1801                 pointMap[e.end()] = minVert;
1802
1803                 // Move shared vertex to mid of edge
1804                 newPoints[minVert] = edgeMids[edgeI];
1805
1806                 // Protect neighbouring faces
1807                 protectNeighbours(surf, e.start(), faceStatus);
1808                 protectNeighbours(surf, e.end(), faceStatus);
1809                 protectNeighbours
1810                 (
1811                     surf,
1812                     oppositeVertex(surf, neighbours[0], edgeI),
1813                     faceStatus
1814                 );
1815                 protectNeighbours
1816                 (
1817                     surf,
1818                     oppositeVertex(surf, neighbours[1], edgeI),
1819                     faceStatus
1820                 );
1821
1822                 // Mark all collapsing faces
1823                 labelList collapseFaces =
1824                     getCollapsedFaces
1825                     (
1826                         surf,
1827                         edgeI
1828                     ).toc();
1829
1830                 forAll(collapseFaces, collapseI)
1831                 {
1832                      faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1833                 }
1834             }
1835         }
1836     }
1837
1838
1839     // Storage for new triangles
1840     List<labelledTri> newTris(surf.size());
1841     label newTriI = 0;
1842
1843     const List<labelledTri>& localFaces = surf.localFaces();
1844
1845
1846     // Get only non-collapsed triangles and renumber vertex labels.
1847     forAll(localFaces, faceI)
1848     {
1849         const labelledTri& f = localFaces[faceI];
1850
1851         const label a = pointMap[f[0]];
1852         const label b = pointMap[f[1]];
1853         const label c = pointMap[f[2]];
1854
1855         if
1856         (
1857             (a != b) && (a != c) && (b != c)
1858          && (faceStatus[faceI] != COLLAPSED)
1859         )
1860         {
1861             // uncollapsed triangle
1862             newTris[newTriI++] = labelledTri(a, b, c, f.region());
1863         }
1864         else
1865         {
1866             //Pout<< "Collapsed triangle " << faceI
1867             //    << " vertices:" << f << endl;
1868         }
1869     }
1870     newTris.setSize(newTriI);
1871
1872
1873
1874     // Pack faces
1875
1876     triSurface tempSurf(newTris, surf.patches(), newPoints);
1877
1878     return
1879         triSurface
1880         (
1881             tempSurf.localFaces(),
1882             tempSurf.patches(),
1883             tempSurf.localPoints()
1884         );
1885 }
1886
1887
1888 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1889
1890 Foam::triSurface Foam::triSurfaceTools::redGreenRefine
1891 (
1892     const triSurface& surf,
1893     const labelList& refineFaces
1894 )
1895 {
1896     List<refineType> refineStatus(surf.size(), NONE);
1897
1898     // Mark & propagate refinement
1899     forAll(refineFaces, refineFaceI)
1900     {
1901         calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus);
1902     }
1903
1904     // Do actual refinement
1905     return doRefine(surf, refineStatus);
1906 }
1907
1908
1909 Foam::triSurface Foam::triSurfaceTools::greenRefine
1910 (
1911     const triSurface& surf,
1912     const labelList& refineEdges
1913 )
1914 {
1915     // Storage for marking faces
1916     List<refineType> refineStatus(surf.size(), NONE);
1917
1918     // Storage for new faces
1919     DynamicList<labelledTri> newFaces(0);
1920
1921     pointField newPoints(surf.localPoints());
1922     newPoints.setSize(surf.nPoints() + surf.nEdges());
1923     label newPointI = surf.nPoints();
1924
1925
1926     // Refine edges
1927     forAll(refineEdges, refineEdgeI)
1928     {
1929         label edgeI = refineEdges[refineEdgeI];
1930
1931         const labelList& myFaces = surf.edgeFaces()[edgeI];
1932
1933         bool neighbourIsRefined= false;
1934
1935         forAll(myFaces, myFaceI)
1936         {
1937             if (refineStatus[myFaces[myFaceI]] != NONE)
1938             {
1939                 neighbourIsRefined =  true;
1940             }
1941         }
1942
1943         // Only refine if none of the faces is refined
1944         if (!neighbourIsRefined)
1945         {
1946             // Refine edge
1947             const edge& e = surf.edges()[edgeI];
1948
1949             point mid =
1950                 0.5
1951               * (
1952                     surf.localPoints()[e.start()]
1953                   + surf.localPoints()[e.end()]
1954                 );
1955
1956             newPoints[newPointI] = mid;
1957
1958             // Refine faces using edge
1959             forAll(myFaces, myFaceI)
1960             {
1961                 // Add faces to newFaces
1962                 greenRefine
1963                 (
1964                     surf,
1965                     myFaces[myFaceI],
1966                     edgeI,
1967                     newPointI,
1968                     newFaces
1969                 );
1970
1971                 // Mark as refined
1972                 refineStatus[myFaces[myFaceI]] = GREEN;
1973             }
1974
1975             newPointI++;
1976         }
1977     }
1978
1979     // Add unrefined faces
1980     forAll(surf.localFaces(), faceI)
1981     {
1982         if (refineStatus[faceI] == NONE)
1983         {
1984             newFaces.append(surf.localFaces()[faceI]);
1985         }
1986     }
1987
1988     newFaces.shrink();
1989     newPoints.setSize(newPointI);
1990
1991     return triSurface(newFaces, surf.patches(), newPoints);
1992 }
1993
1994
1995
1996 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1997 // Geometric helper functions:
1998
1999
2000 // Returns element in edgeIndices with minimum length
2001 Foam::label Foam::triSurfaceTools::minEdge
2002 (
2003     const triSurface& surf,
2004     const labelList& edgeIndices
2005 )
2006 {
2007     scalar minLength = GREAT;
2008     label minIndex = -1;
2009     forAll(edgeIndices, i)
2010     {
2011         const edge& e = surf.edges()[edgeIndices[i]];
2012
2013         scalar length =
2014             mag
2015             (
2016                 surf.localPoints()[e.end()]
2017               - surf.localPoints()[e.start()]
2018             );
2019
2020         if (length < minLength)
2021         {
2022             minLength = length;
2023             minIndex = i;
2024         }
2025     }
2026     return edgeIndices[minIndex];
2027 }
2028
2029
2030 // Returns element in edgeIndices with maximum length
2031 Foam::label Foam::triSurfaceTools::maxEdge
2032 (
2033     const triSurface& surf,
2034     const labelList& edgeIndices
2035 )
2036 {
2037     scalar maxLength = -GREAT;
2038     label maxIndex = -1;
2039     forAll(edgeIndices, i)
2040     {
2041         const edge& e = surf.edges()[edgeIndices[i]];
2042
2043         scalar length =
2044             mag
2045             (
2046                 surf.localPoints()[e.end()]
2047               - surf.localPoints()[e.start()]
2048             );
2049
2050         if (length > maxLength)
2051         {
2052             maxLength = length;
2053             maxIndex = i;
2054         }
2055     }
2056     return edgeIndices[maxIndex];
2057 }
2058
2059
2060 // Merge points and reconstruct surface
2061 Foam::triSurface Foam::triSurfaceTools::mergePoints
2062 (
2063     const triSurface& surf,
2064     const scalar mergeTol
2065 )
2066 {
2067     pointField newPoints(surf.nPoints());
2068
2069     labelList pointMap(surf.nPoints());
2070
2071     bool hasMerged = Foam::mergePoints
2072     (
2073         surf.localPoints(),
2074         mergeTol,
2075         false,
2076         pointMap,
2077         newPoints
2078     );
2079
2080     if (hasMerged)
2081     {
2082         // Pack the triangles
2083
2084         // Storage for new triangles
2085         List<labelledTri> newTriangles(surf.size());
2086         label newTriangleI = 0;
2087
2088         forAll(surf, faceI)
2089         {
2090             const labelledTri& f = surf.localFaces()[faceI];
2091
2092             label newA = pointMap[f[0]];
2093             label newB = pointMap[f[1]];
2094             label newC = pointMap[f[2]];
2095
2096             if ((newA != newB) && (newA != newC) && (newB != newC))
2097             {
2098                 newTriangles[newTriangleI++] =
2099                     labelledTri(newA, newB, newC, f.region());
2100             }
2101         }
2102         newTriangles.setSize(newTriangleI);
2103
2104         return triSurface
2105         (
2106             newTriangles,
2107             surf.patches(),
2108             newPoints
2109         );
2110     }
2111     else
2112     {
2113         return surf;
2114     }
2115 }
2116
2117
2118 // Calculate normal on triangle
2119 Foam::vector Foam::triSurfaceTools::surfaceNormal
2120 (
2121     const triSurface& surf,
2122     const label nearestFaceI,
2123     const point& nearestPt
2124 )
2125 {
2126     const labelledTri& f = surf[nearestFaceI];
2127     const pointField& points = surf.points();
2128
2129     label nearType;
2130     label nearLabel;
2131     triPointRef
2132     (
2133         points[f[0]],
2134         points[f[1]],
2135         points[f[2]]
2136     ).classify(nearestPt, 1E-6, nearType, nearLabel);
2137
2138     if (nearType == triPointRef::NONE)
2139     {
2140         // Nearest to face
2141         return surf.faceNormals()[nearestFaceI];
2142     }
2143     else if (nearType == triPointRef::EDGE)
2144     {
2145         // Nearest to edge. Assume order of faceEdges same as face vertices.
2146         label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
2147
2148         // Calculate edge normal by averaging face normals
2149         const labelList& eFaces = surf.edgeFaces()[edgeI];
2150
2151         vector edgeNormal(vector::zero);
2152
2153         forAll(eFaces, i)
2154         {
2155             edgeNormal += surf.faceNormals()[eFaces[i]];
2156         }
2157         return edgeNormal/(mag(edgeNormal) + VSMALL);
2158     }
2159     else
2160     {
2161         // Nearest to point
2162         const labelledTri& localF = surf.localFaces()[nearestFaceI];
2163         return surf.pointNormals()[localF[nearLabel]];
2164     }
2165 }
2166
2167
2168 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::edgeSide
2169 (
2170     const triSurface& surf,
2171     const point& sample,
2172     const point& nearestPoint,
2173     const label edgeI
2174 )
2175 {
2176     const labelList& eFaces = surf.edgeFaces()[edgeI];
2177
2178     if (eFaces.size() != 2)
2179     {
2180         // Surface not closed.
2181         return UNKNOWN;
2182     }
2183     else
2184     {
2185         const vectorField& faceNormals = surf.faceNormals();
2186
2187         // Compare to bisector. This is actually correct since edge is
2188         // nearest so there is a knife-edge.
2189
2190         vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2191
2192         if (((sample - nearestPoint) & n) > 0)
2193         {
2194             return OUTSIDE;
2195         }
2196         else
2197         {
2198             return INSIDE;
2199         }
2200     }
2201 }
2202
2203
2204 // Calculate normal on triangle
2205 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
2206 (
2207     const triSurface& surf,
2208     const point& sample,
2209     const label nearestFaceI,   // nearest face
2210     const point& nearestPoint   // nearest point on nearest face
2211 )
2212 {
2213     const labelledTri& f = surf[nearestFaceI];
2214     const pointField& points = surf.points();
2215
2216     // Find where point is on triangle. Note tolerance needed. Is relative
2217     // one (used in comparing normalized [0..1] triangle coordinates).
2218     label nearType, nearLabel;
2219     triPointRef
2220     (
2221         points[f[0]],
2222         points[f[1]],
2223         points[f[2]]
2224     ).classify(nearestPoint, 1E-6, nearType, nearLabel);
2225
2226     if (nearType == triPointRef::NONE)
2227     {
2228         // Nearest to face interior. Use faceNormal to determine side
2229         scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
2230
2231         if (c > 0)
2232         {
2233             return OUTSIDE;
2234         }
2235         else
2236         {
2237             return INSIDE;
2238         }
2239     }
2240     else if (nearType == triPointRef::EDGE)
2241     {
2242         // Nearest to edge nearLabel. Note that this can only be a knife-edge
2243         // situation since otherwise the nearest point could never be the edge.
2244
2245         // Get the edge. Assume order of faceEdges same as face vertices.
2246         label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
2247
2248         //if (debug)
2249         //{
2250         //    // Check order of faceEdges same as face vertices.
2251         //    const edge& e = surf.edges()[edgeI];
2252         //    const labelList& meshPoints = surf.meshPoints();
2253         //    const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
2254         //
2255         //    if
2256         //    (
2257         //        meshEdge
2258         //     != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
2259         //    )
2260         //    {
2261         //        FatalErrorIn("triSurfaceTools::surfaceSide")
2262         //            << "Edge:" << edgeI << " local vertices:" << e
2263         //            << " mesh vertices:" << meshEdge
2264         //            << " not at position " << nearLabel
2265         //            << " in face " << f
2266         //            << abort(FatalError);
2267         //    }
2268         //}
2269
2270         return edgeSide(surf, sample, nearestPoint, edgeI);
2271     }
2272     else
2273     {
2274         // Nearest to point. Could use pointNormal here but is not correct.
2275         // Instead determine which edge using point is nearest and use test
2276         // above (nearType == triPointRef::EDGE).
2277
2278
2279         const labelledTri& localF = surf.localFaces()[nearestFaceI];
2280         label nearPointI = localF[nearLabel];
2281
2282         const edgeList& edges = surf.edges();
2283         const pointField& localPoints = surf.localPoints();
2284         const point& base = localPoints[nearPointI];
2285
2286         const labelList& pEdges = surf.pointEdges()[nearPointI];
2287
2288         scalar minDistSqr = Foam::sqr(GREAT);
2289         label minEdgeI = -1;
2290
2291         forAll(pEdges, i)
2292         {
2293             label edgeI = pEdges[i];
2294
2295             const edge& e = edges[edgeI];
2296
2297             label otherPointI = e.otherVertex(nearPointI);
2298
2299             // Get edge normal.
2300             vector eVec(localPoints[otherPointI] - base);
2301             scalar magEVec = mag(eVec);
2302
2303             if (magEVec > VSMALL)
2304             {
2305                 eVec /= magEVec;
2306
2307                 // Get point along vector and determine closest.
2308                 const point perturbPoint = base + eVec;
2309
2310                 scalar distSqr = Foam::magSqr(sample - perturbPoint);
2311
2312                 if (distSqr < minDistSqr)
2313                 {
2314                     minDistSqr = distSqr;
2315                     minEdgeI = edgeI;
2316                 }
2317             }
2318         }
2319
2320         if (minEdgeI == -1)
2321         {
2322             FatalErrorIn("treeDataTriSurface::getSide")
2323                 << "Problem: did not find edge closer than " << minDistSqr
2324                 << abort(FatalError);
2325         }
2326
2327         return edgeSide(surf, sample, nearestPoint, minEdgeI);
2328     }
2329 }
2330
2331
2332 // triangulation of boundaryMesh
2333 Foam::triSurface Foam::triSurfaceTools::triangulate
2334 (
2335     const polyBoundaryMesh& bMesh,
2336     const labelHashSet& includePatches,
2337     const bool verbose
2338 )
2339 {
2340     const polyMesh& mesh = bMesh.mesh();
2341
2342     // Storage for surfaceMesh. Size estimate.
2343     DynamicList<labelledTri> triangles
2344     (
2345         mesh.nFaces() - mesh.nInternalFaces()
2346     );
2347
2348     label newPatchI = 0;
2349
2350     for
2351     (
2352         labelHashSet::const_iterator iter = includePatches.begin();
2353         iter != includePatches.end();
2354         ++iter
2355     )
2356     {
2357         label patchI = iter.key();
2358
2359         const polyPatch& patch = bMesh[patchI];
2360
2361         const pointField& points = patch.points();
2362
2363         label nTriTotal = 0;
2364
2365         forAll(patch, patchFaceI)
2366         {
2367             const face& f = patch[patchFaceI];
2368
2369             faceList triFaces(f.nTriangles(points));
2370
2371             label nTri = 0;
2372
2373             f.triangles(points, nTri, triFaces);
2374
2375             forAll(triFaces, triFaceI)
2376             {
2377                 const face& f = triFaces[triFaceI];
2378
2379                 triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
2380
2381                 nTriTotal++;
2382             }
2383         }
2384
2385         if (verbose)
2386         {
2387             Pout<< patch.name() << " : generated " << nTriTotal
2388                 << " triangles from " << patch.size() << " faces with"
2389                 << " new patchid " << newPatchI << endl;
2390         }
2391
2392         newPatchI++;
2393     }
2394     triangles.shrink();
2395
2396     // Create globally numbered tri surface
2397     triSurface rawSurface(triangles, mesh.points());
2398
2399     // Create locally numbered tri surface
2400     triSurface surface
2401     (
2402         rawSurface.localFaces(),
2403         rawSurface.localPoints()
2404     );
2405
2406     // Add patch names to surface
2407     surface.patches().setSize(newPatchI);
2408
2409     newPatchI = 0;
2410
2411     for
2412     (
2413         labelHashSet::const_iterator iter = includePatches.begin();
2414         iter != includePatches.end();
2415         ++iter
2416     )
2417     {
2418         label patchI = iter.key();
2419
2420         const polyPatch& patch = bMesh[patchI];
2421
2422         surface.patches()[newPatchI].name() = patch.name();
2423
2424         surface.patches()[newPatchI].geometricType() = patch.type();
2425
2426         newPatchI++;
2427     }
2428
2429     return surface;
2430 }
2431
2432
2433 // triangulation of boundaryMesh
2434 Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
2435 (
2436     const polyBoundaryMesh& bMesh,
2437     const labelHashSet& includePatches,
2438     const bool verbose
2439 )
2440 {
2441     const polyMesh& mesh = bMesh.mesh();
2442
2443     // Storage for new points = meshpoints + face centres.
2444     const pointField& points = mesh.points();
2445     const pointField& faceCentres = mesh.faceCentres();
2446
2447     pointField newPoints(points.size() + faceCentres.size());
2448
2449     label newPointI = 0;
2450
2451     forAll(points, pointI)
2452     {
2453         newPoints[newPointI++] = points[pointI];
2454     }
2455     forAll(faceCentres, faceI)
2456     {
2457         newPoints[newPointI++] = faceCentres[faceI];
2458     }
2459
2460
2461     // Count number of faces.
2462     DynamicList<labelledTri> triangles
2463     (
2464         mesh.nFaces() - mesh.nInternalFaces()
2465     );
2466
2467     label newPatchI = 0;
2468
2469     for
2470     (
2471         labelHashSet::const_iterator iter = includePatches.begin();
2472         iter != includePatches.end();
2473         ++iter
2474     )
2475     {
2476         label patchI = iter.key();
2477
2478         const polyPatch& patch = bMesh[patchI];
2479
2480         label nTriTotal = 0;
2481
2482         forAll(patch, patchFaceI)
2483         {
2484             // Face in global coords.
2485             const face& f = patch[patchFaceI];
2486
2487             // Index in newPointI of face centre.
2488             label fc = points.size() + patchFaceI + patch.start();
2489
2490             forAll(f, fp)
2491             {
2492                 label fp1 = (fp + 1) % f.size();
2493
2494                 triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchI));
2495
2496                 nTriTotal++;
2497             }
2498         }
2499
2500         if (verbose)
2501         {
2502             Pout<< patch.name() << " : generated " << nTriTotal
2503                 << " triangles from " << patch.size() << " faces with"
2504                 << " new patchid " << newPatchI << endl;
2505         }
2506
2507         newPatchI++;
2508     }
2509     triangles.shrink();
2510
2511
2512     // Create globally numbered tri surface
2513     triSurface rawSurface(triangles, newPoints);
2514
2515     // Create locally numbered tri surface
2516     triSurface surface
2517     (
2518         rawSurface.localFaces(),
2519         rawSurface.localPoints()
2520     );
2521
2522     // Add patch names to surface
2523     surface.patches().setSize(newPatchI);
2524
2525     newPatchI = 0;
2526
2527     for
2528     (
2529         labelHashSet::const_iterator iter = includePatches.begin();
2530         iter != includePatches.end();
2531         ++iter
2532     )
2533     {
2534         label patchI = iter.key();
2535
2536         const polyPatch& patch = bMesh[patchI];
2537
2538         surface.patches()[newPatchI].name() = patch.name();
2539
2540         surface.patches()[newPatchI].geometricType() = patch.type();
2541
2542         newPatchI++;
2543     }
2544
2545     return surface;
2546 }
2547
2548
2549 Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
2550 {
2551     // Vertices in geompack notation. Note that could probably just use
2552     // pts.begin() if double precision.
2553     List<doubleScalar> geompackVertices(2*pts.size());
2554     label doubleI = 0;
2555     forAll(pts, i)
2556     {
2557         geompackVertices[doubleI++] = pts[i][0];
2558         geompackVertices[doubleI++] = pts[i][1];
2559     }
2560
2561     // Storage for triangles
2562     int m2 = 3;
2563     List<int> triangle_node(m2*3*pts.size());
2564     List<int> triangle_neighbor(m2*3*pts.size());
2565
2566     // Triangulate
2567     int nTris = 0;
2568     dtris2
2569     (
2570         pts.size(),
2571         geompackVertices.begin(),
2572         &nTris,
2573         triangle_node.begin(),
2574         triangle_neighbor.begin()
2575     );
2576
2577     // Trim
2578     triangle_node.setSize(3*nTris);
2579     triangle_neighbor.setSize(3*nTris);
2580
2581     // Convert to triSurface.
2582     List<labelledTri> faces(nTris);
2583
2584     forAll(faces, i)
2585     {
2586         faces[i] = labelledTri
2587         (
2588             triangle_node[3*i]-1,
2589             triangle_node[3*i+1]-1,
2590             triangle_node[3*i+2]-1,
2591             0
2592         );
2593     }
2594
2595     pointField points(pts.size());
2596     forAll(pts, i)
2597     {
2598         points[i][0] = pts[i][0];
2599         points[i][1] = pts[i][1];
2600         points[i][2] = 0.0;
2601     }
2602
2603     return triSurface(faces, points);
2604 }
2605
2606
2607 //// Use CGAL to do Delaunay
2608 //#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
2609 //#include <CGAL/Delaunay_triangulation_2.h>
2610 //#include <CGAL/Triangulation_vertex_base_with_info_2.h>
2611 //#include <cassert>
2612 //
2613 //struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
2614 //
2615 //typedef CGAL::Triangulation_vertex_base_with_info_2<Foam::label, K> Vb;
2616 //typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
2617 //typedef CGAL::Delaunay_triangulation_2<K,Tds> Triangulation;
2618 //
2619 //typedef Triangulation::Vertex_handle Vertex_handle;
2620 //typedef Triangulation::Vertex_iterator Vertex_iterator;
2621 //typedef Triangulation::Face_handle Face_handle;
2622 //typedef Triangulation::Finite_faces_iterator Finite_faces_iterator;
2623 //typedef Triangulation::Point Point;
2624 //Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
2625 //{
2626 //    Triangulation T;
2627 //
2628 //    // Insert 2D vertices; building triangulation
2629 //    forAll(pts, i)
2630 //    {
2631 //        const point& pt = pts[i];
2632 //
2633 //        T.insert(Point(pt[0], pt[1]));
2634 //    }
2635 //
2636 //
2637 //    // Number vertices
2638 //    // ~~~~~~~~~~~~~~~
2639 //
2640 //    label vertI = 0;
2641 //    for
2642 //    (
2643 //        Vertex_iterator it = T.vertices_begin();
2644 //        it != T.vertices_end();
2645 //        ++it
2646 //    )
2647 //    {
2648 //        it->info() = vertI++;
2649 //    }
2650 //
2651 //    // Extract faces
2652 //    // ~~~~~~~~~~~~~
2653 //
2654 //    List<labelledTri> faces(T.number_of_faces());
2655 //
2656 //    label faceI = 0;
2657 //
2658 //    for
2659 //    (
2660 //        Finite_faces_iterator fc = T.finite_faces_begin();
2661 //        fc != T.finite_faces_end();
2662 //        ++fc
2663 //    )
2664 //    {
2665 //        faces[faceI++] = labelledTri
2666 //        (
2667 //            fc->vertex(0)->info(),
2668 //            f[1] = fc->vertex(1)->info(),
2669 //            f[2] = fc->vertex(2)->info()
2670 //        );
2671 //    }
2672 //
2673 //    pointField points(pts.size());
2674 //    forAll(pts, i)
2675 //    {
2676 //        points[i][0] = pts[i][0];
2677 //        points[i][1] = pts[i][1];
2678 //        points[i][2] = 0.0;
2679 //    }
2680 //
2681 //    return triSurface(faces, points);
2682 //}
2683
2684
2685 void Foam::triSurfaceTools::calcInterpolationWeights
2686 (
2687     const triPointRef& tri,
2688     const point& p,
2689     FixedList<scalar, 3>& weights
2690 )
2691 {
2692     // calculate triangle edge vectors and triangle face normal
2693     // the 'i':th edge is opposite node i
2694     FixedList<vector, 3> edge;
2695     edge[0] = tri.c()-tri.b();
2696     edge[1] = tri.a()-tri.c();
2697     edge[2] = tri.b()-tri.a();
2698
2699     vector triangleFaceNormal = edge[1] ^ edge[2];
2700
2701     // calculate edge normal (pointing inwards)
2702     FixedList<vector, 3> normal;
2703     for(label i=0; i<3; i++)
2704     {
2705         normal[i] = triangleFaceNormal ^ edge[i];
2706         normal[i] /= mag(normal[i]) + VSMALL;
2707     }
2708
2709     weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
2710     weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
2711     weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
2712 }
2713
2714
2715 // Calculate weighting factors from samplePts to triangle it is in.
2716 // Uses linear search.
2717 void Foam::triSurfaceTools::calcInterpolationWeights
2718 (
2719     const triSurface& s,
2720     const pointField& samplePts,
2721     List<FixedList<label, 3> >& allVerts,
2722     List<FixedList<scalar, 3> >& allWeights
2723 )
2724 {
2725     allVerts.setSize(samplePts.size());
2726     allWeights.setSize(samplePts.size());
2727
2728     const pointField& points = s.points();
2729
2730     forAll(samplePts, i)
2731     {
2732         const point& samplePt = samplePts[i];
2733
2734
2735         FixedList<label, 3>& verts = allVerts[i];
2736         FixedList<scalar, 3>& weights = allWeights[i];
2737
2738         scalar minDistance = GREAT;
2739
2740         forAll(s, faceI)
2741         {
2742             const labelledTri& f = s[faceI];
2743
2744             triPointRef tri(f.tri(points));
2745
2746             pointHit nearest = tri.nearestPoint(samplePt);
2747
2748             if (nearest.hit())
2749             {
2750                 // samplePt inside triangle
2751                 verts[0] = f[0];
2752                 verts[1] = f[1];
2753                 verts[2] = f[2];
2754
2755                 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2756
2757                 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2758                 //    << " inside triangle:" << faceI
2759                 //    << " verts:" << verts
2760                 //    << " weights:" << weights
2761                 //    << endl;
2762
2763                 break;
2764             }
2765             else if (nearest.distance() < minDistance)
2766             {
2767                 minDistance = nearest.distance();
2768
2769                 // Outside triangle. Store nearest.
2770                 label nearType, nearLabel;
2771                 tri.classify
2772                 (
2773                     nearest.rawPoint(),
2774                     1E-6,               // relative tolerance
2775                     nearType,
2776                     nearLabel
2777                 );
2778
2779                 if (nearType == triPointRef::POINT)
2780                 {
2781                     verts[0] = f[nearLabel];
2782                     weights[0] = 1;
2783                     verts[1] = -1;
2784                     weights[1] = -GREAT;
2785                     verts[2] = -1;
2786                     weights[2] = -GREAT;
2787
2788                     //Pout<< "calcScalingFactors : samplePt:" << samplePt
2789                     //    << " distance:" << nearest.distance()
2790                     //    << " from point:" << points[f[nearLabel]]
2791                     //    << endl;
2792                 }
2793                 else if (nearType == triPointRef::EDGE)
2794                 {
2795                     verts[0] = f[nearLabel];
2796                     verts[1] = f[f.fcIndex(nearLabel)];
2797                     verts[2] = -1;
2798
2799                     const point& p0 = points[verts[0]];
2800                     const point& p1 = points[verts[1]];
2801
2802                     scalar s = min
2803                     (
2804                         1,
2805                         max
2806                         (
2807                             0,
2808                             mag(nearest.rawPoint()-p0)/mag(p1-p0)
2809                         )
2810                     );
2811
2812                     // Interpolate
2813                     weights[0] = 1-s;
2814                     weights[1] = s;
2815                     weights[2] = -GREAT;
2816
2817                     //Pout<< "calcScalingFactors : samplePt:" << samplePt
2818                     //    << " distance:" << nearest.distance()
2819                     //    << " from edge:" << p0 << p1 << " s:" << s
2820                     //    << endl;
2821                 }
2822                 else
2823                 {
2824                     // triangle. Can only happen because of truncation errors.
2825                     verts[0] = f[0];
2826                     verts[1] = f[1];
2827                     verts[2] = f[2];
2828
2829                     calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2830
2831                     break;
2832                 }
2833             }
2834         }
2835     }
2836 }
2837
2838
2839 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2840 // Tracking:
2841
2842
2843 // Test point on surface to see if is on face,edge or point.
2844 Foam::surfaceLocation Foam::triSurfaceTools::classify
2845 (
2846     const triSurface& s,
2847     const label triI,
2848     const point& trianglePoint
2849 )
2850 {
2851     surfaceLocation nearest;
2852
2853     // Nearest point could be on point or edge. Retest.
2854     label index, elemType;
2855     //bool inside =
2856     triPointRef(s[triI].tri(s.points())).classify
2857     (
2858         trianglePoint,
2859         1E-6,
2860         elemType,
2861         index
2862     );
2863
2864     nearest.setPoint(trianglePoint);
2865
2866     if (elemType == triPointRef::NONE)
2867     {
2868         nearest.setHit();
2869         nearest.setIndex(triI);
2870         nearest.elementType() = triPointRef::NONE;
2871     }
2872     else if (elemType == triPointRef::EDGE)
2873     {
2874         nearest.setMiss();
2875         nearest.setIndex(s.faceEdges()[triI][index]);
2876         nearest.elementType() = triPointRef::EDGE;
2877     }
2878     else //if (elemType == triPointRef::POINT)
2879     {
2880         nearest.setMiss();
2881         nearest.setIndex(s.localFaces()[triI][index]);
2882         nearest.elementType() = triPointRef::POINT;
2883     }
2884
2885     return nearest;
2886 }
2887
2888
2889 Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge
2890 (
2891     const triSurface& s,
2892     const surfaceLocation& start,
2893     const surfaceLocation& end,
2894     const plane& cutPlane
2895 )
2896 {
2897     // Start off from starting point
2898     surfaceLocation nearest = start;
2899     nearest.setMiss();
2900
2901     // See if in same triangle as endpoint. If so snap.
2902     snapToEnd(s, end, nearest);
2903
2904     if (!nearest.hit())
2905     {
2906         // Not yet at end point
2907
2908         if (start.elementType() == triPointRef::NONE)
2909         {
2910             // Start point is inside triangle. Trivial cases already handled
2911             // above.
2912
2913             // end point is on edge or point so cross currrent triangle to
2914             // see which edge is cut.
2915
2916             nearest = cutEdge
2917             (
2918                 s,
2919                 start.index(),          // triangle
2920                 -1,                     // excludeEdge
2921                 -1,                     // excludePoint
2922                 start.rawPoint(),
2923                 cutPlane,
2924                 end.rawPoint()
2925             );
2926             nearest.elementType() = triPointRef::EDGE;
2927             nearest.triangle() = start.index();
2928             nearest.setMiss();
2929         }
2930         else if (start.elementType() == triPointRef::EDGE)
2931         {
2932             // Pick connected triangle that is most in direction.
2933             const labelList& eFaces = s.edgeFaces()[start.index()];
2934
2935             nearest = visitFaces
2936             (
2937                 s,
2938                 eFaces,
2939                 start,
2940                 start.index(),      // excludeEdgeI
2941                 -1,                 // excludePointI
2942                 end,
2943                 cutPlane
2944             );
2945         }
2946         else    // start.elementType() == triPointRef::POINT
2947         {
2948             const labelList& pFaces = s.pointFaces()[start.index()];
2949
2950             nearest = visitFaces
2951             (
2952                 s,
2953                 pFaces,
2954                 start,
2955                 -1,                 // excludeEdgeI
2956                 start.index(),      // excludePointI
2957                 end,
2958                 cutPlane
2959             );
2960         }
2961         snapToEnd(s, end, nearest);
2962     }
2963     return nearest;
2964 }
2965
2966
2967 void Foam::triSurfaceTools::track
2968 (
2969     const triSurface& s,
2970     const surfaceLocation& endInfo,
2971     const plane& cutPlane,
2972     surfaceLocation& hitInfo
2973 )
2974 {
2975     //OFstream str("track.obj");
2976     //label vertI = 0;
2977     //meshTools::writeOBJ(str, hitInfo.rawPoint());
2978     //vertI++;
2979
2980     // Track across surface.
2981     while (true)
2982     {
2983         //Pout<< "Tracking from:" << nl
2984         //    << "    " << hitInfo.info()
2985         //    << endl;
2986
2987         hitInfo = trackToEdge
2988         (
2989             s,
2990             hitInfo,
2991             endInfo,
2992             cutPlane
2993         );
2994
2995         //meshTools::writeOBJ(str, hitInfo.rawPoint());
2996         //vertI++;
2997         //str<< "l " << vertI-1 << ' ' << vertI << nl;
2998
2999         //Pout<< "Tracked to:" << nl
3000         //    << "    " << hitInfo.info() << endl;
3001
3002         if (hitInfo.hit() || hitInfo.triangle() == -1)
3003         {
3004             break;
3005         }
3006     }
3007 }
3008
3009
3010 // ************************************************************************* //