66 bool Foam::snappySnapDriver::isFeaturePoint
68 const scalar featureCos,
82 if (isFeatureEdge[pEdges[i]])
86 for (
label j = i+1; j < pEdges.
size(); j++)
88 if (isFeatureEdge[pEdges[j]])
90 const edge& eI = edges[pEdges[i]];
91 const edge& eJ = edges[pEdges[j]];
93 const point&
p = points[pointi];
98 scalar vIMag =
mag(vI);
101 scalar vJMag =
mag(vJ);
107 && ((vI/vIMag & vJ/vJMag) < featureCos)
127 void Foam::snappySnapDriver::smoothAndConstrain
138 for (
label avgIter = 0; avgIter < 20; avgIter++)
159 forAll(pointEdges, pointi)
161 const labelList& pEdges = pointEdges[pointi];
163 label nConstraints = constraints[pointi].
first();
165 if (nConstraints <= 1)
169 label edgeI = pEdges[i];
171 if (isPatchMasterEdge[edgeI])
173 label nbrPointi = edges[edgeI].otherVertex(pointi);
174 if (constraints[nbrPointi].first() >= nConstraints)
176 dispSum[pointi] += disp[nbrPointi];
204 forAll(constraints, pointi)
206 if (dispCount[pointi] > 0)
211 *(disp[pointi] + dispSum[pointi]/dispCount[pointi]);
218 void Foam::snappySnapDriver::calcNearestFace
235 faceSurfaceNormal.
setSize(pp.size());
236 faceSurfaceNormal =
Zero;
237 faceSurfaceGlobalRegion.
setSize(pp.size());
238 faceSurfaceGlobalRegion = -1;
259 label zoneSurfI = zonedSurfaces[i];
261 const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
268 <<
"Problem. Cannot find zone " << faceZoneName
275 isZonedFace[fZone[i]] = 1;
280 forAll(pp.addressing(), i)
282 if (isZonedFace[pp.addressing()[i]])
284 snapSurf[i] = zoneSurfI;
286 meshFaces.append(pp.addressing()[i]);
320 if (hitInfo[hitI].hit())
322 label facei = ppFaces[hitI];
323 faceDisp[facei] = hitInfo[hitI].hitPoint() - fc[hitI];
324 faceSurfaceNormal[facei] = hitNormal[hitI];
334 <<
"Did not find surface near face centre " << fc[hitI]
347 forAll(pp.addressing(), i)
349 if (snapSurf[i] == -1)
352 meshFaces.append(pp.addressing()[i]);
384 if (hitInfo[hitI].hit())
386 label facei = ppFaces[hitI];
387 faceDisp[facei] = hitInfo[hitI].hitPoint() - fc[hitI];
388 faceSurfaceNormal[facei] = hitNormal[hitI];
398 <<
"Did not find surface near face centre " << fc[hitI]
408 faceRotation.
setSize(pp.size());
411 forAll(faceRotation, facei)
414 faceRotation[facei] =
416 ^ faceSurfaceNormal[facei];
424 /
"faceDisp_" +
name(iter) +
".obj",
431 /
"faceRotation_" +
name(iter) +
".obj",
439 void Foam::snappySnapDriver::calcNearestFacePointProperties
446 const labelList& faceSurfaceGlobalRegion,
463 pointFaceSurfNormals.setSize(pp.
nPoints());
464 pointFaceDisp.setSize(pp.
nPoints());
465 pointFaceCentres.setSize(pp.
nPoints());
477 label facei = pFaces[i];
478 if (isMasterFace[facei] && faceSurfaceGlobalRegion[facei] != -1)
485 List<point>& pNormals = pointFaceSurfNormals[pointi];
491 labelList& pFid = pointFacePatchID[pointi];
497 label facei = pFaces[i];
498 label globalRegionI = faceSurfaceGlobalRegion[facei];
500 if (isMasterFace[facei] && globalRegionI != -1)
502 pNormals[nFaces] = faceSurfaceNormal[facei];
503 pDisp[nFaces] = faceDisp[facei];
505 pFid[nFaces] = globalToMasterPatch_[globalRegionI];
528 if (pp.
coupled() || isA<emptyPolyPatch>(pp))
539 forAll(pp.addressing(), i)
541 label meshFacei = pp.addressing()[i];
550 label pointi = boundaryPoints[i];
555 List<point>& pNormals = pointFaceSurfNormals[pointi];
558 labelList& pFid = pointFacePatchID[pointi];
562 label meshFacei = pFaces[i];
584 pointFaceSurfNormals,
620 forAll(pointFaceDisp, pointi)
622 List<point>& pNormals = pointFaceSurfNormals[pointi];
625 labelList& pFid = pointFacePatchID[pointi];
637 void Foam::snappySnapDriver::correctAttraction
649 scalar tang = ((pt-edgePt)&edgeNormal);
654 if (order[0] < order[1])
658 vector attractD = surfacePoints[order[0]]-edgePt;
660 scalar tang2 = (attractD&edgeNormal);
662 attractD -= tang2*edgeNormal;
664 scalar magAttractD =
mag(attractD);
665 scalar fraction = magAttractD/(magAttractD+
mag(edgeOffset));
669 + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
670 edgeOffset = linePt-pt;
685 label patch0 = patchIDs[0];
687 for (
label i = 1; i < patchIDs.
size(); i++)
689 if (patchIDs[i] != patch0)
701 const scalar featureCos,
710 scalar cosAngle = (n&surfaceNormals[j]);
714 (cosAngle >= featureCos)
715 || (cosAngle < (-1+0.001))
734 if (patchIDs.
empty())
740 label patch0 = patchIDs[0];
742 for (
label i = 1; i < patchIDs.
size(); i++)
744 if (patchIDs[i] != patch0)
758 if (surfaceNormals.
size() == 1)
767 forAll(faceToNormalBin, i)
769 if (faceToNormalBin[i] != -1)
771 label& patch = normalToPatch[faceToNormalBin[i]];
777 else if (patch == -2)
781 else if (patch != patchIDs[i])
789 forAll(normalToPatch, normalI)
791 if (normalToPatch[normalI] == -2)
806 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
809 const scalar featureCos,
829 patchAttraction =
Zero;
832 const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
834 const List<point>& pfCentres = pointFaceCentres[pointi];
844 surfacePoints.
clear();
845 surfaceNormals.
clear();
849 faceToNormalBin = -1;
853 const point& fc = pfCentres[i];
854 const vector& fSNormal = pfSurfNormals[i];
855 const vector& fDisp = pfDisp[i];
858 if (
magSqr(fDisp) <
sqr(snapDist[pointi]) &&
mag(fSNormal) > VSMALL)
860 const point pt = fc + fDisp;
863 faceToNormalBin[i] = findNormal
870 if (faceToNormalBin[i] != -1)
878 if (surfacePoints.
size() <= 1)
881 faceToNormalBin[i] = surfaceNormals.
size();
882 surfaceNormals.
append(fSNormal);
884 else if (surfacePoints.
size() == 2)
886 plane pl0(surfacePoints[0], surfaceNormals[0]);
887 plane pl1(surfacePoints[1], surfaceNormals[1]);
889 vector featureNormal = r.dir() /
mag(r.dir());
891 if (
mag(fSNormal&featureNormal) >= 0.001)
895 faceToNormalBin[i] = surfaceNormals.
size();
896 surfaceNormals.
append(fSNormal);
899 else if (surfacePoints.
size() == 3)
903 plane pl0(surfacePoints[0], surfaceNormals[0]);
904 plane pl1(surfacePoints[1], surfaceNormals[1]);
905 plane pl2(surfacePoints[2], surfaceNormals[2]);
909 vector featureNormal = r.dir() /
mag(r.dir());
910 if (
mag(fSNormal&featureNormal) >= 0.001)
912 plane pl3(pt, fSNormal);
915 if (
mag(p012-p013) > snapDist[pointi])
919 faceToNormalBin[i] = surfaceNormals.
size();
920 surfaceNormals.
append(fSNormal);
932 if (surfaceNormals.
size() == 1)
936 ((surfacePoints[0]-pt) & surfaceNormals[0])
950 else if (surfaceNormals.
size() == 2)
952 plane pl0(surfacePoints[0], surfaceNormals[0]);
953 plane pl1(surfacePoints[1], surfaceNormals[1]);
958 vector d = r.refPoint()-pt;
973 else if (surfaceNormals.
size() == 3)
976 plane pl0(surfacePoints[0], surfaceNormals[0]);
977 plane pl1(surfacePoints[1], surfaceNormals[1]);
978 plane pl2(surfacePoints[2], surfaceNormals[2]);
998 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
1001 const bool avoidSnapProblems,
1002 const scalar featureCos,
1025 meshRefiner_.mesh().time().path()
1026 /
"implicitFeatureEdge_" +
name(iter) +
".obj" 1029 Info<<
"Dumping implicit feature-edge direction to " 1030 << feStr().name() <<
endl;
1036 meshRefiner_.mesh().time().path()
1037 /
"implicitFeaturePoint_" +
name(iter) +
".obj" 1040 Info<<
"Dumping implicit feature-point direction to " 1041 << fpStr().name() <<
endl;
1054 featureAttractionUsingReconstruction
1065 pointFaceSurfNormals,
1080 (constraint.
first() > patchConstraints[pointi].
first())
1082 (constraint.
first() == patchConstraints[pointi].
first())
1083 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
1087 patchAttraction[pointi] = attraction;
1088 patchConstraints[pointi] = constraint;
1092 if (patchConstraints[pointi].first() == 2 && feStr.
valid())
1094 feStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1096 else if (patchConstraints[pointi].first() == 3 && fpStr.
valid())
1098 fpStr().write(
linePointRef(pt, pt+patchAttraction[pointi]));
1105 void Foam::snappySnapDriver::stringFeatureEdges
1108 const scalar featureCos,
1147 forAll(pointEdges, pointi)
1149 if (patchConstraints[pointi].first() == 2)
1152 const labelList& pEdges = pointEdges[pointi];
1153 const vector& featVec = patchConstraints[pointi].second();
1157 bool hasPos =
false;
1158 bool hasNeg =
false;
1165 if (patchConstraints[nbrPointi].first() > 1)
1168 const point featPt =
1169 nbrPt + patchAttraction[nbrPointi];
1170 const scalar cosAngle = (featVec & (featPt-pt));
1183 if (!hasPos || !hasNeg)
1189 label bestPosPointi = -1;
1190 scalar minPosDistSqr = GREAT;
1191 label bestNegPointi = -1;
1192 scalar minNegDistSqr = GREAT;
1201 patchConstraints[nbrPointi].first() <= 1
1202 && rawPatchConstraints[nbrPointi].first() > 1
1205 const vector& nbrFeatVec =
1206 rawPatchConstraints[pointi].second();
1208 if (
mag(featVec&nbrFeatVec) > featureCos)
1215 rawPatchAttraction[nbrPointi]
1218 const point featPt =
1220 + rawPatchAttraction[nbrPointi];
1221 const scalar cosAngle =
1222 (featVec & (featPt-pt));
1226 if (!hasPos && d2 < minPosDistSqr)
1229 bestPosPointi = nbrPointi;
1234 if (!hasNeg && d2 < minNegDistSqr)
1237 bestNegPointi = nbrPointi;
1244 if (bestPosPointi != -1)
1254 patchAttraction[bestPosPointi] =
1255 0.5*rawPatchAttraction[bestPosPointi];
1256 patchConstraints[bestPosPointi] =
1257 rawPatchConstraints[bestPosPointi];
1261 if (bestNegPointi != -1)
1271 patchAttraction[bestNegPointi] =
1272 0.5*rawPatchAttraction[bestNegPointi];
1273 patchConstraints[bestNegPointi] =
1274 rawPatchConstraints[bestNegPointi];
1284 Info<<
"Stringing feature edges : changed " << nChanged <<
" points" 1294 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1297 const scalar featureCos,
1319 meshRefiner_.mesh().time().path()
1320 /
"multiPatch_" +
name(iter) +
".obj" 1323 Info<<
"Dumping removed constraints due to same-face" 1324 <<
" multi-patch points to " 1325 << multiPatchStr().name() <<
endl;
1332 forAll(pointFacePatchID, pointi)
1337 pointFacePatchID[pointi],
1338 pointFaceCentres[pointi]
1340 isMultiPatchPoint[pointi] = multiPatchPt.
hit();
1344 forAll(isMultiPatchPoint, pointi)
1346 if (isMultiPatchPoint[pointi])
1350 patchConstraints[pointi].first() <= 1
1351 && rawPatchConstraints[pointi].first() > 1
1354 patchAttraction[pointi] = rawPatchAttraction[pointi];
1355 patchConstraints[pointi] = rawPatchConstraints[pointi];
1378 label nMultiPatchPoints = 0;
1381 label pointi = f[fp];
1384 isMultiPatchPoint[pointi]
1385 && patchConstraints[pointi].first() > 1
1388 nMultiPatchPoints++;
1392 if (nMultiPatchPoints > 0)
1396 label pointi = f[fp];
1399 !isMultiPatchPoint[pointi]
1400 && patchConstraints[pointi].first() > 1
1406 patchAttraction[pointi] =
Zero;
1410 if (multiPatchStr.
valid())
1420 Info<<
"Removing constraints near multi-patch points : changed " 1421 << nChanged <<
" points" <<
endl;
1444 label pointi = f[fp];
1445 if (patchConstraints[pointi].first() >= 2)
1448 if (attractIndices[0] == -1)
1451 attractIndices[0] = fp;
1453 else if (attractIndices[1] == -1)
1457 label fp0 = attractIndices[0];
1466 attractIndices[1] = fp;
1479 if (attractIndices[1] == -1)
1485 return attractIndices;
1489 void Foam::snappySnapDriver::avoidDiagonalAttraction
1492 const scalar featureCos,
1512 if (diag[0] != -1 && diag[1] != -1)
1516 const label i0 = f[diag[0]];
1519 const label i1 = f[diag[1]];
1522 const point mid = 0.5*(pt0+pt1);
1524 const scalar cosAngle =
mag 1526 patchConstraints[i0].second()
1527 & patchConstraints[i1].second()
1536 if (cosAngle > featureCos)
1541 scalar minDistSqr = GREAT;
1544 label pointi = f[fp];
1545 if (patchConstraints[pointi].first() <= 1)
1548 scalar distSqr =
magSqr(mid-pt);
1549 if (distSqr < minDistSqr)
1551 distSqr = minDistSqr;
1558 label minPointi = f[minFp];
1559 patchAttraction[minPointi] =
1561 patchConstraints[minPointi] =
1562 patchConstraints[f[diag[0]]];
1588 Foam::snappySnapDriver::findNearFeatureEdge
1590 const bool isRegionEdge,
1595 const point& estimatedPt,
1633 label featI = nearEdgeFeat[0];
1644 edgeConstraints[featI][nearInfo.
index()].append(c);
1647 patchAttraction[pointi] =
1649 patchConstraints[pointi] =
c;
1656 Foam::snappySnapDriver::findNearFeaturePoint
1658 const bool isRegionPoint,
1663 const point& estimatedPt,
1688 label featI = nearFeat[0];
1694 label featPointi = nearInfo[0].index();
1695 const point& featPt = nearInfo[0].hitPoint();
1696 scalar distSqr =
magSqr(featPt-pt);
1699 label oldPointi = pointAttractor[featI][featPointi];
1701 if (oldPointi != -1)
1713 pointAttractor[featI][featPointi] = pointi;
1718 patchAttraction[pointi] = featPt-pt;
1719 patchConstraints[pointi] =
1723 patchAttraction[oldPointi] =
Zero;
1745 pointAttractor[featI][featPointi] = pointi;
1750 patchAttraction[pointi] = featPt-pt;
1759 void Foam::snappySnapDriver::determineFeatures
1762 const scalar featureCos,
1763 const bool multiRegionFeatureSnap,
1791 featureEdgeStr.
reset 1795 meshRefiner_.mesh().time().path()
1796 /
"featureEdge_" +
name(iter) +
".obj" 1799 Info<<
"Dumping feature-edge sampling to " 1800 << featureEdgeStr().name() <<
endl;
1806 meshRefiner_.mesh().time().path()
1807 /
"missedFeatureEdge_" +
name(iter) +
".obj" 1810 Info<<
"Dumping feature-edges that are too far away to " 1811 << missedEdgeStr().name() <<
endl;
1813 featurePointStr.
reset 1817 meshRefiner_.mesh().time().path()
1818 /
"featurePoint_" +
name(iter) +
".obj" 1821 Info<<
"Dumping feature-point sampling to " 1822 << featurePointStr().name() <<
endl;
1837 featureAttractionUsingReconstruction
1848 pointFaceSurfNormals,
1864 (constraint.
first() > patchConstraints[pointi].
first())
1866 (constraint.
first() == patchConstraints[pointi].
first())
1867 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
1871 patchAttraction[pointi] = attraction;
1872 patchConstraints[pointi] = constraint;
1875 if (patchConstraints[pointi].first() == 1)
1878 if (multiRegionFeatureSnap)
1880 const point estimatedPt(pt + nearestDisp[pointi]);
1886 pointFacePatchID[pointi],
1892 if (multiPatchPt.
hit())
1917 if (featureEdgeStr.
valid())
1919 featureEdgeStr().write
1927 if (missedEdgeStr.
valid())
1929 missedEdgeStr().write
1938 else if (patchConstraints[pointi].first() == 2)
1943 const point estimatedPt(pt + patchAttraction[pointi]);
1948 bool hasSnapped =
false;
1949 if (multiRegionFeatureSnap)
1956 pointFacePatchID[pointi],
1961 if (multiPatchPt.
hit())
1963 if (multiPatchPt.
index() == 0)
1966 nearInfo = findNearFeatureEdge
1986 nearInfo = findNearFeaturePoint
2015 nearInfo = findNearFeatureEdge
2038 patchConstraints[pointi].first() == 3
2039 && featurePointStr.
valid()
2042 featurePointStr().write
2049 patchConstraints[pointi].first() == 2
2050 && featureEdgeStr.
valid()
2053 featureEdgeStr().write
2061 if (missedEdgeStr.
valid())
2063 missedEdgeStr().write
2070 else if (patchConstraints[pointi].first() == 3)
2073 const point estimatedPt(pt + patchAttraction[pointi]);
2077 if (multiRegionFeatureSnap)
2084 pointFacePatchID[pointi],
2089 if (multiPatchPt.
hit())
2092 nearInfo = findNearFeaturePoint
2113 nearInfo = findNearFeaturePoint
2136 nearInfo = findNearFeaturePoint
2157 if (info.
hit() && featurePointStr.
valid())
2159 featurePointStr().write
2170 void Foam::snappySnapDriver::determineBaffleFeatures
2173 const scalar featureCos,
2213 label facei = eFaces[i];
2248 meshRefiner_.mesh().time().path()
2249 /
"baffleEdge_" +
name(iter) +
".obj" 2252 Info<<
nl <<
"Dumping baffle-edges to " 2253 << baffleEdgeStr().name() <<
endl;
2259 label nBaffleEdges = 0;
2266 forAll(edgeFaceNormals, edgeI)
2270 if (efn.
size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2272 isBaffleEdge[edgeI] =
true;
2275 pointStatus[e[0]] = 0;
2276 pointStatus[e[1]] = 0;
2278 if (baffleEdgeStr.
valid())
2289 <<
" baffle edges out of " 2291 <<
" edges." <<
endl;
2315 forAll(pointStatus, pointi)
2319 if (pointStatus[pointi] == 0)
2335 if (!nearInfo.
second().hit())
2341 else if (pointStatus[pointi] == 1)
2353 label featI = nearFeat[0];
2357 label featPointi = nearInfo[0].index();
2358 const point& featPt = nearInfo[0].hitPoint();
2359 scalar distSqr =
magSqr(featPt-pt);
2362 label oldPointi = pointAttractor[featI][featPointi];
2373 pointAttractor[featI][featPointi] = pointi;
2379 patchAttraction[pointi] = featPt-pt;
2380 patchConstraints[pointi] =
2383 if (oldPointi != -1)
2438 void Foam::snappySnapDriver::reverseAttractMeshPoints
2475 bb = bb.extend(rndGen, 1
e-4);
2476 bb.min() -=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2477 bb.max() +=
point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
2487 forAll(rawPatchConstraints, pointi)
2489 if (rawPatchConstraints[pointi].first() >= 2)
2491 isFeatureEdgeOrPoint[pointi] =
true;
2498 <<
" for reverse attraction." <<
endl;
2506 isFeatureEdgeOrPoint,
2511 for (
label nGrow = 0; nGrow < 1; nGrow++)
2513 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
2521 if (isFeatureEdgeOrPoint[f[fp]])
2526 newIsFeatureEdgeOrPoint[f[fp]] =
true;
2533 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
2539 isFeatureEdgeOrPoint,
2547 forAll(isFeatureEdgeOrPoint, pointi)
2549 if (isFeatureEdgeOrPoint[pointi])
2551 attractPoints.append(pointi);
2558 <<
" for reverse attraction." <<
endl;
2573 patchAttraction =
Zero;
2577 forAll(edgeAttractors, featI)
2581 edgeConstraints[featI];
2583 forAll(edgeAttr, featEdgeI)
2589 const point& featPt = attr[i];
2606 patchConstraints[pointi].first() <= 1
2607 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
2610 patchAttraction[pointi] = attraction;
2611 patchConstraints[pointi] = edgeConstr[featEdgeI][i];
2617 <<
"Did not find pp point near " << featPt
2634 forAll(pointAttractor, featI)
2636 const labelList& pointAttr = pointAttractor[featI];
2639 forAll(pointAttr, featPointi)
2641 if (pointAttr[featPointi] != -1)
2643 const point& featPt = features[featI].points()
2661 const point attraction = featPt-pt;
2666 if (patchConstraints[pointi].first() <= 1)
2668 patchAttraction[pointi] = attraction;
2669 patchConstraints[pointi] = pointConstr[featPointi];
2671 else if (patchConstraints[pointi].first() == 2)
2673 patchAttraction[pointi] = attraction;
2674 patchConstraints[pointi] = pointConstr[featPointi];
2676 else if (patchConstraints[pointi].first() == 3)
2682 <
magSqr(patchAttraction[pointi])
2685 patchAttraction[pointi] = attraction;
2686 patchConstraints[pointi] =
2687 pointConstr[featPointi];
2697 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
2700 const bool avoidSnapProblems,
2701 const scalar featureCos,
2702 const bool multiRegionFeatureSnap,
2731 label nFeatEdges = features[featI].edges().
size();
2732 edgeAttractors[featI].setSize(nFeatEdges);
2733 edgeConstraints[featI].setSize(nFeatEdges);
2743 label nFeatPoints = features[featI].points().
size();
2744 pointAttractor[featI].
setSize(nFeatPoints, -1);
2756 multiRegionFeatureSnap,
2762 pointFaceSurfNormals,
2791 determineBaffleFeatures
2816 reverseAttractMeshPoints
2832 rawPatchConstraints,
2845 meshRefiner_.mesh().time().path()
2846 /
"edgeAttractors_" +
name(iter) +
".obj" 2848 Info<<
"Dumping feature-edge attraction to " 2849 << featureEdgeStr.name() <<
endl;
2853 meshRefiner_.mesh().time().path()
2854 /
"pointAttractors_" +
name(iter) +
".obj" 2856 Info<<
"Dumping feature-point attraction to " 2857 << featurePointStr.name() <<
endl;
2859 forAll(patchConstraints, pointi)
2862 const vector& attr = patchAttraction[pointi];
2864 if (patchConstraints[pointi].first() == 2)
2868 else if (patchConstraints[pointi].first() == 3)
2876 if (avoidSnapProblems)
2881 if (multiRegionFeatureSnap)
2883 releasePointsNextToMultiPatch
2895 rawPatchConstraints,
2917 rawPatchConstraints,
2928 avoidDiagonalAttraction
2938 if (debug&meshRefinement::ATTRACTION)
2942 meshRefiner_.mesh().time().path()
2943 /
"patchAttraction_" +
name(iter) +
".obj",
2951 void Foam::snappySnapDriver::preventFaceSqueeze
2954 const scalar featureCos,
2978 label nConstraints = 0;
2981 label pointi = f[fp];
2984 if (patchConstraints[pointi].first() > 1)
2986 points[fp] = pt + patchAttraction[pointi];
2995 if (nConstraints == f.
size())
2998 scalar newArea = singleF.
mag(points);
2999 if (newArea < 0.1*oldArea)
3006 scalar
s =
magSqr(patchAttraction[f[fp]]);
3015 label pointi = f[maxFp];
3017 patchAttraction[pointi] *= 0.5;
3028 const bool avoidSnapProblems,
3030 const scalar featureCos,
3031 const scalar featureAttract,
3043 Info<<
"Overriding displacement on features :" <<
nl 3044 <<
" implicit features : " << implicitFeatureAttraction <<
nl 3045 <<
" explicit features : " << explicitFeatureAttraction <<
nl 3046 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl 3074 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[f[fp]]);
3086 labelList faceSurfaceGlobalRegion(pp.size(), -1);
3096 faceSurfaceGlobalRegion,
3106 calcNearestFacePointProperties
3113 faceSurfaceGlobalRegion,
3115 pointFaceSurfNormals,
3135 patchAttraction =
Zero;
3140 if (implicitFeatureAttraction)
3145 featureAttractionUsingReconstruction
3155 pointFaceSurfNormals,
3165 if (explicitFeatureAttraction)
3173 featureAttractionUsingFeatureEdges
3178 multiRegionFeatureSnap,
3184 pointFaceSurfNormals,
3227 Info<<
"Attraction:" << endl
3229 <<
" avg:" << avgPatchDisp << endl
3230 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
3231 <<
" avg:" << avgPatchAttr <<
endl;
3245 forAll(patchConstraints, pointi)
3247 if (patchConstraints[pointi].first() > 1)
3250 (1.0-featureAttract)*patchDisp[pointi]
3251 + featureAttract*patchAttraction[pointi];
3259 label nMasterPoints = 0;
3264 forAll(patchConstraints, pointi)
3266 if (isPatchMasterPoint[pointi])
3270 if (patchConstraints[pointi].first() == 1)
3274 else if (patchConstraints[pointi].first() == 2)
3278 else if (patchConstraints[pointi].first() == 3)
3289 Info<<
"Feature analysis : total master points:" 3291 <<
" attraction to :" <<
nl 3292 <<
" feature point : " << nPoint <<
nl 3293 <<
" feature edge : " << nEdge <<
nl 3294 <<
" nearest surface : " << nPlanar <<
nl 3295 <<
" rest : " << nMasterPoints-nPoint-nEdge-nPlanar
3312 if (featureAttract < 1-0.001)
3352 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
3363 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
3370 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
3377 (1.0-featureAttract)*smoothedPatchDisp
3378 + featureAttract*tangPatchDisp;
3382 const scalar
relax = featureAttract;
3394 vector(GREAT, GREAT, GREAT)
const labelListList & pointFaces() const
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
cachedRandom rndGen(label(0),-1)
const labelList & patchID() const
Per boundary face label the patch index.
label nPoints() const
Return number of points supporting patch faces.
Switch multiRegionFeatureSnap() const
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Holds (reference to) pointField. Encapsulation of data needed for octree searches. Used for searching for nearest point. No bounding boxes around points. Only overlaps and calcNearest are implemented, rest makes little sense.
A direction and a reference point.
errorManipArg< error, int > exit(error &err, const int errNo=1)
const Type1 & first() const
Return first.
A face is a list of labels corresponding to mesh vertices.
const double e
Elementary charge.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 2-tuple for storing two objects of different types.
bool empty() const
Return true if the UList is empty (ie, size() is zero)
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Switch explicitFeatureSnap() const
bool hit() const
Is there a hit.
const vectorField & faceAreas() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
const labelList & boundaryPoints() const
Return list of boundary points,.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
label otherVertex(const label a) const
Given one vertex, return the other.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
PointIndexHit< point > pointIndexHit
Vector< scalar > vector
A scalar version of the templated Vector.
const indirectPrimitivePatch & patch() const
Reference to patch.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
const vectorField & faceCentres() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
T & first()
Return the first element of the list.
void findNearestRegion(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &hitSurface, labelList &hitRegion) const
Find nearest point on surfaces. Return surface and region on.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
virtual const pointField & points() const
Return raw points.
const labelListList & pointEdges() const
Return point-edge addressing.
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
void operator()(List< T > &x, const List< T > &y) const
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Application of (multi-)patch point contraints.
Type gMaxMagSqr(const UList< Type > &f, const label comm)
Encapsulates queries for features.
void applyConstraint(const vector &cd)
Apply and accumulate the effect of the given constraint direction.
A list of faces which address into the list of points.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
label findZoneID(const word &zoneName) const
Find zone index given a name.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
vectorField pointField
pointField is a vectorField.
ray planeIntersect(const plane &) const
Return the cutting line between this plane and another.
void reset(T *=0)
If object pointer already set, delete object and set to given.
An ordered pair of two objects of type <T> with first() and second() elements.
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
label start() const
Return start label of this patch in the polyMesh face list.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
A class for handling words, derived from string.
static PackedBoolList getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Simple container to keep together snap specific information.
void append(const T &)
Append an element at the end of the list.
void findNearestPoint(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo) const
Find nearest feature point. Sets:
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Pair< label > labelPair
Label pair.
static const label labelMax
List< label > labelList
A List of labels.
const labelListList & pointEdges() const
Accumulates point constraints through successive applications of the applyConstraint function...
point planePlaneIntersect(const plane &, const plane &) const
Return the cutting point between this plane and two other planes.
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Simple random number generator.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
OFstream which keeps track of vertices.
const Field< PointType > & faceCentres() const
Return face centres for patch.
const Type2 & second() const
Return second.
const Point & missPoint() const
Return miss point.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
label nEdges() const
Return number of edges in patch.
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
void setSize(const label)
Reset size of List.
const PtrList< surfaceZonesInfo > & surfZones() const
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
const Field< PointType > & faceNormals() const
Return face normals for patch.
vector point
Point is a vector.
Non-pointer based hierarchical recursive searching.
#define WarningInFunction
Report a warning using Foam::Warning.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
const Type & shapes() const
Reference to shape.
Mesh data needed to do the Finite Volume discretisation.
label index() const
Return index.
A List with indirect addressing.
const dimensionedScalar c
Speed of light in a vacuum.
const Point & hitPoint() const
Return hit point.
const labelList & surfaces() const
Standard boundBox + extra functionality for use in octree.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
dimensioned< scalar > mag(const dimensioned< Type > &)
const faceZoneMesh & faceZones() const
Return face zone mesh.
const labelListList & edgeFaces() const
Return edge-face addressing.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
static const Vector< scalar > zero
A subset of mesh faces organised as a primitive patch.
const labelListList & pointFaces() const
Return point-face addressing.
A patch is a list of labels that address the faces in the global face list.
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
void clear()
Clear the addressed list, i.e. set the size to zero.
virtual const faceList & faces() const
Return raw faces.
fileName path() const
Return path.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
label nInternalFaces() const
A List with indirect addressing.
scalar mag(const pointField &) const
Magnitude of face area.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
void findNearestEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest feature edge. Sets:
void findNearestRegionEdge(const pointField &samples, const scalarField &nearestDistSqr, labelList &nearFeature, List< pointIndexHit > &nearInfo, vectorField &nearNormal) const
Find nearest point on nearest region edge. Sets:
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
Switch implicitFeatureSnap() const
dimensionSet transform(const dimensionSet &)
const Time & time() const
Return the top-level database.
label size() const
Return the number of elements in the UPtrList.