44 Foam::scalar Foam::isoSurfaceCell::isoFraction
81 Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
89 const cell& cFaces = mesh_.cells()[cellI];
91 if (isTet.get(cellI) == 1)
95 const face&
f = mesh_.faces()[cFaces[cFaceI]];
97 for (label fp = 1; fp < f.size() - 1; fp++)
99 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
101 bool aLower = (pointValues[tri[0]] < iso_);
102 bool bLower = (pointValues[tri[1]] < iso_);
103 bool cLower = (pointValues[tri[2]] < iso_);
105 if (aLower == bLower && aLower == cLower)
117 bool cellLower = (cellValues[cellI] < iso_);
120 bool edgeCut =
false;
124 const face& f = mesh_.faces()[cFaces[cFaceI]];
129 if ((pointValues[f[fp]] < iso_) != cellLower)
141 for (label fp = 1; fp < f.size() - 1; fp++)
143 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
149 bool aLower = (pointValues[tri[0]] < iso_);
150 bool bLower = (pointValues[tri[1]] < iso_);
151 bool cLower = (pointValues[tri[2]] < iso_);
153 if (aLower == bLower && aLower == cLower)
174 const labelList& cPoints = mesh_.cellPoints(cellI);
180 if ((pointValues[cPoints[i]] < iso_) != cellLower)
186 if (nPyrCuts == cPoints.size())
203 void Foam::isoSurfaceCell::calcCutTypes
210 cellCutType_.setSize(mesh_.nCells());
212 forAll(mesh_.cells(), cellI)
214 cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
216 if (cellCutType_[cellI] == CUT)
224 Pout<<
"isoSurfaceCell : detected " << nCutCells_
225 <<
" candidate cut cells." <<
endl;
234 const labelledTri& tri0,
235 const labelledTri& tri1
254 label fp0p1 = tri0.fcIndex(fp0);
255 label fp1p1 = tri1.fcIndex(fp1);
256 label fp1m1 = tri1.rcIndex(fp1);
258 if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
260 common[0] = tri0[fp0];
261 common[1] = tri0[fp0p1];
269 Foam::point Foam::isoSurfaceCell::calcCentre(
const triSurface& s)
275 sum += s[i].centre(s.points());
286 DynamicList<labelledTri, 64>& localTris
291 if (localTris.size() == 1)
293 const labelledTri& tri = localTris[0];
294 info.setPoint(tri.centre(localPoints));
297 else if (localTris.size() == 2)
300 const labelledTri& tri0 = localTris[0];
301 const labelledTri& tri1 = localTris[0];
303 labelPair shared = findCommonPoints(tri0, tri1);
311 tri0.centre(localPoints)
312 + tri1.centre(localPoints)
318 else if (localTris.size())
328 localTris.clearStorage();
331 label nZones = surf.markZones
339 info.setPoint(calcCentre(surf));
348 void Foam::isoSurfaceCell::calcSnappedCc
354 DynamicList<point>& snappedPoints,
361 snappedCc.
setSize(mesh_.nCells());
365 DynamicList<point, 64> localPoints(64);
366 DynamicList<labelledTri, 64> localTris(64);
367 Map<label> pointToLocal(64);
369 forAll(mesh_.cells(), cellI)
371 if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
373 scalar cVal = cVals[cellI];
375 const cell& cFaces = mesh_.cells()[cellI];
379 pointToLocal.clear();
386 const face& f = mesh_.faces()[cFaces[cFaceI]];
390 label pointI = f[fp];
392 scalar s = isoFraction(cVal, pVals[pointI]);
394 if (s >= 0.0 && s <= 0.5)
396 if (pointToLocal.insert(pointI, localPoints.size()))
398 localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
404 if (localPoints.size() == 1)
407 snappedCc[cellI] = snappedPoints.size();
408 snappedPoints.append(localPoints[0]);
416 else if (localPoints.size() == 2)
419 snappedCc[cellI] = snappedPoints.size();
420 snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
428 else if (localPoints.size())
433 const face& f = mesh_.faces()[cFaces[cFaceI]];
439 for (label fp = 1; fp < f.size() - 1; fp++)
441 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
448 FixedList<scalar, 3> s(3);
449 s[0] = isoFraction(cVal, pVals[tri[0]]);
450 s[1] = isoFraction(cVal, pVals[tri[1]]);
451 s[2] = isoFraction(cVal, pVals[tri[2]]);
455 (s[0] >= 0.0 && s[0] <= 0.5)
456 && (s[1] >= 0.0 && s[1] <= 0.5)
457 && (s[2] >= 0.0 && s[2] <= 0.5)
464 pointToLocal[tri[0]],
465 pointToLocal[tri[1]],
466 pointToLocal[tri[2]],
480 snappedCc[cellI] = snappedPoints.size();
481 snappedPoints.append(info.hitPoint());
496 void Foam::isoSurfaceCell::genPointTris
503 DynamicList<point, 64>& localTriPoints
509 for (label fp = 1; fp < f.size() - 1; fp++)
511 triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
525 label
b = tri[tri.fcIndex(index)];
526 label c = tri[tri.rcIndex(index)];
529 FixedList<scalar, 3> s(3);
530 s[0] = isoFraction(pointValues[pointI], pointValues[b]);
531 s[1] = isoFraction(pointValues[pointI], pointValues[c]);
532 s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
536 (s[0] >= 0.0 && s[0] <= 0.5)
537 && (s[1] >= 0.0 && s[1] <= 0.5)
538 && (s[2] >= 0.0 && s[2] <= 0.5)
541 localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
542 localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
543 localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
550 void Foam::isoSurfaceCell::genPointTris
555 DynamicList<point, 64>& localTriPoints
559 const cell& cFaces = mesh_.cells()[cellI];
561 FixedList<label, 4> tet;
563 label face0 = cFaces[0];
564 const face& f0 = mesh_.faces()[face0];
566 if (mesh_.faceOwner()[face0] != cellI)
580 const face& f1 = mesh_.faces()[cFaces[1]];
586 if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
597 label i1 = tet.fcIndex(i0);
598 label i2 = tet.fcIndex(i1);
599 label i3 = tet.fcIndex(i2);
602 FixedList<scalar, 3> s(3);
603 s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
604 s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
605 s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
609 (s[0] >= 0.0 && s[0] <= 0.5)
610 && (s[1] >= 0.0 && s[1] <= 0.5)
611 && (s[2] >= 0.0 && s[2] <= 0.5)
614 localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
615 localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
616 localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
621 void Foam::isoSurfaceCell::calcSnappedPoint
628 DynamicList<point>& snappedPoints,
632 const point greatPoint(VGREAT, VGREAT, VGREAT);
633 pointField collapsedPoint(mesh_.nPoints(), greatPoint);
637 DynamicList<point, 64> localTriPoints(100);
640 forAll(mesh_.pointFaces(), pointI)
642 if (isBoundaryPoint.get(pointI) == 1)
653 label faceI = pFaces[i];
657 cellCutType_[mesh_.faceOwner()[faceI]] == CUT
659 mesh_.isInternalFace(faceI)
660 && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
675 localPointCells.clear();
676 localTriPoints.clear();
680 label faceI = pFaces[pFaceI];
681 const face& f = mesh_.faces()[faceI];
682 label own = mesh_.faceOwner()[faceI];
685 if (isTet.get(own) == 1)
687 if (localPointCells.insert(own))
689 genPointTris(pVals, pointI, own, localTriPoints);
694 genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
697 if (mesh_.isInternalFace(faceI))
699 label nei = mesh_.faceNeighbour()[faceI];
701 if (isTet.get(nei) == 1)
703 if (localPointCells.insert(nei))
705 genPointTris(pVals, pointI, nei, localTriPoints);
710 genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
715 if (localTriPoints.size() == 3)
720 collapsedPoint[pointI] =
sum(points)/points.size();
726 else if (localTriPoints.size())
745 label nZones = surf.markZones
753 collapsedPoint[pointI] = calcCentre(surf);
770 snappedPoint.setSize(mesh_.nPoints());
773 forAll(collapsedPoint, pointI)
775 if (collapsedPoint[pointI] != greatPoint)
777 snappedPoint[pointI] = snappedPoints.size();
778 snappedPoints.append(collapsedPoint[pointI]);
788 const bool checkDuplicates,
789 const List<point>& triPoints,
794 label nTris = triPoints.
size()/3;
796 if ((triPoints.size() % 3) != 0)
799 <<
"Problem: number of points " << triPoints.size()
830 <<
"Merged points contain duplicates"
831 <<
" when merging with distance " << mergeDistance_ <<
endl
832 <<
"merged:" << newPoints.size() <<
" re-merged:"
833 << newNewPoints.size()
839 List<labelledTri> tris;
841 DynamicList<labelledTri> dynTris(nTris);
843 DynamicList<label> newToOldTri(nTris);
845 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
849 triPointReverseMap[rawPointI],
850 triPointReverseMap[rawPointI+1],
851 triPointReverseMap[rawPointI+2],
856 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
858 newToOldTri.append(oldTriI);
863 triMap.transfer(newToOldTri);
864 tris.transfer(dynTris);
876 Pout<<
"isoSurfaceCell : merged from " << nTris
877 <<
" down to " << tris.size() <<
" triangles." <<
endl;
883 centres[triI] = tris[triI].centre(newPoints);
899 Pout<<
"isoSurfaceCell : detected "
900 << centres.size()-mergedCentres.size()
901 <<
" duplicate triangles." <<
endl;
908 DynamicList<label> newToOldTri(tris.size());
909 labelList newToMaster(mergedCentres.size(), -1);
912 label mergedI = oldToMerged[triI];
914 if (newToMaster[mergedI] == -1)
916 newToMaster[mergedI] = triI;
917 newToOldTri.append(triMap[triI]);
918 tris[newTriI++] = tris[triI];
922 triMap.transfer(newToOldTri);
923 tris.setSize(newTriI);
932 bool Foam::isoSurfaceCell::validTri(
const triSurface& surf,
const label faceI)
936 const labelledTri& f = surf[faceI];
940 (f[0] < 0) || (f[0] >= surf.points().size())
941 || (f[1] < 0) || (f[1] >= surf.points().size())
942 || (f[2] < 0) || (f[2] >= surf.points().size())
945 WarningIn(
"validTri(const triSurface&, const label)")
946 <<
"triangle " << faceI <<
" vertices " << f
947 <<
" uses point indices outside point range 0.."
948 << surf.points().size()-1 <<
endl;
953 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
955 WarningIn(
"validTri(const triSurface&, const label)")
956 <<
"triangle " << faceI
957 <<
" uses non-unique vertices " << f
964 const labelList& fFaces = surf.faceFaces()[faceI];
970 label nbrFaceI = fFaces[i];
972 if (nbrFaceI <= faceI)
978 const labelledTri& nbrF = surf[nbrFaceI];
982 ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
983 && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
984 && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
987 WarningIn(
"validTri(const triSurface&, const label)")
988 <<
"triangle " << faceI <<
" vertices " << f
989 <<
" has the same vertices as triangle " << nbrFaceI
990 <<
" vertices " << nbrF
1000 void Foam::isoSurfaceCell::calcAddressing
1002 const triSurface& surf,
1003 List<FixedList<label, 3> >& faceEdges,
1006 Map<labelList>& edgeFacesRest
1015 const labelledTri& tri = surf[triI];
1016 edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1017 edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1018 edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1034 Pout<<
"isoSurfaceCell : detected "
1035 << mergedCentres.size()
1036 <<
" edges on " << surf.size() <<
" triangles." <<
endl;
1041 if (surf.size() == 1)
1043 faceEdges.setSize(1);
1044 faceEdges[0][0] = 0;
1045 faceEdges[0][1] = 1;
1046 faceEdges[0][2] = 2;
1047 edgeFace0.setSize(1);
1049 edgeFace1.setSize(1);
1051 edgeFacesRest.clear();
1058 faceEdges.setSize(surf.size());
1062 faceEdges[triI][0] = oldToMerged[edgeI++];
1063 faceEdges[triI][1] = oldToMerged[edgeI++];
1064 faceEdges[triI][2] = oldToMerged[edgeI++];
1069 edgeFace0.setSize(mergedCentres.size());
1071 edgeFace1.setSize(mergedCentres.size());
1073 edgeFacesRest.clear();
1075 forAll(oldToMerged, oldEdgeI)
1077 label triI = oldEdgeI / 3;
1078 label edgeI = oldToMerged[oldEdgeI];
1080 if (edgeFace0[edgeI] == -1)
1082 edgeFace0[edgeI] = triI;
1084 else if (edgeFace1[edgeI] == -1)
1086 edgeFace1[edgeI] = triI;
1097 if (iter != edgeFacesRest.end())
1100 label sz = eFaces.size();
1101 eFaces.setSize(sz+1);
1106 edgeFacesRest.insert(edgeI,
labelList(1, triI));
1113 void Foam::isoSurfaceCell::walkOrientation
1115 const triSurface& surf,
1116 const List<FixedList<label, 3> >& faceEdges,
1119 const label seedTriI,
1124 DynamicList<label> changedFaces(surf.size());
1126 changedFaces.append(seedTriI);
1128 while (changedFaces.size())
1130 DynamicList<label> newChangedFaces(changedFaces.size());
1134 label triI = changedFaces[i];
1135 const labelledTri& tri = surf[triI];
1136 const FixedList<label, 3>& fEdges = faceEdges[triI];
1140 label edgeI = fEdges[fp];
1144 label p1 = tri[tri.fcIndex(fp)];
1148 edgeFace0[edgeI] != triI
1153 if (nbrI != -1 && flipState[nbrI] == -1)
1155 const labelledTri& nbrTri = surf[nbrI];
1159 label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1161 bool sameOrientation = (p1 == nbrP1);
1163 if (flipState[triI] == 0)
1165 flipState[nbrI] = (sameOrientation ? 0 : 1);
1169 flipState[nbrI] = (sameOrientation ? 1 : 0);
1171 newChangedFaces.append(nbrI);
1176 changedFaces.transfer(newChangedFaces);
1181 void Foam::isoSurfaceCell::orientSurface
1184 const List<FixedList<label, 3> >& faceEdges,
1187 const Map<labelList>& edgeFacesRest
1203 seedTriI < surf.size() && flipState[seedTriI] != -1;
1208 if (seedTriI == surf.size())
1215 flipState[seedTriI] = 0;
1232 if (flipState[triI] == 1)
1234 labelledTri tri(surf[triI]);
1236 surf[triI][0] = tri[0];
1237 surf[triI][1] = tri[2];
1238 surf[triI][2] = tri[1];
1240 else if (flipState[triI] == -1)
1244 "isoSurfaceCell::orientSurface(triSurface&, const label)"
1252 bool Foam::isoSurfaceCell::danglingTriangle
1254 const FixedList<label, 3>& fEdges,
1261 if (edgeFace1[fEdges[i]] == -1)
1267 if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1279 Foam::label Foam::isoSurfaceCell::markDanglingTriangles
1281 const List<FixedList<label, 3> >& faceEdges,
1284 const Map<labelList>& edgeFacesRest,
1288 keepTriangles.setSize(faceEdges.size());
1289 keepTriangles =
true;
1291 label nDangling = 0;
1299 label edgeI = iter.key();
1300 const labelList& otherEdgeFaces = iter();
1303 if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1305 keepTriangles[edgeFace0[edgeI]] =
false;
1308 if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1310 keepTriangles[edgeFace1[edgeI]] =
false;
1313 forAll(otherEdgeFaces, i)
1315 label triI = otherEdgeFaces[i];
1316 if (danglingTriangle(faceEdges[triI], edgeFace1))
1318 keepTriangles[triI] =
false;
1329 const triSurface& s,
1337 createWithValues<boolList>
1346 newToOldPoints.setSize(s.points().size());
1347 oldToNewPoints.setSize(s.points().size());
1348 oldToNewPoints = -1;
1352 forAll(include, oldFacei)
1354 if (include[oldFacei])
1357 const labelledTri& tri = s[oldFacei];
1361 label oldPointI = tri[fp];
1363 if (oldToNewPoints[oldPointI] == -1)
1365 oldToNewPoints[oldPointI] = pointI;
1366 newToOldPoints[pointI++] = oldPointI;
1371 newToOldPoints.setSize(pointI);
1376 forAll(newToOldPoints, i)
1378 newPoints[i] = s.points()[newToOldPoints[i]];
1381 List<labelledTri> newTriangles(newToOldFaces.size());
1386 const labelledTri& tri = s[newToOldFaces[i]];
1388 newTriangles[i][0] = oldToNewPoints[tri[0]];
1389 newTriangles[i][1] = oldToNewPoints[tri[1]];
1390 newTriangles[i][2] = oldToNewPoints[tri[2]];
1391 newTriangles[i].region() = tri.region();
1395 return triSurface(newTriangles, s.patches(), newPoints,
true);
1407 const bool regularise,
1408 const scalar mergeTol
1413 mergeDistance_(mergeTol*mesh.
bounds().
mag())
1422 if (tet.
isA(mesh_, cellI))
1424 isTet.set(cellI, 1);
1439 label faceI = pp.
start();
1442 const face& f = mesh_.faces()[faceI++];
1446 isBoundaryPoint.set(f[fp], 1);
1454 calcCutTypes(isTet, cVals, pVals);
1473 snappedCc.
setSize(mesh_.nCells());
1479 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.
size()
1480 <<
" cell centres to intersection." <<
endl;
1484 label nCellSnaps = snappedPoints.
size();
1502 snappedPoint.
setSize(mesh_.nPoints());
1508 Pout<<
"isoSurfaceCell : shifted " << snappedPoints.
size()-nCellSnaps
1509 <<
" vertices to intersection." <<
endl;
1522 mesh_.cellCentres(),
1535 Pout<<
"isoSurfaceCell : generated " << triMeshCells.
size()
1536 <<
" unmerged triangles." <<
endl;
1554 Pout<<
"isoSurfaceCell : generated " << triMap.
size()
1555 <<
" merged triangles." <<
endl;
1558 meshCells_.setSize(triMap.
size());
1561 meshCells_[i] = triMeshCells[triMap[i]];
1566 Pout<<
"isoSurfaceCell : checking " << size()
1567 <<
" triangles for validity." <<
endl;
1572 validTri(*
this, triI);
1587 calcAddressing(*
this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
1591 label nDangling = markDanglingTriangles
1602 Pout<<
"isoSurfaceCell : detected " << nDangling
1603 <<
" dangling triangles." <<
endl;
1626 meshCells_ =
labelField(meshCells_, subsetTriMap);
1630 orientSurface(*
this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);