57 Foam::scalar Foam::isoSurfaceCell::isoFraction
63 const scalar d = s1-s0;
76 const labelledTri& tri0,
77 const labelledTri& tri1
83 label fp1 = tri1.find(tri0[fp0]);
88 fp1 = tri1.find(tri0[fp0]);
96 label fp0p1 = tri0.fcIndex(fp0);
97 label fp1p1 = tri1.fcIndex(fp1);
98 label fp1m1 = tri1.rcIndex(fp1);
100 if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
102 common[0] = tri0[fp0];
103 common[1] = tri0[fp0p1];
110 Foam::point Foam::isoSurfaceCell::calcCentre(
const triSurface&
s)
116 sum +=
s[i].centre(
s.points());
126 DynamicList<labelledTri, 64>& localTris
131 if (localTris.size() == 1)
133 const labelledTri& tri = localTris[0];
134 info.setPoint(tri.centre(localPoints));
137 else if (localTris.size() == 2)
140 const labelledTri& tri0 = localTris[0];
141 const labelledTri& tri1 = localTris[1];
143 labelPair shared = findCommonPoints(tri0, tri1);
147 const vector n0 = tri0.areaNormal(localPoints);
148 const vector n1 = tri1.areaNormal(localPoints);
155 mag(n0) <= ROOTVSMALL
156 ||
mag(n1) <= ROOTVSMALL
164 tri0.centre(localPoints)
165 + tri1.centre(localPoints)
172 else if (localTris.size())
182 localTris.clearStorage();
185 label
nZones = surf.markZones
194 scalar minCos = GREAT;
195 const vector& n0 = surf.faceNormals()[0];
196 for (label i = 1; i < surf.size(); ++i)
198 scalar cosAngle = (n0 & surf.faceNormals()[i]);
199 if (cosAngle < minCos)
207 info.setPoint(calcCentre(surf));
217 void Foam::isoSurfaceCell::calcSnappedCc
223 DynamicList<point>& snappedPoints,
230 snappedCc.setSize(mesh_.nCells());
234 DynamicList<point, 64> localPoints(64);
235 DynamicList<labelledTri, 64> localTris(64);
236 Map<label> pointToLocal(64);
238 forAll(mesh_.cells(), celli)
240 if (cellCutType_[celli] == cutType::CUT)
242 const scalar cVal = cVals[celli];
244 const cell& cFaces = mesh_.cells()[celli];
248 pointToLocal.clear();
253 for (
const label facei : cFaces)
255 const face&
f = mesh_.faces()[facei];
257 for (
const label pointi :
f)
259 scalar
s = isoFraction(cVal, pVals[pointi]);
261 if (
s >= 0.0 &&
s <= 0.5)
263 if (pointToLocal.insert(pointi, localPoints.size()))
265 localPoints.append((1.0-
s)*cc[celli]+
s*pts[pointi]);
271 if (localPoints.size() == 1)
274 snappedCc[celli] = snappedPoints.size();
275 snappedPoints.append(localPoints[0]);
283 else if (localPoints.size() == 2)
286 snappedCc[celli] = snappedPoints.size();
287 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
295 else if (localPoints.size())
298 for (
const label facei : cFaces)
300 const face&
f = mesh_.faces()[facei];
306 const label fp0 = mesh_.tetBasePtIs()[facei];
307 label fp =
f.fcIndex(fp0);
308 for (label i = 2; i <
f.size(); ++i)
310 label nextFp =
f.fcIndex(fp);
314 FixedList<scalar, 3>
s(3);
315 s[0] = isoFraction(cVal, pVals[tri[0]]);
316 s[1] = isoFraction(cVal, pVals[tri[1]]);
317 s[2] = isoFraction(cVal, pVals[tri[2]]);
321 (
s[0] >= 0.0 &&
s[0] <= 0.5)
322 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
323 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
328 (mesh_.faceOwner()[facei] == celli)
329 == (cVal >= pVals[tri[0]])
336 pointToLocal[tri[1]],
337 pointToLocal[tri[0]],
338 pointToLocal[tri[2]],
349 pointToLocal[tri[0]],
350 pointToLocal[tri[1]],
351 pointToLocal[tri[2]],
363 surfPoints.transfer(localPoints);
373 snappedCc[celli] = snappedPoints.size();
374 snappedPoints.append(info.hitPoint());
388 void Foam::isoSurfaceCell::genPointTris
395 DynamicList<point, 64>& localTriPoints
400 const face&
f = mesh_.faces()[facei];
402 const label fp0 = mesh_.tetBasePtIs()[facei];
403 label fp =
f.fcIndex(fp0);
404 for (label i = 2; i <
f.size(); ++i)
406 label nextFp =
f.fcIndex(fp);
409 label index = tri.find(pointi);
417 label
b = tri[tri.fcIndex(index)];
418 label
c = tri[tri.rcIndex(index)];
421 FixedList<scalar, 3>
s(3);
422 s[0] = isoFraction(pointValues[pointi], pointValues[
b]);
423 s[1] = isoFraction(pointValues[pointi], pointValues[
c]);
424 s[2] = isoFraction(pointValues[pointi], cellValues[celli]);
428 (
s[0] >= 0.0 &&
s[0] <= 0.5)
429 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
430 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
433 point p0 = (1.0-
s[0])*pts[pointi] +
s[0]*pts[
b];
434 point p1 = (1.0-
s[1])*pts[pointi] +
s[1]*pts[
c];
435 point p2 = (1.0-
s[2])*pts[pointi] +
s[2]*cc[celli];
439 (mesh_.faceOwner()[facei] == celli)
440 == (pointValues[pointi] > cellValues[celli])
443 localTriPoints.append(
p0);
444 localTriPoints.append(p1);
445 localTriPoints.append(p2);
449 localTriPoints.append(p1);
450 localTriPoints.append(
p0);
451 localTriPoints.append(p2);
460 void Foam::isoSurfaceCell::genPointTris
466 DynamicList<point, 64>& localTriPoints
470 const cell& cFaces = mesh_.cells()[celli];
474 const face&
f = mesh_.faces()[facei];
478 for (
const label cfacei : cFaces)
480 const face& f1 = mesh_.faces()[cfacei];
481 for (
const label p1 : f1)
497 label index =
f.find(pointi);
498 label
b =
f[
f.fcIndex(index)];
499 label
c =
f[
f.rcIndex(index)];
505 FixedList<scalar, 3>
s(3);
506 s[0] = isoFraction(pointValues[pointi], pointValues[
b]);
507 s[1] = isoFraction(pointValues[pointi], pointValues[
c]);
508 s[2] = isoFraction(pointValues[pointi], pointValues[ccPointi]);
512 (
s[0] >= 0.0 &&
s[0] <= 0.5)
513 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
514 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
517 point p0 = (1.0-
s[0])*pts[pointi] +
s[0]*pts[
b];
518 point p1 = (1.0-
s[1])*pts[pointi] +
s[1]*pts[
c];
519 point p2 = (1.0-
s[2])*pts[pointi] +
s[2]*pts[ccPointi];
521 if (mesh_.faceOwner()[facei] != celli)
523 localTriPoints.append(
p0);
524 localTriPoints.append(p1);
525 localTriPoints.append(p2);
529 localTriPoints.append(p1);
530 localTriPoints.append(
p0);
531 localTriPoints.append(p2);
537 void Foam::isoSurfaceCell::calcSnappedPoint
543 DynamicList<point>& snappedPoints,
547 const labelList& faceOwn = mesh_.faceOwner();
548 const labelList& faceNei = mesh_.faceNeighbour();
552 bitSet isBoundaryPoint(mesh_.nPoints());
553 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
555 for (
const polyPatch& pp :
patches)
559 for (
const label facei : pp.range())
561 const face&
f = mesh_.faces()[facei];
563 isBoundaryPoint.
set(
f);
569 const point greatPoint(GREAT, GREAT, GREAT);
571 pointField collapsedPoint(mesh_.nPoints(), greatPoint);
575 DynamicList<point, 64> localTriPoints(100);
578 forAll(mesh_.pointFaces(), pointi)
580 constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
582 if (isBoundaryPoint.test(pointi))
591 for (
const label facei :
pFaces)
595 (cellCutType_[faceOwn[facei]] & realCut) != 0
597 mesh_.isInternalFace(facei)
598 && (cellCutType_[faceNei[facei]] & realCut) != 0
615 localPointCells.clear();
616 localTriPoints.clear();
618 for (
const label facei :
pFaces)
620 const label own = faceOwn[facei];
626 if (localPointCells.insert(own))
628 genPointTris(pVals, pointi, facei, own, localTriPoints);
644 if (mesh_.isInternalFace(facei))
646 const label nei = faceNei[facei];
650 if (localPointCells.insert(nei))
652 genPointTris(pVals, pointi, facei, nei, localTriPoints);
670 if (localTriPoints.size() == 3)
674 points.transfer(localTriPoints);
681 else if (localTriPoints.size())
700 label
nZones = surf.markZones
709 scalar minCos = GREAT;
710 const vector& n0 = surf.faceNormals()[0];
711 for (label i = 1; i < surf.size(); ++i)
713 const vector&
n = surf.faceNormals()[i];
714 scalar cosAngle = (n0 &
n);
715 if (cosAngle < minCos)
722 collapsedPoint[pointi] = calcCentre(surf);
732 minMagSqrEqOp<point>(),
736 snappedPoint.setSize(mesh_.nPoints());
739 forAll(collapsedPoint, pointi)
743 if (
magSqr(collapsedPoint[pointi]) < 0.5*
magSqr(greatPoint))
745 snappedPoint[pointi] = snappedPoints.size();
746 snappedPoints.append(collapsedPoint[pointi]);
754 const bool checkDuplicates,
755 const List<point>& triPoints,
760 label nTris = triPoints.size()/3;
762 if ((triPoints.size() % 3) != 0)
765 <<
"Problem: number of points " << triPoints.size()
782 Pout<<
"isoSurfaceCell : merged from " << triPoints.size()
783 <<
" points down to " << newPoints.size() <<
endl;
799 <<
"Merged points contain duplicates"
800 <<
" when merging with distance " << mergeDistance_ <<
endl
801 <<
"merged:" << newPoints.size() <<
" re-merged:"
802 << newNewPoints.size()
808 List<labelledTri> tris;
810 DynamicList<labelledTri> dynTris(nTris);
812 DynamicList<label> newToOldTri(nTris);
814 for (label oldTriI = 0; oldTriI < nTris; ++oldTriI)
818 triPointReverseMap[rawPointi],
819 triPointReverseMap[rawPointi+1],
820 triPointReverseMap[rawPointi+2],
823 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
825 newToOldTri.append(oldTriI);
832 triMap.transfer(newToOldTri);
833 tris.transfer(dynTris);
845 Pout<<
"isoSurfaceCell : merged from " << nTris
846 <<
" down to " << tris.size() <<
" triangles." <<
endl;
852 centres[triI] = tris[triI].centre(newPoints);
868 Pout<<
"isoSurfaceCell : detected "
869 << centres.size()-mergedCentres.size()
870 <<
" duplicate triangles." <<
endl;
877 DynamicList<label> newToOldTri(tris.size());
878 labelList newToMaster(mergedCentres.size(), -1);
881 label mergedI = oldToMerged[triI];
883 if (newToMaster[mergedI] == -1)
885 newToMaster[mergedI] = triI;
886 newToOldTri.append(triMap[triI]);
887 tris[newTriI++] = tris[triI];
891 triMap.transfer(newToOldTri);
892 tris.setSize(newTriI);
900 void Foam::isoSurfaceCell::calcAddressing
902 const triSurface& surf,
903 List<FixedList<label, 3>>& faceEdges,
906 Map<labelList>& edgeFacesRest
915 const labelledTri& tri = surf[triI];
916 edgeCentres[edgeI++] = 0.5*(
points[tri[0]]+
points[tri[1]]);
917 edgeCentres[edgeI++] = 0.5*(
points[tri[1]]+
points[tri[2]]);
918 edgeCentres[edgeI++] = 0.5*(
points[tri[2]]+
points[tri[0]]);
934 Pout<<
"isoSurfaceCell : detected "
935 << mergedCentres.size()
936 <<
" edges on " << surf.size() <<
" triangles." <<
endl;
946 faceEdges.setSize(surf.size());
950 faceEdges[triI][0] = oldToMerged[edgeI++];
951 faceEdges[triI][1] = oldToMerged[edgeI++];
952 faceEdges[triI][2] = oldToMerged[edgeI++];
957 edgeFace0.setSize(mergedCentres.size());
959 edgeFace1.setSize(mergedCentres.size());
961 edgeFacesRest.clear();
963 forAll(oldToMerged, oldEdgeI)
965 label triI = oldEdgeI / 3;
966 label edgeI = oldToMerged[oldEdgeI];
968 if (edgeFace0[edgeI] == -1)
970 edgeFace0[edgeI] = triI;
972 else if (edgeFace1[edgeI] == -1)
974 edgeFace1[edgeI] = triI;
985 if (iter != edgeFacesRest.end())
988 label sz = eFaces.size();
989 eFaces.setSize(sz+1);
994 edgeFacesRest.insert(edgeI,
labelList(1, triI));
1001 bool Foam::isoSurfaceCell::danglingTriangle
1003 const FixedList<label, 3>& fEdges,
1008 for (
const label edgei : fEdges)
1010 if (edgeFace1[edgei] == -1)
1016 return (nOpen == 1 || nOpen == 2 || nOpen == 3);
1020 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1022 const List<FixedList<label, 3>>& faceEdges,
1025 const Map<labelList>& edgeFacesRest,
1029 keepTriangles.setSize(faceEdges.size());
1030 keepTriangles =
true;
1032 label nDangling = 0;
1040 const label edgeI = iter.key();
1041 const labelList& otherEdgeFaces = iter.val();
1044 if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1046 keepTriangles[edgeFace0[edgeI]] =
false;
1049 if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1051 keepTriangles[edgeFace1[edgeI]] =
false;
1054 for (
const label triI : otherEdgeFaces)
1056 if (danglingTriangle(faceEdges[triI], edgeFace1))
1058 keepTriangles[triI] =
false;
1069 const triSurface&
s,
1077 ListOps::createWithValue<bool>(
s.size(), newToOldFaces,
true,
false)
1080 newToOldPoints.setSize(
s.points().size());
1081 oldToNewPoints.setSize(
s.points().size());
1082 oldToNewPoints = -1;
1086 forAll(include, oldFacei)
1088 if (include[oldFacei])
1091 for (
const label oldPointi :
s[oldFacei])
1093 if (oldToNewPoints[oldPointi] == -1)
1095 oldToNewPoints[oldPointi] = pointi;
1096 newToOldPoints[pointi++] = oldPointi;
1101 newToOldPoints.setSize(pointi);
1106 forAll(newToOldPoints, i)
1108 newPoints[i] =
s.points()[newToOldPoints[i]];
1111 List<labelledTri> newTriangles(newToOldFaces.size());
1116 const labelledTri& tri =
s[newToOldFaces[i]];
1118 newTriangles[i][0] = oldToNewPoints[tri[0]];
1119 newTriangles[i][1] = oldToNewPoints[tri[1]];
1120 newTriangles[i][2] = oldToNewPoints[tri[2]];
1121 newTriangles[i].region() = tri.region();
1125 return triSurface(newTriangles,
s.patches(), newPoints,
true);
1138 const bitSet& ignoreCells
1143 cellCutType_(
mesh.
nCells(), cutType::UNVISITED)
1145 const bool regularise = (params.
filter() != filterType::NONE);
1149 Pout<<
"isoSurfaceCell:" <<
nl
1150 <<
" cell min/max : " <<
minMax(cVals_) <<
nl
1151 <<
" point min/max : " <<
minMax(pVals_) <<
nl
1152 <<
" isoValue : " << iso <<
nl
1153 <<
" filter : " <<
Switch(regularise) <<
nl
1156 <<
" mergeDistance : " << mergeDistance_ <<
nl
1157 <<
" ignoreCells : " << ignoreCells.
count()
1158 <<
" / " << cVals_.size() <<
nl
1163 label nBlockedCells = 0;
1166 nBlockedCells += blockCells(cellCutType_, ignoreCells);
1171 (void)mesh_.tetBasePtIs();
1175 bitSet isTet(mesh_.nCells());
1177 for (label celli = 0; celli < mesh_.nCells(); ++celli)
1187 nCutCells_ = calcCellCuts(cellCutType_);
1191 Pout<<
"isoSurfaceCell : candidate cells cut "
1193 <<
" blocked " << nBlockedCells
1194 <<
" total " << mesh_.nCells() <<
endl;
1199 const auto& fvmesh = dynamicCast<const fvMesh>(
mesh);
1205 "isoSurfaceCell.cutType",
1206 fvmesh.time().timeName(),
1218 forAll(cellCutType_, celli)
1220 debugFld[celli] = cellCutType_[celli];
1223 Pout<<
"Writing cut types:"
1224 << debugField.objectPath() <<
endl;
1247 snappedCc.
setSize(mesh_.nCells());
1253 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.size()
1254 <<
" cell centres to intersection." <<
endl;
1258 label nCellSnaps = snappedPoints.size();
1275 snappedPoint.
setSize(mesh_.nPoints());
1281 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
1282 <<
" vertices to intersection." <<
endl;
1298 mesh_.cellCentres(),
1311 Pout<<
"isoSurfaceCell : generated " << triMeshCells.size()
1312 <<
" unmerged triangles." <<
endl;
1322 if (getClipBounds().valid())
1324 isoSurfacePoint::trimToBox
1330 interpolatedPoints_,
1331 interpolatedOldPoints_,
1332 interpolationWeights_
1334 triMeshCells =
labelField(triMeshCells, trimTriMap);
1340 tmpsurf = stitchTriPoints
1350 Pout<<
"isoSurfaceCell : generated " << triMap.size()
1351 <<
" merged triangles." <<
endl;
1354 if (getClipBounds().valid())
1360 labelList newTriPointMergeMap(nOldPoints, -1);
1361 forAll(trimTriPointMap, trimPointI)
1363 label oldPointI = trimTriPointMap[trimPointI];
1366 label pointI = triPointMergeMap_[trimPointI];
1369 newTriPointMergeMap[oldPointI] = pointI;
1373 triPointMergeMap_.
transfer(newTriPointMergeMap);
1376 meshCells_.setSize(triMap.size());
1379 meshCells_[i] = triMeshCells[triMap[i]];
1386 Pout<<
"isoSurfaceCell : checking " << tmpsurf.size()
1387 <<
" triangles for validity." <<
endl;
1417 label nDangling = markDanglingTriangles
1428 Pout<<
"isoSurfaceCell : detected " << nDangling
1429 <<
" dangling triangles." <<
endl;
1442 tmpsurf = subsetMesh
1449 meshCells_ =
labelField(meshCells_, subsetTriMap);
1469 this->Mesh::transfer(updated);