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