63 const bool aLower = (pointValues[a] < isoValue);
64 const bool bLower = (pointValues[
b] < isoValue);
65 const bool cLower = (pointValues[
c] < isoValue);
67 return !(aLower == bLower && aLower == cLower);
74 Foam::scalar Foam::isoSurfaceCell::isoFraction
80 const scalar d = s1-s0;
91 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
99 if (ignoreCells_.test(celli))
104 const cell& cFaces = mesh_.cells()[celli];
106 if (isTet.test(celli))
108 for (
const label facei : cFaces)
110 const face&
f = mesh_.faces()[facei];
112 for (label fp = 1; fp <
f.size() - 1; ++fp)
114 if (
isTriCut(
f[0],
f[fp],
f[
f.fcIndex(fp)], pointValues, iso_))
124 const bool cellLower = (cellValues[celli] < iso_);
127 bool edgeCut =
false;
129 for (
const label facei : cFaces)
131 const face&
f = mesh_.faces()[facei];
134 for (
const label pointi :
f)
136 if (cellLower != (pointValues[pointi] < iso_))
151 (mesh_.tetBasePtIs()[facei] < 0 ? 0 : mesh_.tetBasePtIs()[facei]);
153 label fp =
f.fcIndex(fp0);
154 for (label i = 2; i <
f.size(); ++i)
156 const label nextFp =
f.fcIndex(fp);
158 if (
isTriCut(
f[fp0],
f[fp],
f[nextFp], pointValues, iso_))
180 const labelList& cPoints = mesh_.cellPoints(celli);
184 for (
const label pointi : cPoints)
186 if ((pointValues[pointi] < iso_) != cellLower)
192 if (nPyrCuts == cPoints.size())
207 void Foam::isoSurfaceCell::calcCutTypes
214 cellCutType_.setSize(mesh_.nCells());
219 (void)mesh_.tetBasePtIs();
221 forAll(cellCutType_, celli)
223 cellCutType_[celli] = calcCutType(isTet, cVals, pVals, celli);
225 if (cellCutType_[celli] == CUT)
233 Pout<<
"isoSurfaceCell : candidate cut cells "
234 << nCutCells_ <<
" / " << mesh_.nCells() <<
endl;
241 const labelledTri& tri0,
242 const labelledTri& tri1
248 label fp1 = tri1.find(tri0[fp0]);
253 fp1 = tri1.find(tri0[fp0]);
261 label fp0p1 = tri0.fcIndex(fp0);
262 label fp1p1 = tri1.fcIndex(fp1);
263 label fp1m1 = tri1.rcIndex(fp1);
265 if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
267 common[0] = tri0[fp0];
268 common[1] = tri0[fp0p1];
275 Foam::point Foam::isoSurfaceCell::calcCentre(
const triSurface&
s)
281 sum +=
s[i].centre(
s.points());
291 DynamicList<labelledTri, 64>& localTris
296 if (localTris.size() == 1)
298 const labelledTri& tri = localTris[0];
299 info.setPoint(tri.centre(localPoints));
302 else if (localTris.size() == 2)
305 const labelledTri& tri0 = localTris[0];
306 const labelledTri& tri1 = localTris[1];
308 labelPair shared = findCommonPoints(tri0, tri1);
312 const vector n0 = tri0.areaNormal(localPoints);
313 const vector n1 = tri1.areaNormal(localPoints);
320 mag(n0) <= ROOTVSMALL
321 ||
mag(n1) <= ROOTVSMALL
329 tri0.centre(localPoints)
330 + tri1.centre(localPoints)
337 else if (localTris.size())
347 localTris.clearStorage();
350 label
nZones = surf.markZones
359 scalar minCos = GREAT;
360 const vector& n0 = surf.faceNormals()[0];
361 for (label i = 1; i < surf.size(); ++i)
363 scalar cosAngle = (n0 & surf.faceNormals()[i]);
364 if (cosAngle < minCos)
372 info.setPoint(calcCentre(surf));
382 void Foam::isoSurfaceCell::calcSnappedCc
388 DynamicList<point>& snappedPoints,
395 snappedCc.setSize(mesh_.nCells());
399 DynamicList<point, 64> localPoints(64);
400 DynamicList<labelledTri, 64> localTris(64);
401 Map<label> pointToLocal(64);
403 forAll(mesh_.cells(), celli)
405 if (cellCutType_[celli] == CUT && !isTet.test(celli))
407 const scalar cVal = cVals[celli];
409 const cell& cFaces = mesh_.cells()[celli];
413 pointToLocal.clear();
418 for (
const label facei : cFaces)
420 const face&
f = mesh_.faces()[facei];
422 for (
const label pointi :
f)
424 scalar
s = isoFraction(cVal, pVals[pointi]);
426 if (
s >= 0.0 &&
s <= 0.5)
428 if (pointToLocal.insert(pointi, localPoints.size()))
430 localPoints.append((1.0-
s)*cc[celli]+
s*pts[pointi]);
436 if (localPoints.size() == 1)
439 snappedCc[celli] = snappedPoints.size();
440 snappedPoints.append(localPoints[0]);
448 else if (localPoints.size() == 2)
451 snappedCc[celli] = snappedPoints.size();
452 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
460 else if (localPoints.size())
463 for (
const label facei : cFaces)
465 const face&
f = mesh_.faces()[facei];
471 const label fp0 = mesh_.tetBasePtIs()[facei];
472 label fp =
f.fcIndex(fp0);
473 for (label i = 2; i <
f.size(); ++i)
475 label nextFp =
f.fcIndex(fp);
479 FixedList<scalar, 3>
s(3);
480 s[0] = isoFraction(cVal, pVals[tri[0]]);
481 s[1] = isoFraction(cVal, pVals[tri[1]]);
482 s[2] = isoFraction(cVal, pVals[tri[2]]);
486 (
s[0] >= 0.0 &&
s[0] <= 0.5)
487 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
488 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
493 (mesh_.faceOwner()[facei] == celli)
494 == (cVal >= pVals[tri[0]])
501 pointToLocal[tri[1]],
502 pointToLocal[tri[0]],
503 pointToLocal[tri[2]],
514 pointToLocal[tri[0]],
515 pointToLocal[tri[1]],
516 pointToLocal[tri[2]],
528 surfPoints.transfer(localPoints);
538 snappedCc[celli] = snappedPoints.size();
539 snappedPoints.append(info.hitPoint());
553 void Foam::isoSurfaceCell::genPointTris
560 DynamicList<point, 64>& localTriPoints
565 const face&
f = mesh_.faces()[facei];
567 const label fp0 = mesh_.tetBasePtIs()[facei];
568 label fp =
f.fcIndex(fp0);
569 for (label i = 2; i <
f.size(); ++i)
571 label nextFp =
f.fcIndex(fp);
574 label index = tri.find(pointi);
582 label
b = tri[tri.fcIndex(index)];
583 label
c = tri[tri.rcIndex(index)];
586 FixedList<scalar, 3>
s(3);
587 s[0] = isoFraction(pointValues[pointi], pointValues[
b]);
588 s[1] = isoFraction(pointValues[pointi], pointValues[
c]);
589 s[2] = isoFraction(pointValues[pointi], cellValues[celli]);
593 (
s[0] >= 0.0 &&
s[0] <= 0.5)
594 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
595 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
598 point p0 = (1.0-
s[0])*pts[pointi] +
s[0]*pts[
b];
599 point p1 = (1.0-
s[1])*pts[pointi] +
s[1]*pts[
c];
600 point p2 = (1.0-
s[2])*pts[pointi] +
s[2]*cc[celli];
604 (mesh_.faceOwner()[facei] == celli)
605 == (pointValues[pointi] > cellValues[celli])
608 localTriPoints.append(
p0);
609 localTriPoints.append(p1);
610 localTriPoints.append(p2);
614 localTriPoints.append(p1);
615 localTriPoints.append(
p0);
616 localTriPoints.append(p2);
625 void Foam::isoSurfaceCell::genPointTris
631 DynamicList<point, 64>& localTriPoints
635 const cell& cFaces = mesh_.cells()[celli];
639 const face&
f = mesh_.faces()[facei];
643 for (
const label cfacei : cFaces)
645 const face& f1 = mesh_.faces()[cfacei];
646 for (
const label p1 : f1)
662 label index =
f.find(pointi);
663 label
b =
f[
f.fcIndex(index)];
664 label
c =
f[
f.rcIndex(index)];
670 FixedList<scalar, 3>
s(3);
671 s[0] = isoFraction(pointValues[pointi], pointValues[
b]);
672 s[1] = isoFraction(pointValues[pointi], pointValues[
c]);
673 s[2] = isoFraction(pointValues[pointi], pointValues[ccPointi]);
677 (
s[0] >= 0.0 &&
s[0] <= 0.5)
678 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
679 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
682 point p0 = (1.0-
s[0])*pts[pointi] +
s[0]*pts[
b];
683 point p1 = (1.0-
s[1])*pts[pointi] +
s[1]*pts[
c];
684 point p2 = (1.0-
s[2])*pts[pointi] +
s[2]*pts[ccPointi];
686 if (mesh_.faceOwner()[facei] != celli)
688 localTriPoints.append(
p0);
689 localTriPoints.append(p1);
690 localTriPoints.append(p2);
694 localTriPoints.append(p1);
695 localTriPoints.append(
p0);
696 localTriPoints.append(p2);
702 void Foam::isoSurfaceCell::calcSnappedPoint
708 DynamicList<point>& snappedPoints,
714 bitSet isBoundaryPoint(mesh_.nPoints());
715 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
716 for (
const polyPatch& pp :
patches)
720 label facei = pp.
start();
723 const face&
f = mesh_.faces()[facei++];
725 isBoundaryPoint.
set(
f);
731 const point greatPoint(GREAT, GREAT, GREAT);
733 pointField collapsedPoint(mesh_.nPoints(), greatPoint);
737 DynamicList<point, 64> localTriPoints(100);
740 forAll(mesh_.pointFaces(), pointi)
742 if (isBoundaryPoint.test(pointi))
751 for (
const label facei :
pFaces)
755 cellCutType_[mesh_.faceOwner()[facei]] == CUT
757 mesh_.isInternalFace(facei)
758 && cellCutType_[mesh_.faceNeighbour()[facei]] == CUT
775 localPointCells.clear();
776 localTriPoints.clear();
778 for (
const label facei :
pFaces)
780 const label own = mesh_.faceOwner()[facei];
786 if (localPointCells.insert(own))
788 genPointTris(pVals, pointi, facei, own, localTriPoints);
804 if (mesh_.isInternalFace(facei))
806 const label nei = mesh_.faceNeighbour()[facei];
810 if (localPointCells.insert(nei))
812 genPointTris(pVals, pointi, facei, nei, localTriPoints);
830 if (localTriPoints.size() == 3)
834 points.transfer(localTriPoints);
841 else if (localTriPoints.size())
860 label
nZones = surf.markZones
869 scalar minCos = GREAT;
870 const vector& n0 = surf.faceNormals()[0];
871 for (label i = 1; i < surf.size(); ++i)
873 const vector&
n = surf.faceNormals()[i];
874 scalar cosAngle = (n0 &
n);
875 if (cosAngle < minCos)
882 collapsedPoint[pointi] = calcCentre(surf);
892 minMagSqrEqOp<point>(),
896 snappedPoint.setSize(mesh_.nPoints());
899 forAll(collapsedPoint, pointi)
903 if (
magSqr(collapsedPoint[pointi]) < 0.5*
magSqr(greatPoint))
905 snappedPoint[pointi] = snappedPoints.size();
906 snappedPoints.append(collapsedPoint[pointi]);
914 const bool checkDuplicates,
915 const List<point>& triPoints,
920 label nTris = triPoints.size()/3;
922 if ((triPoints.size() % 3) != 0)
925 <<
"Problem: number of points " << triPoints.size()
942 Pout<<
"isoSurfaceCell : merged from " << triPoints.size()
943 <<
" points down to " << newPoints.size() <<
endl;
959 <<
"Merged points contain duplicates"
960 <<
" when merging with distance " << mergeDistance_ <<
endl
961 <<
"merged:" << newPoints.size() <<
" re-merged:"
962 << newNewPoints.size()
968 List<labelledTri> tris;
970 DynamicList<labelledTri> dynTris(nTris);
972 DynamicList<label> newToOldTri(nTris);
974 for (label oldTriI = 0; oldTriI < nTris; ++oldTriI)
978 triPointReverseMap[rawPointi],
979 triPointReverseMap[rawPointi+1],
980 triPointReverseMap[rawPointi+2],
983 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
985 newToOldTri.append(oldTriI);
992 triMap.transfer(newToOldTri);
993 tris.transfer(dynTris);
1001 if (checkDuplicates)
1005 Pout<<
"isoSurfaceCell : merged from " << nTris
1006 <<
" down to " << tris.size() <<
" triangles." <<
endl;
1012 centres[triI] = tris[triI].centre(newPoints);
1028 Pout<<
"isoSurfaceCell : detected "
1029 << centres.size()-mergedCentres.size()
1030 <<
" duplicate triangles." <<
endl;
1037 DynamicList<label> newToOldTri(tris.size());
1038 labelList newToMaster(mergedCentres.size(), -1);
1041 label mergedI = oldToMerged[triI];
1043 if (newToMaster[mergedI] == -1)
1045 newToMaster[mergedI] = triI;
1046 newToOldTri.append(triMap[triI]);
1047 tris[newTriI++] = tris[triI];
1051 triMap.transfer(newToOldTri);
1052 tris.setSize(newTriI);
1060 void Foam::isoSurfaceCell::calcAddressing
1062 const triSurface& surf,
1063 List<FixedList<label, 3>>& faceEdges,
1066 Map<labelList>& edgeFacesRest
1075 const labelledTri& tri = surf[triI];
1076 edgeCentres[edgeI++] = 0.5*(
points[tri[0]]+
points[tri[1]]);
1077 edgeCentres[edgeI++] = 0.5*(
points[tri[1]]+
points[tri[2]]);
1078 edgeCentres[edgeI++] = 0.5*(
points[tri[2]]+
points[tri[0]]);
1094 Pout<<
"isoSurfaceCell : detected "
1095 << mergedCentres.size()
1096 <<
" edges on " << surf.size() <<
" triangles." <<
endl;
1106 faceEdges.setSize(surf.size());
1110 faceEdges[triI][0] = oldToMerged[edgeI++];
1111 faceEdges[triI][1] = oldToMerged[edgeI++];
1112 faceEdges[triI][2] = oldToMerged[edgeI++];
1117 edgeFace0.setSize(mergedCentres.size());
1119 edgeFace1.setSize(mergedCentres.size());
1121 edgeFacesRest.clear();
1123 forAll(oldToMerged, oldEdgeI)
1125 label triI = oldEdgeI / 3;
1126 label edgeI = oldToMerged[oldEdgeI];
1128 if (edgeFace0[edgeI] == -1)
1130 edgeFace0[edgeI] = triI;
1132 else if (edgeFace1[edgeI] == -1)
1134 edgeFace1[edgeI] = triI;
1145 if (iter != edgeFacesRest.end())
1148 label sz = eFaces.size();
1149 eFaces.setSize(sz+1);
1154 edgeFacesRest.insert(edgeI,
labelList(1, triI));
1161 bool Foam::isoSurfaceCell::danglingTriangle
1163 const FixedList<label, 3>& fEdges,
1168 for (
const label edgei : fEdges)
1170 if (edgeFace1[edgei] == -1)
1176 return (nOpen == 1 || nOpen == 2 || nOpen == 3);
1180 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1182 const List<FixedList<label, 3>>& faceEdges,
1185 const Map<labelList>& edgeFacesRest,
1189 keepTriangles.setSize(faceEdges.size());
1190 keepTriangles =
true;
1192 label nDangling = 0;
1200 const label edgeI = iter.key();
1201 const labelList& otherEdgeFaces = iter.val();
1204 if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1206 keepTriangles[edgeFace0[edgeI]] =
false;
1209 if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1211 keepTriangles[edgeFace1[edgeI]] =
false;
1214 for (
const label triI : otherEdgeFaces)
1216 if (danglingTriangle(faceEdges[triI], edgeFace1))
1218 keepTriangles[triI] =
false;
1229 const triSurface&
s,
1237 ListOps::createWithValue<bool>(
s.size(), newToOldFaces,
true,
false)
1240 newToOldPoints.setSize(
s.points().size());
1241 oldToNewPoints.setSize(
s.points().size());
1242 oldToNewPoints = -1;
1246 forAll(include, oldFacei)
1248 if (include[oldFacei])
1251 for (
const label oldPointi :
s[oldFacei])
1253 if (oldToNewPoints[oldPointi] == -1)
1255 oldToNewPoints[oldPointi] = pointi;
1256 newToOldPoints[pointi++] = oldPointi;
1261 newToOldPoints.setSize(pointi);
1266 forAll(newToOldPoints, i)
1268 newPoints[i] =
s.points()[newToOldPoints[i]];
1271 List<labelledTri> newTriangles(newToOldFaces.size());
1276 const labelledTri& tri =
s[newToOldFaces[i]];
1278 newTriangles[i][0] = oldToNewPoints[tri[0]];
1279 newTriangles[i][1] = oldToNewPoints[tri[1]];
1280 newTriangles[i][2] = oldToNewPoints[tri[2]];
1281 newTriangles[i].region() = tri.region();
1285 return triSurface(newTriangles,
s.patches(), newPoints,
true);
1299 const scalar mergeTol,
1300 const bitSet& ignoreCells
1309 (filter != filterType::NONE),
1323 const bool regularise,
1325 const scalar mergeTol,
1326 const bitSet& ignoreCells
1332 pVals_(pointValues),
1333 ignoreCells_(ignoreCells),
1338 Pout<<
"isoSurfaceCell::"
1339 <<
" cell min/max : "
1340 <<
min(cVals_) <<
" / "
1341 <<
max(cVals_) <<
nl
1342 <<
" point min/max : "
1343 <<
min(pVals_) <<
" / "
1344 <<
max(pVals_) <<
nl
1345 <<
" isoValue : " << iso <<
nl
1346 <<
" filter : " <<
Switch(regularise) <<
nl
1347 <<
" mergeTol : " << mergeTol <<
nl
1349 <<
" mergeDistance : " << mergeDistance_ <<
nl
1350 <<
" ignoreCells : " << ignoreCells_.count()
1351 <<
" / " << cVals_.size() <<
nl
1356 bitSet isTet(mesh_.nCells());
1362 if (tet.
isA(mesh_, celli))
1371 calcCutTypes(isTet, cellValues, pointValues);
1375 const fvMesh& fvm = dynamicCast<const fvMesh&>(
mesh);
1381 "isoSurfaceCell.cutType",
1394 forAll(cellCutType_, celli)
1396 debugFld[celli] = cellCutType_[celli];
1399 Pout<<
"Writing cut types:"
1400 << debugField.objectPath() <<
endl;
1423 snappedCc.
setSize(mesh_.nCells());
1429 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.size()
1430 <<
" cell centres to intersection." <<
endl;
1434 label nCellSnaps = snappedPoints.size();
1451 snappedPoint.
setSize(mesh_.nPoints());
1457 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
1458 <<
" vertices to intersection." <<
endl;
1474 mesh_.cellCentres(),
1487 Pout<<
"isoSurfaceCell : generated " << triMeshCells.size()
1488 <<
" unmerged triangles." <<
endl;
1498 if (bounds_.valid())
1500 isoSurface::trimToBox
1506 interpolatedPoints_,
1507 interpolatedOldPoints_,
1508 interpolationWeights_
1510 triMeshCells =
labelField(triMeshCells, trimTriMap);
1516 tmpsurf = stitchTriPoints
1526 Pout<<
"isoSurfaceCell : generated " << triMap.size()
1527 <<
" merged triangles." <<
endl;
1530 if (bounds_.valid())
1536 labelList newTriPointMergeMap(nOldPoints, -1);
1537 forAll(trimTriPointMap, trimPointI)
1539 label oldPointI = trimTriPointMap[trimPointI];
1542 label pointI = triPointMergeMap_[trimPointI];
1545 newTriPointMergeMap[oldPointI] = pointI;
1549 triPointMergeMap_.
transfer(newTriPointMergeMap);
1552 meshCells_.setSize(triMap.size());
1555 meshCells_[i] = triMeshCells[triMap[i]];
1562 Pout<<
"isoSurfaceCell : checking " << tmpsurf.size()
1563 <<
" triangles for validity." <<
endl;
1593 label nDangling = markDanglingTriangles
1604 Pout<<
"isoSurfaceCell : detected " << nDangling
1605 <<
" dangling triangles." <<
endl;
1618 tmpsurf = subsetMesh
1625 meshCells_ =
labelField(meshCells_, subsetTriMap);
1645 this->MeshStorage::transfer(updated);