00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "isoSurfaceCell.H"
00027 #include <OpenFOAM/dictionary.H>
00028 #include <OpenFOAM/polyMesh.H>
00029 #include <OpenFOAM/mergePoints.H>
00030 #include <OpenFOAM/tetMatcher.H>
00031 #include <OpenFOAM/syncTools.H>
00032 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00033 #include <triSurface/faceTriangulation.H>
00034
00035
00036
00037 namespace Foam
00038 {
00039 defineTypeNameAndDebug(isoSurfaceCell, 0);
00040 }
00041
00042
00043
00044 Foam::scalar Foam::isoSurfaceCell::isoFraction
00045 (
00046 const scalar s0,
00047 const scalar s1
00048 ) const
00049 {
00050 scalar d = s1-s0;
00051
00052 if (mag(d) > VSMALL)
00053 {
00054 return (iso_-s0)/d;
00055 }
00056 else
00057 {
00058 return -1.0;
00059 }
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
00082 (
00083 const PackedBoolList& isTet,
00084 const scalarField& cellValues,
00085 const scalarField& pointValues,
00086 const label cellI
00087 ) const
00088 {
00089 const cell& cFaces = mesh_.cells()[cellI];
00090
00091 if (isTet.get(cellI) == 1)
00092 {
00093 forAll(cFaces, cFaceI)
00094 {
00095 const face& f = mesh_.faces()[cFaces[cFaceI]];
00096
00097 for (label fp = 1; fp < f.size() - 1; fp++)
00098 {
00099 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00100
00101 bool aLower = (pointValues[tri[0]] < iso_);
00102 bool bLower = (pointValues[tri[1]] < iso_);
00103 bool cLower = (pointValues[tri[2]] < iso_);
00104
00105 if (aLower == bLower && aLower == cLower)
00106 {}
00107 else
00108 {
00109 return CUT;
00110 }
00111 }
00112 }
00113 return NOTCUT;
00114 }
00115 else
00116 {
00117 bool cellLower = (cellValues[cellI] < iso_);
00118
00119
00120 bool edgeCut = false;
00121
00122 forAll(cFaces, cFaceI)
00123 {
00124 const face& f = mesh_.faces()[cFaces[cFaceI]];
00125
00126
00127 forAll(f, fp)
00128 {
00129 if ((pointValues[f[fp]] < iso_) != cellLower)
00130 {
00131 edgeCut = true;
00132 break;
00133 }
00134 }
00135
00136 if (edgeCut)
00137 {
00138 break;
00139 }
00140
00141 for (label fp = 1; fp < f.size() - 1; fp++)
00142 {
00143 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00144
00145
00146
00147
00148
00149 bool aLower = (pointValues[tri[0]] < iso_);
00150 bool bLower = (pointValues[tri[1]] < iso_);
00151 bool cLower = (pointValues[tri[2]] < iso_);
00152
00153 if (aLower == bLower && aLower == cLower)
00154 {}
00155 else
00156 {
00157 edgeCut = true;
00158 break;
00159 }
00160 }
00161
00162 if (edgeCut)
00163 {
00164 break;
00165 }
00166 }
00167
00168 if (edgeCut)
00169 {
00170
00171
00172
00173
00174 const labelList& cPoints = mesh_.cellPoints(cellI);
00175
00176 label nPyrCuts = 0;
00177
00178 forAll(cPoints, i)
00179 {
00180 if ((pointValues[cPoints[i]] < iso_) != cellLower)
00181 {
00182 nPyrCuts++;
00183 }
00184 }
00185
00186 if (nPyrCuts == cPoints.size())
00187 {
00188 return SPHERE;
00189 }
00190 else
00191 {
00192 return CUT;
00193 }
00194 }
00195 else
00196 {
00197 return NOTCUT;
00198 }
00199 }
00200 }
00201
00202
00203 void Foam::isoSurfaceCell::calcCutTypes
00204 (
00205 const PackedBoolList& isTet,
00206 const scalarField& cVals,
00207 const scalarField& pVals
00208 )
00209 {
00210 cellCutType_.setSize(mesh_.nCells());
00211 nCutCells_ = 0;
00212 forAll(mesh_.cells(), cellI)
00213 {
00214 cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
00215
00216 if (cellCutType_[cellI] == CUT)
00217 {
00218 nCutCells_++;
00219 }
00220 }
00221
00222 if (debug)
00223 {
00224 Pout<< "isoSurfaceCell : detected " << nCutCells_
00225 << " candidate cut cells." << endl;
00226 }
00227 }
00228
00229
00230
00231
00232 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
00233 (
00234 const labelledTri& tri0,
00235 const labelledTri& tri1
00236 )
00237 {
00238 labelPair common(-1, -1);
00239
00240 label fp0 = 0;
00241 label fp1 = findIndex(tri1, tri0[fp0]);
00242
00243 if (fp1 == -1)
00244 {
00245 fp0 = 1;
00246 fp1 = findIndex(tri1, tri0[fp0]);
00247 }
00248
00249 if (fp1 != -1)
00250 {
00251
00252
00253
00254 label fp0p1 = tri0.fcIndex(fp0);
00255 label fp1p1 = tri1.fcIndex(fp1);
00256 label fp1m1 = tri1.rcIndex(fp1);
00257
00258 if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
00259 {
00260 common[0] = tri0[fp0];
00261 common[1] = tri0[fp0p1];
00262 }
00263 }
00264 return common;
00265 }
00266
00267
00268
00269 Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
00270 {
00271 vector sum = vector::zero;
00272
00273 forAll(s, i)
00274 {
00275 sum += s[i].centre(s.points());
00276 }
00277 return sum/s.size();
00278 }
00279
00280
00281
00282
00283 Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
00284 (
00285 pointField& localPoints,
00286 DynamicList<labelledTri, 64>& localTris
00287 )
00288 {
00289 pointIndexHit info(false, vector::zero, localTris.size());
00290
00291 if (localTris.size() == 1)
00292 {
00293 const labelledTri& tri = localTris[0];
00294 info.setPoint(tri.centre(localPoints));
00295 info.setHit();
00296 }
00297 else if (localTris.size() == 2)
00298 {
00299
00300 const labelledTri& tri0 = localTris[0];
00301 const labelledTri& tri1 = localTris[0];
00302
00303 labelPair shared = findCommonPoints(tri0, tri1);
00304
00305 if (shared[0] != -1)
00306 {
00307 info.setPoint
00308 (
00309 0.5
00310 * (
00311 tri0.centre(localPoints)
00312 + tri1.centre(localPoints)
00313 )
00314 );
00315 info.setHit();
00316 }
00317 }
00318 else if (localTris.size())
00319 {
00320
00321 triSurface surf
00322 (
00323 localTris,
00324 geometricSurfacePatchList(0),
00325 localPoints,
00326 true
00327 );
00328 localTris.clearStorage();
00329
00330 labelList faceZone;
00331 label nZones = surf.markZones
00332 (
00333 boolList(surf.nEdges(), false),
00334 faceZone
00335 );
00336
00337 if (nZones == 1)
00338 {
00339 info.setPoint(calcCentre(surf));
00340 info.setHit();
00341 }
00342 }
00343
00344 return info;
00345 }
00346
00347
00348 void Foam::isoSurfaceCell::calcSnappedCc
00349 (
00350 const PackedBoolList& isTet,
00351 const scalarField& cVals,
00352 const scalarField& pVals,
00353
00354 DynamicList<point>& snappedPoints,
00355 labelList& snappedCc
00356 ) const
00357 {
00358 const pointField& cc = mesh_.cellCentres();
00359 const pointField& pts = mesh_.points();
00360
00361 snappedCc.setSize(mesh_.nCells());
00362 snappedCc = -1;
00363
00364
00365 DynamicList<point, 64> localPoints(64);
00366 DynamicList<labelledTri, 64> localTris(64);
00367 Map<label> pointToLocal(64);
00368
00369 forAll(mesh_.cells(), cellI)
00370 {
00371 if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
00372 {
00373 scalar cVal = cVals[cellI];
00374
00375 const cell& cFaces = mesh_.cells()[cellI];
00376
00377 localPoints.clear();
00378 localTris.clear();
00379 pointToLocal.clear();
00380
00381
00382
00383
00384 forAll(cFaces, cFaceI)
00385 {
00386 const face& f = mesh_.faces()[cFaces[cFaceI]];
00387
00388 forAll(f, fp)
00389 {
00390 label pointI = f[fp];
00391
00392 scalar s = isoFraction(cVal, pVals[pointI]);
00393
00394 if (s >= 0.0 && s <= 0.5)
00395 {
00396 if (pointToLocal.insert(pointI, localPoints.size()))
00397 {
00398 localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
00399 }
00400 }
00401 }
00402 }
00403
00404 if (localPoints.size() == 1)
00405 {
00406
00407 snappedCc[cellI] = snappedPoints.size();
00408 snappedPoints.append(localPoints[0]);
00409
00410
00411
00412
00413
00414
00415 }
00416 else if (localPoints.size() == 2)
00417 {
00418
00419 snappedCc[cellI] = snappedPoints.size();
00420 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
00421
00422
00423
00424
00425
00426
00427 }
00428 else if (localPoints.size())
00429 {
00430
00431 forAll(cFaces, cFaceI)
00432 {
00433 const face& f = mesh_.faces()[cFaces[cFaceI]];
00434
00435
00436
00437
00438
00439 for (label fp = 1; fp < f.size() - 1; fp++)
00440 {
00441 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00442
00443
00444
00445
00446
00447
00448 FixedList<scalar, 3> s(3);
00449 s[0] = isoFraction(cVal, pVals[tri[0]]);
00450 s[1] = isoFraction(cVal, pVals[tri[1]]);
00451 s[2] = isoFraction(cVal, pVals[tri[2]]);
00452
00453 if
00454 (
00455 (s[0] >= 0.0 && s[0] <= 0.5)
00456 && (s[1] >= 0.0 && s[1] <= 0.5)
00457 && (s[2] >= 0.0 && s[2] <= 0.5)
00458 )
00459 {
00460 localTris.append
00461 (
00462 labelledTri
00463 (
00464 pointToLocal[tri[0]],
00465 pointToLocal[tri[1]],
00466 pointToLocal[tri[2]],
00467 0
00468 )
00469 );
00470 }
00471 }
00472 }
00473
00474 pointField surfPoints;
00475 surfPoints.transfer(localPoints);
00476 pointIndexHit info = collapseSurface(surfPoints, localTris);
00477
00478 if (info.hit())
00479 {
00480 snappedCc[cellI] = snappedPoints.size();
00481 snappedPoints.append(info.hitPoint());
00482
00483
00484
00485
00486
00487
00488 }
00489 }
00490 }
00491 }
00492 }
00493
00494
00495
00496 void Foam::isoSurfaceCell::genPointTris
00497 (
00498 const scalarField& cellValues,
00499 const scalarField& pointValues,
00500 const label pointI,
00501 const face& f,
00502 const label cellI,
00503 DynamicList<point, 64>& localTriPoints
00504 ) const
00505 {
00506 const pointField& cc = mesh_.cellCentres();
00507 const pointField& pts = mesh_.points();
00508
00509 for (label fp = 1; fp < f.size() - 1; fp++)
00510 {
00511 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
00512
00513
00514
00515
00516
00517 label index = findIndex(tri, pointI);
00518
00519 if (index == -1)
00520 {
00521 continue;
00522 }
00523
00524
00525 label b = tri[tri.fcIndex(index)];
00526 label c = tri[tri.rcIndex(index)];
00527
00528
00529 FixedList<scalar, 3> s(3);
00530 s[0] = isoFraction(pointValues[pointI], pointValues[b]);
00531 s[1] = isoFraction(pointValues[pointI], pointValues[c]);
00532 s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
00533
00534 if
00535 (
00536 (s[0] >= 0.0 && s[0] <= 0.5)
00537 && (s[1] >= 0.0 && s[1] <= 0.5)
00538 && (s[2] >= 0.0 && s[2] <= 0.5)
00539 )
00540 {
00541 localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
00542 localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
00543 localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
00544 }
00545 }
00546 }
00547
00548
00549
00550 void Foam::isoSurfaceCell::genPointTris
00551 (
00552 const scalarField& pointValues,
00553 const label pointI,
00554 const label cellI,
00555 DynamicList<point, 64>& localTriPoints
00556 ) const
00557 {
00558 const pointField& pts = mesh_.points();
00559 const cell& cFaces = mesh_.cells()[cellI];
00560
00561 FixedList<label, 4> tet;
00562
00563 label face0 = cFaces[0];
00564 const face& f0 = mesh_.faces()[face0];
00565
00566 if (mesh_.faceOwner()[face0] != cellI)
00567 {
00568 tet[0] = f0[0];
00569 tet[1] = f0[1];
00570 tet[2] = f0[2];
00571 }
00572 else
00573 {
00574 tet[0] = f0[0];
00575 tet[1] = f0[2];
00576 tet[2] = f0[1];
00577 }
00578
00579
00580 const face& f1 = mesh_.faces()[cFaces[1]];
00581
00582 forAll(f1, fp)
00583 {
00584 label p1 = f1[fp];
00585
00586 if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
00587 {
00588 tet[3] = p1;
00589 break;
00590 }
00591 }
00592
00593
00594
00595
00596 label i0 = findIndex(tet, pointI);
00597 label i1 = tet.fcIndex(i0);
00598 label i2 = tet.fcIndex(i1);
00599 label i3 = tet.fcIndex(i2);
00600
00601
00602 FixedList<scalar, 3> s(3);
00603 s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
00604 s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
00605 s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
00606
00607 if
00608 (
00609 (s[0] >= 0.0 && s[0] <= 0.5)
00610 && (s[1] >= 0.0 && s[1] <= 0.5)
00611 && (s[2] >= 0.0 && s[2] <= 0.5)
00612 )
00613 {
00614 localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
00615 localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
00616 localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
00617 }
00618 }
00619
00620
00621 void Foam::isoSurfaceCell::calcSnappedPoint
00622 (
00623 const PackedBoolList& isBoundaryPoint,
00624 const PackedBoolList& isTet,
00625 const scalarField& cVals,
00626 const scalarField& pVals,
00627
00628 DynamicList<point>& snappedPoints,
00629 labelList& snappedPoint
00630 ) const
00631 {
00632 const point greatPoint(VGREAT, VGREAT, VGREAT);
00633 pointField collapsedPoint(mesh_.nPoints(), greatPoint);
00634
00635
00636
00637 DynamicList<point, 64> localTriPoints(100);
00638 labelHashSet localPointCells(100);
00639
00640 forAll(mesh_.pointFaces(), pointI)
00641 {
00642 if (isBoundaryPoint.get(pointI) == 1)
00643 {
00644 continue;
00645 }
00646
00647 const labelList& pFaces = mesh_.pointFaces()[pointI];
00648
00649 bool anyCut = false;
00650
00651 forAll(pFaces, i)
00652 {
00653 label faceI = pFaces[i];
00654
00655 if
00656 (
00657 cellCutType_[mesh_.faceOwner()[faceI]] == CUT
00658 || (
00659 mesh_.isInternalFace(faceI)
00660 && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
00661 )
00662 )
00663 {
00664 anyCut = true;
00665 break;
00666 }
00667 }
00668
00669 if (!anyCut)
00670 {
00671 continue;
00672 }
00673
00674
00675 localPointCells.clear();
00676 localTriPoints.clear();
00677
00678 forAll(pFaces, pFaceI)
00679 {
00680 label faceI = pFaces[pFaceI];
00681 const face& f = mesh_.faces()[faceI];
00682 label own = mesh_.faceOwner()[faceI];
00683
00684
00685 if (isTet.get(own) == 1)
00686 {
00687 if (localPointCells.insert(own))
00688 {
00689 genPointTris(pVals, pointI, own, localTriPoints);
00690 }
00691 }
00692 else
00693 {
00694 genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
00695 }
00696
00697 if (mesh_.isInternalFace(faceI))
00698 {
00699 label nei = mesh_.faceNeighbour()[faceI];
00700
00701 if (isTet.get(nei) == 1)
00702 {
00703 if (localPointCells.insert(nei))
00704 {
00705 genPointTris(pVals, pointI, nei, localTriPoints);
00706 }
00707 }
00708 else
00709 {
00710 genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
00711 }
00712 }
00713 }
00714
00715 if (localTriPoints.size() == 3)
00716 {
00717
00718 pointField points;
00719 points.transfer(localTriPoints);
00720 collapsedPoint[pointI] = sum(points)/points.size();
00721
00722
00723
00724
00725 }
00726 else if (localTriPoints.size())
00727 {
00728
00729
00730
00731 labelList triMap;
00732 labelList triPointReverseMap;
00733 triSurface surf
00734 (
00735 stitchTriPoints
00736 (
00737 false,
00738 localTriPoints,
00739 triPointReverseMap,
00740 triMap
00741 )
00742 );
00743
00744 labelList faceZone;
00745 label nZones = surf.markZones
00746 (
00747 boolList(surf.nEdges(), false),
00748 faceZone
00749 );
00750
00751 if (nZones == 1)
00752 {
00753 collapsedPoint[pointI] = calcCentre(surf);
00754
00755
00756
00757 }
00758 }
00759 }
00760
00761 syncTools::syncPointList
00762 (
00763 mesh_,
00764 collapsedPoint,
00765 minEqOp<point>(),
00766 greatPoint,
00767 true
00768 );
00769
00770 snappedPoint.setSize(mesh_.nPoints());
00771 snappedPoint = -1;
00772
00773 forAll(collapsedPoint, pointI)
00774 {
00775 if (collapsedPoint[pointI] != greatPoint)
00776 {
00777 snappedPoint[pointI] = snappedPoints.size();
00778 snappedPoints.append(collapsedPoint[pointI]);
00779 }
00780 }
00781 }
00782
00783
00784
00785
00786 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
00787 (
00788 const bool checkDuplicates,
00789 const List<point>& triPoints,
00790 labelList& triPointReverseMap,
00791 labelList& triMap
00792 ) const
00793 {
00794 label nTris = triPoints.size()/3;
00795
00796 if ((triPoints.size() % 3) != 0)
00797 {
00798 FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
00799 << "Problem: number of points " << triPoints.size()
00800 << " not a multiple of 3." << abort(FatalError);
00801 }
00802
00803 pointField newPoints;
00804 mergePoints
00805 (
00806 triPoints,
00807 mergeDistance_,
00808 false,
00809 triPointReverseMap,
00810 newPoints
00811 );
00812
00813
00814 if (debug)
00815 {
00816 pointField newNewPoints;
00817 labelList oldToNew;
00818 bool hasMerged = mergePoints
00819 (
00820 newPoints,
00821 mergeDistance_,
00822 true,
00823 oldToNew,
00824 newNewPoints
00825 );
00826
00827 if (hasMerged)
00828 {
00829 FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
00830 << "Merged points contain duplicates"
00831 << " when merging with distance " << mergeDistance_ << endl
00832 << "merged:" << newPoints.size() << " re-merged:"
00833 << newNewPoints.size()
00834 << abort(FatalError);
00835 }
00836 }
00837
00838
00839 List<labelledTri> tris;
00840 {
00841 DynamicList<labelledTri> dynTris(nTris);
00842 label rawPointI = 0;
00843 DynamicList<label> newToOldTri(nTris);
00844
00845 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
00846 {
00847 labelledTri tri
00848 (
00849 triPointReverseMap[rawPointI],
00850 triPointReverseMap[rawPointI+1],
00851 triPointReverseMap[rawPointI+2],
00852 0
00853 );
00854 rawPointI += 3;
00855
00856 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
00857 {
00858 newToOldTri.append(oldTriI);
00859 dynTris.append(tri);
00860 }
00861 }
00862
00863 triMap.transfer(newToOldTri);
00864 tris.transfer(dynTris);
00865 }
00866
00867
00868
00869
00870
00871
00872 if (checkDuplicates)
00873 {
00874 if (debug)
00875 {
00876 Pout<< "isoSurfaceCell : merged from " << nTris
00877 << " down to " << tris.size() << " triangles." << endl;
00878 }
00879
00880 pointField centres(tris.size());
00881 forAll(tris, triI)
00882 {
00883 centres[triI] = tris[triI].centre(newPoints);
00884 }
00885
00886 pointField mergedCentres;
00887 labelList oldToMerged;
00888 bool hasMerged = mergePoints
00889 (
00890 centres,
00891 mergeDistance_,
00892 false,
00893 oldToMerged,
00894 mergedCentres
00895 );
00896
00897 if (debug)
00898 {
00899 Pout<< "isoSurfaceCell : detected "
00900 << centres.size()-mergedCentres.size()
00901 << " duplicate triangles." << endl;
00902 }
00903
00904 if (hasMerged)
00905 {
00906
00907 label newTriI = 0;
00908 DynamicList<label> newToOldTri(tris.size());
00909 labelList newToMaster(mergedCentres.size(), -1);
00910 forAll(tris, triI)
00911 {
00912 label mergedI = oldToMerged[triI];
00913
00914 if (newToMaster[mergedI] == -1)
00915 {
00916 newToMaster[mergedI] = triI;
00917 newToOldTri.append(triMap[triI]);
00918 tris[newTriI++] = tris[triI];
00919 }
00920 }
00921
00922 triMap.transfer(newToOldTri);
00923 tris.setSize(newTriI);
00924 }
00925 }
00926
00927 return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
00928 }
00929
00930
00931
00932 bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
00933 {
00934
00935
00936 const labelledTri& f = surf[faceI];
00937
00938 if
00939 (
00940 (f[0] < 0) || (f[0] >= surf.points().size())
00941 || (f[1] < 0) || (f[1] >= surf.points().size())
00942 || (f[2] < 0) || (f[2] >= surf.points().size())
00943 )
00944 {
00945 WarningIn("validTri(const triSurface&, const label)")
00946 << "triangle " << faceI << " vertices " << f
00947 << " uses point indices outside point range 0.."
00948 << surf.points().size()-1 << endl;
00949
00950 return false;
00951 }
00952
00953 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
00954 {
00955 WarningIn("validTri(const triSurface&, const label)")
00956 << "triangle " << faceI
00957 << " uses non-unique vertices " << f
00958 << endl;
00959 return false;
00960 }
00961
00962
00963
00964 const labelList& fFaces = surf.faceFaces()[faceI];
00965
00966
00967
00968 forAll(fFaces, i)
00969 {
00970 label nbrFaceI = fFaces[i];
00971
00972 if (nbrFaceI <= faceI)
00973 {
00974
00975 continue;
00976 }
00977
00978 const labelledTri& nbrF = surf[nbrFaceI];
00979
00980 if
00981 (
00982 ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
00983 && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
00984 && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
00985 )
00986 {
00987 WarningIn("validTri(const triSurface&, const label)")
00988 << "triangle " << faceI << " vertices " << f
00989 << " has the same vertices as triangle " << nbrFaceI
00990 << " vertices " << nbrF
00991 << endl;
00992
00993 return false;
00994 }
00995 }
00996 return true;
00997 }
00998
00999
01000 void Foam::isoSurfaceCell::calcAddressing
01001 (
01002 const triSurface& surf,
01003 List<FixedList<label, 3> >& faceEdges,
01004 labelList& edgeFace0,
01005 labelList& edgeFace1,
01006 Map<labelList>& edgeFacesRest
01007 ) const
01008 {
01009 const pointField& points = surf.points();
01010
01011 pointField edgeCentres(3*surf.size());
01012 label edgeI = 0;
01013 forAll(surf, triI)
01014 {
01015 const labelledTri& tri = surf[triI];
01016 edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
01017 edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
01018 edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
01019 }
01020
01021 pointField mergedCentres;
01022 labelList oldToMerged;
01023 bool hasMerged = mergePoints
01024 (
01025 edgeCentres,
01026 mergeDistance_,
01027 false,
01028 oldToMerged,
01029 mergedCentres
01030 );
01031
01032 if (debug)
01033 {
01034 Pout<< "isoSurfaceCell : detected "
01035 << mergedCentres.size()
01036 << " edges on " << surf.size() << " triangles." << endl;
01037 }
01038
01039 if (!hasMerged)
01040 {
01041 if (surf.size() == 1)
01042 {
01043 faceEdges.setSize(1);
01044 faceEdges[0][0] = 0;
01045 faceEdges[0][1] = 1;
01046 faceEdges[0][2] = 2;
01047 edgeFace0.setSize(1);
01048 edgeFace0[0] = 0;
01049 edgeFace1.setSize(1);
01050 edgeFace1[0] = -1;
01051 edgeFacesRest.clear();
01052 }
01053 return;
01054 }
01055
01056
01057
01058 faceEdges.setSize(surf.size());
01059 edgeI = 0;
01060 forAll(surf, triI)
01061 {
01062 faceEdges[triI][0] = oldToMerged[edgeI++];
01063 faceEdges[triI][1] = oldToMerged[edgeI++];
01064 faceEdges[triI][2] = oldToMerged[edgeI++];
01065 }
01066
01067
01068
01069 edgeFace0.setSize(mergedCentres.size());
01070 edgeFace0 = -1;
01071 edgeFace1.setSize(mergedCentres.size());
01072 edgeFace1 = -1;
01073 edgeFacesRest.clear();
01074
01075 forAll(oldToMerged, oldEdgeI)
01076 {
01077 label triI = oldEdgeI / 3;
01078 label edgeI = oldToMerged[oldEdgeI];
01079
01080 if (edgeFace0[edgeI] == -1)
01081 {
01082 edgeFace0[edgeI] = triI;
01083 }
01084 else if (edgeFace1[edgeI] == -1)
01085 {
01086 edgeFace1[edgeI] = triI;
01087 }
01088 else
01089 {
01090
01091
01092
01093
01094
01095 Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
01096
01097 if (iter != edgeFacesRest.end())
01098 {
01099 labelList& eFaces = iter();
01100 label sz = eFaces.size();
01101 eFaces.setSize(sz+1);
01102 eFaces[sz] = triI;
01103 }
01104 else
01105 {
01106 edgeFacesRest.insert(edgeI, labelList(1, triI));
01107 }
01108 }
01109 }
01110 }
01111
01112
01113 void Foam::isoSurfaceCell::walkOrientation
01114 (
01115 const triSurface& surf,
01116 const List<FixedList<label, 3> >& faceEdges,
01117 const labelList& edgeFace0,
01118 const labelList& edgeFace1,
01119 const label seedTriI,
01120 labelList& flipState
01121 )
01122 {
01123
01124 DynamicList<label> changedFaces(surf.size());
01125
01126 changedFaces.append(seedTriI);
01127
01128 while (changedFaces.size())
01129 {
01130 DynamicList<label> newChangedFaces(changedFaces.size());
01131
01132 forAll(changedFaces, i)
01133 {
01134 label triI = changedFaces[i];
01135 const labelledTri& tri = surf[triI];
01136 const FixedList<label, 3>& fEdges = faceEdges[triI];
01137
01138 forAll(fEdges, fp)
01139 {
01140 label edgeI = fEdges[fp];
01141
01142
01143 label p0 = tri[fp];
01144 label p1 = tri[tri.fcIndex(fp)];
01145
01146 label nbrI =
01147 (
01148 edgeFace0[edgeI] != triI
01149 ? edgeFace0[edgeI]
01150 : edgeFace1[edgeI]
01151 );
01152
01153 if (nbrI != -1 && flipState[nbrI] == -1)
01154 {
01155 const labelledTri& nbrTri = surf[nbrI];
01156
01157
01158 label nbrFp = findIndex(nbrTri, p0);
01159 label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
01160
01161 bool sameOrientation = (p1 == nbrP1);
01162
01163 if (flipState[triI] == 0)
01164 {
01165 flipState[nbrI] = (sameOrientation ? 0 : 1);
01166 }
01167 else
01168 {
01169 flipState[nbrI] = (sameOrientation ? 1 : 0);
01170 }
01171 newChangedFaces.append(nbrI);
01172 }
01173 }
01174 }
01175
01176 changedFaces.transfer(newChangedFaces);
01177 }
01178 }
01179
01180
01181 void Foam::isoSurfaceCell::orientSurface
01182 (
01183 triSurface& surf,
01184 const List<FixedList<label, 3> >& faceEdges,
01185 const labelList& edgeFace0,
01186 const labelList& edgeFace1,
01187 const Map<labelList>& edgeFacesRest
01188 )
01189 {
01190
01191
01192
01193 labelList flipState(surf.size(), -1);
01194
01195 label seedTriI = 0;
01196
01197 while (true)
01198 {
01199
01200 for
01201 (
01202 ;
01203 seedTriI < surf.size() && flipState[seedTriI] != -1;
01204 seedTriI++
01205 )
01206 {}
01207
01208 if (seedTriI == surf.size())
01209 {
01210 break;
01211 }
01212
01213
01214
01215 flipState[seedTriI] = 0;
01216
01217 walkOrientation
01218 (
01219 surf,
01220 faceEdges,
01221 edgeFace0,
01222 edgeFace1,
01223 seedTriI,
01224 flipState
01225 );
01226 }
01227
01228
01229 surf.clearOut();
01230 forAll(surf, triI)
01231 {
01232 if (flipState[triI] == 1)
01233 {
01234 labelledTri tri(surf[triI]);
01235
01236 surf[triI][0] = tri[0];
01237 surf[triI][1] = tri[2];
01238 surf[triI][2] = tri[1];
01239 }
01240 else if (flipState[triI] == -1)
01241 {
01242 FatalErrorIn
01243 (
01244 "isoSurfaceCell::orientSurface(triSurface&, const label)"
01245 ) << "problem" << abort(FatalError);
01246 }
01247 }
01248 }
01249
01250
01251
01252 bool Foam::isoSurfaceCell::danglingTriangle
01253 (
01254 const FixedList<label, 3>& fEdges,
01255 const labelList& edgeFace1
01256 )
01257 {
01258 label nOpen = 0;
01259 forAll(fEdges, i)
01260 {
01261 if (edgeFace1[fEdges[i]] == -1)
01262 {
01263 nOpen++;
01264 }
01265 }
01266
01267 if (nOpen == 1 || nOpen == 2 || nOpen == 3)
01268 {
01269 return true;
01270 }
01271 else
01272 {
01273 return false;
01274 }
01275 }
01276
01277
01278
01279 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
01280 (
01281 const List<FixedList<label, 3> >& faceEdges,
01282 const labelList& edgeFace0,
01283 const labelList& edgeFace1,
01284 const Map<labelList>& edgeFacesRest,
01285 boolList& keepTriangles
01286 )
01287 {
01288 keepTriangles.setSize(faceEdges.size());
01289 keepTriangles = true;
01290
01291 label nDangling = 0;
01292
01293
01294 forAllConstIter(Map<labelList>, edgeFacesRest, iter)
01295 {
01296
01297
01298
01299 label edgeI = iter.key();
01300 const labelList& otherEdgeFaces = iter();
01301
01302
01303 if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
01304 {
01305 keepTriangles[edgeFace0[edgeI]] = false;
01306 nDangling++;
01307 }
01308 if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
01309 {
01310 keepTriangles[edgeFace1[edgeI]] = false;
01311 nDangling++;
01312 }
01313 forAll(otherEdgeFaces, i)
01314 {
01315 label triI = otherEdgeFaces[i];
01316 if (danglingTriangle(faceEdges[triI], edgeFace1))
01317 {
01318 keepTriangles[triI] = false;
01319 nDangling++;
01320 }
01321 }
01322 }
01323 return nDangling;
01324 }
01325
01326
01327 Foam::triSurface Foam::isoSurfaceCell::subsetMesh
01328 (
01329 const triSurface& s,
01330 const labelList& newToOldFaces,
01331 labelList& oldToNewPoints,
01332 labelList& newToOldPoints
01333 )
01334 {
01335 const boolList include
01336 (
01337 createWithValues<boolList>
01338 (
01339 s.size(),
01340 false,
01341 newToOldFaces,
01342 true
01343 )
01344 );
01345
01346 newToOldPoints.setSize(s.points().size());
01347 oldToNewPoints.setSize(s.points().size());
01348 oldToNewPoints = -1;
01349 {
01350 label pointI = 0;
01351
01352 forAll(include, oldFacei)
01353 {
01354 if (include[oldFacei])
01355 {
01356
01357 const labelledTri& tri = s[oldFacei];
01358
01359 forAll(tri, fp)
01360 {
01361 label oldPointI = tri[fp];
01362
01363 if (oldToNewPoints[oldPointI] == -1)
01364 {
01365 oldToNewPoints[oldPointI] = pointI;
01366 newToOldPoints[pointI++] = oldPointI;
01367 }
01368 }
01369 }
01370 }
01371 newToOldPoints.setSize(pointI);
01372 }
01373
01374
01375 pointField newPoints(newToOldPoints.size());
01376 forAll(newToOldPoints, i)
01377 {
01378 newPoints[i] = s.points()[newToOldPoints[i]];
01379 }
01380
01381 List<labelledTri> newTriangles(newToOldFaces.size());
01382
01383 forAll(newToOldFaces, i)
01384 {
01385
01386 const labelledTri& tri = s[newToOldFaces[i]];
01387
01388 newTriangles[i][0] = oldToNewPoints[tri[0]];
01389 newTriangles[i][1] = oldToNewPoints[tri[1]];
01390 newTriangles[i][2] = oldToNewPoints[tri[2]];
01391 newTriangles[i].region() = tri.region();
01392 }
01393
01394
01395 return triSurface(newTriangles, s.patches(), newPoints, true);
01396 }
01397
01398
01399
01400
01401 Foam::isoSurfaceCell::isoSurfaceCell
01402 (
01403 const polyMesh& mesh,
01404 const scalarField& cVals,
01405 const scalarField& pVals,
01406 const scalar iso,
01407 const bool regularise,
01408 const scalar mergeTol
01409 )
01410 :
01411 mesh_(mesh),
01412 iso_(iso),
01413 mergeDistance_(mergeTol*mesh.bounds().mag())
01414 {
01415
01416 PackedBoolList isTet(mesh_.nCells());
01417 {
01418 tetMatcher tet;
01419
01420 forAll(isTet, cellI)
01421 {
01422 if (tet.isA(mesh_, cellI))
01423 {
01424 isTet.set(cellI, 1);
01425 }
01426 }
01427 }
01428
01429
01430
01431 PackedBoolList isBoundaryPoint(mesh_.nPoints());
01432 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01433 forAll(patches, patchI)
01434 {
01435 const polyPatch& pp = patches[patchI];
01436
01437 if (!pp.coupled())
01438 {
01439 label faceI = pp.start();
01440 forAll(pp, i)
01441 {
01442 const face& f = mesh_.faces()[faceI++];
01443
01444 forAll(f, fp)
01445 {
01446 isBoundaryPoint.set(f[fp], 1);
01447 }
01448 }
01449 }
01450 }
01451
01452
01453
01454 calcCutTypes(isTet, cVals, pVals);
01455
01456 DynamicList<point> snappedPoints(nCutCells_);
01457
01458
01459 labelList snappedCc;
01460 if (regularise)
01461 {
01462 calcSnappedCc
01463 (
01464 isTet,
01465 cVals,
01466 pVals,
01467 snappedPoints,
01468 snappedCc
01469 );
01470 }
01471 else
01472 {
01473 snappedCc.setSize(mesh_.nCells());
01474 snappedCc = -1;
01475 }
01476
01477 if (debug)
01478 {
01479 Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
01480 << " cell centres to intersection." << endl;
01481 }
01482
01483 snappedPoints.shrink();
01484 label nCellSnaps = snappedPoints.size();
01485
01486
01487 labelList snappedPoint;
01488 if (regularise)
01489 {
01490 calcSnappedPoint
01491 (
01492 isBoundaryPoint,
01493 isTet,
01494 cVals,
01495 pVals,
01496 snappedPoints,
01497 snappedPoint
01498 );
01499 }
01500 else
01501 {
01502 snappedPoint.setSize(mesh_.nPoints());
01503 snappedPoint = -1;
01504 }
01505
01506 if (debug)
01507 {
01508 Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
01509 << " vertices to intersection." << endl;
01510 }
01511
01512
01513
01514 DynamicList<point> triPoints(nCutCells_);
01515 DynamicList<label> triMeshCells(nCutCells_);
01516
01517 generateTriPoints
01518 (
01519 cVals,
01520 pVals,
01521
01522 mesh_.cellCentres(),
01523 mesh_.points(),
01524
01525 snappedPoints,
01526 snappedCc,
01527 snappedPoint,
01528
01529 triPoints,
01530 triMeshCells
01531 );
01532
01533 if (debug)
01534 {
01535 Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
01536 << " unmerged triangles." << endl;
01537 }
01538
01539
01540 labelList triMap;
01541 triSurface::operator=
01542 (
01543 stitchTriPoints
01544 (
01545 true,
01546 triPoints,
01547 triPointMergeMap_,
01548 triMap
01549 )
01550 );
01551
01552 if (debug)
01553 {
01554 Pout<< "isoSurfaceCell : generated " << triMap.size()
01555 << " merged triangles." << endl;
01556 }
01557
01558 meshCells_.setSize(triMap.size());
01559 forAll(triMap, i)
01560 {
01561 meshCells_[i] = triMeshCells[triMap[i]];
01562 }
01563
01564 if (debug)
01565 {
01566 Pout<< "isoSurfaceCell : checking " << size()
01567 << " triangles for validity." << endl;
01568
01569 forAll(*this, triI)
01570 {
01571
01572 validTri(*this, triI);
01573 }
01574 }
01575
01576
01577
01578 {
01579 List<FixedList<label, 3> > faceEdges;
01580 labelList edgeFace0, edgeFace1;
01581 Map<labelList> edgeFacesRest;
01582
01583
01584 while (true)
01585 {
01586
01587 calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
01588
01589
01590 boolList keepTriangles;
01591 label nDangling = markDanglingTriangles
01592 (
01593 faceEdges,
01594 edgeFace0,
01595 edgeFace1,
01596 edgeFacesRest,
01597 keepTriangles
01598 );
01599
01600 if (debug)
01601 {
01602 Pout<< "isoSurfaceCell : detected " << nDangling
01603 << " dangling triangles." << endl;
01604 }
01605
01606 if (nDangling == 0)
01607 {
01608 break;
01609 }
01610
01611
01612 labelList subsetTriMap(findIndices(keepTriangles, true));
01613
01614 labelList subsetPointMap;
01615 labelList reversePointMap;
01616 triSurface::operator=
01617 (
01618 subsetMesh
01619 (
01620 *this,
01621 subsetTriMap,
01622 reversePointMap,
01623 subsetPointMap
01624 )
01625 );
01626 meshCells_ = labelField(meshCells_, subsetTriMap);
01627 inplaceRenumber(reversePointMap, triPointMergeMap_);
01628 }
01629
01630 orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
01631 }
01632
01633
01634 }
01635
01636
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01755
01756