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 label globalRegioni = faceSurfaceGlobalRegion[facei];
518 if (isMasterFace[
pp.addressing()[facei]] && globalRegioni != -1)
525 List<point>& pNormals = pointFaceSurfNormals[pointi];
526 pNormals.setSize(nFaces);
527 List<point>& pDisp = pointFaceDisp[pointi];
528 pDisp.setSize(nFaces);
529 List<point>& pFc = pointFaceCentres[pointi];
531 labelList& pFid = pointFacePatchID[pointi];
538 label globalRegioni = faceSurfaceGlobalRegion[facei];
540 if (isMasterFace[
pp.addressing()[facei]] && globalRegioni != -1)
542 pNormals[nFaces] = faceSurfaceNormal[facei];
543 pDisp[nFaces] = faceDisp[facei];
545 pFid[nFaces] = globalToMasterPatch_[globalRegioni];
566 const polyPatch&
pp =
pbm[patchi];
568 if (
pp.coupled() || isA<emptyPolyPatch>(
pp))
572 label meshFacei =
pp.start()+i;
581 label meshFacei =
pp.addressing()[i];
605 const edge&
e = edges[edgei];
606 const labelList& eFaces = edgeFaces[edgei];
608 if (eFaces.size() == 1)
611 isBoundaryPoint.
set(
e[0]);
612 isBoundaryPoint.set(
e[1]);
614 else if (eFaces.size() == 2 && bafflePair[eFaces[0]] == eFaces[1])
617 isBoundaryPoint.set(
e[0]);
618 isBoundaryPoint.set(
e[1]);
632 label patchi =
patchID[bFacei];
641 label pointi = meshToPatchPoint[
f[fp]];
643 if (pointi != -1 && isBoundaryPoint.test(pointi))
645 List<point>& pNormals = pointFaceSurfNormals[pointi];
646 List<point>& pDisp = pointFaceDisp[pointi];
647 List<point>& pFc = pointFaceCentres[pointi];
648 labelList& pFid = pointFacePatchID[pointi];
667 pointFaceSurfNormals,
668 listPlusEqOp<point>(),
677 listPlusEqOp<point>(),
686 forAll(pointFaceCentres, pointi)
690 List<point>& pFc = pointFaceCentres[pointi];
701 listPlusEqOp<point>(),
705 forAll(pointFaceCentres, pointi)
709 List<point>& pFc = pointFaceCentres[pointi];
722 listPlusEqOp<label>(),
730 forAll(pointFaceDisp, pointi)
732 List<point>& pNormals = pointFaceSurfNormals[pointi];
733 List<point>& pDisp = pointFaceDisp[pointi];
734 List<point>& pFc = pointFaceCentres[pointi];
735 labelList& pFid = pointFacePatchID[pointi];
739 pNormals = List<point>(pNormals, visitOrder);
740 pDisp = List<point>(pDisp, visitOrder);
741 pFc = List<point>(pFc, visitOrder);
750 void Foam::snappySnapDriver::correctAttraction
752 const DynamicList<point>& surfacePoints,
753 const DynamicList<label>& surfaceCounts,
762 scalar tang = ((pt-edgePt)&edgeNormal);
766 if (order[0] < order[1])
770 vector attractD = surfacePoints[order[0]]-edgePt;
772 scalar tang2 = (attractD&edgeNormal);
774 attractD -= tang2*edgeNormal;
776 scalar magAttractD =
mag(attractD);
777 scalar fraction = magAttractD/(magAttractD+
mag(edgeOffset));
781 + ((1.0-fraction)*tang2 + fraction*tang)*edgeNormal;
782 edgeOffset = linePt-pt;
791 const List<point>& faceCentres
811 Foam::label Foam::snappySnapDriver::findNormal
813 const scalar featureCos,
815 const DynamicList<vector>& surfaceNormals
822 scalar cosAngle = (
n&surfaceNormals[j]);
826 (cosAngle >= featureCos)
827 || (cosAngle < (-1+0.001))
848 const DynamicList<vector>& surfaceNormals,
876 if (surfaceNormals.size() == 1)
884 labelList normalToPatch(surfaceNormals.size(), -1);
885 forAll(faceToNormalBin, i)
887 if (faceToNormalBin[i] != -1)
889 label&
patch = normalToPatch[faceToNormalBin[i]];
895 else if (
patch == -2)
907 forAll(normalToPatch, normali)
909 if (normalToPatch[normali] == -2)
924 void Foam::snappySnapDriver::writeStats
927 const bitSet& isPatchMasterPoint,
928 const List<pointConstraint>& patchConstraints
931 label nMasterPoints = 0;
936 forAll(patchConstraints, pointi)
938 if (isPatchMasterPoint[pointi])
942 if (patchConstraints[pointi].first() == 1)
946 else if (patchConstraints[pointi].first() == 2)
950 else if (patchConstraints[pointi].first() == 3)
957 reduce(nMasterPoints, sumOp<label>());
958 reduce(nPlanar, sumOp<label>());
959 reduce(nEdge, sumOp<label>());
960 reduce(nPoint, sumOp<label>());
961 Info<<
"total master points :" << nMasterPoints
962 <<
" of which attracted to :" <<
nl 963 <<
" feature point : " << nPoint <<
nl 964 <<
" feature edge : " << nEdge <<
nl 965 <<
" nearest surface : " << nPlanar <<
nl 966 <<
" rest : " << nMasterPoints-nPoint-nEdge-nPlanar
972 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
975 const scalar featureCos,
982 const List<List<point>>& pointFaceSurfNormals,
983 const List<List<point>>& pointFaceDisp,
984 const List<List<point>>& pointFaceCentres,
987 DynamicList<point>& surfacePoints,
988 DynamicList<vector>& surfaceNormals,
992 pointConstraint& patchConstraint
995 patchAttraction =
Zero;
996 patchConstraint = pointConstraint();
998 const List<point>& pfSurfNormals = pointFaceSurfNormals[pointi];
999 const List<point>& pfDisp = pointFaceDisp[pointi];
1000 const List<point>& pfCentres = pointFaceCentres[pointi];
1010 surfacePoints.clear();
1011 surfaceNormals.clear();
1014 faceToNormalBin.setSize(pfDisp.size());
1015 faceToNormalBin = -1;
1019 const point& fc = pfCentres[i];
1020 const vector& fSNormal = pfSurfNormals[i];
1021 const vector& fDisp = pfDisp[i];
1024 if (
magSqr(fDisp) <
sqr(snapDist[pointi]) &&
mag(fSNormal) > VSMALL)
1026 const point pt = fc + fDisp;
1029 faceToNormalBin[i] = findNormal
1036 if (faceToNormalBin[i] != -1)
1044 if (surfacePoints.size() <= 1)
1046 surfacePoints.append(pt);
1047 faceToNormalBin[i] = surfaceNormals.size();
1048 surfaceNormals.append(fSNormal);
1050 else if (surfacePoints.size() == 2)
1052 plane pl0(surfacePoints[0], surfaceNormals[0]);
1053 plane pl1(surfacePoints[1], surfaceNormals[1]);
1054 plane::ray r(pl0.planeIntersect(pl1));
1055 vector featureNormal = r.dir() /
mag(r.dir());
1057 if (
mag(fSNormal&featureNormal) >= 0.001)
1060 surfacePoints.append(pt);
1061 faceToNormalBin[i] = surfaceNormals.size();
1062 surfaceNormals.append(fSNormal);
1065 else if (surfacePoints.size() == 3)
1069 plane pl0(surfacePoints[0], surfaceNormals[0]);
1070 plane pl1(surfacePoints[1], surfaceNormals[1]);
1071 plane pl2(surfacePoints[2], surfaceNormals[2]);
1072 point p012(pl0.planePlaneIntersect(pl1, pl2));
1074 plane::ray r(pl0.planeIntersect(pl1));
1075 vector featureNormal = r.dir() /
mag(r.dir());
1076 if (
mag(fSNormal&featureNormal) >= 0.001)
1078 plane pl3(pt, fSNormal);
1079 point p013(pl0.planePlaneIntersect(pl1, pl3));
1081 if (
mag(p012-p013) > snapDist[pointi])
1084 surfacePoints.append(pt);
1085 faceToNormalBin[i] = surfaceNormals.size();
1086 surfaceNormals.append(fSNormal);
1098 if (surfaceNormals.size() == 1)
1102 ((surfacePoints[0]-pt) & surfaceNormals[0])
1111 patchAttraction = d;
1114 patchConstraint.applyConstraint(surfaceNormals[0]);
1116 else if (surfaceNormals.size() == 2)
1118 plane pl0(surfacePoints[0], surfaceNormals[0]);
1119 plane pl1(surfacePoints[1], surfaceNormals[1]);
1120 plane::ray r(pl0.planeIntersect(pl1));
1124 vector d = r.refPoint()-pt;
1133 patchAttraction = d;
1136 patchConstraint.applyConstraint(surfaceNormals[0]);
1137 patchConstraint.applyConstraint(surfaceNormals[1]);
1139 else if (surfaceNormals.size() == 3)
1142 plane pl0(surfacePoints[0], surfaceNormals[0]);
1143 plane pl1(surfacePoints[1], surfaceNormals[1]);
1144 plane pl2(surfacePoints[2], surfaceNormals[2]);
1145 point cornerPt(pl0.planePlaneIntersect(pl1, pl2));
1146 vector d = cornerPt - pt;
1154 patchAttraction = d;
1157 patchConstraint.applyConstraint(surfaceNormals[0]);
1158 patchConstraint.applyConstraint(surfaceNormals[1]);
1159 patchConstraint.applyConstraint(surfaceNormals[2]);
1165 void Foam::snappySnapDriver::featureAttractionUsingReconstruction
1168 const scalar featureCos,
1174 const List<List<point>>& pointFaceSurfNormals,
1175 const List<List<point>>& pointFaceDisp,
1176 const List<List<point>>& pointFaceCentres,
1180 List<pointConstraint>& patchConstraints
1183 autoPtr<OBJstream> feStr;
1184 autoPtr<OBJstream> fpStr;
1191 meshRefiner_.mesh().time().path()
1192 /
"implicitFeatureEdge_" +
name(iter) +
".obj" 1195 Info<<
"Dumping implicit feature-edge direction to " 1196 << feStr().name() <<
endl;
1202 meshRefiner_.mesh().time().path()
1203 /
"implicitFeaturePoint_" +
name(iter) +
".obj" 1206 Info<<
"Dumping implicit feature-point direction to " 1207 << fpStr().name() <<
endl;
1211 DynamicList<point> surfacePoints(4);
1212 DynamicList<vector> surfaceNormals(4);
1218 pointConstraint constraint;
1220 featureAttractionUsingReconstruction
1231 pointFaceSurfNormals,
1246 (constraint.first() > patchConstraints[pointi].first())
1248 (constraint.first() == patchConstraints[pointi].first())
1249 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
1253 patchAttraction[pointi] = attraction;
1254 patchConstraints[pointi] = constraint;
1258 if (feStr && patchConstraints[pointi].first() == 2)
1260 feStr().writeLine(pt, pt+patchAttraction[pointi]);
1262 else if (fpStr && patchConstraints[pointi].first() == 3)
1264 fpStr().writeLine(pt, pt+patchAttraction[pointi]);
1271 void Foam::snappySnapDriver::stringFeatureEdges
1274 const scalar featureCos,
1280 const List<pointConstraint>& rawPatchConstraints,
1283 List<pointConstraint>& patchConstraints
1313 forAll(pointEdges, pointi)
1315 if (patchConstraints[pointi].first() == 2)
1318 const labelList& pEdges = pointEdges[pointi];
1319 const vector& featVec = patchConstraints[pointi].second();
1323 bool hasPos =
false;
1324 bool hasNeg =
false;
1328 const edge&
e =
pp.
edges()[pEdges[pEdgei]];
1329 label nbrPointi =
e.otherVertex(pointi);
1331 if (patchConstraints[nbrPointi].first() > 1)
1334 const point featPt =
1335 nbrPt + patchAttraction[nbrPointi];
1336 const scalar cosAngle = (featVec & (featPt-pt));
1349 if (!hasPos || !hasNeg)
1355 label bestPosPointi = -1;
1356 scalar minPosDistSqr = GREAT;
1357 label bestNegPointi = -1;
1358 scalar minNegDistSqr = GREAT;
1362 const edge&
e =
pp.
edges()[pEdges[pEdgei]];
1363 label nbrPointi =
e.otherVertex(pointi);
1367 patchConstraints[nbrPointi].first() <= 1
1368 && rawPatchConstraints[nbrPointi].first() > 1
1371 const vector& nbrFeatVec =
1372 rawPatchConstraints[pointi].second();
1374 if (
mag(featVec&nbrFeatVec) > featureCos)
1381 rawPatchAttraction[nbrPointi]
1384 const point featPt =
1386 + rawPatchAttraction[nbrPointi];
1387 const scalar cosAngle =
1388 (featVec & (featPt-pt));
1392 if (!hasPos && d2 < minPosDistSqr)
1395 bestPosPointi = nbrPointi;
1400 if (!hasNeg && d2 < minNegDistSqr)
1403 bestNegPointi = nbrPointi;
1410 if (bestPosPointi != -1)
1420 patchAttraction[bestPosPointi] =
1421 0.5*rawPatchAttraction[bestPosPointi];
1422 patchConstraints[bestPosPointi] =
1423 rawPatchConstraints[bestPosPointi];
1427 if (bestNegPointi != -1)
1437 patchAttraction[bestNegPointi] =
1438 0.5*rawPatchAttraction[bestNegPointi];
1439 patchConstraints[bestNegPointi] =
1440 rawPatchConstraints[bestNegPointi];
1449 reduce(nChanged, sumOp<label>());
1450 Info<<
"Stringing feature edges : changed " << nChanged <<
" points" 1460 void Foam::snappySnapDriver::releasePointsNextToMultiPatch
1463 const scalar featureCos,
1468 const List<List<point>>& pointFaceCentres,
1472 const List<pointConstraint>& rawPatchConstraints,
1475 List<pointConstraint>& patchConstraints
1478 autoPtr<OBJstream> multiPatchStr;
1485 meshRefiner_.mesh().time().path()
1486 /
"multiPatch_" +
name(iter) +
".obj" 1489 Info<<
"Dumping removed constraints due to same-face" 1490 <<
" multi-patch points to " 1491 << multiPatchStr().name() <<
endl;
1496 bitSet isMultiPatchPoint(
pp.size());
1498 forAll(pointFacePatchID, pointi)
1503 pointFacePatchID[pointi],
1504 pointFaceCentres[pointi]
1506 isMultiPatchPoint.set(pointi, multiPatchPt.hit());
1510 forAll(isMultiPatchPoint, pointi)
1512 if (isMultiPatchPoint.test(pointi))
1516 patchConstraints[pointi].first() <= 1
1517 && rawPatchConstraints[pointi].first() > 1
1520 patchAttraction[pointi] = rawPatchAttraction[pointi];
1521 patchConstraints[pointi] = rawPatchConstraints[pointi];
1544 label nMultiPatchPoints = 0;
1547 label pointi =
f[fp];
1550 isMultiPatchPoint.test(pointi)
1551 && patchConstraints[pointi].first() > 1
1554 ++nMultiPatchPoints;
1558 if (nMultiPatchPoints > 0)
1562 label pointi =
f[fp];
1565 !isMultiPatchPoint.test(pointi)
1566 && patchConstraints[pointi].first() > 1
1572 patchAttraction[pointi] =
Zero;
1573 patchConstraints[pointi] = pointConstraint();
1585 reduce(nChanged, sumOp<label>());
1586 Info<<
"Removing constraints near multi-patch points : changed " 1587 << nChanged <<
" points" <<
endl;
1595 const List<pointConstraint>& patchConstraints,
1607 for (label startFp = 0; startFp <
f.
size()-2; startFp++)
1614 endFp <
f.
size() && endFp != minFp;
1620 patchConstraints[
f[startFp]].first() >= 2
1621 && patchConstraints[
f[endFp]].first() >= 2
1624 attractIndices =
labelPair(startFp, endFp);
1630 return attractIndices;
1634 bool Foam::snappySnapDriver::isSplitAlignedWithFeature
1636 const scalar featureCos,
1638 const pointConstraint& pc0,
1640 const pointConstraint& pc1
1644 scalar magD =
mag(d);
1657 if (pc0.first() == 2 &&
mag(d & pc0.second()) > featureCos)
1661 else if (pc1.first() == 2 &&
mag(d & pc1.second()) > featureCos)
1674 bool Foam::snappySnapDriver::isConcave
1680 const scalar concaveCos
1684 scalar magN0 =
mag(n0);
1693 scalar d = (
c1-c0)&n0;
1704 scalar magN1 =
mag(n1);
1712 if ((n0&n1) < concaveCos)
1726 const scalar featureCos,
1727 const scalar concaveCos,
1728 const scalar minAreaRatio,
1731 const List<pointConstraint>& patchConstraints,
1737 DynamicField<point>& points1
1744 if (localF.size() >= 4)
1784 for (label startFp = 0; startFp < localF.size()-2; startFp++)
1786 label minFp = localF.rcIndex(startFp);
1790 label endFp = localF.fcIndex(localF.fcIndex(startFp));
1791 endFp < localF.size() && endFp != minFp;
1795 label startPti = localF[startFp];
1796 label endPti = localF[endFp];
1798 const pointConstraint& startPc = patchConstraints[startPti];
1799 const pointConstraint& endPc = patchConstraints[endPti];
1801 if (startPc.first() >= 2 && endPc.first() >= 2)
1803 if (startPc.first() == 2 || endPc.first() == 2)
1808 point start = localPts[startPti]+patchAttr[startPti];
1809 point end = localPts[endPti]+patchAttr[endPti];
1813 !isSplitAlignedWithFeature
1839 f0.setSize(endFp-startFp+1);
1841 for (label fp = startFp; fp <= endFp; fp++)
1843 f0[i0++] = localF[fp];
1847 const face compact0(
identity(f0.size()));
1850 for (label fp=1; fp < f0.size()-1; fp++)
1857 localPts[f0.last()] + patchAttr[f0.last()]
1863 f1.setSize(localF.size()+2-f0.size());
1869 fp = localF.fcIndex(fp)
1872 f1[i1++] = localF[fp];
1874 f1[i1++] = localF[startFp];
1878 const face compact1(
identity(f1.size()));
1880 points1.append(localPts[f1[0]] + patchAttr[f1[0]]);
1881 for (label fp=1; fp < f1.size()-1; fp++)
1884 points1.append(localPts[
pi] + nearestAttr[
pi]);
1888 localPts[f1.last()] + patchAttr[f1.last()]
1902 compact1.centre(points1),
1903 compact1.areaNormal(points1),
1923 const scalar area0 = f0.mag(localPts);
1924 const scalar area1 = f1.mag(localPts);
1928 area0/area1 >= minAreaRatio
1929 && area1/area0 >= minAreaRatio
1932 attractIndices =
labelPair(startFp, endFp);
1939 return attractIndices;
1943 void Foam::snappySnapDriver::splitDiagonals
1945 const scalar featureCos,
1946 const scalar concaveCos,
1947 const scalar minAreaRatio,
1954 List<pointConstraint>& patchConstraints,
1955 DynamicList<label>& splitFaces,
1956 DynamicList<labelPair>& splits
1962 splitFaces.setCapacity(bFaces.size());
1964 splits.setCapacity(bFaces.size());
1968 DynamicField<point> facePoints0;
1969 DynamicField<point> facePoints1;
1975 findDiagonalAttraction
1996 splitFaces.append(bFaces[facei]);
1997 splits.append(
split);
2009 && patchConstraints[
f[fp]].first() >= 2
2012 patchConstraints[
f[fp]] = pointConstraint
2014 Tuple2<label, vector>
2017 nearestNormal[
f[fp]]
2020 patchAttraction[
f[fp]] = nearestAttraction[
f[fp]];
2028 void Foam::snappySnapDriver::avoidDiagonalAttraction
2031 const scalar featureCos,
2036 List<pointConstraint>& patchConstraints
2051 if (
diag[0] != -1 &&
diag[1] != -1)
2055 const label i0 =
f[
diag[0]];
2058 const label i1 =
f[
diag[1]];
2061 const point mid = 0.5*(pt0+pt1);
2063 const scalar cosAngle =
mag 2065 patchConstraints[i0].second()
2066 & patchConstraints[i1].second()
2075 if (cosAngle > featureCos)
2081 scalar minDistSqr = GREAT;
2084 label pointi =
f[fp];
2085 if (patchConstraints[pointi].first() <= 1)
2088 if (distSqr < minDistSqr)
2096 label minPointi =
f[minFp];
2097 patchAttraction[minPointi] =
2099 patchConstraints[minPointi] = patchConstraints[
f[
diag[0]]];
2125 Foam::snappySnapDriver::findNearFeatureEdge
2127 const bool isRegionEdge,
2132 const point& estimatedPt,
2134 List<List<DynamicList<point>>>& edgeAttractors,
2135 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2137 List<pointConstraint>& patchConstraints
2140 const refinementFeatures& features = meshRefiner_.features();
2143 List<pointIndexHit> nearEdgeInfo;
2148 features.findNearestRegionEdge
2159 features.findNearestEdge
2170 label feati = nearEdgeFeat[0];
2176 edgeAttractors[feati][nearInfo.index()].append(nearInfo.point());
2177 pointConstraint
c(Tuple2<label, vector>(2, nearNormal[0]));
2178 edgeConstraints[feati][nearInfo.index()].append(
c);
2181 patchAttraction[pointi] = nearInfo.point()-
pp.
localPoints()[pointi];
2182 patchConstraints[pointi] =
c;
2184 return Tuple2<label, pointIndexHit>(feati, nearInfo);
2189 Foam::snappySnapDriver::findNearFeaturePoint
2191 const bool isRegionPoint,
2196 const point& estimatedPt,
2199 List<labelList>& pointAttractor,
2200 List<List<pointConstraint>>& pointConstraints,
2202 List<List<DynamicList<point>>>& edgeAttractors,
2203 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2206 List<pointConstraint>& patchConstraints
2209 const refinementFeatures& features = meshRefiner_.features();
2214 List<pointIndexHit> nearInfo;
2215 features.findNearestPoint
2223 label feati = nearFeat[0];
2229 label featPointi = nearInfo[0].index();
2230 const point& featPt = nearInfo[0].hitPoint();
2231 scalar distSqr = featPt.distSqr(pt);
2234 label oldPointi = pointAttractor[feati][featPointi];
2236 if (oldPointi != -1)
2248 pointAttractor[feati][featPointi] = pointi;
2249 pointConstraints[feati][featPointi].first() = 3;
2250 pointConstraints[feati][featPointi].second() =
Zero;
2253 patchAttraction[pointi] = featPt-pt;
2254 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2257 patchAttraction[oldPointi] =
Zero;
2258 patchConstraints[oldPointi] = pointConstraint();
2279 pointAttractor[feati][featPointi] = pointi;
2280 pointConstraints[feati][featPointi].first() = 3;
2281 pointConstraints[feati][featPointi].second() =
Zero;
2284 patchAttraction[pointi] = featPt-pt;
2285 patchConstraints[pointi] = pointConstraints[feati][featPointi];
2289 return Tuple2<label, pointIndexHit>(feati, nearInfo[0]);
2295 void Foam::snappySnapDriver::determineFeatures
2298 const scalar featureCos,
2299 const bool multiRegionFeatureSnap,
2305 const List<List<point>>& pointFaceSurfNormals,
2306 const List<List<point>>& pointFaceDisp,
2307 const List<List<point>>& pointFaceCentres,
2311 List<labelList>& pointAttractor,
2312 List<List<pointConstraint>>& pointConstraints,
2314 List<List<DynamicList<point>>>& edgeAttractors,
2315 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2318 List<pointConstraint>& patchConstraints
2321 autoPtr<OBJstream> featureEdgeStr;
2322 autoPtr<OBJstream> missedEdgeStr;
2323 autoPtr<OBJstream> featurePointStr;
2324 autoPtr<OBJstream> missedMP0Str;
2325 autoPtr<OBJstream> missedMP1Str;
2329 featureEdgeStr.reset
2333 meshRefiner_.mesh().time().path()
2334 /
"featureEdge_" +
name(iter) +
".obj" 2337 Info<<
"Dumping feature-edge sampling to " 2338 << featureEdgeStr().name() <<
endl;
2344 meshRefiner_.mesh().time().path()
2345 /
"missedFeatureEdge_" +
name(iter) +
".obj" 2348 Info<<
"Dumping feature-edges that are too far away to " 2349 << missedEdgeStr().name() <<
endl;
2351 featurePointStr.reset
2355 meshRefiner_.mesh().time().path()
2356 /
"featurePoint_" +
name(iter) +
".obj" 2359 Info<<
"Dumping feature-point sampling to " 2360 << featurePointStr().name() <<
endl;
2366 meshRefiner_.mesh().time().path()
2367 /
"missedFeatureEdgeFromMPEdge_" +
name(iter) +
".obj" 2370 Info<<
"Dumping region-edges that are too far away to " 2371 << missedMP0Str().name() <<
endl;
2377 meshRefiner_.mesh().time().path()
2378 /
"missedFeatureEdgeFromMPPoint_" +
name(iter) +
".obj" 2381 Info<<
"Dumping region-points that are too far away to " 2382 << missedMP1Str().name() <<
endl;
2386 DynamicList<point> surfacePoints(4);
2387 DynamicList<vector> surfaceNormals(4);
2404 pointConstraint constraint;
2406 featureAttractionUsingReconstruction
2417 pointFaceSurfNormals,
2462 (constraint.first() > patchConstraints[pointi].first())
2464 (constraint.first() == patchConstraints[pointi].first())
2465 && (
magSqr(attraction) <
magSqr(patchAttraction[pointi]))
2469 patchAttraction[pointi] = attraction;
2470 patchConstraints[pointi] = constraint;
2473 if (patchConstraints[pointi].first() == 1)
2476 if (multiRegionFeatureSnap)
2478 const point estimatedPt(pt + nearestDisp[pointi]);
2484 pointFacePatchID[pointi],
2490 if (multiPatchPt.hit())
2495 Tuple2<label, pointIndexHit> nearInfo =
2502 multiPatchPt.point(),
2517 featureEdgeStr().writeLine
2528 missedEdgeStr().writeLine
2531 multiPatchPt.point()
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().writeLine(pt, estimatedPt);
2596 nearInfo = findNearFeaturePoint
2625 if (!nearInfo.second().hit())
2627 nearInfo = findNearFeatureEdge
2645 if (missedMP1Str && !nearInfo.second().hit())
2647 missedMP1Str().writeLine(pt, estimatedPt);
2659 nearInfo = findNearFeatureEdge
2683 && patchConstraints[pointi].first() == 3
2686 featurePointStr().writeLine(pt, info.point());
2691 && patchConstraints[pointi].first() == 2
2694 featureEdgeStr().writeLine(pt, info.point());
2701 missedEdgeStr().writeLine(pt, estimatedPt);
2705 else if (patchConstraints[pointi].first() == 3)
2708 const point estimatedPt(pt + patchAttraction[pointi]);
2712 if (multiRegionFeatureSnap)
2719 pointFacePatchID[pointi],
2724 if (multiPatchPt.hit())
2727 nearInfo = findNearFeaturePoint
2748 nearInfo = findNearFeaturePoint
2771 nearInfo = findNearFeaturePoint
2792 if (featurePointStr && info.hit())
2794 featurePointStr().writeLine(pt, info.point());
2813 void Foam::snappySnapDriver::determineBaffleFeatures
2816 const bool baffleFeaturePoints,
2817 const scalar featureCos,
2823 List<labelList>& pointAttractor,
2824 List<List<pointConstraint>>& pointConstraints,
2826 List<List<DynamicList<point>>>& edgeAttractors,
2827 List<List<DynamicList<pointConstraint>>>& edgeConstraints,
2830 List<pointConstraint>& patchConstraints
2833 const fvMesh&
mesh = meshRefiner_.mesh();
2834 const refinementFeatures& features = meshRefiner_.features();
2837 List<List<point>> edgeFaceNormals(
pp.
nEdges());
2843 List<point>& eFc = edgeFaceNormals[edgei];
2847 label facei = eFaces[i];
2864 listPlusEqOp<point>(),
2877 autoPtr<OBJstream> baffleEdgeStr;
2884 meshRefiner_.mesh().time().path()
2885 /
"baffleEdge_" +
name(iter) +
".obj" 2888 Info<<
nl <<
"Dumping baffle-edges to " 2889 << baffleEdgeStr().name() <<
endl;
2895 label nBaffleEdges = 0;
2902 forAll(edgeFaceNormals, edgei)
2904 const List<point>& efn = edgeFaceNormals[edgei];
2906 if (efn.size() == 2 && (efn[0]&efn[1]) < baffleFeatureCos)
2908 isBaffleEdge.set(edgei);
2911 pointStatus[
e[0]] = 0;
2912 pointStatus[
e[1]] = 0;
2921 reduce(nBaffleEdges, sumOp<label>());
2923 Info<<
"Detected " << nBaffleEdges
2924 <<
" baffle edges out of " 2926 <<
" edges." <<
endl;
2951 label nBafflePoints = 0;
2952 forAll(pointStatus, pointi)
2954 if (pointStatus[pointi] != -1)
2959 reduce(nBafflePoints, sumOp<label>());
2962 label nPointAttract = 0;
2963 label nEdgeAttract = 0;
2965 forAll(pointStatus, pointi)
2969 if (pointStatus[pointi] == 0)
2973 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
2992 if (nearInfo.second().hit())
2996 if (baffleFeaturePoints)
2998 nearInfo = findNearFeaturePoint
3018 if (nearInfo.first() != -1)
3026 else if (pointStatus[pointi] == 1)
3029 List<pointIndexHit> nearInfo;
3030 features.findNearestPoint
3038 label feati = nearFeat[0];
3044 label featPointi = nearInfo[0].index();
3045 const point& featPt = nearInfo[0].hitPoint();
3046 scalar distSqr = featPt.distSqr(pt);
3049 label oldPointi = pointAttractor[feati][featPointi];
3060 pointAttractor[feati][featPointi] = pointi;
3061 pointConstraints[feati][featPointi].first() = 3;
3062 pointConstraints[feati][featPointi].second() =
Zero;
3065 patchAttraction[pointi] = featPt-pt;
3066 patchConstraints[pointi] =
3067 pointConstraints[feati][featPointi];
3069 if (oldPointi != -1)
3105 Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
3119 if (nearInfo.first() != -1)
3127 reduce(nPointAttract, sumOp<label>());
3128 reduce(nEdgeAttract, sumOp<label>());
3130 Info<<
"Baffle points : " << nBafflePoints
3131 <<
" of which attracted to :" <<
nl 3132 <<
" feature point : " << nPointAttract <<
nl 3133 <<
" feature edge : " << nEdgeAttract <<
nl 3134 <<
" rest : " << nBafflePoints-nPointAttract-nEdgeAttract
3140 void Foam::snappySnapDriver::reverseAttractMeshPoints
3148 const List<labelList>& pointAttractor,
3149 const List<List<pointConstraint>>& pointConstraints,
3151 const List<List<DynamicList<point>>>& edgeAttractors,
3152 const List<List<DynamicList<pointConstraint>>>& edgeConstraints,
3155 const List<pointConstraint>& rawPatchConstraints,
3159 List<pointConstraint>& patchConstraints
3162 const refinementFeatures& features = meshRefiner_.features();
3177 bb.inflate(
rndGen, 1
e-4, ROOTVSMALL);
3181 DynamicList<label> attractPoints(
pp.
nPoints());
3183 const fvMesh&
mesh = meshRefiner_.mesh();
3187 forAll(rawPatchConstraints, pointi)
3189 if (rawPatchConstraints[pointi].first() >= 2)
3191 isFeatureEdgeOrPoint[pointi] =
true;
3197 <<
" mesh points out of " 3199 <<
" for reverse attraction." <<
endl;
3207 isFeatureEdgeOrPoint,
3212 for (label nGrow = 0; nGrow < 1; nGrow++)
3214 boolList newIsFeatureEdgeOrPoint(isFeatureEdgeOrPoint);
3222 if (isFeatureEdgeOrPoint[
f[fp]])
3227 newIsFeatureEdgeOrPoint[
f[fp]] =
true;
3234 isFeatureEdgeOrPoint = newIsFeatureEdgeOrPoint;
3240 isFeatureEdgeOrPoint,
3248 forAll(isFeatureEdgeOrPoint, pointi)
3250 if (isFeatureEdgeOrPoint[pointi])
3252 attractPoints.append(pointi);
3258 <<
" mesh points out of " 3260 <<
" for reverse attraction." <<
endl;
3264 indexedOctree<treeDataPoint> ppTree
3275 patchAttraction =
Zero;
3277 patchConstraints = pointConstraint();
3279 forAll(edgeAttractors, feati)
3281 const List<DynamicList<point>>& edgeAttr = edgeAttractors[feati];
3282 const List<DynamicList<pointConstraint>>& edgeConstr =
3283 edgeConstraints[feati];
3285 forAll(edgeAttr, featEdgei)
3287 const DynamicList<point>& attr = edgeAttr[featEdgei];
3291 const point& featPt = attr[i];
3301 ppTree.shapes().objectIndex(nearInfo.index());
3303 const point attraction
3306 - ppTree.shapes()[nearInfo.index()]
3313 patchConstraints[pointi].first() <= 1
3314 ||
magSqr(attraction) <
magSqr(patchAttraction[pointi])
3317 patchAttraction[pointi] = attraction;
3318 patchConstraints[pointi] = edgeConstr[featEdgei][i];
3323 static label nWarn = 0;
3328 <<
"Did not find pp point near " << featPt
3334 <<
"Reached warning limit " << nWarn
3335 <<
". Suppressing further warnings." <<
endl;
3354 forAll(pointAttractor, feati)
3356 const labelList& pointAttr = pointAttractor[feati];
3357 const List<pointConstraint>& pointConstr = pointConstraints[feati];
3359 forAll(pointAttr, featPointi)
3361 if (pointAttr[featPointi] != -1)
3363 const point& featPt = features[feati].points()
3378 ppTree.shapes().objectIndex(nearInfo.index());
3380 const point attraction
3383 - ppTree.shapes()[nearInfo.index()]
3389 if (patchConstraints[pointi].first() <= 1)
3391 patchAttraction[pointi] = attraction;
3392 patchConstraints[pointi] = pointConstr[featPointi];
3394 else if (patchConstraints[pointi].first() == 2)
3396 patchAttraction[pointi] = attraction;
3397 patchConstraints[pointi] = pointConstr[featPointi];
3399 else if (patchConstraints[pointi].first() == 3)
3405 <
magSqr(patchAttraction[pointi])
3408 patchAttraction[pointi] = attraction;
3409 patchConstraints[pointi] =
3410 pointConstr[featPointi];
3420 void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
3423 const bool multiRegionFeatureSnap,
3425 const bool detectBaffles,
3426 const bool baffleFeaturePoints,
3428 const bool releasePoints,
3429 const bool stringFeatures,
3430 const bool avoidDiagonal,
3432 const scalar featureCos,
3439 const List<List<point>>& pointFaceSurfNormals,
3440 const List<List<point>>& pointFaceDisp,
3441 const List<List<point>>& pointFaceCentres,
3445 List<pointConstraint>& patchConstraints
3448 const refinementFeatures& features = meshRefiner_.features();
3449 const fvMesh&
mesh = meshRefiner_.mesh();
3451 const bitSet isPatchMasterPoint
3466 List<List<DynamicList<point>>> edgeAttractors(features.size());
3467 List<List<DynamicList<pointConstraint>>> edgeConstraints
3473 label nFeatEdges = features[feati].edges().size();
3474 edgeAttractors[feati].setSize(nFeatEdges);
3475 edgeConstraints[feati].setSize(nFeatEdges);
3481 List<labelList> pointAttractor(features.size());
3482 List<List<pointConstraint>> pointConstraints(features.size());
3485 label nFeatPoints = features[feati].points().size();
3486 pointAttractor[feati].setSize(nFeatPoints, -1);
3487 pointConstraints[feati].setSize(nFeatPoints);
3492 List<pointConstraint> rawPatchConstraints(
pp.
nPoints());
3498 multiRegionFeatureSnap,
3504 pointFaceSurfNormals,
3523 Info<<
"Raw geometric feature analysis : ";
3524 writeStats(
pp, isPatchMasterPoint, rawPatchConstraints);
3540 determineBaffleFeatures
3543 baffleFeaturePoints,
3564 Info<<
"After baffle feature analysis : ";
3565 writeStats(
pp, isPatchMasterPoint, rawPatchConstraints);
3573 reverseAttractMeshPoints
3589 rawPatchConstraints,
3599 Info<<
"Reverse attract feature analysis : ";
3606 OBJstream featureEdgeStr
3608 meshRefiner_.mesh().time().path()
3609 /
"edgeAttractors_" +
name(iter) +
".obj" 3611 Info<<
"Dumping feature-edge attraction to " 3612 << featureEdgeStr.name() <<
endl;
3614 OBJstream featurePointStr
3616 meshRefiner_.mesh().time().path()
3617 /
"pointAttractors_" +
name(iter) +
".obj" 3619 Info<<
"Dumping feature-point attraction to " 3620 << featurePointStr.name() <<
endl;
3622 forAll(patchConstraints, pointi)
3625 const vector& attr = patchAttraction[pointi];
3627 if (patchConstraints[pointi].first() == 2)
3629 featureEdgeStr.writeLine(pt, pt+attr);
3631 else if (patchConstraints[pointi].first() == 3)
3633 featurePointStr.writeLine(pt, pt+attr);
3645 releasePointsNextToMultiPatch
3657 rawPatchConstraints,
3680 rawPatchConstraints,
3693 avoidDiagonalAttraction
3708 meshRefiner_.mesh().time().path()
3709 /
"patchAttraction_" +
name(iter) +
".obj",
3718 void Foam::snappySnapDriver::preventFaceSqueeze
3721 const scalar featureCos,
3728 List<pointConstraint>& patchConstraints
3731 autoPtr<OBJstream> strPtr;
3738 meshRefiner_.mesh().time().path()
3739 /
"faceSqueeze_" +
name(iter) +
".obj" 3742 Info<<
"Dumping faceSqueeze corrections to " 3743 << strPtr().name() <<
endl;
3755 singleF.setSize(
f.
size());
3756 for (label i = 0; i <
f.
size(); i++)
3761 label nConstraints = 0;
3764 label pointi =
f[fp];
3767 if (patchConstraints[pointi].first() > 1)
3769 points[fp] = pt + patchAttraction[pointi];
3778 if (nConstraints ==
f.
size())
3792 scalar
s = pt.distSqr(nextPt);
3801 label pointi =
f.prevLabel(maxFp);
3812 patchAttraction[pointi] = nearestAttraction[pointi];
3816 strPtr().writeLine(pt, pt+patchAttraction[pointi]);
3823 scalar newArea = singleF.mag(
points);
3824 if (newArea < 0.1*oldArea)
3831 scalar
s =
magSqr(patchAttraction[
f[fp]]);
3840 label pointi =
f[maxFp];
3842 patchAttraction[pointi] *= 0.5;
3853 const snapParameters& snapParams,
3854 const bool alignMeshEdges,
3856 const scalar featureCos,
3857 const scalar featureAttract,
3861 motionSmoother& meshMover,
3863 List<pointConstraint>& patchConstraints,
3865 DynamicList<label>& splitFaces,
3866 DynamicList<labelPair>& splits
3876 const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
3877 const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
3878 const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
3880 Info<<
"Overriding displacement on features :" <<
nl 3881 <<
" implicit features : " << implicitFeatureAttraction <<
nl 3882 <<
" explicit features : " << explicitFeatureAttraction <<
nl 3883 <<
" multi-patch features : " << multiRegionFeatureSnap <<
nl 3888 const fvMesh&
mesh = meshRefiner_.mesh();
3892 const bitSet isPatchMasterPoint
3905 List<List<point>> pointFaceSurfNormals;
3906 List<List<point>> pointFaceDisp;
3907 List<List<point>> pointFaceCentres;
3908 List<labelList> pointFacePatchID;
3919 faceSnapDist[facei] =
max(faceSnapDist[facei], snapDist[
f[fp]]);
3941 faceSurfaceGlobalRegion
3951 calcNearestFacePointProperties
3958 faceSurfaceGlobalRegion,
3960 pointFaceSurfNormals,
3979 patchAttraction.
setSize(localPoints.size());
3980 patchAttraction =
Zero;
3982 patchConstraints.setSize(localPoints.size());
3983 patchConstraints = pointConstraint();
3985 if (implicitFeatureAttraction)
3990 featureAttractionUsingReconstruction
3999 pointFaceSurfNormals,
4009 if (explicitFeatureAttraction)
4012 bool releasePoints =
false;
4013 bool stringFeatures =
false;
4014 bool avoidDiagonal =
false;
4017 releasePoints = snapParams.releasePoints();
4018 stringFeatures = snapParams.stringFeatures();
4019 avoidDiagonal = snapParams.avoidDiagonal();
4029 featureAttractionUsingFeatureEdges
4032 multiRegionFeatureSnap,
4034 snapParams.detectBaffles(),
4035 snapParams.baffleFeaturePoints(),
4049 pointFaceSurfNormals,
4059 if (!alignMeshEdges)
4063 degToRad(snapParams.concaveAngle())
4065 const scalar minAreaRatio = snapParams.minAreaRatio();
4067 Info<<
"Experimental: introducing face splits to avoid rotating" 4068 <<
" mesh edges. Splitting faces when" <<
nl 4069 <<
indent <<
"- angle not concave by more than " 4070 << snapParams.concaveAngle() <<
" degrees" <<
nl 4071 <<
indent <<
"- resulting triangles of similar area " 4072 <<
" (ratio within " << minAreaRatio <<
")" <<
nl 4093 Info<<
"Diagonal attraction feature correction : ";
4126 <<
" avg:" << avgPatchDisp <<
endl 4127 <<
" feature : max:" <<
gMaxMagSqr(patchAttraction)
4128 <<
" avg:" << avgPatchAttr <<
endl;
4142 forAll(patchConstraints, pointi)
4144 if (patchConstraints[pointi].first() > 1)
4147 (1.0-featureAttract)*patchDisp[pointi]
4148 + featureAttract*patchAttraction[pointi];
4156 Info<<
"Feature analysis : ";
4172 if (featureAttract < 1-0.001)
4179 const bitSet isPatchMasterEdge
4211 tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
4222 tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
4229 /
"tangPatchDispConstrained_" +
name(iter) +
".obj",
4236 (1.0-featureAttract)*smoothedPatchDisp
4237 + featureAttract*tangPatchDisp;
4241 const scalar
relax = featureAttract;
4252 minMagSqrEqOp<point>(),
4253 vector(GREAT, GREAT, GREAT)
const labelListList & pointEdges() const
Return point-edge addressing.
const labelList patchIDs(pbm.indices(polyPatchNames, true))
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.
fileName path() const
Return path = rootPath/caseName. Same as TimePaths::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 expressions::valueTypeCode::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), works like std::iota() but returning a...
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.
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.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
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)