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];
2882 listPlusEqOp<point>(),
2895 autoPtr<OBJstream> baffleEdgeStr;
2902 meshRefiner_.mesh().time().path()
2903 /
"baffleEdge_" +
name(iter) +
".obj"
2906 Info<<
nl <<
"Dumping baffle-edges to "
2907 << baffleEdgeStr().name() <<
endl;
2912 bitSet isBaffleEdge(pp.nEdges());
2913 label nBaffleEdges = 0;
2918 labelList pointStatus(pp.nPoints(), -1);
2920 forAll(edgeFaceNormals, edgei)
2922 const List<point>& efn = edgeFaceNormals[edgei];
2924 if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2926 isBaffleEdge.set(edgei);
2928 const edge&
e = pp.edges()[edgei];
2929 pointStatus[
e[0]] = 0;
2930 pointStatus[
e[1]] = 0;
2934 const point&
p0 = pp.localPoints()[
e[0]];
2935 const point& p1 = pp.localPoints()[
e[1]];
2941 reduce(nBaffleEdges, sumOp<label>());
2943 Info<<
"Detected " << nBaffleEdges
2944 <<
" baffle edges out of "
2946 <<
" edges." <<
endl;
2971 label nBafflePoints = 0;
2972 forAll(pointStatus, pointi)
2974 if (pointStatus[pointi] != -1)
2979 reduce(nBafflePoints, sumOp<label>());
2982 label nPointAttract = 0;
2983 label nEdgeAttract = 0;
2985 forAll(pointStatus, pointi)
2987 const point& pt = pp.localPoints()[pointi];
2989 if (pointStatus[pointi] == 0)
2993 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3012 if (nearInfo.second().hit())
3016 if (baffleFeaturePoints)
3018 nearInfo = findNearFeaturePoint
3038 if (nearInfo.first() != -1)
3046 else if (pointStatus[pointi] == 1)
3049 List<pointIndexHit> nearInfo;
3050 features.findNearestPoint
3058 label feati = nearFeat[0];
3064 label featPointi = nearInfo[0].index();
3065 const point& featPt = nearInfo[0].hitPoint();
3066 scalar distSqr =
magSqr(featPt-pt);
3069 label oldPointi = pointAttractor[feati][featPointi];
3076 <
magSqr(featPt-pp.localPoints()[oldPointi])
3080 pointAttractor[feati][featPointi] = pointi;
3081 pointConstraints[feati][featPointi].first() = 3;
3082 pointConstraints[feati][featPointi].second() =
Zero;
3085 patchAttraction[pointi] = featPt-pt;
3086 patchConstraints[pointi] =
3087 pointConstraints[feati][featPointi];
3089 if (oldPointi != -1)
3100 pp.localPoints()[oldPointi],
3125 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3139 if (nearInfo.first() != -1)
3147 reduce(nPointAttract, sumOp<label>());
3148 reduce(nEdgeAttract, sumOp<label>());
3150 Info<<
"Baffle points : " << nBafflePoints
3151 <<
" of which attracted to :" <<
nl
3152 <<
" feature point : " << nPointAttract <<
nl
3153 <<
" feature edge : " << nEdgeAttract <<
nl
3154 <<
" rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3160 void Foam::snappySnapDriver::reverseAttractMeshPoints
3168 const List<labelList>& pointAttractor,
3169 const List<List<pointConstraint>>& pointConstraints,
3171 const List<List<DynamicList<point>>>& edgeAttractors,
3172 const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3175 const List<pointConstraint>& rawPatchConstraints,
3179 List<pointConstraint>& patchConstraints
3182 const refinementFeatures& features = meshRefiner_.features();
3190 treeBoundBox bb(pp.localPoints());
3203 DynamicList<label> attractPoints(pp.nPoints());
3205 const fvMesh&
mesh = meshRefiner_.mesh();
3207 boolList isFeatureEdgeOrPoint(pp.nPoints(),
false);
3209 forAll(rawPatchConstraints, pointi)
3211 if (rawPatchConstraints[pointi].first() >= 2)
3213 isFeatureEdgeOrPoint[pointi] =
true;
3219 <<
" mesh points out of "
3221 <<
" for reverse attraction." <<
endl;
3229 isFeatureEdgeOrPoint,
3234 for (label nGrow = 0; nGrow < 1; nGrow++)
3236 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3238 forAll(pp.localFaces(), facei)
3240 const face&
f = pp.localFaces()[facei];
3244 if (isFeatureEdgeOrPoint[
f[fp]])
3249 newIsFeatureEdgeOrPoint[
f[fp]] =
true;
3256 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3262 isFeatureEdgeOrPoint,
3270 forAll(isFeatureEdgeOrPoint, pointi)
3272 if (isFeatureEdgeOrPoint[pointi])
3274 attractPoints.append(pointi);
3280 <<
" mesh points out of "
3282 <<
" for reverse attraction." <<
endl;
3286 indexedOctree<treeDataPoint> ppTree
3288 treeDataPoint(pp.localPoints(), attractPoints),
3296 patchAttraction.setSize(pp.nPoints());
3297 patchAttraction =
Zero;
3298 patchConstraints.setSize(pp.nPoints());
3299 patchConstraints = pointConstraint();
3301 forAll(edgeAttractors, feati)
3303 const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3304 const List<DynamicList<pointConstraint>>& edgeConstr =
3305 edgeConstraints[feati];
3307 forAll(edgeAttr, featEdgei)
3309 const DynamicList<point>& attr = edgeAttr[featEdgei];
3313 const point& featPt = attr[i];
3323 ppTree.shapes().pointLabels()[nearInfo.index()];
3324 const point attraction = featPt-pp.localPoints()[pointi];
3330 patchConstraints[pointi].first() <= 1
3331 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
3334 patchAttraction[pointi] = attraction;
3335 patchConstraints[pointi] = edgeConstr[featEdgei][i];
3340 static label nWarn = 0;
3345 <<
"Did not find pp point near " << featPt
3351 <<
"Reached warning limit " << nWarn
3352 <<
". Suppressing further warnings." <<
endl;
3371 forAll(pointAttractor, feati)
3373 const labelList& pointAttr = pointAttractor[feati];
3374 const List<pointConstraint>& pointConstr = pointConstraints[feati];
3376 forAll(pointAttr, featPointi)
3378 if (pointAttr[featPointi] != -1)
3380 const point& featPt = features[feati].points()
3395 ppTree.shapes().pointLabels()[nearInfo.index()];
3397 const point& pt = pp.localPoints()[pointi];
3398 const point attraction = featPt-pt;
3403 if (patchConstraints[pointi].first() <= 1)
3405 patchAttraction[pointi] = attraction;
3406 patchConstraints[pointi] = pointConstr[featPointi];
3408 else if (patchConstraints[pointi].first() == 2)
3410 patchAttraction[pointi] = attraction;
3411 patchConstraints[pointi] = pointConstr[featPointi];
3413 else if (patchConstraints[pointi].first() == 3)
3419 <
magSqr(patchAttraction[pointi])
3422 patchAttraction[pointi] = attraction;
3423 patchConstraints[pointi] =
3424 pointConstr[featPointi];
3434 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3437 const bool multiRegionFeatureSnap,
3439 const bool detectBaffles,
3440 const bool baffleFeaturePoints,
3442 const bool releasePoints,
3443 const bool stringFeatures,
3444 const bool avoidDiagonal,
3446 const scalar featureCos,
3453 const List<List<point>>& pointFaceSurfNormals,
3454 const List<List<point>>& pointFaceDisp,
3455 const List<List<point>>& pointFaceCentres,
3459 List<pointConstraint>& patchConstraints
3462 const refinementFeatures& features = meshRefiner_.features();
3463 const fvMesh&
mesh = meshRefiner_.mesh();
3465 const bitSet isPatchMasterPoint
3480 List<List<DynamicList<point>>> edgeAttractors(features.size());
3481 List<List<DynamicList<pointConstraint>>> edgeConstraints
3487 label nFeatEdges = features[feati].edges().size();
3488 edgeAttractors[feati].setSize(nFeatEdges);
3489 edgeConstraints[feati].setSize(nFeatEdges);
3495 List<labelList> pointAttractor(features.size());
3496 List<List<pointConstraint>> pointConstraints(features.size());
3499 label nFeatPoints = features[feati].points().size();
3500 pointAttractor[feati].setSize(nFeatPoints, -1);
3501 pointConstraints[feati].setSize(nFeatPoints);
3506 List<pointConstraint> rawPatchConstraints(pp.nPoints());
3512 multiRegionFeatureSnap,
3518 pointFaceSurfNormals,
3537 Info<<
"Raw geometric feature analysis : ";
3538 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3554 determineBaffleFeatures
3557 baffleFeaturePoints,
3578 Info<<
"After baffle feature analysis : ";
3579 writeStats(pp, isPatchMasterPoint, rawPatchConstraints);
3587 reverseAttractMeshPoints
3603 rawPatchConstraints,
3613 Info<<
"Reverse attract feature analysis : ";
3614 writeStats(pp, isPatchMasterPoint, patchConstraints);
3620 OBJstream featureEdgeStr
3622 meshRefiner_.mesh().time().path()
3623 /
"edgeAttractors_" +
name(iter) +
".obj"
3625 Info<<
"Dumping feature-edge attraction to "
3626 << featureEdgeStr.name() <<
endl;
3628 OBJstream featurePointStr
3630 meshRefiner_.mesh().time().path()
3631 /
"pointAttractors_" +
name(iter) +
".obj"
3633 Info<<
"Dumping feature-point attraction to "
3634 << featurePointStr.name() <<
endl;
3636 forAll(patchConstraints, pointi)
3638 const point& pt = pp.localPoints()[pointi];
3639 const vector& attr = patchAttraction[pointi];
3641 if (patchConstraints[pointi].first() == 2)
3645 else if (patchConstraints[pointi].first() == 3)
3659 releasePointsNextToMultiPatch
3671 rawPatchConstraints,
3694 rawPatchConstraints,
3707 avoidDiagonalAttraction
3722 meshRefiner_.mesh().time().path()
3723 /
"patchAttraction_" +
name(iter) +
".obj",
3725 pp.localPoints() + patchAttraction
3732 void Foam::snappySnapDriver::preventFaceSqueeze
3735 const scalar featureCos,
3742 List<pointConstraint>& patchConstraints
3745 autoPtr<OBJstream> strPtr;
3752 meshRefiner_.mesh().time().path()
3753 /
"faceSqueeze_" +
name(iter) +
".obj"
3756 Info<<
"Dumping faceSqueeze corrections to "
3757 << strPtr().name() <<
endl;
3762 forAll(pp.localFaces(), facei)
3764 const face&
f = pp.localFaces()[facei];
3766 if (
f.size() !=
points.size())
3769 singleF.setSize(
f.size());
3770 for (label i = 0; i <
f.size(); i++)
3775 label nConstraints = 0;
3778 label pointi =
f[fp];
3779 const point& pt = pp.localPoints()[pointi];
3781 if (patchConstraints[pointi].first() > 1)
3783 points[fp] = pt + patchAttraction[pointi];
3792 if (nConstraints ==
f.size())
3803 const point& pt = pp.localPoints()[
f[fp]];
3804 const point& nextPt = pp.localPoints()[
f.nextLabel(fp)];
3815 label pointi =
f.prevLabel(maxFp);
3819 const point& pt = pp.localPoints()[pointi];
3826 patchAttraction[pointi] = nearestAttraction[pointi];
3839 scalar oldArea =
f.mag(pp.localPoints());
3840 scalar newArea = singleF.mag(
points);
3841 if (newArea < 0.1*oldArea)
3848 scalar
s =
magSqr(patchAttraction[
f[fp]]);
3857 label pointi =
f[maxFp];
3859 patchAttraction[pointi] *= 0.5;
3870 const snapParameters& snapParams,
3871 const bool alignMeshEdges,
3873 const scalar featureCos,
3874 const scalar featureAttract,
3878 motionSmoother& meshMover,
3880 List<pointConstraint>& patchConstraints,
3882 DynamicList<label>& splitFaces,
3883 DynamicList<labelPair>& splits
3893 const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3894 const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3895 const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3897 Info<<
"Overriding displacement on features :" <<
nl
3898 <<
" implicit features : " << implicitFeatureAttraction <<
nl
3899 <<
" explicit features : " << explicitFeatureAttraction <<
nl
3900 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl
3904 const pointField& localPoints = pp.localPoints();
3905 const fvMesh&
mesh = meshRefiner_.mesh();
3909 const bitSet isPatchMasterPoint
3922 List<List<point>> pointFaceSurfNormals;
3923 List<List<point>> pointFaceDisp;
3924 List<List<point>> pointFaceCentres;
3925 List<labelList> pointFacePatchID;
3931 forAll(pp.localFaces(), facei)
3933 const face&
f = pp.localFaces()[facei];
3936 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[
f[fp]]);
3948 labelList faceSurfaceGlobalRegion(pp.size(), -1);
3958 faceSurfaceGlobalRegion
3968 calcNearestFacePointProperties
3975 faceSurfaceGlobalRegion,
3977 pointFaceSurfNormals,
3996 patchAttraction.setSize(localPoints.size());
3997 patchAttraction =
Zero;
3999 patchConstraints.setSize(localPoints.size());
4000 patchConstraints = pointConstraint();
4002 if (implicitFeatureAttraction)
4007 featureAttractionUsingReconstruction
4016 pointFaceSurfNormals,
4026 if (explicitFeatureAttraction)
4029 bool releasePoints =
false;
4030 bool stringFeatures =
false;
4031 bool avoidDiagonal =
false;
4034 releasePoints = snapParams.releasePoints();
4035 stringFeatures = snapParams.stringFeatures();
4036 avoidDiagonal = snapParams.avoidDiagonal();
4046 featureAttractionUsingFeatureEdges
4049 multiRegionFeatureSnap,
4051 snapParams.detectBaffles(),
4052 snapParams.baffleFeaturePoints(),
4066 pointFaceSurfNormals,
4076 if (!alignMeshEdges)
4080 degToRad(snapParams.concaveAngle())
4082 const scalar minAreaRatio = snapParams.minAreaRatio();
4084 Info<<
"Experimental: introducing face splits to avoid rotating"
4085 <<
" mesh edges. Splitting faces when" <<
nl
4086 <<
indent <<
"- angle not concave by more than "
4087 << snapParams.concaveAngle() <<
" degrees" <<
nl
4088 <<
indent <<
"- resulting triangles of similar area "
4089 <<
" (ratio within " << minAreaRatio <<
")" <<
nl
4110 Info<<
"Diagonal attraction feature correction : ";
4111 writeStats(pp, isPatchMasterPoint, patchConstraints);
4143 <<
" avg:" << avgPatchDisp <<
endl
4144 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
4145 <<
" avg:" << avgPatchAttr <<
endl;
4159 forAll(patchConstraints, pointi)
4161 if (patchConstraints[pointi].first() > 1)
4164 (1.0-featureAttract)*patchDisp[pointi]
4165 + featureAttract*patchAttraction[pointi];
4173 Info<<
"Feature analysis : ";
4174 writeStats(pp, isPatchMasterPoint, patchConstraints);
4189 if (featureAttract < 1-0.001)
4196 const bitSet isPatchMasterEdge
4228 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4239 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4246 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
4248 pp.localPoints() + tangPatchDisp
4253 (1.0-featureAttract)*smoothedPatchDisp
4254 + featureAttract*tangPatchDisp;
4258 const scalar
relax = featureAttract;
4269 minMagSqrEqOp<point>(),
4270 vector(GREAT, GREAT, GREAT)