55 void Foam::polyMesh::calcDirections()
const 57 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
65 bool hasEmptyPatches =
false;
66 bool hasWedgePatches =
false;
74 if (isA<emptyPolyPatch>(pp))
82 hasEmptyPatches =
true;
86 else if (isA<wedgePolyPatch>(pp))
88 const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>(pp);
92 (void)wpp.faceNormals();
96 hasWedgePatches =
true;
97 wedgeDirVec +=
cmptMag(wpp.centreNormal());
105 reduce(emptyDirVec, sumOp<vector>());
107 emptyDirVec.normalise();
109 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
111 if (emptyDirVec[cmpt] > 1
e-6)
113 solutionD_[cmpt] = -1;
117 solutionD_[cmpt] = 1;
125 geometricD_ = solutionD_;
129 reduce(wedgeDirVec, sumOp<vector>());
131 wedgeDirVec.normalise();
133 for (
direction cmpt=0; cmpt<vector::nComponents; cmpt++)
135 if (wedgeDirVec[cmpt] > 1
e-6)
137 geometricD_[cmpt] = -1;
141 geometricD_[cmpt] = 1;
171 Foam::polyMesh::polyMesh(
const IOobject&
io,
const bool doInit)
180 time().findInstance(meshDir(),
"points"),
192 time().findInstance(meshDir(),
"faces"),
223 clearedPrimitives_(false),
247 tetBasePtIsPtr_(readTetBasePtIs()),
248 cellTreePtr_(nullptr),
306 globalMeshDataPtr_(nullptr),
308 topoChanging_(false),
309 storeOldCellCentres_(false),
311 oldPointsPtr_(nullptr),
312 oldCellCentresPtr_(nullptr)
343 <<
"Missing mesh boundary on one or more domains" <<
endl;
349 <<
"No points in mesh" <<
endl;
354 <<
"No cells in mesh" <<
endl;
373 boundary_.updateMesh();
376 boundary_.calcGeometry();
385 Foam::polyMesh::polyMesh
449 clearedPrimitives_(false),
464 bounds_(points_, syncPar),
468 tetBasePtIsPtr_(nullptr),
469 cellTreePtr_(nullptr),
512 globalMeshDataPtr_(nullptr),
514 topoChanging_(false),
515 storeOldCellCentres_(false),
517 oldPointsPtr_(nullptr),
518 oldCellCentresPtr_(nullptr)
523 const face& curFace = faces_[facei];
525 if (
min(curFace) < 0 ||
max(curFace) > points_.
size())
528 <<
"Face " << facei <<
"contains vertex labels out of range: " 529 << curFace <<
" Max point index = " << points_.
size()
539 Foam::polyMesh::polyMesh
602 clearedPrimitives_(false),
617 bounds_(points_, syncPar),
621 tetBasePtIsPtr_(nullptr),
622 cellTreePtr_(nullptr),
665 globalMeshDataPtr_(nullptr),
667 topoChanging_(false),
668 storeOldCellCentres_(false),
670 oldPointsPtr_(nullptr),
671 oldCellCentresPtr_(nullptr)
676 const face& curFace = faces_[facei];
678 if (
min(curFace) < 0 ||
max(curFace) > points_.
size())
681 <<
"Face " << facei <<
"contains vertex labels out of range: " 682 << curFace <<
" Max point index = " << points_.
size()
693 const cell& curCell = cLst[celli];
695 if (
min(curCell) < 0 ||
max(curCell) > faces_.
size())
698 <<
"Cell " << celli <<
"contains face labels out of range: " 699 << curCell <<
" Max face index = " << faces_.
size()
723 const bool validBoundary
727 clearAddressing(
true);
733 points_.transfer(*
points);
734 bounds_ =
boundBox(points_, validBoundary);
739 faces_.transfer(*faces);
744 owner_.transfer(*owner);
749 neighbour_.transfer(*neighbour);
756 boundary_[patchi] = polyPatch
773 const face& curFace = faces_[facei];
775 if (
min(curFace) < 0 ||
max(curFace) > points_.size())
778 <<
"Face " << facei <<
" contains vertex labels out of range: " 779 << curFace <<
" Max point index = " << points_.size()
797 boundary_.updateMesh();
800 boundary_.calcGeometry();
806 <<
"No points or no cells in mesh" <<
endl;
835 return parent().dbDir();
844 return dbDir()/meshSubDir;
856 return points_.instance();
862 return faces_.instance();
868 if (geometricD_.x() == 0)
885 if (solutionD_.x() == 0)
902 if (!tetBasePtIsPtr_)
907 <<
"Forcing storage of base points." 911 tetBasePtIsPtr_.reset
929 return *tetBasePtIsPtr_;
940 treeBoundBox overallBb(
points());
941 overallBb.inflate(
rndGen, 1
e-4, ROOTVSMALL);
945 new indexedOctree<treeDataCell>
961 return *cellTreePtr_;
968 const bool validBoundary
971 if (boundaryMesh().size())
974 <<
"boundary already exists" 982 boundary_.transfer(plist);
988 globalMeshDataPtr_.clear();
993 boundary_.updateMesh();
996 boundary_.calcGeometry();
998 boundary_.checkDefinition();
1005 PtrList<pointZone>&& pz,
1006 PtrList<faceZone>&& fz,
1007 PtrList<cellZone>&& cz
1010 if (pointZones_.size() || faceZones_.size() || cellZones_.size())
1013 <<
"point, face or cell zone already exists" 1020 pointZones_.clear();
1021 pointZones_.transfer(pz);
1029 faceZones_.transfer(fz);
1037 cellZones_.transfer(cz);
1045 const List<polyPatch*>&
p,
1046 const bool validBoundary
1052 addPatches(plist, validBoundary);
1075 if (clearedPrimitives_)
1078 <<
"points deallocated" 1088 return io.upToDate(points_);
1094 io.eventNo() = points_.eventNo()+1;
1100 if (clearedPrimitives_)
1103 <<
"faces deallocated" 1137 oldPointsPtr_.reset(
new pointField(points_));
1138 curMotionTimeIndex_ = time().timeIndex();
1141 return *oldPointsPtr_;
1147 storeOldCellCentres_ =
true;
1151 return cellCentres();
1154 if (!oldCellCentresPtr_)
1159 return *oldCellCentresPtr_;
1166 <<
"Moving points for time " << time().value()
1167 <<
" index " << time().timeIndex() <<
endl;
1169 if (newPoints.size() != points_.size())
1172 <<
"Size of newPoints " << newPoints.size()
1173 <<
" does not correspond to current mesh points size " 1182 if (curMotionTimeIndex_ != time().
timeIndex())
1186 Info<<
"tmp<scalarField> polyMesh::movePoints(const pointField&) : " 1187 <<
" Storing current points for time " << time().value()
1188 <<
" index " << time().timeIndex() <<
endl;
1191 if (storeOldCellCentres_)
1193 oldCellCentresPtr_.clear();
1194 oldCellCentresPtr_.reset(
new pointField(cellCentres()));
1198 oldPointsPtr_.clear();
1199 oldPointsPtr_.reset(
new pointField(points_));
1200 curMotionTimeIndex_ = time().timeIndex();
1203 points_ = newPoints;
1205 bool moveError =
false;
1209 if (checkMeshMotion(points_,
true))
1214 <<
"Moving the mesh with given points will " 1215 <<
"invalidate the mesh." <<
nl 1216 <<
"Mesh motion should not be executed." <<
endl;
1221 points_.instance() = time().timeName();
1222 points_.eventNo() = getEvent();
1224 if (tetBasePtIsPtr_)
1227 tetBasePtIsPtr_->instance() = time().timeName();
1228 tetBasePtIsPtr_->eventNo() = getEvent();
1244 if (globalMeshDataPtr_)
1246 globalMeshDataPtr_->movePoints(points_);
1251 bounds_ = boundBox(points_);
1252 boundary_.movePoints(points_);
1254 pointZones_.movePoints(points_);
1255 faceZones_.movePoints(points_);
1256 cellZones_.movePoints(points_);
1267 cellTreePtr_.clear();
1277 meshObject::movePoints<polyMesh>(*this);
1278 meshObject::movePoints<pointMesh>(*this);
1280 const_cast<Time&
>(time()).functionObjects().movePoints(*
this);
1283 if (
debug && moveError)
1294 curMotionTimeIndex_ = 0;
1295 oldPointsPtr_.clear();
1296 oldCellCentresPtr_.clear();
1302 if (!globalMeshDataPtr_)
1306 Pout<<
"polyMesh::globalData() const : " 1307 <<
"Constructing parallelData from processor topology" 1314 return *globalMeshDataPtr_;
1332 fileName meshFilesPath = thisDb().time().
path()/instanceDir/meshDir();
1334 rm(meshFilesPath/
"points");
1335 rm(meshFilesPath/
"faces");
1336 rm(meshFilesPath/
"owner");
1337 rm(meshFilesPath/
"neighbour");
1338 rm(meshFilesPath/
"cells");
1339 rm(meshFilesPath/
"boundary");
1340 rm(meshFilesPath/
"pointZones");
1341 rm(meshFilesPath/
"faceZones");
1342 rm(meshFilesPath/
"cellZones");
1343 rm(meshFilesPath/
"meshModifiers");
1344 rm(meshFilesPath/
"parallelData");
1349 rmDir(meshFilesPath/
"sets");
1356 removeFiles(instance());
1375 celli =
tree.findInside(
p);
1380 findTetFacePt(celli,
p, tetFacei, tetPti);
1393 const polyMesh&
mesh = *
this;
1396 tetFacei = tet.face();
1397 tetPti = tet.tetPt();
1405 const cellDecomposition decompMode
1416 case FACE_CENTRE_TRIS:
1420 const cell& cFaces =
cells()[celli];
1424 label facei = cFaces[cFacei];
1425 const face&
f = faces_[facei];
1426 const point& fc = faceCentres()[facei];
1427 bool isOwn = (owner_[facei] == celli);
1437 nextPointi =
f.nextLabel(fp);
1441 pointi =
f.nextLabel(fp);
1452 vector proj =
p - faceTri.centre();
1454 if ((faceTri.areaNormal() & proj) > 0)
1464 case FACE_DIAG_TRIS:
1468 const cell& cFaces =
cells()[celli];
1472 label facei = cFaces[cFacei];
1473 const face&
f = faces_[facei];
1475 for (label tetPti = 1; tetPti <
f.
size() - 1; tetPti++)
1478 tetIndices faceTetIs(celli, facei, tetPti);
1482 vector proj =
p - faceTri.centre();
1484 if ((faceTri.areaNormal() & proj) > 0)
1500 findTetFacePt(celli,
p, tetFacei, tetPti);
1502 return tetFacei != -1;
1514 const cellDecomposition decompMode
1520 && (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
1529 (void)tetBasePtIs();
1537 if (decompMode == CELL_TETS)
1546 findCellFacePt(
p, celli, tetFacei, tetPti);
1561 (void)tetBasePtIs();
1565 label celli = findNearestCell(
p);
1568 if (pointInCell(
p, celli, decompMode))
1576 for (label celli = 0; celli < nCells(); celli++)
1578 if (pointInCell(
p, celli, decompMode))
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
void size(const label n)
Older name for setAddressableSize.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
A class for handling file names.
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
Find a suitable base point for each face for decomposition into tets.
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual void movePoints(const pointField &)
Move points.
const fileName & facesInstance() const
Return the current instance directory for faces.
A face is a list of labels corresponding to mesh vertices.
label nPoints() const noexcept
Number of mesh points.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
const word & regionName() const
The mesh region name or word::null if polyMesh::defaultRegion.
virtual void setUpToDatePoints(regIOobject &io) const
Set io to be up-to-date with points.
const word & name() const noexcept
Return the object name.
Cell-face mesh analysis engine.
void resetMotion() const
Reset motion.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
constexpr char nl
The newline '\n' character (0x0a)
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
void resetPrimitives(autoPtr< pointField > &&points, autoPtr< faceList > &&faces, autoPtr< labelList > &&owner, autoPtr< labelList > &&neighbour, const labelUList &patchSizes, const labelUList &patchStarts, const bool validBoundary=true)
Reset mesh primitive data. Assumes all patch info correct.
A bounding box defined in terms of min/max extrema points.
const cellList & cells() const
bool pointInCell(const point &p, label celli) const
Return true if the point is in the cell.
virtual const fileName & dbDir() const
Override the objectRegistry dbDir for a single-region case.
static std::string path(const std::string &str)
Return directory path name (part before last /)
Ignore writing from objectRegistry::writeObject()
bool pointInCell(const point &p, label celli, const cellDecomposition=CELL_TETS) const
Test if point p is in the celli.
void movePoints(const pointField &p, const pointField &oldP)
Move points.
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
bool hasHeaderClass() const noexcept
True if headerClassName() is non-empty (after reading)
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
bool isDir(const fileName &name, const bool followLink=true)
Does the name exist as a DIRECTORY in the file system?
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
const fileName & pointsInstance() const
Return the current instance directory for points.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
void removeFiles() const
Remove all files from mesh instance()
A class for handling words, derived from Foam::string.
const Time & time() const noexcept
Return time registry.
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
#define DebugInFunction
Report an information message using Foam::Info.
static word defaultRegion
Return the default region name.
Tree tree(triangles.begin(), triangles.end())
void addZones(PtrList< pointZone > &&pz, PtrList< faceZone > &&fz, PtrList< cellZone > &&cz)
Add mesh zones.
virtual const pointField & oldPoints() const
Return old points (mesh motion)
virtual const labelList & faceOwner() const
Return face owner.
label comm() const noexcept
Return communicator used for parallel communication.
static const word null
An empty word.
const globalMeshData & globalData() const
Return parallel info.
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
virtual const faceList & faces() const
Return raw faces.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
bool rmDir(const fileName &directory, const bool silent=false, const bool emptyOnly=false)
Remove a directory and its contents recursively,.
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross-product operators.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
CompactIOList< cell, label > cellCompactIOList
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
virtual bool upToDatePoints(const regIOobject &io) const
Return true if io is up-to-date with points.
bool empty() const noexcept
True if the list is empty (ie, size() is zero)
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
virtual const pointField & oldCellCentres() const
Return old cellCentres (mesh motion)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
IOobject(const IOobject &)=default
Copy construct.
fileName meshDir() const
Return the local mesh directory (dbDir()/meshSubDir)
Automatically write from objectRegistry::writeObject()
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
triangle< point, const point & > triPointRef
A triangle using referred points.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
virtual ~polyMesh()
Destructor.
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
messageStream Info
Information stream (stdout output on master, null elsewhere)
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
SubField< vector > subField
Declare type of subField.
virtual const fileName & dbDir() const
Local directory path of this objectRegistry relative to the time.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
Mesh consisting of general polyhedral cells.
A subset of mesh faces organised as a primitive patch.
Registry of regIOobjects.
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
virtual bool write(const bool valid=true) const
Write using setting from DB.
void addPatches(polyPatchList &plist, const bool validBoundary=true)
Add boundary patches.
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
Defines the attributes of an object for which implicit objectRegistry management is supported...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< cell > cellList
A List of cells.
Inter-processor communications stream.
IOList< label > labelIOList
Label container classes.
bool rm(const fileName &file)
Remove a file (or its gz equivalent), returning true if successful.
#define InfoInFunction
Report an information message using Foam::Info.
static constexpr const zero Zero
Global zero (0)