71 bool Foam::isoSurface::noTransform(
const tensor& tt)
const
74 (
mag(tt.xx()-1) < mergeDistance_)
75 && (
mag(tt.yy()-1) < mergeDistance_)
76 && (
mag(tt.zz()-1) < mergeDistance_)
77 && (
mag(tt.xy()) < mergeDistance_)
78 && (
mag(tt.xz()) < mergeDistance_)
79 && (
mag(tt.yx()) < mergeDistance_)
80 && (
mag(tt.yz()) < mergeDistance_)
81 && (
mag(tt.zx()) < mergeDistance_)
82 && (
mag(tt.zy()) < mergeDistance_);
86 bool Foam::isoSurface::collocatedPatch(
const polyPatch& pp)
88 const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
89 return cpp.parallel() && !cpp.separated();
95 const coupledPolyPatch& pp
99 bitSet collocated(pp.size());
101 if (isA<processorPolyPatch>(pp))
103 if (collocatedPatch(pp))
111 else if (isA<cyclicPolyPatch>(pp))
113 if (collocatedPatch(pp))
124 <<
"Unhandled coupledPolyPatch type " << pp.type()
131 void Foam::isoSurface::syncUnseparatedPoints
134 const point& nullValue
139 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
148 isA<processorPolyPatch>(
p)
150 && collocatedPatch(
p)
153 const processorPolyPatch& pp =
154 refCast<const processorPolyPatch>(
p);
156 const labelList& meshPts = pp.meshPoints();
157 const labelList& nbrPts = pp.neighbPoints();
163 const label nbrPointi = nbrPts[pointi];
164 patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
181 isA<processorPolyPatch>(
p)
183 && collocatedPatch(
p)
186 const processorPolyPatch& pp =
187 refCast<const processorPolyPatch>(
p);
198 fromNbr >> nbrPatchInfo;
201 const labelList& meshPts = pp.meshPoints();
205 const label meshPointi = meshPts[pointi];
208 pointValues[meshPointi],
219 if (isA<cyclicPolyPatch>(
p))
221 const cyclicPolyPatch& cycPatch =
222 refCast<const cyclicPolyPatch>(
p);
224 if (cycPatch.owner() && collocatedPatch(cycPatch))
228 const edgeList& coupledPoints = cycPatch.coupledPoints();
229 const labelList& meshPts = cycPatch.meshPoints();
230 const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
231 const labelList& nbrMeshPoints = nbrPatch.meshPoints();
238 const edge&
e = coupledPoints[i];
239 half0Values[i] = pointValues[meshPts[
e[0]]];
240 half1Values[i] = pointValues[nbrMeshPoints[
e[1]]];
245 const edge&
e = coupledPoints[i];
246 const label
p0 = meshPts[
e[0]];
247 const label p1 = nbrMeshPoints[
e[1]];
249 minEqOp<point>()(pointValues[
p0], half1Values[i]);
250 minEqOp<point>()(pointValues[p1], half0Values[i]);
257 const globalMeshData& pd = mesh_.globalData();
259 if (pd.nGlobalPoints() > 0)
262 pointField sharedPts(pd.nGlobalPoints(), nullValue);
264 forAll(pd.sharedPointLabels(), i)
266 const label meshPointi = pd.sharedPointLabels()[i];
268 sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointi];
277 forAll(pd.sharedPointLabels(), i)
279 const label meshPointi = pd.sharedPointLabels()[i];
280 pointValues[meshPointi] = sharedPts[pd.sharedPointAddr()[i]];
286 Foam::scalar Foam::isoSurface::isoFraction
292 const scalar d = s1-s0;
303 bool Foam::isoSurface::isEdgeOfFaceCut
313 const bool fpLower = (pVals[
f[fp]] < iso_);
318 || fpLower != neiLower
319 || fpLower != (pVals[
f[
f.fcIndex(fp)]] < iso_)
330 void Foam::isoSurface::getNeighbour
341 const labelList& own = mesh_.faceOwner();
342 const labelList& nei = mesh_.faceNeighbour();
344 if (mesh_.isInternalFace(facei))
346 const label nbr = (own[facei] == celli ? nei[facei] : own[facei]);
347 nbrValue = cVals[nbr];
348 nbrPoint = meshC[nbr];
352 const label bFacei = facei-mesh_.nInternalFaces();
353 const label patchi = boundaryRegion[bFacei];
354 const polyPatch& pp = mesh_.boundaryMesh()[patchi];
355 const label patchFacei = facei-pp.start();
357 nbrValue = cVals.boundaryField()[patchi][patchFacei];
358 nbrPoint = meshC.boundaryField()[patchi][patchFacei];
363 void Foam::isoSurface::calcCutTypes
371 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
372 const labelList& own = mesh_.faceOwner();
373 const labelList& nei = mesh_.faceNeighbour();
375 faceCutType_.
setSize(mesh_.nFaces());
376 faceCutType_ = NOTCUT;
378 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
380 const scalar& ownValue = cVals[own[facei]];
381 const bool ownLower = (ownValue < iso_);
396 const bool neiLower = (nbrValue < iso_);
398 if (ownLower != neiLower)
400 faceCutType_[facei] = CUT;
406 const face
f = mesh_.faces()[facei];
408 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
410 faceCutType_[facei] = CUT;
415 for (
const polyPatch& pp :
patches)
417 label facei = pp.start();
421 const scalar& ownValue = cVals[own[facei]];
422 const bool ownLower = (ownValue < iso_);
437 const bool neiLower = (nbrValue < iso_);
439 if (ownLower != neiLower)
441 faceCutType_[facei] = CUT;
447 const face
f = mesh_.faces()[facei];
449 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
451 faceCutType_[facei] = CUT;
460 cellCutType_.
setSize(mesh_.nCells());
461 cellCutType_ = NOTCUT;
466 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
468 if (faceCutType_[facei] == NOTCUT)
473 if (cellCutType_[own[facei]] == NOTCUT)
475 cellCutType_[own[facei]] = CUT;
478 if (cellCutType_[nei[facei]] == NOTCUT)
480 cellCutType_[nei[facei]] = CUT;
488 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
490 if (faceCutType_[facei] == NOTCUT)
495 if (cellCutType_[own[facei]] == NOTCUT)
497 cellCutType_[own[facei]] = CUT;
504 Pout<<
"isoSurface : candidate cut cells "
505 << nCutCells_ <<
" / " << mesh_.nCells() <<
endl;
510 Foam::point Foam::isoSurface::calcCentre(
const triSurface&
s)
518 sum +=
s[i].centre(
s.points());
524 void Foam::isoSurface::calcSnappedCc
531 DynamicList<point>& snappedPoints,
538 snappedCc.setSize(mesh_.nCells());
542 DynamicList<point, 64> localTriPoints(64);
544 forAll(mesh_.cells(), celli)
546 if (cellCutType_[celli] == CUT)
548 scalar cVal = cVals[celli];
550 const cell& cFaces = mesh_.cells()[celli];
552 localTriPoints.clear();
561 label facei = cFaces[cFacei];
577 FixedList<scalar, 3>
s;
578 FixedList<point, 3> pt;
581 s[2] = isoFraction(cVal, nbrValue);
582 pt[2] = (1.0-
s[2])*cc[celli] +
s[2]*nbrPoint;
584 const face&
f = mesh_.faces()[cFaces[cFacei]];
590 s[0] = isoFraction(cVal, pVals[
p0]);
591 pt[0] = (1.0-
s[0])*cc[celli] +
s[0]*pts[
p0];
594 label p1 =
f[
f.fcIndex(fp)];
595 s[1] = isoFraction(cVal, pVals[p1]);
596 pt[1] = (1.0-
s[1])*cc[celli] +
s[1]*pts[p1];
600 (
s[0] >= 0.0 &&
s[0] <= 0.5)
601 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
602 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
605 localTriPoints.append(pt[0]);
606 localTriPoints.append(pt[1]);
607 localTriPoints.append(pt[2]);
614 if (
s[i] >= 0.0 &&
s[i] <= 0.5)
616 otherPointSum += pt[i];
624 if (localTriPoints.size() == 0)
630 snappedCc[celli] = snappedPoints.size();
631 snappedPoints.append(otherPointSum/nOther);
638 else if (localTriPoints.size() == 3)
642 points.transfer(localTriPoints);
643 snappedCc[celli] = snappedPoints.size();
669 label
nZones = surf.markZones
677 snappedCc[celli] = snappedPoints.size();
678 snappedPoints.append(calcCentre(surf));
689 void Foam::isoSurface::calcSnappedPoint
691 const bitSet& isBoundaryPoint,
697 DynamicList<point>& snappedPoints,
707 DynamicList<point, 64> localTriPoints(100);
709 forAll(mesh_.pointFaces(), pointi)
711 if (isBoundaryPoint.test(pointi))
720 for (
const label facei :
pFaces)
722 if (faceCutType_[facei] == CUT)
735 localTriPoints.clear();
739 for (
const label facei :
pFaces)
743 const face&
f = mesh_.faces()[facei];
744 label own = mesh_.faceOwner()[facei];
761 FixedList<scalar, 4>
s;
762 FixedList<point, 4> pt;
764 label fp =
f.find(pointi);
765 s[0] = isoFraction(pVals[pointi], cVals[own]);
766 pt[0] = (1.0-
s[0])*pts[pointi] +
s[0]*cc[own];
768 s[1] = isoFraction(pVals[pointi], nbrValue);
769 pt[1] = (1.0-
s[1])*pts[pointi] +
s[1]*nbrPoint;
771 label nextPointi =
f[
f.fcIndex(fp)];
772 s[2] = isoFraction(pVals[pointi], pVals[nextPointi]);
773 pt[2] = (1.0-
s[2])*pts[pointi] +
s[2]*pts[nextPointi];
775 label prevPointi =
f[
f.rcIndex(fp)];
776 s[3] = isoFraction(pVals[pointi], pVals[prevPointi]);
777 pt[3] = (1.0-
s[3])*pts[pointi] +
s[3]*pts[prevPointi];
781 (
s[0] >= 0.0 &&
s[0] <= 0.5)
782 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
783 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
786 localTriPoints.append(pt[0]);
787 localTriPoints.append(pt[1]);
788 localTriPoints.append(pt[2]);
792 (
s[0] >= 0.0 &&
s[0] <= 0.5)
793 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
794 && (
s[3] >= 0.0 &&
s[3] <= 0.5)
797 localTriPoints.append(pt[3]);
798 localTriPoints.append(pt[0]);
799 localTriPoints.append(pt[1]);
805 if (
s[i] >= 0.0 &&
s[i] <= 0.5)
807 otherPointSum += pt[i];
813 if (localTriPoints.size() == 0)
819 collapsedPoint[pointi] = otherPointSum/nOther;
822 else if (localTriPoints.size() == 3)
826 points.transfer(localTriPoints);
848 label
nZones = surf.markZones
856 collapsedPoint[pointi] = calcCentre(surf);
863 syncUnseparatedPoints(collapsedPoint,
point::max);
866 snappedPoint.setSize(mesh_.nPoints());
869 forAll(collapsedPoint, pointi)
873 snappedPoint[pointi] = snappedPoints.size();
874 snappedPoints.append(collapsedPoint[pointi]);
882 const bool checkDuplicates,
883 const List<point>& triPoints,
888 label nTris = triPoints.size()/3;
890 if ((triPoints.size() % 3) != 0)
893 <<
"Problem: number of points " << triPoints.size()
910 Pout<<
"isoSurface : merged from " << triPoints.size()
911 <<
" down to " << newPoints.size() <<
" unique points." <<
endl;
927 <<
"Merged points contain duplicates"
928 <<
" when merging with distance " << mergeDistance_ <<
endl
929 <<
"merged:" << newPoints.size() <<
" re-merged:"
930 << newNewPoints.size()
936 List<labelledTri> tris;
938 DynamicList<labelledTri> dynTris(nTris);
940 DynamicList<label> newToOldTri(nTris);
942 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
946 triPointReverseMap[rawPointi],
947 triPointReverseMap[rawPointi+1],
948 triPointReverseMap[rawPointi+2],
953 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
955 newToOldTri.append(oldTriI);
960 triMap.transfer(newToOldTri);
961 tris.transfer(dynTris);
976 DynamicList<label> newToOldTri(tris.size());
980 const labelledTri& tri = tris[triI];
989 label nbrTriI =
pFaces[i];
991 if (nbrTriI > triI && (tris[nbrTriI] == tri))
1004 label newTriI = newToOldTri.size();
1005 newToOldTri.append(triMap[triI]);
1006 tris[newTriI] = tris[triI];
1010 triMap.transfer(newToOldTri);
1011 tris.setSize(triMap.size());
1015 Pout<<
"isoSurface : merged from " << nTris
1016 <<
" down to " << tris.size() <<
" unique triangles." <<
endl;
1026 const labelledTri&
f = surf[facei];
1027 const labelList& fFaces = surf.faceFaces()[facei];
1031 label nbrFacei = fFaces[i];
1033 if (nbrFacei <= facei)
1039 const labelledTri& nbrF = surf[nbrFacei];
1045 <<
" triangle " << facei <<
" vertices " <<
f
1046 <<
" fc:" <<
f.centre(surf.points())
1047 <<
" has the same vertices as triangle " << nbrFacei
1048 <<
" vertices " << nbrF
1049 <<
" fc:" << nbrF.centre(surf.points())
1061 void Foam::isoSurface::trimToPlanes
1063 const PtrList<plane>& planes,
1065 DynamicList<point>& newTriPoints
1069 DynamicList<triPoints> insideTrisA;
1070 storeOp insideOpA(insideTrisA);
1073 DynamicList<triPoints> insideTrisB;
1074 storeOp insideOpB(insideTrisB);
1076 triPointRef::dummyOp dop;
1079 insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
1086 const plane& pln = planes[faceI];
1090 insideOpB.tris_.clear();
1091 forAll(insideOpA.tris_, i)
1093 const triPoints& tri = insideOpA.tris_[i];
1094 triPointRef(tri).sliceWithPlane(pln, insideOpB, dop);
1099 insideOpA.tris_.clear();
1100 forAll(insideOpB.tris_, i)
1102 const triPoints& tri = insideOpB.tris_[i];
1103 triPointRef(tri).sliceWithPlane(pln, insideOpA, dop);
1113 forAll(insideOpA.tris_, i)
1115 const triPoints& tri = insideOpA.tris_[i];
1116 newTriPoints.append(tri[0]);
1117 newTriPoints.append(tri[1]);
1118 newTriPoints.append(tri[2]);
1123 forAll(insideOpB.tris_, i)
1125 const triPoints& tri = insideOpB.tris_[i];
1126 newTriPoints.append(tri[0]);
1127 newTriPoints.append(tri[1]);
1128 newTriPoints.append(tri[2]);
1134 void Foam::isoSurface::trimToBox
1136 const treeBoundBox& bb,
1137 DynamicList<point>& triPoints,
1138 DynamicList<label>& triMap
1143 Pout<<
"isoSurface : trimming to " << bb <<
endl;
1151 planes.set(faceI,
new plane(bb.faceCentre(faceI), -
n));
1154 label nTris = triPoints.size()/3;
1156 DynamicList<point> newTriPoints(triPoints.size()/16);
1157 triMap.setCapacity(nTris/16);
1160 for (label triI = 0; triI < nTris; triI++)
1162 const point&
p0 = triPoints[vertI++];
1163 const point& p1 = triPoints[vertI++];
1164 const point& p2 = triPoints[vertI++];
1166 label oldNPoints = newTriPoints.
size();
1174 label nCells = (newTriPoints.size()-oldNPoints)/3;
1175 for (label i = 0; i < nCells; i++)
1177 triMap.append(triI);
1183 Pout<<
"isoSurface : trimmed from " << nTris
1184 <<
" down to " << triMap.size()
1185 <<
" triangles." <<
endl;
1188 triPoints.transfer(newTriPoints);
1192 void Foam::isoSurface::trimToBox
1194 const treeBoundBox& bb,
1195 DynamicList<point>& triPoints,
1196 DynamicList<label>& triMap,
1199 List<FixedList<label, 3>>& interpolatedOldPoints,
1200 List<FixedList<scalar, 3>>& interpolationWeights
1203 const List<point> oldTriPoints(triPoints);
1206 trimToBox(bb, triPoints, triMap);
1213 label sz = oldTriPoints.size()/100;
1214 DynamicList<label> dynInterpolatedPoints(sz);
1215 DynamicList<FixedList<label, 3>> dynInterpolatedOldPoints(sz);
1216 DynamicList<FixedList<scalar, 3>> dynInterpolationWeights(sz);
1219 triPointMap.setSize(triPoints.size());
1222 label oldTriI = triMap[triI];
1225 for (label i = 0; i < 3; i++)
1227 label pointI = 3*triI+i;
1228 const point& pt = triPoints[pointI];
1231 label matchPointI = -1;
1232 for (label j = 0; j < 3; j++)
1234 label oldPointI = 3*oldTriI+j;
1235 if (pt == oldTriPoints[oldPointI])
1237 matchPointI = oldPointI;
1242 triPointMap[pointI] = matchPointI;
1245 if (matchPointI == -1)
1247 dynInterpolatedPoints.append(pointI);
1249 FixedList<label, 3> oldPoints
1251 {3*oldTriI, 3*oldTriI+1, 3*oldTriI+2}
1253 dynInterpolatedOldPoints.append(oldPoints);
1257 FixedList<scalar, 3> weights({bary.a(), bary.b(), bary.c()});
1259 dynInterpolationWeights.append(weights);
1264 interpolatedPoints.transfer(dynInterpolatedPoints);
1265 interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
1266 interpolationWeights.transfer(dynInterpolationWeights);
1272 const triSurface&
s,
1280 ListOps::createWithValue<bool>(
s.size(), newToOldFaces,
true,
false)
1283 newToOldPoints.setSize(
s.points().size());
1284 oldToNewPoints.setSize(
s.points().size());
1285 oldToNewPoints = -1;
1289 forAll(include, oldFacei)
1291 if (include[oldFacei])
1294 const labelledTri& tri =
s[oldFacei];
1298 label oldPointi = tri[fp];
1300 if (oldToNewPoints[oldPointi] == -1)
1302 oldToNewPoints[oldPointi] = pointi;
1303 newToOldPoints[pointi++] = oldPointi;
1308 newToOldPoints.setSize(pointi);
1313 forAll(newToOldPoints, i)
1315 newPoints[i] =
s.points()[newToOldPoints[i]];
1318 List<labelledTri> newTriangles(newToOldFaces.size());
1323 const labelledTri& tri =
s[newToOldFaces[i]];
1325 newTriangles[i][0] = oldToNewPoints[tri[0]];
1326 newTriangles[i][1] = oldToNewPoints[tri[1]];
1327 newTriangles[i][2] = oldToNewPoints[tri[2]];
1328 newTriangles[i].region() = tri.region();
1332 return triSurface(newTriangles,
s.patches(), newPoints,
true);
1345 const scalar mergeTol
1353 (filter != filterType::NONE),
1365 const bool regularise,
1367 const scalar mergeTol
1371 mesh_(cellValues.mesh()),
1372 pVals_(pointValues),
1373 regularise_(regularise),
1374 mergeDistance_(mergeTol*mesh_.bounds().mag())
1378 Pout<<
"isoSurface:" <<
nl
1379 <<
" isoField : " << cellValues.name() <<
nl
1380 <<
" cell min/max : "
1383 <<
" point min/max : "
1384 <<
min(pVals_) <<
" / "
1385 <<
max(pVals_) <<
nl
1386 <<
" isoValue : " << iso <<
nl
1387 <<
" filter : " <<
Switch(regularise_) <<
nl
1388 <<
" mergeTol : " << mergeTol <<
nl
1400 cValsPtr_.reset(adaptPatchFields(cellValues).ptr());
1413 mesh_.pointsInstance(),
1422 mesh_.cellCentres(),
1434 meshC.boundaryField()[patchi]
1439 collocatedFaces(refCast<const coupledPolyPatch>(pp))
1444 if (!isCollocated[i])
1446 pfld[i] = mesh_.faceCentres()[pp.
start()+i];
1450 else if (isA<emptyPolyPatch>(pp))
1452 auto& bfld =
const_cast<slicedVolVectorField::Boundary&
>
1454 meshC.boundaryField()
1458 bfld.set(patchi,
nullptr);
1466 mesh_.boundary()[patchi],
1472 bfld[patchi] = pp.
patchSlice(mesh_.faceCentres());
1485 label facei = pp.
start();
1501 const fvMesh& fvm = mesh_;
1507 "isoSurface.cutType",
1520 forAll(cellCutType_, celli)
1522 debugFld[celli] = cellCutType_[celli];
1525 Pout<<
"Writing cut types:"
1526 << debugField.objectPath() <<
endl;
1551 snappedCc.
setSize(mesh_.nCells());
1559 Pout<<
"isoSurface : shifted " << snappedPoints.size()
1560 <<
" cell centres to intersection." <<
endl;
1563 label nCellSnaps = snappedPoints.size();
1571 bitSet isBoundaryPoint(mesh_.nPoints());
1581 refCast<const coupledPolyPatch>(pp);
1583 bitSet isCollocated(collocatedFaces(cpp));
1587 if (!isCollocated[i])
1589 const face&
f = mesh_.faces()[cpp.
start()+i];
1591 isBoundaryPoint.
set(
f);
1599 const face&
f = mesh_.faces()[pp.start()+i];
1601 isBoundaryPoint.
set(
f);
1620 snappedPoint.
setSize(mesh_.nPoints());
1626 Pout<<
"isoSurface : shifted " << snappedPoints.size()-nCellSnaps
1627 <<
" vertices to intersection." <<
endl;
1656 Pout<<
"isoSurface : generated " << triMeshCells.size()
1658 <<
" unmerged points." <<
endl;
1667 if (bounds_.valid())
1675 interpolatedPoints_,
1676 interpolatedOldPoints_,
1677 interpolationWeights_
1679 triMeshCells =
labelField(triMeshCells, trimTriMap);
1685 tmpsurf = stitchTriPoints
1695 Pout<<
"isoSurface : generated " << triMap.size()
1696 <<
" merged triangles." <<
endl;
1700 if (bounds_.valid())
1706 labelList newTriPointMergeMap(nOldPoints, -1);
1707 forAll(trimTriPointMap, trimPointI)
1709 label oldPointI = trimTriPointMap[trimPointI];
1712 label pointI = triPointMergeMap_[trimPointI];
1715 newTriPointMergeMap[oldPointI] = pointI;
1719 triPointMergeMap_.
transfer(newTriPointMergeMap);
1722 meshCells_.setSize(triMap.size());
1725 meshCells_[i] = triMeshCells[triMap[i]];
1731 Pout<<
"isoSurface : checking " << tmpsurf.size()
1732 <<
" triangles for validity." <<
endl;
1740 Pout<<
"Dumping surface to " << stlFile <<
endl;
1741 tmpsurf.
write(stlFile);
1759 this->MeshStorage::transfer(updated);