58 for (
direction i = 0; i < vector::nComponents; i++)
88 return less(values_[a], values_[
b]);
134 label nMasterChanged = 0;
155 bitSet refinedInternalFace(nInternalFaces);
159 for (label faceI = 0; faceI < nInternalFaces; faceI++)
161 label oldOwn = map.
cellMap()[faceOwner[faceI]];
162 label oldNei = map.
cellMap()[faceNeighbour[faceI]];
166 oldOwn >= 0 && !oldRefineCell.test(oldOwn)
167 && oldNei >= 0 && !oldRefineCell.test(oldNei)
174 refinedInternalFace.
set(faceI);
187 label faceI = pp.
start();
191 label oldOwn = map.
cellMap()[faceOwner[faceI]];
193 if (oldOwn >= 0 && !oldRefineCell.test(oldOwn))
199 refinedBoundaryFace[faceI-nInternalFaces] =
true;
219 forAll(refinedInternalFace, faceI)
221 if (refinedInternalFace.test(faceI))
223 const cell& ownFaces =
cells[faceOwner[faceI]];
226 changedFace[ownFaces[ownI]] =
true;
228 const cell& neiFaces =
cells[faceNeighbour[faceI]];
231 changedFace[neiFaces[neiI]] =
true;
236 forAll(refinedBoundaryFace, i)
238 if (refinedBoundaryFace[i])
240 const cell& ownFaces =
cells[faceOwner[i+nInternalFaces]];
243 changedFace[ownFaces[ownI]] =
true;
264 forAll(changedFace, faceI)
266 if (changedFace[faceI] && isMasterFace.
test(faceI))
276 Pout<<
"getChangedFaces : Detected " 277 <<
" local:" << changedFaces.
size()
282 faceSet changedFacesSet(
mesh,
"changedFaces", changedFaces);
283 Pout<<
"getChangedFaces : Writing " << changedFaces.
size()
284 <<
" changed faces to faceSet " << changedFacesSet.
name()
286 changedFacesSet.
write();
296 bool Foam::meshRefinement::markForRefine
298 const label markValue,
299 const label nAllowRefine,
307 cellValue = markValue;
311 return nRefine <= nAllowRefine;
315 void Foam::meshRefinement::markFeatureCellLevel
337 Cloud<trackedParticle> startPointCloud
341 IDLList<trackedParticle>()
349 for (
const point& keepPoint : keepPoints)
351 const label celli = mesh_.cellTree().findInside(keepPoint);
359 const edgeMesh& featureMesh = features_[feati];
360 const label featureLevel = features_.levels()[feati][0];
365 label nRegions = featureMesh.regions(edgeRegion);
368 bitSet regionVisited(nRegions);
374 forAll(pointEdges, pointi)
376 if (pointEdges[pointi].size() != 2)
380 Pout<<
"Adding particle from point:" << pointi
381 <<
" coord:" << featureMesh.points()[pointi]
382 <<
" since number of emanating edges:" 383 << pointEdges[pointi].size()
388 startPointCloud.addParticle
395 featureMesh.points()[pointi],
404 if (pointEdges[pointi].size() > 0)
406 label e0 = pointEdges[pointi][0];
407 label regioni = edgeRegion[e0];
408 regionVisited.set(regioni);
416 forAll(featureMesh.edges(), edgei)
418 if (regionVisited.set(edgeRegion[edgei]))
420 const edge&
e = featureMesh.edges()[edgei];
421 label pointi =
e.start();
424 Pout<<
"Adding particle from point:" << pointi
425 <<
" coord:" << featureMesh.points()[pointi]
426 <<
" on circular region:" << edgeRegion[edgei]
431 startPointCloud.addParticle
438 featureMesh.points()[pointi],
453 maxFeatureLevel =
labelList(mesh_.nCells(), -1);
456 List<bitSet> featureEdgeVisited(features_.size());
460 featureEdgeVisited[featI].setSize(features_[featI].edges().size());
461 featureEdgeVisited[featI] =
false;
465 trackedParticle::trackingData td
476 scalar maxTrackLen = 2.0*mesh_.bounds().mag();
480 Pout<<
"Tracking " << startPointCloud.size()
481 <<
" particles over distance " << maxTrackLen
482 <<
" to find the starting cell" <<
endl;
484 startPointCloud.move(startPointCloud, td, maxTrackLen);
488 maxFeatureLevel = -1;
491 featureEdgeVisited[featI] =
false;
495 Cloud<trackedParticle> cloud
499 IDLList<trackedParticle>()
504 Pout<<
"Constructing cloud for cell marking" <<
endl;
507 for (
const trackedParticle& startTp : startPointCloud)
509 const label featI = startTp.i();
510 const label pointI = startTp.j();
512 const edgeMesh& featureMesh = features_[featI];
513 const labelList& pEdges = featureMesh.pointEdges()[pointI];
518 label edgeI = pEdges[pEdgeI];
520 if (featureEdgeVisited[featI].
set(edgeI))
525 const edge&
e = featureMesh.edges()[edgeI];
526 label otherPointi =
e.otherVertex(pointI);
528 trackedParticle* tp(
new trackedParticle(startTp));
529 tp->start() = tp->position();
530 tp->end() = featureMesh.points()[otherPointi];
531 tp->j() = otherPointi;
536 Pout<<
"Adding particle for point:" << pointI
537 <<
" coord:" << tp->position()
538 <<
" feature:" << featI
539 <<
" to track to:" << tp->end()
543 cloud.addParticle(tp);
548 startPointCloud.clear();
556 Pout<<
"Tracking " << cloud.size()
557 <<
" particles over distance " << maxTrackLen
558 <<
" to mark cells" <<
endl;
560 cloud.move(cloud, td, maxTrackLen);
563 for (trackedParticle& tp : cloud)
565 const label featI = tp.i();
566 const label pointI = tp.j();
568 const edgeMesh& featureMesh = features_[featI];
569 const labelList& pEdges = featureMesh.pointEdges()[pointI];
574 bool keepParticle =
false;
578 label edgeI = pEdges[i];
580 if (featureEdgeVisited[featI].
set(edgeI))
585 const edge&
e = featureMesh.edges()[edgeI];
586 label otherPointi =
e.otherVertex(pointI);
588 tp.start() = tp.position();
589 tp.end() = featureMesh.points()[otherPointi];
590 tp.j() = otherPointi;
601 cloud.deleteParticle(tp);
608 Pout<<
"Remaining particles " << cloud.size() <<
endl;
636 Foam::label Foam::meshRefinement::markFeatureRefinement
639 const label nAllowRefine,
647 markFeatureCellLevel(keepPoints, maxFeatureLevel);
652 const labelList& cellLevel = meshCutter_.cellLevel();
654 label oldNRefine = nRefine;
656 forAll(maxFeatureLevel, cellI)
658 if (maxFeatureLevel[cellI] > cellLevel[cellI])
684 Info<<
"Reached refinement limit." <<
endl;
687 return returnReduce(nRefine-oldNRefine, sumOp<label>());
692 Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
694 const label nAllowRefine,
700 const labelList& cellLevel = meshCutter_.cellLevel();
701 const pointField& cellCentres = mesh_.cellCentres();
704 if (features_.maxDistance() <= 0.0)
709 label oldNRefine = nRefine;
713 labelList testLevels(cellLevel.size()-nRefine);
718 if (refineCell[cellI] == -1)
720 testCc[testI] = cellCentres[cellI];
721 testLevels[testI] = cellLevel[cellI];
728 features_.findHigherLevel(testCc, testLevels, maxLevel);
735 if (refineCell[cellI] == -1)
737 if (maxLevel[testI] > testLevels[testI])
739 bool reachedLimit = !markForRefine
751 Pout<<
"Stopped refining internal cells" 752 <<
" since reaching my cell limit of " 753 << mesh_.nCells()+7*nRefine <<
endl;
768 Info<<
"Reached refinement limit." <<
endl;
771 return returnReduce(nRefine-oldNRefine, sumOp<label>());
776 Foam::label Foam::meshRefinement::markInternalRefinement
778 const label nAllowRefine,
784 const labelList& cellLevel = meshCutter_.cellLevel();
785 const pointField& cellCentres = mesh_.cellCentres();
787 label oldNRefine = nRefine;
791 labelList testLevels(cellLevel.size()-nRefine);
796 if (refineCell[cellI] == -1)
798 testCc[testI] = cellCentres[cellI];
799 testLevels[testI] = cellLevel[cellI];
806 shells_.findHigherLevel(testCc, testLevels, maxLevel);
813 if (refineCell[cellI] == -1)
815 if (maxLevel[testI] > testLevels[testI])
817 bool reachedLimit = !markForRefine
829 Pout<<
"Stopped refining internal cells" 830 <<
" since reaching my cell limit of " 831 << mesh_.nCells()+7*nRefine <<
endl;
846 Info<<
"Reached refinement limit." <<
endl;
849 return returnReduce(nRefine-oldNRefine, sumOp<label>());
853 Foam::label Foam::meshRefinement::unmarkInternalRefinement
859 const labelList& cellLevel = meshCutter_.cellLevel();
860 const pointField& cellCentres = mesh_.cellCentres();
862 label oldNRefine = nRefine;
871 if (refineCell[cellI] >= 0)
873 testCc[testI] = cellCentres[cellI];
874 testLevels[testI] = cellLevel[cellI];
881 limitShells_.findLevel(testCc, testLevels, levelShell);
888 if (refineCell[cellI] >= 0)
890 if (levelShell[testI] != -1)
892 refineCell[cellI] = -1;
899 return returnReduce(oldNRefine-nRefine, sumOp<label>());
914 const labelList& surfIndex = surfaceIndex();
921 if (surfIndex[faceI] != -1)
923 label own = mesh_.faceOwner()[faceI];
925 if (mesh_.isInternalFace(faceI))
927 label nei = mesh_.faceNeighbour()[faceI];
929 if (refineCell[own] == -1 || refineCell[nei] == -1)
931 testFaces[nTest++] = faceI;
936 const label bFacei = faceI - mesh_.nInternalFaces();
937 if (refineCell[own] == -1 || boundaryRefineCell[bFacei] == -1)
939 testFaces[nTest++] = faceI;
944 testFaces.setSize(nTest);
951 Foam::label Foam::meshRefinement::markSurfaceRefinement
953 const label nAllowRefine,
961 const labelList& cellLevel = meshCutter_.cellLevel();
963 label oldNRefine = nRefine;
971 labelList testFaces(getRefineCandidateFaces(refineCell));
996 surfaces_.findHigherIntersection
1013 label faceI = testFaces[i];
1015 label surfI = surfaceHit[i];
1025 label own = mesh_.faceOwner()[faceI];
1027 if (surfaceMinLevel[i] > cellLevel[own])
1045 if (mesh_.isInternalFace(faceI))
1047 label nei = mesh_.faceNeighbour()[faceI];
1048 if (surfaceMinLevel[i] > cellLevel[nei])
1075 Info<<
"Reached refinement limit." <<
endl;
1078 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1083 Foam::label Foam::meshRefinement::countMatches
1085 const List<point>& normals1,
1086 const List<point>& normals2,
1094 const vector& n1 = normals1[i];
1098 const vector& n2 = normals2[j];
1123 bool Foam::meshRefinement::highCurvature
1125 const scalar minCosAngle,
1126 const scalar lengthScale,
1133 const scalar cosAngle = (n0&n1);
1135 if (cosAngle < minCosAngle)
1140 else if (cosAngle > 1-1
e-6)
1145 else if (lengthScale > SMALL)
1150 const plane pl0(
p0, (n0 ^ axis));
1151 const scalar r1 = pl0.normalIntersect(p1, n1);
1155 const plane pl1(p1, (n1 ^ axis));
1156 const scalar r0 = pl1.normalIntersect(
p0, n0);
1161 const scalar radius = 0.5*(
mag(r1)+
mag(r0));
1163 if (radius < lengthScale)
1181 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
1183 const scalar curvature,
1184 const label nAllowRefine,
1192 const labelList& cellLevel = meshCutter_.cellLevel();
1193 const pointField& cellCentres = mesh_.cellCentres();
1195 label oldNRefine = nRefine;
1209 labelList testFaces(getRefineCandidateFaces(refineCell));
1219 label faceI = testFaces[i];
1221 label own = mesh_.faceOwner()[faceI];
1223 if (mesh_.isInternalFace(faceI))
1225 label nei = mesh_.faceNeighbour()[faceI];
1227 start[i] = cellCentres[own];
1228 end[i] = cellCentres[nei];
1229 minLevel[i] =
min(cellLevel[own], cellLevel[nei]);
1233 label bFaceI = faceI - mesh_.nInternalFaces();
1235 start[i] = cellCentres[own];
1236 end[i] = neiCc[bFaceI];
1238 if (!isMasterFace[faceI])
1240 std::swap(start[i],
end[i]);
1243 minLevel[i] =
min(cellLevel[own], neiLevel[bFaceI]);
1256 const bool hasCurvatureLevels = (
max(surfaces_.maxCurvatureLevel()) > 0);
1262 List<pointList> cellSurfLocations(mesh_.nCells());
1263 List<vectorList> cellSurfNormals(mesh_.nCells());
1267 List<pointList> surfaceLocation;
1269 List<vectorList> surfaceNormal;
1273 surfaces_.findAllIntersections
1280 labelList(surfaces_.maxLevel().size(), 0),
1281 surfaces_.maxLevel(),
1293 forAll(surfaceNormal, pointI)
1295 pointList& pLocations = surfaceLocation[pointI];
1296 vectorList& pNormals = surfaceNormal[pointI];
1297 labelList& pLevel = surfaceLevel[pointI];
1299 sortedOrder(pNormals, visitOrder, normalLess(pLocations));
1301 pLocations = List<point>(pLocations, visitOrder);
1302 pNormals = List<vector>(pNormals, visitOrder);
1377 label faceI = testFaces[i];
1378 label own = mesh_.faceOwner()[faceI];
1380 const pointList& fPoints = surfaceLocation[i];
1381 const vectorList& fNormals = surfaceNormal[i];
1382 const labelList& fLevels = surfaceLevel[i];
1386 if (fLevels[hitI] > cellLevel[own])
1388 cellSurfLevels[own].append(fLevels[hitI]);
1389 cellSurfLocations[own].append(fPoints[hitI]);
1390 cellSurfNormals[own].append(fNormals[hitI]);
1393 if (mesh_.isInternalFace(faceI))
1395 label nei = mesh_.faceNeighbour()[faceI];
1396 if (fLevels[hitI] > cellLevel[nei])
1398 cellSurfLevels[nei].append(fLevels[hitI]);
1399 cellSurfLocations[nei].append(fPoints[hitI]);
1400 cellSurfNormals[nei].append(fNormals[hitI]);
1414 forAll(cellSurfNormals, cellI)
1416 const vectorList& normals = cellSurfNormals[cellI];
1420 nNormals += normals.
size();
1423 reduce(nSet, sumOp<label>());
1424 reduce(nNormals, sumOp<label>());
1425 Info<<
"markSurfaceCurvatureRefinement :" 1426 <<
" cells:" << mesh_.globalData().nTotalCells()
1427 <<
" of which with normals:" << nSet
1428 <<
" ; total normals stored:" << nNormals
1434 bool reachedLimit =
false;
1443 !reachedLimit && cellI < cellSurfLocations.size();
1448 const vectorList& normals = cellSurfNormals[cellI];
1449 const labelList& levels = cellSurfLevels[cellI];
1452 const scalar cellSize =
1453 meshCutter_.level0EdgeLength()/
pow(2.0, cellLevel[cellI]);
1456 for (label i = 0; !reachedLimit && i < normals.size(); i++)
1458 for (label j = i+1; !reachedLimit && j < normals.size(); j++)
1469 (hasCurvatureLevels ? cellSize : scalar(0)),
1477 const label maxLevel =
max(levels[i], levels[j]);
1479 if (cellLevel[cellI] < maxLevel)
1494 Pout<<
"Stopped refining since reaching my cell" 1495 <<
" limit of " << mesh_.nCells()+7*nRefine
1498 reachedLimit =
true;
1518 !reachedLimit && faceI < mesh_.nInternalFaces();
1522 label own = mesh_.faceOwner()[faceI];
1523 label nei = mesh_.faceNeighbour()[faceI];
1525 const pointList& ownPoints = cellSurfLocations[own];
1526 const vectorList& ownNormals = cellSurfNormals[own];
1527 const labelList& ownLevels = cellSurfLevels[own];
1528 const pointList& neiPoints = cellSurfLocations[nei];
1529 const vectorList& neiNormals = cellSurfNormals[nei];
1530 const labelList& neiLevels = cellSurfLevels[nei];
1538 countMatches(ownNormals, neiNormals)
1539 == ownNormals.
size();
1542 countMatches(neiNormals, ownNormals)
1543 == neiNormals.size();
1546 if (!ownIsSubset && !neiIsSubset)
1549 const scalar cellSize =
1550 meshCutter_.level0EdgeLength()
1551 /
pow(2.0,
min(cellLevel[own], cellLevel[nei]));
1554 for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1556 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1565 (hasCurvatureLevels ? cellSize : scalar(0)),
1574 if (cellLevel[own] < ownLevels[i])
1589 Pout<<
"Stopped refining since reaching" 1590 <<
" my cell limit of " 1591 << mesh_.nCells()+7*nRefine <<
endl;
1593 reachedLimit =
true;
1597 if (cellLevel[nei] < neiLevels[j])
1612 Pout<<
"Stopped refining since reaching" 1613 <<
" my cell limit of " 1614 << mesh_.nCells()+7*nRefine <<
endl;
1616 reachedLimit =
true;
1630 List<vectorList> neiSurfacePoints;
1632 List<vectorList> neiSurfaceNormals;
1638 label faceI = mesh_.nInternalFaces();
1639 !reachedLimit && faceI < mesh_.nFaces();
1643 label own = mesh_.faceOwner()[faceI];
1644 label bFaceI = faceI - mesh_.nInternalFaces();
1646 const pointList& ownPoints = cellSurfLocations[own];
1647 const vectorList& ownNormals = cellSurfNormals[own];
1648 const labelList& ownLevels = cellSurfLevels[own];
1650 const pointList& neiPoints = neiSurfacePoints[bFaceI];
1651 const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
1659 countMatches(ownNormals, neiNormals)
1660 == ownNormals.
size();
1663 countMatches(neiNormals, ownNormals)
1664 == neiNormals.size();
1667 if (!ownIsSubset && !neiIsSubset)
1670 const scalar cellSize =
1671 meshCutter_.level0EdgeLength()
1672 /
pow(2.0,
min(cellLevel[own], neiLevel[bFaceI]));
1675 for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1677 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1686 (hasCurvatureLevels ? cellSize : scalar(0)),
1694 if (cellLevel[own] < ownLevels[i])
1709 Pout<<
"Stopped refining since reaching" 1710 <<
" my cell limit of " 1711 << mesh_.nCells()+7*nRefine
1714 reachedLimit =
true;
1731 Info<<
"Reached refinement limit." <<
endl;
1734 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1740 const scalar planarCos,
1750 vector d = point1-point0;
1751 scalar magD =
mag(d);
1753 if (magD > mergeDistance())
1755 scalar cosAngle = (normal0 & normal1);
1758 if (cosAngle < (-1+planarCos))
1761 avg = 0.5*(normal0-normal1);
1763 else if (cosAngle > (1-planarCos))
1765 avg = 0.5*(normal0+normal1);
1768 if (avg != vector::zero)
1773 if (
mag(avg&d) > mergeDistance())
1787 const scalar planarCos,
1801 vector d = point1-point0;
1802 scalar magD =
mag(d);
1804 if (magD > mergeDistance())
1806 scalar cosAngle = (normal0 & normal1);
1809 if (cosAngle < (-1+planarCos))
1812 avgNormal = 0.5*(normal0-normal1);
1814 else if (cosAngle > (1-planarCos))
1816 avgNormal = 0.5*(normal0+normal1);
1819 if (avgNormal != vector::zero)
1821 avgNormal /=
mag(avgNormal);
1843 bool Foam::meshRefinement::checkProximity
1845 const scalar planarCos,
1846 const label nAllowRefine,
1848 const label surfaceLevel,
1849 const vector& surfaceLocation,
1850 const vector& surfaceNormal,
1854 label& cellMaxLevel,
1862 const labelList& cellLevel = meshCutter_.cellLevel();
1865 if (surfaceLevel > cellLevel[cellI])
1867 if (cellMaxLevel == -1)
1870 cellMaxLevel = surfaceLevel;
1871 cellMaxLocation = surfaceLocation;
1872 cellMaxNormal = surfaceNormal;
1881 bool closeSurfaces = isNormalGap
1894 if (surfaceLevel > cellMaxLevel)
1896 cellMaxLevel = surfaceLevel;
1897 cellMaxLocation = surfaceLocation;
1898 cellMaxNormal = surfaceNormal;
1912 return markForRefine
1928 Foam::label Foam::meshRefinement::markProximityRefinement
1930 const scalar planarCos,
1935 const label nAllowRefine,
1943 const labelList& cellLevel = meshCutter_.cellLevel();
1945 label oldNRefine = nRefine;
1956 const labelList testFaces(getRefineCandidateFaces(refineCell));
1977 labelList cellMaxLevel(mesh_.nCells(), -1);
1983 List<vectorList> surfaceLocation;
1984 List<vectorList> surfaceNormal;
1988 surfaces_.findAllIntersections
2012 label faceI = testFaces[i];
2013 label own = mesh_.faceOwner()[faceI];
2015 const labelList& fLevels = surfaceLevel[i];
2016 const vectorList& fPoints = surfaceLocation[i];
2017 const vectorList& fNormals = surfaceNormal[i];
2032 cellMaxLocation[own],
2040 if (mesh_.isInternalFace(faceI))
2042 label nei = mesh_.faceNeighbour()[faceI];
2057 cellMaxLocation[nei],
2071 labelList neiBndMaxLevel(mesh_.nBoundaryFaces());
2072 pointField neiBndMaxLocation(mesh_.nBoundaryFaces());
2073 vectorField neiBndMaxNormal(mesh_.nBoundaryFaces());
2075 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2077 label bFaceI = faceI-mesh_.nInternalFaces();
2078 label own = mesh_.faceOwner()[faceI];
2080 neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
2081 neiBndMaxLocation[bFaceI] = cellMaxLocation[own];
2082 neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
2091 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2093 label own = mesh_.faceOwner()[faceI];
2094 label nei = mesh_.faceNeighbour()[faceI];
2096 if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
2105 cellMaxLocation[own],
2108 cellMaxLocation[nei],
2114 if (cellLevel[own] < cellMaxLevel[own])
2129 Pout<<
"Stopped refining since reaching my cell" 2130 <<
" limit of " << mesh_.nCells()+7*nRefine
2137 if (cellLevel[nei] < cellMaxLevel[nei])
2152 Pout<<
"Stopped refining since reaching my cell" 2153 <<
" limit of " << mesh_.nCells()+7*nRefine
2163 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2165 label own = mesh_.faceOwner()[faceI];
2166 label bFaceI = faceI - mesh_.nInternalFaces();
2168 if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
2177 cellMaxLocation[own],
2180 neiBndMaxLocation[bFaceI],
2181 neiBndMaxNormal[bFaceI]
2198 Pout<<
"Stopped refining since reaching my cell" 2199 <<
" limit of " << mesh_.nCells()+7*nRefine
2214 Info<<
"Reached refinement limit." <<
endl;
2217 return returnReduce(nRefine-oldNRefine, sumOp<label>());
2231 const scalar curvature,
2232 const scalar planarAngle,
2234 const bool featureRefinement,
2235 const bool featureDistanceRefinement,
2236 const bool internalRefinement,
2237 const bool surfaceRefinement,
2238 const bool curvatureRefinement,
2239 const bool smallFeatureRefinement,
2240 const bool gapRefinement,
2241 const bool bigGapRefinement,
2242 const bool spreadGapSize,
2243 const label maxGlobalCells,
2244 const label maxLocalCells
2247 label totNCells = mesh_.globalData().nTotalCells();
2251 if (totNCells >= maxGlobalCells)
2253 Info<<
"No cells marked for refinement since reached limit " 2254 << maxGlobalCells <<
'.' <<
endl;
2288 labelList neiLevel(mesh_.nBoundaryFaces());
2290 calcNeighbourData(neiLevel, neiCc);
2299 if (featureRefinement)
2301 label nFeatures = markFeatureRefinement
2310 Info<<
"Marked for refinement due to explicit features " 2311 <<
": " << nFeatures <<
" cells." <<
endl;
2317 if (featureDistanceRefinement)
2319 label nShell = markInternalDistanceToFeatureRefinement
2326 Info<<
"Marked for refinement due to distance to explicit features " 2327 ": " << nShell <<
" cells." <<
endl;
2333 if (internalRefinement)
2335 label nShell = markInternalRefinement
2342 Info<<
"Marked for refinement due to refinement shells " 2343 <<
": " << nShell <<
" cells." <<
endl;
2349 if (surfaceRefinement)
2351 label nSurf = markSurfaceRefinement
2360 Info<<
"Marked for refinement due to surface intersection " 2361 <<
": " << nSurf <<
" cells." <<
endl;
2369 &&
max(shells_.maxGapLevel()) > 0
2372 label nGapSurf = markSurfaceGapRefinement
2382 Info<<
"Marked for refinement due to surface intersection" 2384 <<
": " << nGapSurf <<
" cells." <<
endl;
2394 && (curvature >= -1 && curvature <= 1)
2395 && (surfaces_.minLevel() != surfaces_.maxLevel())
2398 label nCurv = markSurfaceCurvatureRefinement
2408 Info<<
"Marked for refinement due to curvature/regions " 2409 <<
": " << nCurv <<
" cells." <<
endl;
2416 if (smallFeatureRefinement)
2423 &&
max(shells_.maxGapLevel()) > 0
2426 nGap = markSmallFeatureRefinement
2437 Info<<
"Marked for refinement due to close opposite surfaces " 2438 <<
": " << nGap <<
" cells." <<
endl;
2441 if (
max(surfaces_.maxCurvatureLevel()) > 0)
2443 nCurv = markSurfaceFieldRefinement
2452 Info<<
"Marked for refinement" 2453 <<
" due to curvature " 2454 <<
": " << nCurv <<
" cells." <<
endl;
2462 const labelList& surfaceGapLevel = surfaces_.gapLevel();
2467 && (planarCos >= -1 && planarCos <= 1)
2468 && (
max(surfaceGapLevel) > -1)
2471 Info<<
"Specified gap level : " <<
max(surfaceGapLevel)
2472 <<
", planar angle " << planarAngle <<
endl;
2474 label nGap = markProximityRefinement
2488 Info<<
"Marked for refinement due to close opposite surfaces " 2489 <<
": " << nGap <<
" cells." <<
endl;
2499 && (planarCos >= -1 && planarCos <= 1)
2500 &&
max(shells_.maxGapLevel()) > 0
2507 label nGap = markInternalGapRefinement
2518 Info<<
"Marked for refinement due to opposite surfaces" 2520 <<
": " << nGap <<
" cells." <<
endl;
2527 if (limitShells_.shells().size())
2529 label nUnmarked = unmarkInternalRefinement(
refineCell, nRefine);
2532 Info<<
"Unmarked for refinement due to limit shells" 2533 <<
" : " << nUnmarked <<
" cells." <<
endl;
2542 cellsToRefine.
setSize(nRefine);
2549 cellsToRefine[nRefine++] = cellI;
2554 return cellsToRefine;
2567 meshCutter_.setRefinement(cellsToRefine, meshMod);
2571 mesh_.moving(
false);
2578 mesh_.updateMesh(map);
2595 updateMesh(map, getChangedFaces(map, cellsToRefine));
2607 decompositionMethod& decomposer,
2608 fvMeshDistribute& distributor,
2610 const scalar maxLoadUnbalance
2614 refine(cellsToRefine);
2618 Pout<<
"Writing refined but unbalanced " << msg
2626 Pout<<
"Dumped debug data in = " 2627 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2633 Info<<
"Refined mesh in = " 2634 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2635 printMeshInfo(
debug,
"After refinement " + msg);
2645 scalar nIdealCells =
2646 mesh_.globalData().nTotalCells()
2651 mag(1.0-mesh_.nCells()/nIdealCells),
2655 if (unbalance <= maxLoadUnbalance)
2657 Info<<
"Skipping balancing since max unbalance " << unbalance
2658 <<
" is less than allowable " << maxLoadUnbalance
2674 Info<<
"Balanced mesh in = " 2675 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2677 printMeshInfo(
debug,
"After balancing " + msg);
2682 Pout<<
"Writing balanced " << msg
2687 writeType(writeLevel() | WRITEMESH),
2690 Pout<<
"Dumped debug data in = " 2691 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2708 decompositionMethod& decomposer,
2709 fvMeshDistribute& distributor,
2711 const scalar maxLoadUnbalance
2714 labelList cellsToRefine(initCellsToRefine);
2749 scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
2750 scalar nIdealNewCells =
2755 mag(1.0-nNewCells/nIdealNewCells),
2759 if (unbalance <= maxLoadUnbalance)
2761 Info<<
"Skipping balancing since max unbalance " << unbalance
2762 <<
" is less than allowable " << maxLoadUnbalance
2770 cellWeights[cellsToRefine[i]] += 7;
2783 distMap().distributeCellIndices(cellsToRefine);
2785 Info<<
"Balanced mesh in = " 2786 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2801 printMeshInfo(
debug,
"After balancing " + msg);
2805 Pout<<
"Writing balanced " << msg
2813 Pout<<
"Dumped debug data in = " 2814 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2825 refine(cellsToRefine);
2829 Pout<<
"Writing refined " << msg
2834 writeType(writeLevel() | WRITEMESH),
2837 Pout<<
"Dumped debug data in = " 2838 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2844 Info<<
"Refined mesh in = " 2845 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2859 printMeshInfo(
debug,
"After refinement " + msg);
2867 const label maxGlobalCells,
2868 const label maxLocalCells,
2873 const labelList& cellLevel = meshCutter_.cellLevel();
2874 const pointField& cellCentres = mesh_.cellCentres();
2876 label totNCells = mesh_.globalData().nTotalCells();
2880 if (totNCells >= maxGlobalCells)
2882 Info<<
"No cells marked for refinement since reached limit " 2883 << maxGlobalCells <<
'.' <<
endl;
2897 shells_.findDirectionalLevel
2907 forAll(insideShell, celli)
2909 if (insideShell[celli] >= 0)
2911 bool reachedLimit = !markForRefine
2923 Pout<<
"Stopped refining cells" 2924 <<
" since reaching my cell limit of " 2925 << mesh_.nCells()+nAllowRefine <<
endl;
2936 label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
2939 Info<<
"Unmarked for refinement due to limit shells" 2940 <<
" : " << nUnmarked <<
" cells." <<
endl;
2949 cellsToRefine.setSize(nRefine);
2952 forAll(refineCell, cellI)
2954 if (refineCell[cellI] != -1)
2956 cellsToRefine[nRefine++] = cellI;
2961 return cellsToRefine;
2975 List<refineCell> refCells(cellsToRefine.size());
2978 refCells[i] = refineCell(cellsToRefine[i], refDir);
2982 hexCellLooper cellWalker(mesh_);
2985 cellCuts cuts(mesh_, cellWalker, refCells);
2993 meshRefiner.setRefinement(cuts, meshMod);
2997 mesh_.moving(
false);
3002 mesh_.updateMesh(*morphMap);
3005 if (morphMap().hasMotionPoints())
3007 mesh_.movePoints(morphMap().preMotionPoints());
3019 meshRefiner.updateMesh(*morphMap);
3022 updateMesh(*morphMap, getChangedFaces(*morphMap, cellsToRefine));
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
List< labelList > labelListList
A List of labelList.
void size(const label n)
Older name for setAddressableSize.
label nOldCells() const
Number of old cells.
const labelList & cellMap() const
Old cell map.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
autoPtr< mapDistributePolyMesh > balanceAndRefine(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance)
Balance before refining some cells.
virtual const labelList & faceNeighbour() const
Return face neighbour.
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
const word & name() const noexcept
Return the object name.
bool isNormalGap(const scalar planarCos, const label level0, const vector &point0, const vector &normal0, const label level1, const vector &point1, const vector &normal1) const
Is local topology a small gap normal to the test vector.
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
A traits class, which is primarily used for primitives.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
const cellList & cells() const
labelList refineCandidates(const pointField &keepPoints, const scalar curvature, const scalar planarAngle, const bool featureRefinement, const bool featureDistanceRefinement, const bool internalRefinement, const bool surfaceRefinement, const bool curvatureRefinement, const bool smallFeatureRefinement, const bool gapRefinement, const bool bigGapRefinement, const bool spreadGapSize, const label maxGlobalCells, const label maxLocalCells) const
Calculate list of cells to refine.
List< point > pointList
A List of points.
label nFaces() const noexcept
Number of mesh faces.
autoPtr< mapPolyMesh > directionalRefine(const string &msg, const direction cmpt, const labelList &cellsToRefine)
Directionally refine in direction cmpt.
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.
Description of cuts across cells.
static bool less(const vector &x, const vector &y)
To compare normals.
bool hasMotionPoints() const
Has valid preMotionPoints?
List< vector > vectorList
A List of vectors.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const fvMesh & mesh() const
Reference to mesh.
#define forAll(list, i)
Loop across all elements in list.
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator) is 1 for serial run.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
dimensionedScalar cos(const dimensionedScalar &ds)
autoPtr< mapDistributePolyMesh > refineAndBalance(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance)
Refine some cells and rebalance.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void clear()
Clear the list, i.e. set size to zero.
label checkData(const fvMesh &mesh, const instantList &timeDirs, wordList &objectNames)
Check if fields are good to use (available at all times)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelList & faceOwner() const
Return face owner.
const globalMeshData & globalData() const
Return parallel info.
label nInternalFaces() const noexcept
Number of internal faces.
Container with cells to refine. Refinement given as single direction.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
labelList directionalRefineCandidates(const label maxGlobalCells, const label maxLocalCells, const labelList ¤tLevel, const direction dir) const
Calculate list of cells to directionally refine.
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
autoPtr< mapPolyMesh > refine(const labelList &cellsToRefine)
Refine some cells.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
UIndirectList< label > labelUIndList
UIndirectList of labels.
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
A cell is defined as a list of faces with extra functionality.
Direct mesh changes based on v1.3 polyTopoChange syntax.
const polyMesh & mesh() const
Return polyMesh.
label start() const
Return start label of this patch in the polyMesh face list.
normalLess(const vectorList &values)
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.
messageStream Info
Information stream (stdout output on master, null elsewhere)
label nTotalFaces() const noexcept
Return total number of faces in decomposed mesh. Not.
writeType
Enumeration for what to write. Used as a bit-pattern.
Field< vector > vectorField
Specialisation of Field<T> for vector.
debugType
Enumeration for what to debug. Used as a bit-pattern.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
Mesh consisting of general polyhedral cells.
bool isGap(const scalar, const vector &, const vector &, const vector &, const vector &) const
Is local topology a small gap?
List< label > labelList
A List of labels.
A patch is a list of labels that address the faces in the global face list.
virtual bool write(const bool valid=true) const
Write using setting from DB.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
const pointField & preMotionPoints() const
Pre-motion point positions.
bool operator()(const label a, const label b) const
static constexpr const zero Zero
Global zero (0)