src/sampling/sampledSurface/isoSurface/isoSurfaceCell.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "isoSurfaceCell.H"
  28 #include "dictionary.H"
  29 #include "polyMesh.H"
  30 #include "mergePoints.H"
  31 #include "tetMatcher.H"
  32 #include "syncTools.H"
  33 #include "addToRunTimeSelectionTable.H"
  34 #include "faceTriangulation.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(isoSurfaceCell, 0);
  41 }
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 Foam::scalar Foam::isoSurfaceCell::isoFraction
  46 (
  47     const scalar s0,
  48     const scalar s1
  49 ) const
  50 {
  51     scalar d = s1-s0;
  52
  53     if (mag(d) > VSMALL)
  54     {
  55         return (iso_-s0)/d;
  56     }
  57     else
  58     {
  59         return -1.0;
  60     }
  61 }
  62
  63
  64 //Foam::List<Foam::triFace> Foam::isoSurfaceCell::triangulate(const face& f)
  65 // const
  66 //{
  67 //    faceList triFaces(f.nTriangles(mesh_.points()));
  68 //    label triFaceI = 0;
  69 //    f.triangles(mesh_.points(), triFaceI, triFaces);
  70 //
  71 //    List<triFace> tris(triFaces.size());
  72 //    forAll(triFaces, i)
  73 //    {
  74 //        tris[i][0] = triFaces[i][0];
  75 //        tris[i][1] = triFaces[i][1];
  76 //        tris[i][2] = triFaces[i][2];
  77 //    }
  78 //    return tris;
  79 //}
  80
  81
  82 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
  83 (
  84     const PackedBoolList& isTet,
  85     const scalarField& cellValues,
  86     const scalarField& pointValues,
  87     const label cellI
  88 ) const
  89 {
  90     const cell& cFaces = mesh_.cells()[cellI];
  91
  92     if (isTet.get(cellI) == 1)
  93     {
  94         forAll(cFaces, cFaceI)
  95         {
  96             const face& f = mesh_.faces()[cFaces[cFaceI]];
  97
  98             for (label fp = 1; fp < f.size() - 1; fp++)
  99             {
 100                 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
 101
 102                 bool aLower = (pointValues[tri[0]] < iso_);
 103                 bool bLower = (pointValues[tri[1]] < iso_);
 104                 bool cLower = (pointValues[tri[2]] < iso_);
 105
 106                 if (aLower == bLower && aLower == cLower)
 107                 {}
 108                 else
 109                 {
 110                     return CUT;
 111                 }
 112             }
 113         }
 114         return NOTCUT;
 115     }
 116     else
 117     {
 118         bool cellLower = (cellValues[cellI] < iso_);
 119
 120         // First check if there is any cut in cell
 121         bool edgeCut = false;
 122
 123         forAll(cFaces, cFaceI)
 124         {
 125             const face& f = mesh_.faces()[cFaces[cFaceI]];
 126
 127             // Check pyramids cut
 128             forAll(f, fp)
 129             {
 130                 if ((pointValues[f[fp]] < iso_) != cellLower)
 131                 {
 132                     edgeCut = true;
 133                     break;
 134                 }
 135             }
 136
 137             if (edgeCut)
 138             {
 139                 break;
 140             }
 141
 142             for (label fp = 1; fp < f.size() - 1; fp++)
 143             {
 144                 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
 145             //List<triFace> tris(triangulate(f));
 146             //forAll(tris, i)
 147             //{
 148             //    const triFace& tri = tris[i];
 149
 150                 bool aLower = (pointValues[tri[0]] < iso_);
 151                 bool bLower = (pointValues[tri[1]] < iso_);
 152                 bool cLower = (pointValues[tri[2]] < iso_);
 153
 154                 if (aLower == bLower && aLower == cLower)
 155                 {}
 156                 else
 157                 {
 158                     edgeCut = true;
 159                     break;
 160                 }
 161             }
 162
 163             if (edgeCut)
 164             {
 165                 break;
 166             }
 167         }
 168
 169         if (edgeCut)
 170         {
 171             // Count actual cuts (expensive since addressing needed)
 172             // Note: not needed if you don't want to preserve maxima/minima
 173             // centred around cellcentre. In that case just always return CUT
 174
 175             const labelList& cPoints = mesh_.cellPoints(cellI);
 176
 177             label nPyrCuts = 0;
 178
 179             forAll(cPoints, i)
 180             {
 181                 if ((pointValues[cPoints[i]] < iso_) != cellLower)
 182                 {
 183                     nPyrCuts++;
 184                 }
 185             }
 186
 187             if (nPyrCuts == cPoints.size())
 188             {
 189                 return SPHERE;
 190             }
 191             else
 192             {
 193                 return CUT;
 194             }
 195         }
 196         else
 197         {
 198             return NOTCUT;
 199         }
 200     }
 201 }
 202
 203
 204 void Foam::isoSurfaceCell::calcCutTypes
 205 (
 206     const PackedBoolList& isTet,
 207     const scalarField& cVals,
 208     const scalarField& pVals
 209 )
 210 {
 211     cellCutType_.setSize(mesh_.nCells());
 212     nCutCells_ = 0;
 213     forAll(mesh_.cells(), cellI)
 214     {
 215         cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
 216
 217         if (cellCutType_[cellI] == CUT)
 218         {
 219             nCutCells_++;
 220         }
 221     }
 222
 223     if (debug)
 224     {
 225         Pout<< "isoSurfaceCell : detected " << nCutCells_
 226             << " candidate cut cells." << endl;
 227     }
 228 }
 229
 230
 231
 232 // Return the two common points between two triangles
 233 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
 234 (
 235     const labelledTri& tri0,
 236     const labelledTri& tri1
 237 )
 238 {
 239     labelPair common(-1, -1);
 240
 241     label fp0 = 0;
 242     label fp1 = findIndex(tri1, tri0[fp0]);
 243
 244     if (fp1 == -1)
 245     {
 246         fp0 = 1;
 247         fp1 = findIndex(tri1, tri0[fp0]);
 248     }
 249
 250     if (fp1 != -1)
 251     {
 252         // So tri0[fp0] is tri1[fp1]
 253
 254         // Find next common point
 255         label fp0p1 = tri0.fcIndex(fp0);
 256         label fp1p1 = tri1.fcIndex(fp1);
 257         label fp1m1 = tri1.rcIndex(fp1);
 258
 259         if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
 260         {
 261             common[0] = tri0[fp0];
 262             common[1] = tri0[fp0p1];
 263         }
 264     }
 265     return common;
 266 }
 267
 268
 269 // Caculate centre of surface.
 270 Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
 271 {
 272     vector sum = vector::zero;
 273
 274     forAll(s, i)
 275     {
 276         sum += s[i].centre(s.points());
 277     }
 278     return sum/s.size();
 279 }
 280
 281
 282 // Replace surface (localPoints, localTris) with single point. Returns
 283 // point. Destructs arguments.
 284 Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
 285 (
 286     pointField& localPoints,
 287     DynamicList<labelledTri, 64>& localTris
 288 )
 289 {
 290     pointIndexHit info(false, vector::zero, localTris.size());
 291
 292     if (localTris.size() == 1)
 293     {
 294         const labelledTri& tri = localTris[0];
 295         info.setPoint(tri.centre(localPoints));
 296         info.setHit();
 297     }
 298     else if (localTris.size() == 2)
 299     {
 300         // Check if the two triangles share an edge.
 301         const labelledTri& tri0 = localTris[0];
 302         const labelledTri& tri1 = localTris[0];
 303
 304         labelPair shared = findCommonPoints(tri0, tri1);
 305
 306         if (shared[0] != -1)
 307         {
 308             info.setPoint
 309             (
 310                 0.5
 311               * (
 312                     tri0.centre(localPoints)
 313                   + tri1.centre(localPoints)
 314                 )
 315             );
 316             info.setHit();
 317         }
 318     }
 319     else if (localTris.size())
 320     {
 321         // Check if single region. Rare situation.
 322         triSurface surf
 323         (
 324             localTris,
 325             geometricSurfacePatchList(0),
 326             localPoints,
 327             true
 328         );
 329         localTris.clearStorage();
 330
 331         labelList faceZone;
 332         label nZones = surf.markZones
 333         (
 334             boolList(surf.nEdges(), false),
 335             faceZone
 336         );
 337
 338         if (nZones == 1)
 339         {
 340             info.setPoint(calcCentre(surf));
 341             info.setHit();
 342         }
 343     }
 344
 345     return info;
 346 }
 347
 348
 349 void Foam::isoSurfaceCell::calcSnappedCc
 350 (
 351     const PackedBoolList& isTet,
 352     const scalarField& cVals,
 353     const scalarField& pVals,
 354
 355     DynamicList<point>& snappedPoints,
 356     labelList& snappedCc
 357 ) const
 358 {
 359     const pointField& cc = mesh_.cellCentres();
 360     const pointField& pts = mesh_.points();
 361
 362     snappedCc.setSize(mesh_.nCells());
 363     snappedCc = -1;
 364
 365     // Work arrays
 366     DynamicList<point, 64> localPoints(64);
 367     DynamicList<labelledTri, 64> localTris(64);
 368     Map<label> pointToLocal(64);
 369
 370     forAll(mesh_.cells(), cellI)
 371     {
 372         if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
 373         {
 374             scalar cVal = cVals[cellI];
 375
 376             const cell& cFaces = mesh_.cells()[cellI];
 377
 378             localPoints.clear();
 379             localTris.clear();
 380             pointToLocal.clear();
 381
 382             // Create points for all intersections close to cell centre
 383             // (i.e. from pyramid edges)
 384
 385             forAll(cFaces, cFaceI)
 386             {
 387                 const face& f = mesh_.faces()[cFaces[cFaceI]];
 388
 389                 forAll(f, fp)
 390                 {
 391                     label pointI = f[fp];
 392
 393                     scalar s = isoFraction(cVal, pVals[pointI]);
 394
 395                     if (s >= 0.0 && s <= 0.5)
 396                     {
 397                         if (pointToLocal.insert(pointI, localPoints.size()))
 398                         {
 399                             localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
 400                         }
 401                     }
 402                 }
 403             }
 404
 405             if (localPoints.size() == 1)
 406             {
 407                 // No need for any analysis.
 408                 snappedCc[cellI] = snappedPoints.size();
 409                 snappedPoints.append(localPoints[0]);
 410
 411                 //Pout<< "cell:" << cellI
 412                 //    << " at " << mesh_.cellCentres()[cellI]
 413                 //    << " collapsing " << localPoints
 414                 //    << " intersections down to "
 415                 //    << snappedPoints[snappedCc[cellI]] << endl;
 416             }
 417             else if (localPoints.size() == 2)
 418             {
 419                 //? No need for any analysis.???
 420                 snappedCc[cellI] = snappedPoints.size();
 421                 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
 422
 423                 //Pout<< "cell:" << cellI
 424                 //    << " at " << mesh_.cellCentres()[cellI]
 425                 //    << " collapsing " << localPoints
 426                 //    << " intersections down to "
 427                 //    << snappedPoints[snappedCc[cellI]] << endl;
 428             }
 429             else if (localPoints.size())
 430             {
 431                 // Need to analyse
 432                 forAll(cFaces, cFaceI)
 433                 {
 434                     const face& f = mesh_.faces()[cFaces[cFaceI]];
 435
 436                     // Do a tetrahedrisation. Each face to cc becomes pyr.
 437                     // Each pyr gets split into tets by diagonalisation
 438                     // of face.
 439
 440                     for (label fp = 1; fp < f.size() - 1; fp++)
 441                     {
 442                         triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
 443                     //List<triFace> tris(triangulate(f));
 444                     //forAll(tris, i)
 445                     //{
 446                     //    const triFace& tri = tris[i];
 447
 448                         // Get fractions for the three edges to cell centre
 449                         FixedList<scalar, 3> s(3);
 450                         s[0] = isoFraction(cVal, pVals[tri[0]]);
 451                         s[1] = isoFraction(cVal, pVals[tri[1]]);
 452                         s[2] = isoFraction(cVal, pVals[tri[2]]);
 453
 454                         if
 455                         (
 456                             (s[0] >= 0.0 && s[0] <= 0.5)
 457                          && (s[1] >= 0.0 && s[1] <= 0.5)
 458                          && (s[2] >= 0.0 && s[2] <= 0.5)
 459                         )
 460                         {
 461                             localTris.append
 462                             (
 463                                 labelledTri
 464                                 (
 465                                     pointToLocal[tri[0]],
 466                                     pointToLocal[tri[1]],
 467                                     pointToLocal[tri[2]],
 468                                     0
 469                                 )
 470                             );
 471                         }
 472                     }
 473                 }
 474
 475                 pointField surfPoints;
 476                 surfPoints.transfer(localPoints);
 477                 pointIndexHit info = collapseSurface(surfPoints, localTris);
 478
 479                 if (info.hit())
 480                 {
 481                     snappedCc[cellI] = snappedPoints.size();
 482                     snappedPoints.append(info.hitPoint());
 483
 484                     //Pout<< "cell:" << cellI
 485                     //    << " at " << mesh_.cellCentres()[cellI]
 486                     //    << " collapsing " << surfPoints
 487                     //    << " intersections down to "
 488                     //    << snappedPoints[snappedCc[cellI]] << endl;
 489                 }
 490             }
 491         }
 492     }
 493 }
 494
 495
 496 // Generate triangles for face connected to pointI
 497 void Foam::isoSurfaceCell::genPointTris
 498 (
 499     const scalarField& cellValues,
 500     const scalarField& pointValues,
 501     const label pointI,
 502     const face& f,
 503     const label cellI,
 504     DynamicList<point, 64>& localTriPoints
 505 ) const
 506 {
 507     const pointField& cc = mesh_.cellCentres();
 508     const pointField& pts = mesh_.points();
 509
 510     for (label fp = 1; fp < f.size() - 1; fp++)
 511     {
 512         triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
 513     //List<triFace> tris(triangulate(f));
 514     //forAll(tris, i)
 515     //{
 516     //    const triFace& tri = tris[i];
 517
 518         label index = findIndex(tri, pointI);
 519
 520         if (index == -1)
 521         {
 522             continue;
 523         }
 524
 525         // Tet between index..index-1, index..index+1, index..cc
 526         label b = tri[tri.fcIndex(index)];
 527         label c = tri[tri.rcIndex(index)];
 528
 529         // Get fractions for the three edges emanating from point
 530         FixedList<scalar, 3> s(3);
 531         s[0] = isoFraction(pointValues[pointI], pointValues[b]);
 532         s[1] = isoFraction(pointValues[pointI], pointValues[c]);
 533         s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
 534
 535         if
 536         (
 537             (s[0] >= 0.0 && s[0] <= 0.5)
 538          && (s[1] >= 0.0 && s[1] <= 0.5)
 539          && (s[2] >= 0.0 && s[2] <= 0.5)
 540         )
 541         {
 542             localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
 543             localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
 544             localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
 545         }
 546     }
 547 }
 548
 549
 550 // Generate triangle for tet connected to pointI
 551 void Foam::isoSurfaceCell::genPointTris
 552 (
 553     const scalarField& pointValues,
 554     const label pointI,
 555     const label cellI,
 556     DynamicList<point, 64>& localTriPoints
 557 ) const
 558 {
 559     const pointField& pts = mesh_.points();
 560     const cell& cFaces = mesh_.cells()[cellI];
 561
 562     FixedList<label, 4> tet;
 563
 564     label face0 = cFaces[0];
 565     const face& f0 = mesh_.faces()[face0];
 566
 567     if (mesh_.faceOwner()[face0] != cellI)
 568     {
 569         tet[0] = f0[0];
 570         tet[1] = f0[1];
 571         tet[2] = f0[2];
 572     }
 573     else
 574     {
 575         tet[0] = f0[0];
 576         tet[1] = f0[2];
 577         tet[2] = f0[1];
 578     }
 579
 580     // Find the point on the next face that is not on f0
 581     const face& f1 = mesh_.faces()[cFaces[1]];
 582
 583     forAll(f1, fp)
 584     {
 585         label p1 = f1[fp];
 586
 587         if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
 588         {
 589             tet[3] = p1;
 590             break;
 591         }
 592     }
 593
 594
 595     // Get the index of pointI
 596
 597     label i0 = findIndex(tet, pointI);
 598     label i1 = tet.fcIndex(i0);
 599     label i2 = tet.fcIndex(i1);
 600     label i3 = tet.fcIndex(i2);
 601
 602     // Get fractions for the three edges emanating from point
 603     FixedList<scalar, 3> s(3);
 604     s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
 605     s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
 606     s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
 607
 608     if
 609     (
 610         (s[0] >= 0.0 && s[0] <= 0.5)
 611      && (s[1] >= 0.0 && s[1] <= 0.5)
 612      && (s[2] >= 0.0 && s[2] <= 0.5)
 613     )
 614     {
 615         localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
 616         localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
 617         localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
 618     }
 619 }
 620
 621
 622 void Foam::isoSurfaceCell::calcSnappedPoint
 623 (
 624     const PackedBoolList& isBoundaryPoint,
 625     const PackedBoolList& isTet,
 626     const scalarField& cVals,
 627     const scalarField& pVals,
 628
 629     DynamicList<point>& snappedPoints,
 630     labelList& snappedPoint
 631 ) const
 632 {
 633     const point greatPoint(VGREAT, VGREAT, VGREAT);
 634     pointField collapsedPoint(mesh_.nPoints(), greatPoint);
 635
 636
 637     // Work arrays
 638     DynamicList<point, 64> localTriPoints(100);
 639     labelHashSet localPointCells(100);
 640
 641     forAll(mesh_.pointFaces(), pointI)
 642     {
 643         if (isBoundaryPoint.get(pointI) == 1)
 644         {
 645             continue;
 646         }
 647
 648         const labelList& pFaces = mesh_.pointFaces()[pointI];
 649
 650         bool anyCut = false;
 651
 652         forAll(pFaces, i)
 653         {
 654             label faceI = pFaces[i];
 655
 656             if
 657             (
 658                 cellCutType_[mesh_.faceOwner()[faceI]] == CUT
 659              || (
 660                     mesh_.isInternalFace(faceI)
 661                  && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
 662                 )
 663             )
 664             {
 665                 anyCut = true;
 666                 break;
 667             }
 668         }
 669
 670         if (!anyCut)
 671         {
 672             continue;
 673         }
 674
 675
 676         localPointCells.clear();
 677         localTriPoints.clear();
 678
 679         forAll(pFaces, pFaceI)
 680         {
 681             label faceI = pFaces[pFaceI];
 682             const face& f = mesh_.faces()[faceI];
 683             label own = mesh_.faceOwner()[faceI];
 684
 685             // Triangulate around f[0] on owner side
 686             if (isTet.get(own) == 1)
 687             {
 688                 if (localPointCells.insert(own))
 689                 {
 690                     genPointTris(pVals, pointI, own, localTriPoints);
 691                 }
 692             }
 693             else
 694             {
 695                 genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
 696             }
 697
 698             if (mesh_.isInternalFace(faceI))
 699             {
 700                 label nei = mesh_.faceNeighbour()[faceI];
 701
 702                 if (isTet.get(nei) == 1)
 703                 {
 704                     if (localPointCells.insert(nei))
 705                     {
 706                         genPointTris(pVals, pointI, nei, localTriPoints);
 707                     }
 708                 }
 709                 else
 710                 {
 711                     genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
 712                 }
 713             }
 714         }
 715
 716         if (localTriPoints.size() == 3)
 717         {
 718             // Single triangle. No need for any analysis. Average points.
 719             pointField points;
 720             points.transfer(localTriPoints);
 721             collapsedPoint[pointI] = sum(points)/points.size();
 722
 723             //Pout<< "    point:" << pointI
 724             //    << " replacing coord:" << mesh_.points()[pointI]
 725             //    << " by average:" << collapsedPoint[pointI] << endl;
 726         }
 727         else if (localTriPoints.size())
 728         {
 729             // Convert points into triSurface.
 730
 731             // Merge points and compact out non-valid triangles
 732             labelList triMap;               // merged to unmerged triangle
 733             labelList triPointReverseMap;   // unmerged to merged point
 734             triSurface surf
 735             (
 736                 stitchTriPoints
 737                 (
 738                     false,                  // do not check for duplicate tris
 739                     localTriPoints,
 740                     triPointReverseMap,
 741                     triMap
 742                 )
 743             );
 744
 745             labelList faceZone;
 746             label nZones = surf.markZones
 747             (
 748                 boolList(surf.nEdges(), false),
 749                 faceZone
 750             );
 751
 752             if (nZones == 1)
 753             {
 754                 collapsedPoint[pointI] = calcCentre(surf);
 755                 //Pout<< "    point:" << pointI << " nZones:" << nZones
 756                 //    << " replacing coord:" << mesh_.points()[pointI]
 757                 //    << " by average:" << collapsedPoint[pointI] << endl;
 758             }
 759         }
 760     }
 761
 762     syncTools::syncPointList
 763     (
 764         mesh_,
 765         collapsedPoint,
 766         minEqOp<point>(),
 767         greatPoint,
 768         true                // are coordinates so separate
 769     );
 770
 771     snappedPoint.setSize(mesh_.nPoints());
 772     snappedPoint = -1;
 773
 774     forAll(collapsedPoint, pointI)
 775     {
 776         if (collapsedPoint[pointI] != greatPoint)
 777         {
 778             snappedPoint[pointI] = snappedPoints.size();
 779             snappedPoints.append(collapsedPoint[pointI]);
 780         }
 781     }
 782 }
 783
 784
 785
 786
 787 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
 788 (
 789     const bool checkDuplicates,
 790     const List<point>& triPoints,
 791     labelList& triPointReverseMap,  // unmerged to merged point
 792     labelList& triMap               // merged to unmerged triangle
 793 ) const
 794 {
 795     label nTris = triPoints.size()/3;
 796
 797     if ((triPoints.size() % 3) != 0)
 798     {
 799         FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
 800             << "Problem: number of points " << triPoints.size()
 801             << " not a multiple of 3." << abort(FatalError);
 802     }
 803
 804     pointField newPoints;
 805     mergePoints
 806     (
 807         triPoints,
 808         mergeDistance_,
 809         false,
 810         triPointReverseMap,
 811         newPoints
 812     );
 813
 814     // Check that enough merged.
 815     if (debug)
 816     {
 817         pointField newNewPoints;
 818         labelList oldToNew;
 819         bool hasMerged = mergePoints
 820         (
 821             newPoints,
 822             mergeDistance_,
 823             true,
 824             oldToNew,
 825             newNewPoints
 826         );
 827
 828         if (hasMerged)
 829         {
 830             FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
 831                 << "Merged points contain duplicates"
 832                 << " when merging with distance " << mergeDistance_ << endl
 833                 << "merged:" << newPoints.size() << " re-merged:"
 834                 << newNewPoints.size()
 835                 << abort(FatalError);
 836         }
 837     }
 838
 839
 840     List<labelledTri> tris;
 841     {
 842         DynamicList<labelledTri> dynTris(nTris);
 843         label rawPointI = 0;
 844         DynamicList<label> newToOldTri(nTris);
 845
 846         for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
 847         {
 848             labelledTri tri
 849             (
 850                 triPointReverseMap[rawPointI],
 851                 triPointReverseMap[rawPointI+1],
 852                 triPointReverseMap[rawPointI+2],
 853                 0
 854             );
 855             rawPointI += 3;
 856
 857             if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
 858             {
 859                 newToOldTri.append(oldTriI);
 860                 dynTris.append(tri);
 861             }
 862         }
 863
 864         triMap.transfer(newToOldTri);
 865         tris.transfer(dynTris);
 866     }
 867
 868
 869     // Use face centres to determine 'flat hole' situation (see RMT paper).
 870     // Two unconnected triangles get connected because (some of) the edges
 871     // separating them get collapsed. Below only checks for duplicate triangles,
 872     // not non-manifold edge connectivity.
 873     if (checkDuplicates)
 874     {
 875         if (debug)
 876         {
 877             Pout<< "isoSurfaceCell : merged from " << nTris
 878                 << " down to " << tris.size() << " triangles." << endl;
 879         }
 880
 881         pointField centres(tris.size());
 882         forAll(tris, triI)
 883         {
 884             centres[triI] = tris[triI].centre(newPoints);
 885         }
 886
 887         pointField mergedCentres;
 888         labelList oldToMerged;
 889         bool hasMerged = mergePoints
 890         (
 891             centres,
 892             mergeDistance_,
 893             false,
 894             oldToMerged,
 895             mergedCentres
 896         );
 897
 898         if (debug)
 899         {
 900             Pout<< "isoSurfaceCell : detected "
 901                 << centres.size()-mergedCentres.size()
 902                 << " duplicate triangles." << endl;
 903         }
 904
 905         if (hasMerged)
 906         {
 907             // Filter out duplicates.
 908             label newTriI = 0;
 909             DynamicList<label> newToOldTri(tris.size());
 910             labelList newToMaster(mergedCentres.size(), -1);
 911             forAll(tris, triI)
 912             {
 913                 label mergedI = oldToMerged[triI];
 914
 915                 if (newToMaster[mergedI] == -1)
 916                 {
 917                     newToMaster[mergedI] = triI;
 918                     newToOldTri.append(triMap[triI]);
 919                     tris[newTriI++] = tris[triI];
 920                 }
 921             }
 922
 923             triMap.transfer(newToOldTri);
 924             tris.setSize(newTriI);
 925         }
 926     }
 927
 928     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
 929 }
 930
 931
 932 // Does face use valid vertices?
 933 bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
 934 {
 935     // Simple check on indices ok.
 936
 937     const labelledTri& f = surf[faceI];
 938
 939     if
 940     (
 941         (f[0] < 0) || (f[0] >= surf.points().size())
 942      || (f[1] < 0) || (f[1] >= surf.points().size())
 943      || (f[2] < 0) || (f[2] >= surf.points().size())
 944     )
 945     {
 946         WarningIn("validTri(const triSurface&, const label)")
 947             << "triangle " << faceI << " vertices " << f
 948             << " uses point indices outside point range 0.."
 949             << surf.points().size()-1 << endl;
 950
 951         return false;
 952     }
 953
 954     if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
 955     {
 956         WarningIn("validTri(const triSurface&, const label)")
 957             << "triangle " << faceI
 958             << " uses non-unique vertices " << f
 959             << endl;
 960         return false;
 961     }
 962
 963     // duplicate triangle check
 964
 965     const labelList& fFaces = surf.faceFaces()[faceI];
 966
 967     // Check if faceNeighbours use same points as this face.
 968     // Note: discards normal information - sides of baffle are merged.
 969     forAll(fFaces, i)
 970     {
 971         label nbrFaceI = fFaces[i];
 972
 973         if (nbrFaceI <= faceI)
 974         {
 975             // lower numbered faces already checked
 976             continue;
 977         }
 978
 979         const labelledTri& nbrF = surf[nbrFaceI];
 980
 981         if
 982         (
 983             ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
 984          && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
 985          && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
 986         )
 987         {
 988             WarningIn("validTri(const triSurface&, const label)")
 989                 << "triangle " << faceI << " vertices " << f
 990                 << " has the same vertices as triangle " << nbrFaceI
 991                 << " vertices " << nbrF
 992                 << endl;
 993
 994             return false;
 995         }
 996     }
 997     return true;
 998 }
 999
1000
1001 void Foam::isoSurfaceCell::calcAddressing
1002 (
1003     const triSurface& surf,
1004     List<FixedList<label, 3> >& faceEdges,
1005     labelList& edgeFace0,
1006     labelList& edgeFace1,
1007     Map<labelList>& edgeFacesRest
1008 ) const
1009 {
1010     const pointField& points = surf.points();
1011
1012     pointField edgeCentres(3*surf.size());
1013     label edgeI = 0;
1014     forAll(surf, triI)
1015     {
1016         const labelledTri& tri = surf[triI];
1017         edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1018         edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1019         edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1020     }
1021
1022     pointField mergedCentres;
1023     labelList oldToMerged;
1024     bool hasMerged = mergePoints
1025     (
1026         edgeCentres,
1027         mergeDistance_,
1028         false,
1029         oldToMerged,
1030         mergedCentres
1031     );
1032
1033     if (debug)
1034     {
1035         Pout<< "isoSurfaceCell : detected "
1036             << mergedCentres.size()
1037             << " edges on " << surf.size() << " triangles." << endl;
1038     }
1039
1040     if (!hasMerged)
1041     {
1042         if (surf.size() == 1)
1043         {
1044             faceEdges.setSize(1);
1045             faceEdges[0][0] = 0;
1046             faceEdges[0][1] = 1;
1047             faceEdges[0][2] = 2;
1048             edgeFace0.setSize(1);
1049             edgeFace0[0] = 0;
1050             edgeFace1.setSize(1);
1051             edgeFace1[0] = -1;
1052             edgeFacesRest.clear();
1053         }
1054         return;
1055     }
1056
1057
1058     // Determine faceEdges
1059     faceEdges.setSize(surf.size());
1060     edgeI = 0;
1061     forAll(surf, triI)
1062     {
1063         faceEdges[triI][0] = oldToMerged[edgeI++];
1064         faceEdges[triI][1] = oldToMerged[edgeI++];
1065         faceEdges[triI][2] = oldToMerged[edgeI++];
1066     }
1067
1068
1069     // Determine edgeFaces
1070     edgeFace0.setSize(mergedCentres.size());
1071     edgeFace0 = -1;
1072     edgeFace1.setSize(mergedCentres.size());
1073     edgeFace1 = -1;
1074     edgeFacesRest.clear();
1075
1076     forAll(oldToMerged, oldEdgeI)
1077     {
1078         label triI = oldEdgeI / 3;
1079         label edgeI = oldToMerged[oldEdgeI];
1080
1081         if (edgeFace0[edgeI] == -1)
1082         {
1083             edgeFace0[edgeI] = triI;
1084         }
1085         else if (edgeFace1[edgeI] == -1)
1086         {
1087             edgeFace1[edgeI] = triI;
1088         }
1089         else
1090         {
1091             //WarningIn("orientSurface(triSurface&)")
1092             //    << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
1093             //    << " used by more than two triangles: " << edgeFace0[edgeI]
1094             //    << ", "
1095             //    << edgeFace1[edgeI] << " and " << triI << endl;
1096             Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
1097
1098             if (iter != edgeFacesRest.end())
1099             {
1100                 labelList& eFaces = iter();
1101                 label sz = eFaces.size();
1102                 eFaces.setSize(sz+1);
1103                 eFaces[sz] = triI;
1104             }
1105             else
1106             {
1107                 edgeFacesRest.insert(edgeI, labelList(1, triI));
1108             }
1109         }
1110     }
1111 }
1112
1113
1114 void Foam::isoSurfaceCell::walkOrientation
1115 (
1116     const triSurface& surf,
1117     const List<FixedList<label, 3> >& faceEdges,
1118     const labelList& edgeFace0,
1119     const labelList& edgeFace1,
1120     const label seedTriI,
1121     labelList& flipState
1122 )
1123 {
1124     // Do walk for consistent orientation.
1125     DynamicList<label> changedFaces(surf.size());
1126
1127     changedFaces.append(seedTriI);
1128
1129     while (changedFaces.size())
1130     {
1131         DynamicList<label> newChangedFaces(changedFaces.size());
1132
1133         forAll(changedFaces, i)
1134         {
1135             label triI = changedFaces[i];
1136             const labelledTri& tri = surf[triI];
1137             const FixedList<label, 3>& fEdges = faceEdges[triI];
1138
1139             forAll(fEdges, fp)
1140             {
1141                 label edgeI = fEdges[fp];
1142
1143                 // my points:
1144                 label p0 = tri[fp];
1145                 label p1 = tri[tri.fcIndex(fp)];
1146
1147                 label nbrI =
1148                 (
1149                     edgeFace0[edgeI] != triI
1150                   ? edgeFace0[edgeI]
1151                   : edgeFace1[edgeI]
1152                 );
1153
1154                 if (nbrI != -1 && flipState[nbrI] == -1)
1155                 {
1156                     const labelledTri& nbrTri = surf[nbrI];
1157
1158                     // nbr points
1159                     label nbrFp = findIndex(nbrTri, p0);
1160                     label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1161
1162                     bool sameOrientation = (p1 == nbrP1);
1163
1164                     if (flipState[triI] == 0)
1165                     {
1166                         flipState[nbrI] = (sameOrientation ? 0 : 1);
1167                     }
1168                     else
1169                     {
1170                         flipState[nbrI] = (sameOrientation ? 1 : 0);
1171                     }
1172                     newChangedFaces.append(nbrI);
1173                 }
1174             }
1175         }
1176
1177         changedFaces.transfer(newChangedFaces);
1178     }
1179 }
1180
1181
1182 void Foam::isoSurfaceCell::orientSurface
1183 (
1184     triSurface& surf,
1185     const List<FixedList<label, 3> >& faceEdges,
1186     const labelList& edgeFace0,
1187     const labelList& edgeFace1,
1188     const Map<labelList>& edgeFacesRest
1189 )
1190 {
1191     // -1 : unvisited
1192     //  0 : leave as is
1193     //  1 : flip
1194     labelList flipState(surf.size(), -1);
1195
1196     label seedTriI = 0;
1197
1198     while (true)
1199     {
1200         // Find first unvisited triangle
1201         for
1202         (
1203             ;
1204             seedTriI < surf.size() && flipState[seedTriI] != -1;
1205             seedTriI++
1206         )
1207         {}
1208
1209         if (seedTriI == surf.size())
1210         {
1211             break;
1212         }
1213
1214         // Note: Determine orientation of seedTriI?
1215         // for now assume it is ok
1216         flipState[seedTriI] = 0;
1217
1218         walkOrientation
1219         (
1220             surf,
1221             faceEdges,
1222             edgeFace0,
1223             edgeFace1,
1224             seedTriI,
1225             flipState
1226         );
1227     }
1228
1229     // Do actual flipping
1230     surf.clearOut();
1231     forAll(surf, triI)
1232     {
1233         if (flipState[triI] == 1)
1234         {
1235             labelledTri tri(surf[triI]);
1236
1237             surf[triI][0] = tri[0];
1238             surf[triI][1] = tri[2];
1239             surf[triI][2] = tri[1];
1240         }
1241         else if (flipState[triI] == -1)
1242         {
1243             FatalErrorIn
1244             (
1245                 "isoSurfaceCell::orientSurface(triSurface&, const label)"
1246             )   << "problem" << abort(FatalError);
1247         }
1248     }
1249 }
1250
1251
1252 // Checks if triangle is connected through edgeI only.
1253 bool Foam::isoSurfaceCell::danglingTriangle
1254 (
1255     const FixedList<label, 3>& fEdges,
1256     const labelList& edgeFace1
1257 )
1258 {
1259     label nOpen = 0;
1260     forAll(fEdges, i)
1261     {
1262         if (edgeFace1[fEdges[i]] == -1)
1263         {
1264             nOpen++;
1265         }
1266     }
1267
1268     if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1269     {
1270         return true;
1271     }
1272     else
1273     {
1274         return false;
1275     }
1276 }
1277
1278
1279 // Mark triangles to keep. Returns number of dangling triangles.
1280 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1281 (
1282     const List<FixedList<label, 3> >& faceEdges,
1283     const labelList& edgeFace0,
1284     const labelList& edgeFace1,
1285     const Map<labelList>& edgeFacesRest,
1286     boolList& keepTriangles
1287 )
1288 {
1289     keepTriangles.setSize(faceEdges.size());
1290     keepTriangles = true;
1291
1292     label nDangling = 0;
1293
1294     // Remove any dangling triangles
1295     forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
1296     {
1297         // These are all the non-manifold edges. Filter out all triangles
1298         // with only one connected edge (= this edge)
1299
1300         label edgeI = iter.key();
1301         const labelList& otherEdgeFaces = iter();
1302
1303         // Remove all dangling triangles
1304         if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1305         {
1306             keepTriangles[edgeFace0[edgeI]] = false;
1307             nDangling++;
1308         }
1309         if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1310         {
1311             keepTriangles[edgeFace1[edgeI]] = false;
1312             nDangling++;
1313         }
1314         forAll(otherEdgeFaces, i)
1315         {
1316             label triI = otherEdgeFaces[i];
1317             if (danglingTriangle(faceEdges[triI], edgeFace1))
1318             {
1319                 keepTriangles[triI] = false;
1320                 nDangling++;
1321             }
1322         }
1323     }
1324     return nDangling;
1325 }
1326
1327
1328 Foam::triSurface Foam::isoSurfaceCell::subsetMesh
1329 (
1330     const triSurface& s,
1331     const labelList& newToOldFaces,
1332     labelList& oldToNewPoints,
1333     labelList& newToOldPoints
1334 )
1335 {
1336     const boolList include
1337     (
1338         createWithValues<boolList>
1339         (
1340             s.size(),
1341             false,
1342             newToOldFaces,
1343             true
1344         )
1345     );
1346
1347     newToOldPoints.setSize(s.points().size());
1348     oldToNewPoints.setSize(s.points().size());
1349     oldToNewPoints = -1;
1350     {
1351         label pointI = 0;
1352
1353         forAll(include, oldFacei)
1354         {
1355             if (include[oldFacei])
1356             {
1357                 // Renumber labels for triangle
1358                 const labelledTri& tri = s[oldFacei];
1359
1360                 forAll(tri, fp)
1361                 {
1362                     label oldPointI = tri[fp];
1363
1364                     if (oldToNewPoints[oldPointI] == -1)
1365                     {
1366                         oldToNewPoints[oldPointI] = pointI;
1367                         newToOldPoints[pointI++] = oldPointI;
1368                     }
1369                 }
1370             }
1371         }
1372         newToOldPoints.setSize(pointI);
1373     }
1374
1375     // Extract points
1376     pointField newPoints(newToOldPoints.size());
1377     forAll(newToOldPoints, i)
1378     {
1379         newPoints[i] = s.points()[newToOldPoints[i]];
1380     }
1381     // Extract faces
1382     List<labelledTri> newTriangles(newToOldFaces.size());
1383
1384     forAll(newToOldFaces, i)
1385     {
1386         // Get old vertex labels
1387         const labelledTri& tri = s[newToOldFaces[i]];
1388
1389         newTriangles[i][0] = oldToNewPoints[tri[0]];
1390         newTriangles[i][1] = oldToNewPoints[tri[1]];
1391         newTriangles[i][2] = oldToNewPoints[tri[2]];
1392         newTriangles[i].region() = tri.region();
1393     }
1394
1395     // Reuse storage.
1396     return triSurface(newTriangles, s.patches(), newPoints, true);
1397 }
1398
1399
1400 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1401
1402 Foam::isoSurfaceCell::isoSurfaceCell
1403 (
1404     const polyMesh& mesh,
1405     const scalarField& cVals,
1406     const scalarField& pVals,
1407     const scalar iso,
1408     const bool regularise,
1409     const scalar mergeTol
1410 )
1411 :
1412     mesh_(mesh),
1413     iso_(iso),
1414     mergeDistance_(mergeTol*mesh.bounds().mag())
1415 {
1416     // Determine if cell is tet
1417     PackedBoolList isTet(mesh_.nCells());
1418     {
1419         tetMatcher tet;
1420
1421         forAll(isTet, cellI)
1422         {
1423             if (tet.isA(mesh_, cellI))
1424             {
1425                 isTet.set(cellI, 1);
1426             }
1427         }
1428     }
1429
1430     // Determine if point is on boundary. Points on boundaries are never
1431     // snapped. Coupled boundaries are handled explicitly so not marked here.
1432     PackedBoolList isBoundaryPoint(mesh_.nPoints());
1433     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1434     forAll(patches, patchI)
1435     {
1436         const polyPatch& pp = patches[patchI];
1437
1438         if (!pp.coupled())
1439         {
1440             label faceI = pp.start();
1441             forAll(pp, i)
1442             {
1443                 const face& f = mesh_.faces()[faceI++];
1444
1445                 forAll(f, fp)
1446                 {
1447                     isBoundaryPoint.set(f[fp], 1);
1448                 }
1449             }
1450         }
1451     }
1452
1453
1454     // Determine if any cut through cell
1455     calcCutTypes(isTet, cVals, pVals);
1456
1457     DynamicList<point> snappedPoints(nCutCells_);
1458
1459     // Per cc -1 or a point inside snappedPoints.
1460     labelList snappedCc;
1461     if (regularise)
1462     {
1463         calcSnappedCc
1464         (
1465             isTet,
1466             cVals,
1467             pVals,
1468             snappedPoints,
1469             snappedCc
1470         );
1471     }
1472     else
1473     {
1474         snappedCc.setSize(mesh_.nCells());
1475         snappedCc = -1;
1476     }
1477
1478     if (debug)
1479     {
1480         Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
1481             << " cell centres to intersection." << endl;
1482     }
1483
1484     snappedPoints.shrink();
1485     label nCellSnaps = snappedPoints.size();
1486
1487     // Per point -1 or a point inside snappedPoints.
1488     labelList snappedPoint;
1489     if (regularise)
1490     {
1491         calcSnappedPoint
1492         (
1493             isBoundaryPoint,
1494             isTet,
1495             cVals,
1496             pVals,
1497             snappedPoints,
1498             snappedPoint
1499         );
1500     }
1501     else
1502     {
1503         snappedPoint.setSize(mesh_.nPoints());
1504         snappedPoint = -1;
1505     }
1506
1507     if (debug)
1508     {
1509         Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
1510             << " vertices to intersection." << endl;
1511     }
1512
1513
1514
1515     DynamicList<point> triPoints(nCutCells_);
1516     DynamicList<label> triMeshCells(nCutCells_);
1517
1518     generateTriPoints
1519     (
1520         cVals,
1521         pVals,
1522
1523         mesh_.cellCentres(),
1524         mesh_.points(),
1525
1526         snappedPoints,
1527         snappedCc,
1528         snappedPoint,
1529
1530         triPoints,
1531         triMeshCells
1532     );
1533
1534     if (debug)
1535     {
1536         Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
1537             << " unmerged triangles." << endl;
1538     }
1539
1540     // Merge points and compact out non-valid triangles
1541     labelList triMap;           // merged to unmerged triangle
1542     triSurface::operator=
1543     (
1544         stitchTriPoints
1545         (
1546             true,               // check for duplicate tris
1547             triPoints,
1548             triPointMergeMap_,  // unmerged to merged point
1549             triMap
1550         )
1551     );
1552
1553     if (debug)
1554     {
1555         Pout<< "isoSurfaceCell : generated " << triMap.size()
1556             << " merged triangles." << endl;
1557     }
1558
1559     meshCells_.setSize(triMap.size());
1560     forAll(triMap, i)
1561     {
1562         meshCells_[i] = triMeshCells[triMap[i]];
1563     }
1564
1565     if (debug)
1566     {
1567         Pout<< "isoSurfaceCell : checking " << size()
1568             << " triangles for validity." << endl;
1569
1570         forAll(*this, triI)
1571         {
1572             // Copied from surfaceCheck
1573             validTri(*this, triI);
1574         }
1575     }
1576
1577
1578 //if (false)
1579 {
1580     List<FixedList<label, 3> > faceEdges;
1581     labelList edgeFace0, edgeFace1;
1582     Map<labelList> edgeFacesRest;
1583
1584
1585     while (true)
1586     {
1587         // Calculate addressing
1588         calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
1589
1590         // See if any dangling triangles
1591         boolList keepTriangles;
1592         label nDangling = markDanglingTriangles
1593         (
1594             faceEdges,
1595             edgeFace0,
1596             edgeFace1,
1597             edgeFacesRest,
1598             keepTriangles
1599         );
1600
1601         if (debug)
1602         {
1603             Pout<< "isoSurfaceCell : detected " << nDangling
1604                 << " dangling triangles." << endl;
1605         }
1606
1607         if (nDangling == 0)
1608         {
1609             break;
1610         }
1611
1612         // Create face map (new to old)
1613         labelList subsetTriMap(findIndices(keepTriangles, true));
1614
1615         labelList subsetPointMap;
1616         labelList reversePointMap;
1617         triSurface::operator=
1618         (
1619             subsetMesh
1620             (
1621                 *this,
1622                 subsetTriMap,
1623                 reversePointMap,
1624                 subsetPointMap
1625             )
1626         );
1627         meshCells_ = labelField(meshCells_, subsetTriMap);
1628         inplaceRenumber(reversePointMap, triPointMergeMap_);
1629     }
1630
1631     orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
1632 }
1633
1634     //combineCellTriangles();
1635 }
1636
1637
1638 ////XXXXXXX
1639 //// Experimental retriangulation of triangles per cell. Problem is that
1640 //// -it is very expensive   -only gets rid of internal points, not of boundary
1641 //// ones so limited benefit (e.g. 60 v.s. 88 triangles)
1642 //void Foam::isoSurfaceCell::combineCellTriangles()
1643 //{
1644 //    if (size())
1645 //    {
1646 //        DynamicList<labelledTri> newTris(size());
1647 //        DynamicList<label> newTriToCell(size());
1648 //
1649 //        label startTriI = 0;
1650 //
1651 //        DynamicList<labelledTri> tris;
1652 //
1653 //        for (label triI = 1; triI <= meshCells_.size(); triI++)
1654 //        {
1655 //            if
1656 //            (
1657 //                triI == meshCells_.size()
1658 //             || meshCells_[triI] != meshCells_[startTriI]
1659 //            )
1660 //            {
1661 //                label nTris = triI-startTriI;
1662 //
1663 //                if (nTris == 1)
1664 //                {
1665 //                    newTris.append(operator[](startTriI));
1666 //                    newTriToCell.append(meshCells_[startTriI]);
1667 //                }
1668 //                else
1669 //                {
1670 //                    // Collect from startTriI to triI in a triSurface
1671 //                    tris.clear();
1672 //                    for (label i = startTriI; i < triI; i++)
1673 //                    {
1674 //                        tris.append(operator[](i));
1675 //                    }
1676 //                    triSurface cellTris(tris, patches(), points());
1677 //                    tris.clear();
1678 //
1679 //                    // Get outside
1680 //                    const labelListList& loops = cellTris.edgeLoops();
1681 //
1682 //                    forAll(loops, i)
1683 //                    {
1684 //                        // Do proper triangulation of loop
1685 //                        face loop(renumber(cellTris.meshPoints(), loops[i]));
1686 //
1687 //                        faceTriangulation faceTris
1688 //                        (
1689 //                            points(),
1690 //                            loop,
1691 //                            true
1692 //                        );
1693 //
1694 //                        // Copy into newTris
1695 //                        forAll(faceTris, faceTriI)
1696 //                        {
1697 //                            const triFace& tri = faceTris[faceTriI];
1698 //
1699 //                            newTris.append
1700 //                            (
1701 //                                labelledTri
1702 //                                (
1703 //                                    tri[0],
1704 //                                    tri[1],
1705 //                                    tri[2],
1706 //                                    operator[](startTriI).region()
1707 //                                )
1708 //                            );
1709 //                            newTriToCell.append(meshCells_[startTriI]);
1710 //                        }
1711 //                    }
1712 //                }
1713 //
1714 //                startTriI = triI;
1715 //            }
1716 //        }
1717 //        newTris.shrink();
1718 //        newTriToCell.shrink();
1719 //
1720 //        // Compact
1721 //        pointField newPoints(points().size());
1722 //        label newPointI = 0;
1723 //        labelList oldToNewPoint(points().size(), -1);
1724 //
1725 //        forAll(newTris, i)
1726 //        {
1727 //            labelledTri& tri = newTris[i];
1728 //            forAll(tri, j)
1729 //            {
1730 //                label pointI = tri[j];
1731 //
1732 //                if (oldToNewPoint[pointI] == -1)
1733 //                {
1734 //                    oldToNewPoint[pointI] = newPointI;
1735 //                    newPoints[newPointI++] = points()[pointI];
1736 //                }
1737 //                tri[j] = oldToNewPoint[pointI];
1738 //            }
1739 //        }
1740 //        newPoints.setSize(newPointI);
1741 //
1742 //        triSurface::operator=
1743 //        (
1744 //            triSurface
1745 //            (
1746 //                newTris,
1747 //                patches(),
1748 //                newPoints,
1749 //                true
1750 //            )
1751 //        );
1752 //        meshCells_.transfer(newTriToCell);
1753 //    }
1754 //}
1755 ////XXXXXXX
1756
1757 // ************************************************************************* //