70static const scalar edgeTol = 1
e-3;
75 Info<<
"Writing " << msg <<
" (" <<
cells.size() <<
") to cellSet "
90 if (
mag(normal &
vector(1, 0, 0)) > 1-edgeTol)
94 else if (
mag(normal &
vector(0, 1, 0)) > 1-edgeTol)
98 else if (
mag(normal &
vector(0, 0, 1)) > 1-edgeTol)
117 scalar maxX = -GREAT;
121 scalar maxY = -GREAT;
125 scalar maxZ = -GREAT;
128 scalar minOther = GREAT;
129 scalar maxOther = -GREAT;
135 const edge&
e = edges[edgei];
139 scalar eMag =
mag(eVec);
143 if (
mag(eVec &
x) > 1-edgeTol)
145 minX =
min(minX, eMag);
146 maxX =
max(maxX, eMag);
149 else if (
mag(eVec &
y) > 1-edgeTol)
151 minY =
min(minY, eMag);
152 maxY =
max(maxY, eMag);
155 else if (
mag(eVec & z) > 1-edgeTol)
157 minZ =
min(minZ, eMag);
158 maxZ =
max(maxZ, eMag);
163 minOther =
min(minOther, eMag);
164 maxOther =
max(maxOther, eMag);
169 <<
"Mesh edge statistics:" <<
nl
170 <<
" x aligned : number:" << nX <<
"\tminLen:" << minX
171 <<
"\tmaxLen:" << maxX <<
nl
172 <<
" y aligned : number:" << nY <<
"\tminLen:" << minY
173 <<
"\tmaxLen:" << maxY <<
nl
174 <<
" z aligned : number:" << nZ <<
"\tminLen:" << minZ
175 <<
"\tmaxLen:" << maxZ <<
nl
176 <<
" other : number:" <<
mesh.nEdges() - nX - nY - nZ
177 <<
"\tminLen:" << minOther
178 <<
"\tmaxLen:" << maxOther <<
nl <<
endl;
180 if (excludeCmpt == 0)
182 return min(minY,
min(minZ, minOther));
184 else if (excludeCmpt == 1)
186 return min(minX,
min(minZ, minOther));
188 else if (excludeCmpt == 2)
190 return min(minX,
min(minY, minOther));
194 return min(minX,
min(minY,
min(minZ, minOther)));
202 label patchi =
mesh.boundaryMesh().findPatchID(patchName);
218 mesh.nInternalFaces(),
221 emptyPolyPatch::typeName
238 mesh.removeBoundary();
239 mesh.addPatches(newPatches);
241 Info<<
"Created patch oldInternalFaces at " << patchi <<
endl;
245 Info<<
"Reusing patch oldInternalFaces at " << patchi <<
endl;
254void selectCurvatureCells
259 const scalar nearDist,
260 const scalar curvature,
286 cutSource.applyToSet(topoSetSource::ADD,
cells);
295 const bool selectInside,
296 const bool selectOutside,
306 for (
const label celli : cutCells)
312 const label facei = cFaces[i];
314 if (
mesh.isInternalFace(facei))
316 label nbr =
mesh.faceOwner()[facei];
320 nbr =
mesh.faceNeighbour()[facei];
323 if (selectInside && inside.
found(nbr))
325 addCutFaces.insert(nbr);
327 else if (selectOutside && outside.
found(nbr))
329 addCutFaces.insert(nbr);
335 Info<<
" Selected an additional " << addCutFaces.size()
336 <<
" neighbours of cutCells to refine" <<
endl;
338 for (
const label facei : addCutFaces)
340 cutCells.insert(facei);
349bool limitRefinementLevel
352 const label limitDiff,
361 if (!excludeCells.
found(celli))
367 label nbr = cCells[i];
369 if (!excludeCells.
found(nbr))
371 if (refLevel[celli] - refLevel[nbr] >= limitDiff)
374 <<
"Level difference between neighbouring cells "
375 << celli <<
" and " << nbr
376 <<
" greater than or equal to " << limitDiff <<
endl
377 <<
"refLevels:" << refLevel[celli] <<
' '
378 << refLevel[nbr] <<
abort(FatalError);
388 for (
const label celli : cutCells)
395 const label nbr = cCells[i];
397 if (!excludeCells.
found(nbr) && !cutCells.found(nbr))
399 if (refLevel[celli] + 1 - refLevel[nbr] >= limitDiff)
401 addCutCells.insert(nbr);
407 if (addCutCells.size())
411 Info<<
"Added an additional " << addCutCells.size() <<
" cells"
412 <<
" to satisfy 1:" << limitDiff <<
" refinement level"
415 for (
const label celli : addCutCells)
417 cutCells.insert(celli);
422 Info<<
"Added no additional cells"
423 <<
" to satisfy 1:" << limitDiff <<
" refinement level"
440 label oldCells =
mesh.nCells();
456 for (label celli = oldCells; celli <
mesh.nCells(); celli++)
463 forAll(addedCells, oldCelli)
465 const labelList& added = addedCells[oldCelli];
471 label masterLevel = ++refLevel[oldCelli];
475 refLevel[added[i]] = masterLevel;
486 const label writeMesh,
497 Info<<
"Mesh has:" <<
mesh.nCells() <<
" cells."
498 <<
" Removing:" << cellLabels.size() <<
" cells" <<
endl;
501 labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
504 cellRemover.setRefinement
519 if (morphMap().hasMotionPoints())
521 mesh.movePoints(morphMap().preMotionPoints());
525 cellRemover.updateMesh(morphMap());
528 const labelList& cellMap = morphMap().cellMap();
534 newRefLevel[i] = refLevel[cellMap[i]];
542 Info<<
"Writing refined mesh to time " <<
runTime.timeName() <<
nl
545 IOstream::defaultPrecision(
max(10u, IOstream::defaultPrecision()));
568 const bool selectCut,
569 const bool selectInside,
570 const bool selectOutside,
572 const label nCutLayers,
581 surfaceSets::getSurfaceSets
599 selected.
addSet(cutCells);
610 Info<<
"Determined cell status:" <<
endl
611 <<
" inside :" << inside.
size() <<
endl
612 <<
" outside :" << outside.
size() <<
endl
613 <<
" cutCells:" << cutCells.
size() <<
endl
614 <<
" selected:" << selected.
size() <<
endl
617 writeSet(inside,
"inside cells");
618 writeSet(outside,
"outside cells");
619 writeSet(cutCells,
"cut cells");
620 writeSet(selected,
"selected cells");
625int main(
int argc,
char *argv[])
629 "Refine cells near to a surface"
631 argList::noParallel();
638 label newPatchi = addPatch(
mesh,
"oldInternalFaces");
645 Info<<
"Checking for motionProperties\n" <<
endl;
652 IOobject::MUST_READ_IF_MODIFIED,
661 Info<<
"Reading " <<
runTime.constant() /
"motionProperties"
666 if (motionProperties.get<
bool>(
"twoDMotion"))
677 "snappyRefineMeshDict",
680 IOobject::MUST_READ_IF_MODIFIED,
687 const label nCutLayers(refineDict.
get<label>(
"nCutLayers"));
688 const label cellLimit(refineDict.
get<label>(
"cellLimit"));
689 const bool selectCut(refineDict.
get<
bool>(
"selectCut"));
690 const bool selectInside(refineDict.
get<
bool>(
"selectInside"));
691 const bool selectOutside(refineDict.
get<
bool>(
"selectOutside"));
692 const bool selectHanging(refineDict.
get<
bool>(
"selectHanging"));
694 const scalar minEdgeLen(refineDict.
get<scalar>(
"minEdgeLen"));
695 const scalar maxEdgeLen(refineDict.
get<scalar>(
"maxEdgeLen"));
696 const scalar curvature(refineDict.
get<scalar>(
"curvature"));
697 const scalar curvDist(refineDict.
get<scalar>(
"curvatureDistance"));
699 const label refinementLimit(refineDict.
get<label>(
"splitLevel"));
700 const bool writeMesh(refineDict.
get<
bool>(
"writeMesh"));
702 Info<<
"Cells to be used for meshing (0=false, 1=true):" <<
nl
703 <<
" cells cut by surface : " << selectCut <<
nl
704 <<
" cells inside of surface : " << selectInside <<
nl
705 <<
" cells outside of surface : " << selectOutside <<
nl
706 <<
" hanging cells : " << selectHanging <<
nl
710 if (nCutLayers > 0 && selectInside)
713 <<
"Illogical settings : Both nCutLayers>0 and selectInside true."
715 <<
"This would mean that inside cells get removed but should be"
716 <<
" included in final mesh" <<
endl;
729 forAll(outsidePts, outsidei)
731 const point& outsidePoint = outsidePts[outsidei];
733 if (queryMesh.findCell(outsidePoint, -1,
false) == -1)
736 <<
"outsidePoint " << outsidePoint
737 <<
" is not inside any cell"
751 polyMesh::defaultRegion,
753 IOobject::READ_IF_PRESENT,
759 label maxLevel =
max(refLevel);
763 Info<<
"Read existing refinement level from file "
765 <<
" min level : " <<
min(refLevel) <<
nl
766 <<
" max level : " << maxLevel <<
nl
771 Info<<
"Created zero refinement level in file "
780 direction normalDir(getNormalDir(correct2DPtr));
781 scalar meshMinEdgeLen = getEdgeStats(
mesh, normalDir);
783 while (meshMinEdgeLen > minEdgeLen)
810 Info<<
" Selected " << cutCells.
name() <<
" with "
811 << cutCells.
size() <<
" cells" <<
endl;
813 if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
830 Info<<
" Selected " << curveCells.name() <<
" with "
831 << curveCells.size() <<
" cells" <<
endl;
850 cutCells.
subset(curveCells);
852 Info<<
" Removed cells not on surface curvature. Selected "
861 subsetMesh(
mesh, writeMesh, newPatchi, inside, cutCells, refLevel);
880 Info<<
" Current cells : " <<
mesh.nCells() <<
nl
881 <<
" Selected for refinement :" << cutCells.
size() <<
nl
884 if (cutCells.
empty())
886 Info<<
"Stopping refining since 0 cells selected to be refined ..."
891 if ((
mesh.nCells() + 8*cutCells.
size()) > cellLimit)
893 Info<<
"Stopping refining since cell limit reached ..." <<
nl
894 <<
"Would refine from " <<
mesh.nCells()
895 <<
" to " <<
mesh.nCells() + 8*cutCells.
size() <<
" cells."
908 Info<<
" After refinement:" <<
mesh.nCells() <<
nl
913 Info<<
" Writing refined mesh to time " <<
runTime.timeName()
916 IOstream::defaultPrecision(
max(10u, IOstream::defaultPrecision()));
922 meshMinEdgeLen = getEdgeStats(
mesh, normalDir);
958 Info<<
"Detected " << hanging.
size() <<
" hanging cells"
959 <<
" (cells with all points on"
960 <<
" outside of cellSet 'selected').\nThese would probably result"
961 <<
" in flattened cells when snapping the mesh to the surface"
964 Info<<
"Refining " << hanging.
size() <<
" hanging cells" <<
nl
981 doRefinement(
mesh, refineDict, hanging, refLevel);
983 Info<<
"Writing refined mesh to time " <<
runTime.timeName() <<
nl
987 IOstream::defaultPrecision(
max(10u, IOstream::defaultPrecision()));
994 Info<<
"Writing refined mesh to time " <<
runTime.timeName() <<
nl
998 IOstream::defaultPrecision(
max(10u, IOstream::defaultPrecision()));
List< Key > toc() const
The table of contents (the keys) in unsorted order.
bool empty() const noexcept
True if the hash table is empty.
bool found(const Key &key) const
Return true if hashed entry is found in table.
label size() const noexcept
The number of elements in table.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const word & name() const noexcept
Return the object name.
fileName objectPath() const
The complete path + object name.
void transfer(List< T > &list)
void setSize(const label n)
Alias for resize()
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
void size(const label n)
Older name for setAddressableSize.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
A bounding box defined in terms of min/max extrema points.
A collection of cell labels.
virtual void updateMesh(const mapPolyMesh &morphMap)
Update any stored data for new labels.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Empty front and back plane patch. Used for 2-D geometries.
A class for handling file names.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
Does multiple pass refinement to refine cells in multiple directions.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
virtual autoPtr< polyPatch > clone(const labelList &faceCells) const
Construct and return a clone, setting faceCells.
label start() const
Return start label of this patch in the polyMesh face list.
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
virtual bool write(const bool valid=true) const
Write using setting from DB.
Given list of cells to remove, insert all the topology changes.
string & expand(const bool allowEmpty=false)
A topoSetCellSource to select cells based on relation to a surface given by an external file.
virtual void addSet(const topoSet &set)
Add elements present in set.
virtual void subset(const topoSet &set)
Subset contents. Only elements present in both sets remain.
Helper class to search on triSurface.
const triSurface & surface() const
Return reference to the surface.
Triangulated surface description with patch information.
Class applies a two-dimensional correction to mesh motion point field.
const vector & planeNormal() const
Return plane normal.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.