64 static const scalar edgeTol = 1e-3;
67 void writeSet(
const cellSet& cells,
const string& msg)
69 Info<<
"Writing " << msg <<
" (" << cells.
size() <<
") to cellSet " 84 if (
mag(normal &
vector(1, 0, 0)) > 1-edgeTol)
88 else if (
mag(normal &
vector(0, 1, 0)) > 1-edgeTol)
92 else if (
mag(normal &
vector(0, 0, 1)) > 1-edgeTol)
111 scalar maxX = -GREAT;
115 scalar maxY = -GREAT;
119 scalar maxZ = -GREAT;
122 scalar minOther = GREAT;
123 scalar maxOther = -GREAT;
129 const edge& e = edges[edgeI];
133 scalar eMag =
mag(eVec);
137 if (
mag(eVec & x) > 1-edgeTol)
139 minX =
min(minX, eMag);
140 maxX =
max(maxX, eMag);
143 else if (
mag(eVec & y) > 1-edgeTol)
145 minY =
min(minY, eMag);
146 maxY =
max(maxY, eMag);
149 else if (
mag(eVec & z) > 1-edgeTol)
151 minZ =
min(minZ, eMag);
152 maxZ =
max(maxZ, eMag);
157 minOther =
min(minOther, eMag);
158 maxOther =
max(maxOther, eMag);
163 <<
"Mesh edge statistics:" <<
nl 164 <<
" x aligned : number:" << nX <<
"\tminLen:" << minX
165 <<
"\tmaxLen:" << maxX <<
nl 166 <<
" y aligned : number:" << nY <<
"\tminLen:" << minY
167 <<
"\tmaxLen:" << maxY <<
nl 168 <<
" z aligned : number:" << nZ <<
"\tminLen:" << minZ
169 <<
"\tmaxLen:" << maxZ <<
nl 170 <<
" other : number:" << mesh.
nEdges() - nX - nY - nZ
171 <<
"\tminLen:" << minOther
172 <<
"\tmaxLen:" << maxOther <<
nl <<
endl;
174 if (excludeCmpt == 0)
176 return min(minY,
min(minZ, minOther));
178 else if (excludeCmpt == 1)
180 return min(minX,
min(minZ, minOther));
182 else if (excludeCmpt == 2)
184 return min(minX,
min(minY, minOther));
188 return min(minX,
min(minY,
min(minZ, minOther)));
215 emptyPolyPatch::typeName
235 Info<<
"Created patch oldInternalFaces at " << patchi <<
endl;
239 Info<<
"Reusing patch oldInternalFaces at " << patchi <<
endl;
248 void selectCurvatureCells
253 const scalar nearDist,
254 const scalar curvature,
286 void addCutNeighbours
289 const bool selectInside,
290 const bool selectOutside,
302 const label celli = iter.key();
307 const label facei = cFaces[i];
318 if (selectInside && inside.
found(nbr))
320 addCutFaces.insert(nbr);
322 else if (selectOutside && outside.
found(nbr))
324 addCutFaces.insert(nbr);
330 Info<<
" Selected an additional " << addCutFaces.size()
331 <<
" neighbours of cutCells to refine" <<
endl;
335 cutCells.
insert(iter.key());
344 bool limitRefinementLevel
347 const label limitDiff,
356 if (!excludeCells.
found(celli))
362 label nbr = cCells[i];
364 if (!excludeCells.
found(nbr))
366 if (refLevel[celli] - refLevel[nbr] >= limitDiff)
369 <<
"Level difference between neighbouring cells " 370 << celli <<
" and " << nbr
371 <<
" greater than or equal to " << limitDiff << endl
372 <<
"refLevels:" << refLevel[celli] <<
' ' 386 const label celli = iter.key();
391 const label nbr = cCells[i];
393 if (!excludeCells.
found(nbr) && !cutCells.
found(nbr))
395 if (refLevel[celli] + 1 - refLevel[nbr] >= limitDiff)
397 addCutCells.insert(nbr);
403 if (addCutCells.size())
407 Info<<
"Added an additional " << addCutCells.size() <<
" cells" 408 <<
" to satisfy 1:" << limitDiff <<
" refinement level" 413 cutCells.
insert(iter.key());
419 Info<<
"Added no additional cells" 420 <<
" to satisfy 1:" << limitDiff <<
" refinement level" 454 for (
label celli = oldCells; celli < mesh.
nCells(); celli++)
461 forAll(addedCells, oldCelli)
463 const labelList& added = addedCells[oldCelli];
469 label masterLevel = ++refLevel[oldCelli];
473 refLevel[added[i]] = masterLevel;
484 const label writeMesh,
495 Info<<
"Mesh has:" << mesh.
nCells() <<
" cells." 496 <<
" Removing:" << cellLabels.size() <<
" cells" <<
endl;
499 labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
502 cellRemover.setRefinement
517 if (morphMap().hasMotionPoints())
519 mesh.
movePoints(morphMap().preMotionPoints());
523 cellRemover.updateMesh(morphMap());
526 const labelList& cellMap = morphMap().cellMap();
532 newRefLevel[i] = refLevel[cellMap[i]];
540 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl 566 const bool selectCut,
567 const bool selectInside,
568 const bool selectOutside,
570 const label nCutLayers,
597 selected.
addSet(cutCells);
608 Info<<
"Determined cell status:" << endl
609 <<
" inside :" << inside.size() << endl
610 <<
" outside :" << outside.size() << endl
611 <<
" cutCells:" << cutCells.
size() << endl
612 <<
" selected:" << selected.
size() << endl
615 writeSet(inside,
"inside cells");
616 writeSet(outside,
"outside cells");
617 writeSet(cutCells,
"cut cells");
618 writeSet(selected,
"selected cells");
623 int main(
int argc,
char *argv[])
632 label newPatchi = addPatch(mesh,
"oldInternalFaces");
639 Info<<
"Checking for motionProperties\n" <<
endl;
653 if (motionObj.headerOk())
655 Info<<
"Reading " << runTime.
constant() /
"motionProperties" 660 Switch twoDMotion(motionProperties.lookup(
"twoDMotion"));
664 Info<<
"Correcting for 2D motion" << endl <<
endl;
673 "autoRefineMeshDict",
698 Info<<
"Cells to be used for meshing (0=false, 1=true):" <<
nl 699 <<
" cells cut by surface : " << selectCut <<
nl 700 <<
" cells inside of surface : " << selectInside <<
nl 701 <<
" cells outside of surface : " << selectOutside <<
nl 702 <<
" hanging cells : " << selectHanging <<
nl 706 if (nCutLayers > 0 && selectInside)
709 <<
"Illogical settings : Both nCutLayers>0 and selectInside true." 711 <<
"This would mean that inside cells get removed but should be" 712 <<
" included in final mesh" <<
endl;
725 forAll(outsidePts, outsideI)
727 const point& outsidePoint = outsidePts[outsideI];
729 if (queryMesh.findCell(outsidePoint, -1,
false) == -1)
732 <<
"outsidePoint " << outsidePoint
733 <<
" is not inside any cell" 759 Info<<
"Read existing refinement level from file " 761 <<
" min level : " <<
min(refLevel) <<
nl 762 <<
" max level : " << maxLevel <<
nl 767 Info<<
"Created zero refinement level in file " 776 direction normalDir(getNormalDir(correct2DPtr));
777 scalar meshMinEdgeLen = getEdgeStats(mesh, normalDir);
779 while (meshMinEdgeLen > minEdgeLen)
782 cellSet inside(mesh,
"inside", mesh.nCells()/10);
783 cellSet outside(mesh,
"outside", mesh.nCells()/10);
784 cellSet cutCells(mesh,
"cutCells", mesh.nCells()/10);
785 cellSet selected(mesh,
"selected", mesh.nCells()/10);
806 Info<<
" Selected " << cutCells.
name() <<
" with " 807 << cutCells.
size() <<
" cells" <<
endl;
809 if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
814 cellSet curveCells(mesh,
"curveCells", mesh.nCells()/10);
826 Info<<
" Selected " << curveCells.name() <<
" with " 827 << curveCells.size() <<
" cells" <<
endl;
846 cutCells.
subset(curveCells);
848 Info<<
" Removed cells not on surface curvature. Selected " 857 subsetMesh(mesh, writeMesh, newPatchi, inside, cutCells, refLevel);
876 Info<<
" Current cells : " << mesh.nCells() <<
nl 877 <<
" Selected for refinement :" << cutCells.
size() <<
nl 880 if (cutCells.
empty())
882 Info<<
"Stopping refining since 0 cells selected to be refined ..." 887 if ((mesh.nCells() + 8*cutCells.
size()) > cellLimit)
889 Info<<
"Stopping refining since cell limit reached ..." <<
nl 890 <<
"Would refine from " << mesh.nCells()
891 <<
" to " << mesh.nCells() + 8*cutCells.
size() <<
" cells." 904 Info<<
" After refinement:" << mesh.nCells() <<
nl 909 Info<<
" Writing refined mesh to time " << runTime.
timeName()
918 meshMinEdgeLen = getEdgeStats(mesh, normalDir);
925 cellSet inside(mesh,
"inside", mesh.nCells()/10);
926 cellSet outside(mesh,
"outside", mesh.nCells()/10);
927 cellSet cutCells(mesh,
"cutCells", mesh.nCells()/10);
928 cellSet selected(mesh,
"selected", mesh.nCells()/10);
954 Info<<
"Detected " << hanging.
size() <<
" hanging cells" 955 <<
" (cells with all points on" 956 <<
" outside of cellSet 'selected').\nThese would probably result" 957 <<
" in flattened cells when snapping the mesh to the surface" 960 Info<<
"Refining " << hanging.
size() <<
" hanging cells" <<
nl 977 doRefinement(mesh, refineDict, hanging, refLevel);
979 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl 990 Info<<
"Writing refined mesh to time " << runTime.
timeName() <<
nl const Time & time() const
Return time.
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
const labelListList & cellCells() const
fileName objectPath() const
Return complete path + object name.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
virtual autoPtr< polyPatch > clone(const polyBoundaryMesh &bm) const
Construct and return a clone, resetting the boundary mesh.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const triSurface & surface() const
Return reference to the surface.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Cell-face mesh analysis engine.
const vector & planeNormal() const
Return plane normal.
static unsigned int defaultPrecision()
Return the default precision.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Does multiple pass refinement to refine cells in multiple directions.
void size(const label)
Override size to be inconsistent with allocated storage.
Given list of cells to remove insert all the topology changes.
static word defaultRegion
Return the default region name.
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual void updateMesh(const mapPolyMesh &morphMap)
Update any stored data for new labels.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
virtual const pointField & points() const =0
Return mesh points.
static void noParallel()
Remove the parallel options.
A bounding box defined in terms of the points at its extremities.
label size() const
Return number of elements in table.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Vector< scalar > vector
A scalar version of the templated Vector.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
bool insert(const Key &key)
Insert a new entry.
const cellList & cells() const
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
static labelHashSet getHangingCells(const primitiveMesh &mesh, const labelHashSet &internalCells)
Get cells using points on 'outside' only.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Helper class to search on triSurface.
vectorField pointField
pointField is a vectorField.
virtual void addSet(const topoSet &set)
Add elements present in set.
label start() const
Return start label of this patch in the polyMesh face list.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Class applies a two-dimensional correction to mesh motion point field.
A class for handling words, derived from string.
const word & constant() const
Return constant name.
const vectorField & cellCentres() const
virtual void subset(const topoSet &set)
Subset contents. Only elements present in both sets remain.
List< label > labelList
A List of labels.
const fileName & local() const
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
bool empty() const
Return true if the hash table is empty.
label readLabel(Istream &is)
List< Key > toc() const
Return the table of contents.
static void getSurfaceSets(const polyMesh &mesh, const fileName &surfName, const triSurface &surf, const triSurfaceSearch &querySurf, const pointField &outsidePts, const label nCutLayers, labelHashSet &inside, labelHashSet &outside, labelHashSet &cut)
Divide cells into cut,inside and outside.
bool found(const Key &) const
Return true if hashedEntry is found in table.
void removeBoundary()
Remove boundary patches.
Empty front and back plane patch. Used for 2-D geometries.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void addPatches(const List< polyPatch * > &, const bool validBoundary=true)
Add boundary patches.
A topoSetSource to select cells based on relation to surface.
void setSize(const label)
Reset size of List.
virtual const labelList & faceNeighbour() const
Return face neighbour.
#define WarningInFunction
Report a warning using Foam::Warning.
string & expand(const bool allowEmpty=false)
Expand initial tildes and all occurences of environment variables.
A collection of cell labels.
Direct mesh changes based on v1.3 polyTopoChange syntax.
vector vec(const pointField &) const
Return the vector (end - start)
virtual bool write() const
Write using setting from DB.
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Mesh consisting of general polyhedral cells.
const word & system() const
Return system name.
virtual const labelList & faceOwner() const
Return face owner.
A patch is a list of labels that address the faces in the global face list.
Triangulated surface description with patch information.
label findPatchID(const word &patchName) const
Find patch index given a name.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const fileName & instance() const
label nInternalFaces() const
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
const word & name() const
Return name.
label size() const
Return the number of elements in the UPtrList.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.