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>(),
683 pointField localPoints(pp.points(), pp.meshPoints());
684 forAll(pointFaceCentres, pointi)
686 const point& pt = pp.points()[pp.meshPoints()[pointi]];
688 List<point>& pFc = pointFaceCentres[pointi];
699 listPlusEqOp<point>(),
703 forAll(pointFaceCentres, pointi)
705 const point& pt = pp.points()[pp.meshPoints()[pointi]];
707 List<point>& pFc = pointFaceCentres[pointi];
720 listPlusEqOp<label>(),
728 forAll(pointFaceDisp, pointi)
730 List<point>& pNormals = pointFaceSurfNormals[pointi];
731 List<point>& pDisp = pointFaceDisp[pointi];
732 List<point>& pFc = pointFaceCentres[pointi];
733 labelList& pFid = pointFacePatchID[pointi];
737 pNormals = List<point>(pNormals, visitOrder);
738 pDisp = List<point>(pDisp, visitOrder);
739 pFc = List<point>(pFc, visitOrder);
748 void Foam::snappySnapDriver::correctAttraction
750 const DynamicList<point>& surfacePoints,
751 const DynamicList<label>& surfaceCounts,
760 scalar tang = ((pt-edgePt)&edgeNormal);
764 if (order[0] < order[1])
768 vector attractD = surfacePoints[order[0]]-edgePt;
770 scalar tang2 = (attractD&edgeNormal);
772 attractD -= tang2*edgeNormal;
774 scalar magAttractD =
mag(attractD);
775 scalar fraction = magAttractD/(magAttractD+
mag(edgeOffset));
779 + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
780 edgeOffset = linePt-pt;
789 const List<point>& faceCentres
795 label patch0 = patchIDs[0];
797 for (label i = 1; i < patchIDs.size(); i++)
799 if (patchIDs[i] != patch0)
809 Foam::label Foam::snappySnapDriver::findNormal
811 const scalar featureCos,
813 const DynamicList<vector>& surfaceNormals
820 scalar cosAngle = (
n&surfaceNormals[j]);
824 (cosAngle >= featureCos)
825 || (cosAngle < (-1+0.001))
846 const DynamicList<vector>& surfaceNormals,
850 if (patchIDs.empty())
856 label patch0 = patchIDs[0];
858 for (label i = 1; i < patchIDs.size(); i++)
860 if (patchIDs[i] != patch0)
874 if (surfaceNormals.size() == 1)
882 labelList normalToPatch(surfaceNormals.size(), -1);
883 forAll(faceToNormalBin, i)
885 if (faceToNormalBin[i] != -1)
887 label&
patch = normalToPatch[faceToNormalBin[i]];
893 else if (
patch == -2)
897 else if (
patch != patchIDs[i])
905 forAll(normalToPatch, normali)
907 if (normalToPatch[normali] == -2)
922 void Foam::snappySnapDriver::writeStats
925 const bitSet& isPatchMasterPoint,
926 const List<pointConstraint>& patchConstraints
929 label nMasterPoints = 0;
934 forAll(patchConstraints, pointi)
936 if (isPatchMasterPoint[pointi])
940 if (patchConstraints[pointi].first() == 1)
944 else if (patchConstraints[pointi].first() == 2)
948 else if (patchConstraints[pointi].first() == 3)
955 reduce(nMasterPoints, sumOp<label>());
956 reduce(nPlanar, sumOp<label>());
957 reduce(nEdge, sumOp<label>());
958 reduce(nPoint, sumOp<label>());
959 Info<<
"total master points :" << nMasterPoints
960 <<
" of which attracted to :" <<
nl
961 <<
" feature point : " << nPoint <<
nl
962 <<
" feature edge : " << nEdge <<
nl
963 <<
" nearest surface : " << nPlanar <<
nl
964 <<
" rest : " << nMasterPoints-nPoint-nEdge-nPlanar
970 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
973 const scalar featureCos,
980 const List<List<point>>& pointFaceSurfNormals,
981 const List<List<point>>& pointFaceDisp,
982 const List<List<point>>& pointFaceCentres,
985 DynamicList<point>& surfacePoints,
986 DynamicList<vector>& surfaceNormals,
990 pointConstraint& patchConstraint
993 patchAttraction =
Zero;
994 patchConstraint = pointConstraint();
996 const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
997 const List<point>& pfDisp = pointFaceDisp[pointi];
998 const List<point>& pfCentres = pointFaceCentres[pointi];
1008 surfacePoints.clear();
1009 surfaceNormals.clear();
1012 faceToNormalBin.setSize(pfDisp.size());
1013 faceToNormalBin = -1;
1017 const point& fc = pfCentres[i];
1018 const vector& fSNormal = pfSurfNormals[i];
1019 const vector& fDisp = pfDisp[i];
1022 if (
magSqr(fDisp) <
sqr(snapDist[pointi]) &&
mag(fSNormal) > VSMALL)
1024 const point pt = fc + fDisp;
1027 faceToNormalBin[i] = findNormal
1034 if (faceToNormalBin[i] != -1)
1042 if (surfacePoints.size() <= 1)
1044 surfacePoints.append(pt);
1045 faceToNormalBin[i] = surfaceNormals.size();
1046 surfaceNormals.append(fSNormal);
1048 else if (surfacePoints.size() == 2)
1050 plane pl0(surfacePoints[0], surfaceNormals[0]);
1051 plane pl1(surfacePoints[1], surfaceNormals[1]);
1052 plane::ray r(pl0.planeIntersect(pl1));
1053 vector featureNormal = r.dir() /
mag(r.dir());
1055 if (
mag(fSNormal&featureNormal) >= 0.001)
1058 surfacePoints.append(pt);
1059 faceToNormalBin[i] = surfaceNormals.size();
1060 surfaceNormals.append(fSNormal);
1063 else if (surfacePoints.size() == 3)
1067 plane pl0(surfacePoints[0], surfaceNormals[0]);
1068 plane pl1(surfacePoints[1], surfaceNormals[1]);
1069 plane pl2(surfacePoints[2], surfaceNormals[2]);
1070 point p012(pl0.planePlaneIntersect(pl1, pl2));
1072 plane::ray r(pl0.planeIntersect(pl1));
1073 vector featureNormal = r.dir() /
mag(r.dir());
1074 if (
mag(fSNormal&featureNormal) >= 0.001)
1076 plane pl3(pt, fSNormal);
1077 point p013(pl0.planePlaneIntersect(pl1, pl3));
1079 if (
mag(p012-p013) > snapDist[pointi])
1082 surfacePoints.append(pt);
1083 faceToNormalBin[i] = surfaceNormals.
size();
1084 surfaceNormals.append(fSNormal);
1093 const point& pt = pp.localPoints()[pointi];
1096 if (surfaceNormals.size() == 1)
1100 ((surfacePoints[0]-pt) & surfaceNormals[0])
1109 patchAttraction = d;
1112 patchConstraint.applyConstraint(surfaceNormals[0]);
1114 else if (surfaceNormals.size() == 2)
1116 plane pl0(surfacePoints[0], surfaceNormals[0]);
1117 plane pl1(surfacePoints[1], surfaceNormals[1]);
1118 plane::ray r(pl0.planeIntersect(pl1));
1122 vector d = r.refPoint()-pt;
1131 patchAttraction = d;
1134 patchConstraint.applyConstraint(surfaceNormals[0]);
1135 patchConstraint.applyConstraint(surfaceNormals[1]);
1137 else if (surfaceNormals.size() == 3)
1140 plane pl0(surfacePoints[0], surfaceNormals[0]);
1141 plane pl1(surfacePoints[1], surfaceNormals[1]);
1142 plane pl2(surfacePoints[2], surfaceNormals[2]);
1143 point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
1144 vector d = cornerPt - pt;
1152 patchAttraction = d;
1155 patchConstraint.applyConstraint(surfaceNormals[0]);
1156 patchConstraint.applyConstraint(surfaceNormals[1]);
1157 patchConstraint.applyConstraint(surfaceNormals[2]);
1163 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
1166 const scalar featureCos,
1172 const List<List<point>>& pointFaceSurfNormals,
1173 const List<List<point>>& pointFaceDisp,
1174 const List<List<point>>& pointFaceCentres,
1178 List<pointConstraint>& patchConstraints
1181 autoPtr<OBJstream> feStr;
1182 autoPtr<OBJstream> fpStr;
1189 meshRefiner_.mesh().time().path()
1190 /
"implicitFeatureEdge_" +
name(iter) +
".obj"
1193 Info<<
"Dumping implicit feature-edge direction to "
1194 << feStr().name() <<
endl;
1200 meshRefiner_.mesh().time().path()
1201 /
"implicitFeaturePoint_" +
name(iter) +
".obj"
1204 Info<<
"Dumping implicit feature-point direction to "
1205 << fpStr().name() <<
endl;
1209 DynamicList<point> surfacePoints(4);
1210 DynamicList<vector> surfaceNormals(4);
1213 forAll(pp.localPoints(), pointi)
1216 pointConstraint constraint;
1218 featureAttractionUsingReconstruction
1229 pointFaceSurfNormals,
1244 (constraint.first() > patchConstraints[pointi].first())
1246 (constraint.first() == patchConstraints[pointi].first())
1247 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
1251 patchAttraction[pointi] = attraction;
1252 patchConstraints[pointi] = constraint;
1254 const point& pt = pp.localPoints()[pointi];
1256 if (feStr && patchConstraints[pointi].first() == 2)
1258 feStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1260 else if (fpStr && patchConstraints[pointi].first() == 3)
1262 fpStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1269 void Foam::snappySnapDriver::stringFeatureEdges
1272 const scalar featureCos,
1278 const List<pointConstraint>& rawPatchConstraints,
1281 List<pointConstraint>& patchConstraints
1311 forAll(pointEdges, pointi)
1313 if (patchConstraints[pointi].first() == 2)
1315 const point& pt = pp.localPoints()[pointi];
1316 const labelList& pEdges = pointEdges[pointi];
1317 const vector& featVec = patchConstraints[pointi].second();
1321 bool hasPos =
false;
1322 bool hasNeg =
false;
1326 const edge&
e = pp.edges()[pEdges[pEdgei]];
1327 label nbrPointi =
e.otherVertex(pointi);
1329 if (patchConstraints[nbrPointi].first() > 1)
1331 const point& nbrPt = pp.localPoints()[nbrPointi];
1332 const point featPt =
1333 nbrPt + patchAttraction[nbrPointi];
1334 const scalar cosAngle = (featVec & (featPt-pt));
1347 if (!hasPos || !hasNeg)
1353 label bestPosPointi = -1;
1354 scalar minPosDistSqr = GREAT;
1355 label bestNegPointi = -1;
1356 scalar minNegDistSqr = GREAT;
1360 const edge&
e = pp.edges()[pEdges[pEdgei]];
1361 label nbrPointi =
e.otherVertex(pointi);
1365 patchConstraints[nbrPointi].first() <= 1
1366 && rawPatchConstraints[nbrPointi].first() > 1
1369 const vector& nbrFeatVec =
1370 rawPatchConstraints[pointi].second();
1372 if (
mag(featVec&nbrFeatVec) > featureCos)
1379 rawPatchAttraction[nbrPointi]
1382 const point featPt =
1383 pp.localPoints()[nbrPointi]
1384 + rawPatchAttraction[nbrPointi];
1385 const scalar cosAngle =
1386 (featVec & (featPt-pt));
1390 if (!hasPos && d2 < minPosDistSqr)
1393 bestPosPointi = nbrPointi;
1398 if (!hasNeg && d2 < minNegDistSqr)
1401 bestNegPointi = nbrPointi;
1408 if (bestPosPointi != -1)
1418 patchAttraction[bestPosPointi] =
1419 0.5*rawPatchAttraction[bestPosPointi];
1420 patchConstraints[bestPosPointi] =
1421 rawPatchConstraints[bestPosPointi];
1425 if (bestNegPointi != -1)
1435 patchAttraction[bestNegPointi] =
1436 0.5*rawPatchAttraction[bestNegPointi];
1437 patchConstraints[bestNegPointi] =
1438 rawPatchConstraints[bestNegPointi];
1447 reduce(nChanged, sumOp<label>());
1448 Info<<
"Stringing feature edges : changed " << nChanged <<
" points"
1458 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1461 const scalar featureCos,
1466 const List<List<point>>& pointFaceCentres,
1470 const List<pointConstraint>& rawPatchConstraints,
1473 List<pointConstraint>& patchConstraints
1476 autoPtr<OBJstream> multiPatchStr;
1483 meshRefiner_.mesh().time().path()
1484 /
"multiPatch_" +
name(iter) +
".obj"
1487 Info<<
"Dumping removed constraints due to same-face"
1488 <<
" multi-patch points to "
1489 << multiPatchStr().name() <<
endl;
1494 bitSet isMultiPatchPoint(pp.size());
1496 forAll(pointFacePatchID, pointi)
1500 pp.localPoints()[pointi],
1501 pointFacePatchID[pointi],
1502 pointFaceCentres[pointi]
1504 isMultiPatchPoint.set(pointi, multiPatchPt.hit());
1508 forAll(isMultiPatchPoint, pointi)
1510 if (isMultiPatchPoint.test(pointi))
1514 patchConstraints[pointi].first() <= 1
1515 && rawPatchConstraints[pointi].first() > 1
1518 patchAttraction[pointi] = rawPatchAttraction[pointi];
1519 patchConstraints[pointi] = rawPatchConstraints[pointi];
1538 forAll(pp.localFaces(), facei)
1540 const face&
f = pp.localFaces()[facei];
1542 label nMultiPatchPoints = 0;
1545 label pointi =
f[fp];
1548 isMultiPatchPoint.test(pointi)
1549 && patchConstraints[pointi].first() > 1
1552 ++nMultiPatchPoints;
1556 if (nMultiPatchPoints > 0)
1560 label pointi =
f[fp];
1563 !isMultiPatchPoint.test(pointi)
1564 && patchConstraints[pointi].first() > 1
1570 patchAttraction[pointi] =
Zero;
1571 patchConstraints[pointi] = pointConstraint();
1576 multiPatchStr().write(pp.localPoints()[pointi]);
1583 reduce(nChanged, sumOp<label>());
1584 Info<<
"Removing constraints near multi-patch points : changed "
1585 << nChanged <<
" points" <<
endl;
1593 const List<pointConstraint>& patchConstraints,
1597 const face&
f = pp.localFaces()[facei];
1605 for (label startFp = 0; startFp <
f.size()-2; startFp++)
1607 label minFp =
f.rcIndex(startFp);
1611 label endFp =
f.fcIndex(
f.fcIndex(startFp));
1612 endFp <
f.size() && endFp != minFp;
1618 patchConstraints[
f[startFp]].first() >= 2
1619 && patchConstraints[
f[endFp]].first() >= 2
1622 attractIndices =
labelPair(startFp, endFp);
1628 return attractIndices;
1632 bool Foam::snappySnapDriver::isSplitAlignedWithFeature
1634 const scalar featureCos,
1636 const pointConstraint& pc0,
1638 const pointConstraint& pc1
1642 scalar magD =
mag(d);
1655 if (pc0.first() == 2 &&
mag(d & pc0.second()) > featureCos)
1659 else if (pc1.first() == 2 &&
mag(d & pc1.second()) > featureCos)
1672 bool Foam::snappySnapDriver::isConcave
1678 const scalar concaveCos
1682 scalar magN0 =
mag(n0);
1691 scalar d = (
c1-c0)&n0;
1702 scalar magN1 =
mag(n1);
1710 if ((n0&n1) < concaveCos)
1724 const scalar featureCos,
1725 const scalar concaveCos,
1726 const scalar minAreaRatio,
1729 const List<pointConstraint>& patchConstraints,
1735 DynamicField<point>& points1
1738 const face& localF = pp.localFaces()[facei];
1742 if (localF.size() >= 4)
1744 const pointField& localPts = pp.localPoints();
1782 for (label startFp = 0; startFp < localF.size()-2; startFp++)
1784 label minFp = localF.rcIndex(startFp);
1788 label endFp = localF.fcIndex(localF.fcIndex(startFp));
1789 endFp < localF.size() && endFp != minFp;
1793 label startPti = localF[startFp];
1794 label endPti = localF[endFp];
1796 const pointConstraint& startPc = patchConstraints[startPti];
1797 const pointConstraint& endPc = patchConstraints[endPti];
1799 if (startPc.first() >= 2 && endPc.first() >= 2)
1801 if (startPc.first() == 2 || endPc.first() == 2)
1806 point start = localPts[startPti]+patchAttr[startPti];
1807 point end = localPts[endPti]+patchAttr[endPti];
1811 !isSplitAlignedWithFeature
1837 f0.setSize(endFp-startFp+1);
1839 for (label fp = startFp; fp <= endFp; fp++)
1841 f0[i0++] = localF[fp];
1845 const face compact0(
identity(f0.size()));
1847 points0.append(localPts[f0[0]] + patchAttr[f0[0]]);
1848 for (label fp=1; fp < f0.size()-1; fp++)
1855 localPts[f0.last()] + patchAttr[f0.last()]
1861 f1.setSize(localF.size()+2-f0.size());
1867 fp = localF.fcIndex(fp)
1870 f1[i1++] = localF[fp];
1872 f1[i1++] = localF[startFp];
1876 const face compact1(
identity(f1.size()));
1878 points1.append(localPts[f1[0]] + patchAttr[f1[0]]);
1879 for (label fp=1; fp < f1.size()-1; fp++)
1882 points1.append(localPts[
pi] + nearestAttr[
pi]);
1886 localPts[f1.last()] + patchAttr[f1.last()]
1900 compact1.centre(points1),
1901 compact1.areaNormal(points1),
1921 const scalar area0 = f0.mag(localPts);
1922 const scalar area1 = f1.mag(localPts);
1926 area0/area1 >= minAreaRatio
1927 && area1/area0 >= minAreaRatio
1930 attractIndices =
labelPair(startFp, endFp);
1937 return attractIndices;
1941 void Foam::snappySnapDriver::splitDiagonals
1943 const scalar featureCos,
1944 const scalar concaveCos,
1945 const scalar minAreaRatio,
1952 List<pointConstraint>& patchConstraints,
1953 DynamicList<label>& splitFaces,
1954 DynamicList<labelPair>& splits
1957 const labelList& bFaces = pp.addressing();
1960 splitFaces.setCapacity(bFaces.size());
1962 splits.setCapacity(bFaces.size());
1966 DynamicField<point> facePoints0;
1967 DynamicField<point> facePoints1;
1973 findDiagonalAttraction
1994 splitFaces.append(bFaces[facei]);
1995 splits.append(
split);
1997 const face&
f = pp.localFaces()[facei];
2007 && patchConstraints[
f[fp]].first() >= 2
2010 patchConstraints[
f[fp]] = pointConstraint
2012 Tuple2<label, vector>
2015 nearestNormal[
f[fp]]
2018 patchAttraction[
f[fp]] = nearestAttraction[
f[fp]];
2026 void Foam::snappySnapDriver::avoidDiagonalAttraction
2029 const scalar featureCos,
2034 List<pointConstraint>& patchConstraints
2037 forAll(pp.localFaces(), facei)
2039 const face&
f = pp.localFaces()[facei];
2049 if (
diag[0] != -1 &&
diag[1] != -1)
2053 const label i0 =
f[
diag[0]];
2055 pp.localPoints()[i0]+patchAttraction[i0];
2056 const label i1 =
f[
diag[1]];
2058 pp.localPoints()[i1]+patchAttraction[i1];
2059 const point mid = 0.5*(pt0+pt1);
2061 const scalar cosAngle =
mag
2063 patchConstraints[i0].second()
2064 & patchConstraints[i1].second()
2073 if (cosAngle > featureCos)
2079 scalar minDistSqr = GREAT;
2082 label pointi =
f[fp];
2083 if (patchConstraints[pointi].first() <= 1)
2085 const point& pt = pp.localPoints()[pointi];
2086 scalar distSqr =
magSqr(mid-pt);
2087 if (distSqr < minDistSqr)
2095 label minPointi =
f[minFp];
2096 patchAttraction[minPointi] =
2097 mid-pp.localPoints()[minPointi];
2098 patchConstraints[minPointi] = patchConstraints[
f[
diag[0]]];
2124 Foam::snappySnapDriver::findNearFeatureEdge
2126 const bool isRegionEdge,
2131 const point& estimatedPt,
2133 List<List<DynamicList<point>>>& edgeAttractors,
2134 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2136 List<pointConstraint>& patchConstraints
2139 const refinementFeatures& features = meshRefiner_.features();
2142 List<pointIndexHit> nearEdgeInfo;
2147 features.findNearestRegionEdge
2158 features.findNearestEdge
2169 label feati = nearEdgeFeat[0];
2175 edgeAttractors[feati][nearInfo.index()].append
2179 pointConstraint
c(Tuple2<label, vector>(2, nearNormal[0]));
2180 edgeConstraints[feati][nearInfo.index()].append(
c);
2183 patchAttraction[pointi] = nearInfo.hitPoint()-pp.localPoints()[pointi];
2184 patchConstraints[pointi] =
c;
2186 return Tuple2<label, pointIndexHit>(feati, nearInfo);
2191 Foam::snappySnapDriver::findNearFeaturePoint
2193 const bool isRegionPoint,
2198 const point& estimatedPt,
2201 List<labelList>& pointAttractor,
2202 List<List<pointConstraint>>& pointConstraints,
2204 List<List<DynamicList<point>>>& edgeAttractors,
2205 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2208 List<pointConstraint>& patchConstraints
2211 const refinementFeatures& features = meshRefiner_.features();
2216 List<pointIndexHit> nearInfo;
2217 features.findNearestPoint
2225 label feati = nearFeat[0];
2229 const point& pt = pp.localPoints()[pointi];
2231 label featPointi = nearInfo[0].index();
2232 const point& featPt = nearInfo[0].hitPoint();
2233 scalar distSqr =
magSqr(featPt-pt);
2236 label oldPointi = pointAttractor[feati][featPointi];
2238 if (oldPointi != -1)
2241 if (distSqr >=
magSqr(featPt-pp.localPoints()[oldPointi]))
2250 pointAttractor[feati][featPointi] = pointi;
2251 pointConstraints[feati][featPointi].first() = 3;
2252 pointConstraints[feati][featPointi].second() =
Zero;
2255 patchAttraction[pointi] = featPt-pt;
2256 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2259 patchAttraction[oldPointi] =
Zero;
2260 patchConstraints[oldPointi] = pointConstraint();
2269 pp.localPoints()[oldPointi],
2281 pointAttractor[feati][featPointi] = pointi;
2282 pointConstraints[feati][featPointi].first() = 3;
2283 pointConstraints[feati][featPointi].second() =
Zero;
2286 patchAttraction[pointi] = featPt-pt;
2287 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2291 return Tuple2<label, pointIndexHit>(feati, nearInfo[0]);
2297 void Foam::snappySnapDriver::determineFeatures
2300 const scalar featureCos,
2301 const bool multiRegionFeatureSnap,
2307 const List<List<point>>& pointFaceSurfNormals,
2308 const List<List<point>>& pointFaceDisp,
2309 const List<List<point>>& pointFaceCentres,
2313 List<labelList>& pointAttractor,
2314 List<List<pointConstraint>>& pointConstraints,
2316 List<List<DynamicList<point>>>& edgeAttractors,
2317 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2320 List<pointConstraint>& patchConstraints
2323 autoPtr<OBJstream> featureEdgeStr;
2324 autoPtr<OBJstream> missedEdgeStr;
2325 autoPtr<OBJstream> featurePointStr;
2326 autoPtr<OBJstream> missedMP0Str;
2327 autoPtr<OBJstream> missedMP1Str;
2331 featureEdgeStr.reset
2335 meshRefiner_.mesh().time().path()
2336 /
"featureEdge_" +
name(iter) +
".obj"
2339 Info<<
"Dumping feature-edge sampling to "
2340 << featureEdgeStr().name() <<
endl;
2346 meshRefiner_.mesh().time().path()
2347 /
"missedFeatureEdge_" +
name(iter) +
".obj"
2350 Info<<
"Dumping feature-edges that are too far away to "
2351 << missedEdgeStr().name() <<
endl;
2353 featurePointStr.reset
2357 meshRefiner_.mesh().time().path()
2358 /
"featurePoint_" +
name(iter) +
".obj"
2361 Info<<
"Dumping feature-point sampling to "
2362 << featurePointStr().name() <<
endl;
2368 meshRefiner_.mesh().time().path()
2369 /
"missedFeatureEdgeFromMPEdge_" +
name(iter) +
".obj"
2372 Info<<
"Dumping region-edges that are too far away to "
2373 << missedMP0Str().name() <<
endl;
2379 meshRefiner_.mesh().time().path()
2380 /
"missedFeatureEdgeFromMPPoint_" +
name(iter) +
".obj"
2383 Info<<
"Dumping region-points that are too far away to "
2384 << missedMP1Str().name() <<
endl;
2388 DynamicList<point> surfacePoints(4);
2389 DynamicList<vector> surfaceNormals(4);
2392 forAll(pp.localPoints(), pointi)
2394 const point& pt = pp.localPoints()[pointi];
2406 pointConstraint constraint;
2408 featureAttractionUsingReconstruction
2419 pointFaceSurfNormals,
2464 (constraint.first() > patchConstraints[pointi].first())
2466 (constraint.first() == patchConstraints[pointi].first())
2467 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
2471 patchAttraction[pointi] = attraction;
2472 patchConstraints[pointi] = constraint;
2475 if (patchConstraints[pointi].first() == 1)
2478 if (multiRegionFeatureSnap)
2480 const point estimatedPt(pt + nearestDisp[pointi]);
2486 pointFacePatchID[pointi],
2492 if (multiPatchPt.hit())
2497 Tuple2<label, pointIndexHit> nearInfo =
2504 multiPatchPt.hitPoint(),
2519 featureEdgeStr().write
2529 missedEdgeStr().write
2538 else if (patchConstraints[pointi].first() == 2)
2543 const point estimatedPt(pt + patchAttraction[pointi]);
2548 bool hasSnapped =
false;
2549 if (multiRegionFeatureSnap)
2556 pointFacePatchID[pointi],
2561 if (multiPatchPt.hit())
2563 if (multiPatchPt.index() == 0)
2566 nearInfo = findNearFeatureEdge
2583 if (missedMP0Str && !nearInfo.second().hit())
2585 missedMP0Str().write
2599 nearInfo = findNearFeaturePoint
2628 if (!nearInfo.second().hit())
2630 nearInfo = findNearFeatureEdge
2648 if (missedMP1Str && !nearInfo.second().hit())
2650 missedMP1Str().write
2665 nearInfo = findNearFeatureEdge
2689 && patchConstraints[pointi].first() == 3
2692 featurePointStr().write
2700 && patchConstraints[pointi].first() == 2
2703 featureEdgeStr().write
2713 missedEdgeStr().write
2720 else if (patchConstraints[pointi].first() == 3)
2723 const point estimatedPt(pt + patchAttraction[pointi]);
2727 if (multiRegionFeatureSnap)
2734 pointFacePatchID[pointi],
2739 if (multiPatchPt.hit())
2742 nearInfo = findNearFeaturePoint
2763 nearInfo = findNearFeaturePoint
2786 nearInfo = findNearFeaturePoint
2807 if (featurePointStr && info.hit())
2809 featurePointStr().write
2831 void Foam::snappySnapDriver::determineBaffleFeatures
2834 const bool baffleFeaturePoints,
2835 const scalar featureCos,
2841 List<labelList>& pointAttractor,
2842 List<List<pointConstraint>>& pointConstraints,
2844 List<List<DynamicList<point>>>& edgeAttractors,
2845 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2848 List<pointConstraint>& patchConstraints
2851 const fvMesh&
mesh = meshRefiner_.mesh();
2852 const refinementFeatures& features = meshRefiner_.features();
2855 List<List<point>> edgeFaceNormals(pp.nEdges());
2858 forAll(pp.edgeFaces(), edgei)
2860 const labelList& eFaces = pp.edgeFaces()[edgei];
2861 List<point>& eFc = edgeFaceNormals[edgei];
2865 label facei = eFaces[i];
2866 eFc[i] = pp.faceNormals()[facei];
2881 listPlusEqOp<point>(),
2893 autoPtr<OBJstream> baffleEdgeStr;
2900 meshRefiner_.mesh().time().path()
2901 /
"baffleEdge_" +
name(iter) +
".obj"
2904 Info<<
nl <<
"Dumping baffle-edges to "
2905 << baffleEdgeStr().name() <<
endl;
2910 bitSet isBaffleEdge(pp.nEdges());
2911 label nBaffleEdges = 0;
2916 labelList pointStatus(pp.nPoints(), -1);
2918 forAll(edgeFaceNormals, edgei)
2920 const List<point>& efn = edgeFaceNormals[edgei];
2922 if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2924 isBaffleEdge.set(edgei);
2926 const edge&
e = pp.edges()[edgei];
2927 pointStatus[
e[0]] = 0;
2928 pointStatus[
e[1]] = 0;
2932 const point&
p0 = pp.localPoints()[
e[0]];
2933 const point& p1 = pp.localPoints()[
e[1]];
2939 reduce(nBaffleEdges, sumOp<label>());
2941 Info<<
"Detected " << nBaffleEdges
2942 <<
" baffle edges out of "
2944 <<
" edges." <<
endl;
2969 label nBafflePoints = 0;
2970 forAll(pointStatus, pointi)
2972 if (pointStatus[pointi] != -1)
2977 reduce(nBafflePoints, sumOp<label>());
2980 label nPointAttract = 0;
2981 label nEdgeAttract = 0;
2983 forAll(pointStatus, pointi)
2985 const point& pt = pp.localPoints()[pointi];
2987 if (pointStatus[pointi] == 0)
2991 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3010 if (nearInfo.second().hit())
3014 if (baffleFeaturePoints)
3016 nearInfo = findNearFeaturePoint
3036 if (nearInfo.first() != -1)
3044 else if (pointStatus[pointi] == 1)
3047 List<pointIndexHit> nearInfo;
3048 features.findNearestPoint
3056 label feati = nearFeat[0];
3062 label featPointi = nearInfo[0].index();
3063 const point& featPt = nearInfo[0].hitPoint();
3064 scalar distSqr =
magSqr(featPt-pt);
3067 label oldPointi = pointAttractor[feati][featPointi];
3074 <
magSqr(featPt-pp.localPoints()[oldPointi])
3078 pointAttractor[feati][featPointi] = pointi;
3079 pointConstraints[feati][featPointi].first() = 3;
3080 pointConstraints[feati][featPointi].second() =
Zero;
3083 patchAttraction[pointi] = featPt-pt;
3084 patchConstraints[pointi] =
3085 pointConstraints[feati][featPointi];
3087 if (oldPointi != -1)
3098 pp.localPoints()[oldPointi],
3123 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3137 if (nearInfo.first() != -1)
3145 reduce(nPointAttract, sumOp<label>());
3146 reduce(nEdgeAttract, sumOp<label>());
3148 Info<<
"Baffle points : " << nBafflePoints
3149 <<
" of which attracted to :" <<
nl
3150 <<
" feature point : " << nPointAttract <<
nl
3151 <<
" feature edge : " << nEdgeAttract <<
nl
3152 <<
" rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3158 void Foam::snappySnapDriver::reverseAttractMeshPoints
3166 const List<labelList>& pointAttractor,
3167 const List<List<pointConstraint>>& pointConstraints,
3169 const List<List<DynamicList<point>>>& edgeAttractors,
3170 const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3173 const List<pointConstraint>& rawPatchConstraints,
3177 List<pointConstraint>& patchConstraints
3180 const refinementFeatures& features = meshRefiner_.features();
3188 treeBoundBox bb(pp.localPoints());
3201 DynamicList<label> attractPoints(pp.nPoints());
3203 const fvMesh&
mesh = meshRefiner_.mesh();
3205 boolList isFeatureEdgeOrPoint(pp.nPoints(),
false);
3207 forAll(rawPatchConstraints, pointi)
3209 if (rawPatchConstraints[pointi].first() >= 2)
3211 isFeatureEdgeOrPoint[pointi] =
true;
3217 <<
" mesh points out of "
3219 <<
" for reverse attraction." <<
endl;
3227 isFeatureEdgeOrPoint,
3232 for (label nGrow = 0; nGrow < 1; nGrow++)
3234 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3236 forAll(pp.localFaces(), facei)
3238 const face&
f = pp.localFaces()[facei];
3242 if (isFeatureEdgeOrPoint[
f[fp]])
3247 newIsFeatureEdgeOrPoint[
f[fp]] =
true;
3254 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3260 isFeatureEdgeOrPoint,
3268 forAll(isFeatureEdgeOrPoint, pointi)
3270 if (isFeatureEdgeOrPoint[pointi])
3272 attractPoints.append(pointi);
3278 <<
" mesh points out of "
3280 <<
" for reverse attraction." <<
endl;
3284 indexedOctree<treeDataPoint> ppTree
3286 treeDataPoint(pp.localPoints(), attractPoints),
3294 patchAttraction.setSize(pp.nPoints());
3295 patchAttraction =
Zero;
3296 patchConstraints.setSize(pp.nPoints());
3297 patchConstraints = pointConstraint();
3299 forAll(edgeAttractors, feati)
3301 const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3302 const List<DynamicList<pointConstraint>>& edgeConstr =
3303 edgeConstraints[feati];
3305 forAll(edgeAttr, featEdgei)
3307 const DynamicList<point>& attr = edgeAttr[featEdgei];
3311 const point& featPt = attr[i];
3321 ppTree.shapes().pointLabels()[nearInfo.index()];
3322 const point attraction = featPt-pp.localPoints()[pointi];
3328 patchConstraints[pointi].first() <= 1
3329 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
3332 patchAttraction[pointi] = attraction;
3333 patchConstraints[pointi] = edgeConstr[featEdgei][i];
3338 static label nWarn = 0;
3343 <<
"Did not find pp point near " << featPt
3349 <<
"Reached warning limit " << nWarn
3350 <<
". Suppressing further warnings." <<
endl;
3369 forAll(pointAttractor, feati)
3371 const labelList& pointAttr = pointAttractor[feati];
3372 const List<pointConstraint>& pointConstr = pointConstraints[feati];
3374 forAll(pointAttr, featPointi)
3376 if (pointAttr[featPointi] != -1)
3378 const point& featPt = features[feati].points()
3393 ppTree.shapes().pointLabels()[nearInfo.index()];
3395 const point& pt = pp.localPoints()[pointi];
3396 const point attraction = featPt-pt;
3401 if (patchConstraints[pointi].first() <= 1)
3403 patchAttraction[pointi] = attraction;
3404 patchConstraints[pointi] = pointConstr[featPointi];
3406 else if (patchConstraints[pointi].first() == 2)
3408 patchAttraction[pointi] = attraction;
3409 patchConstraints[pointi] = pointConstr[featPointi];
3411 else if (patchConstraints[pointi].first() == 3)
3417 <
magSqr(patchAttraction[pointi])
3420 patchAttraction[pointi] = attraction;
3421 patchConstraints[pointi] =
3422 pointConstr[featPointi];
3432 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3435 const bool multiRegionFeatureSnap,
3437 const bool detectBaffles,
3438 const bool baffleFeaturePoints,
3440 const bool releasePoints,
3441 const bool stringFeatures,
3442 const bool avoidDiagonal,
3444 const scalar featureCos,
3451 const List<List<point>>& pointFaceSurfNormals,
3452 const List<List<point>>& pointFaceDisp,
3453 const List<List<point>>& pointFaceCentres,
3457 List<pointConstraint>& patchConstraints
3460 const refinementFeatures& features = meshRefiner_.features();
3461 const fvMesh&
mesh = meshRefiner_.mesh();
3463 const bitSet isPatchMasterPoint
3478 List<List<DynamicList<point>>> edgeAttractors(features.size());
3479 List<List<DynamicList<pointConstraint>>> edgeConstraints
3485 label nFeatEdges = features[feati].edges().size();
3486 edgeAttractors[feati].setSize(nFeatEdges);
3487 edgeConstraints[feati].setSize(nFeatEdges);
3493 List<labelList> pointAttractor(features.size());
3494 List<List<pointConstraint>> pointConstraints(features.size());
3497 label nFeatPoints = features[feati].points().size();
3498 pointAttractor[feati].setSize(nFeatPoints, -1);
3499 pointConstraints[feati].setSize(nFeatPoints);
3504 List<pointConstraint> rawPatchConstraints(pp.nPoints());
3510 multiRegionFeatureSnap,
3516 pointFaceSurfNormals,
3535 Info<<
"Raw geometric feature analysis : ";
3536 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3552 determineBaffleFeatures
3555 baffleFeaturePoints,
3576 Info<<
"After baffle feature analysis : ";
3577 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3585 reverseAttractMeshPoints
3601 rawPatchConstraints,
3611 Info<<
"Reverse attract feature analysis : ";
3612 writeStats(pp, isPatchMasterPoint, patchConstraints);
3618 OBJstream featureEdgeStr
3620 meshRefiner_.mesh().time().path()
3621 /
"edgeAttractors_" +
name(iter) +
".obj"
3623 Info<<
"Dumping feature-edge attraction to "
3624 << featureEdgeStr.name() <<
endl;
3626 OBJstream featurePointStr
3628 meshRefiner_.mesh().time().path()
3629 /
"pointAttractors_" +
name(iter) +
".obj"
3631 Info<<
"Dumping feature-point attraction to "
3632 << featurePointStr.name() <<
endl;
3634 forAll(patchConstraints, pointi)
3636 const point& pt = pp.localPoints()[pointi];
3637 const vector& attr = patchAttraction[pointi];
3639 if (patchConstraints[pointi].first() == 2)
3643 else if (patchConstraints[pointi].first() == 3)
3657 releasePointsNextToMultiPatch
3669 rawPatchConstraints,
3692 rawPatchConstraints,
3705 avoidDiagonalAttraction
3720 meshRefiner_.mesh().time().path()
3721 /
"patchAttraction_" +
name(iter) +
".obj",
3723 pp.localPoints() + patchAttraction
3730 void Foam::snappySnapDriver::preventFaceSqueeze
3733 const scalar featureCos,
3740 List<pointConstraint>& patchConstraints
3743 autoPtr<OBJstream> strPtr;
3750 meshRefiner_.mesh().time().path()
3751 /
"faceSqueeze_" +
name(iter) +
".obj"
3754 Info<<
"Dumping faceSqueeze corrections to "
3755 << strPtr().name() <<
endl;
3760 forAll(pp.localFaces(), facei)
3762 const face&
f = pp.localFaces()[facei];
3764 if (
f.size() !=
points.size())
3767 singleF.setSize(
f.size());
3768 for (label i = 0; i <
f.size(); i++)
3773 label nConstraints = 0;
3776 label pointi =
f[fp];
3777 const point& pt = pp.localPoints()[pointi];
3779 if (patchConstraints[pointi].first() > 1)
3781 points[fp] = pt + patchAttraction[pointi];
3790 if (nConstraints ==
f.size())
3801 const point& pt = pp.localPoints()[
f[fp]];
3802 const point& nextPt = pp.localPoints()[
f.nextLabel(fp)];
3813 label pointi =
f.prevLabel(maxFp);
3817 const point& pt = pp.localPoints()[pointi];
3824 patchAttraction[pointi] = nearestAttraction[pointi];
3837 scalar oldArea =
f.mag(pp.localPoints());
3838 scalar newArea = singleF.mag(
points);
3839 if (newArea < 0.1*oldArea)
3846 scalar
s =
magSqr(patchAttraction[
f[fp]]);
3855 label pointi =
f[maxFp];
3857 patchAttraction[pointi] *= 0.5;
3868 const snapParameters& snapParams,
3869 const bool alignMeshEdges,
3871 const scalar featureCos,
3872 const scalar featureAttract,
3876 motionSmoother& meshMover,
3878 List<pointConstraint>& patchConstraints,
3880 DynamicList<label>& splitFaces,
3881 DynamicList<labelPair>& splits
3891 const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3892 const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3893 const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3895 Info<<
"Overriding displacement on features :" <<
nl
3896 <<
" implicit features : " << implicitFeatureAttraction <<
nl
3897 <<
" explicit features : " << explicitFeatureAttraction <<
nl
3898 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl
3902 const pointField& localPoints = pp.localPoints();
3903 const fvMesh&
mesh = meshRefiner_.mesh();
3907 const bitSet isPatchMasterPoint
3920 List<List<point>> pointFaceSurfNormals;
3921 List<List<point>> pointFaceDisp;
3922 List<List<point>> pointFaceCentres;
3923 List<labelList> pointFacePatchID;
3929 forAll(pp.localFaces(), facei)
3931 const face&
f = pp.localFaces()[facei];
3934 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[
f[fp]]);
3946 labelList faceSurfaceGlobalRegion(pp.size(), -1);
3956 faceSurfaceGlobalRegion
3966 calcNearestFacePointProperties
3973 faceSurfaceGlobalRegion,
3975 pointFaceSurfNormals,
3994 patchAttraction.setSize(localPoints.size());
3995 patchAttraction =
Zero;
3997 patchConstraints.setSize(localPoints.size());
3998 patchConstraints = pointConstraint();
4000 if (implicitFeatureAttraction)
4005 featureAttractionUsingReconstruction
4014 pointFaceSurfNormals,
4024 if (explicitFeatureAttraction)
4027 bool releasePoints =
false;
4028 bool stringFeatures =
false;
4029 bool avoidDiagonal =
false;
4032 releasePoints = snapParams.releasePoints();
4033 stringFeatures = snapParams.stringFeatures();
4034 avoidDiagonal = snapParams.avoidDiagonal();
4044 featureAttractionUsingFeatureEdges
4047 multiRegionFeatureSnap,
4049 snapParams.detectBaffles(),
4050 snapParams.baffleFeaturePoints(),
4064 pointFaceSurfNormals,
4074 if (!alignMeshEdges)
4078 degToRad(snapParams.concaveAngle())
4080 const scalar minAreaRatio = snapParams.minAreaRatio();
4082 Info<<
"Experimental: introducing face splits to avoid rotating"
4083 <<
" mesh edges. Splitting faces when" <<
nl
4084 <<
indent <<
"- angle not concave by more than "
4085 << snapParams.concaveAngle() <<
" degrees" <<
nl
4086 <<
indent <<
"- resulting triangles of similar area "
4087 <<
" (ratio within " << minAreaRatio <<
")" <<
nl
4108 Info<<
"Diagonal attraction feature correction : ";
4109 writeStats(pp, isPatchMasterPoint, patchConstraints);
4141 <<
" avg:" << avgPatchDisp <<
endl
4142 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
4143 <<
" avg:" << avgPatchAttr <<
endl;
4157 forAll(patchConstraints, pointi)
4159 if (patchConstraints[pointi].first() > 1)
4162 (1.0-featureAttract)*patchDisp[pointi]
4163 + featureAttract*patchAttraction[pointi];
4171 Info<<
"Feature analysis : ";
4172 writeStats(pp, isPatchMasterPoint, patchConstraints);
4187 if (featureAttract < 1-0.001)
4194 const bitSet isPatchMasterEdge
4226 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4237 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4244 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
4246 pp.localPoints() + tangPatchDisp
4251 (1.0-featureAttract)*smoothedPatchDisp
4252 + featureAttract*tangPatchDisp;
4256 const scalar
relax = featureAttract;
4267 minMagSqrEqOp<point>(),
4268 vector(GREAT, GREAT, GREAT)