Go to the documentation of this file.
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 =
105 mesh.cellCentres()[
mesh.faceOwner()[facei]],
111 if (
mesh.isInternalFace(facei))
113 const scalar neiQuality =
117 mesh.cellCentres()[
mesh.faceNeighbour()[facei]],
123 if (neiQuality < ownQuality)
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
181 mesh.cellCentres()[
mesh.faceNeighbour()[fI]],
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);
239 label nInternalFaces =
mesh.nInternalFaces();
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]];
253 syncTools::swapBoundaryFacePositions(
mesh, neighbourCellCentres);
260 mesh.nFaces() - nInternalFaces,
266 label fI = nInternalFaces, bFI = 0;
271 label patchi =
mesh.boundaryMesh().patchID()[bFI];
276 refCast<const coupledPolyPatch>(
patches[patchi]);
280 boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
284 neighbourCellCentres[bFI],
293 boundaryFaceTetBasePtIs[bFI] = -2;
298 boundaryFaceTetBasePtIs[bFI] = findBasePoint
311 syncTools::syncBoundaryFaceList
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 "
331 << fI <<
" at " <<
mesh.faceCentres()[fI]
335 if (bFTetBasePtI < 1)
343 label patchi =
mesh.boundaryMesh().patchID()[bFI];
348 refCast<const coupledPolyPatch>(
patches[patchi]);
374 bFTetBasePtI =
mesh.faces()[fI].size() - bFTetBasePtI;
401 for (label facei =
mesh.nInternalFaces(); facei <
mesh.nFaces(); ++facei)
403 neiCc[facei -
mesh.nInternalFaces()] = cc[own[facei]];
406 syncTools::swapBoundaryFacePositions(
mesh, neiCc);
412 label nErrorTets = 0;
416 const face&
f = fcs[facei];
423 p[
f.nextLabel(fPtI)],
440 if (
mesh.isInternalFace(facei))
443 const face&
f = fcs[facei];
450 p[
f.nextLabel(fPtI)],
467 if (findSharedBasePoint(
mesh, facei, tol, report) == -1)
479 label patchi =
patches.patchID()[facei -
mesh.nInternalFaces()];
489 neiCc[facei -
mesh.nInternalFaces()],
505 if (findBasePoint(
mesh, facei, tol, report) == -1)
524 Info<<
" ***Error in face tets: "
525 << nErrorTets <<
" faces with low quality or negative volume "
526 <<
"decomposition tets." <<
endl;
552 label nTets =
f.size() - 2;
556 for (label tetPti = 1; tetPti <
f.size() - 1; ++tetPti)
558 faceTets[tetPti - 1] =
tetIndices(celli, facei, tetPti);
574 const cell& thisCell = pCells[celli];
578 for (
const label facei : thisCell)
580 nTets +=
pFaces[facei].size() - 2;
585 for (
const label facei : thisCell)
587 cellTets.
append(faceTetIndices(
mesh, facei, celli));
604 const cell& thisCell = pCells[celli];
609 for (
const label facei : thisCell)
613 for (label tetPti = 1; tetPti <
f.size() - 1; ++tetPti)
621 tetContainingPt = faceTetIs;
626 if (tetContainingPt.
cell() != -1)
632 return tetContainingPt;
656 forAll(tetBasePtIs, facei)
658 if (tetBasePtIs[facei] == -1)
660 problemCells.set(faceOwn[facei]);
662 if (
mesh.isInternalFace(facei))
664 problemCells.set(faceNei[facei]);
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()
698 <<
" at:" <<
mesh.points()[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])
720 || (
mesh.isInternalFace(facei) && problemCells.test(faceNei[facei]))
723 const face&
f = faces[facei];
727 const label fp0 = tetBasePtIs[facei] < 0 ? 0 : tetBasePtIs[facei];
731 !problemPoints.test(
f.rcValue(fp0))
732 && !problemPoints.test(
f.fcValue(fp0))
740 scalar maxQ = -GREAT;
746 !problemPoints.test(
f.rcValue(fp))
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;
static bool checkFaceTets(const polyMesh &mesh, scalar tol=minTetQuality, const bool report=false, labelHashSet *setPtr=nullptr)
Check face-decomposition tet volume.
bool inside(const point &pt) const
Return true if point is inside tetrahedron.
static const scalar minTetQuality
Minimum tetrahedron quality.
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
label cell() const
Return the cell.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
A List obtained as a section of another List.
Ostream & endl(Ostream &os)
Add newline and flush stream.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
virtual bool owner() const =0
Does this side own the patch ?
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Mesh consisting of general polyhedral cells.
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
#define forAll(list, i)
Loop across all elements in list.
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.
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
messageStream Info
Information stream (stdout output on master, null elsewhere)
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
reduce(hasMovingMesh, orOp< bool >())
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
errorManip< error > abort(error &err)
errorManipArg< error, int > exit(error &err, const int errNo=1)
static labelList adjustTetBasePtIs(const polyMesh &mesh, const bool report=false)
Return an adjusted list of tet base points.
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
forAllConstIters(mixture.phases(), phase)
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
const polyBoundaryMesh & patches
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
A face is a list of labels corresponding to mesh vertices.
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
static scalar minQuality(const polyMesh &mesh, const point &cC, const label fI, const bool isOwner, const label faceBasePtI)
A cell is defined as a list of faces with extra functionality.
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]
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.