src/sampling/sampledSurface/isoSurface/isoSurface.C

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