49 const mapPolyMesh& map,
53 const polyMesh&
mesh = map.mesh();
65 const labelList& faceOwner = mesh.faceOwner();
66 const labelList& faceNeighbour = mesh.faceNeighbour();
68 const label nInternalFaces = mesh.nInternalFaces();
73 forAll(oldCellsToRefine, i)
75 oldRefineCell.set(oldCellsToRefine[i], 1u);
83 for (label faceI = 0; faceI < nInternalFaces; faceI++)
85 label oldOwn = map.cellMap()[faceOwner[faceI]];
86 label oldNei = map.cellMap()[faceNeighbour[faceI]];
91 && oldRefineCell.get(oldOwn) == 0u
93 && oldRefineCell.get(oldNei) == 0u
100 refinedInternalFace.set(faceI, 1u);
107 boolList refinedBoundaryFace(mesh.nFaces()-nInternalFaces,
false);
109 forAll(mesh.boundaryMesh(), patchI)
111 const polyPatch& pp = mesh.boundaryMesh()[patchI];
113 label faceI = pp.start();
117 label oldOwn = map.cellMap()[faceOwner[faceI]];
119 if (oldOwn >= 0 && oldRefineCell.get(oldOwn) == 0u)
125 refinedBoundaryFace[faceI-nInternalFaces] =
true;
144 boolList changedFace(mesh.nFaces(),
false);
146 forAll(refinedInternalFace, faceI)
148 if (refinedInternalFace.get(faceI) == 1u)
150 const cell& ownFaces = cells[faceOwner[faceI]];
153 changedFace[ownFaces[ownI]] =
true;
155 const cell& neiFaces = cells[faceNeighbour[faceI]];
158 changedFace[neiFaces[neiI]] =
true;
163 forAll(refinedBoundaryFace, i)
165 if (refinedBoundaryFace[i])
167 const cell& ownFaces = cells[faceOwner[i+nInternalFaces]];
170 changedFace[ownFaces[ownI]] =
true;
189 forAll(changedFace, faceI)
191 if (changedFace[faceI])
197 changedFaces.setSize(nChanged);
200 forAll(changedFace, faceI)
202 if (changedFace[faceI])
204 changedFaces[nChanged++] = faceI;
209 label nChangedFaces = changedFaces.size();
210 reduce(nChangedFaces, sumOp<label>());
214 Pout<<
"getChangedFaces : Detected "
215 <<
" local:" << changedFaces.size()
216 <<
" global:" << nChangedFaces
217 <<
" changed faces out of " << mesh.globalData().nTotalFaces()
220 faceSet changedFacesSet(mesh,
"changedFaces", changedFaces);
221 Pout<<
"getChangedFaces : Writing " << changedFaces.size()
222 <<
" changed faces to faceSet " << changedFacesSet.
name()
224 changedFacesSet.
write();
234 bool Foam::meshRefinement::markForRefine
236 const label markValue,
237 const label nAllowRefine,
245 cellValue = markValue;
249 return nRefine <= nAllowRefine;
254 Foam::label Foam::meshRefinement::markFeatureRefinement
256 const point& keepPoint,
257 const PtrList<featureEdgeMesh>& featureMeshes,
259 const label nAllowRefine,
281 Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
284 label cellI = mesh_.findCell(keepPoint);
288 forAll(featureMeshes, featI)
290 const featureEdgeMesh& featureMesh = featureMeshes[featI];
293 forAll(pointEdges, pointI)
295 if (pointEdges[pointI].size() != 2)
299 Pout<<
"Adding particle from point:" << pointI
300 <<
" coord:" << featureMesh.points()[pointI]
301 <<
" pEdges:" << pointEdges[pointI]
313 featureMesh.points()[pointI],
314 featureLevels[featI],
326 labelList maxFeatureLevel(mesh_.nCells(), -1);
329 trackedParticle::trackData td(cloud, maxFeatureLevel);
335 maxFeatureLevel = -1;
338 List<PackedBoolList> featureEdgeVisited(featureMeshes.size());
340 forAll(featureMeshes, featI)
342 featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size());
343 featureEdgeVisited[featI] = 0u;
348 label nParticles = 0;
351 forAllIter(Cloud<trackedParticle>, cloud, iter)
353 trackedParticle& tp = iter();
355 label featI = tp.i();
356 label pointI = tp.j();
358 const featureEdgeMesh& featureMesh = featureMeshes[featI];
359 const labelList& pEdges = featureMesh.pointEdges()[pointI];
364 bool keepParticle =
false;
368 label edgeI = pEdges[i];
370 if (featureEdgeVisited[featI].
set(edgeI, 1u))
375 const edge&
e = featureMesh.edges()[edgeI];
376 label otherPointI = e.otherVertex(pointI);
378 tp.
end() = featureMesh.points()[otherPointI];
379 tp.j() = otherPointI;
389 cloud.deleteParticle(tp);
398 reduce(nParticles, sumOp<label>());
412 const labelList& cellLevel = meshCutter_.cellLevel();
414 label oldNRefine = nRefine;
416 forAll(maxFeatureLevel, cellI)
418 if (maxFeatureLevel[cellI] > cellLevel[cellI])
444 Info<<
"Reached refinement limit." <<
endl;
447 return returnReduce(nRefine-oldNRefine, sumOp<label>());
452 Foam::label Foam::meshRefinement::markInternalRefinement
454 const label nAllowRefine,
460 const labelList& cellLevel = meshCutter_.cellLevel();
461 const pointField& cellCentres = mesh_.cellCentres();
463 label oldNRefine = nRefine;
467 labelList testLevels(cellLevel.size()-nRefine);
472 if (refineCell[cellI] == -1)
474 testCc[testI] = cellCentres[cellI];
475 testLevels[testI] = cellLevel[cellI];
482 shells_.findHigherLevel(testCc, testLevels, maxLevel);
489 if (refineCell[cellI] == -1)
491 if (maxLevel[testI] > testLevels[testI])
493 bool reachedLimit = !markForRefine
505 Pout<<
"Stopped refining internal cells"
506 <<
" since reaching my cell limit of "
507 << mesh_.nCells()+7*nRefine <<
endl;
522 Info<<
"Reached refinement limit." <<
endl;
525 return returnReduce(nRefine-oldNRefine, sumOp<label>());
540 forAll(surfaceIndex_, faceI)
542 if (surfaceIndex_[faceI] != -1)
544 label own = mesh_.faceOwner()[faceI];
546 if (mesh_.isInternalFace(faceI))
548 label nei = mesh_.faceNeighbour()[faceI];
550 if (refineCell[own] == -1 || refineCell[nei] == -1)
552 testFaces[nTest++] = faceI;
557 if (refineCell[own] == -1)
559 testFaces[nTest++] = faceI;
564 testFaces.setSize(nTest);
571 Foam::label Foam::meshRefinement::markSurfaceRefinement
573 const label nAllowRefine,
581 const labelList& cellLevel = meshCutter_.cellLevel();
582 const pointField& cellCentres = mesh_.cellCentres();
584 label oldNRefine = nRefine;
592 labelList testFaces(getRefineCandidateFaces(refineCell));
603 label faceI = testFaces[i];
605 label own = mesh_.faceOwner()[faceI];
607 if (mesh_.isInternalFace(faceI))
609 label nei = mesh_.faceNeighbour()[faceI];
611 start[i] = cellCentres[own];
612 end[i] = cellCentres[nei];
613 minLevel[i] =
min(cellLevel[own], cellLevel[nei]);
617 label bFaceI = faceI - mesh_.nInternalFaces();
619 start[i] = cellCentres[own];
620 end[i] = neiCc[bFaceI];
621 minLevel[i] =
min(cellLevel[own], neiLevel[bFaceI]);
631 surfaces_.findHigherIntersection
647 label faceI = testFaces[i];
649 label surfI = surfaceHit[i];
659 label own = mesh_.faceOwner()[faceI];
661 if (surfaceMinLevel[i] > cellLevel[own])
679 if (mesh_.isInternalFace(faceI))
681 label nei = mesh_.faceNeighbour()[faceI];
682 if (surfaceMinLevel[i] > cellLevel[nei])
709 Info<<
"Reached refinement limit." <<
endl;
712 return returnReduce(nRefine-oldNRefine, sumOp<label>());
720 bool Foam::meshRefinement::checkCurvature
722 const scalar curvature,
723 const label nAllowRefine,
725 const label surfaceLevel,
726 const vector& surfaceNormal,
737 const labelList& cellLevel = meshCutter_.cellLevel();
740 if (surfaceLevel > cellLevel[cellI])
742 if (cellMaxLevel == -1)
745 cellMaxLevel = surfaceLevel;
746 cellMaxNormal = surfaceNormal;
751 if ((cellMaxNormal & surfaceNormal) < curvature)
764 if (surfaceLevel > cellMaxLevel)
766 cellMaxLevel = surfaceLevel;
767 cellMaxNormal = surfaceNormal;
778 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
780 const scalar curvature,
781 const label nAllowRefine,
789 const labelList& cellLevel = meshCutter_.cellLevel();
790 const pointField& cellCentres = mesh_.cellCentres();
792 label oldNRefine = nRefine;
803 labelList testFaces(getRefineCandidateFaces(refineCell));
812 label faceI = testFaces[i];
814 label own = mesh_.faceOwner()[faceI];
816 if (mesh_.isInternalFace(faceI))
818 label nei = mesh_.faceNeighbour()[faceI];
820 start[i] = cellCentres[own];
821 end[i] = cellCentres[nei];
822 minLevel[i] =
min(cellLevel[own], cellLevel[nei]);
826 label bFaceI = faceI - mesh_.nInternalFaces();
828 start[i] = cellCentres[own];
829 end[i] = neiCc[bFaceI];
830 minLevel[i] =
min(cellLevel[own], neiLevel[bFaceI]);
837 labelList cellMaxLevel(mesh_.nCells(), -1);
842 List<vectorList> surfaceNormal;
846 surfaces_.findAllHigherIntersections
863 label faceI = testFaces[i];
864 label own = mesh_.faceOwner()[faceI];
866 const vectorList& fNormals = surfaceNormal[i];
867 const labelList& fLevels = surfaceLevel[i];
888 if (mesh_.isInternalFace(faceI))
890 label nei = mesh_.faceNeighbour()[faceI];
917 labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
918 vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
920 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
922 label bFaceI = faceI-mesh_.nInternalFaces();
923 label own = mesh_.faceOwner()[faceI];
925 neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
926 neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
934 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
936 label own = mesh_.faceOwner()[faceI];
937 label nei = mesh_.faceNeighbour()[faceI];
939 if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
942 if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
945 if (cellLevel[own] < cellMaxLevel[own])
960 Pout<<
"Stopped refining since reaching my cell"
961 <<
" limit of " << mesh_.nCells()+7*nRefine
968 if (cellLevel[nei] < cellMaxLevel[nei])
983 Pout<<
"Stopped refining since reaching my cell"
984 <<
" limit of " << mesh_.nCells()+7*nRefine
994 for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
996 label own = mesh_.faceOwner()[faceI];
997 label bFaceI = faceI - mesh_.nInternalFaces();
999 if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
1002 if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
1017 Pout<<
"Stopped refining since reaching my cell"
1018 <<
" limit of " << mesh_.nCells()+7*nRefine
1033 Info<<
"Reached refinement limit." <<
endl;
1036 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1049 const point& keepPoint,
1050 const scalar curvature,
1055 const bool featureRefinement,
1056 const bool internalRefinement,
1057 const bool surfaceRefinement,
1058 const bool curvatureRefinement,
1059 const label maxGlobalCells,
1060 const label maxLocalCells
1063 label totNCells = mesh_.globalData().nTotalCells();
1067 if (totNCells >= maxGlobalCells)
1069 Info<<
"No cells marked for refinement since reached limit "
1070 << maxGlobalCells <<
'.' <<
endl;
1099 labelList refineCell(mesh_.nCells(), -1);
1104 labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1105 pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1106 calcNeighbourData(neiLevel, neiCc);
1113 if (featureRefinement)
1115 label nFeatures = markFeatureRefinement
1126 Info<<
"Marked for refinement due to explicit features : "
1127 << nFeatures <<
" cells." <<
endl;
1133 if (internalRefinement)
1135 label nShell = markInternalRefinement
1142 Info<<
"Marked for refinement due to refinement shells : "
1143 << nShell <<
" cells." <<
endl;
1149 if (surfaceRefinement)
1151 label nSurf = markSurfaceRefinement
1160 Info<<
"Marked for refinement due to surface intersection : "
1161 << nSurf <<
" cells." <<
endl;
1170 && (curvature >= -1 && curvature <= 1)
1171 && (surfaces_.minLevel() != surfaces_.maxLevel())
1174 label nCurv = markSurfaceCurvatureRefinement
1184 Info<<
"Marked for refinement due to curvature/regions : "
1185 << nCurv <<
" cells." <<
endl;
1191 cellsToRefine.
setSize(nRefine);
1194 forAll(refineCell, cellI)
1196 if (refineCell[cellI] != -1)
1198 cellsToRefine[nRefine++] = cellI;
1203 return cellsToRefine;
1216 meshCutter_.setRefinement(cellsToRefine, meshMod);
1222 mesh_.updateMesh(map);
1225 if (map().hasMotionPoints())
1227 mesh_.movePoints(map().preMotionPoints());
1237 mesh_.setInstance(oldInstance());
1241 updateMesh(map, getChangedFaces(map, cellsToRefine));
1256 const scalar maxLoadUnbalance
1260 refine(cellsToRefine);
1264 Pout<<
"Writing refined but unbalanced " << msg
1272 Pout<<
"Dumped debug data in = "
1273 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1279 Info<<
"Refined mesh in = "
1280 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1281 printMeshInfo(debug,
"After refinement " + msg);
1291 scalar nIdealCells =
1292 mesh_.globalData().nTotalCells()
1297 mag(1.0-mesh_.nCells()/nIdealCells),
1301 if (unbalance <= maxLoadUnbalance)
1303 Info<<
"Skipping balancing since max unbalance " << unbalance
1304 <<
" is less than allowable " << maxLoadUnbalance
1320 Info<<
"Balanced mesh in = "
1321 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1323 printMeshInfo(debug,
"After balancing " + msg);
1328 Pout<<
"Writing balanced " << msg
1335 Pout<<
"Dumped debug data in = "
1336 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1356 const scalar maxLoadUnbalance
1359 labelList cellsToRefine(initCellsToRefine);
1394 scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.
size());
1395 scalar nIdealNewCells =
1400 mag(1.0-nNewCells/nIdealNewCells),
1404 if (unbalance <= maxLoadUnbalance)
1406 Info<<
"Skipping balancing since max unbalance " << unbalance
1407 <<
" is less than allowable " << maxLoadUnbalance
1415 cellWeights[cellsToRefine[i]] += 7;
1428 distMap().distributeCellIndices(cellsToRefine);
1430 Info<<
"Balanced mesh in = "
1431 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1446 printMeshInfo(debug,
"After balancing " + msg);
1450 Pout<<
"Writing balanced " << msg
1457 Pout<<
"Dumped debug data in = "
1458 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1469 refine(cellsToRefine);
1473 Pout<<
"Writing refined " << msg
1481 Pout<<
"Dumped debug data in = "
1482 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1488 Info<<
"Refined mesh in = "
1489 << mesh_.time().cpuTimeIncrement() <<
" s" <<
endl;
1503 printMeshInfo(debug,
"After refinement " + msg);