73 list.append(std::move(tri));
81 const auto* cpp = isA<coupledPolyPatch>(pp);
82 return bool(cpp) && cpp->parallel() && !cpp->separated();
88 static bool isCollocatedPatch(
const coupledPolyPatch& pp)
90 if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
96 <<
"Unhandled coupledPolyPatch type " << pp.type() <<
nl
108 bool Foam::isoSurfacePoint::noTransform(
const tensor& tt)
const
111 (
mag(tt.xx()-1) < mergeDistance_)
112 && (
mag(tt.yy()-1) < mergeDistance_)
113 && (
mag(tt.zz()-1) < mergeDistance_)
114 && (
mag(tt.xy()) < mergeDistance_)
115 && (
mag(tt.xz()) < mergeDistance_)
116 && (
mag(tt.yx()) < mergeDistance_)
117 && (
mag(tt.yz()) < mergeDistance_)
118 && (
mag(tt.zx()) < mergeDistance_)
119 && (
mag(tt.zy()) < mergeDistance_);
125 const coupledPolyPatch& pp
129 bitSet collocated(pp.size());
131 if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
142 <<
"Unhandled coupledPolyPatch type " << pp.type()
149 void Foam::isoSurfacePoint::syncUnseparatedPoints
152 const point& nullValue
157 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
166 isA<processorPolyPatch>(
p)
171 const processorPolyPatch& pp =
172 refCast<const processorPolyPatch>(
p);
174 const labelList& meshPts = pp.meshPoints();
175 const labelList& nbrPts = pp.neighbPoints();
181 const label nbrPointi = nbrPts[pointi];
182 patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
199 isA<processorPolyPatch>(
p)
204 const processorPolyPatch& pp =
205 refCast<const processorPolyPatch>(
p);
216 fromNbr >> nbrPatchInfo;
219 const labelList& meshPts = pp.meshPoints();
223 const label meshPointi = meshPts[pointi];
226 pointValues[meshPointi],
237 if (isA<cyclicPolyPatch>(
p))
239 const cyclicPolyPatch& cycPatch =
240 refCast<const cyclicPolyPatch>(
p);
246 const edgeList& coupledPoints = cycPatch.coupledPoints();
247 const labelList& meshPts = cycPatch.meshPoints();
248 const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
249 const labelList& nbrMeshPoints = nbrPatch.meshPoints();
256 const edge&
e = coupledPoints[i];
257 half0Values[i] = pointValues[meshPts[
e[0]]];
258 half1Values[i] = pointValues[nbrMeshPoints[
e[1]]];
263 const edge&
e = coupledPoints[i];
264 const label
p0 = meshPts[
e[0]];
265 const label p1 = nbrMeshPoints[
e[1]];
267 minEqOp<point>()(pointValues[
p0], half1Values[i]);
268 minEqOp<point>()(pointValues[p1], half0Values[i]);
275 const globalMeshData& pd = mesh_.globalData();
277 if (pd.nGlobalPoints() > 0)
280 pointField sharedPts(pd.nGlobalPoints(), nullValue);
282 forAll(pd.sharedPointLabels(), i)
284 const label meshPointi = pd.sharedPointLabels()[i];
286 sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointi];
295 forAll(pd.sharedPointLabels(), i)
297 const label meshPointi = pd.sharedPointLabels()[i];
298 pointValues[meshPointi] = sharedPts[pd.sharedPointAddr()[i]];
304 Foam::scalar Foam::isoSurfacePoint::isoFraction
310 const scalar d = s1-s0;
321 bool Foam::isoSurfacePoint::isEdgeOfFaceCut
331 const bool fpLower = (pVals[
f[fp]] < iso_);
336 || fpLower != neiLower
337 || fpLower != (pVals[
f[
f.fcIndex(fp)]] < iso_)
348 void Foam::isoSurfacePoint::getNeighbour
359 const labelList& own = mesh_.faceOwner();
360 const labelList& nei = mesh_.faceNeighbour();
362 if (mesh_.isInternalFace(facei))
364 const label nbr = (own[facei] == celli ? nei[facei] : own[facei]);
365 nbrValue = cVals[nbr];
366 nbrPoint = meshC[nbr];
370 const label bFacei = facei-mesh_.nInternalFaces();
371 const label patchi = boundaryRegion[bFacei];
372 const polyPatch& pp = mesh_.boundaryMesh()[patchi];
373 const label patchFacei = facei-pp.start();
375 nbrValue = cVals.boundaryField()[patchi][patchFacei];
376 nbrPoint = meshC.boundaryField()[patchi][patchFacei];
381 void Foam::isoSurfacePoint::calcCutTypes
389 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
390 const labelList& own = mesh_.faceOwner();
391 const labelList& nei = mesh_.faceNeighbour();
393 faceCutType_.
resize(mesh_.nFaces());
394 faceCutType_ = cutType::NOTCUT;
396 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
398 const scalar ownValue = cVals[own[facei]];
399 const bool ownLower = (ownValue < iso_);
414 const bool neiLower = (nbrValue < iso_);
416 if (ownLower != neiLower)
418 faceCutType_[facei] = cutType::CUT;
424 const face&
f = mesh_.faces()[facei];
426 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
428 faceCutType_[facei] = cutType::CUT;
433 for (
const polyPatch& pp :
patches)
435 for (
const label facei : pp.range())
437 const scalar ownValue = cVals[own[facei]];
438 const bool ownLower = (ownValue < iso_);
453 const bool neiLower = (nbrValue < iso_);
455 if (ownLower != neiLower)
457 faceCutType_[facei] = cutType::CUT;
463 const face&
f = mesh_.faces()[facei];
465 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
467 faceCutType_[facei] = cutType::CUT;
474 cellCutType_.
resize(mesh_.nCells());
475 cellCutType_ = cutType::NOTCUT;
480 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
482 if (faceCutType_[facei] == cutType::NOTCUT)
487 if (cellCutType_[own[facei]] == cutType::NOTCUT)
489 cellCutType_[own[facei]] = cutType::CUT;
492 if (cellCutType_[nei[facei]] == cutType::NOTCUT)
494 cellCutType_[nei[facei]] = cutType::CUT;
502 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
504 if (faceCutType_[facei] == cutType::NOTCUT)
509 if (cellCutType_[own[facei]] == cutType::NOTCUT)
511 cellCutType_[own[facei]] = cutType::CUT;
518 Pout<<
"isoSurfacePoint : candidate cut cells "
519 << nCutCells_ <<
" / " << mesh_.nCells() <<
endl;
524 Foam::point Foam::isoSurfacePoint::calcCentre(
const triSurface&
s)
532 sum +=
s[i].centre(
s.points());
538 void Foam::isoSurfacePoint::calcSnappedCc
545 DynamicList<point>& snappedPoints,
552 snappedCc.setSize(mesh_.nCells());
556 DynamicList<point, 64> localTriPoints(64);
558 forAll(mesh_.cells(), celli)
560 if (cellCutType_[celli] == cutType::CUT)
562 const scalar cVal = cVals[celli];
564 localTriPoints.clear();
571 for (
const label facei : mesh_.cells()[celli])
573 const face&
f = mesh_.faces()[facei];
589 FixedList<scalar, 3>
s;
590 FixedList<point, 3> pt;
593 s[2] = isoFraction(cVal, nbrValue);
594 pt[2] = (1.0-
s[2])*cc[celli] +
s[2]*nbrPoint;
600 s[0] = isoFraction(cVal, pVals[
p0]);
601 pt[0] = (1.0-
s[0])*cc[celli] +
s[0]*pts[
p0];
604 label p1 =
f[
f.fcIndex(fp)];
605 s[1] = isoFraction(cVal, pVals[p1]);
606 pt[1] = (1.0-
s[1])*cc[celli] +
s[1]*pts[p1];
610 (
s[0] >= 0.0 &&
s[0] <= 0.5)
611 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
612 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
615 localTriPoints.append(pt[0]);
616 localTriPoints.append(pt[1]);
617 localTriPoints.append(pt[2]);
624 if (
s[i] >= 0.0 &&
s[i] <= 0.5)
626 otherPointSum += pt[i];
634 if (localTriPoints.size() == 0)
640 snappedCc[celli] = snappedPoints.size();
641 snappedPoints.append(otherPointSum/nOther);
648 else if (localTriPoints.size() == 3)
651 snappedCc[celli] = snappedPoints.size();
652 snappedPoints.append(
sum(localTriPoints)/3);
653 localTriPoints.clear();
678 label
nZones = surf.markZones
686 snappedCc[celli] = snappedPoints.size();
687 snappedPoints.append(calcCentre(surf));
698 void Foam::isoSurfacePoint::calcSnappedPoint
700 const bitSet& isBoundaryPoint,
706 DynamicList<point>& snappedPoints,
716 DynamicList<point, 64> localTriPoints(100);
718 forAll(mesh_.pointFaces(), pointi)
720 if (isBoundaryPoint.test(pointi))
729 for (
const label facei :
pFaces)
731 if (faceCutType_[facei] == cutType::CUT)
744 localTriPoints.clear();
748 for (
const label facei :
pFaces)
752 const face&
f = mesh_.faces()[facei];
753 const label own = mesh_.faceOwner()[facei];
770 FixedList<scalar, 4>
s;
771 FixedList<point, 4> pt;
773 label fp =
f.find(pointi);
774 s[0] = isoFraction(pVals[pointi], cVals[own]);
775 pt[0] = (1.0-
s[0])*pts[pointi] +
s[0]*cc[own];
777 s[1] = isoFraction(pVals[pointi], nbrValue);
778 pt[1] = (1.0-
s[1])*pts[pointi] +
s[1]*nbrPoint;
780 label nextPointi =
f[
f.fcIndex(fp)];
781 s[2] = isoFraction(pVals[pointi], pVals[nextPointi]);
782 pt[2] = (1.0-
s[2])*pts[pointi] +
s[2]*pts[nextPointi];
784 label prevPointi =
f[
f.rcIndex(fp)];
785 s[3] = isoFraction(pVals[pointi], pVals[prevPointi]);
786 pt[3] = (1.0-
s[3])*pts[pointi] +
s[3]*pts[prevPointi];
790 (
s[0] >= 0.0 &&
s[0] <= 0.5)
791 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
792 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
795 localTriPoints.append(pt[0]);
796 localTriPoints.append(pt[1]);
797 localTriPoints.append(pt[2]);
801 (
s[0] >= 0.0 &&
s[0] <= 0.5)
802 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
803 && (
s[3] >= 0.0 &&
s[3] <= 0.5)
806 localTriPoints.append(pt[3]);
807 localTriPoints.append(pt[0]);
808 localTriPoints.append(pt[1]);
814 if (
s[i] >= 0.0 &&
s[i] <= 0.5)
816 otherPointSum += pt[i];
822 if (localTriPoints.size() == 0)
828 collapsedPoint[pointi] = otherPointSum/nOther;
831 else if (localTriPoints.size() == 3)
835 points.transfer(localTriPoints);
857 label
nZones = surf.markZones
865 collapsedPoint[pointi] = calcCentre(surf);
872 syncUnseparatedPoints(collapsedPoint,
point::max);
875 snappedPoint.setSize(mesh_.nPoints());
878 forAll(collapsedPoint, pointi)
882 snappedPoint[pointi] = snappedPoints.size();
883 snappedPoints.append(collapsedPoint[pointi]);
891 const bool checkDuplicates,
892 const List<point>& triPoints,
897 label nTris = triPoints.size()/3;
899 if ((triPoints.size() % 3) != 0)
902 <<
"Problem: number of points " << triPoints.size()
919 Pout<<
"isoSurfacePoint : merged from " << triPoints.size()
920 <<
" down to " << newPoints.size() <<
" unique points." <<
endl;
936 <<
"Merged points contain duplicates"
937 <<
" when merging with distance " << mergeDistance_ <<
endl
938 <<
"merged:" << newPoints.size() <<
" re-merged:"
939 << newNewPoints.size()
945 List<labelledTri> tris;
947 DynamicList<labelledTri> dynTris(nTris);
949 DynamicList<label> newToOldTri(nTris);
951 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
955 triPointReverseMap[rawPointi],
956 triPointReverseMap[rawPointi+1],
957 triPointReverseMap[rawPointi+2],
962 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
964 newToOldTri.append(oldTriI);
969 triMap.transfer(newToOldTri);
970 tris.transfer(dynTris);
985 DynamicList<label> newToOldTri(tris.size());
989 const labelledTri& tri = tris[triI];
998 label nbrTriI =
pFaces[i];
1000 if (nbrTriI > triI && (tris[nbrTriI] == tri))
1013 label newTriI = newToOldTri.size();
1014 newToOldTri.append(triMap[triI]);
1015 tris[newTriI] = tris[triI];
1019 triMap.transfer(newToOldTri);
1020 tris.setSize(triMap.size());
1024 Pout<<
"isoSurfacePoint : merged from " << nTris
1025 <<
" down to " << tris.size() <<
" unique triangles." <<
endl;
1035 const labelledTri&
f = surf[facei];
1036 const labelList& fFaces = surf.faceFaces()[facei];
1040 label nbrFacei = fFaces[i];
1042 if (nbrFacei <= facei)
1048 const labelledTri& nbrF = surf[nbrFacei];
1054 <<
" triangle " << facei <<
" vertices " <<
f
1055 <<
" fc:" <<
f.centre(surf.points())
1056 <<
" has the same vertices as triangle " << nbrFacei
1057 <<
" vertices " << nbrF
1058 <<
" fc:" << nbrF.centre(surf.points())
1070 void Foam::isoSurfacePoint::trimToPlanes
1072 const PtrList<plane>& planes,
1074 DynamicList<point>& newTriPoints
1078 DynamicList<triPoints> insideTrisA;
1079 storeOp insideOpA(insideTrisA);
1082 DynamicList<triPoints> insideTrisB;
1083 storeOp insideOpB(insideTrisB);
1085 triPointRef::dummyOp dop;
1088 insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
1095 const plane& pln = planes[faceI];
1099 insideTrisB.clear();
1100 for (
const triPoints& tri : insideTrisA)
1102 triPointRef(tri).sliceWithPlane(pln, insideOpB, dop);
1107 insideTrisA.clear();
1108 for (
const triPoints& tri : insideTrisB)
1110 triPointRef(tri).sliceWithPlane(pln, insideOpA, dop);
1117 DynamicList<triPoints>& newTris = (useA ? insideTrisA : insideTrisB);
1119 newTriPoints.reserve(newTriPoints.size() + 3*newTris.size());
1122 for (
const triPoints& tri : newTris)
1124 newTriPoints.append(tri[0]);
1125 newTriPoints.append(tri[1]);
1126 newTriPoints.append(tri[2]);
1131 void Foam::isoSurfacePoint::trimToBox
1133 const treeBoundBox& bb,
1134 DynamicList<point>& triPoints,
1135 DynamicList<label>& triMap
1140 Pout<<
"isoSurfacePoint : trimming to " << bb <<
endl;
1148 planes.set(faceI,
new plane(bb.faceCentre(faceI), -
n));
1151 label nTris = triPoints.size()/3;
1153 DynamicList<point> newTriPoints(triPoints.size()/16);
1154 triMap.setCapacity(nTris/16);
1157 for (label triI = 0; triI < nTris; triI++)
1159 const point&
p0 = triPoints[vertI++];
1160 const point& p1 = triPoints[vertI++];
1161 const point& p2 = triPoints[vertI++];
1163 label oldNPoints = newTriPoints.
size();
1171 label nCells = (newTriPoints.size()-oldNPoints)/3;
1172 for (label i = 0; i < nCells; i++)
1174 triMap.append(triI);
1180 Pout<<
"isoSurfacePoint : trimmed from " << nTris
1181 <<
" down to " << triMap.size()
1182 <<
" triangles." <<
endl;
1185 triPoints.transfer(newTriPoints);
1189 void Foam::isoSurfacePoint::trimToBox
1191 const treeBoundBox& bb,
1192 DynamicList<point>& triPoints,
1193 DynamicList<label>& triMap,
1196 List<FixedList<label, 3>>& interpolatedOldPoints,
1197 List<FixedList<scalar, 3>>& interpolationWeights
1200 const List<point> oldTriPoints(triPoints);
1203 trimToBox(bb, triPoints, triMap);
1210 label sz = oldTriPoints.size()/100;
1211 DynamicList<label> dynInterpolatedPoints(sz);
1212 DynamicList<FixedList<label, 3>> dynInterpolatedOldPoints(sz);
1213 DynamicList<FixedList<scalar, 3>> dynInterpolationWeights(sz);
1216 triPointMap.setSize(triPoints.size());
1219 label oldTriI = triMap[triI];
1222 for (label i = 0; i < 3; i++)
1224 label pointI = 3*triI+i;
1225 const point& pt = triPoints[pointI];
1228 label matchPointI = -1;
1229 for (label j = 0; j < 3; j++)
1231 label oldPointI = 3*oldTriI+j;
1232 if (pt == oldTriPoints[oldPointI])
1234 matchPointI = oldPointI;
1239 triPointMap[pointI] = matchPointI;
1242 if (matchPointI == -1)
1244 dynInterpolatedPoints.append(pointI);
1246 FixedList<label, 3> oldPoints
1248 {3*oldTriI, 3*oldTriI+1, 3*oldTriI+2}
1250 dynInterpolatedOldPoints.append(oldPoints);
1254 FixedList<scalar, 3> weights({bary.a(), bary.b(), bary.c()});
1256 dynInterpolationWeights.append(weights);
1261 interpolatedPoints.transfer(dynInterpolatedPoints);
1262 interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
1263 interpolationWeights.transfer(dynInterpolationWeights);
1269 const triSurface&
s,
1277 ListOps::createWithValue<bool>(
s.size(), newToOldFaces,
true,
false)
1280 newToOldPoints.setSize(
s.points().size());
1281 oldToNewPoints.setSize(
s.points().size());
1282 oldToNewPoints = -1;
1286 forAll(include, oldFacei)
1288 if (include[oldFacei])
1291 const labelledTri& tri =
s[oldFacei];
1295 label oldPointi = tri[fp];
1297 if (oldToNewPoints[oldPointi] == -1)
1299 oldToNewPoints[oldPointi] = pointi;
1300 newToOldPoints[pointi++] = oldPointi;
1305 newToOldPoints.setSize(pointi);
1310 forAll(newToOldPoints, i)
1312 newPoints[i] =
s.points()[newToOldPoints[i]];
1315 List<labelledTri> newTriangles(newToOldFaces.size());
1320 const labelledTri& tri =
s[newToOldFaces[i]];
1322 newTriangles[i][0] = oldToNewPoints[tri[0]];
1323 newTriangles[i][1] = oldToNewPoints[tri[1]];
1324 newTriangles[i][2] = oldToNewPoints[tri[2]];
1325 newTriangles[i].region() = tri.region();
1329 return triSurface(newTriangles,
s.patches(), newPoints,
true);
1353 mergeDistance_(params.
mergeTol()*mesh_.bounds().mag())
1355 const bool regularise = (params.
filter() != filterType::NONE);
1356 const fvMesh& fvmesh = cellValues.mesh();
1360 Pout<<
"isoSurfacePoint:" <<
nl
1361 <<
" isoField : " << cellValues.name() <<
nl
1362 <<
" cell min/max : " <<
minMax(cVals_) <<
nl
1363 <<
" point min/max : " <<
minMax(pVals_) <<
nl
1364 <<
" isoValue : " << iso <<
nl
1365 <<
" filter : " <<
Switch(regularise) <<
nl
1377 cValsPtr_.reset(adaptPatchFields(cellValues).ptr());
1411 meshC.boundaryField()[patchi]
1416 collocatedFaces(refCast<const coupledPolyPatch>(pp))
1421 if (!isCollocated[i])
1423 pfld[i] = mesh_.faceCentres()[pp.
start()+i];
1427 else if (isA<emptyPolyPatch>(pp))
1429 auto& bfld =
const_cast<slicedVolVectorField::Boundary&
>
1431 meshC.boundaryField()
1435 bfld.set(patchi,
nullptr);
1449 bfld[patchi] = pp.
patchSlice(mesh_.faceCentres());
1472 "isoSurfacePoint.cutType",
1485 forAll(cellCutType_, celli)
1487 debugFld[celli] = cellCutType_[celli];
1490 Pout<<
"Writing cut types:"
1491 << debugField.objectPath() <<
endl;
1516 snappedCc.
setSize(mesh_.nCells());
1524 Pout<<
"isoSurfacePoint : shifted " << snappedPoints.size()
1525 <<
" cell centres to intersection." <<
endl;
1528 label nCellSnaps = snappedPoints.size();
1536 bitSet isBoundaryPoint(mesh_.nPoints());
1546 refCast<const coupledPolyPatch>(pp);
1548 bitSet isCollocated(collocatedFaces(cpp));
1552 if (!isCollocated[i])
1554 const face&
f = mesh_.faces()[cpp.
start()+i];
1556 isBoundaryPoint.
set(
f);
1564 const face&
f = mesh_.faces()[pp.start()+i];
1566 isBoundaryPoint.
set(
f);
1585 snappedPoint.
setSize(mesh_.nPoints());
1591 Pout<<
"isoSurfacePoint : shifted " << snappedPoints.size()-nCellSnaps
1592 <<
" vertices to intersection." <<
endl;
1621 Pout<<
"isoSurfacePoint : generated " << triMeshCells.size()
1623 <<
" unmerged points." <<
endl;
1632 if (getClipBounds().valid())
1640 interpolatedPoints_,
1641 interpolatedOldPoints_,
1642 interpolationWeights_
1644 triMeshCells =
labelField(triMeshCells, trimTriMap);
1650 tmpsurf = stitchTriPoints
1660 Pout<<
"isoSurfacePoint : generated " << triMap.size()
1661 <<
" merged triangles." <<
endl;
1665 if (getClipBounds().valid())
1671 labelList newTriPointMergeMap(nOldPoints, -1);
1672 forAll(trimTriPointMap, trimPointI)
1674 label oldPointI = trimTriPointMap[trimPointI];
1677 label pointI = triPointMergeMap_[trimPointI];
1680 newTriPointMergeMap[oldPointI] = pointI;
1684 triPointMergeMap_.
transfer(newTriPointMergeMap);
1687 meshCells_.setSize(triMap.size());
1690 meshCells_[i] = triMeshCells[triMap[i]];
1696 Pout<<
"isoSurfacePoint : checking " << tmpsurf.size()
1697 <<
" triangles for validity." <<
endl;
1705 Pout<<
"Dumping surface to " << stlFile <<
endl;
1706 tmpsurf.
write(stlFile);
1724 this->Mesh::transfer(updated);