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-2009 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:" << UIndirectList<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:" << UIndirectList<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     // points in face in same order as edge
 465     bool edgeOrder = (f[f.fcIndex(startIndex)] == e.end());
 466
 467     // Get faces using edge in sorted order. (sorted such that walking
 468     // around them in anti-clockwise order corresponds to edge vector
 469     // acc. to right-hand rule)
 470     const labelList& eFaces = surf.sortedEdgeFaces()[sharedEdgeI];
 471
 472     // Get position of face in sorted edge faces
 473     label faceIndex = findIndex(eFaces, firstFaceI);
 474
 475     if (edgeOrder)
 476     {
 477         // Get face before firstFaceI
 478         return eFaces[eFaces.rcIndex(faceIndex)];
 479     }
 480     else
 481     {
 482         // Get face after firstFaceI
 483         return eFaces[eFaces.fcIndex(faceIndex)];
 484     }
 485 }
 486
 487
 488 // Calculate (inward pointing) normals on edges shared by faces in faceToEdge and
 489 // averages them to pointNormals.
 490 void calcPointVecs
 491 (
 492     const triSurface& surf,
 493     const Map<label>& faceToEdge,
 494     const Map<label>& faceToPoint,
 495     vectorField& borderPointVec
 496 )
 497 {
 498     const labelListList& sortedEdgeFaces = surf.sortedEdgeFaces();
 499     const edgeList& edges = surf.edges();
 500     const pointField& points = surf.localPoints();
 501
 502     boolList edgeDone(surf.nEdges(), false);
 503
 504     for
 505     (
 506         Map<label>::const_iterator iter = faceToEdge.begin();
 507         iter != faceToEdge.end();
 508         ++iter
 509     )
 510     {
 511         label edgeI = iter();
 512
 513         if (!edgeDone[edgeI])
 514         {
 515             edgeDone[edgeI] = true;
 516
 517             // Get the two connected faces in sorted order
 518             // Note: should have stored this when walking ...
 519
 520             label face0I = -1;
 521             label face1I = -1;
 522
 523             const labelList& eFaces = sortedEdgeFaces[edgeI];
 524
 525             forAll(eFaces, i)
 526             {
 527                 label faceI = eFaces[i];
 528
 529                 if (faceToEdge.found(faceI))
 530                 {
 531                     if (face0I == -1)
 532                     {
 533                         face0I = faceI;
 534                     }
 535                     else if (face1I == -1)
 536                     {
 537                         face1I = faceI;
 538
 539                         break;
 540                     }
 541                 }
 542             }
 543
 544             if (face0I == -1 && face1I == -1)
 545             {
 546                 Info<< "Writing surface to errorSurf.ftr" << endl;
 547
 548                 surf.write("errorSurf.ftr");
 549
 550                 FatalErrorIn("calcPointVecs")
 551                     << "Cannot find two faces using border edge " << edgeI
 552                     << " verts:" << edges[edgeI]
 553                     << " eFaces:" << eFaces << endl
 554                     << "face0I:" << face0I
 555                     << " face1I:" << face1I << nl
 556                     << "faceToEdge:" << faceToEdge << nl
 557                     << "faceToPoint:" << faceToPoint
 558                     << "Written surface to errorSurf.ftr"
 559                     << abort(FatalError);
 560             }
 561
 562             // Now we have edge and the two faces in counter-clockwise order
 563             // as seen from edge vector. Calculate normal.
 564             const edge& e = edges[edgeI];
 565             vector eVec = e.vec(points);
 566
 567             // Determine vector as average of the vectors in the two faces.
 568             // If there is only one face available use only one vector.
 569             vector midVec(vector::zero);
 570
 571             if (face0I != -1)
 572             {
 573                 label v0 = triSurfaceTools::oppositeVertex(surf, face0I, edgeI);
 574                 vector e0 = (points[v0] - points[e.start()]) ^ eVec;
 575                 e0 /= mag(e0);
 576                 midVec = e0;
 577             }
 578
 579             if (face1I != -1)
 580             {
 581                 label v1 = triSurfaceTools::oppositeVertex(surf, face1I, edgeI);
 582                 vector e1 = (points[e.start()] - points[v1]) ^ eVec;
 583                 e1 /= mag(e1);
 584                 midVec += e1;
 585             }
 586
 587             scalar magMidVec = mag(midVec);
 588
 589             if (magMidVec > SMALL)
 590             {
 591                 midVec /= magMidVec;
 592
 593                 // Average to pointVec
 594                 borderPointVec[e.start()] += midVec;
 595                 borderPointVec[e.end()] += midVec;
 596             }
 597         }
 598     }
 599 }
 600
 601
 602 // Renumbers vertices (of triangles in faceToEdge) of which the pointMap is
 603 // not -1.
 604 void renumberFaces
 605 (
 606     const triSurface& surf,
 607     const labelList& pointMap,
 608     const Map<label>& faceToEdge,
 609     List<labelledTri>& newTris
 610 )
 611 {
 612     for
 613     (
 614         Map<label>::const_iterator iter = faceToEdge.begin();
 615         iter != faceToEdge.end();
 616         ++iter
 617     )
 618     {
 619         label faceI = iter.key();
 620
 621         const labelledTri& f = surf.localFaces()[faceI];
 622
 623         forAll(f, fp)
 624         {
 625             if (pointMap[f[fp]] != -1)
 626             {
 627                 newTris[faceI][fp] = pointMap[f[fp]];
 628             }
 629         }
 630     }
 631 }
 632
 633
 634 // Split all borderEdges that don't have borderPoint. Return true if split
 635 // anything.
 636 bool splitBorderEdges
 637 (
 638     triSurface& surf,
 639     const boolList& borderEdge,
 640     const labelList& borderPoint
 641 )
 642 {
 643     labelList edgesToBeSplit(surf.nEdges());
 644     label nSplit = 0;
 645
 646     forAll(borderEdge, edgeI)
 647     {
 648         if (borderEdge[edgeI])
 649         {
 650             const edge& e = surf.edges()[edgeI];
 651
 652             if (borderPoint[e.start()] == -1 && borderPoint[e.end()] == -1)
 653             {
 654                 // None of the points of the edge is borderPoint. Split edge
 655                 // to introduce border point.
 656                 edgesToBeSplit[nSplit++] = edgeI;
 657             }
 658         }
 659     }
 660     edgesToBeSplit.setSize(nSplit);
 661
 662     if (nSplit > 0)
 663     {
 664         Info<< "Splitting surface along " << nSplit << " borderEdges that don't"
 665             << " neighbour other borderEdges" << nl << endl;
 666
 667         surf = triSurfaceTools::greenRefine(surf, edgesToBeSplit);
 668
 669         return true;
 670     }
 671     else
 672     {
 673         Info<< "No edges to be split" <<nl << endl;
 674
 675         return false;
 676     }
 677 }
 678
 679
 680 // Main program:
 681
 682 int main(int argc, char *argv[])
 683 {
 684     argList::noParallel();
 685     argList::validArgs.clear();
 686
 687     argList::validArgs.append("surface file");
 688     argList::validArgs.append("output surface file");
 689     argList::validOptions.insert("debug", "");
 690
 691     argList args(argc, argv);
 692
 693     fileName inSurfName(args.additionalArgs()[0]);
 694     fileName outSurfName(args.additionalArgs()[1]);
 695     bool debug = args.optionFound("debug");
 696
 697
 698     Info<< "Reading surface from " << inSurfName << endl;
 699
 700     triSurface surf(inSurfName);
 701
 702     // Make sure sortedEdgeFaces is calculated correctly
 703     testSortedEdgeFaces(surf);
 704
 705     // Get all quad connected edges. These are seen as borders when walking.
 706     boolList borderEdge(surf.nEdges(), false);
 707     markBorderEdges(debug, surf, borderEdge);
 708
 709     // Points on two sides connected to borderEdges are called borderPoints and
 710     // will be duplicated. borderPoint contains label of newly introduced vertex.
 711     labelList borderPoint(surf.nPoints(), -1);
 712     markBorderPoints(debug, surf, borderEdge, borderPoint);
 713
 714     // Split edges where there would be no borderPoint to duplicate.
 715     splitBorderEdges(surf, borderEdge, borderPoint);
 716
 717
 718     Info<< "Writing split surface to " << outSurfName << nl << endl;
 719     surf.write(outSurfName);
 720     Info<< "Finished writing surface to " << outSurfName << nl << endl;
 721
 722
 723     // Last iteration values.
 724     label nOldBorderEdges = -1;
 725     label nOldBorderPoints = -1;
 726
 727     label iteration = 0;
 728
 729     do
 730     {
 731         // Redo borderEdge/borderPoint calculation.
 732         boolList borderEdge(surf.nEdges(), false);
 733
 734         label nBorderEdges = markBorderEdges(debug, surf, borderEdge);
 735
 736         if (nBorderEdges == 0)
 737         {
 738             Info<< "Found no border edges. Exiting." << nl << nl;
 739
 740             break;
 741         }
 742
 743         // Label of newly introduced duplicate.
 744         labelList borderPoint(surf.nPoints(), -1);
 745
 746         label nBorderPoints =
 747             markBorderPoints
 748             (
 749                 debug,
 750                 surf,
 751                 borderEdge,
 752                 borderPoint
 753             );
 754
 755         if (nBorderPoints == 0)
 756         {
 757             Info<< "Found no border points. Exiting." << nl << nl;
 758
 759             break;
 760         }
 761
 762         Info<< "Found:\n"
 763             << "    border edges :" << nBorderEdges << nl
 764             << "    border points:" << nBorderPoints << nl
 765             << endl;
 766
 767         if
 768         (
 769             nBorderPoints == nOldBorderPoints
 770          && nBorderEdges == nOldBorderEdges
 771         )
 772         {
 773             Info<< "Stopping since number of border edges and point is same"
 774                 << " as in previous iteration" << nl << endl;
 775
 776             break;
 777         }
 778
 779         //
 780         // Define splitLine as a series of connected borderEdges. Find start
 781         // of one (as edge and point on edge)
 782         //
 783
 784         label startEdgeI = -1;
 785         label startPointI = -1;
 786
 787         forAll(borderEdge, edgeI)
 788         {
 789             if (borderEdge[edgeI])
 790             {
 791                 const edge& e = surf.edges()[edgeI];
 792
 793                 if ((borderPoint[e[0]] != -1) && (borderPoint[e[1]] == -1))
 794                 {
 795                     startEdgeI = edgeI;
 796                     startPointI = e[0];
 797
 798                     break;
 799                 }
 800                 else if ((borderPoint[e[0]] == -1) && (borderPoint[e[1]] != -1))
 801                 {
 802                     startEdgeI = edgeI;
 803                     startPointI = e[1];
 804
 805                     break;
 806                 }
 807             }
 808         }
 809
 810         if (startEdgeI == -1)
 811         {
 812             Info<< "Cannot find starting point of splitLine\n" << endl;
 813
 814             break;
 815         }
 816
 817         // Pick any face using edge to start from.
 818         const labelList& eFaces = surf.edgeFaces()[startEdgeI];
 819
 820         label firstFaceI = eFaces[0];
 821
 822         // Find second face which is from same surface i.e. has outwards
 823         // pointing normal as well (actually bit more complex than this)
 824         label secondFaceI = sharedFace(surf, firstFaceI, startEdgeI);
 825
 826         Info<< "Starting local walk from:" << endl
 827             << "    edge :" << startEdgeI << endl
 828             << "    point:" << startPointI << endl
 829             << "    face0:" << firstFaceI << endl
 830             << "    face1:" << secondFaceI << endl
 831             << endl;
 832
 833         // From face on border edge to edge.
 834         Map<label> faceToEdge(2*nBorderEdges);
 835         // From face connected to border point (but not border edge) to point.
 836         Map<label> faceToPoint(2*nBorderPoints);
 837
 838         faceToEdge.insert(firstFaceI, startEdgeI);
 839
 840         walkSplitLine
 841         (
 842             surf,
 843             borderEdge,
 844             borderPoint,
 845
 846             firstFaceI,
 847             startEdgeI,
 848             startPointI,
 849
 850             faceToEdge,
 851             faceToPoint
 852         );
 853
 854         faceToEdge.insert(secondFaceI, startEdgeI);
 855
 856         walkSplitLine
 857         (
 858             surf,
 859             borderEdge,
 860             borderPoint,
 861
 862             secondFaceI,
 863             startEdgeI,
 864             startPointI,
 865
 866             faceToEdge,
 867             faceToPoint
 868         );
 869
 870         Info<< "Finished local walk and visited" << nl
 871             << "    border edges                :" << faceToEdge.size() << nl
 872             << "    border points(but not edges):" << faceToPoint.size() << nl
 873             << endl;
 874
 875         if (debug)
 876         {
 877             dumpFaces("faceToEdge.obj", surf, faceToEdge);
 878             dumpFaces("faceToPoint.obj", surf, faceToPoint);
 879         }
 880
 881         //
 882         // Create coordinates for borderPoints by duplicating the existing
 883         // point and then slightly shifting it inwards. To determine the
 884         // inwards direction get the average normal of both connectedFaces on
 885         // the edge and then interpolate this to the (border)point.
 886         //
 887
 888         vectorField borderPointVec(surf.nPoints(), vector(GREAT, GREAT, GREAT));
 889
 890         calcPointVecs(surf, faceToEdge, faceToPoint, borderPointVec);
 891
 892
 893         // New position. Start off from copy of old points.
 894         pointField newPoints(surf.localPoints());
 895         newPoints.setSize(newPoints.size() + nBorderPoints);
 896
 897         forAll(borderPoint, pointI)
 898         {
 899             label newPointI = borderPoint[pointI];
 900
 901             if (newPointI != -1)
 902             {
 903                 scalar minLen = minEdgeLen(surf, pointI);
 904
 905                 vector n = borderPointVec[pointI];
 906                 n /= mag(n);
 907
 908                 newPoints[newPointI] = newPoints[pointI] + 0.1 * minLen * n;
 909             }
 910         }
 911
 912
 913         //
 914         // Renumber all faces in connectedFaces
 915         //
 916
 917         // Start off from copy of faces.
 918         List<labelledTri> newTris(surf.size());
 919
 920         forAll(surf, faceI)
 921         {
 922             newTris[faceI] = surf.localFaces()[faceI];
 923
 924             newTris[faceI].region() = surf[faceI].region();
 925         }
 926
 927         // Renumber all faces in faceToEdge
 928         renumberFaces(surf, borderPoint, faceToEdge, newTris);
 929         // Renumber all faces in faceToPoint
 930         renumberFaces(surf, borderPoint, faceToPoint, newTris);
 931
 932
 933         // Check if faces use unmoved points.
 934         forAll(newTris, faceI)
 935         {
 936             const labelledTri& f = newTris[faceI];
 937
 938             forAll(f, fp)
 939             {
 940                 const point& pt = newPoints[f[fp]];
 941
 942                 if (mag(pt) >= GREAT/2)
 943                 {
 944                     Info<< "newTri:" << faceI << " verts:" << f
 945                         << " vert:" << f[fp] << " point:" << pt << endl;
 946                 }
 947             }
 948         }
 949
 950
 951         surf = triSurface(newTris, surf.patches(), newPoints);
 952
 953         if (debug || (iteration != 0 && (iteration % 20) == 0))
 954         {
 955             Info<< "Writing surface to " << outSurfName << nl << endl;
 956
 957             surf.write(outSurfName);
 958
 959             Info<< "Finished writing surface to " << outSurfName << nl << endl;
 960         }
 961
 962         // Save prev iteration values.
 963         nOldBorderEdges = nBorderEdges;
 964         nOldBorderPoints = nBorderPoints;
 965
 966         iteration++;
 967     }
 968     while (true);
 969
 970     Info<< "Writing final surface to " << outSurfName << nl << endl;
 971
 972     surf.write(outSurfName);
 973
 974     Info << "End\n" << endl;
 975
 976     return 0;
 977 }
 978
 979
 980 // ************************************************************************* //