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