applications/utilities/surface/surfaceClean/collapseBase.C

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