80 static const scalar edgeTol = 1
E-3;
85 Info<<
"Writing " << msg <<
" (" << cells.
size() <<
") to cellSet "
100 if (
mag(normal &
vector(1, 0, 0)) > 1-edgeTol)
104 else if (
mag(normal &
vector(0, 1, 0)) > 1-edgeTol)
108 else if (
mag(normal &
vector(0, 0, 1)) > 1-edgeTol)
127 scalar maxX = -GREAT;
131 scalar maxY = -GREAT;
135 scalar maxZ = -GREAT;
138 scalar minOther = GREAT;
139 scalar maxOther = -GREAT;
145 const edge&
e = edges[edgeI];
149 scalar eMag =
mag(eVec);
153 if (
mag(eVec & x) > 1-edgeTol)
155 minX =
min(minX, eMag);
156 maxX =
max(maxX, eMag);
159 else if (
mag(eVec &
y) > 1-edgeTol)
161 minY =
min(minY, eMag);
162 maxY =
max(maxY, eMag);
165 else if (
mag(eVec & z) > 1-edgeTol)
167 minZ =
min(minZ, eMag);
168 maxZ =
max(maxZ, eMag);
173 minOther =
min(minOther, eMag);
174 maxOther =
max(maxOther, eMag);
179 <<
"Mesh edge statistics:" <<
nl
180 <<
" x aligned : number:" << nX <<
"\tminLen:" << minX
181 <<
"\tmaxLen:" << maxX <<
nl
182 <<
" y aligned : number:" << nY <<
"\tminLen:" << minY
183 <<
"\tmaxLen:" << maxY <<
nl
184 <<
" z aligned : number:" << nZ <<
"\tminLen:" << minZ
185 <<
"\tmaxLen:" << maxZ <<
nl
186 <<
" other : number:" << mesh.
nEdges() - nX - nY - nZ
187 <<
"\tminLen:" << minOther
188 <<
"\tmaxLen:" << maxOther <<
nl <<
endl;
190 if (excludeCmpt == 0)
192 return min(minY,
min(minZ, minOther));
194 else if (excludeCmpt == 1)
196 return min(minX,
min(minZ, minOther));
198 else if (excludeCmpt == 2)
200 return min(minX,
min(minY, minOther));
204 return min(minX,
min(minY,
min(minZ, minOther)));
250 Info<<
"Created patch oldInternalFaces at " << patchI <<
endl;
254 Info<<
"Reusing patch oldInternalFaces at " << patchI <<
endl;
263 void selectCurvatureCells
268 const scalar nearDist,
269 const scalar curvature,
300 void addCutNeighbours
303 const bool selectInside,
304 const bool selectOutside,
317 iter != cutCells.
end();
321 label cellI = iter.key();
327 label faceI = cFaces[i];
338 if (selectInside && inside.
found(nbr))
340 addCutFaces.insert(nbr);
342 else if (selectOutside && outside.
found(nbr))
344 addCutFaces.insert(nbr);
350 Info<<
" Selected an additional " << addCutFaces.size()
351 <<
" neighbours of cutCells to refine" <<
endl;
356 iter != addCutFaces.end();
360 cutCells.
insert(iter.key());
369 bool limitRefinementLevel
372 const label limitDiff,
381 if (!excludeCells.
found(cellI))
387 label nbr = cCells[i];
389 if (!excludeCells.
found(nbr))
391 if (refLevel[cellI] - refLevel[nbr] >= limitDiff)
394 <<
"Level difference between neighbouring cells "
395 << cellI <<
" and " << nbr
396 <<
" greater than or equal to " << limitDiff << endl
397 <<
"refLevels:" << refLevel[cellI] <<
' '
411 iter != cutCells.
end();
416 label cellI = iter.key();
422 label nbr = cCells[i];
424 if (!excludeCells.
found(nbr) && !cutCells.
found(nbr))
426 if (refLevel[cellI] + 1 - refLevel[nbr] >= limitDiff)
428 addCutCells.insert(nbr);
434 if (addCutCells.size())
438 Info<<
"Added an additional " << addCutCells.size() <<
" cells"
439 <<
" to satisfy 1:" << limitDiff <<
" refinement level"
444 iter != addCutCells.end();
448 cutCells.
insert(iter.key());
454 Info<<
"Added no additional cells"
455 <<
" to satisfy 1:" << limitDiff <<
" refinement level"
473 label oldCells = mesh.
nCells();
489 for (label cellI = oldCells; cellI < mesh.
nCells(); cellI++)
496 forAll(addedCells, oldCellI)
498 const labelList& added = addedCells[oldCellI];
504 label masterLevel = ++refLevel[oldCellI];
508 refLevel[added[i]] = masterLevel;
519 const label writeMesh,
530 Info<<
"Mesh has:" << mesh.
nCells() <<
" cells."
531 <<
" Removing:" << cellLabels.size() <<
" cells" <<
endl;
534 labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
537 cellRemover.setRefinement
552 if (morphMap().hasMotionPoints())
554 mesh.
movePoints(morphMap().preMotionPoints());
558 cellRemover.updateMesh(morphMap());
561 const labelList& cellMap = morphMap().cellMap();
567 newRefLevel[i] = refLevel[cellMap[i]];
575 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl
601 const bool selectCut,
602 const bool selectInside,
603 const bool selectOutside,
605 const label nCutLayers,
632 selected.
addSet(cutCells);
643 Info<<
"Determined cell status:" << endl
644 <<
" inside :" << inside.size() << endl
645 <<
" outside :" << outside.size() << endl
646 <<
" cutCells:" << cutCells.
size() << endl
647 <<
" selected:" << selected.
size() << endl
650 writeSet(inside,
"inside cells");
651 writeSet(outside,
"outside cells");
652 writeSet(cutCells,
"cut cells");
653 writeSet(selected,
"selected cells");
659 int main(
int argc,
char *argv[])
668 label newPatchI = addPatch(mesh,
"oldInternalFaces");
675 Info<<
"Checking for motionProperties\n" <<
endl;
689 if (motionObj.headerOk())
691 Info<<
"Reading " << runTime.
constant() /
"motionProperties"
696 Switch twoDMotion(motionProperties.lookup(
"twoDMotion"));
700 Info<<
"Correcting for 2D motion" << endl <<
endl;
709 "autoRefineMeshDict",
719 label nCutLayers(readLabel(refineDict.
lookup(
"nCutLayers")));
720 label cellLimit(readLabel(refineDict.
lookup(
"cellLimit")));
731 label refinementLimit(readLabel(refineDict.
lookup(
"splitLevel")));
734 Info<<
"Cells to be used for meshing (0=false, 1=true):" <<
nl
735 <<
" cells cut by surface : " << selectCut <<
nl
736 <<
" cells inside of surface : " << selectInside <<
nl
737 <<
" cells outside of surface : " << selectOutside <<
nl
738 <<
" hanging cells : " << selectHanging <<
nl
742 if (nCutLayers > 0 && selectInside)
745 <<
"Illogical settings : Both nCutLayers>0 and selectInside true."
747 <<
"This would mean that inside cells get removed but should be"
748 <<
" included in final mesh" << endl;
761 forAll(outsidePts, outsideI)
763 const point& outsidePoint = outsidePts[outsideI];
765 if (queryMesh.findCell(outsidePoint, -1,
false) == -1)
768 <<
"outsidePoint " << outsidePoint
769 <<
" is not inside any cell"
791 label maxLevel =
max(refLevel);
795 Info<<
"Read existing refinement level from file "
797 <<
" min level : " <<
min(refLevel) <<
nl
798 <<
" max level : " << maxLevel <<
nl
803 Info<<
"Created zero refinement level in file "
812 direction normalDir(getNormalDir(correct2DPtr));
813 scalar meshMinEdgeLen = getEdgeStats(mesh, normalDir);
815 while (meshMinEdgeLen > minEdgeLen)
818 cellSet inside(mesh,
"inside", mesh.nCells()/10);
819 cellSet outside(mesh,
"outside", mesh.nCells()/10);
820 cellSet cutCells(mesh,
"cutCells", mesh.nCells()/10);
821 cellSet selected(mesh,
"selected", mesh.nCells()/10);
842 Info<<
" Selected " << cutCells.
name() <<
" with "
843 << cutCells.
size() <<
" cells" <<
endl;
845 if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
850 cellSet curveCells(mesh,
"curveCells", mesh.nCells()/10);
862 Info<<
" Selected " << curveCells.name() <<
" with "
863 << curveCells.size() <<
" cells" <<
endl;
882 cutCells.
subset(curveCells);
884 Info<<
" Removed cells not on surface curvature. Selected "
893 subsetMesh(mesh, writeMesh, newPatchI, inside, cutCells, refLevel);
912 Info<<
" Current cells : " << mesh.nCells() <<
nl
913 <<
" Selected for refinement :" << cutCells.
size() <<
nl
916 if (cutCells.
empty())
918 Info<<
"Stopping refining since 0 cells selected to be refined ..."
923 if ((mesh.nCells() + 8*cutCells.
size()) > cellLimit)
925 Info<<
"Stopping refining since cell limit reached ..." <<
nl
926 <<
"Would refine from " << mesh.nCells()
927 <<
" to " << mesh.nCells() + 8*cutCells.
size() <<
" cells."
940 Info<<
" After refinement:" << mesh.nCells() <<
nl
945 Info<<
" Writing refined mesh to time " << runTime.
timeName()
954 meshMinEdgeLen = getEdgeStats(mesh, normalDir);
961 cellSet inside(mesh,
"inside", mesh.nCells()/10);
962 cellSet outside(mesh,
"outside", mesh.nCells()/10);
963 cellSet cutCells(mesh,
"cutCells", mesh.nCells()/10);
964 cellSet selected(mesh,
"selected", mesh.nCells()/10);
990 Info<<
"Detected " << hanging.
size() <<
" hanging cells"
991 <<
" (cells with all points on"
992 <<
" outside of cellSet 'selected').\nThese would probably result"
993 <<
" in flattened cells when snapping the mesh to the surface"
996 Info<<
"Refining " << hanging.
size() <<
" hanging cells" <<
nl
1002 limitRefinementLevel
1013 doRefinement(mesh, refineDict, hanging, refLevel);
1015 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl
1024 else if (!writeMesh)
1026 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl