61 const bool aLower = (pointValues[a] < isoValue);
62 const bool bLower = (pointValues[
b] < isoValue);
63 const bool cLower = (pointValues[
c] < isoValue);
65 return !(aLower == bLower && aLower == cLower);
72 Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
78 if (ignoreCells_.
test(celli))
83 const cell& cFaces = mesh_.
cells()[celli];
87 for (
const label facei : cFaces)
98 const face&
f = mesh_.
faces()[facei];
100 for (label fp = 1; fp <
f.size() - 1; ++fp)
112 const bool cellLower = (cVals_[celli] <
iso_);
115 bool edgeCut =
false;
117 for (
const label facei : cFaces)
128 const face&
f = mesh_.
faces()[facei];
131 for (
const label pointi :
f)
133 if ((pVals_[pointi] <
iso_) != cellLower)
147 const label fp0 = (tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei]);
149 label fp =
f.fcIndex(fp0);
150 for (label i = 2; i <
f.size(); ++i)
152 const label nextFp =
f.fcIndex(fp);
179 for (
const label pointi : cPoints)
181 if ((pVals_[pointi] <
iso_) != cellLower)
187 if (nPyrCuts == cPoints.size())
201 Foam::label Foam::isoSurfaceTopo::calcCutTypes
204 List<cellCutType>& cellCutTypes
207 cellCutTypes.setSize(mesh_.nCells());
210 forAll(cellCutTypes, celli)
212 cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
214 if (cellCutTypes[celli] == CUT)
222 Pout<<
"isoSurfaceCell : candidate cut cells "
223 << nCutCells <<
" / " << mesh_.nCells() <<
endl;
229 Foam::scalar Foam::isoSurfaceTopo::minTetQ
232 const label faceBasePtI
238 mesh_.cellCentres()[mesh_.faceOwner()[facei]],
244 if (mesh_.isInternalFace(facei))
252 mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
263 void Foam::isoSurfaceTopo::fixTetBasePtIs()
267 const faceList& faces = mesh_.faces();
268 const labelList& faceOwner = mesh_.faceOwner();
269 const labelList& faceNeighbour = mesh_.faceNeighbour();
273 tetBasePtIs_ = mesh_.tetBasePtIs();
277 bitSet problemCells(
cells.size());
279 forAll(tetBasePtIs_, facei)
281 if (tetBasePtIs_[facei] == -1)
283 problemCells.set(faceOwner[facei]);
285 if (mesh_.isInternalFace(facei))
287 problemCells.set(faceNeighbour[facei]);
295 bitSet problemPoints(mesh_.points().size());
300 Map<label> pointCount;
303 for (
const label celli : problemCells)
307 for (
const label facei :
cells[celli])
309 for (
const label pointi : faces[facei])
311 ++pointCount(pointi);
320 <<
"point:" << iter.key()
321 <<
" at:" << mesh_.points()[iter.key()]
322 <<
" only used by one face" <<
nl
325 else if (iter.val() == 2)
327 problemPoints.set(iter.key());
338 forAll(tetBasePtIs_, facei)
342 problemCells[faceOwner[facei]]
343 || (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
346 const face&
f = faces[facei];
350 const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
354 !problemPoints[
f[
f.rcIndex(fp0)]]
355 && !problemPoints[
f[
f.fcIndex(fp0)]]
363 scalar maxQ = -GREAT;
369 !problemPoints[
f[
f.rcIndex(fp)]]
370 && !problemPoints[
f[
f.fcIndex(fp)]]
373 const scalar q = minTetQ(facei, fp);
385 tetBasePtIs_[facei] = maxFp;
393 tetBasePtIs_[facei] = 0;
403 Pout<<
"isoSurface : adapted starting point of triangulation on "
404 << nAdapted <<
" faces." <<
endl;
411 Foam::label Foam::isoSurfaceTopo::generatePoint
414 const bool edgeIsDiag,
417 DynamicList<edge>& pointToVerts,
418 DynamicList<label>& pointToFace,
419 DynamicList<bool>& pointFromDiag,
420 EdgeMap<label>& vertsToPoint
423 const auto edgeFnd = vertsToPoint.cfind(
vertices);
426 return edgeFnd.val();
430 label pointi = pointToVerts.size();
433 pointToFace.append(facei);
434 pointFromDiag.append(edgeIsDiag);
435 vertsToPoint.insert(
vertices, pointi);
440 void Foam::isoSurfaceTopo::generateTriPoints
443 const FixedList<scalar, 4>&
s,
444 const FixedList<point, 4>&
p,
445 const FixedList<label, 4>& pIndex,
446 const FixedList<bool, 6>& edgeIsDiag,
448 DynamicList<edge>& pointToVerts,
449 DynamicList<label>& pointToFace,
450 DynamicList<bool>& pointFromDiag,
452 EdgeMap<label>& vertsToPoint,
453 DynamicList<label>& verts
490 edge(pIndex[0], pIndex[1]),
491 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
500 edge(pIndex[0], pIndex[2]),
501 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
510 edge(pIndex[0], pIndex[3]),
511 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
515 if (triIndex == 0x0E)
518 label sz = verts.size();
519 Swap(verts[sz-2], verts[sz-1]);
533 edge(pIndex[1], pIndex[0]),
534 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
543 edge(pIndex[1], pIndex[3]),
544 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
553 edge(pIndex[1], pIndex[2]),
554 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
558 if (triIndex == 0x0D)
561 label sz = verts.size();
562 Swap(verts[sz-2], verts[sz-1]);
576 edge(pIndex[0], pIndex[2]),
577 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
586 edge(pIndex[1], pIndex[3]),
587 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
597 edge(pIndex[0], pIndex[3]),
598 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
610 edge(pIndex[1], pIndex[2]),
611 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
616 if (triIndex == 0x0C)
619 label sz = verts.size();
620 Swap(verts[sz-5], verts[sz-4]);
621 Swap(verts[sz-2], verts[sz-1]);
635 edge(pIndex[2], pIndex[0]),
636 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
645 edge(pIndex[2], pIndex[1]),
646 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
655 edge(pIndex[2], pIndex[3]),
656 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
660 if (triIndex == 0x0B)
663 label sz = verts.size();
664 Swap(verts[sz-2], verts[sz-1]);
678 edge(pIndex[0], pIndex[1]),
679 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
688 edge(pIndex[2], pIndex[3]),
689 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
701 edge(pIndex[0], pIndex[3]),
702 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
712 edge(pIndex[1], pIndex[2]),
713 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
718 if (triIndex == 0x0A)
721 label sz = verts.size();
722 Swap(verts[sz-5], verts[sz-4]);
723 Swap(verts[sz-2], verts[sz-1]);
737 edge(pIndex[0], pIndex[1]),
738 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
747 edge(pIndex[2], pIndex[3]),
748 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
759 edge(pIndex[1], pIndex[3]),
760 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
772 edge(pIndex[0], pIndex[2]),
773 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
777 if (triIndex == 0x09)
780 label sz = verts.size();
781 Swap(verts[sz-5], verts[sz-4]);
782 Swap(verts[sz-2], verts[sz-1]);
796 edge(pIndex[3], pIndex[0]),
797 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
806 edge(pIndex[3], pIndex[2]),
807 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
816 edge(pIndex[3], pIndex[1]),
817 pointToVerts, pointToFace, pointFromDiag, vertsToPoint
820 if (triIndex == 0x07)
823 label sz = verts.size();
824 Swap(verts[sz-2], verts[sz-1]);
832 void Foam::isoSurfaceTopo::generateTriPoints
837 DynamicList<edge>& pointToVerts,
838 DynamicList<label>& pointToFace,
839 DynamicList<bool>& pointFromDiag,
841 EdgeMap<label>& vertsToPoint,
842 DynamicList<label>& verts,
843 DynamicList<label>& faceLabels
846 const faceList& faces = mesh_.faces();
847 const labelList& faceOwner = mesh_.faceOwner();
849 const cell& cFaces = mesh_.cells()[celli];
856 label facei = cFaces[0];
857 const face& f0 = faces[facei];
861 const face& f1 = faces[cFaces[1]];
862 label oppositeI = -1;
866 if (!f0.found(oppositeI))
876 FixedList<bool, 6> edgeIsDiag(
false);
878 if (faceOwner[facei] == celli)
891 label startTrii = verts.size();
925 label nTris = (verts.size()-startTrii)/3;
926 for (label i = 0; i < nTris; ++i)
928 faceLabels.append(facei);
933 const pointField& cellCentres = mesh_.cellCentres();
935 for (
const label facei : cFaces)
939 !mesh_.isInternalFace(facei)
940 && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
946 const face&
f = faces[facei];
948 label fp0 = tetBasePtIs_[facei];
950 label startTrii = verts.size();
958 label fp =
f.fcIndex(fp0);
959 for (label i = 2; i <
f.size(); ++i)
961 label nextFp =
f.fcIndex(fp);
963 FixedList<bool, 6> edgeIsDiag(
false);
967 label p2 =
f[nextFp];
968 if (faceOwner[facei] == celli)
971 if (i != 2) edgeIsDiag[1] =
true;
972 if (i !=
f.size()-1) edgeIsDiag[0] =
true;
976 if (i != 2) edgeIsDiag[0] =
true;
977 if (i !=
f.size()-1) edgeIsDiag[1] =
true;
1010 mesh_.nPoints()+celli
1024 label nTris = (verts.size()-startTrii)/3;
1025 for (label i = 0; i < nTris; ++i)
1027 faceLabels.append(facei);
1034 void Foam::isoSurfaceTopo::triangulateOutside
1036 const bool filterDiag,
1042 DynamicList<face>& compactFaces,
1043 DynamicList<label>& compactCellIDs
1058 if (loop.size() > 2)
1060 compactFaces.append(face(0));
1061 face&
f = compactFaces.last();
1067 label pointi =
mp[loop[i]];
1068 if (filterDiag && pointFromDiag[pointi])
1070 label prevPointi =
mp[loop[loop.fcIndex(i)]];
1073 pointFromDiag[prevPointi]
1074 && (pointToFace[pointi] != pointToFace[prevPointi])
1099 label pointi =
mp[loop[i]];
1103 compactCellIDs.
append(cellID);
1111 const bool filterDiag,
1112 const MeshStorage&
s,
1116 DynamicList<label>& pointCompactMap,
1117 DynamicList<label>& compactCellIDs
1122 pointCompactMap.clear();
1123 compactCellIDs.clear();
1125 DynamicList<face> compactFaces(
s.size()/8);
1127 for (label celli = 0; celli < start.size()-1; ++celli)
1131 label nTris = start[celli+1]-start[celli];
1135 const SubList<face> cellFaces(
s, nTris, start[celli]);
1157 DynamicField<point> compactPoints(
points.size());
1158 pointCompactMap.
clear();
1160 for (face&
f : compactFaces)
1164 label pointi =
f[fp];
1165 label compacti = oldToCompact[pointi];
1168 compacti = compactPoints.size();
1169 oldToCompact[pointi] = compacti;
1171 pointCompactMap.append(pointi);
1178 MeshStorage cpSurface
1180 std::move(compactPoints),
1181 std::move(compactFaces),
1199 const bitSet& ignoreCells
1206 ignoreCells_(ignoreCells)
1210 Pout<<
"isoSurfaceTopo::"
1211 <<
" cell min/max : "
1212 <<
min(cVals_) <<
" / "
1213 <<
max(cVals_) <<
nl
1214 <<
" point min/max : "
1215 <<
min(pVals_) <<
" / "
1216 <<
max(pVals_) <<
nl
1217 <<
" isoValue : " << iso <<
nl
1221 <<
" ignoreCells : " << ignoreCells_.count()
1222 <<
" / " << cVals_.size() <<
nl
1234 if (isA<cyclicACMIPolyPatch>(pp))
1236 ignoreBoundaryFaces_.
setSize(mesh_.nBoundaryFaces());
1239 label facei = pp.
start()+i;
1240 ignoreBoundaryFaces_.set(facei-mesh_.nInternalFaces());
1253 const label nCutCells = calcCutTypes(tet, cellCutTypes);
1275 for (label celli = 0; celli < mesh_.nCells(); ++celli)
1277 startTri[celli] = faceLabels.size();
1278 if (cellCutTypes[celli] != NOTCUT)
1283 tet.
isA(mesh_, celli),
1294 for (label i = startTri[celli]; i < faceLabels.size(); ++i)
1296 cellLabels.
append(celli);
1300 startTri[mesh_.nCells()] = faceLabels.size();
1303 pointToVerts_.
transfer(pointToVerts);
1304 meshCells_.transfer(cellLabels);
1305 pointToFace_.transfer(pointToFace);
1311 mesh_.cellCentres(),
1318 faceList allTris(faceLabels.size());
1320 for (
face& allTri : allTris)
1323 allTri[0] = verts[verti++];
1324 allTri[1] = verts[verti++];
1325 allTri[2] = verts[verti++];
1338 MeshStorage::transfer(updated);
1351 Pout<<
"isoSurfaceTopo : generated " << size() <<
" faces." <<
endl;
1355 if (filter != filterType::NONE)
1361 MeshStorage::operator=
1365 (filter == filterType::DIAGCELL ?
true :
false),
1378 meshCells_.transfer(compactCellIDs);
1382 Pout<<
"isoSurfaceTopo :"
1383 <<
" after removing cell centre and face-diag triangles : "
1388 if (filter == filterType::DIAGCELL)
1413 if (!ignoreCells_.empty())
1415 const labelList& faceOwner = mesh_.faceOwner();
1416 const labelList& faceNeighbour = mesh_.faceNeighbour();
1427 if (ignoreCells_.test(faceOwner[facei])) ++
n;
1428 if (ignoreCells_.test(faceNeighbour[facei])) ++
n;
1434 isBoundaryPoint.set(
mesh.
faces()[facei]);
1455 const labelList& eFaces = edgeFaces[edgei];
1456 if (eFaces.size() == 1)
1460 const edge&
e = edges()[edgei];
1461 const edge& verts0 = pointToVerts_[
mp[
e[0]]];
1462 const edge& verts1 = pointToVerts_[
mp[
e[1]]];
1465 isBoundaryPoint[verts0[0]]
1466 && isBoundaryPoint[verts0[1]]
1467 && isBoundaryPoint[verts1[0]]
1468 && isBoundaryPoint[verts1[1]]
1483 Pout<<
"isoSurfaceTopo :"
1485 <<
" faces since on open edges" <<
endl;
1501 MeshStorage::subsetMesh
1508 MeshStorage::transfer(filteredSurf);