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
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;
496 Cloud<trackedParticle> cloud(mesh_,
Foam::zero{},
"featureCloud");
500 Pout<<
"Constructing cloud for cell marking" <<
endl;
503 for (
const trackedParticle& startTp : startPointCloud)
505 const label featI = startTp.i();
506 const label pointI = startTp.j();
508 const edgeMesh& featureMesh = features_[featI];
509 const labelList& pEdges = featureMesh.pointEdges()[pointI];
514 label edgeI = pEdges[pEdgeI];
516 if (featureEdgeVisited[featI].
set(edgeI))
521 const edge&
e = featureMesh.edges()[edgeI];
522 label otherPointi =
e.otherVertex(pointI);
524 trackedParticle* tp(
new trackedParticle(startTp));
525 tp->start() = tp->position();
526 tp->end() = featureMesh.points()[otherPointi];
527 tp->j() = otherPointi;
532 Pout<<
"Adding particle for point:" << pointI
533 <<
" coord:" << tp->position()
534 <<
" feature:" << featI
535 <<
" to track to:" << tp->end()
539 cloud.addParticle(tp);
544 startPointCloud.clear();
552 Pout<<
"Tracking " << cloud.size()
553 <<
" particles over distance " << maxTrackLen
554 <<
" to mark cells" <<
endl;
556 cloud.move(cloud,
td, maxTrackLen);
559 for (trackedParticle& tp : cloud)
561 const label featI = tp.i();
562 const label pointI = tp.j();
564 const edgeMesh& featureMesh = features_[featI];
565 const labelList& pEdges = featureMesh.pointEdges()[pointI];
570 bool keepParticle =
false;
574 label edgeI = pEdges[i];
576 if (featureEdgeVisited[featI].
set(edgeI))
581 const edge&
e = featureMesh.edges()[edgeI];
582 label otherPointi =
e.otherVertex(pointI);
584 tp.start() = tp.position();
585 tp.end() = featureMesh.points()[otherPointi];
586 tp.j() = otherPointi;
597 cloud.deleteParticle(tp);
604 Pout<<
"Remaining particles " << cloud.size() <<
endl;
632 Foam::label Foam::meshRefinement::markFeatureRefinement
635 const label nAllowRefine,
643 markFeatureCellLevel(keepPoints, maxFeatureLevel);
648 const labelList& cellLevel = meshCutter_.cellLevel();
652 forAll(maxFeatureLevel, cellI)
654 if (maxFeatureLevel[cellI] > cellLevel[cellI])
680 Info<<
"Reached refinement limit." <<
endl;
688 Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
690 const label nAllowRefine,
696 const labelList& cellLevel = meshCutter_.cellLevel();
697 const pointField& cellCentres = mesh_.cellCentres();
700 if (features_.maxDistance() <= 0.0)
709 labelList testLevels(cellLevel.size()-nRefine);
714 if (refineCell[cellI] == -1)
716 testCc[testI] = cellCentres[cellI];
717 testLevels[testI] = cellLevel[cellI];
724 features_.findHigherLevel(testCc, testLevels, maxLevel);
731 if (refineCell[cellI] == -1)
733 if (maxLevel[testI] > testLevels[testI])
735 bool reachedLimit = !markForRefine
747 Pout<<
"Stopped refining internal cells" 748 <<
" since reaching my cell limit of " 749 << mesh_.nCells()+7*nRefine <<
endl;
764 Info<<
"Reached refinement limit." <<
endl;
772 Foam::label Foam::meshRefinement::markInternalRefinement
774 const label nAllowRefine,
780 const labelList& cellLevel = meshCutter_.cellLevel();
781 const pointField& cellCentres = mesh_.cellCentres();
787 labelList testLevels(cellLevel.size()-nRefine);
792 if (refineCell[cellI] == -1)
794 testCc[testI] = cellCentres[cellI];
795 testLevels[testI] = cellLevel[cellI];
802 shells_.findHigherLevel(testCc, testLevels, maxLevel);
809 if (refineCell[cellI] == -1)
811 if (maxLevel[testI] > testLevels[testI])
813 bool reachedLimit = !markForRefine
825 Pout<<
"Stopped refining internal cells" 826 <<
" since reaching my cell limit of " 827 << mesh_.nCells()+7*nRefine <<
endl;
842 Info<<
"Reached refinement limit." <<
endl;
849 Foam::label Foam::meshRefinement::unmarkInternalRefinement
855 const labelList& cellLevel = meshCutter_.cellLevel();
856 const pointField& cellCentres = mesh_.cellCentres();
867 if (refineCell[cellI] >= 0)
869 testCc[testI] = cellCentres[cellI];
870 testLevels[testI] = cellLevel[cellI];
877 limitShells_.findLevel(testCc, testLevels, levelShell);
884 if (refineCell[cellI] >= 0)
886 if (levelShell[testI] != -1)
888 refineCell[cellI] = -1;
910 const labelList& surfIndex = surfaceIndex();
917 if (surfIndex[faceI] != -1)
919 label own = mesh_.faceOwner()[faceI];
921 if (mesh_.isInternalFace(faceI))
923 label nei = mesh_.faceNeighbour()[faceI];
925 if (refineCell[own] == -1 || refineCell[nei] == -1)
927 testFaces[nTest++] = faceI;
932 const label bFacei = faceI - mesh_.nInternalFaces();
933 if (refineCell[own] == -1 || boundaryRefineCell[bFacei] == -1)
935 testFaces[nTest++] = faceI;
940 testFaces.setSize(nTest);
947 Foam::label Foam::meshRefinement::markSurfaceRefinement
949 const label nAllowRefine,
957 const labelList& cellLevel = meshCutter_.cellLevel();
967 labelList testFaces(getRefineCandidateFaces(refineCell));
992 surfaces_.findHigherIntersection
1009 label faceI = testFaces[i];
1011 label surfI = surfaceHit[i];
1021 label own = mesh_.faceOwner()[faceI];
1023 if (surfaceMinLevel[i] > cellLevel[own])
1041 if (mesh_.isInternalFace(faceI))
1043 label nei = mesh_.faceNeighbour()[faceI];
1044 if (surfaceMinLevel[i] > cellLevel[nei])
1071 Info<<
"Reached refinement limit." <<
endl;
1079 Foam::label Foam::meshRefinement::countMatches
1081 const List<point>& normals1,
1082 const List<point>& normals2,
1090 const vector& n1 = normals1[i];
1094 const vector& n2 = normals2[j];
1119 bool Foam::meshRefinement::highCurvature
1121 const scalar minCosAngle,
1122 const scalar lengthScale,
1129 const scalar cosAngle = (n0&n1);
1131 if (cosAngle < minCosAngle)
1136 else if (cosAngle > 1-1
e-6)
1141 else if (lengthScale > SMALL)
1146 const plane pl0(
p0, (n0 ^ axis));
1147 const scalar r1 = pl0.normalIntersect(p1, n1);
1151 const plane pl1(p1, (n1 ^ axis));
1152 const scalar r0 = pl1.normalIntersect(
p0, n0);
1157 const scalar radius = 0.5*(
mag(r1)+
mag(r0));
1159 if (radius < lengthScale)
1177 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
1179 const scalar curvature,
1180 const label nAllowRefine,
1188 const labelList& cellLevel = meshCutter_.cellLevel();
1189 const pointField& cellCentres = mesh_.cellCentres();
1205 labelList testFaces(getRefineCandidateFaces(refineCell));
1215 label faceI = testFaces[i];
1217 label own = mesh_.faceOwner()[faceI];
1219 if (mesh_.isInternalFace(faceI))
1221 label nei = mesh_.faceNeighbour()[faceI];
1223 start[i] = cellCentres[own];
1224 end[i] = cellCentres[nei];
1225 minLevel[i] =
min(cellLevel[own], cellLevel[nei]);
1229 label bFaceI = faceI - mesh_.nInternalFaces();
1231 start[i] = cellCentres[own];
1232 end[i] = neiCc[bFaceI];
1234 if (!isMasterFace[faceI])
1236 std::swap(start[i],
end[i]);
1239 minLevel[i] =
min(cellLevel[own], neiLevel[bFaceI]);
1252 const bool hasCurvatureLevels = (
max(surfaces_.maxCurvatureLevel()) > 0);
1258 List<pointList> cellSurfLocations(mesh_.nCells());
1259 List<vectorList> cellSurfNormals(mesh_.nCells());
1263 List<pointList> surfaceLocation;
1265 List<vectorList> surfaceNormal;
1269 surfaces_.findAllIntersections
1276 labelList(surfaces_.maxLevel().size(), 0),
1277 surfaces_.maxLevel(),
1289 forAll(surfaceNormal, pointI)
1291 pointList& pLocations = surfaceLocation[pointI];
1292 vectorList& pNormals = surfaceNormal[pointI];
1293 labelList& pLevel = surfaceLevel[pointI];
1295 sortedOrder(pNormals, visitOrder, normalLess(pLocations));
1297 pLocations = List<point>(pLocations, visitOrder);
1298 pNormals = List<vector>(pNormals, visitOrder);
1373 label faceI = testFaces[i];
1374 label own = mesh_.faceOwner()[faceI];
1376 const pointList& fPoints = surfaceLocation[i];
1377 const vectorList& fNormals = surfaceNormal[i];
1378 const labelList& fLevels = surfaceLevel[i];
1382 if (fLevels[hitI] > cellLevel[own])
1384 cellSurfLevels[own].append(fLevels[hitI]);
1385 cellSurfLocations[own].append(fPoints[hitI]);
1386 cellSurfNormals[own].append(fNormals[hitI]);
1389 if (mesh_.isInternalFace(faceI))
1391 label nei = mesh_.faceNeighbour()[faceI];
1392 if (fLevels[hitI] > cellLevel[nei])
1394 cellSurfLevels[nei].append(fLevels[hitI]);
1395 cellSurfLocations[nei].append(fPoints[hitI]);
1396 cellSurfNormals[nei].append(fNormals[hitI]);
1410 forAll(cellSurfNormals, cellI)
1412 const vectorList& normals = cellSurfNormals[cellI];
1416 nNormals += normals.
size();
1419 reduce(nSet, sumOp<label>());
1420 reduce(nNormals, sumOp<label>());
1421 Info<<
"markSurfaceCurvatureRefinement :" 1422 <<
" cells:" << mesh_.globalData().nTotalCells()
1423 <<
" of which with normals:" << nSet
1424 <<
" ; total normals stored:" << nNormals
1430 bool reachedLimit =
false;
1439 !reachedLimit && cellI < cellSurfLocations.size();
1444 const vectorList& normals = cellSurfNormals[cellI];
1445 const labelList& levels = cellSurfLevels[cellI];
1448 const scalar cellSize =
1449 meshCutter_.level0EdgeLength()/
pow(2.0, cellLevel[cellI]);
1452 for (label i = 0; !reachedLimit && i < normals.size(); i++)
1454 for (label j = i+1; !reachedLimit && j < normals.size(); j++)
1465 (hasCurvatureLevels ? cellSize : scalar(0)),
1473 const label maxLevel =
max(levels[i], levels[j]);
1475 if (cellLevel[cellI] < maxLevel)
1490 Pout<<
"Stopped refining since reaching my cell" 1491 <<
" limit of " << mesh_.nCells()+7*nRefine
1494 reachedLimit =
true;
1514 !reachedLimit && faceI < mesh_.nInternalFaces();
1518 label own = mesh_.faceOwner()[faceI];
1519 label nei = mesh_.faceNeighbour()[faceI];
1521 const pointList& ownPoints = cellSurfLocations[own];
1522 const vectorList& ownNormals = cellSurfNormals[own];
1523 const labelList& ownLevels = cellSurfLevels[own];
1524 const pointList& neiPoints = cellSurfLocations[nei];
1525 const vectorList& neiNormals = cellSurfNormals[nei];
1526 const labelList& neiLevels = cellSurfLevels[nei];
1534 countMatches(ownNormals, neiNormals)
1535 == ownNormals.
size();
1538 countMatches(neiNormals, ownNormals)
1539 == neiNormals.size();
1542 if (!ownIsSubset && !neiIsSubset)
1545 const scalar cellSize =
1546 meshCutter_.level0EdgeLength()
1547 /
pow(2.0,
min(cellLevel[own], cellLevel[nei]));
1550 for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1552 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1561 (hasCurvatureLevels ? cellSize : scalar(0)),
1570 if (cellLevel[own] < ownLevels[i])
1585 Pout<<
"Stopped refining since reaching" 1586 <<
" my cell limit of " 1587 << mesh_.nCells()+7*nRefine <<
endl;
1589 reachedLimit =
true;
1593 if (cellLevel[nei] < neiLevels[j])
1608 Pout<<
"Stopped refining since reaching" 1609 <<
" my cell limit of " 1610 << mesh_.nCells()+7*nRefine <<
endl;
1612 reachedLimit =
true;
1626 List<vectorList> neiSurfacePoints;
1628 List<vectorList> neiSurfaceNormals;
1634 label faceI = mesh_.nInternalFaces();
1635 !reachedLimit && faceI < mesh_.nFaces();
1639 label own = mesh_.faceOwner()[faceI];
1640 label bFaceI = faceI - mesh_.nInternalFaces();
1642 const pointList& ownPoints = cellSurfLocations[own];
1643 const vectorList& ownNormals = cellSurfNormals[own];
1644 const labelList& ownLevels = cellSurfLevels[own];
1646 const pointList& neiPoints = neiSurfacePoints[bFaceI];
1647 const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
1655 countMatches(ownNormals, neiNormals)
1656 == ownNormals.
size();
1659 countMatches(neiNormals, ownNormals)
1660 == neiNormals.size();
1663 if (!ownIsSubset && !neiIsSubset)
1666 const scalar cellSize =
1667 meshCutter_.level0EdgeLength()
1668 /
pow(2.0,
min(cellLevel[own], neiLevel[bFaceI]));
1671 for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1673 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1682 (hasCurvatureLevels ? cellSize : scalar(0)),
1690 if (cellLevel[own] < ownLevels[i])
1705 Pout<<
"Stopped refining since reaching" 1706 <<
" my cell limit of " 1707 << mesh_.nCells()+7*nRefine
1710 reachedLimit =
true;
1727 Info<<
"Reached refinement limit." <<
endl;
1736 const scalar planarCos,
1746 vector d = point1-point0;
1747 scalar magD =
mag(d);
1749 if (magD > mergeDistance())
1751 scalar cosAngle = (normal0 & normal1);
1754 if (cosAngle < (-1+planarCos))
1757 avg = 0.5*(normal0-normal1);
1759 else if (cosAngle > (1-planarCos))
1761 avg = 0.5*(normal0+normal1);
1764 if (avg != vector::zero)
1769 if (
mag(avg&d) > mergeDistance())
1783 const scalar planarCos,
1797 vector d = point1-point0;
1798 scalar magD =
mag(d);
1800 if (magD > mergeDistance())
1802 scalar cosAngle = (normal0 & normal1);
1805 if (cosAngle < (-1+planarCos))
1808 avgNormal = 0.5*(normal0-normal1);
1810 else if (cosAngle > (1-planarCos))
1812 avgNormal = 0.5*(normal0+normal1);
1815 if (avgNormal != vector::zero)
1817 avgNormal /=
mag(avgNormal);
1839 bool Foam::meshRefinement::checkProximity
1841 const scalar planarCos,
1842 const label nAllowRefine,
1844 const label surfaceLevel,
1845 const vector& surfaceLocation,
1846 const vector& surfaceNormal,
1850 label& cellMaxLevel,
1858 const labelList& cellLevel = meshCutter_.cellLevel();
1861 if (surfaceLevel > cellLevel[cellI])
1863 if (cellMaxLevel == -1)
1866 cellMaxLevel = surfaceLevel;
1867 cellMaxLocation = surfaceLocation;
1868 cellMaxNormal = surfaceNormal;
1877 bool closeSurfaces = isNormalGap
1890 if (surfaceLevel > cellMaxLevel)
1892 cellMaxLevel = surfaceLevel;
1893 cellMaxLocation = surfaceLocation;
1894 cellMaxNormal = surfaceNormal;
1908 return markForRefine
1924 Foam::label Foam::meshRefinement::markProximityRefinement
1926 const scalar planarCos,
1931 const label nAllowRefine,
1939 const labelList& cellLevel = meshCutter_.cellLevel();
1952 const labelList testFaces(getRefineCandidateFaces(refineCell));
1973 labelList cellMaxLevel(mesh_.nCells(), -1);
1979 List<vectorList> surfaceLocation;
1980 List<vectorList> surfaceNormal;
1984 surfaces_.findAllIntersections
2008 label faceI = testFaces[i];
2009 label own = mesh_.faceOwner()[faceI];
2011 const labelList& fLevels = surfaceLevel[i];
2012 const vectorList& fPoints = surfaceLocation[i];
2013 const vectorList& fNormals = surfaceNormal[i];
2028 cellMaxLocation[own],
2036 if (mesh_.isInternalFace(faceI))
2038 label nei = mesh_.faceNeighbour()[faceI];
2053 cellMaxLocation[nei],
2067 labelList neiBndMaxLevel(mesh_.nBoundaryFaces());
2068 pointField neiBndMaxLocation(mesh_.nBoundaryFaces());
2069 vectorField neiBndMaxNormal(mesh_.nBoundaryFaces());
2071 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2073 label bFaceI = faceI-mesh_.nInternalFaces();
2074 label own = mesh_.faceOwner()[faceI];
2076 neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
2077 neiBndMaxLocation[bFaceI] = cellMaxLocation[own];
2078 neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
2087 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2089 label own = mesh_.faceOwner()[faceI];
2090 label nei = mesh_.faceNeighbour()[faceI];
2092 if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
2101 cellMaxLocation[own],
2104 cellMaxLocation[nei],
2110 if (cellLevel[own] < cellMaxLevel[own])
2125 Pout<<
"Stopped refining since reaching my cell" 2126 <<
" limit of " << mesh_.nCells()+7*nRefine
2133 if (cellLevel[nei] < cellMaxLevel[nei])
2148 Pout<<
"Stopped refining since reaching my cell" 2149 <<
" limit of " << mesh_.nCells()+7*nRefine
2159 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
2161 label own = mesh_.faceOwner()[faceI];
2162 label bFaceI = faceI - mesh_.nInternalFaces();
2164 if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
2173 cellMaxLocation[own],
2176 neiBndMaxLocation[bFaceI],
2177 neiBndMaxNormal[bFaceI]
2194 Pout<<
"Stopped refining since reaching my cell" 2195 <<
" limit of " << mesh_.nCells()+7*nRefine
2210 Info<<
"Reached refinement limit." <<
endl;
2227 const scalar curvature,
2228 const scalar planarAngle,
2230 const bool featureRefinement,
2231 const bool featureDistanceRefinement,
2232 const bool internalRefinement,
2233 const bool surfaceRefinement,
2234 const bool curvatureRefinement,
2235 const bool smallFeatureRefinement,
2236 const bool gapRefinement,
2237 const bool bigGapRefinement,
2238 const bool spreadGapSize,
2239 const label maxGlobalCells,
2240 const label maxLocalCells
2243 label totNCells = mesh_.globalData().nTotalCells();
2247 if (totNCells >= maxGlobalCells)
2249 Info<<
"No cells marked for refinement since reached limit " 2250 << maxGlobalCells <<
'.' <<
endl;
2284 labelList neiLevel(mesh_.nBoundaryFaces());
2286 calcNeighbourData(neiLevel, neiCc);
2295 if (featureRefinement)
2297 label nFeatures = markFeatureRefinement
2306 Info<<
"Marked for refinement due to explicit features " 2307 <<
": " << nFeatures <<
" cells." <<
endl;
2313 if (featureDistanceRefinement)
2315 label nShell = markInternalDistanceToFeatureRefinement
2322 Info<<
"Marked for refinement due to distance to explicit features " 2323 ": " << nShell <<
" cells." <<
endl;
2329 if (internalRefinement)
2331 label nShell = markInternalRefinement
2338 Info<<
"Marked for refinement due to refinement shells " 2339 <<
": " << nShell <<
" cells." <<
endl;
2345 if (surfaceRefinement)
2347 label nSurf = markSurfaceRefinement
2356 Info<<
"Marked for refinement due to surface intersection " 2357 <<
": " << nSurf <<
" cells." <<
endl;
2365 &&
max(shells_.maxGapLevel()) > 0
2368 label nGapSurf = markSurfaceGapRefinement
2378 Info<<
"Marked for refinement due to surface intersection" 2380 <<
": " << nGapSurf <<
" cells." <<
endl;
2390 && (curvature >= -1 && curvature <= 1)
2391 && (surfaces_.minLevel() != surfaces_.maxLevel())
2394 label nCurv = markSurfaceCurvatureRefinement
2404 Info<<
"Marked for refinement due to curvature/regions " 2405 <<
": " << nCurv <<
" cells." <<
endl;
2412 if (smallFeatureRefinement)
2419 &&
max(shells_.maxGapLevel()) > 0
2422 nGap = markSmallFeatureRefinement
2433 Info<<
"Marked for refinement due to close opposite surfaces " 2434 <<
": " << nGap <<
" cells." <<
endl;
2437 if (
max(surfaces_.maxCurvatureLevel()) > 0)
2439 nCurv = markSurfaceFieldRefinement
2448 Info<<
"Marked for refinement" 2449 <<
" due to curvature " 2450 <<
": " << nCurv <<
" cells." <<
endl;
2458 const labelList& surfaceGapLevel = surfaces_.gapLevel();
2463 && (planarCos >= -1 && planarCos <= 1)
2464 && (
max(surfaceGapLevel) > -1)
2467 Info<<
"Specified gap level : " <<
max(surfaceGapLevel)
2468 <<
", planar angle " << planarAngle <<
endl;
2470 label nGap = markProximityRefinement
2484 Info<<
"Marked for refinement due to close opposite surfaces " 2485 <<
": " << nGap <<
" cells." <<
endl;
2495 && (planarCos >= -1 && planarCos <= 1)
2496 &&
max(shells_.maxGapLevel()) > 0
2503 label nGap = markInternalGapRefinement
2514 Info<<
"Marked for refinement due to opposite surfaces" 2516 <<
": " << nGap <<
" cells." <<
endl;
2523 if (limitShells_.shells().size())
2525 label nUnmarked = unmarkInternalRefinement(
refineCell, nRefine);
2528 Info<<
"Unmarked for refinement due to limit shells" 2529 <<
" : " << nUnmarked <<
" cells." <<
endl;
2538 cellsToRefine.
setSize(nRefine);
2545 cellsToRefine[nRefine++] = cellI;
2550 return cellsToRefine;
2563 meshCutter_.setRefinement(cellsToRefine, meshMod);
2567 mesh_.moving(
false);
2574 mesh_.updateMesh(map);
2591 updateMesh(map, getChangedFaces(map, cellsToRefine));
2601 decompositionMethod& decomposer,
2602 fvMeshDistribute& distributor,
2604 const scalar maxLoadUnbalance,
2605 const label maxCellUnbalance
2608 autoPtr<mapDistributePolyMesh> distMap;
2614 const scalar nNewCells =
2615 scalar(mesh_.nCells() + 7*cellsToRefine.size());
2616 const scalar nNewCellsAll =
returnReduce(nNewCells, sumOp<scalar>());
2620 mag(1.0-nNewCells/nIdealNewCells),
2626 const scalar nNewCellsOnly = scalar(7*cellsToRefine.size());
2628 const label maxNewCells =
2631 const label maxDeltaCells =
2641 (maxNewCells <= maxCellUnbalance)
2642 && (maxDeltaCells <= maxCellUnbalance)
2645 Info<<
"Skipping balancing since trigger value not reached:" <<
"\n" 2646 <<
" Trigger cell count: " << maxCellUnbalance <<
nl 2647 <<
" Max new cell count in proc: " << maxNewCells <<
nl 2648 <<
" Max difference between new cells and balanced: " 2649 << maxDeltaCells <<
nl 2650 <<
" Max load unbalance " << maxLoadUnbalance
2655 if (unbalance <= maxLoadUnbalance)
2657 Info<<
"Skipping balancing since max unbalance " << unbalance
2658 <<
" is less than allowable " << maxLoadUnbalance
2663 Info<<
"Balancing since max unbalance " << unbalance
2664 <<
" is larger than allowable " << maxLoadUnbalance
2670 cellWeights[cellsToRefine[i]] += 7;
2683 distMap().distributeCellIndices(cellsToRefine);
2685 Info<<
"Balanced mesh in = " 2686 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2688 printMeshInfo(
debug,
"After balancing " + msg,
true);
2693 Pout<<
"Writing balanced " << msg
2698 writeType(writeLevel() | WRITEMESH),
2701 Pout<<
"Dumped debug data in = " 2702 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2721 decompositionMethod& decomposer,
2722 fvMeshDistribute& distributor,
2724 const scalar maxLoadUnbalance,
2725 const label maxCellUnbalance
2731 refine(cellsToRefine);
2735 Pout<<
"Writing refined but unbalanced " << msg
2743 Pout<<
"Dumped debug data in = " 2744 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2750 Info<<
"Refined mesh in = " 2751 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2752 printMeshInfo(
debug,
"After refinement " + msg,
true);
2760 auto distMap = balance
2779 decompositionMethod& decomposer,
2780 fvMeshDistribute& distributor,
2782 const scalar maxLoadUnbalance,
2783 const label maxCellUnbalance
2786 labelList cellsToRefine(initCellsToRefine);
2816 auto distMap = balance
2831 refine(cellsToRefine);
2835 Pout<<
"Writing refined " << msg
2843 Pout<<
"Dumped debug data in = " 2844 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2850 Info<<
"Refined mesh in = " 2851 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
2865 printMeshInfo(
debug,
"After refinement " + msg,
true);
2873 const label maxGlobalCells,
2874 const label maxLocalCells,
2879 const labelList& cellLevel = meshCutter_.cellLevel();
2880 const pointField& cellCentres = mesh_.cellCentres();
2882 label totNCells = mesh_.globalData().nTotalCells();
2886 if (totNCells >= maxGlobalCells)
2888 Info<<
"No cells marked for refinement since reached limit " 2889 << maxGlobalCells <<
'.' <<
endl;
2903 shells_.findDirectionalLevel
2913 forAll(insideShell, celli)
2915 if (insideShell[celli] >= 0)
2917 bool reachedLimit = !markForRefine
2929 Pout<<
"Stopped refining cells" 2930 <<
" since reaching my cell limit of " 2931 << mesh_.nCells()+nAllowRefine <<
endl;
2942 label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
2945 Info<<
"Unmarked for refinement due to limit shells" 2946 <<
" : " << nUnmarked <<
" cells." <<
endl;
2955 cellsToRefine.setSize(nRefine);
2958 forAll(refineCell, cellI)
2960 if (refineCell[cellI] != -1)
2962 cellsToRefine[nRefine++] = cellI;
2967 return cellsToRefine;
2981 List<refineCell> refCells(cellsToRefine.size());
2984 refCells[i] = refineCell(cellsToRefine[i], refDir);
2988 hexCellLooper cellWalker(mesh_);
2991 cellCuts cuts(mesh_, cellWalker, refCells);
2999 meshRefiner.setRefinement(cuts, meshMod);
3003 mesh_.moving(
false);
3008 mesh_.updateMesh(*morphMap);
3011 if (morphMap().hasMotionPoints())
3013 mesh_.movePoints(morphMap().preMotionPoints());
3025 meshRefiner.updateMesh(*morphMap);
3028 updateMesh(*morphMap, getChangedFaces(*morphMap, cellsToRefine));
void size(const label n)
Older name for setAddressableSize.
const polyMesh & mesh() const noexcept
Return polyMesh.
autoPtr< mapDistributePolyMesh > refineAndBalance(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance, const label maxCellUnbalance)
Refine some cells and rebalance.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
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.
constexpr char nl
The newline '\n' character (0x0a)
UIndirectList< label > labelUIndList
UIndirectList of labels.
Ostream & endl(Ostream &os)
Add newline and flush stream.
A traits class, which is primarily used for primitives and vector-space.
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
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
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< vector > vectorList
List of vector.
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.
List< labelList > labelListList
List of labelList.
Description of cuts across cells.
static bool less(const vector &x, const vector &y)
To compare normals.
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). It 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)
const labelList & cellMap() const noexcept
Old cell map.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
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.
virtual const labelList & faceOwner() const
Return face owner.
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
const globalMeshData & globalData() const
Return parallel info (demand-driven)
label nInternalFaces() const noexcept
Number of internal faces.
autoPtr< mapDistributePolyMesh > balance(const bool keepZoneFaces, const bool keepBaffles, const scalarField &cellWeights, decompositionMethod &decomposer, fvMeshDistribute &distributor)
Redecompose according to cell count.
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 for True value at specified position, never auto-vivify entries.
autoPtr< mapDistributePolyMesh > balanceAndRefine(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance, const label maxCellUnbalance)
Balance before refining some cells.
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.
label nOldCells() const noexcept
Number of old cells.
bool hasMotionPoints() const noexcept
Has valid preMotionPoints?
autoPtr< mapPolyMesh > refine(const labelList &cellsToRefine)
Refine some cells.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
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.
normalLess(const vectorList &values)
messageStream Info
Information stream (stdout output on master, null elsewhere)
label nTotalFaces() const noexcept
Total global number of mesh faces. Not compensated for duplicate faces!
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
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.
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< point > pointList
List of point.
List< label > labelList
A List of labels.
A patch is a list of labels that address the faces in the global face list.
label checkData(const objectRegistry &obr, const instantList &timeDirs, wordList &objectNames, const fileName &local=fileName::null)
Check if fields are good to use (available at all times)
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.
const pointField & preMotionPoints() const noexcept
Pre-motion point positions.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
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())
bool operator()(const label a, const label b) const
static constexpr const zero Zero
Global zero (0)