59 x.setSize(sz+
y.size());
71 bool Foam::snappySnapDriver::isFeaturePoint
73 const scalar featureCos,
75 const bitSet& isFeatureEdge,
81 const labelList& pEdges = pp.pointEdges()[pointi];
87 if (isFeatureEdge[pEdges[i]])
91 for (label j = i+1; j < pEdges.size(); j++)
93 if (isFeatureEdge[pEdges[j]])
95 const edge& ei = edges[pEdges[i]];
96 const edge& ej = edges[pEdges[j]];
103 scalar viMag =
mag(vi);
106 scalar vjMag =
mag(vj);
112 && ((vi/viMag & vj/vjMag) < featureCos)
132 void Foam::snappySnapDriver::smoothAndConstrain
134 const bitSet& isPatchMasterEdge,
137 const List<pointConstraint>& constraints,
141 const fvMesh&
mesh = meshRefiner_.mesh();
143 for (label avgIter = 0; avgIter < 20; avgIter++)
164 forAll(pointEdges, pointi)
166 const labelList& pEdges = pointEdges[pointi];
168 label nConstraints = constraints[pointi].first();
170 if (nConstraints <= 1)
174 label edgei = pEdges[i];
176 if (isPatchMasterEdge[edgei])
178 label nbrPointi = edges[edgei].otherVertex(pointi);
179 if (constraints[nbrPointi].first() >= nConstraints)
181 dispSum[pointi] += disp[nbrPointi];
209 forAll(constraints, pointi)
211 if (dispCount[pointi] > 0)
216 *(disp[pointi] + dispSum[pointi]/dispCount[pointi]);
223 void Foam::snappySnapDriver::calcNearestFace
234 const fvMesh&
mesh = meshRefiner_.mesh();
235 const refinementSurfaces& surfaces = meshRefiner_.surfaces();
238 faceDisp.setSize(pp.size());
240 faceSurfaceNormal.setSize(pp.size());
241 faceSurfaceNormal =
Zero;
242 faceSurfaceGlobalRegion.setSize(pp.size());
243 faceSurfaceGlobalRegion = -1;
260 const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
264 label zoneSurfi = zonedSurfaces[i];
266 const wordList& faceZoneNames = surfZones[zoneSurfi].faceZoneNames();
269 DynamicList<label> ppFaces;
270 DynamicList<label> meshFaces;
271 forAll(faceZoneNames, fzi)
273 const word& faceZoneName = faceZoneNames[fzi];
278 <<
"Problem. Cannot find zone " << faceZoneName
282 const bitSet isZonedFace(
mesh.
nFaces(), fZone);
284 ppFaces.reserve(ppFaces.capacity()+fZone.size());
285 meshFaces.reserve(meshFaces.capacity()+fZone.size());
287 forAll(pp.addressing(), i)
289 if (isZonedFace[pp.addressing()[i]])
291 snapSurf[i] = zoneSurfi;
293 meshFaces.append(pp.addressing()[i]);
306 IndirectList<face>(
mesh.
faces(), meshFaces),
311 List<pointIndexHit> hitInfo;
315 surfaces.findNearestRegion
328 if (hitInfo[hiti].hit())
330 label facei = ppFaces[hiti];
331 faceDisp[facei] = hitInfo[hiti].hitPoint() - fc[hiti];
332 faceSurfaceNormal[facei] = hitNormal[hiti];
333 faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
341 static label nWarn = 0;
346 <<
"Did not find surface near face centre " << fc[hiti]
352 <<
"Reached warning limit " << nWarn
353 <<
". Suppressing further warnings." <<
endl;
365 DynamicList<label> ppFaces(pp.size());
366 DynamicList<label> meshFaces(pp.size());
367 forAll(pp.addressing(), i)
369 if (snapSurf[i] == -1)
372 meshFaces.append(pp.addressing()[i]);
382 IndirectList<face>(
mesh.
faces(), meshFaces),
387 List<pointIndexHit> hitInfo;
391 surfaces.findNearestRegion
404 if (hitInfo[hiti].hit())
406 label facei = ppFaces[hiti];
407 faceDisp[facei] = hitInfo[hiti].hitPoint() - fc[hiti];
408 faceSurfaceNormal[facei] = hitNormal[hiti];
409 faceSurfaceGlobalRegion[facei] = surfaces.globalRegion
417 static label nWarn = 0;
422 <<
"Did not find surface near face centre " << fc[hiti]
429 <<
"Reached warning limit " << nWarn
430 <<
". Suppressing further warnings." <<
endl;
477 void Foam::snappySnapDriver::calcNearestFacePointProperties
484 const labelList& faceSurfaceGlobalRegion,
486 List<List<point>>& pointFaceSurfNormals,
487 List<List<point>>& pointFaceDisp,
488 List<List<point>>& pointFaceCentres,
489 List<labelList>& pointFacePatchID
492 const fvMesh&
mesh = meshRefiner_.mesh();
501 pointFaceSurfNormals.setSize(pp.nPoints());
502 pointFaceDisp.setSize(pp.nPoints());
503 pointFaceCentres.setSize(pp.nPoints());
504 pointFacePatchID.setSize(pp.nPoints());
507 forAll(pp.pointFaces(), pointi)
516 if (isMasterFace[facei] && faceSurfaceGlobalRegion[facei] != -1)
523 List<point>& pNormals = pointFaceSurfNormals[pointi];
524 pNormals.setSize(nFaces);
525 List<point>& pDisp = pointFaceDisp[pointi];
526 pDisp.setSize(nFaces);
527 List<point>& pFc = pointFaceCentres[pointi];
529 labelList& pFid = pointFacePatchID[pointi];
536 label globalRegioni = faceSurfaceGlobalRegion[facei];
538 if (isMasterFace[facei] && globalRegioni != -1)
540 pNormals[nFaces] = faceSurfaceNormal[facei];
541 pDisp[nFaces] = faceDisp[facei];
542 pFc[nFaces] = pp.faceCentres()[facei];
543 pFid[nFaces] = globalToMasterPatch_[globalRegioni];
564 const polyPatch& pp = pbm[patchi];
566 if (pp.coupled() || isA<emptyPolyPatch>(pp))
570 label meshFacei = pp.start()+i;
577 forAll(pp.addressing(), i)
579 label meshFacei = pp.addressing()[i];
596 bitSet isBoundaryPoint(pp.nPoints());
603 const edge&
e = edges[edgei];
604 const labelList& eFaces = edgeFaces[edgei];
606 if (eFaces.size() == 1)
609 isBoundaryPoint.
set(
e[0]);
610 isBoundaryPoint.set(
e[1]);
612 else if (eFaces.size() == 2 && bafflePair[eFaces[0]] == eFaces[1])
615 isBoundaryPoint.set(
e[0]);
616 isBoundaryPoint.set(
e[1]);
623 forAll(pp.meshPoints(), pointi)
625 meshToPatchPoint[pp.meshPoints()[pointi]] = pointi;
630 label patchi =
patchID[bFacei];
639 label pointi = meshToPatchPoint[
f[fp]];
641 if (pointi != -1 && isBoundaryPoint.test(pointi))
643 List<point>& pNormals = pointFaceSurfNormals[pointi];
644 List<point>& pDisp = pointFaceDisp[pointi];
645 List<point>& pFc = pointFaceCentres[pointi];
646 labelList& pFid = pointFacePatchID[pointi];
651 pNormals.append(fn/
mag(fn));
665 pointFaceSurfNormals,
666 listPlusEqOp<point>(),
675 listPlusEqOp<point>(),
684 listPlusEqOp<point>(),
686 mapDistribute::transformPosition()
693 listPlusEqOp<label>(),
701 forAll(pointFaceDisp, pointi)
703 List<point>& pNormals = pointFaceSurfNormals[pointi];
704 List<point>& pDisp = pointFaceDisp[pointi];
705 List<point>& pFc = pointFaceCentres[pointi];
706 labelList& pFid = pointFacePatchID[pointi];
710 pNormals = List<point>(pNormals, visitOrder);
711 pDisp = List<point>(pDisp, visitOrder);
712 pFc = List<point>(pFc, visitOrder);
721 void Foam::snappySnapDriver::correctAttraction
723 const DynamicList<point>& surfacePoints,
724 const DynamicList<label>& surfaceCounts,
733 scalar tang = ((pt-edgePt)&edgeNormal);
737 if (order[0] < order[1])
741 vector attractD = surfacePoints[order[0]]-edgePt;
743 scalar tang2 = (attractD&edgeNormal);
745 attractD -= tang2*edgeNormal;
747 scalar magAttractD =
mag(attractD);
748 scalar fraction = magAttractD/(magAttractD+
mag(edgeOffset));
752 + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
753 edgeOffset = linePt-pt;
762 const List<point>& faceCentres
768 label patch0 = patchIDs[0];
770 for (label i = 1; i < patchIDs.size(); i++)
772 if (patchIDs[i] != patch0)
782 Foam::label Foam::snappySnapDriver::findNormal
784 const scalar featureCos,
786 const DynamicList<vector>& surfaceNormals
793 scalar cosAngle = (
n&surfaceNormals[j]);
797 (cosAngle >= featureCos)
798 || (cosAngle < (-1+0.001))
819 const DynamicList<vector>& surfaceNormals,
823 if (patchIDs.empty())
829 label patch0 = patchIDs[0];
831 for (label i = 1; i < patchIDs.size(); i++)
833 if (patchIDs[i] != patch0)
847 if (surfaceNormals.size() == 1)
855 labelList normalToPatch(surfaceNormals.size(), -1);
856 forAll(faceToNormalBin, i)
858 if (faceToNormalBin[i] != -1)
860 label&
patch = normalToPatch[faceToNormalBin[i]];
866 else if (
patch == -2)
870 else if (
patch != patchIDs[i])
878 forAll(normalToPatch, normali)
880 if (normalToPatch[normali] == -2)
895 void Foam::snappySnapDriver::writeStats
898 const bitSet& isPatchMasterPoint,
899 const List<pointConstraint>& patchConstraints
902 label nMasterPoints = 0;
907 forAll(patchConstraints, pointi)
909 if (isPatchMasterPoint[pointi])
913 if (patchConstraints[pointi].first() == 1)
917 else if (patchConstraints[pointi].first() == 2)
921 else if (patchConstraints[pointi].first() == 3)
928 reduce(nMasterPoints, sumOp<label>());
929 reduce(nPlanar, sumOp<label>());
930 reduce(nEdge, sumOp<label>());
931 reduce(nPoint, sumOp<label>());
932 Info<<
"total master points :" << nMasterPoints
933 <<
" of which attracted to :" <<
nl
934 <<
" feature point : " << nPoint <<
nl
935 <<
" feature edge : " << nEdge <<
nl
936 <<
" nearest surface : " << nPlanar <<
nl
937 <<
" rest : " << nMasterPoints-nPoint-nEdge-nPlanar
943 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
946 const scalar featureCos,
953 const List<List<point>>& pointFaceSurfNormals,
954 const List<List<point>>& pointFaceDisp,
955 const List<List<point>>& pointFaceCentres,
958 DynamicList<point>& surfacePoints,
959 DynamicList<vector>& surfaceNormals,
963 pointConstraint& patchConstraint
966 patchAttraction =
Zero;
967 patchConstraint = pointConstraint();
969 const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
970 const List<point>& pfDisp = pointFaceDisp[pointi];
971 const List<point>& pfCentres = pointFaceCentres[pointi];
981 surfacePoints.clear();
982 surfaceNormals.clear();
985 faceToNormalBin.setSize(pfDisp.size());
986 faceToNormalBin = -1;
990 const point& fc = pfCentres[i];
991 const vector& fSNormal = pfSurfNormals[i];
992 const vector& fDisp = pfDisp[i];
995 if (
magSqr(fDisp) <
sqr(snapDist[pointi]) &&
mag(fSNormal) > VSMALL)
997 const point pt = fc + fDisp;
1000 faceToNormalBin[i] = findNormal
1007 if (faceToNormalBin[i] != -1)
1015 if (surfacePoints.size() <= 1)
1017 surfacePoints.append(pt);
1018 faceToNormalBin[i] = surfaceNormals.size();
1019 surfaceNormals.append(fSNormal);
1021 else if (surfacePoints.size() == 2)
1023 plane pl0(surfacePoints[0], surfaceNormals[0]);
1024 plane pl1(surfacePoints[1], surfaceNormals[1]);
1025 plane::ray r(pl0.planeIntersect(pl1));
1026 vector featureNormal = r.dir() /
mag(r.dir());
1028 if (
mag(fSNormal&featureNormal) >= 0.001)
1031 surfacePoints.append(pt);
1032 faceToNormalBin[i] = surfaceNormals.size();
1033 surfaceNormals.append(fSNormal);
1036 else if (surfacePoints.size() == 3)
1040 plane pl0(surfacePoints[0], surfaceNormals[0]);
1041 plane pl1(surfacePoints[1], surfaceNormals[1]);
1042 plane pl2(surfacePoints[2], surfaceNormals[2]);
1043 point p012(pl0.planePlaneIntersect(pl1, pl2));
1045 plane::ray r(pl0.planeIntersect(pl1));
1046 vector featureNormal = r.dir() /
mag(r.dir());
1047 if (
mag(fSNormal&featureNormal) >= 0.001)
1049 plane pl3(pt, fSNormal);
1050 point p013(pl0.planePlaneIntersect(pl1, pl3));
1052 if (
mag(p012-p013) > snapDist[pointi])
1055 surfacePoints.append(pt);
1056 faceToNormalBin[i] = surfaceNormals.
size();
1057 surfaceNormals.append(fSNormal);
1066 const point& pt = pp.localPoints()[pointi];
1069 if (surfaceNormals.size() == 1)
1073 ((surfacePoints[0]-pt) & surfaceNormals[0])
1082 patchAttraction = d;
1085 patchConstraint.applyConstraint(surfaceNormals[0]);
1087 else if (surfaceNormals.size() == 2)
1089 plane pl0(surfacePoints[0], surfaceNormals[0]);
1090 plane pl1(surfacePoints[1], surfaceNormals[1]);
1091 plane::ray r(pl0.planeIntersect(pl1));
1095 vector d = r.refPoint()-pt;
1104 patchAttraction = d;
1107 patchConstraint.applyConstraint(surfaceNormals[0]);
1108 patchConstraint.applyConstraint(surfaceNormals[1]);
1110 else if (surfaceNormals.size() == 3)
1113 plane pl0(surfacePoints[0], surfaceNormals[0]);
1114 plane pl1(surfacePoints[1], surfaceNormals[1]);
1115 plane pl2(surfacePoints[2], surfaceNormals[2]);
1116 point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
1117 vector d = cornerPt - pt;
1125 patchAttraction = d;
1128 patchConstraint.applyConstraint(surfaceNormals[0]);
1129 patchConstraint.applyConstraint(surfaceNormals[1]);
1130 patchConstraint.applyConstraint(surfaceNormals[2]);
1136 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
1139 const scalar featureCos,
1145 const List<List<point>>& pointFaceSurfNormals,
1146 const List<List<point>>& pointFaceDisp,
1147 const List<List<point>>& pointFaceCentres,
1151 List<pointConstraint>& patchConstraints
1154 autoPtr<OBJstream> feStr;
1155 autoPtr<OBJstream> fpStr;
1162 meshRefiner_.mesh().time().path()
1163 /
"implicitFeatureEdge_" +
name(iter) +
".obj"
1166 Info<<
"Dumping implicit feature-edge direction to "
1167 << feStr().name() <<
endl;
1173 meshRefiner_.mesh().time().path()
1174 /
"implicitFeaturePoint_" +
name(iter) +
".obj"
1177 Info<<
"Dumping implicit feature-point direction to "
1178 << fpStr().name() <<
endl;
1182 DynamicList<point> surfacePoints(4);
1183 DynamicList<vector> surfaceNormals(4);
1186 forAll(pp.localPoints(), pointi)
1189 pointConstraint constraint;
1191 featureAttractionUsingReconstruction
1202 pointFaceSurfNormals,
1217 (constraint.first() > patchConstraints[pointi].first())
1219 (constraint.first() == patchConstraints[pointi].first())
1220 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
1224 patchAttraction[pointi] = attraction;
1225 patchConstraints[pointi] = constraint;
1227 const point& pt = pp.localPoints()[pointi];
1229 if (patchConstraints[pointi].first() == 2 && feStr.valid())
1231 feStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1233 else if (patchConstraints[pointi].first() == 3 && fpStr.valid())
1235 fpStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1242 void Foam::snappySnapDriver::stringFeatureEdges
1245 const scalar featureCos,
1251 const List<pointConstraint>& rawPatchConstraints,
1254 List<pointConstraint>& patchConstraints
1284 forAll(pointEdges, pointi)
1286 if (patchConstraints[pointi].first() == 2)
1288 const point& pt = pp.localPoints()[pointi];
1289 const labelList& pEdges = pointEdges[pointi];
1290 const vector& featVec = patchConstraints[pointi].second();
1294 bool hasPos =
false;
1295 bool hasNeg =
false;
1299 const edge&
e = pp.edges()[pEdges[pEdgei]];
1300 label nbrPointi =
e.otherVertex(pointi);
1302 if (patchConstraints[nbrPointi].first() > 1)
1304 const point& nbrPt = pp.localPoints()[nbrPointi];
1305 const point featPt =
1306 nbrPt + patchAttraction[nbrPointi];
1307 const scalar cosAngle = (featVec & (featPt-pt));
1320 if (!hasPos || !hasNeg)
1326 label bestPosPointi = -1;
1327 scalar minPosDistSqr = GREAT;
1328 label bestNegPointi = -1;
1329 scalar minNegDistSqr = GREAT;
1333 const edge&
e = pp.edges()[pEdges[pEdgei]];
1334 label nbrPointi =
e.otherVertex(pointi);
1338 patchConstraints[nbrPointi].first() <= 1
1339 && rawPatchConstraints[nbrPointi].first() > 1
1342 const vector& nbrFeatVec =
1343 rawPatchConstraints[pointi].second();
1345 if (
mag(featVec&nbrFeatVec) > featureCos)
1352 rawPatchAttraction[nbrPointi]
1355 const point featPt =
1356 pp.localPoints()[nbrPointi]
1357 + rawPatchAttraction[nbrPointi];
1358 const scalar cosAngle =
1359 (featVec & (featPt-pt));
1363 if (!hasPos && d2 < minPosDistSqr)
1366 bestPosPointi = nbrPointi;
1371 if (!hasNeg && d2 < minNegDistSqr)
1374 bestNegPointi = nbrPointi;
1381 if (bestPosPointi != -1)
1391 patchAttraction[bestPosPointi] =
1392 0.5*rawPatchAttraction[bestPosPointi];
1393 patchConstraints[bestPosPointi] =
1394 rawPatchConstraints[bestPosPointi];
1398 if (bestNegPointi != -1)
1408 patchAttraction[bestNegPointi] =
1409 0.5*rawPatchAttraction[bestNegPointi];
1410 patchConstraints[bestNegPointi] =
1411 rawPatchConstraints[bestNegPointi];
1420 reduce(nChanged, sumOp<label>());
1421 Info<<
"Stringing feature edges : changed " << nChanged <<
" points"
1431 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1434 const scalar featureCos,
1439 const List<List<point>>& pointFaceCentres,
1443 const List<pointConstraint>& rawPatchConstraints,
1446 List<pointConstraint>& patchConstraints
1449 autoPtr<OBJstream> multiPatchStr;
1456 meshRefiner_.mesh().time().path()
1457 /
"multiPatch_" +
name(iter) +
".obj"
1460 Info<<
"Dumping removed constraints due to same-face"
1461 <<
" multi-patch points to "
1462 << multiPatchStr().name() <<
endl;
1467 bitSet isMultiPatchPoint(pp.size());
1469 forAll(pointFacePatchID, pointi)
1473 pp.localPoints()[pointi],
1474 pointFacePatchID[pointi],
1475 pointFaceCentres[pointi]
1477 isMultiPatchPoint.set(pointi, multiPatchPt.hit());
1481 forAll(isMultiPatchPoint, pointi)
1483 if (isMultiPatchPoint.test(pointi))
1487 patchConstraints[pointi].first() <= 1
1488 && rawPatchConstraints[pointi].first() > 1
1491 patchAttraction[pointi] = rawPatchAttraction[pointi];
1492 patchConstraints[pointi] = rawPatchConstraints[pointi];
1511 forAll(pp.localFaces(), facei)
1513 const face&
f = pp.localFaces()[facei];
1515 label nMultiPatchPoints = 0;
1518 label pointi =
f[fp];
1521 isMultiPatchPoint.test(pointi)
1522 && patchConstraints[pointi].first() > 1
1525 ++nMultiPatchPoints;
1529 if (nMultiPatchPoints > 0)
1533 label pointi =
f[fp];
1536 !isMultiPatchPoint.test(pointi)
1537 && patchConstraints[pointi].first() > 1
1543 patchAttraction[pointi] =
Zero;
1544 patchConstraints[pointi] = pointConstraint();
1547 if (multiPatchStr.valid())
1549 multiPatchStr().write(pp.localPoints()[pointi]);
1556 reduce(nChanged, sumOp<label>());
1557 Info<<
"Removing constraints near multi-patch points : changed "
1558 << nChanged <<
" points" <<
endl;
1566 const List<pointConstraint>& patchConstraints,
1570 const face&
f = pp.localFaces()[facei];
1578 for (label startFp = 0; startFp <
f.size()-2; startFp++)
1580 label minFp =
f.rcIndex(startFp);
1584 label endFp =
f.fcIndex(
f.fcIndex(startFp));
1585 endFp <
f.size() && endFp != minFp;
1591 patchConstraints[
f[startFp]].first() >= 2
1592 && patchConstraints[
f[endFp]].first() >= 2
1595 attractIndices =
labelPair(startFp, endFp);
1601 return attractIndices;
1605 bool Foam::snappySnapDriver::isSplitAlignedWithFeature
1607 const scalar featureCos,
1609 const pointConstraint& pc0,
1611 const pointConstraint& pc1
1615 scalar magD =
mag(d);
1628 if (pc0.first() == 2 &&
mag(d & pc0.second()) > featureCos)
1632 else if (pc1.first() == 2 &&
mag(d & pc1.second()) > featureCos)
1645 bool Foam::snappySnapDriver::isConcave
1651 const scalar concaveCos
1655 scalar magN0 =
mag(n0);
1664 scalar d = (
c1-c0)&n0;
1675 scalar magN1 =
mag(n1);
1683 if ((n0&n1) < concaveCos)
1697 const scalar featureCos,
1698 const scalar concaveCos,
1699 const scalar minAreaRatio,
1702 const List<pointConstraint>& patchConstraints,
1708 DynamicField<point>& points1
1711 const face& localF = pp.localFaces()[facei];
1715 if (localF.size() >= 4)
1717 const pointField& localPts = pp.localPoints();
1755 for (label startFp = 0; startFp < localF.size()-2; startFp++)
1757 label minFp = localF.rcIndex(startFp);
1761 label endFp = localF.fcIndex(localF.fcIndex(startFp));
1762 endFp < localF.size() && endFp != minFp;
1766 label startPti = localF[startFp];
1767 label endPti = localF[endFp];
1769 const pointConstraint& startPc = patchConstraints[startPti];
1770 const pointConstraint& endPc = patchConstraints[endPti];
1772 if (startPc.first() >= 2 && endPc.first() >= 2)
1774 if (startPc.first() == 2 || endPc.first() == 2)
1779 point start = localPts[startPti]+patchAttr[startPti];
1780 point end = localPts[endPti]+patchAttr[endPti];
1784 !isSplitAlignedWithFeature
1810 f0.setSize(endFp-startFp+1);
1812 for (label fp = startFp; fp <= endFp; fp++)
1814 f0[i0++] = localF[fp];
1818 const face compact0(
identity(f0.size()));
1820 points0.append(localPts[f0[0]] + patchAttr[f0[0]]);
1821 for (label fp=1; fp < f0.size()-1; fp++)
1828 localPts[f0.last()] + patchAttr[f0.last()]
1834 f1.setSize(localF.size()+2-f0.size());
1840 fp = localF.fcIndex(fp)
1843 f1[i1++] = localF[fp];
1845 f1[i1++] = localF[startFp];
1849 const face compact1(
identity(f1.size()));
1851 points1.append(localPts[f1[0]] + patchAttr[f1[0]]);
1852 for (label fp=1; fp < f1.size()-1; fp++)
1855 points1.append(localPts[
pi] + nearestAttr[
pi]);
1859 localPts[f1.last()] + patchAttr[f1.last()]
1873 compact1.centre(points1),
1874 compact1.areaNormal(points1),
1894 const scalar area0 = f0.mag(localPts);
1895 const scalar area1 = f1.mag(localPts);
1899 area0/area1 >= minAreaRatio
1900 && area1/area0 >= minAreaRatio
1903 attractIndices =
labelPair(startFp, endFp);
1910 return attractIndices;
1914 void Foam::snappySnapDriver::splitDiagonals
1916 const scalar featureCos,
1917 const scalar concaveCos,
1918 const scalar minAreaRatio,
1925 List<pointConstraint>& patchConstraints,
1926 DynamicList<label>& splitFaces,
1927 DynamicList<labelPair>& splits
1930 const labelList& bFaces = pp.addressing();
1933 splitFaces.setCapacity(bFaces.size());
1935 splits.setCapacity(bFaces.size());
1939 DynamicField<point> facePoints0;
1940 DynamicField<point> facePoints1;
1946 findDiagonalAttraction
1967 splitFaces.append(bFaces[facei]);
1968 splits.append(
split);
1970 const face&
f = pp.localFaces()[facei];
1980 && patchConstraints[
f[fp]].first() >= 2
1983 patchConstraints[
f[fp]] = pointConstraint
1985 Tuple2<label, vector>
1988 nearestNormal[
f[fp]]
1991 patchAttraction[
f[fp]] = nearestAttraction[
f[fp]];
1999 void Foam::snappySnapDriver::avoidDiagonalAttraction
2002 const scalar featureCos,
2007 List<pointConstraint>& patchConstraints
2010 forAll(pp.localFaces(), facei)
2012 const face&
f = pp.localFaces()[facei];
2022 if (
diag[0] != -1 &&
diag[1] != -1)
2026 const label i0 =
f[
diag[0]];
2028 pp.localPoints()[i0]+patchAttraction[i0];
2029 const label i1 =
f[
diag[1]];
2031 pp.localPoints()[i1]+patchAttraction[i1];
2032 const point mid = 0.5*(pt0+pt1);
2034 const scalar cosAngle =
mag
2036 patchConstraints[i0].second()
2037 & patchConstraints[i1].second()
2046 if (cosAngle > featureCos)
2052 scalar minDistSqr = GREAT;
2055 label pointi =
f[fp];
2056 if (patchConstraints[pointi].first() <= 1)
2058 const point& pt = pp.localPoints()[pointi];
2059 scalar distSqr =
magSqr(mid-pt);
2060 if (distSqr < minDistSqr)
2068 label minPointi =
f[minFp];
2069 patchAttraction[minPointi] =
2070 mid-pp.localPoints()[minPointi];
2071 patchConstraints[minPointi] = patchConstraints[
f[
diag[0]]];
2097 Foam::snappySnapDriver::findNearFeatureEdge
2099 const bool isRegionEdge,
2104 const point& estimatedPt,
2106 List<List<DynamicList<point>>>& edgeAttractors,
2107 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2109 List<pointConstraint>& patchConstraints
2112 const refinementFeatures& features = meshRefiner_.features();
2115 List<pointIndexHit> nearEdgeInfo;
2120 features.findNearestRegionEdge
2131 features.findNearestEdge
2142 label feati = nearEdgeFeat[0];
2148 edgeAttractors[feati][
nearInfo.index()].append
2152 pointConstraint
c(Tuple2<label, vector>(2, nearNormal[0]));
2153 edgeConstraints[feati][
nearInfo.index()].append(
c);
2156 patchAttraction[pointi] =
nearInfo.hitPoint()-pp.localPoints()[pointi];
2157 patchConstraints[pointi] =
c;
2159 return Tuple2<label, pointIndexHit>(feati,
nearInfo);
2164 Foam::snappySnapDriver::findNearFeaturePoint
2166 const bool isRegionPoint,
2171 const point& estimatedPt,
2174 List<labelList>& pointAttractor,
2175 List<List<pointConstraint>>& pointConstraints,
2177 List<List<DynamicList<point>>>& edgeAttractors,
2178 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2181 List<pointConstraint>& patchConstraints
2184 const refinementFeatures& features = meshRefiner_.features();
2190 features.findNearestPoint
2198 label feati = nearFeat[0];
2202 const point& pt = pp.localPoints()[pointi];
2204 label featPointi =
nearInfo[0].index();
2206 scalar distSqr =
magSqr(featPt-pt);
2209 label oldPointi = pointAttractor[feati][featPointi];
2211 if (oldPointi != -1)
2214 if (distSqr >=
magSqr(featPt-pp.localPoints()[oldPointi]))
2223 pointAttractor[feati][featPointi] = pointi;
2224 pointConstraints[feati][featPointi].first() = 3;
2225 pointConstraints[feati][featPointi].second() =
Zero;
2228 patchAttraction[pointi] = featPt-pt;
2229 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2232 patchAttraction[oldPointi] =
Zero;
2233 patchConstraints[oldPointi] = pointConstraint();
2242 pp.localPoints()[oldPointi],
2254 pointAttractor[feati][featPointi] = pointi;
2255 pointConstraints[feati][featPointi].first() = 3;
2256 pointConstraints[feati][featPointi].second() =
Zero;
2259 patchAttraction[pointi] = featPt-pt;
2260 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2264 return Tuple2<label, pointIndexHit>(feati,
nearInfo[0]);
2270 void Foam::snappySnapDriver::determineFeatures
2273 const scalar featureCos,
2274 const bool multiRegionFeatureSnap,
2280 const List<List<point>>& pointFaceSurfNormals,
2281 const List<List<point>>& pointFaceDisp,
2282 const List<List<point>>& pointFaceCentres,
2286 List<labelList>& pointAttractor,
2287 List<List<pointConstraint>>& pointConstraints,
2289 List<List<DynamicList<point>>>& edgeAttractors,
2290 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2293 List<pointConstraint>& patchConstraints
2296 autoPtr<OBJstream> featureEdgeStr;
2297 autoPtr<OBJstream> missedEdgeStr;
2298 autoPtr<OBJstream> featurePointStr;
2299 autoPtr<OBJstream> missedMP0Str;
2300 autoPtr<OBJstream> missedMP1Str;
2304 featureEdgeStr.reset
2308 meshRefiner_.mesh().time().path()
2309 /
"featureEdge_" +
name(iter) +
".obj"
2312 Info<<
"Dumping feature-edge sampling to "
2313 << featureEdgeStr().name() <<
endl;
2319 meshRefiner_.mesh().time().path()
2320 /
"missedFeatureEdge_" +
name(iter) +
".obj"
2323 Info<<
"Dumping feature-edges that are too far away to "
2324 << missedEdgeStr().name() <<
endl;
2326 featurePointStr.reset
2330 meshRefiner_.mesh().time().path()
2331 /
"featurePoint_" +
name(iter) +
".obj"
2334 Info<<
"Dumping feature-point sampling to "
2335 << featurePointStr().name() <<
endl;
2341 meshRefiner_.mesh().time().path()
2342 /
"missedFeatureEdgeFromMPEdge_" +
name(iter) +
".obj"
2345 Info<<
"Dumping region-edges that are too far away to "
2346 << missedMP0Str().name() <<
endl;
2352 meshRefiner_.mesh().time().path()
2353 /
"missedFeatureEdgeFromMPPoint_" +
name(iter) +
".obj"
2356 Info<<
"Dumping region-points that are too far away to "
2357 << missedMP1Str().name() <<
endl;
2361 DynamicList<point> surfacePoints(4);
2362 DynamicList<vector> surfaceNormals(4);
2365 forAll(pp.localPoints(), pointi)
2367 const point& pt = pp.localPoints()[pointi];
2379 pointConstraint constraint;
2381 featureAttractionUsingReconstruction
2392 pointFaceSurfNormals,
2437 (constraint.first() > patchConstraints[pointi].first())
2439 (constraint.first() == patchConstraints[pointi].first())
2440 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
2444 patchAttraction[pointi] = attraction;
2445 patchConstraints[pointi] = constraint;
2448 if (patchConstraints[pointi].first() == 1)
2451 if (multiRegionFeatureSnap)
2453 const point estimatedPt(pt + nearestDisp[pointi]);
2459 pointFacePatchID[pointi],
2465 if (multiPatchPt.hit())
2470 Tuple2<label, pointIndexHit>
nearInfo =
2477 multiPatchPt.hitPoint(),
2490 if (featureEdgeStr.valid())
2492 featureEdgeStr().write
2500 if (missedEdgeStr.valid())
2502 missedEdgeStr().write
2511 else if (patchConstraints[pointi].first() == 2)
2516 const point estimatedPt(pt + patchAttraction[pointi]);
2521 bool hasSnapped =
false;
2522 if (multiRegionFeatureSnap)
2529 pointFacePatchID[pointi],
2534 if (multiPatchPt.hit())
2536 if (multiPatchPt.index() == 0)
2558 missedMP0Str.valid()
2562 missedMP0Str().write
2627 missedMP1Str.valid()
2631 missedMP1Str().write
2669 patchConstraints[pointi].first() == 3
2670 && featurePointStr.valid()
2673 featurePointStr().write
2680 patchConstraints[pointi].first() == 2
2681 && featureEdgeStr.valid()
2684 featureEdgeStr().write
2692 if (missedEdgeStr.valid())
2694 missedEdgeStr().write
2701 else if (patchConstraints[pointi].first() == 3)
2704 const point estimatedPt(pt + patchAttraction[pointi]);
2708 if (multiRegionFeatureSnap)
2715 pointFacePatchID[pointi],
2720 if (multiPatchPt.hit())
2788 if (info.hit() && featurePointStr.valid())
2790 featurePointStr().write
2812 void Foam::snappySnapDriver::determineBaffleFeatures
2815 const bool baffleFeaturePoints,
2816 const scalar featureCos,
2822 List<labelList>& pointAttractor,
2823 List<List<pointConstraint>>& pointConstraints,
2825 List<List<DynamicList<point>>>& edgeAttractors,
2826 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2829 List<pointConstraint>& patchConstraints
2832 const fvMesh&
mesh = meshRefiner_.mesh();
2833 const refinementFeatures& features = meshRefiner_.features();
2836 List<List<point>> edgeFaceNormals(pp.nEdges());
2839 forAll(pp.edgeFaces(), edgei)
2841 const labelList& eFaces = pp.edgeFaces()[edgei];
2842 List<point>& eFc = edgeFaceNormals[edgei];
2846 label facei = eFaces[i];
2847 eFc[i] = pp.faceNormals()[facei];
2862 listPlusEqOp<point>(),
2874 autoPtr<OBJstream> baffleEdgeStr;
2881 meshRefiner_.mesh().time().path()
2882 /
"baffleEdge_" +
name(iter) +
".obj"
2885 Info<<
nl <<
"Dumping baffle-edges to "
2886 << baffleEdgeStr().name() <<
endl;
2891 bitSet isBaffleEdge(pp.nEdges());
2892 label nBaffleEdges = 0;
2897 labelList pointStatus(pp.nPoints(), -1);
2899 forAll(edgeFaceNormals, edgei)
2901 const List<point>& efn = edgeFaceNormals[edgei];
2903 if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2905 isBaffleEdge.set(edgei);
2907 const edge&
e = pp.edges()[edgei];
2908 pointStatus[
e[0]] = 0;
2909 pointStatus[
e[1]] = 0;
2911 if (baffleEdgeStr.valid())
2913 const point&
p0 = pp.localPoints()[
e[0]];
2914 const point& p1 = pp.localPoints()[
e[1]];
2920 reduce(nBaffleEdges, sumOp<label>());
2922 Info<<
"Detected " << nBaffleEdges
2923 <<
" baffle edges out of "
2925 <<
" edges." <<
endl;
2950 label nBafflePoints = 0;
2951 forAll(pointStatus, pointi)
2953 if (pointStatus[pointi] != -1)
2958 reduce(nBafflePoints, sumOp<label>());
2961 label nPointAttract = 0;
2962 label nEdgeAttract = 0;
2964 forAll(pointStatus, pointi)
2966 const point& pt = pp.localPoints()[pointi];
2968 if (pointStatus[pointi] == 0)
2972 Tuple2<label, pointIndexHit>
nearInfo = findNearFeatureEdge
2995 if (baffleFeaturePoints)
3025 else if (pointStatus[pointi] == 1)
3029 features.findNearestPoint
3037 label feati = nearFeat[0];
3043 label featPointi =
nearInfo[0].index();
3045 scalar distSqr =
magSqr(featPt-pt);
3048 label oldPointi = pointAttractor[feati][featPointi];
3055 <
magSqr(featPt-pp.localPoints()[oldPointi])
3059 pointAttractor[feati][featPointi] = pointi;
3060 pointConstraints[feati][featPointi].first() = 3;
3061 pointConstraints[feati][featPointi].second() =
Zero;
3064 patchAttraction[pointi] = featPt-pt;
3065 patchConstraints[pointi] =
3066 pointConstraints[feati][featPointi];
3068 if (oldPointi != -1)
3079 pp.localPoints()[oldPointi],
3104 Tuple2<label, pointIndexHit>
nearInfo = findNearFeatureEdge
3126 reduce(nPointAttract, sumOp<label>());
3127 reduce(nEdgeAttract, sumOp<label>());
3129 Info<<
"Baffle points : " << nBafflePoints
3130 <<
" of which attracted to :" <<
nl
3131 <<
" feature point : " << nPointAttract <<
nl
3132 <<
" feature edge : " << nEdgeAttract <<
nl
3133 <<
" rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3139 void Foam::snappySnapDriver::reverseAttractMeshPoints
3147 const List<labelList>& pointAttractor,
3148 const List<List<pointConstraint>>& pointConstraints,
3150 const List<List<DynamicList<point>>>& edgeAttractors,
3151 const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3154 const List<pointConstraint>& rawPatchConstraints,
3158 List<pointConstraint>& patchConstraints
3161 const refinementFeatures& features = meshRefiner_.features();
3169 treeBoundBox bb(pp.localPoints());
3182 DynamicList<label> attractPoints(pp.nPoints());
3184 const fvMesh&
mesh = meshRefiner_.mesh();
3186 boolList isFeatureEdgeOrPoint(pp.nPoints(),
false);
3188 forAll(rawPatchConstraints, pointi)
3190 if (rawPatchConstraints[pointi].first() >= 2)
3192 isFeatureEdgeOrPoint[pointi] =
true;
3198 <<
" mesh points out of "
3200 <<
" for reverse attraction." <<
endl;
3208 isFeatureEdgeOrPoint,
3213 for (label nGrow = 0; nGrow < 1; nGrow++)
3215 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3217 forAll(pp.localFaces(), facei)
3219 const face&
f = pp.localFaces()[facei];
3223 if (isFeatureEdgeOrPoint[
f[fp]])
3228 newIsFeatureEdgeOrPoint[
f[fp]] =
true;
3235 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3241 isFeatureEdgeOrPoint,
3249 forAll(isFeatureEdgeOrPoint, pointi)
3251 if (isFeatureEdgeOrPoint[pointi])
3253 attractPoints.append(pointi);
3259 <<
" mesh points out of "
3261 <<
" for reverse attraction." <<
endl;
3265 indexedOctree<treeDataPoint> ppTree
3267 treeDataPoint(pp.localPoints(), attractPoints),
3275 patchAttraction.setSize(pp.nPoints());
3276 patchAttraction =
Zero;
3277 patchConstraints.setSize(pp.nPoints());
3278 patchConstraints = pointConstraint();
3280 forAll(edgeAttractors, feati)
3282 const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3283 const List<DynamicList<pointConstraint>>& edgeConstr =
3284 edgeConstraints[feati];
3286 forAll(edgeAttr, featEdgei)
3288 const DynamicList<point>& attr = edgeAttr[featEdgei];
3292 const point& featPt = attr[i];
3302 ppTree.shapes().pointLabels()[
nearInfo.index()];
3303 const point attraction = featPt-pp.localPoints()[pointi];
3309 patchConstraints[pointi].first() <= 1
3310 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
3313 patchAttraction[pointi] = attraction;
3314 patchConstraints[pointi] = edgeConstr[featEdgei][i];
3319 static label nWarn = 0;
3324 <<
"Did not find pp point near " << featPt
3330 <<
"Reached warning limit " << nWarn
3331 <<
". Suppressing further warnings." <<
endl;
3350 forAll(pointAttractor, feati)
3352 const labelList& pointAttr = pointAttractor[feati];
3353 const List<pointConstraint>& pointConstr = pointConstraints[feati];
3355 forAll(pointAttr, featPointi)
3357 if (pointAttr[featPointi] != -1)
3359 const point& featPt = features[feati].points()
3374 ppTree.shapes().pointLabels()[
nearInfo.index()];
3376 const point& pt = pp.localPoints()[pointi];
3377 const point attraction = featPt-pt;
3382 if (patchConstraints[pointi].first() <= 1)
3384 patchAttraction[pointi] = attraction;
3385 patchConstraints[pointi] = pointConstr[featPointi];
3387 else if (patchConstraints[pointi].first() == 2)
3389 patchAttraction[pointi] = attraction;
3390 patchConstraints[pointi] = pointConstr[featPointi];
3392 else if (patchConstraints[pointi].first() == 3)
3398 <
magSqr(patchAttraction[pointi])
3401 patchAttraction[pointi] = attraction;
3402 patchConstraints[pointi] =
3403 pointConstr[featPointi];
3413 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3416 const bool multiRegionFeatureSnap,
3418 const bool detectBaffles,
3419 const bool baffleFeaturePoints,
3421 const bool releasePoints,
3422 const bool stringFeatures,
3423 const bool avoidDiagonal,
3425 const scalar featureCos,
3432 const List<List<point>>& pointFaceSurfNormals,
3433 const List<List<point>>& pointFaceDisp,
3434 const List<List<point>>& pointFaceCentres,
3438 List<pointConstraint>& patchConstraints
3441 const refinementFeatures& features = meshRefiner_.features();
3442 const fvMesh&
mesh = meshRefiner_.mesh();
3444 const bitSet isPatchMasterPoint
3459 List<List<DynamicList<point>>> edgeAttractors(features.size());
3460 List<List<DynamicList<pointConstraint>>> edgeConstraints
3466 label nFeatEdges = features[feati].edges().size();
3467 edgeAttractors[feati].setSize(nFeatEdges);
3468 edgeConstraints[feati].setSize(nFeatEdges);
3474 List<labelList> pointAttractor(features.size());
3475 List<List<pointConstraint>> pointConstraints(features.size());
3478 label nFeatPoints = features[feati].points().size();
3479 pointAttractor[feati].setSize(nFeatPoints, -1);
3480 pointConstraints[feati].setSize(nFeatPoints);
3485 List<pointConstraint> rawPatchConstraints(pp.nPoints());
3491 multiRegionFeatureSnap,
3497 pointFaceSurfNormals,
3516 Info<<
"Raw geometric feature analysis : ";
3517 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3533 determineBaffleFeatures
3536 baffleFeaturePoints,
3557 Info<<
"After baffle feature analysis : ";
3558 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3566 reverseAttractMeshPoints
3582 rawPatchConstraints,
3592 Info<<
"Reverse attract feature analysis : ";
3593 writeStats(pp, isPatchMasterPoint, patchConstraints);
3599 OBJstream featureEdgeStr
3601 meshRefiner_.mesh().time().path()
3602 /
"edgeAttractors_" +
name(iter) +
".obj"
3604 Info<<
"Dumping feature-edge attraction to "
3605 << featureEdgeStr.name() <<
endl;
3607 OBJstream featurePointStr
3609 meshRefiner_.mesh().time().path()
3610 /
"pointAttractors_" +
name(iter) +
".obj"
3612 Info<<
"Dumping feature-point attraction to "
3613 << featurePointStr.name() <<
endl;
3615 forAll(patchConstraints, pointi)
3617 const point& pt = pp.localPoints()[pointi];
3618 const vector& attr = patchAttraction[pointi];
3620 if (patchConstraints[pointi].first() == 2)
3624 else if (patchConstraints[pointi].first() == 3)
3638 releasePointsNextToMultiPatch
3650 rawPatchConstraints,
3673 rawPatchConstraints,
3686 avoidDiagonalAttraction
3701 meshRefiner_.mesh().time().path()
3702 /
"patchAttraction_" +
name(iter) +
".obj",
3704 pp.localPoints() + patchAttraction
3711 void Foam::snappySnapDriver::preventFaceSqueeze
3714 const scalar featureCos,
3721 List<pointConstraint>& patchConstraints
3724 autoPtr<OBJstream> strPtr;
3731 meshRefiner_.mesh().time().path()
3732 /
"faceSqueeze_" +
name(iter) +
".obj"
3735 Info<<
"Dumping faceSqueeze corrections to "
3736 << strPtr().name() <<
endl;
3741 forAll(pp.localFaces(), facei)
3743 const face&
f = pp.localFaces()[facei];
3745 if (
f.size() !=
points.size())
3748 singleF.setSize(
f.size());
3749 for (label i = 0; i <
f.size(); i++)
3754 label nConstraints = 0;
3757 label pointi =
f[fp];
3758 const point& pt = pp.localPoints()[pointi];
3760 if (patchConstraints[pointi].first() > 1)
3762 points[fp] = pt + patchAttraction[pointi];
3771 if (nConstraints ==
f.size())
3782 const point& pt = pp.localPoints()[
f[fp]];
3783 const point& nextPt = pp.localPoints()[
f.nextLabel(fp)];
3794 label pointi =
f.prevLabel(maxFp);
3798 const point& pt = pp.localPoints()[pointi];
3805 patchAttraction[pointi] = nearestAttraction[pointi];
3818 scalar oldArea =
f.mag(pp.localPoints());
3819 scalar newArea = singleF.mag(
points);
3820 if (newArea < 0.1*oldArea)
3827 scalar
s =
magSqr(patchAttraction[
f[fp]]);
3836 label pointi =
f[maxFp];
3838 patchAttraction[pointi] *= 0.5;
3849 const snapParameters& snapParams,
3850 const bool alignMeshEdges,
3852 const scalar featureCos,
3853 const scalar featureAttract,
3857 motionSmoother& meshMover,
3859 List<pointConstraint>& patchConstraints,
3861 DynamicList<label>& splitFaces,
3862 DynamicList<labelPair>& splits
3872 const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3873 const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3874 const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3876 Info<<
"Overriding displacement on features :" <<
nl
3877 <<
" implicit features : " << implicitFeatureAttraction <<
nl
3878 <<
" explicit features : " << explicitFeatureAttraction <<
nl
3879 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl
3883 const pointField& localPoints = pp.localPoints();
3884 const fvMesh&
mesh = meshRefiner_.mesh();
3888 const bitSet isPatchMasterPoint
3901 List<List<point>> pointFaceSurfNormals;
3902 List<List<point>> pointFaceDisp;
3903 List<List<point>> pointFaceCentres;
3904 List<labelList> pointFacePatchID;
3910 forAll(pp.localFaces(), facei)
3912 const face&
f = pp.localFaces()[facei];
3915 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[
f[fp]]);
3927 labelList faceSurfaceGlobalRegion(pp.size(), -1);
3937 faceSurfaceGlobalRegion
3947 calcNearestFacePointProperties
3954 faceSurfaceGlobalRegion,
3956 pointFaceSurfNormals,
3975 patchAttraction.setSize(localPoints.size());
3976 patchAttraction =
Zero;
3978 patchConstraints.setSize(localPoints.size());
3979 patchConstraints = pointConstraint();
3981 if (implicitFeatureAttraction)
3986 featureAttractionUsingReconstruction
3995 pointFaceSurfNormals,
4005 if (explicitFeatureAttraction)
4008 bool releasePoints =
false;
4009 bool stringFeatures =
false;
4010 bool avoidDiagonal =
false;
4013 releasePoints = snapParams.releasePoints();
4014 stringFeatures = snapParams.stringFeatures();
4015 avoidDiagonal = snapParams.avoidDiagonal();
4025 featureAttractionUsingFeatureEdges
4028 multiRegionFeatureSnap,
4030 snapParams.detectBaffles(),
4031 snapParams.baffleFeaturePoints(),
4045 pointFaceSurfNormals,
4055 if (!alignMeshEdges)
4059 degToRad(snapParams.concaveAngle())
4061 const scalar minAreaRatio = snapParams.minAreaRatio();
4063 Info<<
"Experimental: introducing face splits to avoid rotating"
4064 <<
" mesh edges. Splitting faces when" <<
nl
4065 <<
indent <<
"- angle not concave by more than "
4066 << snapParams.concaveAngle() <<
" degrees" <<
nl
4067 <<
indent <<
"- resulting triangles of similar area "
4068 <<
" (ratio within " << minAreaRatio <<
")" <<
nl
4089 Info<<
"Diagonal attraction feature correction : ";
4090 writeStats(pp, isPatchMasterPoint, patchConstraints);
4122 <<
" avg:" << avgPatchDisp <<
endl
4123 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
4124 <<
" avg:" << avgPatchAttr <<
endl;
4138 forAll(patchConstraints, pointi)
4140 if (patchConstraints[pointi].first() > 1)
4143 (1.0-featureAttract)*patchDisp[pointi]
4144 + featureAttract*patchAttraction[pointi];
4152 Info<<
"Feature analysis : ";
4153 writeStats(pp, isPatchMasterPoint, patchConstraints);
4168 if (featureAttract < 1-0.001)
4175 const bitSet isPatchMasterEdge
4207 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4218 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4225 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
4227 pp.localPoints() + tangPatchDisp
4232 (1.0-featureAttract)*smoothedPatchDisp
4233 + featureAttract*tangPatchDisp;
4237 const scalar
relax = featureAttract;
4248 minMagSqrEqOp<point>(),
4249 vector(GREAT, GREAT, GREAT)