applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.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     Takes multiply connected surface and tries to split surface at
  27     multiply connected edges by duplicating points. Introduces concept of
  28     - borderEdge. Edge with 4 faces connected to it.
  29     - borderPoint. Point connected to exactly 2 borderEdges.
  30     - borderLine. Connected list of borderEdges.
  31
  32     By duplicating borderPoints this will split 'borderLines'. As a
  33     preprocessing step it can detect borderEdges without any borderPoints
  34     and explicitly split these triangles.
  35
  36     The problems in this algorithm are:
  37     - determining which two (of the four) faces form a surface. Done by walking
  38       face-edge-face while keeping and edge or point on the borderEdge
  39       borderPoint.
  40     - determining the outwards pointing normal to be used to slightly offset the
  41       duplicated point.
  42
  43     Uses sortedEdgeFaces quite a bit.
  44
  45     Is tested on simple borderLines resulting from extracting a surface
  46     from a hex mesh. Will quite possibly go wrong on more complicated border
  47     lines (i.e. ones forming a loop).
  48
  49     Dumps surface every so often since might take a long time to complete.
  50
  51 \*---------------------------------------------------------------------------*/
  52
  53 #include "argList.H"
  54 #include "triSurface.H"
  55 #include "OFstream.H"
  56 #include "ListOps.H"
  57 #include "triSurfaceTools.H"
  58
  59 using namespace Foam;
  60
  61 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  62
  63 void writeOBJ(Ostream& os, const pointField& pts)
  64 {
  65     forAll(pts, i)
  66     {
  67         const point& pt = pts[i];
  68
  69         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
  70     }
  71 }
  72
  73
  74 void dumpPoints(const triSurface& surf, const labelList& borderPoint)
  75 {
  76     fileName fName("borderPoints.obj");
  77
  78     Info<< "Dumping borderPoints as Lightwave .obj file to " << fName
  79         << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
  80
  81     OFstream os(fName);
  82
  83     forAll(borderPoint, pointI)
  84     {
  85         if (borderPoint[pointI] != -1)
  86         {
  87             const point& pt = surf.localPoints()[pointI];
  88
  89             os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
  90         }
  91     }
  92 }
  93
  94
  95 void dumpEdges(const triSurface& surf, const boolList& borderEdge)
  96 {
  97     fileName fName("borderEdges.obj");
  98
  99     Info<< "Dumping borderEdges as Lightwave .obj file to " << fName
 100         << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
 101
 102     OFstream os(fName);
 103
 104     writeOBJ(os, surf.localPoints());
 105
 106     forAll(borderEdge, edgeI)
 107     {
 108         if (borderEdge[edgeI])
 109         {
 110             const edge& e = surf.edges()[edgeI];
 111
 112             os << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
 113         }
 114     }
 115 }
 116
 117
 118 void dumpFaces
 119 (
 120     const fileName& fName,
 121     const triSurface& surf,
 122     const Map<label>& connectedFaces
 123 )
 124 {
 125     Info<< "Dumping connectedFaces as Lightwave .obj file to " << fName
 126         << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
 127
 128     OFstream os(fName);
 129
 130     for
 131     (
 132         Map<label>::const_iterator iter = connectedFaces.begin();
 133         iter != connectedFaces.end();
 134         ++iter
 135     )
 136     {
 137         const labelledTri& f = surf.localFaces()[iter.key()];
 138
 139         point ctr(f.centre(surf.localPoints()));
 140
 141         os << "v " << ctr.x() << ' ' << ctr.y() << ' ' << ctr.z() << endl;
 142     }
 143 }
 144
 145
 146 void testSortedEdgeFaces(const triSurface& surf)
 147 {
 148     const labelListList& edgeFaces = surf.edgeFaces();
 149     const labelListList& sortedEdgeFaces = surf.sortedEdgeFaces();
 150
 151     forAll(edgeFaces, edgeI)
 152     {
 153         const labelList& myFaces = edgeFaces[edgeI];
 154         const labelList& sortMyFaces = sortedEdgeFaces[edgeI];
 155
 156         forAll(myFaces, i)
 157         {
 158             if (findIndex(sortMyFaces, myFaces[i]) == -1)
 159             {
 160                 FatalErrorIn("testSortedEdgeFaces") << abort(FatalError);
 161             }
 162         }
 163         forAll(sortMyFaces, i)
 164         {
 165             if (findIndex(myFaces, sortMyFaces[i]) == -1)
 166             {
 167                 FatalErrorIn("testSortedEdgeFaces") << abort(FatalError);
 168             }
 169         }
 170     }
 171 }
 172
 173
 174 // Mark off all border edges. Return number of border edges.
 175 label markBorderEdges
 176 (
 177     const bool debug,
 178     const triSurface& surf,
 179     boolList& borderEdge
 180 )
 181 {
 182     label nBorderEdges = 0;
 183
 184     const labelListList& edgeFaces = surf.edgeFaces();
 185
 186     forAll(edgeFaces, edgeI)
 187     {
 188         if (edgeFaces[edgeI].size() == 4)
 189         {
 190             borderEdge[edgeI] = true;
 191
 192             nBorderEdges++;
 193         }
 194     }
 195
 196     if (debug)
 197     {
 198         dumpEdges(surf, borderEdge);
 199     }
 200
 201     return nBorderEdges;
 202 }
 203
 204
 205 // Mark off all border points. Return number of border points. Border points
 206 // marked by setting value to newly introduced point.
 207 label markBorderPoints
 208 (
 209     const bool debug,
 210     const triSurface& surf,
 211     const boolList& borderEdge,
 212     labelList& borderPoint
 213 )
 214 {
 215     label nPoints = surf.nPoints();
 216
 217     const labelListList& pointEdges = surf.pointEdges();
 218
 219     forAll(pointEdges, pointI)
 220     {
 221         const labelList& pEdges = pointEdges[pointI];
 222
 223         label nBorderEdges = 0;
 224
 225         forAll(pEdges, i)
 226         {
 227             if (borderEdge[pEdges[i]])
 228             {
 229                 nBorderEdges++;
 230             }
 231         }
 232
 233         if (nBorderEdges == 2 && borderPoint[pointI] == -1)
 234         {
 235             borderPoint[pointI] = nPoints++;
 236         }
 237     }
 238
 239     label nBorderPoints = nPoints - surf.nPoints();
 240
 241     if (debug)
 242     {
 243         dumpPoints(surf, borderPoint);
 244     }
 245
 246     return nBorderPoints;
 247 }
 248
 249
 250 // Get minumum length of edges connected to pointI
 251 // Serves to get some local length scale.
 252 scalar minEdgeLen(const triSurface& surf, const label pointI)
 253 {
 254     const pointField& points = surf.localPoints();
 255
 256     const labelList& pEdges = surf.pointEdges()[pointI];
 257
 258     scalar minLen = GREAT;
 259
 260     forAll(pEdges, i)
 261     {
 262         label edgeI = pEdges[i];
 263
 264         scalar len = surf.edges()[edgeI].mag(points);
 265
 266         if (len < minLen)
 267         {
 268             minLen = len;
 269         }
 270     }
 271     return minLen;
 272 }
 273
 274
 275 // Find edge among edgeLabels with endpoints v0,v1
 276 label findEdge
 277 (
 278     const triSurface& surf,
 279     const labelList& edgeLabels,
 280     const label v0,
 281     const label v1
 282 )
 283 {
 284     forAll(edgeLabels, i)
 285     {
 286         label edgeI = edgeLabels[i];
 287
 288         const edge& e = surf.edges()[edgeI];
 289
 290         if
 291         (
 292             (
 293                 e.start() == v0
 294              && e.end() == v1
 295             )
 296          || (
 297                 e.start() == v1
 298              && e.end() == v0
 299             )
 300         )
 301         {
 302             return edgeI;
 303         }
 304     }
 305
 306
 307     FatalErrorIn("findEdge") << "Cannot find edge with labels " << v0
 308         << ' ' << v1 << " in candidates " << edgeLabels
 309         << " with vertices:" << IndirectList<edge>(surf.edges(), edgeLabels)()
 310         << abort(FatalError);
 311
 312     return -1;
 313 }
 314
 315
 316 // Get the other edge connected to pointI on faceI.
 317 label otherEdge
 318 (
 319     const triSurface& surf,
 320     const label faceI,
 321     const label otherEdgeI,
 322     const label pointI
 323 )
 324 {
 325     const labelList& fEdges = surf.faceEdges()[faceI];
 326
 327     forAll(fEdges, i)
 328     {
 329         label edgeI = fEdges[i];
 330
 331         const edge& e = surf.edges()[edgeI];
 332
 333         if
 334         (
 335             edgeI != otherEdgeI
 336          && (
 337                 e.start() == pointI
 338              || e.end() == pointI
 339             )
 340         )
 341         {
 342             return edgeI;
 343         }
 344     }
 345
 346     FatalErrorIn("otherEdge") << "Cannot find other edge on face " << faceI
 347         << " verts:" << surf.localPoints()[faceI]
 348         << " connected to point " << pointI
 349         << " faceEdges:" << IndirectList<edge>(surf.edges(), fEdges)()
 350         << abort(FatalError);
 351
 352     return -1;
 353 }
 354
 355
 356 // Starting from startPoint on startEdge on startFace walk along border
 357 // and insert faces along the way. Walk keeps always one point or one edge
 358 // on the border line.
 359 void walkSplitLine
 360 (
 361     const triSurface& surf,
 362     const boolList& borderEdge,
 363     const labelList& borderPoint,
 364
 365     const label startFaceI,
 366     const label startEdgeI,     // is border edge
 367     const label startPointI,    // is border point
 368
 369     Map<label>& faceToEdge,
 370     Map<label>& faceToPoint
 371 )
 372 {
 373     label faceI = startFaceI;
 374     label edgeI = startEdgeI;
 375     label pointI = startPointI;
 376
 377     do
 378     {
 379         //
 380         // Stick to pointI and walk face-edge-face until back on border edge.
 381         //
 382
 383         do
 384         {
 385             // Cross face to next edge.
 386             edgeI = otherEdge(surf, faceI, edgeI, pointI);
 387
 388             if (borderEdge[edgeI])
 389             {
 390                 if (!faceToEdge.insert(faceI, edgeI))
 391                 {
 392                     // Was already visited.
 393                     return;
 394                 }
 395                 else
 396                 {
 397                     // First visit to this borderEdge. We're back on borderline.
 398                     break;
 399                 }
 400             }
 401             else if (!faceToPoint.insert(faceI, pointI))
 402             {
 403                 // Was already visited.
 404                 return;
 405             }
 406
 407             // Cross edge to other face
 408             const labelList& eFaces = surf.edgeFaces()[edgeI];
 409
 410             if (eFaces.size() != 2)
 411             {
 412                 FatalErrorIn("walkSplitPoint") << abort(FatalError);
 413             }
 414
 415             if (eFaces[0] == faceI)
 416             {
 417                 faceI = eFaces[1];
 418             }
 419             else if (eFaces[1] == faceI)
 420             {
 421                 faceI = eFaces[0];
 422             }
 423             else
 424             {
 425                 FatalErrorIn("walkSplitPoint") << abort(FatalError);
 426             }
 427         }
 428         while (true);
 429
 430
 431         //
 432         // Back on border edge. Cross to other point on edge.
 433         //
 434
 435         pointI = surf.edges()[edgeI].otherVertex(pointI);
 436
 437         if (borderPoint[pointI] == -1)
 438         {
 439             return;
 440         }
 441
 442     }
 443     while (true);
 444 }
 445
 446
 447 // Find second face which is from same surface i.e. has points on the
 448 // shared edge in reverse order.
 449 label sharedFace
 450 (
 451     const triSurface& surf,
 452     const label firstFaceI,
 453     const label sharedEdgeI
 454 )
 455 {
 456     // Find ordering of face points in edge.
 457
 458     const edge& e = surf.edges()[sharedEdgeI];
 459
 460     const labelledTri& f = surf.localFaces()[firstFaceI];
 461
 462     label startIndex = findIndex(f, e.start());
 463
 464     bool edgeOrder;
 465
 466     if (f[(startIndex + 1) % f.size()] == e.end())
 467     {
 468         // points in face in same order as edge
 469         edgeOrder = true;
 470     }
 471     else
 472     {
 473         // points in face in reverse order as edge
 474         edgeOrder = false;
 475     }
 476
 477     // Get faces using edge in sorted order. (sorted such that walking
 478     // around them in anti-clockwise order corresponds to edge vector
 479     // acc. to right-hand rule)
 480     const labelList& eFaces = surf.sortedEdgeFaces()[sharedEdgeI];
 481
 482     // Get position of face in sorted edge faces
 483     label faceIndex = findIndex(eFaces, firstFaceI);
 484
 485     if (edgeOrder)
 486     {
 487         // Get face before firstFaceI
 488         if (faceIndex == 0)
 489         {
 490             return eFaces[eFaces.size() - 1];
 491         }
 492         else
 493         {
 494             return eFaces[faceIndex - 1];
 495         }
 496     }
 497     else
 498     {
 499         // Get face after firstFaceI
 500         return eFaces[(faceIndex+1) % eFaces.size()];
 501     }
 502 }
 503
 504
 505 // Calculate (inward pointing) normals on edges shared by faces in faceToEdge and
 506 // averages them to pointNormals.
 507 void calcPointVecs
 508 (
 509     const triSurface& surf,
 510     const Map<label>& faceToEdge,
 511     const Map<label>& faceToPoint,
 512     vectorField& borderPointVec
 513 )
 514 {
 515     const labelListList& sortedEdgeFaces = surf.sortedEdgeFaces();
 516     const edgeList& edges = surf.edges();
 517     const pointField& points = surf.localPoints();
 518
 519     boolList edgeDone(surf.nEdges(), false);
 520
 521     for
 522     (
 523         Map<label>::const_iterator iter = faceToEdge.begin();
 524         iter != faceToEdge.end();
 525         ++iter
 526     )
 527     {
 528         label edgeI = iter();
 529
 530         if (!edgeDone[edgeI])
 531         {
 532             edgeDone[edgeI] = true;
 533
 534             // Get the two connected faces in sorted order
 535             // Note: should have stored this when walking ...
 536
 537             label face0I = -1;
 538             label face1I = -1;
 539
 540             const labelList& eFaces = sortedEdgeFaces[edgeI];
 541
 542             forAll(eFaces, i)
 543             {
 544                 label faceI = eFaces[i];
 545
 546                 if (faceToEdge.found(faceI))
 547                 {
 548                     if (face0I == -1)
 549                     {
 550                         face0I = faceI;
 551                     }
 552                     else if (face1I == -1)
 553                     {
 554                         face1I = faceI;
 555
 556                         break;
 557                     }
 558                 }
 559             }
 560
 561             if (face0I == -1 && face1I == -1)
 562             {
 563                 Info<< "Writing surface to errorSurf.ftr" << endl;
 564
 565                 surf.write("errorSurf.ftr");
 566
 567                 FatalErrorIn("calcPointVecs")
 568                     << "Cannot find two faces using border edge " << edgeI
 569                     << " verts:" << edges[edgeI]
 570                     << " eFaces:" << eFaces << endl
 571                     << "face0I:" << face0I
 572                     << " face1I:" << face1I << nl
 573                     << "faceToEdge:" << faceToEdge << nl
 574                     << "faceToPoint:" << faceToPoint
 575                     << "Written surface to errorSurf.ftr"
 576                     << abort(FatalError);
 577             }
 578
 579             // Now we have edge and the two faces in counter-clockwise order
 580             // as seen from edge vector. Calculate normal.
 581             const edge& e = edges[edgeI];
 582             vector eVec = e.vec(points);
 583
 584             // Determine vector as average of the vectors in the two faces.
 585             // If there is only one face available use only one vector.
 586             vector midVec(vector::zero);
 587
 588             if (face0I != -1)
 589             {
 590                 label v0 = triSurfaceTools::oppositeVertex(surf, face0I, edgeI);
 591                 vector e0 = (points[v0] - points[e.start()]) ^ eVec;
 592                 e0 /= mag(e0);
 593                 midVec = e0;
 594             }
 595
 596             if (face1I != -1)
 597             {
 598                 label v1 = triSurfaceTools::oppositeVertex(surf, face1I, edgeI);
 599                 vector e1 = (points[e.start()] - points[v1]) ^ eVec;
 600                 e1 /= mag(e1);
 601                 midVec += e1;
 602             }
 603
 604             scalar magMidVec = mag(midVec);
 605
 606             if (magMidVec > SMALL)
 607             {
 608                 midVec /= magMidVec;
 609
 610                 // Average to pointVec
 611                 borderPointVec[e.start()] += midVec;
 612                 borderPointVec[e.end()] += midVec;
 613             }
 614         }
 615     }
 616 }
 617
 618
 619 // Renumbers vertices (of triangles in faceToEdge) of which the pointMap is
 620 // not -1.
 621 void renumberFaces
 622 (
 623     const triSurface& surf,
 624     const labelList& pointMap,
 625     const Map<label>& faceToEdge,
 626     List<labelledTri>& newTris
 627 )
 628 {
 629     for
 630     (
 631         Map<label>::const_iterator iter = faceToEdge.begin();
 632         iter != faceToEdge.end();
 633         ++iter
 634     )
 635     {
 636         label faceI = iter.key();
 637
 638         const labelledTri& f = surf.localFaces()[faceI];
 639
 640         forAll(f, fp)
 641         {
 642             if (pointMap[f[fp]] != -1)
 643             {
 644                 newTris[faceI][fp] = pointMap[f[fp]];
 645             }
 646         }
 647     }
 648 }
 649
 650
 651 // Split all borderEdges that don't have borderPoint. Return true if split
 652 // anything.
 653 bool splitBorderEdges
 654 (
 655     triSurface& surf,
 656     const boolList& borderEdge,
 657     const labelList& borderPoint
 658 )
 659 {
 660     labelList edgesToBeSplit(surf.nEdges());
 661     label nSplit = 0;
 662
 663     forAll(borderEdge, edgeI)
 664     {
 665         if (borderEdge[edgeI])
 666         {
 667             const edge& e = surf.edges()[edgeI];
 668
 669             if (borderPoint[e.start()] == -1 && borderPoint[e.end()] == -1)
 670             {
 671                 // None of the points of the edge is borderPoint. Split edge
 672                 // to introduce border point.
 673                 edgesToBeSplit[nSplit++] = edgeI;
 674             }
 675         }
 676     }
 677     edgesToBeSplit.setSize(nSplit);
 678
 679     if (nSplit > 0)
 680     {
 681         Info<< "Splitting surface along " << nSplit << " borderEdges that don't"
 682             << " neighbour other borderEdges" << nl << endl;
 683
 684         surf = triSurfaceTools::greenRefine(surf, edgesToBeSplit);
 685
 686         return true;
 687     }
 688     else
 689     {
 690         Info<< "No edges to be split" <<nl << endl;
 691
 692         return false;
 693     }
 694 }
 695
 696
 697 // Main program:
 698
 699 int main(int argc, char *argv[])
 700 {
 701     argList::noParallel();
 702     argList::validArgs.clear();
 703
 704     argList::validArgs.append("surface file");
 705     argList::validArgs.append("output surface file");
 706     argList::validOptions.insert("debug", "");
 707
 708     argList args(argc, argv);
 709
 710     fileName inSurfName(args.additionalArgs()[0]);
 711     fileName outSurfName(args.additionalArgs()[1]);
 712     bool debug = args.options().found("debug");
 713
 714
 715     Info<< "Reading surface from " << inSurfName << endl;
 716
 717     triSurface surf(inSurfName);
 718
 719     // Make sure sortedEdgeFaces is calculated correctly
 720     testSortedEdgeFaces(surf);
 721
 722     // Get all quad connected edges. These are seen as borders when walking.
 723     boolList borderEdge(surf.nEdges(), false);
 724     markBorderEdges(debug, surf, borderEdge);
 725
 726     // Points on two sides connected to borderEdges are called borderPoints and
 727     // will be duplicated. borderPoint contains label of newly introduced vertex.
 728     labelList borderPoint(surf.nPoints(), -1);
 729     markBorderPoints(debug, surf, borderEdge, borderPoint);
 730
 731     // Split edges where there would be no borderPoint to duplicate.
 732     splitBorderEdges(surf, borderEdge, borderPoint);
 733
 734
 735     Info<< "Writing split surface to " << outSurfName << nl << endl;
 736     surf.write(outSurfName);
 737     Info<< "Finished writing surface to " << outSurfName << nl << endl;
 738
 739
 740     // Last iteration values.
 741     label nOldBorderEdges = -1;
 742     label nOldBorderPoints = -1;
 743
 744     label iteration = 0;
 745
 746     do
 747     {
 748         // Redo borderEdge/borderPoint calculation.
 749         boolList borderEdge(surf.nEdges(), false);
 750
 751         label nBorderEdges = markBorderEdges(debug, surf, borderEdge);
 752
 753         if (nBorderEdges == 0)
 754         {
 755             Info<< "Found no border edges. Exiting." << nl << nl;
 756
 757             break;
 758         }
 759
 760         // Label of newly introduced duplicate.
 761         labelList borderPoint(surf.nPoints(), -1);
 762
 763         label nBorderPoints =
 764             markBorderPoints
 765             (
 766                 debug,
 767                 surf,
 768                 borderEdge,
 769                 borderPoint
 770             );
 771
 772         if (nBorderPoints == 0)
 773         {
 774             Info<< "Found no border points. Exiting." << nl << nl;
 775
 776             break;
 777         }
 778
 779         Info<< "Found:\n"
 780             << "    border edges :" << nBorderEdges << nl
 781             << "    border points:" << nBorderPoints << nl
 782             << endl;
 783
 784         if
 785         (
 786             nBorderPoints == nOldBorderPoints
 787          && nBorderEdges == nOldBorderEdges
 788         )
 789         {
 790             Info<< "Stopping since number of border edges and point is same"
 791                 << " as in previous iteration" << nl << endl;
 792
 793             break;
 794         }
 795
 796         //
 797         // Define splitLine as a series of connected borderEdges. Find start
 798         // of one (as edge and point on edge)
 799         //
 800
 801         label startEdgeI = -1;
 802         label startPointI = -1;
 803
 804         forAll(borderEdge, edgeI)
 805         {
 806             if (borderEdge[edgeI])
 807             {
 808                 const edge& e = surf.edges()[edgeI];
 809
 810                 if ((borderPoint[e[0]] != -1) && (borderPoint[e[1]] == -1))
 811                 {
 812                     startEdgeI = edgeI;
 813                     startPointI = e[0];
 814
 815                     break;
 816                 }
 817                 else if ((borderPoint[e[0]] == -1) && (borderPoint[e[1]] != -1))
 818                 {
 819                     startEdgeI = edgeI;
 820                     startPointI = e[1];
 821
 822                     break;
 823                 }
 824             }
 825         }
 826
 827         if (startEdgeI == -1)
 828         {
 829             Info<< "Cannot find starting point of splitLine\n" << endl;
 830
 831             break;
 832         }
 833
 834         // Pick any face using edge to start from.
 835         const labelList& eFaces = surf.edgeFaces()[startEdgeI];
 836
 837         label firstFaceI = eFaces[0];
 838
 839         // Find second face which is from same surface i.e. has outwards
 840         // pointing normal as well (actually bit more complex than this)
 841         label secondFaceI = sharedFace(surf, firstFaceI, startEdgeI);
 842
 843         Info<< "Starting local walk from:" << endl
 844             << "    edge :" << startEdgeI << endl
 845             << "    point:" << startPointI << endl
 846             << "    face0:" << firstFaceI << endl
 847             << "    face1:" << secondFaceI << endl
 848             << endl;
 849
 850         // From face on border edge to edge.
 851         Map<label> faceToEdge(2*nBorderEdges);
 852         // From face connected to border point (but not border edge) to point.
 853         Map<label> faceToPoint(2*nBorderPoints);
 854
 855         faceToEdge.insert(firstFaceI, startEdgeI);
 856
 857         walkSplitLine
 858         (
 859             surf,
 860             borderEdge,
 861             borderPoint,
 862
 863             firstFaceI,
 864             startEdgeI,
 865             startPointI,
 866
 867             faceToEdge,
 868             faceToPoint
 869         );
 870
 871         faceToEdge.insert(secondFaceI, startEdgeI);
 872
 873         walkSplitLine
 874         (
 875             surf,
 876             borderEdge,
 877             borderPoint,
 878
 879             secondFaceI,
 880             startEdgeI,
 881             startPointI,
 882
 883             faceToEdge,
 884             faceToPoint
 885         );
 886
 887         Info<< "Finished local walk and visited" << nl
 888             << "    border edges                :" << faceToEdge.size() << nl
 889             << "    border points(but not edges):" << faceToPoint.size() << nl
 890             << endl;
 891
 892         if (debug)
 893         {
 894             dumpFaces("faceToEdge.obj", surf, faceToEdge);
 895             dumpFaces("faceToPoint.obj", surf, faceToPoint);
 896         }
 897
 898         //
 899         // Create coordinates for borderPoints by duplicating the existing
 900         // point and then slightly shifting it inwards. To determine the
 901         // inwards direction get the average normal of both connectedFaces on
 902         // the edge and then interpolate this to the (border)point.
 903         //
 904
 905         vectorField borderPointVec(surf.nPoints(), vector(GREAT, GREAT, GREAT));
 906
 907         calcPointVecs(surf, faceToEdge, faceToPoint, borderPointVec);
 908
 909
 910         // New position. Start off from copy of old points.
 911         pointField newPoints(surf.localPoints());
 912         newPoints.setSize(newPoints.size() + nBorderPoints);
 913
 914         forAll(borderPoint, pointI)
 915         {
 916             label newPointI = borderPoint[pointI];
 917
 918             if (newPointI != -1)
 919             {
 920                 scalar minLen = minEdgeLen(surf, pointI);
 921
 922                 vector n = borderPointVec[pointI];
 923                 n /= mag(n);
 924
 925                 newPoints[newPointI] = newPoints[pointI] + 0.1 * minLen * n;
 926             }
 927         }
 928
 929
 930         //
 931         // Renumber all faces in connectedFaces
 932         //
 933
 934         // Start off from copy of faces.
 935         List<labelledTri> newTris(surf.size());
 936
 937         forAll(surf, faceI)
 938         {
 939             newTris[faceI] = surf.localFaces()[faceI];
 940
 941             newTris[faceI].region() = surf[faceI].region();
 942         }
 943
 944         // Renumber all faces in faceToEdge
 945         renumberFaces(surf, borderPoint, faceToEdge, newTris);
 946         // Renumber all faces in faceToPoint
 947         renumberFaces(surf, borderPoint, faceToPoint, newTris);
 948
 949
 950         // Check if faces use unmoved points.
 951         forAll(newTris, faceI)
 952         {
 953             const labelledTri& f = newTris[faceI];
 954
 955             forAll(f, fp)
 956             {
 957                 const point& pt = newPoints[f[fp]];
 958
 959                 if (mag(pt) >= GREAT/2)
 960                 {
 961                     Info<< "newTri:" << faceI << " verts:" << f
 962                         << " vert:" << f[fp] << " point:" << pt << endl;
 963                 }
 964             }
 965         }
 966
 967
 968         surf = triSurface(newTris, surf.patches(), newPoints);
 969
 970         if (debug || (iteration != 0 && (iteration % 20) == 0))
 971         {
 972             Info<< "Writing surface to " << outSurfName << nl << endl;
 973
 974             surf.write(outSurfName);
 975
 976             Info<< "Finished writing surface to " << outSurfName << nl << endl;
 977         }
 978
 979         // Save prev iteration values.
 980         nOldBorderEdges = nBorderEdges;
 981         nOldBorderPoints = nBorderPoints;
 982
 983         iteration++;
 984     }
 985     while (true);
 986
 987     Info<< "Writing final surface to " << outSurfName << nl << endl;
 988
 989     surf.write(outSurfName);
 990
 991     Info << "End\n" << endl;
 992
 993     return 0;
 994 }
 995
 996
 997 // ************************************************************************* //