src/sampling/sampledSurface/isoSurface/isoSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "isoSurface.H"
  28 #include "fvMesh.H"
  29 #include "mergePoints.H"
  30 #include "syncTools.H"
  31 #include "addToRunTimeSelectionTable.H"
  32 #include "slicedVolFields.H"
  33 #include "volFields.H"
  34 #include "surfaceFields.H"
  35 #include "OFstream.H"
  36 #include "meshTools.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(isoSurface, 0);
  43 }
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 bool Foam::isoSurface::noTransform(const tensor& tt) const
  48 {
  49     return
  50         (mag(tt.xx()-1) < mergeDistance_)
  51      && (mag(tt.yy()-1) < mergeDistance_)
  52      && (mag(tt.zz()-1) < mergeDistance_)
  53      && (mag(tt.xy()) < mergeDistance_)
  54      && (mag(tt.xz()) < mergeDistance_)
  55      && (mag(tt.yx()) < mergeDistance_)
  56      && (mag(tt.yz()) < mergeDistance_)
  57      && (mag(tt.zx()) < mergeDistance_)
  58      && (mag(tt.zy()) < mergeDistance_);
  59 }
  60
  61
  62 bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
  63 {
  64     const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
  65
  66     return cpp.parallel() && !cpp.separated();
  67 }
  68
  69
  70 // Calculates per face whether couple is collocated.
  71 Foam::PackedBoolList Foam::isoSurface::collocatedFaces
  72 (
  73     const coupledPolyPatch& pp
  74 ) const
  75 {
  76     // Initialise to false
  77     PackedBoolList collocated(pp.size());
  78
  79     const vectorField& separation = pp.separation();
  80     const tensorField& forwardT = pp.forwardT();
  81
  82     if (forwardT.size() == 0)
  83     {
  84         // Parallel.
  85         if (separation.size() == 0)
  86         {
  87             collocated = 1u;
  88         }
  89         else if (separation.size() == 1)
  90         {
  91             // Fully separate. Do not synchronise.
  92         }
  93         else
  94         {
  95             // Per face separation.
  96             forAll(pp, faceI)
  97             {
  98                 if (mag(separation[faceI]) < mergeDistance_)
  99                 {
 100                     collocated[faceI] = 1u;
 101                 }
 102             }
 103         }
 104     }
 105     else if (forwardT.size() == 1)
 106     {
 107         // Fully transformed.
 108     }
 109     else
 110     {
 111         // Per face transformation.
 112         forAll(pp, faceI)
 113         {
 114             if (noTransform(forwardT[faceI]))
 115             {
 116                 collocated[faceI] = 1u;
 117             }
 118         }
 119     }
 120     return collocated;
 121 }
 122
 123
 124 // Insert the data for local point patchPointI into patch local values
 125 // and/or into the shared values field.
 126 void Foam::isoSurface::insertPointData
 127 (
 128     const processorPolyPatch& pp,
 129     const Map<label>& meshToShared,
 130     const pointField& pointValues,
 131     const label patchPointI,
 132     pointField& patchValues,
 133     pointField& sharedValues
 134 ) const
 135 {
 136     label meshPointI = pp.meshPoints()[patchPointI];
 137
 138     // Store in local field
 139     label nbrPointI = pp.neighbPoints()[patchPointI];
 140     if (nbrPointI >= 0 && nbrPointI < patchValues.size())
 141     {
 142         minEqOp<point>()(patchValues[nbrPointI], pointValues[meshPointI]);
 143     }
 144
 145     // Store in shared field
 146     Map<label>::const_iterator iter = meshToShared.find(meshPointI);
 147     if (iter != meshToShared.end())
 148     {
 149         minEqOp<point>()(sharedValues[iter()], pointValues[meshPointI]);
 150     }
 151 }
 152
 153
 154 void Foam::isoSurface::syncUnseparatedPoints
 155 (
 156     pointField& pointValues,
 157     const point& nullValue
 158 ) const
 159 {
 160     // Until syncPointList handles separated coupled patches with multiple
 161     // transforms do our own synchronisation of non-separated patches only
 162     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 163     const globalMeshData& pd = mesh_.globalData();
 164     const labelList& sharedPtAddr = pd.sharedPointAddr();
 165     const labelList& sharedPtLabels = pd.sharedPointLabels();
 166
 167     // Create map from meshPoint to globalShared index.
 168     Map<label> meshToShared(2*sharedPtLabels.size());
 169     forAll(sharedPtLabels, i)
 170     {
 171         meshToShared.insert(sharedPtLabels[i], sharedPtAddr[i]);
 172     }
 173
 174     // Values on shared points.
 175     pointField sharedInfo(pd.nGlobalPoints(), nullValue);
 176
 177
 178     if (Pstream::parRun())
 179     {
 180         // Send
 181         forAll(patches, patchI)
 182         {
 183             if
 184             (
 185                 isA<processorPolyPatch>(patches[patchI])
 186              && patches[patchI].nPoints() > 0
 187             )
 188             {
 189                 const processorPolyPatch& pp =
 190                     refCast<const processorPolyPatch>(patches[patchI]);
 191                 const labelList& meshPts = pp.meshPoints();
 192
 193                 pointField patchInfo(meshPts.size(), nullValue);
 194
 195                 PackedBoolList isCollocated(collocatedFaces(pp));
 196
 197                 forAll(isCollocated, faceI)
 198                 {
 199                     if (isCollocated[faceI])
 200                     {
 201                         const face& f = pp.localFaces()[faceI];
 202
 203                         forAll(f, fp)
 204                         {
 205                             label pointI = f[fp];
 206
 207                             insertPointData
 208                             (
 209                                 pp,
 210                                 meshToShared,
 211                                 pointValues,
 212                                 pointI,
 213                                 patchInfo,
 214                                 sharedInfo
 215                             );
 216                         }
 217                     }
 218                 }
 219
 220                 OPstream toNbr(Pstream::blocking, pp.neighbProcNo());
 221                 toNbr << patchInfo;
 222             }
 223         }
 224
 225         // Receive and combine.
 226
 227         forAll(patches, patchI)
 228         {
 229             if
 230             (
 231                 isA<processorPolyPatch>(patches[patchI])
 232              && patches[patchI].nPoints() > 0
 233             )
 234             {
 235                 const processorPolyPatch& pp =
 236                     refCast<const processorPolyPatch>(patches[patchI]);
 237
 238                 pointField nbrPatchInfo(pp.nPoints());
 239                 {
 240                     // We do not know the number of points on the other side
 241                     // so cannot use Pstream::read.
 242                     IPstream fromNbr(Pstream::blocking, pp.neighbProcNo());
 243                     fromNbr >> nbrPatchInfo;
 244                 }
 245
 246                 // Null any value which is not on neighbouring processor
 247                 nbrPatchInfo.setSize(pp.nPoints(), nullValue);
 248
 249                 const labelList& meshPts = pp.meshPoints();
 250
 251                 forAll(meshPts, pointI)
 252                 {
 253                     label meshPointI = meshPts[pointI];
 254                     minEqOp<point>()
 255                     (
 256                         pointValues[meshPointI],
 257                         nbrPatchInfo[pointI]
 258                     );
 259                 }
 260             }
 261         }
 262     }
 263
 264
 265     // Don't do cyclics for now. Are almost always separated anyway.
 266
 267
 268     // Shared points
 269
 270     // Combine on master.
 271     Pstream::listCombineGather(sharedInfo, minEqOp<point>());
 272     Pstream::listCombineScatter(sharedInfo);
 273
 274     // Now we will all have the same information. Merge it back with
 275     // my local information. (Note assignment and not combine operator)
 276     forAll(sharedPtLabels, i)
 277     {
 278         label meshPointI = sharedPtLabels[i];
 279         pointValues[meshPointI] = sharedInfo[sharedPtAddr[i]];
 280     }
 281 }
 282
 283
 284 Foam::scalar Foam::isoSurface::isoFraction
 285 (
 286     const scalar s0,
 287     const scalar s1
 288 ) const
 289 {
 290     scalar d = s1-s0;
 291
 292     if (mag(d) > VSMALL)
 293     {
 294         return (iso_-s0)/d;
 295     }
 296     else
 297     {
 298         return -1.0;
 299     }
 300 }
 301
 302
 303 bool Foam::isoSurface::isEdgeOfFaceCut
 304 (
 305     const scalarField& pVals,
 306     const face& f,
 307     const bool ownLower,
 308     const bool neiLower
 309 ) const
 310 {
 311     forAll(f, fp)
 312     {
 313         bool fpLower = (pVals[f[fp]] < iso_);
 314         if
 315         (
 316             (fpLower != ownLower)
 317          || (fpLower != neiLower)
 318          || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
 319         )
 320         {
 321             return true;
 322         }
 323     }
 324     return false;
 325 }
 326
 327
 328 // Get neighbour value and position.
 329 void Foam::isoSurface::getNeighbour
 330 (
 331     const labelList& boundaryRegion,
 332     const volVectorField& meshC,
 333     const volScalarField& cVals,
 334     const label cellI,
 335     const label faceI,
 336     scalar& nbrValue,
 337     point& nbrPoint
 338 ) const
 339 {
 340     const labelList& own = mesh_.faceOwner();
 341     const labelList& nei = mesh_.faceNeighbour();
 342
 343     if (mesh_.isInternalFace(faceI))
 344     {
 345         label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
 346         nbrValue = cVals[nbr];
 347         nbrPoint = meshC[nbr];
 348     }
 349     else
 350     {
 351         label bFaceI = faceI-mesh_.nInternalFaces();
 352         label patchI = boundaryRegion[bFaceI];
 353         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
 354         label patchFaceI = faceI-pp.start();
 355
 356         nbrValue = cVals.boundaryField()[patchI][patchFaceI];
 357         nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
 358     }
 359 }
 360
 361
 362 // Determine for every face/cell whether it (possibly) generates triangles.
 363 void Foam::isoSurface::calcCutTypes
 364 (
 365     const labelList& boundaryRegion,
 366     const volVectorField& meshC,
 367     const volScalarField& cVals,
 368     const scalarField& pVals
 369 )
 370 {
 371     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 372     const labelList& own = mesh_.faceOwner();
 373     const labelList& nei = mesh_.faceNeighbour();
 374
 375     faceCutType_.setSize(mesh_.nFaces());
 376     faceCutType_ = NOTCUT;
 377
 378     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 379     {
 380         // CC edge.
 381         bool ownLower = (cVals[own[faceI]] < iso_);
 382
 383         scalar nbrValue;
 384         point nbrPoint;
 385         getNeighbour
 386         (
 387             boundaryRegion,
 388             meshC,
 389             cVals,
 390             own[faceI],
 391             faceI,
 392             nbrValue,
 393             nbrPoint
 394         );
 395
 396         bool neiLower = (nbrValue < iso_);
 397
 398         if (ownLower != neiLower)
 399         {
 400             faceCutType_[faceI] = CUT;
 401         }
 402         else
 403         {
 404             // See if any mesh edge is cut by looping over all the edges of the
 405             // face.
 406             const face f = mesh_.faces()[faceI];
 407
 408             if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
 409             {
 410                 faceCutType_[faceI] = CUT;
 411             }
 412         }
 413     }
 414
 415     forAll(patches, patchI)
 416     {
 417         const polyPatch& pp = patches[patchI];
 418
 419         label faceI = pp.start();
 420
 421         forAll(pp, i)
 422         {
 423             bool ownLower = (cVals[own[faceI]] < iso_);
 424
 425             scalar nbrValue;
 426             point nbrPoint;
 427             getNeighbour
 428             (
 429                 boundaryRegion,
 430                 meshC,
 431                 cVals,
 432                 own[faceI],
 433                 faceI,
 434                 nbrValue,
 435                 nbrPoint
 436             );
 437
 438             bool neiLower = (nbrValue < iso_);
 439
 440             if (ownLower != neiLower)
 441             {
 442                 faceCutType_[faceI] = CUT;
 443             }
 444             else
 445             {
 446                 // Mesh edge.
 447                 const face f = mesh_.faces()[faceI];
 448
 449                 if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
 450                 {
 451                     faceCutType_[faceI] = CUT;
 452                 }
 453             }
 454
 455             faceI++;
 456         }
 457     }
 458
 459
 460
 461     nCutCells_ = 0;
 462     cellCutType_.setSize(mesh_.nCells());
 463     cellCutType_ = NOTCUT;
 464
 465     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 466     {
 467         if (faceCutType_[faceI] != NOTCUT)
 468         {
 469             if (cellCutType_[own[faceI]] == NOTCUT)
 470             {
 471                 cellCutType_[own[faceI]] = CUT;
 472                 nCutCells_++;
 473             }
 474             if (cellCutType_[nei[faceI]] == NOTCUT)
 475             {
 476                 cellCutType_[nei[faceI]] = CUT;
 477                 nCutCells_++;
 478             }
 479         }
 480     }
 481     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 482     {
 483         if (faceCutType_[faceI] != NOTCUT)
 484         {
 485             if (cellCutType_[own[faceI]] == NOTCUT)
 486             {
 487                 cellCutType_[own[faceI]] = CUT;
 488                 nCutCells_++;
 489             }
 490         }
 491     }
 492
 493     if (debug)
 494     {
 495         Pout<< "isoSurface : detected " << nCutCells_
 496             << " candidate cut cells (out of " << mesh_.nCells()
 497             << ")." << endl;
 498     }
 499 }
 500
 501
 502 // Return the two common points between two triangles
 503 Foam::labelPair Foam::isoSurface::findCommonPoints
 504 (
 505     const labelledTri& tri0,
 506     const labelledTri& tri1
 507 )
 508 {
 509     labelPair common(-1, -1);
 510
 511     label fp0 = 0;
 512     label fp1 = findIndex(tri1, tri0[fp0]);
 513
 514     if (fp1 == -1)
 515     {
 516         fp0 = 1;
 517         fp1 = findIndex(tri1, tri0[fp0]);
 518     }
 519
 520     if (fp1 != -1)
 521     {
 522         // So tri0[fp0] is tri1[fp1]
 523
 524         // Find next common point
 525         label fp0p1 = tri0.fcIndex(fp0);
 526         label fp1p1 = tri1.fcIndex(fp1);
 527         label fp1m1 = tri1.rcIndex(fp1);
 528
 529         if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
 530         {
 531             common[0] = tri0[fp0];
 532             common[1] = tri0[fp0p1];
 533         }
 534     }
 535     return common;
 536 }
 537
 538
 539 // Caculate centre of surface.
 540 Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 541 {
 542     vector sum = vector::zero;
 543
 544     forAll(s, i)
 545     {
 546         sum += s[i].centre(s.points());
 547     }
 548     return sum/s.size();
 549 }
 550
 551
 552 // Replace surface (localPoints, localTris) with single point. Returns
 553 // point. Destructs arguments.
 554 Foam::pointIndexHit Foam::isoSurface::collapseSurface
 555 (
 556     pointField& localPoints,
 557     DynamicList<labelledTri, 64>& localTris
 558 )
 559 {
 560     pointIndexHit info(false, vector::zero, localTris.size());
 561
 562     if (localTris.size() == 1)
 563     {
 564         const labelledTri& tri = localTris[0];
 565         info.setPoint(tri.centre(localPoints));
 566         info.setHit();
 567     }
 568     else if (localTris.size() == 2)
 569     {
 570         // Check if the two triangles share an edge.
 571         const labelledTri& tri0 = localTris[0];
 572         const labelledTri& tri1 = localTris[0];
 573
 574         labelPair shared = findCommonPoints(tri0, tri1);
 575
 576         if (shared[0] != -1)
 577         {
 578             info.setPoint
 579             (
 580                 0.5
 581               * (
 582                     tri0.centre(localPoints)
 583                   + tri1.centre(localPoints)
 584                 )
 585             );
 586             info.setHit();
 587         }
 588     }
 589     else if (localTris.size())
 590     {
 591         // Check if single region. Rare situation.
 592         triSurface surf
 593         (
 594             localTris,
 595             geometricSurfacePatchList(0),
 596             localPoints,
 597             true
 598         );
 599         localTris.clearStorage();
 600
 601         labelList faceZone;
 602         label nZones = surf.markZones
 603         (
 604             boolList(surf.nEdges(), false),
 605             faceZone
 606         );
 607
 608         if (nZones == 1)
 609         {
 610             info.setPoint(calcCentre(surf));
 611             info.setHit();
 612         }
 613     }
 614
 615     return info;
 616 }
 617
 618
 619 // Determine per cell centre whether all the intersections get collapsed
 620 // to a single point
 621 void Foam::isoSurface::calcSnappedCc
 622 (
 623     const labelList& boundaryRegion,
 624     const volVectorField& meshC,
 625     const volScalarField& cVals,
 626     const scalarField& pVals,
 627
 628     DynamicList<point>& snappedPoints,
 629     labelList& snappedCc
 630 ) const
 631 {
 632     const pointField& pts = mesh_.points();
 633     const pointField& cc = mesh_.cellCentres();
 634
 635     snappedCc.setSize(mesh_.nCells());
 636     snappedCc = -1;
 637
 638     // Work arrays
 639     DynamicList<point, 64> localTriPoints(64);
 640
 641     forAll(mesh_.cells(), cellI)
 642     {
 643         if (cellCutType_[cellI] == CUT)
 644         {
 645             scalar cVal = cVals[cellI];
 646
 647             const cell& cFaces = mesh_.cells()[cellI];
 648
 649             localTriPoints.clear();
 650             label nOther = 0;
 651             point otherPointSum = vector::zero;
 652
 653             // Create points for all intersections close to cell centre
 654             // (i.e. from pyramid edges)
 655
 656             forAll(cFaces, cFaceI)
 657             {
 658                 label faceI = cFaces[cFaceI];
 659
 660                 scalar nbrValue;
 661                 point nbrPoint;
 662                 getNeighbour
 663                 (
 664                     boundaryRegion,
 665                     meshC,
 666                     cVals,
 667                     cellI,
 668                     faceI,
 669                     nbrValue,
 670                     nbrPoint
 671                 );
 672
 673                 // Calculate intersection points of edges to cell centre
 674                 FixedList<scalar, 3> s;
 675                 FixedList<point, 3> pt;
 676
 677                 // From cc to neighbour cc.
 678                 s[2] = isoFraction(cVal, nbrValue);
 679                 pt[2] = (1.0-s[2])*cc[cellI] + s[2]*nbrPoint;
 680
 681                 const face& f = mesh_.faces()[cFaces[cFaceI]];
 682
 683                 forAll(f, fp)
 684                 {
 685                     // From cc to fp
 686                     label p0 = f[fp];
 687                     s[0] = isoFraction(cVal, pVals[p0]);
 688                     pt[0] = (1.0-s[0])*cc[cellI] + s[0]*pts[p0];
 689
 690                     // From cc to fp+1
 691                     label p1 = f[f.fcIndex(fp)];
 692                     s[1] = isoFraction(cVal, pVals[p1]);
 693                     pt[1] = (1.0-s[1])*cc[cellI] + s[1]*pts[p1];
 694
 695                     if
 696                     (
 697                         (s[0] >= 0.0 && s[0] <= 0.5)
 698                      && (s[1] >= 0.0 && s[1] <= 0.5)
 699                      && (s[2] >= 0.0 && s[2] <= 0.5)
 700                     )
 701                     {
 702                         localTriPoints.append(pt[0]);
 703                         localTriPoints.append(pt[1]);
 704                         localTriPoints.append(pt[2]);
 705                     }
 706                     else
 707                     {
 708                         // Get average of all other points
 709                         forAll(s, i)
 710                         {
 711                             if (s[i] >= 0.0 && s[i] <= 0.5)
 712                             {
 713                                 otherPointSum += pt[i];
 714                                 nOther++;
 715                             }
 716                         }
 717                     }
 718                 }
 719             }
 720
 721             if (localTriPoints.size() == 0)
 722             {
 723                 // No complete triangles. Get average of all intersection
 724                 // points.
 725                 if (nOther > 0)
 726                 {
 727                     snappedCc[cellI] = snappedPoints.size();
 728                     snappedPoints.append(otherPointSum/nOther);
 729
 730                     //Pout<< "    point:" << pointI
 731                     //    << " replacing coord:" << mesh_.points()[pointI]
 732                     //    << " by average:" << collapsedPoint[pointI] << endl;
 733                 }
 734             }
 735             else if (localTriPoints.size() == 3)
 736             {
 737                 // Single triangle. No need for any analysis. Average points.
 738                 pointField points;
 739                 points.transfer(localTriPoints);
 740                 snappedCc[cellI] = snappedPoints.size();
 741                 snappedPoints.append(sum(points)/points.size());
 742
 743                 //Pout<< "    point:" << pointI
 744                 //    << " replacing coord:" << mesh_.points()[pointI]
 745                 //    << " by average:" << collapsedPoint[pointI] << endl;
 746             }
 747             else
 748             {
 749                 // Convert points into triSurface.
 750
 751                 // Merge points and compact out non-valid triangles
 752                 labelList triMap;               // merged to unmerged triangle
 753                 labelList triPointReverseMap;   // unmerged to merged point
 754                 triSurface surf
 755                 (
 756                     stitchTriPoints
 757                     (
 758                         false,              // do not check for duplicate tris
 759                         localTriPoints,
 760                         triPointReverseMap,
 761                         triMap
 762                     )
 763                 );
 764
 765                 labelList faceZone;
 766                 label nZones = surf.markZones
 767                 (
 768                     boolList(surf.nEdges(), false),
 769                     faceZone
 770                 );
 771
 772                 if (nZones == 1)
 773                 {
 774                     snappedCc[cellI] = snappedPoints.size();
 775                     snappedPoints.append(calcCentre(surf));
 776                     //Pout<< "    point:" << pointI << " nZones:" << nZones
 777                     //    << " replacing coord:" << mesh_.points()[pointI]
 778                     //    << " by average:" << collapsedPoint[pointI] << endl;
 779                 }
 780             }
 781         }
 782     }
 783 }
 784
 785
 786 // Determine per meshpoint whether all the intersections get collapsed
 787 // to a single point
 788 void Foam::isoSurface::calcSnappedPoint
 789 (
 790     const PackedBoolList& isBoundaryPoint,
 791     const labelList& boundaryRegion,
 792     const volVectorField& meshC,
 793     const volScalarField& cVals,
 794     const scalarField& pVals,
 795
 796     DynamicList<point>& snappedPoints,
 797     labelList& snappedPoint
 798 ) const
 799 {
 800     const pointField& pts = mesh_.points();
 801     const pointField& cc = mesh_.cellCentres();
 802
 803
 804     const point greatPoint(VGREAT, VGREAT, VGREAT);
 805     pointField collapsedPoint(mesh_.nPoints(), greatPoint);
 806
 807
 808     // Work arrays
 809     DynamicList<point, 64> localTriPoints(100);
 810
 811     forAll(mesh_.pointFaces(), pointI)
 812     {
 813         if (isBoundaryPoint.get(pointI) == 1)
 814         {
 815             continue;
 816         }
 817
 818         const labelList& pFaces = mesh_.pointFaces()[pointI];
 819
 820         bool anyCut = false;
 821
 822         forAll(pFaces, i)
 823         {
 824             label faceI = pFaces[i];
 825
 826             if (faceCutType_[faceI] == CUT)
 827             {
 828                 anyCut = true;
 829                 break;
 830             }
 831         }
 832
 833         if (!anyCut)
 834         {
 835             continue;
 836         }
 837
 838
 839         localTriPoints.clear();
 840         label nOther = 0;
 841         point otherPointSum = vector::zero;
 842
 843         forAll(pFaces, pFaceI)
 844         {
 845             // Create points for all intersections close to point
 846             // (i.e. from pyramid edges)
 847
 848             label faceI = pFaces[pFaceI];
 849             const face& f = mesh_.faces()[faceI];
 850             label own = mesh_.faceOwner()[faceI];
 851
 852             // Get neighbour value
 853             scalar nbrValue;
 854             point nbrPoint;
 855             getNeighbour
 856             (
 857                 boundaryRegion,
 858                 meshC,
 859                 cVals,
 860                 own,
 861                 faceI,
 862                 nbrValue,
 863                 nbrPoint
 864             );
 865
 866             // Calculate intersection points of edges emanating from point
 867             FixedList<scalar, 4> s;
 868             FixedList<point, 4> pt;
 869
 870             label fp = findIndex(f, pointI);
 871             s[0] = isoFraction(pVals[pointI], cVals[own]);
 872             pt[0] = (1.0-s[0])*pts[pointI] + s[0]*cc[own];
 873
 874             s[1] = isoFraction(pVals[pointI], nbrValue);
 875             pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint;
 876
 877             label nextPointI = f[f.fcIndex(fp)];
 878             s[2] = isoFraction(pVals[pointI], pVals[nextPointI]);
 879             pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI];
 880
 881             label prevPointI = f[f.rcIndex(fp)];
 882             s[3] = isoFraction(pVals[pointI], pVals[prevPointI]);
 883             pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI];
 884
 885             if
 886             (
 887                 (s[0] >= 0.0 && s[0] <= 0.5)
 888              && (s[1] >= 0.0 && s[1] <= 0.5)
 889              && (s[2] >= 0.0 && s[2] <= 0.5)
 890             )
 891             {
 892                 localTriPoints.append(pt[0]);
 893                 localTriPoints.append(pt[1]);
 894                 localTriPoints.append(pt[2]);
 895             }
 896             if
 897             (
 898                 (s[0] >= 0.0 && s[0] <= 0.5)
 899              && (s[1] >= 0.0 && s[1] <= 0.5)
 900              && (s[3] >= 0.0 && s[3] <= 0.5)
 901             )
 902             {
 903                 localTriPoints.append(pt[3]);
 904                 localTriPoints.append(pt[0]);
 905                 localTriPoints.append(pt[1]);
 906             }
 907
 908             // Get average of points as fallback
 909             forAll(s, i)
 910             {
 911                 if (s[i] >= 0.0 && s[i] <= 0.5)
 912                 {
 913                     otherPointSum += pt[i];
 914                     nOther++;
 915                 }
 916             }
 917         }
 918
 919         if (localTriPoints.size() == 0)
 920         {
 921             // No complete triangles. Get average of all intersection
 922             // points.
 923             if (nOther > 0)
 924             {
 925                 collapsedPoint[pointI] = otherPointSum/nOther;
 926             }
 927         }
 928         else if (localTriPoints.size() == 3)
 929         {
 930             // Single triangle. No need for any analysis. Average points.
 931             pointField points;
 932             points.transfer(localTriPoints);
 933             collapsedPoint[pointI] = sum(points)/points.size();
 934         }
 935         else
 936         {
 937             // Convert points into triSurface.
 938
 939             // Merge points and compact out non-valid triangles
 940             labelList triMap;               // merged to unmerged triangle
 941             labelList triPointReverseMap;   // unmerged to merged point
 942             triSurface surf
 943             (
 944                 stitchTriPoints
 945                 (
 946                     false,                  // do not check for duplicate tris
 947                     localTriPoints,
 948                     triPointReverseMap,
 949                     triMap
 950                 )
 951             );
 952
 953             labelList faceZone;
 954             label nZones = surf.markZones
 955             (
 956                 boolList(surf.nEdges(), false),
 957                 faceZone
 958             );
 959
 960             if (nZones == 1)
 961             {
 962                 collapsedPoint[pointI] = calcCentre(surf);
 963             }
 964         }
 965     }
 966
 967
 968     // Synchronise snap location
 969     syncUnseparatedPoints(collapsedPoint, greatPoint);
 970
 971
 972     snappedPoint.setSize(mesh_.nPoints());
 973     snappedPoint = -1;
 974
 975     forAll(collapsedPoint, pointI)
 976     {
 977         if (collapsedPoint[pointI] != greatPoint)
 978         {
 979             snappedPoint[pointI] = snappedPoints.size();
 980             snappedPoints.append(collapsedPoint[pointI]);
 981         }
 982     }
 983 }
 984
 985
 986 Foam::triSurface Foam::isoSurface::stitchTriPoints
 987 (
 988     const bool checkDuplicates,
 989     const List<point>& triPoints,
 990     labelList& triPointReverseMap,  // unmerged to merged point
 991     labelList& triMap               // merged to unmerged triangle
 992 ) const
 993 {
 994     label nTris = triPoints.size()/3;
 995
 996     if ((triPoints.size() % 3) != 0)
 997     {
 998         FatalErrorIn("isoSurface::stitchTriPoints(..)")
 999             << "Problem: number of points " << triPoints.size()
1000             << " not a multiple of 3." << abort(FatalError);
1001     }
1002
1003     pointField newPoints;
1004     mergePoints
1005     (
1006         triPoints,
1007         mergeDistance_,
1008         false,
1009         triPointReverseMap,
1010         newPoints
1011     );
1012
1013     // Check that enough merged.
1014     if (debug)
1015     {
1016         pointField newNewPoints;
1017         labelList oldToNew;
1018         bool hasMerged = mergePoints
1019         (
1020             newPoints,
1021             mergeDistance_,
1022             true,
1023             oldToNew,
1024             newNewPoints
1025         );
1026
1027         if (hasMerged)
1028         {
1029             FatalErrorIn("isoSurface::stitchTriPoints(..)")
1030                 << "Merged points contain duplicates"
1031                 << " when merging with distance " << mergeDistance_ << endl
1032                 << "merged:" << newPoints.size() << " re-merged:"
1033                 << newNewPoints.size()
1034                 << abort(FatalError);
1035         }
1036     }
1037
1038
1039     List<labelledTri> tris;
1040     {
1041         DynamicList<labelledTri> dynTris(nTris);
1042         label rawPointI = 0;
1043         DynamicList<label> newToOldTri(nTris);
1044
1045         for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
1046         {
1047             labelledTri tri
1048             (
1049                 triPointReverseMap[rawPointI],
1050                 triPointReverseMap[rawPointI+1],
1051                 triPointReverseMap[rawPointI+2],
1052                 0
1053             );
1054             rawPointI += 3;
1055
1056             if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
1057             {
1058                 newToOldTri.append(oldTriI);
1059                 dynTris.append(tri);
1060             }
1061         }
1062
1063         triMap.transfer(newToOldTri);
1064         tris.transfer(dynTris);
1065     }
1066
1067
1068
1069     // Determine 'flat hole' situation (see RMT paper).
1070     // Two unconnected triangles get connected because (some of) the edges
1071     // separating them get collapsed. Below only checks for duplicate triangles,
1072     // not non-manifold edge connectivity.
1073     if (checkDuplicates)
1074     {
1075         labelListList pointFaces;
1076         invertManyToMany(newPoints.size(), tris, pointFaces);
1077
1078         // Filter out duplicates.
1079         DynamicList<label> newToOldTri(tris.size());
1080
1081         forAll(tris, triI)
1082         {
1083             const labelledTri& tri = tris[triI];
1084             const labelList& pFaces = pointFaces[tri[0]];
1085
1086             // Find the maximum of any duplicates. Maximum since the tris
1087             // below triI
1088             // get overwritten so we cannot use them in a comparison.
1089             label dupTriI = -1;
1090             forAll(pFaces, i)
1091             {
1092                 label nbrTriI = pFaces[i];
1093
1094                 if (nbrTriI > triI && (tris[nbrTriI] == tri))
1095                 {
1096                     //Pout<< "Duplicate : " << triI << " verts:" << tri
1097                     //    << " to " << nbrTriI << " verts:" << tris[nbrTriI]
1098                     //    << endl;
1099                     dupTriI = nbrTriI;
1100                     break;
1101                 }
1102             }
1103
1104             if (dupTriI == -1)
1105             {
1106                 // There is no (higher numbered) duplicate triangle
1107                 label newTriI = newToOldTri.size();
1108                 newToOldTri.append(triI);
1109                 tris[newTriI] = tris[triI];
1110             }
1111         }
1112
1113         triMap.transfer(newToOldTri);
1114         tris.setSize(triMap.size());
1115
1116         if (debug)
1117         {
1118             Pout<< "isoSurface : merged from " << nTris
1119                 << " down to " << tris.size() << " unique triangles." << endl;
1120         }
1121
1122
1123         if (debug)
1124         {
1125             triSurface surf(tris, geometricSurfacePatchList(0), newPoints);
1126
1127             forAll(surf, faceI)
1128             {
1129                 const labelledTri& f = surf[faceI];
1130                 const labelList& fFaces = surf.faceFaces()[faceI];
1131
1132                 forAll(fFaces, i)
1133                 {
1134                     label nbrFaceI = fFaces[i];
1135
1136                     if (nbrFaceI <= faceI)
1137                     {
1138                         // lower numbered faces already checked
1139                         continue;
1140                     }
1141
1142                     const labelledTri& nbrF = surf[nbrFaceI];
1143
1144                     if (f == nbrF)
1145                     {
1146                         FatalErrorIn("validTri(const triSurface&, const label)")
1147                             << "Check : "
1148                             << " triangle " << faceI << " vertices " << f
1149                             << " fc:" << f.centre(surf.points())
1150                             << " has the same vertices as triangle " << nbrFaceI
1151                             << " vertices " << nbrF
1152                             << " fc:" << nbrF.centre(surf.points())
1153                             << abort(FatalError);
1154                     }
1155                 }
1156             }
1157         }
1158     }
1159
1160     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
1161 }
1162
1163
1164 // Does face use valid vertices?
1165 bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
1166 {
1167     // Simple check on indices ok.
1168
1169     const labelledTri& f = surf[faceI];
1170
1171     if
1172     (
1173         (f[0] < 0) || (f[0] >= surf.points().size())
1174      || (f[1] < 0) || (f[1] >= surf.points().size())
1175      || (f[2] < 0) || (f[2] >= surf.points().size())
1176     )
1177     {
1178         WarningIn("validTri(const triSurface&, const label)")
1179             << "triangle " << faceI << " vertices " << f
1180             << " uses point indices outside point range 0.."
1181             << surf.points().size()-1 << endl;
1182
1183         return false;
1184     }
1185
1186     if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
1187     {
1188         WarningIn("validTri(const triSurface&, const label)")
1189             << "triangle " << faceI
1190             << " uses non-unique vertices " << f
1191             << endl;
1192         return false;
1193     }
1194
1195     // duplicate triangle check
1196
1197     const labelList& fFaces = surf.faceFaces()[faceI];
1198
1199     // Check if faceNeighbours use same points as this face.
1200     // Note: discards normal information - sides of baffle are merged.
1201     forAll(fFaces, i)
1202     {
1203         label nbrFaceI = fFaces[i];
1204
1205         if (nbrFaceI <= faceI)
1206         {
1207             // lower numbered faces already checked
1208             continue;
1209         }
1210
1211         const labelledTri& nbrF = surf[nbrFaceI];
1212
1213         if
1214         (
1215             ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
1216          && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
1217          && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
1218         )
1219         {
1220             WarningIn("validTri(const triSurface&, const label)")
1221                 << "triangle " << faceI << " vertices " << f
1222                 << " fc:" << f.centre(surf.points())
1223                 << " has the same vertices as triangle " << nbrFaceI
1224                 << " vertices " << nbrF
1225                 << " fc:" << nbrF.centre(surf.points())
1226                 << endl;
1227
1228             return false;
1229         }
1230     }
1231     return true;
1232 }
1233
1234
1235 void Foam::isoSurface::calcAddressing
1236 (
1237     const triSurface& surf,
1238     List<FixedList<label, 3> >& faceEdges,
1239     labelList& edgeFace0,
1240     labelList& edgeFace1,
1241     Map<labelList>& edgeFacesRest
1242 ) const
1243 {
1244     const pointField& points = surf.points();
1245
1246     pointField edgeCentres(3*surf.size());
1247     label edgeI = 0;
1248     forAll(surf, triI)
1249     {
1250         const labelledTri& tri = surf[triI];
1251         edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1252         edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1253         edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1254     }
1255
1256     pointField mergedCentres;
1257     labelList oldToMerged;
1258     bool hasMerged = mergePoints
1259     (
1260         edgeCentres,
1261         mergeDistance_,
1262         false,
1263         oldToMerged,
1264         mergedCentres
1265     );
1266
1267     if (debug)
1268     {
1269         Pout<< "isoSurface : detected "
1270             << mergedCentres.size()
1271             << " geometric edges on " << surf.size() << " triangles." << endl;
1272     }
1273
1274     if (!hasMerged)
1275     {
1276         if (surf.size() == 1)
1277         {
1278             faceEdges.setSize(1);
1279             faceEdges[0][0] = 0;
1280             faceEdges[0][1] = 1;
1281             faceEdges[0][2] = 2;
1282             edgeFace0.setSize(1);
1283             edgeFace0[0] = 0;
1284             edgeFace1.setSize(1);
1285             edgeFace1[0] = -1;
1286             edgeFacesRest.clear();
1287         }
1288         return;
1289     }
1290
1291
1292     // Determine faceEdges
1293     faceEdges.setSize(surf.size());
1294     edgeI = 0;
1295     forAll(surf, triI)
1296     {
1297         faceEdges[triI][0] = oldToMerged[edgeI++];
1298         faceEdges[triI][1] = oldToMerged[edgeI++];
1299         faceEdges[triI][2] = oldToMerged[edgeI++];
1300     }
1301
1302
1303     // Determine edgeFaces
1304     edgeFace0.setSize(mergedCentres.size());
1305     edgeFace0 = -1;
1306     edgeFace1.setSize(mergedCentres.size());
1307     edgeFace1 = -1;
1308     edgeFacesRest.clear();
1309
1310     // Overflow edge faces for geometric shared edges that turned
1311     // out to be different anyway.
1312     EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
1313
1314     forAll(oldToMerged, oldEdgeI)
1315     {
1316         label triI = oldEdgeI / 3;
1317         label edgeI = oldToMerged[oldEdgeI];
1318
1319         if (edgeFace0[edgeI] == -1)
1320         {
1321             // First triangle for edge
1322             edgeFace0[edgeI] = triI;
1323         }
1324         else
1325         {
1326             //- Check that the two triangles actually topologically
1327             //  share an edge
1328             const labelledTri& prevTri = surf[edgeFace0[edgeI]];
1329             const labelledTri& tri = surf[triI];
1330
1331             label fp = oldEdgeI % 3;
1332
1333             edge e(tri[fp], tri[tri.fcIndex(fp)]);
1334
1335             label prevTriIndex = -1;
1336
1337             forAll(prevTri, i)
1338             {
1339                 if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
1340                 {
1341                     prevTriIndex = i;
1342                     break;
1343                 }
1344             }
1345
1346             if (prevTriIndex == -1)
1347             {
1348                 // Different edge. Store for later.
1349                 EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
1350
1351                 if (iter != extraEdgeFaces.end())
1352                 {
1353                     labelList& eFaces = iter();
1354                     label sz = eFaces.size();
1355                     eFaces.setSize(sz+1);
1356                     eFaces[sz] = triI;
1357                 }
1358                 else
1359                 {
1360                     extraEdgeFaces.insert(e, labelList(1, triI));
1361                 }
1362             }
1363             else if (edgeFace1[edgeI] == -1)
1364             {
1365                 edgeFace1[edgeI] = triI;
1366             }
1367             else
1368             {
1369                 //WarningIn("orientSurface(triSurface&)")
1370                 //    << "Edge " << edgeI << " with centre "
1371                 //    << mergedCentres[edgeI]
1372                 //    << " used by more than two triangles: "
1373                 //    << edgeFace0[edgeI] << ", "
1374                 //    << edgeFace1[edgeI] << " and " << triI << endl;
1375                 Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
1376
1377                 if (iter != edgeFacesRest.end())
1378                 {
1379                     labelList& eFaces = iter();
1380                     label sz = eFaces.size();
1381                     eFaces.setSize(sz+1);
1382                     eFaces[sz] = triI;
1383                 }
1384                 else
1385                 {
1386                     edgeFacesRest.insert(edgeI, labelList(1, triI));
1387                 }
1388             }
1389         }
1390     }
1391
1392     // Add extraEdgeFaces
1393     edgeI = edgeFace0.size();
1394
1395     edgeFace0.setSize(edgeI + extraEdgeFaces.size());
1396     edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
1397
1398     forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
1399     {
1400         const labelList& eFaces = iter();
1401
1402         // The current edge will become edgeI. Replace all occurrences in
1403         // faceEdges
1404         forAll(eFaces, i)
1405         {
1406             label triI = eFaces[i];
1407             const labelledTri& tri = surf[triI];
1408
1409             FixedList<label, 3>& fEdges = faceEdges[triI];
1410             forAll(tri, fp)
1411             {
1412                 edge e(tri[fp], tri[tri.fcIndex(fp)]);
1413
1414                 if (e == iter.key())
1415                 {
1416                     fEdges[fp] = edgeI;
1417                     break;
1418                 }
1419             }
1420         }
1421
1422
1423         // Add face to edgeFaces
1424
1425         edgeFace0[edgeI] = eFaces[0];
1426
1427         if (eFaces.size() >= 2)
1428         {
1429             edgeFace1[edgeI] = eFaces[1];
1430
1431             if (eFaces.size() > 2)
1432             {
1433                 edgeFacesRest.insert
1434                 (
1435                     edgeI,
1436                     SubList<label>(eFaces, eFaces.size()-2, 2)
1437                 );
1438             }
1439         }
1440
1441         edgeI++;
1442     }
1443 }
1444
1445
1446 void Foam::isoSurface::walkOrientation
1447 (
1448     const triSurface& surf,
1449     const List<FixedList<label, 3> >& faceEdges,
1450     const labelList& edgeFace0,
1451     const labelList& edgeFace1,
1452     const label seedTriI,
1453     labelList& flipState
1454 )
1455 {
1456     // Do walk for consistent orientation.
1457     DynamicList<label> changedFaces(surf.size());
1458
1459     changedFaces.append(seedTriI);
1460
1461     while (changedFaces.size())
1462     {
1463         DynamicList<label> newChangedFaces(changedFaces.size());
1464
1465         forAll(changedFaces, i)
1466         {
1467             label triI = changedFaces[i];
1468             const labelledTri& tri = surf[triI];
1469             const FixedList<label, 3>& fEdges = faceEdges[triI];
1470
1471             forAll(fEdges, fp)
1472             {
1473                 label edgeI = fEdges[fp];
1474
1475                 // my points:
1476                 label p0 = tri[fp];
1477                 label p1 = tri[tri.fcIndex(fp)];
1478
1479                 label nbrI =
1480                 (
1481                     edgeFace0[edgeI] != triI
1482                   ? edgeFace0[edgeI]
1483                   : edgeFace1[edgeI]
1484                 );
1485
1486                 if (nbrI != -1 && flipState[nbrI] == -1)
1487                 {
1488                     const labelledTri& nbrTri = surf[nbrI];
1489
1490                     // nbr points
1491                     label nbrFp = findIndex(nbrTri, p0);
1492
1493                     if (nbrFp == -1)
1494                     {
1495                         FatalErrorIn("isoSurface::walkOrientation(..)")
1496                             << "triI:" << triI
1497                             << " tri:" << tri
1498                             << " p0:" << p0
1499                             << " p1:" << p1
1500                             << " fEdges:" << fEdges
1501                             << " edgeI:" << edgeI
1502                             << " edgeFace0:" << edgeFace0[edgeI]
1503                             << " edgeFace1:" << edgeFace1[edgeI]
1504                             << " nbrI:" << nbrI
1505                             << " nbrTri:" << nbrTri
1506                             << abort(FatalError);
1507                     }
1508
1509
1510                     label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1511
1512                     bool sameOrientation = (p1 == nbrP1);
1513
1514                     if (flipState[triI] == 0)
1515                     {
1516                         flipState[nbrI] = (sameOrientation ? 0 : 1);
1517                     }
1518                     else
1519                     {
1520                         flipState[nbrI] = (sameOrientation ? 1 : 0);
1521                     }
1522                     newChangedFaces.append(nbrI);
1523                 }
1524             }
1525         }
1526
1527         changedFaces.transfer(newChangedFaces);
1528     }
1529 }
1530
1531
1532 void Foam::isoSurface::orientSurface
1533 (
1534     triSurface& surf,
1535     const List<FixedList<label, 3> >& faceEdges,
1536     const labelList& edgeFace0,
1537     const labelList& edgeFace1,
1538     const Map<labelList>& edgeFacesRest
1539 )
1540 {
1541     // -1 : unvisited
1542     //  0 : leave as is
1543     //  1 : flip
1544     labelList flipState(surf.size(), -1);
1545
1546     label seedTriI = 0;
1547
1548     while (true)
1549     {
1550         // Find first unvisited triangle
1551         for
1552         (
1553             ;
1554             seedTriI < surf.size() && flipState[seedTriI] != -1;
1555             seedTriI++
1556         )
1557         {}
1558
1559         if (seedTriI == surf.size())
1560         {
1561             break;
1562         }
1563
1564         // Note: Determine orientation of seedTriI?
1565         // for now assume it is ok
1566         flipState[seedTriI] = 0;
1567
1568         walkOrientation
1569         (
1570             surf,
1571             faceEdges,
1572             edgeFace0,
1573             edgeFace1,
1574             seedTriI,
1575             flipState
1576         );
1577     }
1578
1579     // Do actual flipping
1580     surf.clearOut();
1581     forAll(surf, triI)
1582     {
1583         if (flipState[triI] == 1)
1584         {
1585             labelledTri tri(surf[triI]);
1586
1587             surf[triI][0] = tri[0];
1588             surf[triI][1] = tri[2];
1589             surf[triI][2] = tri[1];
1590         }
1591         else if (flipState[triI] == -1)
1592         {
1593             FatalErrorIn
1594             (
1595                 "isoSurface::orientSurface(triSurface&, const label)"
1596             )   << "problem" << abort(FatalError);
1597         }
1598     }
1599 }
1600
1601
1602 // Checks if triangle is connected through edgeI only.
1603 bool Foam::isoSurface::danglingTriangle
1604 (
1605     const FixedList<label, 3>& fEdges,
1606     const labelList& edgeFace1
1607 )
1608 {
1609     label nOpen = 0;
1610     forAll(fEdges, i)
1611     {
1612         if (edgeFace1[fEdges[i]] == -1)
1613         {
1614             nOpen++;
1615         }
1616     }
1617
1618     if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1619     {
1620         return true;
1621     }
1622     else
1623     {
1624         return false;
1625     }
1626 }
1627
1628
1629 // Mark triangles to keep. Returns number of dangling triangles.
1630 Foam::label Foam::isoSurface::markDanglingTriangles
1631 (
1632     const List<FixedList<label, 3> >& faceEdges,
1633     const labelList& edgeFace0,
1634     const labelList& edgeFace1,
1635     const Map<labelList>& edgeFacesRest,
1636     boolList& keepTriangles
1637 )
1638 {
1639     keepTriangles.setSize(faceEdges.size());
1640     keepTriangles = true;
1641
1642     label nDangling = 0;
1643
1644     // Remove any dangling triangles
1645     forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
1646     {
1647         // These are all the non-manifold edges. Filter out all triangles
1648         // with only one connected edge (= this edge)
1649
1650         label edgeI = iter.key();
1651         const labelList& otherEdgeFaces = iter();
1652
1653         // Remove all dangling triangles
1654         if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1655         {
1656             keepTriangles[edgeFace0[edgeI]] = false;
1657             nDangling++;
1658         }
1659         if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1660         {
1661             keepTriangles[edgeFace1[edgeI]] = false;
1662             nDangling++;
1663         }
1664         forAll(otherEdgeFaces, i)
1665         {
1666             label triI = otherEdgeFaces[i];
1667             if (danglingTriangle(faceEdges[triI], edgeFace1))
1668             {
1669                 keepTriangles[triI] = false;
1670                 nDangling++;
1671             }
1672         }
1673     }
1674     return nDangling;
1675 }
1676
1677
1678 Foam::triSurface Foam::isoSurface::subsetMesh
1679 (
1680     const triSurface& s,
1681     const labelList& newToOldFaces,
1682     labelList& oldToNewPoints,
1683     labelList& newToOldPoints
1684 )
1685 {
1686     const boolList include
1687     (
1688         createWithValues<boolList>
1689         (
1690             s.size(),
1691             false,
1692             newToOldFaces,
1693             true
1694         )
1695     );
1696
1697     newToOldPoints.setSize(s.points().size());
1698     oldToNewPoints.setSize(s.points().size());
1699     oldToNewPoints = -1;
1700     {
1701         label pointI = 0;
1702
1703         forAll(include, oldFacei)
1704         {
1705             if (include[oldFacei])
1706             {
1707                 // Renumber labels for triangle
1708                 const labelledTri& tri = s[oldFacei];
1709
1710                 forAll(tri, fp)
1711                 {
1712                     label oldPointI = tri[fp];
1713
1714                     if (oldToNewPoints[oldPointI] == -1)
1715                     {
1716                         oldToNewPoints[oldPointI] = pointI;
1717                         newToOldPoints[pointI++] = oldPointI;
1718                     }
1719                 }
1720             }
1721         }
1722         newToOldPoints.setSize(pointI);
1723     }
1724
1725     // Extract points
1726     pointField newPoints(newToOldPoints.size());
1727     forAll(newToOldPoints, i)
1728     {
1729         newPoints[i] = s.points()[newToOldPoints[i]];
1730     }
1731     // Extract faces
1732     List<labelledTri> newTriangles(newToOldFaces.size());
1733
1734     forAll(newToOldFaces, i)
1735     {
1736         // Get old vertex labels
1737         const labelledTri& tri = s[newToOldFaces[i]];
1738
1739         newTriangles[i][0] = oldToNewPoints[tri[0]];
1740         newTriangles[i][1] = oldToNewPoints[tri[1]];
1741         newTriangles[i][2] = oldToNewPoints[tri[2]];
1742         newTriangles[i].region() = tri.region();
1743     }
1744
1745     // Reuse storage.
1746     return triSurface(newTriangles, s.patches(), newPoints, true);
1747 }
1748
1749
1750 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1751
1752 Foam::isoSurface::isoSurface
1753 (
1754     const volScalarField& cVals,
1755     const scalarField& pVals,
1756     const scalar iso,
1757     const bool regularise,
1758     const scalar mergeTol
1759 )
1760 :
1761     mesh_(cVals.mesh()),
1762     pVals_(pVals),
1763     iso_(iso),
1764     regularise_(regularise),
1765     mergeDistance_(mergeTol*mesh_.bounds().mag())
1766 {
1767     if (debug)
1768     {
1769         Pout<< "isoSurface:" << nl
1770             << "    isoField      : " << cVals.name() << nl
1771             << "    cell min/max  : "
1772             << min(cVals.internalField()) << " / "
1773             << max(cVals.internalField()) << nl
1774             << "    point min/max : "
1775             << min(pVals_) << " / "
1776             << max(pVals_) << nl
1777             << "    isoValue      : " << iso << nl
1778             << "    regularise    : " << regularise_ << nl
1779             << "    mergeTol      : " << mergeTol << nl
1780             << endl;
1781     }
1782
1783     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1784
1785
1786     // Rewrite input field
1787     // ~~~~~~~~~~~~~~~~~~~
1788     // Rewrite input volScalarField to have interpolated values
1789     // on separated patches.
1790
1791     cValsPtr_.reset(adaptPatchFields(cVals).ptr());
1792
1793
1794     // Construct cell centres field consistent with cVals
1795     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1796     // Generate field to interpolate. This is identical to the mesh.C()
1797     // except on separated coupled patches and on empty patches.
1798
1799     slicedVolVectorField meshC
1800     (
1801         IOobject
1802         (
1803             "C",
1804             mesh_.pointsInstance(),
1805             mesh_.meshSubDir,
1806             mesh_,
1807             IOobject::NO_READ,
1808             IOobject::NO_WRITE,
1809             false
1810         ),
1811         mesh_,
1812         dimLength,
1813         mesh_.cellCentres(),
1814         mesh_.faceCentres()
1815     );
1816     forAll(patches, patchI)
1817     {
1818         const polyPatch& pp = patches[patchI];
1819
1820         // Adapt separated coupled (proc and cyclic) patches
1821         if (isA<coupledPolyPatch>(pp) && !collocatedPatch(pp))
1822         {
1823             fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
1824             (
1825                 meshC.boundaryField()[patchI]
1826             );
1827
1828             PackedBoolList isCollocated
1829             (
1830                 collocatedFaces(refCast<const coupledPolyPatch>(pp))
1831             );
1832
1833             forAll(isCollocated, i)
1834             {
1835                 if (!isCollocated[i])
1836                 {
1837                     pfld[i] = mesh_.faceCentres()[pp.start()+i];
1838                 }
1839             }
1840         }
1841         else if (isA<emptyPolyPatch>(pp))
1842         {
1843             typedef slicedVolVectorField::GeometricBoundaryField bType;
1844
1845             bType& bfld = const_cast<bType&>(meshC.boundaryField());
1846
1847             // Clear old value. Cannot resize it since is a slice.
1848             bfld.set(patchI, NULL);
1849
1850             // Set new value we can change
1851             bfld.set
1852             (
1853                 patchI,
1854                 new calculatedFvPatchField<vector>
1855                 (
1856                     mesh_.boundary()[patchI],
1857                     meshC
1858                 )
1859             );
1860
1861             // Change to face centres
1862             bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
1863         }
1864     }
1865
1866
1867
1868     // Pre-calculate patch-per-face to avoid whichPatch call.
1869     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1870
1871     forAll(patches, patchI)
1872     {
1873         const polyPatch& pp = patches[patchI];
1874
1875         label faceI = pp.start();
1876
1877         forAll(pp, i)
1878         {
1879             boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
1880             faceI++;
1881         }
1882     }
1883
1884
1885
1886     // Determine if any cut through face/cell
1887     calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
1888
1889
1890     DynamicList<point> snappedPoints(nCutCells_);
1891
1892     // Per cc -1 or a point inside snappedPoints.
1893     labelList snappedCc;
1894     if (regularise_)
1895     {
1896         calcSnappedCc
1897         (
1898             boundaryRegion,
1899             meshC,
1900             cValsPtr_(),
1901             pVals_,
1902
1903             snappedPoints,
1904             snappedCc
1905         );
1906     }
1907     else
1908     {
1909         snappedCc.setSize(mesh_.nCells());
1910         snappedCc = -1;
1911     }
1912
1913
1914
1915     if (debug)
1916     {
1917         Pout<< "isoSurface : shifted " << snappedPoints.size()
1918             << " cell centres to intersection." << endl;
1919     }
1920
1921     label nCellSnaps = snappedPoints.size();
1922
1923
1924     // Per point -1 or a point inside snappedPoints.
1925     labelList snappedPoint;
1926     if (regularise_)
1927     {
1928         // Determine if point is on boundary.
1929         PackedBoolList isBoundaryPoint(mesh_.nPoints());
1930
1931         forAll(patches, patchI)
1932         {
1933             // Mark all boundary points that are not physically coupled
1934             // (so anything but collocated coupled patches)
1935
1936             if (patches[patchI].coupled())
1937             {
1938                 if (!collocatedPatch(patches[patchI]))
1939                 {
1940                     const coupledPolyPatch& cpp =
1941                         refCast<const coupledPolyPatch>
1942                         (
1943                             patches[patchI]
1944                         );
1945
1946                     PackedBoolList isCollocated(collocatedFaces(cpp));
1947
1948                     forAll(isCollocated, i)
1949                     {
1950                         if (!isCollocated[i])
1951                         {
1952                             const face& f = mesh_.faces()[cpp.start()+i];
1953
1954                             forAll(f, fp)
1955                             {
1956                                 isBoundaryPoint.set(f[fp], 1);
1957                             }
1958                         }
1959                     }
1960                 }
1961             }
1962             else
1963             {
1964                 const polyPatch& pp = patches[patchI];
1965
1966                 forAll(pp, i)
1967                 {
1968                     const face& f = mesh_.faces()[pp.start()+i];
1969
1970                     forAll(f, fp)
1971                     {
1972                         isBoundaryPoint.set(f[fp], 1);
1973                     }
1974                 }
1975             }
1976         }
1977
1978         calcSnappedPoint
1979         (
1980             isBoundaryPoint,
1981             boundaryRegion,
1982             meshC,
1983             cValsPtr_(),
1984             pVals_,
1985
1986             snappedPoints,
1987             snappedPoint
1988         );
1989     }
1990     else
1991     {
1992         snappedPoint.setSize(mesh_.nPoints());
1993         snappedPoint = -1;
1994     }
1995
1996     if (debug)
1997     {
1998         Pout<< "isoSurface : shifted " << snappedPoints.size()-nCellSnaps
1999             << " vertices to intersection." << endl;
2000     }
2001
2002
2003
2004     DynamicList<point> triPoints(nCutCells_);
2005     DynamicList<label> triMeshCells(nCutCells_);
2006
2007     generateTriPoints
2008     (
2009         cValsPtr_(),
2010         pVals_,
2011
2012         meshC,
2013         mesh_.points(),
2014
2015         snappedPoints,
2016         snappedCc,
2017         snappedPoint,
2018
2019         triPoints,
2020         triMeshCells
2021     );
2022
2023     if (debug)
2024     {
2025         Pout<< "isoSurface : generated " << triMeshCells.size()
2026             << " unmerged triangles from " << triPoints.size()
2027             << " unmerged points." << endl;
2028     }
2029
2030
2031     // Merge points and compact out non-valid triangles
2032     labelList triMap;           // merged to unmerged triangle
2033     triSurface::operator=
2034     (
2035         stitchTriPoints
2036         (
2037             true,               // check for duplicate tris
2038             triPoints,
2039             triPointMergeMap_,  // unmerged to merged point
2040             triMap
2041         )
2042     );
2043
2044     if (debug)
2045     {
2046         Pout<< "isoSurface : generated " << triMap.size()
2047             << " merged triangles." << endl;
2048     }
2049
2050     meshCells_.setSize(triMap.size());
2051     forAll(triMap, i)
2052     {
2053         meshCells_[i] = triMeshCells[triMap[i]];
2054     }
2055
2056     if (debug)
2057     {
2058         Pout<< "isoSurface : checking " << size()
2059             << " triangles for validity." << endl;
2060
2061         forAll(*this, triI)
2062         {
2063             // Copied from surfaceCheck
2064             validTri(*this, triI);
2065         }
2066     }
2067
2068
2069 //if (false)
2070 {
2071     List<FixedList<label, 3> > faceEdges;
2072     labelList edgeFace0, edgeFace1;
2073     Map<labelList> edgeFacesRest;
2074
2075
2076     while (true)
2077     {
2078         // Calculate addressing
2079         calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
2080
2081         // See if any dangling triangles
2082         boolList keepTriangles;
2083         label nDangling = markDanglingTriangles
2084         (
2085             faceEdges,
2086             edgeFace0,
2087             edgeFace1,
2088             edgeFacesRest,
2089             keepTriangles
2090         );
2091
2092         if (debug)
2093         {
2094             Pout<< "isoSurface : detected " << nDangling
2095                 << " dangling triangles." << endl;
2096         }
2097
2098         if (nDangling == 0)
2099         {
2100             break;
2101         }
2102
2103         // Create face map (new to old)
2104         labelList subsetTriMap(findIndices(keepTriangles, true));
2105
2106         labelList subsetPointMap;
2107         labelList reversePointMap;
2108         triSurface::operator=
2109         (
2110             subsetMesh
2111             (
2112                 *this,
2113                 subsetTriMap,
2114                 reversePointMap,
2115                 subsetPointMap
2116             )
2117         );
2118         meshCells_ = labelField(meshCells_, subsetTriMap);
2119         inplaceRenumber(reversePointMap, triPointMergeMap_);
2120     }
2121
2122     orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
2123 }
2124
2125
2126     if (debug)
2127     {
2128         fileName stlFile = mesh_.time().path() + ".stl";
2129         Pout<< "Dumping surface to " << stlFile << endl;
2130         triSurface::write(stlFile);
2131     }
2132 }
2133
2134
2135 // ************************************************************************* //