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