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