59 x.setSize(sz+
y.size());
71 bool Foam::snappySnapDriver::isFeaturePoint
73 const scalar featureCos,
75 const bitSet& isFeatureEdge,
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());
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].point() - 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());
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].point() - 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());
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];
543 pFid[nFaces] = globalToMasterPatch_[globalRegioni];
564 const polyPatch&
pp =
pbm[patchi];
566 if (
pp.coupled() || isA<emptyPolyPatch>(
pp))
570 label meshFacei =
pp.start()+i;
579 label meshFacei =
pp.addressing()[i];
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]);
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];
665 pointFaceSurfNormals,
666 listPlusEqOp<point>(),
675 listPlusEqOp<point>(),
684 forAll(pointFaceCentres, pointi)
688 List<point>& pFc = pointFaceCentres[pointi];
699 listPlusEqOp<point>(),
703 forAll(pointFaceCentres, 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
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,
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)
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);
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);
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;
1256 if (feStr && patchConstraints[pointi].first() == 2)
1258 feStr().writeLine(pt, pt+patchAttraction[pointi]);
1260 else if (fpStr && patchConstraints[pointi].first() == 3)
1262 fpStr().writeLine(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)
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)
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 =
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)
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];
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();
1583 reduce(nChanged, sumOp<label>());
1584 Info<<
"Removing constraints near multi-patch points : changed " 1585 << nChanged <<
" points" <<
endl;
1593 const List<pointConstraint>& patchConstraints,
1605 for (label startFp = 0; startFp <
f.
size()-2; 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
1742 if (localF.size() >= 4)
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()));
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
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);
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
2049 if (
diag[0] != -1 &&
diag[1] != -1)
2053 const label i0 =
f[
diag[0]];
2056 const label i1 =
f[
diag[1]];
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)
2086 if (distSqr < minDistSqr)
2094 label minPointi =
f[minFp];
2095 patchAttraction[minPointi] =
2097 patchConstraints[minPointi] = patchConstraints[
f[
diag[0]]];
2123 Foam::snappySnapDriver::findNearFeatureEdge
2125 const bool isRegionEdge,
2130 const point& estimatedPt,
2132 List<List<DynamicList<point>>>& edgeAttractors,
2133 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2135 List<pointConstraint>& patchConstraints
2138 const refinementFeatures& features = meshRefiner_.features();
2141 List<pointIndexHit> nearEdgeInfo;
2146 features.findNearestRegionEdge
2157 features.findNearestEdge
2168 label feati = nearEdgeFeat[0];
2174 edgeAttractors[feati][nearInfo.index()].append(nearInfo.point());
2175 pointConstraint
c(Tuple2<label, vector>(2, nearNormal[0]));
2176 edgeConstraints[feati][nearInfo.index()].append(
c);
2179 patchAttraction[pointi] = nearInfo.point()-
pp.
localPoints()[pointi];
2180 patchConstraints[pointi] =
c;
2182 return Tuple2<label, pointIndexHit>(feati, nearInfo);
2187 Foam::snappySnapDriver::findNearFeaturePoint
2189 const bool isRegionPoint,
2194 const point& estimatedPt,
2197 List<labelList>& pointAttractor,
2198 List<List<pointConstraint>>& pointConstraints,
2200 List<List<DynamicList<point>>>& edgeAttractors,
2201 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2204 List<pointConstraint>& patchConstraints
2207 const refinementFeatures& features = meshRefiner_.features();
2212 List<pointIndexHit> nearInfo;
2213 features.findNearestPoint
2221 label feati = nearFeat[0];
2227 label featPointi = nearInfo[0].index();
2228 const point& featPt = nearInfo[0].hitPoint();
2229 scalar distSqr = featPt.distSqr(pt);
2232 label oldPointi = pointAttractor[feati][featPointi];
2234 if (oldPointi != -1)
2246 pointAttractor[feati][featPointi] = pointi;
2247 pointConstraints[feati][featPointi].first() = 3;
2248 pointConstraints[feati][featPointi].second() =
Zero;
2251 patchAttraction[pointi] = featPt-pt;
2252 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2255 patchAttraction[oldPointi] =
Zero;
2256 patchConstraints[oldPointi] = pointConstraint();
2277 pointAttractor[feati][featPointi] = pointi;
2278 pointConstraints[feati][featPointi].first() = 3;
2279 pointConstraints[feati][featPointi].second() =
Zero;
2282 patchAttraction[pointi] = featPt-pt;
2283 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2287 return Tuple2<label, pointIndexHit>(feati, nearInfo[0]);
2293 void Foam::snappySnapDriver::determineFeatures
2296 const scalar featureCos,
2297 const bool multiRegionFeatureSnap,
2303 const List<List<point>>& pointFaceSurfNormals,
2304 const List<List<point>>& pointFaceDisp,
2305 const List<List<point>>& pointFaceCentres,
2309 List<labelList>& pointAttractor,
2310 List<List<pointConstraint>>& pointConstraints,
2312 List<List<DynamicList<point>>>& edgeAttractors,
2313 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2316 List<pointConstraint>& patchConstraints
2319 autoPtr<OBJstream> featureEdgeStr;
2320 autoPtr<OBJstream> missedEdgeStr;
2321 autoPtr<OBJstream> featurePointStr;
2322 autoPtr<OBJstream> missedMP0Str;
2323 autoPtr<OBJstream> missedMP1Str;
2327 featureEdgeStr.reset
2331 meshRefiner_.mesh().time().path()
2332 /
"featureEdge_" +
name(iter) +
".obj" 2335 Info<<
"Dumping feature-edge sampling to " 2336 << featureEdgeStr().name() <<
endl;
2342 meshRefiner_.mesh().time().path()
2343 /
"missedFeatureEdge_" +
name(iter) +
".obj" 2346 Info<<
"Dumping feature-edges that are too far away to " 2347 << missedEdgeStr().name() <<
endl;
2349 featurePointStr.reset
2353 meshRefiner_.mesh().time().path()
2354 /
"featurePoint_" +
name(iter) +
".obj" 2357 Info<<
"Dumping feature-point sampling to " 2358 << featurePointStr().name() <<
endl;
2364 meshRefiner_.mesh().time().path()
2365 /
"missedFeatureEdgeFromMPEdge_" +
name(iter) +
".obj" 2368 Info<<
"Dumping region-edges that are too far away to " 2369 << missedMP0Str().name() <<
endl;
2375 meshRefiner_.mesh().time().path()
2376 /
"missedFeatureEdgeFromMPPoint_" +
name(iter) +
".obj" 2379 Info<<
"Dumping region-points that are too far away to " 2380 << missedMP1Str().name() <<
endl;
2384 DynamicList<point> surfacePoints(4);
2385 DynamicList<vector> surfaceNormals(4);
2402 pointConstraint constraint;
2404 featureAttractionUsingReconstruction
2415 pointFaceSurfNormals,
2460 (constraint.first() > patchConstraints[pointi].first())
2462 (constraint.first() == patchConstraints[pointi].first())
2463 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
2467 patchAttraction[pointi] = attraction;
2468 patchConstraints[pointi] = constraint;
2471 if (patchConstraints[pointi].first() == 1)
2474 if (multiRegionFeatureSnap)
2476 const point estimatedPt(pt + nearestDisp[pointi]);
2482 pointFacePatchID[pointi],
2488 if (multiPatchPt.hit())
2493 Tuple2<label, pointIndexHit> nearInfo =
2500 multiPatchPt.point(),
2515 featureEdgeStr().writeLine
2526 missedEdgeStr().writeLine
2529 multiPatchPt.point()
2536 else if (patchConstraints[pointi].first() == 2)
2541 const point estimatedPt(pt + patchAttraction[pointi]);
2546 bool hasSnapped =
false;
2547 if (multiRegionFeatureSnap)
2554 pointFacePatchID[pointi],
2559 if (multiPatchPt.hit())
2561 if (multiPatchPt.index() == 0)
2564 nearInfo = findNearFeatureEdge
2581 if (missedMP0Str && !nearInfo.second().hit())
2583 missedMP0Str().writeLine(pt, estimatedPt);
2594 nearInfo = findNearFeaturePoint
2623 if (!nearInfo.second().hit())
2625 nearInfo = findNearFeatureEdge
2643 if (missedMP1Str && !nearInfo.second().hit())
2645 missedMP1Str().writeLine(pt, estimatedPt);
2657 nearInfo = findNearFeatureEdge
2681 && patchConstraints[pointi].first() == 3
2684 featurePointStr().writeLine(pt, info.point());
2689 && patchConstraints[pointi].first() == 2
2692 featureEdgeStr().writeLine(pt, info.point());
2699 missedEdgeStr().writeLine(pt, estimatedPt);
2703 else if (patchConstraints[pointi].first() == 3)
2706 const point estimatedPt(pt + patchAttraction[pointi]);
2710 if (multiRegionFeatureSnap)
2717 pointFacePatchID[pointi],
2722 if (multiPatchPt.hit())
2725 nearInfo = findNearFeaturePoint
2746 nearInfo = findNearFeaturePoint
2769 nearInfo = findNearFeaturePoint
2790 if (featurePointStr && info.hit())
2792 featurePointStr().writeLine(pt, info.point());
2811 void Foam::snappySnapDriver::determineBaffleFeatures
2814 const bool baffleFeaturePoints,
2815 const scalar featureCos,
2821 List<labelList>& pointAttractor,
2822 List<List<pointConstraint>>& pointConstraints,
2824 List<List<DynamicList<point>>>& edgeAttractors,
2825 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2828 List<pointConstraint>& patchConstraints
2831 const fvMesh&
mesh = meshRefiner_.mesh();
2832 const refinementFeatures& features = meshRefiner_.features();
2835 List<List<point>> edgeFaceNormals(
pp.
nEdges());
2841 List<point>& eFc = edgeFaceNormals[edgei];
2845 label facei = eFaces[i];
2862 listPlusEqOp<point>(),
2875 autoPtr<OBJstream> baffleEdgeStr;
2882 meshRefiner_.mesh().time().path()
2883 /
"baffleEdge_" +
name(iter) +
".obj" 2886 Info<<
nl <<
"Dumping baffle-edges to " 2887 << baffleEdgeStr().name() <<
endl;
2893 label nBaffleEdges = 0;
2900 forAll(edgeFaceNormals, edgei)
2902 const List<point>& efn = edgeFaceNormals[edgei];
2904 if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2906 isBaffleEdge.set(edgei);
2909 pointStatus[
e[0]] = 0;
2910 pointStatus[
e[1]] = 0;
2919 reduce(nBaffleEdges, sumOp<label>());
2921 Info<<
"Detected " << nBaffleEdges
2922 <<
" baffle edges out of " 2924 <<
" edges." <<
endl;
2949 label nBafflePoints = 0;
2950 forAll(pointStatus, pointi)
2952 if (pointStatus[pointi] != -1)
2957 reduce(nBafflePoints, sumOp<label>());
2960 label nPointAttract = 0;
2961 label nEdgeAttract = 0;
2963 forAll(pointStatus, pointi)
2967 if (pointStatus[pointi] == 0)
2971 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
2990 if (nearInfo.second().hit())
2994 if (baffleFeaturePoints)
2996 nearInfo = findNearFeaturePoint
3016 if (nearInfo.first() != -1)
3024 else if (pointStatus[pointi] == 1)
3027 List<pointIndexHit> nearInfo;
3028 features.findNearestPoint
3036 label feati = nearFeat[0];
3042 label featPointi = nearInfo[0].index();
3043 const point& featPt = nearInfo[0].hitPoint();
3044 scalar distSqr = featPt.distSqr(pt);
3047 label oldPointi = pointAttractor[feati][featPointi];
3058 pointAttractor[feati][featPointi] = pointi;
3059 pointConstraints[feati][featPointi].first() = 3;
3060 pointConstraints[feati][featPointi].second() =
Zero;
3063 patchAttraction[pointi] = featPt-pt;
3064 patchConstraints[pointi] =
3065 pointConstraints[feati][featPointi];
3067 if (oldPointi != -1)
3103 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3117 if (nearInfo.first() != -1)
3125 reduce(nPointAttract, sumOp<label>());
3126 reduce(nEdgeAttract, sumOp<label>());
3128 Info<<
"Baffle points : " << nBafflePoints
3129 <<
" of which attracted to :" <<
nl 3130 <<
" feature point : " << nPointAttract <<
nl 3131 <<
" feature edge : " << nEdgeAttract <<
nl 3132 <<
" rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3138 void Foam::snappySnapDriver::reverseAttractMeshPoints
3146 const List<labelList>& pointAttractor,
3147 const List<List<pointConstraint>>& pointConstraints,
3149 const List<List<DynamicList<point>>>& edgeAttractors,
3150 const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3153 const List<pointConstraint>& rawPatchConstraints,
3157 List<pointConstraint>& patchConstraints
3160 const refinementFeatures& features = meshRefiner_.features();
3175 bb.inflate(
rndGen, 1
e-4, ROOTVSMALL);
3179 DynamicList<label> attractPoints(
pp.
nPoints());
3181 const fvMesh&
mesh = meshRefiner_.mesh();
3185 forAll(rawPatchConstraints, pointi)
3187 if (rawPatchConstraints[pointi].first() >= 2)
3189 isFeatureEdgeOrPoint[pointi] =
true;
3195 <<
" mesh points out of " 3197 <<
" for reverse attraction." <<
endl;
3205 isFeatureEdgeOrPoint,
3210 for (label nGrow = 0; nGrow < 1; nGrow++)
3212 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3220 if (isFeatureEdgeOrPoint[
f[fp]])
3225 newIsFeatureEdgeOrPoint[
f[fp]] =
true;
3232 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3238 isFeatureEdgeOrPoint,
3246 forAll(isFeatureEdgeOrPoint, pointi)
3248 if (isFeatureEdgeOrPoint[pointi])
3250 attractPoints.append(pointi);
3256 <<
" mesh points out of " 3258 <<
" for reverse attraction." <<
endl;
3262 indexedOctree<treeDataPoint> ppTree
3273 patchAttraction =
Zero;
3275 patchConstraints = pointConstraint();
3277 forAll(edgeAttractors, feati)
3279 const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3280 const List<DynamicList<pointConstraint>>& edgeConstr =
3281 edgeConstraints[feati];
3283 forAll(edgeAttr, featEdgei)
3285 const DynamicList<point>& attr = edgeAttr[featEdgei];
3289 const point& featPt = attr[i];
3299 ppTree.shapes().objectIndex(nearInfo.index());
3301 const point attraction
3304 - ppTree.shapes()[nearInfo.index()]
3311 patchConstraints[pointi].first() <= 1
3312 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
3315 patchAttraction[pointi] = attraction;
3316 patchConstraints[pointi] = edgeConstr[featEdgei][i];
3321 static label nWarn = 0;
3326 <<
"Did not find pp point near " << featPt
3332 <<
"Reached warning limit " << nWarn
3333 <<
". Suppressing further warnings." <<
endl;
3352 forAll(pointAttractor, feati)
3354 const labelList& pointAttr = pointAttractor[feati];
3355 const List<pointConstraint>& pointConstr = pointConstraints[feati];
3357 forAll(pointAttr, featPointi)
3359 if (pointAttr[featPointi] != -1)
3361 const point& featPt = features[feati].points()
3376 ppTree.shapes().objectIndex(nearInfo.index());
3378 const point attraction
3381 - ppTree.shapes()[nearInfo.index()]
3387 if (patchConstraints[pointi].first() <= 1)
3389 patchAttraction[pointi] = attraction;
3390 patchConstraints[pointi] = pointConstr[featPointi];
3392 else if (patchConstraints[pointi].first() == 2)
3394 patchAttraction[pointi] = attraction;
3395 patchConstraints[pointi] = pointConstr[featPointi];
3397 else if (patchConstraints[pointi].first() == 3)
3403 <
magSqr(patchAttraction[pointi])
3406 patchAttraction[pointi] = attraction;
3407 patchConstraints[pointi] =
3408 pointConstr[featPointi];
3418 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3421 const bool multiRegionFeatureSnap,
3423 const bool detectBaffles,
3424 const bool baffleFeaturePoints,
3426 const bool releasePoints,
3427 const bool stringFeatures,
3428 const bool avoidDiagonal,
3430 const scalar featureCos,
3437 const List<List<point>>& pointFaceSurfNormals,
3438 const List<List<point>>& pointFaceDisp,
3439 const List<List<point>>& pointFaceCentres,
3443 List<pointConstraint>& patchConstraints
3446 const refinementFeatures& features = meshRefiner_.features();
3447 const fvMesh&
mesh = meshRefiner_.mesh();
3449 const bitSet isPatchMasterPoint
3464 List<List<DynamicList<point>>> edgeAttractors(features.size());
3465 List<List<DynamicList<pointConstraint>>> edgeConstraints
3471 label nFeatEdges = features[feati].edges().size();
3472 edgeAttractors[feati].setSize(nFeatEdges);
3473 edgeConstraints[feati].setSize(nFeatEdges);
3479 List<labelList> pointAttractor(features.size());
3480 List<List<pointConstraint>> pointConstraints(features.size());
3483 label nFeatPoints = features[feati].points().size();
3484 pointAttractor[feati].setSize(nFeatPoints, -1);
3485 pointConstraints[feati].setSize(nFeatPoints);
3490 List<pointConstraint> rawPatchConstraints(
pp.
nPoints());
3496 multiRegionFeatureSnap,
3502 pointFaceSurfNormals,
3521 Info<<
"Raw geometric feature analysis : ";
3522 writeStats(
pp, isPatchMasterPoint, rawPatchConstraints);
3538 determineBaffleFeatures
3541 baffleFeaturePoints,
3562 Info<<
"After baffle feature analysis : ";
3563 writeStats(
pp, isPatchMasterPoint, rawPatchConstraints);
3571 reverseAttractMeshPoints
3587 rawPatchConstraints,
3597 Info<<
"Reverse attract feature analysis : ";
3604 OBJstream featureEdgeStr
3606 meshRefiner_.mesh().time().path()
3607 /
"edgeAttractors_" +
name(iter) +
".obj" 3609 Info<<
"Dumping feature-edge attraction to " 3610 << featureEdgeStr.name() <<
endl;
3612 OBJstream featurePointStr
3614 meshRefiner_.mesh().time().path()
3615 /
"pointAttractors_" +
name(iter) +
".obj" 3617 Info<<
"Dumping feature-point attraction to " 3618 << featurePointStr.name() <<
endl;
3620 forAll(patchConstraints, pointi)
3623 const vector& attr = patchAttraction[pointi];
3625 if (patchConstraints[pointi].first() == 2)
3627 featureEdgeStr.writeLine(pt, pt+attr);
3629 else if (patchConstraints[pointi].first() == 3)
3631 featurePointStr.writeLine(pt, pt+attr);
3643 releasePointsNextToMultiPatch
3655 rawPatchConstraints,
3678 rawPatchConstraints,
3691 avoidDiagonalAttraction
3706 meshRefiner_.mesh().time().path()
3707 /
"patchAttraction_" +
name(iter) +
".obj",
3716 void Foam::snappySnapDriver::preventFaceSqueeze
3719 const scalar featureCos,
3726 List<pointConstraint>& patchConstraints
3729 autoPtr<OBJstream> strPtr;
3736 meshRefiner_.mesh().time().path()
3737 /
"faceSqueeze_" +
name(iter) +
".obj" 3740 Info<<
"Dumping faceSqueeze corrections to " 3741 << strPtr().name() <<
endl;
3753 singleF.setSize(
f.
size());
3754 for (label i = 0; i <
f.
size(); i++)
3759 label nConstraints = 0;
3762 label pointi =
f[fp];
3765 if (patchConstraints[pointi].first() > 1)
3767 points[fp] = pt + patchAttraction[pointi];
3776 if (nConstraints ==
f.
size())
3790 scalar
s = pt.distSqr(nextPt);
3799 label pointi =
f.prevLabel(maxFp);
3810 patchAttraction[pointi] = nearestAttraction[pointi];
3814 strPtr().writeLine(pt, pt+patchAttraction[pointi]);
3821 scalar newArea = singleF.mag(
points);
3822 if (newArea < 0.1*oldArea)
3829 scalar
s =
magSqr(patchAttraction[
f[fp]]);
3838 label pointi =
f[maxFp];
3840 patchAttraction[pointi] *= 0.5;
3851 const snapParameters& snapParams,
3852 const bool alignMeshEdges,
3854 const scalar featureCos,
3855 const scalar featureAttract,
3859 motionSmoother& meshMover,
3861 List<pointConstraint>& patchConstraints,
3863 DynamicList<label>& splitFaces,
3864 DynamicList<labelPair>& splits
3874 const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3875 const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3876 const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3878 Info<<
"Overriding displacement on features :" <<
nl 3879 <<
" implicit features : " << implicitFeatureAttraction <<
nl 3880 <<
" explicit features : " << explicitFeatureAttraction <<
nl 3881 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl 3886 const fvMesh&
mesh = meshRefiner_.mesh();
3890 const bitSet isPatchMasterPoint
3903 List<List<point>> pointFaceSurfNormals;
3904 List<List<point>> pointFaceDisp;
3905 List<List<point>> pointFaceCentres;
3906 List<labelList> pointFacePatchID;
3917 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[
f[fp]]);
3939 faceSurfaceGlobalRegion
3949 calcNearestFacePointProperties
3956 faceSurfaceGlobalRegion,
3958 pointFaceSurfNormals,
3977 patchAttraction.
setSize(localPoints.size());
3978 patchAttraction =
Zero;
3980 patchConstraints.setSize(localPoints.size());
3981 patchConstraints = pointConstraint();
3983 if (implicitFeatureAttraction)
3988 featureAttractionUsingReconstruction
3997 pointFaceSurfNormals,
4007 if (explicitFeatureAttraction)
4010 bool releasePoints =
false;
4011 bool stringFeatures =
false;
4012 bool avoidDiagonal =
false;
4015 releasePoints = snapParams.releasePoints();
4016 stringFeatures = snapParams.stringFeatures();
4017 avoidDiagonal = snapParams.avoidDiagonal();
4027 featureAttractionUsingFeatureEdges
4030 multiRegionFeatureSnap,
4032 snapParams.detectBaffles(),
4033 snapParams.baffleFeaturePoints(),
4047 pointFaceSurfNormals,
4057 if (!alignMeshEdges)
4061 degToRad(snapParams.concaveAngle())
4063 const scalar minAreaRatio = snapParams.minAreaRatio();
4065 Info<<
"Experimental: introducing face splits to avoid rotating" 4066 <<
" mesh edges. Splitting faces when" <<
nl 4067 <<
indent <<
"- angle not concave by more than " 4068 << snapParams.concaveAngle() <<
" degrees" <<
nl 4069 <<
indent <<
"- resulting triangles of similar area " 4070 <<
" (ratio within " << minAreaRatio <<
")" <<
nl 4091 Info<<
"Diagonal attraction feature correction : ";
4124 <<
" avg:" << avgPatchDisp <<
endl 4125 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
4126 <<
" avg:" << avgPatchAttr <<
endl;
4140 forAll(patchConstraints, pointi)
4142 if (patchConstraints[pointi].first() > 1)
4145 (1.0-featureAttract)*patchDisp[pointi]
4146 + featureAttract*patchAttraction[pointi];
4154 Info<<
"Feature analysis : ";
4170 if (featureAttract < 1-0.001)
4177 const bitSet isPatchMasterEdge
4209 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4220 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4227 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
4234 (1.0-featureAttract)*smoothedPatchDisp
4235 + featureAttract*tangPatchDisp;
4239 const scalar
relax = featureAttract;
4250 minMagSqrEqOp<point>(),
4251 vector(GREAT, GREAT, GREAT)
const labelListList & pointEdges() const
Return point-edge addressing.
const polyBoundaryMesh & pbm
label nPoints() const
Number of points supporting patch faces.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
void size(const label n)
Older name for setAddressableSize.
const labelList patchIDs(pbm.patchSet(polyPatchNames, false, true).sortedToc())
fileName path() const
Return path.
Ostream & indent(Ostream &os)
Indent stream.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const Field< point_type > & localPoints() const
Return pointField of points in patch.
label nPoints() const noexcept
Number of mesh points.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 2-tuple for storing two objects of dissimilar types. The container is similar in purpose to std::pa...
void append(const T &val)
Append an element at the end of the list.
List< edge > edgeList
List of edge.
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
const labelListList & pointEdges() const
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
bool empty() const noexcept
True if List is empty (ie, size() is zero)
UIndirectList< label > labelUIndList
UIndirectList of labels.
T & first()
Access first element of the list, position [0].
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
const labelList & patchID() const
Per boundary face label the patch index.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
static bitSet getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
const Time & time() const
Return the top-level database.
label nFaces() const noexcept
Number of mesh faces.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field...
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
Type gMaxMagSqr(const UList< Type > &f, const label comm)
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
A list of faces which address into the list of points.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
void clear()
Clear the list, i.e. set size to zero.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Field< point_type > & faceCentres() const
Return face centres for patch.
constexpr scalar pi(M_PI)
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
label nInternalFaces() const noexcept
Number of internal faces.
const Field< point_type > & points() const noexcept
Return reference to global points.
virtual const faceList & faces() const
Return raw faces.
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
const labelListList & pointFaces() const
Return point-face addressing.
label nEdges() const
Number of edges in patch.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
Pair< label > labelPair
A pair of labels.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
label otherVertex(const label vertex) const
Given one vertex label, return the other one.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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]
static bool split(const std::string &line, std::string &key, std::string &val)
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
List< word > wordList
List of word.
pointField points0(pointIOField(IOobject("points", mesh.time().constant(), polyMesh::meshSubDir, mesh, IOobject::MUST_READ, IOobject::NO_WRITE, IOobject::NO_REGISTER)))
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using cell addressing.
static bitSet getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
const vectorField & faceAreas() const
const dimensionedScalar c
Speed of light in a vacuum.
const std::string patch
OpenFOAM patch number as a std::string.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
static labelList findDuplicateFaces(const primitiveMesh &, const labelList &)
Helper routine to find baffles (two boundary faces using the.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
static T gAverage(const bitSet &isMasterElem, const UList< T > &values)
Helper: calculate average.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
void operator()(List< T > &x, const List< T > &y) const
List< bool > boolList
A List of bools.
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
static constexpr const zero Zero
Global zero (0)