46 bool Foam::isoSurface::noTransform(
const tensor& tt)
const
49 (
mag(tt.xx()-1) < mergeDistance_)
50 && (
mag(tt.yy()-1) < mergeDistance_)
51 && (
mag(tt.zz()-1) < mergeDistance_)
52 && (
mag(tt.xy()) < mergeDistance_)
53 && (
mag(tt.xz()) < mergeDistance_)
54 && (
mag(tt.yx()) < mergeDistance_)
55 && (
mag(tt.yz()) < mergeDistance_)
56 && (
mag(tt.zx()) < mergeDistance_)
57 && (
mag(tt.zy()) < mergeDistance_);
61 bool Foam::isoSurface::collocatedPatch(
const polyPatch& pp)
63 const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
65 return cpp.parallel() && !cpp.separated();
72 const coupledPolyPatch& pp
81 if (forwardT.size() == 0)
84 if (separation.size() == 0)
88 else if (separation.size() == 1)
97 if (
mag(separation[faceI]) < mergeDistance_)
99 collocated[faceI] = 1u;
104 else if (forwardT.size() == 1)
113 if (noTransform(forwardT[faceI]))
115 collocated[faceI] = 1u;
125 void Foam::isoSurface::insertPointData
127 const processorPolyPatch& pp,
128 const Map<label>& meshToShared,
130 const label patchPointI,
135 label meshPointI = pp.meshPoints()[patchPointI];
138 label nbrPointI = pp.neighbPoints()[patchPointI];
139 if (nbrPointI >= 0 && nbrPointI < patchValues.size())
141 minEqOp<point>()(patchValues[nbrPointI], pointValues[meshPointI]);
146 if (iter != meshToShared.end())
148 minEqOp<point>()(sharedValues[iter()], pointValues[meshPointI]);
153 void Foam::isoSurface::syncUnseparatedPoints
156 const point& nullValue
161 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
162 const globalMeshData& pd = mesh_.globalData();
163 const labelList& sharedPtAddr = pd.sharedPointAddr();
164 const labelList& sharedPtLabels = pd.sharedPointLabels();
167 Map<label> meshToShared(2*sharedPtLabels.size());
170 meshToShared.insert(sharedPtLabels[i], sharedPtAddr[i]);
174 pointField sharedInfo(pd.nGlobalPoints(), nullValue);
184 isA<processorPolyPatch>(patches[patchI])
185 && patches[patchI].
nPoints() > 0
188 const processorPolyPatch& pp =
189 refCast<const processorPolyPatch>(patches[patchI]);
190 const labelList& meshPts = pp.meshPoints();
192 pointField patchInfo(meshPts.size(), nullValue);
196 forAll(isCollocated, faceI)
198 if (isCollocated[faceI])
200 const face&
f = pp.localFaces()[faceI];
204 label pointI = f[fp];
230 isA<processorPolyPatch>(patches[patchI])
231 && patches[patchI].
nPoints() > 0
234 const processorPolyPatch& pp =
235 refCast<const processorPolyPatch>(patches[patchI]);
242 fromNbr >> nbrPatchInfo;
246 nbrPatchInfo.
setSize(pp.nPoints(), nullValue);
248 const labelList& meshPts = pp.meshPoints();
252 label meshPointI = meshPts[pointI];
255 pointValues[meshPointI],
277 label meshPointI = sharedPtLabels[i];
278 pointValues[meshPointI] = sharedInfo[sharedPtAddr[i]];
283 Foam::scalar Foam::isoSurface::isoFraction
302 bool Foam::isoSurface::isEdgeOfFaceCut
312 bool fpLower = (pVals[f[fp]] < iso_);
315 (fpLower != ownLower)
316 || (fpLower != neiLower)
317 || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
328 void Foam::isoSurface::getNeighbour
339 const labelList& own = mesh_.faceOwner();
340 const labelList& nei = mesh_.faceNeighbour();
342 if (mesh_.isInternalFace(faceI))
344 label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
345 nbrValue = cVals[nbr];
346 nbrPoint = meshC[nbr];
350 label bFaceI = faceI-mesh_.nInternalFaces();
351 label patchI = boundaryRegion[bFaceI];
352 const polyPatch& pp = mesh_.boundaryMesh()[patchI];
353 label patchFaceI = faceI-pp.start();
355 nbrValue = cVals.boundaryField()[patchI][patchFaceI];
356 nbrPoint = meshC.boundaryField()[patchI][patchFaceI];
362 void Foam::isoSurface::calcCutTypes
370 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
371 const labelList& own = mesh_.faceOwner();
372 const labelList& nei = mesh_.faceNeighbour();
374 faceCutType_.
setSize(mesh_.nFaces());
375 faceCutType_ = NOTCUT;
377 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
380 bool ownLower = (cVals[own[faceI]] < iso_);
395 bool neiLower = (nbrValue < iso_);
397 if (ownLower != neiLower)
399 faceCutType_[faceI] = CUT;
405 const face f = mesh_.faces()[faceI];
407 if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
409 faceCutType_[faceI] = CUT;
416 const polyPatch& pp = patches[patchI];
418 label faceI = pp.start();
422 bool ownLower = (cVals[own[faceI]] < iso_);
437 bool neiLower = (nbrValue < iso_);
439 if (ownLower != neiLower)
441 faceCutType_[faceI] = CUT;
446 const face f = mesh_.faces()[faceI];
448 if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
450 faceCutType_[faceI] = CUT;
461 cellCutType_.setSize(mesh_.nCells());
462 cellCutType_ = NOTCUT;
464 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
466 if (faceCutType_[faceI] != NOTCUT)
468 if (cellCutType_[own[faceI]] == NOTCUT)
470 cellCutType_[own[faceI]] = CUT;
473 if (cellCutType_[nei[faceI]] == NOTCUT)
475 cellCutType_[nei[faceI]] = CUT;
480 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
482 if (faceCutType_[faceI] != NOTCUT)
484 if (cellCutType_[own[faceI]] == NOTCUT)
486 cellCutType_[own[faceI]] = CUT;
494 Pout<<
"isoSurface : detected " << nCutCells_
495 <<
" candidate cut cells (out of " << mesh_.nCells()
504 const labelledTri& tri0,
505 const labelledTri& tri1
524 label fp0p1 = tri0.fcIndex(fp0);
525 label fp1p1 = tri1.fcIndex(fp1);
526 label fp1m1 = tri1.rcIndex(fp1);
528 if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
530 common[0] = tri0[fp0];
531 common[1] = tri0[fp0p1];
539 Foam::point Foam::isoSurface::calcCentre(
const triSurface& s)
545 sum += s[i].centre(s.points());
556 DynamicList<labelledTri, 64>& localTris
561 if (localTris.size() == 1)
563 const labelledTri& tri = localTris[0];
564 info.setPoint(tri.centre(localPoints));
567 else if (localTris.size() == 2)
570 const labelledTri& tri0 = localTris[0];
571 const labelledTri& tri1 = localTris[0];
573 labelPair shared = findCommonPoints(tri0, tri1);
581 tri0.centre(localPoints)
582 + tri1.centre(localPoints)
588 else if (localTris.size())
598 localTris.clearStorage();
601 label nZones = surf.markZones
609 info.setPoint(calcCentre(surf));
620 void Foam::isoSurface::calcSnappedCc
627 DynamicList<point>& snappedPoints,
634 snappedCc.
setSize(mesh_.nCells());
638 DynamicList<point, 64> localTriPoints(64);
640 forAll(mesh_.cells(), cellI)
642 if (cellCutType_[cellI] == CUT)
644 scalar cVal = cVals[cellI];
646 const cell& cFaces = mesh_.cells()[cellI];
648 localTriPoints.clear();
657 label faceI = cFaces[cFaceI];
673 FixedList<scalar, 3> s;
674 FixedList<point, 3> pt;
677 s[2] = isoFraction(cVal, nbrValue);
678 pt[2] = (1.0-s[2])*cc[cellI] + s[2]*nbrPoint;
680 const face& f = mesh_.faces()[cFaces[cFaceI]];
686 s[0] = isoFraction(cVal, pVals[p0]);
687 pt[0] = (1.0-s[0])*cc[cellI] + s[0]*pts[p0];
690 label p1 = f[f.fcIndex(fp)];
691 s[1] = isoFraction(cVal, pVals[p1]);
692 pt[1] = (1.0-s[1])*cc[cellI] + s[1]*pts[p1];
696 (s[0] >= 0.0 && s[0] <= 0.5)
697 && (s[1] >= 0.0 && s[1] <= 0.5)
698 && (s[2] >= 0.0 && s[2] <= 0.5)
701 localTriPoints.append(pt[0]);
702 localTriPoints.append(pt[1]);
703 localTriPoints.append(pt[2]);
710 if (s[i] >= 0.0 && s[i] <= 0.5)
712 otherPointSum += pt[i];
720 if (localTriPoints.size() == 0)
726 snappedCc[cellI] = snappedPoints.size();
727 snappedPoints.append(otherPointSum/nOther);
734 else if (localTriPoints.size() == 3)
739 snappedCc[cellI] = snappedPoints.size();
740 snappedPoints.append(
sum(points)/points.size());
765 label nZones = surf.markZones
773 snappedCc[cellI] = snappedPoints.
size();
774 snappedPoints.append(calcCentre(surf));
787 void Foam::isoSurface::calcSnappedPoint
795 DynamicList<point>& snappedPoints,
803 const point greatPoint(VGREAT, VGREAT, VGREAT);
804 pointField collapsedPoint(mesh_.nPoints(), greatPoint);
808 DynamicList<point, 64> localTriPoints(100);
810 forAll(mesh_.pointFaces(), pointI)
812 if (isBoundaryPoint.get(pointI) == 1)
823 label faceI = pFaces[i];
825 if (faceCutType_[faceI] == CUT)
838 localTriPoints.clear();
847 label faceI = pFaces[pFaceI];
848 const face& f = mesh_.faces()[faceI];
849 label own = mesh_.faceOwner()[faceI];
866 FixedList<scalar, 4> s;
867 FixedList<point, 4> pt;
870 s[0] = isoFraction(pVals[pointI], cVals[own]);
871 pt[0] = (1.0-s[0])*pts[pointI] + s[0]*cc[own];
873 s[1] = isoFraction(pVals[pointI], nbrValue);
874 pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint;
876 label nextPointI = f[f.fcIndex(fp)];
877 s[2] = isoFraction(pVals[pointI], pVals[nextPointI]);
878 pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI];
880 label prevPointI = f[f.rcIndex(fp)];
881 s[3] = isoFraction(pVals[pointI], pVals[prevPointI]);
882 pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI];
886 (s[0] >= 0.0 && s[0] <= 0.5)
887 && (s[1] >= 0.0 && s[1] <= 0.5)
888 && (s[2] >= 0.0 && s[2] <= 0.5)
891 localTriPoints.append(pt[0]);
892 localTriPoints.append(pt[1]);
893 localTriPoints.append(pt[2]);
897 (s[0] >= 0.0 && s[0] <= 0.5)
898 && (s[1] >= 0.0 && s[1] <= 0.5)
899 && (s[3] >= 0.0 && s[3] <= 0.5)
902 localTriPoints.append(pt[3]);
903 localTriPoints.append(pt[0]);
904 localTriPoints.append(pt[1]);
910 if (s[i] >= 0.0 && s[i] <= 0.5)
912 otherPointSum += pt[i];
918 if (localTriPoints.size() == 0)
924 collapsedPoint[pointI] = otherPointSum/nOther;
927 else if (localTriPoints.size() == 3)
932 collapsedPoint[pointI] =
sum(points)/points.size();
953 label nZones = surf.markZones
961 collapsedPoint[pointI] = calcCentre(surf);
968 syncUnseparatedPoints(collapsedPoint, greatPoint);
971 snappedPoint.setSize(mesh_.nPoints());
974 forAll(collapsedPoint, pointI)
976 if (collapsedPoint[pointI] != greatPoint)
978 snappedPoint[pointI] = snappedPoints.size();
979 snappedPoints.append(collapsedPoint[pointI]);
987 const bool checkDuplicates,
988 const List<point>& triPoints,
993 label nTris = triPoints.
size()/3;
995 if ((triPoints.size() % 3) != 0)
998 <<
"Problem: number of points " << triPoints.size()
1029 <<
"Merged points contain duplicates"
1030 <<
" when merging with distance " << mergeDistance_ <<
endl
1031 <<
"merged:" << newPoints.size() <<
" re-merged:"
1032 << newNewPoints.size()
1038 List<labelledTri> tris;
1040 DynamicList<labelledTri> dynTris(nTris);
1041 label rawPointI = 0;
1042 DynamicList<label> newToOldTri(nTris);
1044 for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
1048 triPointReverseMap[rawPointI],
1049 triPointReverseMap[rawPointI+1],
1050 triPointReverseMap[rawPointI+2],
1055 if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
1057 newToOldTri.append(oldTriI);
1058 dynTris.append(tri);
1062 triMap.transfer(newToOldTri);
1063 tris.transfer(dynTris);
1072 if (checkDuplicates)
1078 DynamicList<label> newToOldTri(tris.size());
1082 const labelledTri& tri = tris[triI];
1083 const labelList& pFaces = pointFaces[tri[0]];
1091 label nbrTriI = pFaces[i];
1093 if (nbrTriI > triI && (tris[nbrTriI] == tri))
1106 label newTriI = newToOldTri.size();
1107 newToOldTri.append(triI);
1108 tris[newTriI] = tris[triI];
1112 triMap.transfer(newToOldTri);
1113 tris.setSize(triMap.size());
1117 Pout<<
"isoSurface : merged from " << nTris
1118 <<
" down to " << tris.size() <<
" unique triangles." <<
endl;
1128 const labelledTri& f = surf[faceI];
1129 const labelList& fFaces = surf.faceFaces()[faceI];
1133 label nbrFaceI = fFaces[i];
1135 if (nbrFaceI <= faceI)
1141 const labelledTri& nbrF = surf[nbrFaceI];
1145 FatalErrorIn(
"validTri(const triSurface&, const label)")
1147 <<
" triangle " << faceI <<
" vertices " << f
1148 <<
" fc:" << f.centre(surf.points())
1149 <<
" has the same vertices as triangle " << nbrFaceI
1150 <<
" vertices " << nbrF
1151 <<
" fc:" << nbrF.centre(surf.points())
1164 bool Foam::isoSurface::validTri(
const triSurface& surf,
const label faceI)
1168 const labelledTri& f = surf[faceI];
1172 (f[0] < 0) || (f[0] >= surf.points().size())
1173 || (f[1] < 0) || (f[1] >= surf.points().size())
1174 || (f[2] < 0) || (f[2] >= surf.points().size())
1177 WarningIn(
"validTri(const triSurface&, const label)")
1178 <<
"triangle " << faceI <<
" vertices " << f
1179 <<
" uses point indices outside point range 0.."
1180 << surf.points().size()-1 <<
endl;
1185 if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
1187 WarningIn(
"validTri(const triSurface&, const label)")
1188 <<
"triangle " << faceI
1189 <<
" uses non-unique vertices " << f
1196 const labelList& fFaces = surf.faceFaces()[faceI];
1202 label nbrFaceI = fFaces[i];
1204 if (nbrFaceI <= faceI)
1210 const labelledTri& nbrF = surf[nbrFaceI];
1214 ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
1215 && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
1216 && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
1219 WarningIn(
"validTri(const triSurface&, const label)")
1220 <<
"triangle " << faceI <<
" vertices " << f
1221 <<
" fc:" << f.centre(surf.points())
1222 <<
" has the same vertices as triangle " << nbrFaceI
1223 <<
" vertices " << nbrF
1224 <<
" fc:" << nbrF.centre(surf.points())
1234 void Foam::isoSurface::calcAddressing
1236 const triSurface& surf,
1237 List<FixedList<label, 3> >& faceEdges,
1240 Map<labelList>& edgeFacesRest
1249 const labelledTri& tri = surf[triI];
1250 edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
1251 edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
1252 edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
1268 Pout<<
"isoSurface : detected "
1269 << mergedCentres.size()
1270 <<
" geometric edges on " << surf.size() <<
" triangles." <<
endl;
1275 if (surf.size() == 1)
1277 faceEdges.setSize(1);
1278 faceEdges[0][0] = 0;
1279 faceEdges[0][1] = 1;
1280 faceEdges[0][2] = 2;
1281 edgeFace0.setSize(1);
1283 edgeFace1.setSize(1);
1285 edgeFacesRest.clear();
1292 faceEdges.setSize(surf.size());
1296 faceEdges[triI][0] = oldToMerged[edgeI++];
1297 faceEdges[triI][1] = oldToMerged[edgeI++];
1298 faceEdges[triI][2] = oldToMerged[edgeI++];
1303 edgeFace0.setSize(mergedCentres.size());
1305 edgeFace1.setSize(mergedCentres.size());
1307 edgeFacesRest.clear();
1311 EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
1313 forAll(oldToMerged, oldEdgeI)
1315 label triI = oldEdgeI / 3;
1316 label edgeI = oldToMerged[oldEdgeI];
1318 if (edgeFace0[edgeI] == -1)
1321 edgeFace0[edgeI] = triI;
1327 const labelledTri& prevTri = surf[edgeFace0[edgeI]];
1328 const labelledTri& tri = surf[triI];
1330 label fp = oldEdgeI % 3;
1332 edge
e(tri[fp], tri[tri.fcIndex(fp)]);
1334 label prevTriIndex = -1;
1338 if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) ==
e)
1345 if (prevTriIndex == -1)
1348 EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(
e);
1350 if (iter != extraEdgeFaces.end())
1353 label sz = eFaces.size();
1354 eFaces.setSize(sz+1);
1359 extraEdgeFaces.insert(
e,
labelList(1, triI));
1362 else if (edgeFace1[edgeI] == -1)
1364 edgeFace1[edgeI] = triI;
1376 if (iter != edgeFacesRest.end())
1379 label sz = eFaces.size();
1380 eFaces.setSize(sz+1);
1385 edgeFacesRest.insert(edgeI,
labelList(1, triI));
1392 edgeI = edgeFace0.size();
1394 edgeFace0.setSize(edgeI + extraEdgeFaces.size());
1395 edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
1405 label triI = eFaces[i];
1406 const labelledTri& tri = surf[triI];
1408 FixedList<label, 3>& fEdges = faceEdges[triI];
1411 edge
e(tri[fp], tri[tri.fcIndex(fp)]);
1413 if (
e == iter.key())
1424 edgeFace0[edgeI] = eFaces[0];
1426 if (eFaces.size() >= 2)
1428 edgeFace1[edgeI] = eFaces[1];
1430 if (eFaces.size() > 2)
1432 edgeFacesRest.insert
1435 SubList<label>(eFaces, eFaces.size()-2, 2)
1445 void Foam::isoSurface::walkOrientation
1447 const triSurface& surf,
1448 const List<FixedList<label, 3> >& faceEdges,
1451 const label seedTriI,
1456 DynamicList<label> changedFaces(surf.size());
1458 changedFaces.append(seedTriI);
1460 while (changedFaces.size())
1462 DynamicList<label> newChangedFaces(changedFaces.size());
1466 label triI = changedFaces[i];
1467 const labelledTri& tri = surf[triI];
1468 const FixedList<label, 3>& fEdges = faceEdges[triI];
1472 label edgeI = fEdges[fp];
1476 label p1 = tri[tri.fcIndex(fp)];
1480 edgeFace0[edgeI] != triI
1485 if (nbrI != -1 && flipState[nbrI] == -1)
1487 const labelledTri& nbrTri = surf[nbrI];
1499 <<
" fEdges:" << fEdges
1500 <<
" edgeI:" << edgeI
1501 <<
" edgeFace0:" << edgeFace0[edgeI]
1502 <<
" edgeFace1:" << edgeFace1[edgeI]
1504 <<
" nbrTri:" << nbrTri
1509 label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
1511 bool sameOrientation = (p1 == nbrP1);
1513 if (flipState[triI] == 0)
1515 flipState[nbrI] = (sameOrientation ? 0 : 1);
1519 flipState[nbrI] = (sameOrientation ? 1 : 0);
1521 newChangedFaces.append(nbrI);
1526 changedFaces.transfer(newChangedFaces);
1531 void Foam::isoSurface::orientSurface
1534 const List<FixedList<label, 3> >& faceEdges,
1537 const Map<labelList>& edgeFacesRest
1553 seedTriI < surf.size() && flipState[seedTriI] != -1;
1558 if (seedTriI == surf.size())
1565 flipState[seedTriI] = 0;
1582 if (flipState[triI] == 1)
1584 labelledTri tri(surf[triI]);
1586 surf[triI][0] = tri[0];
1587 surf[triI][1] = tri[2];
1588 surf[triI][2] = tri[1];
1590 else if (flipState[triI] == -1)
1594 "isoSurface::orientSurface(triSurface&, const label)"
1602 bool Foam::isoSurface::danglingTriangle
1604 const FixedList<label, 3>& fEdges,
1611 if (edgeFace1[fEdges[i]] == -1)
1617 if (nOpen == 1 || nOpen == 2 || nOpen == 3)
1629 Foam::label Foam::isoSurface::markDanglingTriangles
1631 const List<FixedList<label, 3> >& faceEdges,
1634 const Map<labelList>& edgeFacesRest,
1638 keepTriangles.setSize(faceEdges.size());
1639 keepTriangles =
true;
1641 label nDangling = 0;
1649 label edgeI = iter.key();
1650 const labelList& otherEdgeFaces = iter();
1653 if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
1655 keepTriangles[edgeFace0[edgeI]] =
false;
1658 if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
1660 keepTriangles[edgeFace1[edgeI]] =
false;
1663 forAll(otherEdgeFaces, i)
1665 label triI = otherEdgeFaces[i];
1666 if (danglingTriangle(faceEdges[triI], edgeFace1))
1668 keepTriangles[triI] =
false;
1679 const triSurface& s,
1687 createWithValues<boolList>
1696 newToOldPoints.setSize(s.points().size());
1697 oldToNewPoints.setSize(s.points().size());
1698 oldToNewPoints = -1;
1702 forAll(include, oldFacei)
1704 if (include[oldFacei])
1707 const labelledTri& tri = s[oldFacei];
1711 label oldPointI = tri[fp];
1713 if (oldToNewPoints[oldPointI] == -1)
1715 oldToNewPoints[oldPointI] = pointI;
1716 newToOldPoints[pointI++] = oldPointI;
1721 newToOldPoints.setSize(pointI);
1726 forAll(newToOldPoints, i)
1728 newPoints[i] = s.points()[newToOldPoints[i]];
1731 List<labelledTri> newTriangles(newToOldFaces.size());
1736 const labelledTri& tri = s[newToOldFaces[i]];
1738 newTriangles[i][0] = oldToNewPoints[tri[0]];
1739 newTriangles[i][1] = oldToNewPoints[tri[1]];
1740 newTriangles[i][2] = oldToNewPoints[tri[2]];
1741 newTriangles[i].region() = tri.region();
1745 return triSurface(newTriangles, s.patches(), newPoints,
true);
1756 const bool regularise,
1757 const scalar mergeTol
1760 mesh_(cVals.
mesh()),
1763 regularise_(regularise),
1764 mergeDistance_(mergeTol*mesh_.bounds().mag())
1768 Pout<<
"isoSurface:" <<
nl
1769 <<
" isoField : " << cVals.
name() <<
nl
1770 <<
" cell min/max : "
1773 <<
" point min/max : "
1774 <<
min(pVals_) <<
" / "
1775 <<
max(pVals_) <<
nl
1776 <<
" isoValue : " << iso <<
nl
1777 <<
" regularise : " << regularise_ <<
nl
1778 <<
" mergeTol : " << mergeTol <<
nl
1790 cValsPtr_.reset(adaptPatchFields(cVals).ptr());
1803 mesh_.pointsInstance(),
1812 mesh_.cellCentres(),
1820 if (isA<coupledPolyPatch>(pp) && !collocatedPatch(pp))
1824 meshC.boundaryField()[patchI]
1829 collocatedFaces(refCast<const coupledPolyPatch>(pp))
1834 if (!isCollocated[i])
1836 pfld[i] = mesh_.faceCentres()[pp.
start()+i];
1840 else if (isA<emptyPolyPatch>(pp))
1844 bType& bfld =
const_cast<bType&
>(meshC.boundaryField());
1847 bfld.
set(patchI, NULL);
1855 mesh_.boundary()[patchI],
1861 bfld[patchI] = pp.
patchSlice(mesh_.faceCentres());
1868 labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1874 label faceI = pp.
start();
1878 boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
1886 calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
1908 snappedCc.
setSize(mesh_.nCells());
1916 Pout<<
"isoSurface : shifted " << snappedPoints.
size()
1917 <<
" cell centres to intersection." <<
endl;
1920 label nCellSnaps = snappedPoints.
size();
1935 if (patches[patchI].coupled())
1937 if (!collocatedPatch(patches[patchI]))
1940 refCast<const coupledPolyPatch>
1949 if (!isCollocated[i])
1951 const face& f = mesh_.faces()[cpp.
start()+i];
1955 isBoundaryPoint.set(f[fp], 1);
1967 const face& f = mesh_.faces()[pp.
start()+i];
1971 isBoundaryPoint.set(f[fp], 1);
1991 snappedPoint.
setSize(mesh_.nPoints());
1997 Pout<<
"isoSurface : shifted " << snappedPoints.
size()-nCellSnaps
1998 <<
" vertices to intersection." <<
endl;
2024 Pout<<
"isoSurface : generated " << triMeshCells.
size()
2025 <<
" unmerged triangles from " << triPoints.size()
2026 <<
" unmerged points." <<
endl;
2045 Pout<<
"isoSurface : generated " << triMap.
size()
2046 <<
" merged triangles." <<
endl;
2049 meshCells_.setSize(triMap.
size());
2052 meshCells_[i] = triMeshCells[triMap[i]];
2057 Pout<<
"isoSurface : checking " << size()
2058 <<
" triangles for validity." <<
endl;
2063 validTri(*
this, triI);
2078 calcAddressing(*
this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
2082 label nDangling = markDanglingTriangles
2093 Pout<<
"isoSurface : detected " << nDangling
2094 <<
" dangling triangles." <<
endl;
2117 meshCells_ =
labelField(meshCells_, subsetTriMap);
2121 orientSurface(*
this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
2128 Pout<<
"Dumping surface to " << stlFile <<
endl;
2129 triSurface::write(stlFile);