applications/utilities/surface/surfaceClean/collapseBase.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 \*---------------------------------------------------------------------------*/
  25
  26 #include "collapseBase.H"
  27 #include <meshTools/triSurfaceTools.H>
  28 #include <OpenFOAM/argList.H>
  29 #include <OpenFOAM/OFstream.H>
  30 #include <OpenFOAM/SubList.H>
  31 #include <OpenFOAM/labelPair.H>
  32 #include <meshTools/meshTools.H>
  33 #include <OpenFOAM/OSspecific.H>
  34
  35 using namespace Foam;
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 // Dump collapse region to .obj file
  40 static void writeRegionOBJ
  41 (
  42     const triSurface& surf,
  43     const label regionI,
  44     const labelList& collapseRegion,
  45     const labelList& outsideVerts
  46 )
  47 {
  48     fileName dir("regions");
  49
  50     mkDir(dir);
  51     fileName regionName(dir / "region_" + name(regionI) + ".obj");
  52
  53     Pout<< "Dumping region " << regionI << " to file " << regionName << endl;
  54
  55     boolList include(surf.size(), false);
  56
  57     forAll(collapseRegion, faceI)
  58     {
  59         if (collapseRegion[faceI] == regionI)
  60         {
  61             include[faceI] = true;
  62         }
  63     }
  64
  65     labelList pointMap, faceMap;
  66
  67     triSurface regionSurf(surf.subsetMesh(include, pointMap, faceMap));
  68
  69     //Pout<< "Region " << regionI << " surface:" << nl;
  70     //regionSurf.writeStats(Pout);
  71
  72     regionSurf.write(regionName);
  73
  74
  75     // Dump corresponding outside vertices.
  76     fileName pointsName(dir / "regionPoints_" + name(regionI) + ".obj");
  77
  78     Pout<< "Dumping region " << regionI << " points to file " << pointsName
  79         << endl;
  80
  81     OFstream str(pointsName);
  82
  83     forAll(outsideVerts, i)
  84     {
  85         meshTools::writeOBJ(str, surf.localPoints()[outsideVerts[i]]);
  86     }
  87 }
  88
  89
  90 // Split triangle into multiple triangles because edge e being split
  91 // into multiple edges.
  92 static void splitTri
  93 (
  94     const labelledTri& f,
  95     const edge& e,
  96     const labelList& splitPoints,
  97     DynamicList<labelledTri>& tris
  98 )
  99 {
 100     label oldNTris = tris.size();
 101
 102     label fp = findIndex(f, e[0]);
 103     label fp1 = (fp+1)%3;
 104     label fp2 = (fp1+1)%3;
 105
 106     if (f[fp1] == e[1])
 107     {
 108         // Split triangle along fp to fp1
 109         tris.append(labelledTri(f[fp2], f[fp], splitPoints[0], f.region()));
 110
 111         for (label i = 1; i < splitPoints.size(); i++)
 112         {
 113             tris.append
 114             (
 115                 labelledTri
 116                 (
 117                     f[fp2],
 118                     splitPoints[i-1],
 119                     splitPoints[i],
 120                     f.region()
 121                 )
 122             );
 123         }
 124
 125         tris.append
 126         (
 127             labelledTri
 128             (
 129                 f[fp2],
 130                 splitPoints[splitPoints.size()-1],
 131                 f[fp1],
 132                 f.region()
 133             )
 134         );
 135     }
 136     else if (f[fp2] == e[1])
 137     {
 138         // Split triangle along fp2 to fp. (Reverse order of splitPoints)
 139
 140         tris.append
 141         (
 142             labelledTri
 143             (
 144                 f[fp1],
 145                 f[fp2],
 146                 splitPoints[splitPoints.size()-1],
 147                 f.region()
 148             )
 149         );
 150
 151         for (label i = splitPoints.size()-1; i > 0; --i)
 152         {
 153             tris.append
 154             (
 155                 labelledTri
 156                 (
 157                     f[fp1],
 158                     splitPoints[i],
 159                     splitPoints[i-1],
 160                     f.region()
 161                 )
 162             );
 163         }
 164
 165         tris.append
 166         (
 167             labelledTri
 168             (
 169                 f[fp1],
 170                 splitPoints[0],
 171                 f[fp],
 172                 f.region()
 173             )
 174         );
 175     }
 176     else
 177     {
 178         FatalErrorIn("splitTri")
 179             << "Edge " << e << " not part of triangle " << f
 180             << " fp:" << fp
 181             << " fp1:" << fp1
 182             << " fp2:" << fp2
 183             << abort(FatalError);
 184     }
 185
 186     Pout<< "Split face " << f << " along edge " << e
 187         << " into triangles:" << endl;
 188
 189     for (label i = oldNTris; i < tris.size(); i++)
 190     {
 191         Pout<< "   " << tris[i] << nl;
 192     }
 193 }
 194
 195
 196 // Insert scalar into sortedVerts/sortedWeights so the weights are in
 197 // incrementing order.
 198 static bool insertSorted
 199 (
 200     const label vertI,
 201     const scalar weight,
 202
 203     labelList& sortedVerts,
 204     scalarField& sortedWeights
 205 )
 206 {
 207     if (findIndex(sortedVerts, vertI) != -1)
 208     {
 209         FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
 210             << " which is already in list of sorted vertices "
 211             << sortedVerts << abort(FatalError);
 212     }
 213
 214     if (weight <= 0 || weight >= 1)
 215     {
 216         FatalErrorIn("insertSorted") << "Trying to insert vertex " << vertI
 217             << " with illegal weight " << weight
 218             << " into list of sorted vertices "
 219             << sortedVerts << abort(FatalError);
 220     }
 221
 222
 223     label insertI = sortedVerts.size();
 224
 225     forAll(sortedVerts, sortedI)
 226     {
 227         scalar w = sortedWeights[sortedI];
 228
 229         if (mag(w - weight) < SMALL)
 230         {
 231             WarningIn("insertSorted")
 232                 << "Trying to insert weight " << weight << " which is close to"
 233                 << " existing weight " << w << " in " << sortedWeights
 234                 << endl;
 235         }
 236
 237         if (w > weight)
 238         {
 239             // Start inserting before sortedI.
 240             insertI = sortedI;
 241             break;
 242         }
 243     }
 244
 245
 246     label sz = sortedWeights.size();
 247
 248     sortedWeights.setSize(sz + 1);
 249     sortedVerts.setSize(sz + 1);
 250
 251     // Leave everything up to (not including) insertI intact.
 252
 253     // Make space by copying from insertI up.
 254     for (label i = sz-1; i >= insertI; --i)
 255     {
 256         sortedWeights[i+1] = sortedWeights[i];
 257         sortedVerts[i+1] = sortedVerts[i];
 258     }
 259     sortedWeights[insertI] = weight;
 260     sortedVerts[insertI] = vertI;
 261
 262     return true;
 263 }
 264
 265
 266 // Mark all faces that are going to be collapsed.
 267 // faceToEdge: per face -1 or the base edge of the face.
 268 static void markCollapsedFaces
 269 (
 270     const triSurface& surf,
 271     const scalar minLen,
 272     labelList& faceToEdge
 273 )
 274 {
 275     faceToEdge.setSize(surf.size());
 276     faceToEdge = -1;
 277
 278     const pointField& localPoints = surf.localPoints();
 279     const labelListList& edgeFaces = surf.edgeFaces();
 280
 281     forAll(edgeFaces, edgeI)
 282     {
 283         const edge& e = surf.edges()[edgeI];
 284
 285         const labelList& eFaces = surf.edgeFaces()[edgeI];
 286
 287         forAll(eFaces, i)
 288         {
 289             label faceI = eFaces[i];
 290
 291             // Check distance of vertex to edge.
 292             label opposite0 =
 293                 triSurfaceTools::oppositeVertex
 294                 (
 295                     surf,
 296                     faceI,
 297                     edgeI
 298                 );
 299
 300             pointHit pHit =
 301                 e.line(localPoints).nearestDist
 302                 (
 303                     localPoints[opposite0]
 304                 );
 305
 306             if (pHit.hit() && pHit.distance() < minLen)
 307             {
 308                 // Remove faceI and split all other faces using this
 309                 // edge. This is done by 'replacing' the edgeI with the
 310                 // opposite0 vertex
 311                 Pout<< "Splitting face " << faceI << " since distance "
 312                     << pHit.distance()
 313                     << " from vertex " << opposite0
 314                     << " to edge " << edgeI
 315                     << "  points "
 316                     << localPoints[e[0]]
 317                     << localPoints[e[1]]
 318                     << " is too small" << endl;
 319
 320                 // Mark face as being collapsed
 321                 if (faceToEdge[faceI] != -1)
 322                 {
 323                     FatalErrorIn("markCollapsedFaces")
 324                         << "Cannot collapse face " << faceI << " since "
 325                         << " is marked to be collapsed both to edge "
 326                         << faceToEdge[faceI] << " and " << edgeI
 327                         << abort(FatalError);
 328                 }
 329
 330                 faceToEdge[faceI] = edgeI;
 331             }
 332         }
 333     }
 334 }
 335
 336
 337 // Recurse through collapsed faces marking all of them with regionI (in
 338 // collapseRegion)
 339 static void markRegion
 340 (
 341     const triSurface& surf,
 342     const labelList& faceToEdge,
 343     const label regionI,
 344     const label faceI,
 345     labelList& collapseRegion
 346 )
 347 {
 348     if (faceToEdge[faceI] == -1 || collapseRegion[faceI] != -1)
 349     {
 350         FatalErrorIn("markRegion")
 351             << "Problem : crossed into uncollapsed/regionized face"
 352             << abort(FatalError);
 353     }
 354
 355     collapseRegion[faceI] = regionI;
 356
 357     // Recurse across edges to collapsed neighbours
 358
 359     const labelList& fEdges = surf.faceEdges()[faceI];
 360
 361     forAll(fEdges, fEdgeI)
 362     {
 363         label edgeI = fEdges[fEdgeI];
 364
 365         const labelList& eFaces = surf.edgeFaces()[edgeI];
 366
 367         forAll(eFaces, i)
 368         {
 369             label nbrFaceI = eFaces[i];
 370
 371             if (faceToEdge[nbrFaceI] != -1)
 372             {
 373                 if (collapseRegion[nbrFaceI] == -1)
 374                 {
 375                     markRegion
 376                     (
 377                         surf,
 378                         faceToEdge,
 379                         regionI,
 380                         nbrFaceI,
 381                         collapseRegion
 382                     );
 383                 }
 384                 else if (collapseRegion[nbrFaceI] != regionI)
 385                 {
 386                     FatalErrorIn("markRegion")
 387                         << "Edge:" << edgeI << " between face " << faceI
 388                         << " with region " << regionI
 389                         << " and face " << nbrFaceI
 390                         << " with region " << collapseRegion[nbrFaceI]
 391                         << endl;
 392                 }
 393             }
 394         }
 395     }
 396 }
 397
 398
 399 // Mark every face with region (in collapseRegion) (or -1).
 400 // Return number of regions.
 401 static label markRegions
 402 (
 403     const triSurface& surf,
 404     const labelList& faceToEdge,
 405     labelList& collapseRegion
 406 )
 407 {
 408     label regionI = 0;
 409
 410     forAll(faceToEdge, faceI)
 411     {
 412         if (collapseRegion[faceI] == -1 && faceToEdge[faceI] != -1)
 413         {
 414             Pout<< "markRegions : Marking region:" << regionI
 415                 << " starting from face " << faceI << endl;
 416
 417             // Collapsed face. Mark connected region with current region number
 418             markRegion(surf, faceToEdge, regionI++, faceI, collapseRegion);
 419         }
 420     }
 421     return regionI;
 422 }
 423
 424
 425 // Type of region.
 426 // -1  : edge inbetween uncollapsed faces.
 427 // -2  : edge inbetween collapsed faces
 428 // >=0 : edge inbetween uncollapsed and collapsed region. Returns region.
 429 static label edgeType
 430 (
 431     const triSurface& surf,
 432     const labelList& collapseRegion,
 433     const label edgeI
 434 )
 435 {
 436     const labelList& eFaces = surf.edgeFaces()[edgeI];
 437
 438     // Detect if edge is inbetween collapseRegion and non-collapse face
 439     bool usesUncollapsed = false;
 440     label usesRegion = -1;
 441
 442     forAll(eFaces, i)
 443     {
 444         label faceI = eFaces[i];
 445
 446         label region = collapseRegion[faceI];
 447
 448         if (region == -1)
 449         {
 450             usesUncollapsed = true;
 451         }
 452         else if (usesRegion == -1)
 453         {
 454             usesRegion = region;
 455         }
 456         else if (usesRegion != region)
 457         {
 458             FatalErrorIn("edgeRegion") << abort(FatalError);
 459         }
 460         else
 461         {
 462             // Equal regions.
 463         }
 464     }
 465
 466     if (usesUncollapsed)
 467     {
 468         if (usesRegion == -1)
 469         {
 470             // uncollapsed faces only.
 471             return -1;
 472         }
 473         else
 474         {
 475             // between uncollapsed and collapsed.
 476             return usesRegion;
 477         }
 478     }
 479     else
 480     {
 481         if (usesRegion == -1)
 482         {
 483             FatalErrorIn("edgeRegion") << abort(FatalError);
 484             return -2;
 485         }
 486         else
 487         {
 488             return -2;
 489         }
 490     }
 491 }
 492
 493
 494 // Get points on outside edge of region (= outside points)
 495 static labelListList getOutsideVerts
 496 (
 497     const triSurface& surf,
 498     const labelList& collapseRegion,
 499     const label nRegions
 500 )
 501 {
 502     const labelListList& edgeFaces = surf.edgeFaces();
 503
 504     // Per region all the outside vertices.
 505     labelListList outsideVerts(nRegions);
 506
 507     forAll(edgeFaces, edgeI)
 508     {
 509         // Detect if edge is inbetween collapseRegion and non-collapse face
 510         label regionI = edgeType(surf, collapseRegion, edgeI);
 511
 512         if (regionI >= 0)
 513         {
 514             // Edge borders both uncollapsed face and collapsed face on region
 515             // usesRegion.
 516
 517             const edge& e = surf.edges()[edgeI];
 518
 519             labelList& regionVerts = outsideVerts[regionI];
 520
 521             // Add both edge points to regionVerts.
 522             forAll(e, eI)
 523             {
 524                 label v = e[eI];
 525
 526                 if (findIndex(regionVerts, v) == -1)
 527                 {
 528                     label sz = regionVerts.size();
 529                     regionVerts.setSize(sz+1);
 530                     regionVerts[sz] = v;
 531                 }
 532             }
 533         }
 534     }
 535
 536     return outsideVerts;
 537 }
 538
 539
 540 // n^2 search for furthest removed point pair.
 541 static labelPair getSpanPoints
 542 (
 543     const triSurface& surf,
 544     const labelList& outsideVerts
 545 )
 546 {
 547     const pointField& localPoints = surf.localPoints();
 548
 549     scalar maxDist = -GREAT;
 550     labelPair maxPair;
 551
 552     forAll(outsideVerts, i)
 553     {
 554         label v0 = outsideVerts[i];
 555
 556         for (label j = i+1; j < outsideVerts.size(); j++)
 557         {
 558             label v1 = outsideVerts[j];
 559
 560             scalar d = mag(localPoints[v0] - localPoints[v1]);
 561
 562             if (d > maxDist)
 563             {
 564                 maxDist = d;
 565                 maxPair[0] = v0;
 566                 maxPair[1] = v1;
 567             }
 568         }
 569     }
 570
 571     return maxPair;
 572 }
 573
 574
 575 // Project all non-span points onto the span edge.
 576 static void projectNonSpanPoints
 577 (
 578     const triSurface& surf,
 579     const labelList& outsideVerts,
 580     const labelPair& spanPair,
 581     labelList& sortedVertices,
 582     scalarField& sortedWeights
 583 )
 584 {
 585     const point& p0 = surf.localPoints()[spanPair[0]];
 586     const point& p1 = surf.localPoints()[spanPair[1]];
 587
 588     forAll(outsideVerts, i)
 589     {
 590         label v = outsideVerts[i];
 591
 592         if (v != spanPair[0] && v != spanPair[1])
 593         {
 594             // Is a non-span point. Project onto spanning edge.
 595
 596             pointHit pHit =
 597                 linePointRef(p0, p1).nearestDist
 598                 (
 599                     surf.localPoints()[v]
 600                 );
 601
 602             if (!pHit.hit())
 603             {
 604                 FatalErrorIn("projectNonSpanPoints")
 605                     << abort(FatalError);
 606             }
 607
 608             scalar w = mag(pHit.hitPoint() - p0) / mag(p1 - p0);
 609
 610             insertSorted(v, w, sortedVertices, sortedWeights);
 611         }
 612     }
 613 }
 614
 615
 616 // Slice part of the orderVertices (and optionally reverse) for this edge.
 617 static void getSplitVerts
 618 (
 619     const triSurface& surf,
 620     const label regionI,
 621     const labelPair& spanPoints,
 622     const labelList& orderedVerts,
 623     const scalarField& orderedWeights,
 624     const label edgeI,
 625
 626     labelList& splitVerts,
 627     scalarField& splitWeights
 628 )
 629 {
 630     const edge& e = surf.edges()[edgeI];
 631     const label sz = orderedVerts.size();
 632
 633     if (e[0] == spanPoints[0])
 634     {
 635         // Edge in same order as spanPoints&orderedVerts. Keep order.
 636
 637         if (e[1] == spanPoints[1])
 638         {
 639             // Copy all.
 640             splitVerts = orderedVerts;
 641             splitWeights = orderedWeights;
 642         }
 643         else
 644         {
 645             // Copy upto (but not including) e[1]
 646             label i1 = findIndex(orderedVerts, e[1]);
 647             splitVerts = SubList<label>(orderedVerts, i1, 0);
 648             splitWeights = SubList<scalar>(orderedWeights, i1, 0);
 649         }
 650     }
 651     else if (e[0] == spanPoints[1])
 652     {
 653         // Reverse.
 654
 655         if (e[1] == spanPoints[0])
 656         {
 657             // Copy all.
 658             splitVerts = orderedVerts;
 659             reverse(splitVerts);
 660             splitWeights = orderedWeights;
 661             reverse(splitWeights);
 662         }
 663         else
 664         {
 665             // Copy downto (but not including) e[1]
 666
 667             label i1 = findIndex(orderedVerts, e[1]);
 668             splitVerts = SubList<label>(orderedVerts, sz-(i1+1), i1+1);
 669             reverse(splitVerts);
 670             splitWeights = SubList<scalar>(orderedWeights, sz-(i1+1), i1+1);
 671             reverse(splitWeights);
 672         }
 673     }
 674     else if (e[1] == spanPoints[0])
 675     {
 676         // Reverse.
 677
 678         // Copy upto (but not including) e[0]
 679
 680         label i0 = findIndex(orderedVerts, e[0]);
 681         splitVerts = SubList<label>(orderedVerts, i0, 0);
 682         reverse(splitVerts);
 683         splitWeights = SubList<scalar>(orderedWeights, i0, 0);
 684         reverse(splitWeights);
 685     }
 686     else if (e[1] == spanPoints[1])
 687     {
 688         // Copy from (but not including) e[0] to end
 689
 690         label i0 = findIndex(orderedVerts, e[0]);
 691         splitVerts = SubList<label>(orderedVerts, sz-(i0+1), i0+1);
 692         splitWeights = SubList<scalar>(orderedWeights, sz-(i0+1), i0+1);
 693     }
 694     else
 695     {
 696         label i0 = findIndex(orderedVerts, e[0]);
 697         label i1 = findIndex(orderedVerts, e[1]);
 698
 699         if (i0 == -1 || i1 == -1)
 700         {
 701             FatalErrorIn("getSplitVerts")
 702                 << "Did not find edge in projected vertices." << nl
 703                 << "region:" << regionI << nl
 704                 << "spanPoints:" << spanPoints
 705                 << "  coords:" << surf.localPoints()[spanPoints[0]]
 706                 << surf.localPoints()[spanPoints[1]] << nl
 707                 << "edge:" << edgeI
 708                 << "  verts:" << e
 709                 << "  coords:" << surf.localPoints()[e[0]]
 710                 << surf.localPoints()[e[1]] << nl
 711                 << "orderedVerts:" << orderedVerts << nl
 712                 << abort(FatalError);
 713         }
 714
 715         if (i0 < i1)
 716         {
 717             splitVerts = SubList<label>(orderedVerts, i1-i0-1, i0+1);
 718             splitWeights = SubList<scalar>(orderedWeights, i1-i0-1, i0+1);
 719         }
 720         else
 721         {
 722             splitVerts = SubList<label>(orderedVerts, i0-i1-1, i1+1);
 723             reverse(splitVerts);
 724             splitWeights = SubList<scalar>(orderedWeights, i0-i1-1, i1+1);
 725             reverse(splitWeights);
 726         }
 727     }
 728 }
 729
 730
 731 label collapseBase(triSurface& surf, const scalar minLen)
 732 {
 733     label nTotalSplit = 0;
 734
 735     label iter = 0;
 736
 737     while (true)
 738     {
 739         // Detect faces to collapse
 740         // ~~~~~~~~~~~~~~~~~~~~~~~~
 741
 742         // -1 or edge the face is collapsed onto.
 743         labelList faceToEdge(surf.size(), -1);
 744
 745         // Calculate faceToEdge (face collapses)
 746         markCollapsedFaces(surf, minLen, faceToEdge);
 747
 748
 749         // Find regions of connected collapsed faces
 750         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 751
 752         // per face -1 or region
 753         labelList collapseRegion(surf.size(), -1);
 754
 755         label nRegions = markRegions(surf, faceToEdge, collapseRegion);
 756
 757         Pout<< "Detected " << nRegions << " regions of faces to be collapsed"
 758             << nl << endl;
 759
 760         // Pick up all vertices on outside of region
 761         labelListList outsideVerts
 762         (
 763             getOutsideVerts(surf, collapseRegion, nRegions)
 764         );
 765
 766         // For all regions determine maximum distance between points
 767         List<labelPair> spanPoints(nRegions);
 768         labelListList orderedVertices(nRegions);
 769         List<scalarField> orderedWeights(nRegions);
 770
 771         forAll(spanPoints, regionI)
 772         {
 773             spanPoints[regionI] = getSpanPoints(surf, outsideVerts[regionI]);
 774
 775             Pout<< "For region " << regionI << " found extrema at points "
 776                 << surf.localPoints()[spanPoints[regionI][0]]
 777                 << surf.localPoints()[spanPoints[regionI][1]]
 778                 << endl;
 779
 780             // Project all non-span points onto the span edge.
 781             projectNonSpanPoints
 782             (
 783                 surf,
 784                 outsideVerts[regionI],
 785                 spanPoints[regionI],
 786                 orderedVertices[regionI],
 787                 orderedWeights[regionI]
 788             );
 789
 790             Pout<< "For region:" << regionI
 791                 << " span:" << spanPoints[regionI]
 792                 << " orderedVerts:" << orderedVertices[regionI]
 793                 << " orderedWeights:" << orderedWeights[regionI]
 794                 << endl;
 795
 796             writeRegionOBJ
 797             (
 798                 surf,
 799                 regionI,
 800                 collapseRegion,
 801                 outsideVerts[regionI]
 802             );
 803
 804             Pout<< endl;
 805         }
 806
 807
 808
 809         // Actually split the edges
 810         // ~~~~~~~~~~~~~~~~~~~~~~~~
 811
 812
 813         const List<labelledTri>& localFaces = surf.localFaces();
 814         const edgeList& edges = surf.edges();
 815
 816         label nSplit = 0;
 817
 818         // Storage for new triangles.
 819         DynamicList<labelledTri> newTris(surf.size());
 820
 821         // Whether face has been dealt with (either copied/split or deleted)
 822         boolList faceHandled(surf.size(), false);
 823
 824
 825         forAll(edges, edgeI)
 826         {
 827             const edge& e = edges[edgeI];
 828
 829             // Detect if edge is inbetween collapseRegion and non-collapse face
 830             label regionI = edgeType(surf, collapseRegion, edgeI);
 831
 832             if (regionI == -2)
 833             {
 834                 // inbetween collapsed faces. nothing needs to be done.
 835             }
 836             else if (regionI == -1)
 837             {
 838                 // edge inbetween uncollapsed faces. Handle these later on.
 839             }
 840             else
 841             {
 842                 // some faces around edge are collapsed.
 843
 844                 // Find additional set of points on edge to be used to split
 845                 // the remaining faces.
 846
 847                 labelList splitVerts;
 848                 scalarField splitWeights;
 849                 getSplitVerts
 850                 (
 851                     surf,
 852                     regionI,
 853                     spanPoints[regionI],
 854                     orderedVertices[regionI],
 855                     orderedWeights[regionI],
 856                     edgeI,
 857
 858                     splitVerts,
 859                     splitWeights
 860                 );
 861
 862                 if (splitVerts.size())
 863                 {
 864                     // Split edge using splitVerts. All non-collapsed triangles
 865                     // using edge will get split.
 866
 867
 868                     {
 869                         const pointField& localPoints = surf.localPoints();
 870                         Pout<< "edge " << edgeI << ' ' << e
 871                             << "  points "
 872                             << localPoints[e[0]] << ' ' << localPoints[e[1]]
 873                             << " split into edges with extra points:"
 874                             << endl;
 875                         forAll(splitVerts, i)
 876                         {
 877                             Pout<< "    " << splitVerts[i] << " weight "
 878                                 << splitWeights[i] << nl;
 879                         }
 880                     }
 881
 882                     const labelList& eFaces = surf.edgeFaces()[edgeI];
 883
 884                     forAll(eFaces, i)
 885                     {
 886                         label faceI = eFaces[i];
 887
 888                         if (!faceHandled[faceI] && faceToEdge[faceI] == -1)
 889                         {
 890                             // Split face to use vertices.
 891                             splitTri
 892                             (
 893                                 localFaces[faceI],
 894                                 e,
 895                                 splitVerts,
 896                                 newTris
 897                             );
 898
 899                             faceHandled[faceI] = true;
 900
 901                             nSplit++;
 902                         }
 903                     }
 904                 }
 905             }
 906         }
 907
 908         // Copy all unsplit faces
 909         forAll(faceHandled, faceI)
 910         {
 911             if (!faceHandled[faceI] && faceToEdge[faceI] == -1)
 912             {
 913                 newTris.append(localFaces[faceI]);
 914             }
 915         }
 916
 917         Pout<< "collapseBase : splitting " << nSplit << " triangles"
 918             << endl;
 919
 920         nTotalSplit += nSplit;
 921
 922         if (nSplit == 0)
 923         {
 924             break;
 925         }
 926
 927         // Pack the triangles
 928         newTris.shrink();
 929
 930         Pout<< "Resetting surface from " << surf.size() << " to "
 931             << newTris.size() << " triangles" << endl;
 932         surf = triSurface(newTris, surf.patches(), surf.localPoints());
 933
 934         {
 935             fileName fName("bla" + name(iter) + ".obj");
 936             Pout<< "Writing surf to " << fName << endl;
 937             surf.write(fName);
 938         }
 939
 940         iter++;
 941     }
 942
 943     // Remove any unused vertices
 944     surf = triSurface(surf.localFaces(), surf.patches(), surf.localPoints());
 945
 946     return nTotalSplit;
 947 }
 948
 949
 950 // ************************ vim: set sw=4 sts=4 et: ************************ //