73 list.append(std::move(tri));
81 const auto* cpp = isA<coupledPolyPatch>(pp);
82 return 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();
161 const labelList& procPatches = mesh_.globalData().processorPatches();
164 for (
const label patchi : procPatches)
166 const polyPatch& pp =
patches[patchi];
167 const auto& procPatch = refCast<const processorPolyPatch>(pp);
171 const labelList& meshPts = procPatch.meshPoints();
172 const labelList& nbrPts = procPatch.neighbPoints();
178 const label nbrPointi = nbrPts[pointi];
179 patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
185 procPatch.neighbProcNo()
192 for (
const label patchi : procPatches)
194 const polyPatch& pp =
patches[patchi];
195 const auto& procPatch = refCast<const processorPolyPatch>(pp);
206 procPatch.neighbProcNo()
208 fromNbr >> nbrPatchInfo;
211 const labelList& meshPts = procPatch.meshPoints();
215 const label meshPointi = meshPts[pointi];
218 pointValues[meshPointi],
227 for (
const polyPatch& pp :
patches)
229 const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
235 const auto& cycPatch = *cpp;
236 const auto& nbrPatch = cycPatch.neighbPatch();
238 const edgeList& coupledPoints = cycPatch.coupledPoints();
239 const labelList& meshPts = cycPatch.meshPoints();
240 const labelList& nbrMeshPoints = nbrPatch.meshPoints();
247 const edge&
e = coupledPoints[i];
248 half0Values[i] = pointValues[meshPts[
e[0]]];
249 half1Values[i] = pointValues[nbrMeshPoints[
e[1]]];
254 const edge&
e = coupledPoints[i];
255 const label
p0 = meshPts[
e[0]];
256 const label p1 = nbrMeshPoints[
e[1]];
258 minEqOp<point>()(pointValues[
p0], half1Values[i]);
259 minEqOp<point>()(pointValues[p1], half0Values[i]);
265 const globalMeshData& pd = mesh_.globalData();
267 if (pd.nGlobalPoints() > 0)
270 pointField sharedPts(pd.nGlobalPoints(), nullValue);
272 forAll(pd.sharedPointLabels(), i)
274 const label meshPointi = pd.sharedPointLabels()[i];
276 sharedPts[pd.sharedPointAddr()[i]] = pointValues[meshPointi];
285 forAll(pd.sharedPointLabels(), i)
287 const label meshPointi = pd.sharedPointLabels()[i];
288 pointValues[meshPointi] = sharedPts[pd.sharedPointAddr()[i]];
294 Foam::scalar Foam::isoSurfacePoint::isoFraction
300 const scalar d = s1-s0;
311 bool Foam::isoSurfacePoint::isEdgeOfFaceCut
321 const bool fpLower = (pVals[
f[fp]] < iso_);
326 || fpLower != neiLower
327 || fpLower != (pVals[
f[
f.fcIndex(fp)]] < iso_)
338 void Foam::isoSurfacePoint::getNeighbour
349 const labelList& own = mesh_.faceOwner();
350 const labelList& nei = mesh_.faceNeighbour();
352 if (mesh_.isInternalFace(facei))
354 const label nbr = (own[facei] == celli ? nei[facei] : own[facei]);
355 nbrValue = cVals[nbr];
356 nbrPoint = meshC[nbr];
360 const label bFacei = facei-mesh_.nInternalFaces();
361 const label patchi = boundaryRegion[bFacei];
362 const polyPatch& pp = mesh_.boundaryMesh()[patchi];
363 const label patchFacei = facei-pp.start();
365 nbrValue = cVals.boundaryField()[patchi][patchFacei];
366 nbrPoint = meshC.boundaryField()[patchi][patchFacei];
371 void Foam::isoSurfacePoint::calcCutTypes
379 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
380 const labelList& own = mesh_.faceOwner();
381 const labelList& nei = mesh_.faceNeighbour();
383 faceCutType_.
resize(mesh_.nFaces());
384 faceCutType_ = cutType::NOTCUT;
386 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
388 const scalar ownValue = cVals[own[facei]];
389 const bool ownLower = (ownValue < iso_);
404 const bool neiLower = (nbrValue < iso_);
406 if (ownLower != neiLower)
408 faceCutType_[facei] = cutType::CUT;
414 const face&
f = mesh_.faces()[facei];
416 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
418 faceCutType_[facei] = cutType::CUT;
423 for (
const polyPatch& pp :
patches)
425 for (
const label facei : pp.range())
427 const scalar ownValue = cVals[own[facei]];
428 const bool ownLower = (ownValue < iso_);
443 const bool neiLower = (nbrValue < iso_);
445 if (ownLower != neiLower)
447 faceCutType_[facei] = cutType::CUT;
453 const face&
f = mesh_.faces()[facei];
455 if (isEdgeOfFaceCut(pVals,
f, ownLower, neiLower))
457 faceCutType_[facei] = cutType::CUT;
464 cellCutType_.
resize(mesh_.nCells());
465 cellCutType_ = cutType::NOTCUT;
470 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
472 if (faceCutType_[facei] == cutType::NOTCUT)
477 if (cellCutType_[own[facei]] == cutType::NOTCUT)
479 cellCutType_[own[facei]] = cutType::CUT;
482 if (cellCutType_[nei[facei]] == cutType::NOTCUT)
484 cellCutType_[nei[facei]] = cutType::CUT;
492 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
494 if (faceCutType_[facei] == cutType::NOTCUT)
499 if (cellCutType_[own[facei]] == cutType::NOTCUT)
501 cellCutType_[own[facei]] = cutType::CUT;
508 Pout<<
"isoSurfacePoint : candidate cut cells "
509 << nCutCells_ <<
" / " << mesh_.nCells() <<
endl;
514 Foam::point Foam::isoSurfacePoint::calcCentre(
const triSurface&
s)
522 sum +=
s[i].centre(
s.points());
528 void Foam::isoSurfacePoint::calcSnappedCc
535 DynamicList<point>& snappedPoints,
542 snappedCc.setSize(mesh_.nCells());
546 DynamicList<point, 64> localTriPoints(64);
548 forAll(mesh_.cells(), celli)
550 if (cellCutType_[celli] == cutType::CUT)
552 const scalar cVal = cVals[celli];
554 localTriPoints.clear();
561 for (
const label facei : mesh_.cells()[celli])
563 const face&
f = mesh_.faces()[facei];
579 FixedList<scalar, 3>
s;
580 FixedList<point, 3> pt;
583 s[2] = isoFraction(cVal, nbrValue);
584 pt[2] = (1.0-
s[2])*cc[celli] +
s[2]*nbrPoint;
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)
641 snappedCc[celli] = snappedPoints.size();
642 snappedPoints.append(
sum(localTriPoints)/3);
643 localTriPoints.clear();
668 label
nZones = surf.markZones
676 snappedCc[celli] = snappedPoints.size();
677 snappedPoints.append(calcCentre(surf));
688 void Foam::isoSurfacePoint::calcSnappedPoint
690 const bitSet& isBoundaryPoint,
696 DynamicList<point>& snappedPoints,
706 DynamicList<point, 64> localTriPoints(100);
708 forAll(mesh_.pointFaces(), pointi)
710 if (isBoundaryPoint.test(pointi))
719 for (
const label facei :
pFaces)
721 if (faceCutType_[facei] == cutType::CUT)
734 localTriPoints.clear();
738 for (
const label facei :
pFaces)
742 const face&
f = mesh_.faces()[facei];
743 const label own = mesh_.faceOwner()[facei];
760 FixedList<scalar, 4>
s;
761 FixedList<point, 4> pt;
763 label fp =
f.find(pointi);
764 s[0] = isoFraction(pVals[pointi], cVals[own]);
765 pt[0] = (1.0-
s[0])*pts[pointi] +
s[0]*cc[own];
767 s[1] = isoFraction(pVals[pointi], nbrValue);
768 pt[1] = (1.0-
s[1])*pts[pointi] +
s[1]*nbrPoint;
770 label nextPointi =
f[
f.fcIndex(fp)];
771 s[2] = isoFraction(pVals[pointi], pVals[nextPointi]);
772 pt[2] = (1.0-
s[2])*pts[pointi] +
s[2]*pts[nextPointi];
774 label prevPointi =
f[
f.rcIndex(fp)];
775 s[3] = isoFraction(pVals[pointi], pVals[prevPointi]);
776 pt[3] = (1.0-
s[3])*pts[pointi] +
s[3]*pts[prevPointi];
780 (
s[0] >= 0.0 &&
s[0] <= 0.5)
781 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
782 && (
s[2] >= 0.0 &&
s[2] <= 0.5)
785 localTriPoints.append(pt[0]);
786 localTriPoints.append(pt[1]);
787 localTriPoints.append(pt[2]);
791 (
s[0] >= 0.0 &&
s[0] <= 0.5)
792 && (
s[1] >= 0.0 &&
s[1] <= 0.5)
793 && (
s[3] >= 0.0 &&
s[3] <= 0.5)
796 localTriPoints.append(pt[3]);
797 localTriPoints.append(pt[0]);
798 localTriPoints.append(pt[1]);
804 if (
s[i] >= 0.0 &&
s[i] <= 0.5)
806 otherPointSum += pt[i];
812 if (localTriPoints.size() == 0)
818 collapsedPoint[pointi] = otherPointSum/nOther;
821 else if (localTriPoints.size() == 3)
825 points.transfer(localTriPoints);
847 label
nZones = surf.markZones
855 collapsedPoint[pointi] = calcCentre(surf);
862 syncUnseparatedPoints(collapsedPoint,
point::max);
865 snappedPoint.setSize(mesh_.nPoints());
868 forAll(collapsedPoint, pointi)
872 snappedPoint[pointi] = snappedPoints.size();
873 snappedPoints.append(collapsedPoint[pointi]);
881 const bool checkDuplicates,
882 const List<point>& triPoints,
887 label nTris = triPoints.size()/3;
889 if ((triPoints.size() % 3) != 0)
892 <<
"Problem: number of points " << triPoints.size()
909 Pout<<
"isoSurfacePoint : merged from " << triPoints.size()
910 <<
" down to " << newPoints.size() <<
" unique points." <<
endl;
926 <<
"Merged points contain duplicates"
927 <<
" when merging with distance " << mergeDistance_ <<
endl
928 <<
"merged:" << newPoints.size() <<
" re-merged:"
929 << newNewPoints.size()
935 List<labelledTri> tris;
937 DynamicList<labelledTri> dynTris(nTris);
939 DynamicList<label> newToOldTri(nTris);
941 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
945 triPointReverseMap[rawPointi],
946 triPointReverseMap[rawPointi+1],
947 triPointReverseMap[rawPointi+2],
952 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
954 newToOldTri.append(oldTriI);
959 triMap.transfer(newToOldTri);
960 tris.transfer(dynTris);
975 DynamicList<label> newToOldTri(tris.size());
979 const labelledTri& tri = tris[triI];
988 label nbrTriI =
pFaces[i];
990 if (nbrTriI > triI && (tris[nbrTriI] == tri))
1003 label newTriI = newToOldTri.size();
1004 newToOldTri.append(triMap[triI]);
1005 tris[newTriI] = tris[triI];
1009 triMap.transfer(newToOldTri);
1010 tris.setSize(triMap.size());
1014 Pout<<
"isoSurfacePoint : merged from " << nTris
1015 <<
" down to " << tris.size() <<
" unique triangles." <<
endl;
1025 const labelledTri&
f = surf[facei];
1026 const labelList& fFaces = surf.faceFaces()[facei];
1030 label nbrFacei = fFaces[i];
1032 if (nbrFacei <= facei)
1038 const labelledTri& nbrF = surf[nbrFacei];
1044 <<
" triangle " << facei <<
" vertices " <<
f
1045 <<
" fc:" <<
f.centre(surf.points())
1046 <<
" has the same vertices as triangle " << nbrFacei
1047 <<
" vertices " << nbrF
1048 <<
" fc:" << nbrF.centre(surf.points())
1060 void Foam::isoSurfacePoint::trimToPlanes
1062 const PtrList<plane>& planes,
1064 DynamicList<point>& newTriPoints
1068 DynamicList<triPoints> insideTrisA;
1069 storeOp insideOpA(insideTrisA);
1072 DynamicList<triPoints> insideTrisB;
1073 storeOp insideOpB(insideTrisB);
1075 triPointRef::dummyOp dop;
1078 insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
1085 const plane& pln = planes[faceI];
1089 insideTrisB.clear();
1090 for (
const triPoints& tri : insideTrisA)
1092 triPointRef(tri).sliceWithPlane(pln, insideOpB, dop);
1097 insideTrisA.clear();
1098 for (
const triPoints& tri : insideTrisB)
1100 triPointRef(tri).sliceWithPlane(pln, insideOpA, dop);
1107 DynamicList<triPoints>& newTris = (useA ? insideTrisA : insideTrisB);
1109 newTriPoints.reserve(newTriPoints.size() + 3*newTris.size());
1112 for (
const triPoints& tri : newTris)
1114 newTriPoints.append(tri[0]);
1115 newTriPoints.append(tri[1]);
1116 newTriPoints.append(tri[2]);
1121 void Foam::isoSurfacePoint::trimToBox
1123 const treeBoundBox& bb,
1124 DynamicList<point>& triPoints,
1125 DynamicList<label>& triMap
1130 Pout<<
"isoSurfacePoint : trimming to " << bb <<
endl;
1138 planes.set(faceI,
new plane(bb.faceCentre(faceI), -
n));
1141 label nTris = triPoints.size()/3;
1143 DynamicList<point> newTriPoints(triPoints.size()/16);
1144 triMap.setCapacity(nTris/16);
1147 for (label triI = 0; triI < nTris; triI++)
1149 const point&
p0 = triPoints[vertI++];
1150 const point& p1 = triPoints[vertI++];
1151 const point& p2 = triPoints[vertI++];
1153 label oldNPoints = newTriPoints.
size();
1161 label nCells = (newTriPoints.size()-oldNPoints)/3;
1162 for (label i = 0; i < nCells; i++)
1164 triMap.append(triI);
1170 Pout<<
"isoSurfacePoint : trimmed from " << nTris
1171 <<
" down to " << triMap.size()
1172 <<
" triangles." <<
endl;
1175 triPoints.transfer(newTriPoints);
1179 void Foam::isoSurfacePoint::trimToBox
1181 const treeBoundBox& bb,
1182 DynamicList<point>& triPoints,
1183 DynamicList<label>& triMap,
1186 List<FixedList<label, 3>>& interpolatedOldPoints,
1187 List<FixedList<scalar, 3>>& interpolationWeights
1190 const List<point> oldTriPoints(triPoints);
1193 trimToBox(bb, triPoints, triMap);
1200 label sz = oldTriPoints.size()/100;
1201 DynamicList<label> dynInterpolatedPoints(sz);
1202 DynamicList<FixedList<label, 3>> dynInterpolatedOldPoints(sz);
1203 DynamicList<FixedList<scalar, 3>> dynInterpolationWeights(sz);
1206 triPointMap.setSize(triPoints.size());
1209 label oldTriI = triMap[triI];
1212 for (label i = 0; i < 3; i++)
1214 label pointI = 3*triI+i;
1215 const point& pt = triPoints[pointI];
1218 label matchPointI = -1;
1219 for (label j = 0; j < 3; j++)
1221 label oldPointI = 3*oldTriI+j;
1222 if (pt == oldTriPoints[oldPointI])
1224 matchPointI = oldPointI;
1229 triPointMap[pointI] = matchPointI;
1232 if (matchPointI == -1)
1234 dynInterpolatedPoints.append(pointI);
1236 FixedList<label, 3> oldPoints
1238 {3*oldTriI, 3*oldTriI+1, 3*oldTriI+2}
1240 dynInterpolatedOldPoints.append(oldPoints);
1244 FixedList<scalar, 3> weights({bary.a(), bary.b(), bary.c()});
1246 dynInterpolationWeights.append(weights);
1251 interpolatedPoints.transfer(dynInterpolatedPoints);
1252 interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
1253 interpolationWeights.transfer(dynInterpolationWeights);
1259 const triSurface&
s,
1267 ListOps::createWithValue<bool>(
s.size(), newToOldFaces,
true,
false)
1270 newToOldPoints.setSize(
s.points().size());
1271 oldToNewPoints.setSize(
s.points().size());
1272 oldToNewPoints = -1;
1276 forAll(include, oldFacei)
1278 if (include[oldFacei])
1281 const labelledTri& tri =
s[oldFacei];
1285 label oldPointi = tri[fp];
1287 if (oldToNewPoints[oldPointi] == -1)
1289 oldToNewPoints[oldPointi] = pointi;
1290 newToOldPoints[pointi++] = oldPointi;
1295 newToOldPoints.setSize(pointi);
1300 forAll(newToOldPoints, i)
1302 newPoints[i] =
s.points()[newToOldPoints[i]];
1305 List<labelledTri> newTriangles(newToOldFaces.size());
1310 const labelledTri& tri =
s[newToOldFaces[i]];
1312 newTriangles[i][0] = oldToNewPoints[tri[0]];
1313 newTriangles[i][1] = oldToNewPoints[tri[1]];
1314 newTriangles[i][2] = oldToNewPoints[tri[2]];
1315 newTriangles[i].region() = tri.region();
1319 return triSurface(newTriangles,
s.patches(), newPoints,
true);
1343 mergeDistance_(params.
mergeTol()*mesh_.bounds().mag())
1345 const bool regularise = (params.
filter() != filterType::NONE);
1346 const fvMesh& fvmesh = cellValues.mesh();
1350 Pout<<
"isoSurfacePoint:" <<
nl
1351 <<
" isoField : " << cellValues.name() <<
nl
1352 <<
" cell min/max : " <<
minMax(cVals_) <<
nl
1353 <<
" point min/max : " <<
minMax(pVals_) <<
nl
1354 <<
" isoValue : " << iso <<
nl
1355 <<
" filter : " <<
Switch(regularise) <<
nl
1367 cValsPtr_.reset(adaptPatchFields(cellValues).ptr());
1401 meshC.boundaryField()[patchi]
1406 collocatedFaces(refCast<const coupledPolyPatch>(pp))
1411 if (!isCollocated[i])
1413 pfld[i] = mesh_.faceCentres()[pp.
start()+i];
1417 else if (isA<emptyPolyPatch>(pp))
1419 auto& bfld =
const_cast<slicedVolVectorField::Boundary&
>
1421 meshC.boundaryField()
1425 bfld.set(patchi,
nullptr);
1439 bfld[patchi] = pp.
patchSlice(mesh_.faceCentres());
1462 "isoSurfacePoint.cutType",
1475 forAll(cellCutType_, celli)
1477 debugFld[celli] = cellCutType_[celli];
1480 Pout<<
"Writing cut types:"
1481 << debugField.objectPath() <<
endl;
1506 snappedCc.
setSize(mesh_.nCells());
1514 Pout<<
"isoSurfacePoint : shifted " << snappedPoints.size()
1515 <<
" cell centres to intersection." <<
endl;
1518 label nCellSnaps = snappedPoints.size();
1526 bitSet isBoundaryPoint(mesh_.nPoints());
1536 refCast<const coupledPolyPatch>(pp);
1538 bitSet isCollocated(collocatedFaces(cpp));
1542 if (!isCollocated[i])
1544 const face&
f = mesh_.faces()[cpp.
start()+i];
1546 isBoundaryPoint.
set(
f);
1554 const face&
f = mesh_.faces()[pp.start()+i];
1556 isBoundaryPoint.
set(
f);
1575 snappedPoint.
setSize(mesh_.nPoints());
1581 Pout<<
"isoSurfacePoint : shifted " << snappedPoints.size()-nCellSnaps
1582 <<
" vertices to intersection." <<
endl;
1611 Pout<<
"isoSurfacePoint : generated " << triMeshCells.size()
1613 <<
" unmerged points." <<
endl;
1622 if (getClipBounds().valid())
1630 interpolatedPoints_,
1631 interpolatedOldPoints_,
1632 interpolationWeights_
1634 triMeshCells =
labelField(triMeshCells, trimTriMap);
1640 tmpsurf = stitchTriPoints
1650 Pout<<
"isoSurfacePoint : generated " << triMap.size()
1651 <<
" merged triangles." <<
endl;
1655 if (getClipBounds().valid())
1661 labelList newTriPointMergeMap(nOldPoints, -1);
1662 forAll(trimTriPointMap, trimPointI)
1664 label oldPointI = trimTriPointMap[trimPointI];
1667 label pointI = triPointMergeMap_[trimPointI];
1670 newTriPointMergeMap[oldPointI] = pointI;
1674 triPointMergeMap_.
transfer(newTriPointMergeMap);
1677 meshCells_.setSize(triMap.size());
1680 meshCells_[i] = triMeshCells[triMap[i]];
1686 Pout<<
"isoSurfacePoint : checking " << tmpsurf.size()
1687 <<
" triangles for validity." <<
endl;
1695 Pout<<
"Dumping surface to " << stlFile <<
endl;
1696 tmpsurf.
write(stlFile);
1714 this->Mesh::transfer(updated);