47 const label faceBasePtI
55 const point& tetBasePt = pPts[
f[faceBasePtI]];
57 scalar thisBaseMinTetQuality = VGREAT;
59 for (label tetPtI = 1; tetPtI <
f.size() - 1; tetPtI++)
61 label facePtI = (tetPtI + faceBasePtI) %
f.size();
62 label otherFacePtI =
f.fcIndex(facePtI);
70 ptBI =
f[otherFacePtI];
74 ptAI =
f[otherFacePtI];
78 const point& pA = pPts[ptAI];
79 const point& pB = pPts[ptBI];
83 scalar tetQuality = tet.quality();
85 if (tetQuality < thisBaseMinTetQuality)
87 thisBaseMinTetQuality = tetQuality;
90 return thisBaseMinTetQuality;
98 const label faceBasePtI
101 const scalar ownQuality =
113 const scalar neiQuality =
123 if (neiQuality < ownQuality)
135 const polyMesh&
mesh,
148 label oCI = pOwner[fI];
150 const point& oCc = pC[oCI];
154 scalar ownQuality = minQuality(
mesh, oCc, fI,
true, faceBasePtI);
155 scalar neiQuality = minQuality(
mesh, nCc, fI,
false, faceBasePtI);
157 if (
min(ownQuality, neiQuality) > tol)
177 return findSharedBasePoint
202 label cI = pOwner[fI];
204 bool own = (pOwner[fI] == cI);
206 const point& cC = pC[cI];
210 scalar quality = minQuality(
mesh, cC, fI, own, faceBasePtI);
241 for (label fI = 0; fI < nInternalFaces; ++fI)
243 tetBasePtIs[fI] = findSharedBasePoint(
mesh, fI, tol, report);
248 for (label facei = nInternalFaces; facei <
mesh.
nFaces(); ++facei)
250 neighbourCellCentres[facei - nInternalFaces] = pC[pOwner[facei]];
257 SubList<label> boundaryFaceTetBasePtIs
266 label fI = nInternalFaces, bFI = 0;
275 const coupledPolyPatch& cpp =
276 refCast<const coupledPolyPatch>(
patches[patchi]);
280 boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
284 neighbourCellCentres[bFI],
293 boundaryFaceTetBasePtIs[bFI] = -2;
298 boundaryFaceTetBasePtIs[bFI] = findBasePoint
314 boundaryFaceTetBasePtIs,
320 label fI = nInternalFaces, bFI = 0;
325 label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
327 if (bFTetBasePtI == -2)
330 <<
"Coupled face base point exchange failure for face " 335 if (bFTetBasePtI < 1)
347 const coupledPolyPatch& cpp =
348 refCast<const coupledPolyPatch>(
patches[patchi]);
385 const polyMesh&
mesh,
412 label nErrorTets = 0;
416 const face&
f = fcs[facei];
423 p[
f.nextLabel(fPtI)],
432 setPtr->insert(facei);
443 const face&
f = fcs[facei];
450 p[
f.nextLabel(fPtI)],
459 setPtr->insert(facei);
467 if (findSharedBasePoint(
mesh, facei, tol, report) == -1)
471 setPtr->insert(facei);
497 setPtr->insert(facei);
505 if (findBasePoint(
mesh, facei, tol, report) == -1)
509 setPtr->insert(facei);
518 reduce(nErrorTets, sumOp<label>());
524 Info<<
" ***Error in face tets: " 525 << nErrorTets <<
" faces with low quality or negative volume " 526 <<
"decomposition tets." <<
endl;
543 const polyMesh&
mesh,
552 label nTets =
f.
size() - 2;
554 List<tetIndices> faceTets(nTets);
556 for (label tetPti = 1; tetPti <
f.
size() - 1; ++tetPti)
558 faceTets[tetPti - 1] = tetIndices(celli, facei, tetPti);
567 const polyMesh&
mesh,
574 const cell& thisCell = pCells[celli];
578 for (
const label facei : thisCell)
580 nTets +=
pFaces[facei].size() - 2;
583 DynamicList<tetIndices> cellTets(nTets);
585 for (
const label facei : thisCell)
587 cellTets.append(faceTetIndices(
mesh, facei, celli));
596 const polyMesh&
mesh,
604 const cell& thisCell = pCells[celli];
606 tetIndices tetContainingPt;
609 for (
const label facei : thisCell)
613 for (label tetPti = 1; tetPti <
f.
size() - 1; ++tetPti)
616 tetIndices faceTetIs(celli, facei, tetPti);
619 if (faceTetIs.tet(
mesh).inside(pt))
621 tetContainingPt = faceTetIs;
626 if (tetContainingPt.cell() != -1)
632 return tetContainingPt;
638 const polyMesh&
mesh,
656 forAll(tetBasePtIs, facei)
658 if (tetBasePtIs[facei] == -1)
660 problemCells.set(faceOwn[facei]);
664 problemCells.set(faceNei[facei]);
677 Map<label> pointCount;
680 for (
const label celli : problemCells)
684 for (
const label facei :
cells[celli])
686 for (
const label pointi : faces[facei])
688 ++pointCount(pointi);
697 <<
"point:" << iter.key()
699 <<
" only used by one face" <<
nl 702 else if (iter.val() == 2)
704 problemPoints.set(iter.key());
715 forAll(tetBasePtIs, facei)
719 problemCells.test(faceOwn[facei])
723 const face&
f = faces[facei];
727 const label fp0 = tetBasePtIs[facei] < 0 ? 0 : tetBasePtIs[facei];
732 && !problemPoints.test(
f.
fcValue(fp0))
740 scalar maxQ = -GREAT;
747 && !problemPoints.test(
f.
fcValue(fp))
750 const scalar q = minQuality(
mesh, facei, fp);
762 tetBasePtIs[facei] = maxFp;
770 tetBasePtIs[facei] = 0;
781 Pout<<
"Adapted starting point of triangulation on " 782 << nAdapted <<
" faces." <<
endl;
void size(const label n)
Older name for setAddressableSize.
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
Find a suitable base point for each face for decomposition into tets.
List< cell > cellList
List of cell.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
static scalar minQuality(const polyMesh &mesh, const point &cC, const label fI, const bool isOwner, const label faceBasePtI)
Given a face and cc and starting index for triangulation determine the worst tet quality.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const labelList & patchID() const
Per boundary face label the patch index.
const cellList & cells() const
label nFaces() const noexcept
Number of mesh faces.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
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.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< face > faceList
List of faces.
static bool checkFaceTets(const polyMesh &mesh, scalar tol=minTetQuality, const bool report=false, labelHashSet *setPtr=nullptr)
Check face-decomposition tet volume.
vectorField pointField
pointField is a vectorField.
const T & fcValue(const label i) const
Return forward circular value (ie, next value in the list)
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
static labelList adjustTetBasePtIs(const polyMesh &mesh, const bool report=false)
Return an adjusted list of tet base points.
static label findSharedBasePoint(const polyMesh &mesh, label fI, const point &nCc, scalar tol, bool report=false)
Find the first suitable base point to use for a minimum.
virtual const labelList & faceOwner() const
Return face owner.
label nInternalFaces() const noexcept
Number of internal faces.
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.
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
const vectorField & faceCentres() const
vector point
Point is a vector.
const polyBoundaryMesh & patches
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
const T & rcValue(const label i) const
Return reverse circular value (ie, previous value in the list)
static const scalar minTetQuality
Minimum tetrahedron quality.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
forAllConstIters(mixture.phases(), phase)
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.