45 for (
label i = 0; i < 4; i++)
60 if ((lambda > 0.0) && (lambda < 1.0))
74 const label tetPlaneBasePtI,
85 return movingTetLambda
99 const point& base = pPts[tetPlaneBasePtI];
101 scalar lambdaNumerator = (base - from) & n;
102 scalar lambdaDenominator = (to - from) & n;
109 if (
mag(lambdaDenominator) < tol)
111 if (
mag(lambdaNumerator) < tol)
121 if (
mag((to - from)) < tol/
mag(n))
131 lambdaDenominator =
sign(lambdaDenominator)*SMALL;
136 return lambdaNumerator/lambdaDenominator;
146 const label tetPlaneBasePtI,
148 const label tetFaceI,
154 const pointField& oldPPts = mesh_.oldPoints();
157 const point&
b = pPts[tetPlaneBasePtI];
162 const point& b00 = oldPPts[tetPlaneBasePtI];
166 point b0 = b00 + stepFraction_*(b - b00);
172 tetIndices tetIs(cellI, tetFaceI, tetPtI, mesh_);
180 point tet0PtA = tet00.
a() + stepFraction_*(tet.
a() - tet00.
a());
181 point tet0PtB = tet00.
b() + stepFraction_*(tet.
b() - tet00.
b());
182 point tet0PtC = tet00.
c() + stepFraction_*(tet.
c() - tet00.
c());
183 point tet0PtD = tet00.
d() + stepFraction_*(tet.
d() - tet00.
d());
186 tetPointRef tet0(tet0PtA, tet0PtB, tet0PtC, tet0PtD);
224 scalar lambdaNumerator = 0;
225 scalar lambdaDenominator = 0;
234 scalar a = (dP - dB) & dN;
235 scalar b = ((dP - dB) & n0) + (dS & dN);
242 scalar discriminant =
sqr(b) - 4.0*a*
c;
244 if (discriminant < 0)
278 lambdaNumerator = -
c;
279 lambdaDenominator =
b;
286 lambdaNumerator = -(dS & n0);
287 lambdaDenominator = ((dP - dB) & n0);
290 if (
mag(lambdaDenominator) < tol)
292 if (
mag(lambdaNumerator) < tol)
302 if (
mag((to - from)) < tol/
mag(n))
311 lambdaDenominator =
sign(lambdaDenominator)*SMALL;
316 return lambdaNumerator/lambdaDenominator;
332 if (tetBasePtI == -1)
336 "inline void Foam::particle::tetNeighbour(label triI)" 338 <<
"No base point for face " <<
tetFaceI_ <<
", " << f
339 <<
", produces a valid tet decomposition." 453 "Foam::particle::tetNeighbour(label triI)" 455 <<
"Tet tri face index error, can only be 0..3, supplied " 484 label fI = thisCell[cFI];
499 else if (f == pFaces[fI])
531 if (tetBasePtI == -1)
536 "Foam::particle::crossEdgeConnectedFace" 538 "const label& cellI," 544 <<
"No base point for face " << fI <<
", " << f
545 <<
", produces a decomposition that has a minimum " 546 <<
"volume greater than tolerance." 551 eIndex -= tetBasePtI;
555 eIndex = (eIndex + otherFace.
size()) % otherFace.
size();
564 else if (eIndex == otherFace.
size() - 1)
568 tetPtI = otherFace.
size() - 2;
589 WarningIn(
"particle::getNewParticleID() const")
590 <<
"Particle counter has overflowed. This might cause problems" 591 <<
" when reconstructing particle tracks." <<
endl;
702 <<
"cell, tetFace and tetPt search failure at position " 736 <<
" cell, tetFace and tetPt search failure at " 738 <<
" for requested cell " << oldCellI <<
nl 739 <<
" If this is a restart or " 740 "reconstruction/decomposition etc. it is likely that" 741 " the write precision is not sufficient.\n" 742 " Either increase 'writePrecision' or " 743 "set 'writeFormat' to 'binary'" 782 }
while (
tetFaceI_ < 0 && iterNo <= trap);
787 <<
"cell, tetFace and tetPt search failure at position " 793 WarningIn(
"void Foam::particle::initCellFacePt()")
795 <<
" to " << newPosition
799 <<
" (A fraction of " 801 <<
" of the distance to the cell centre)" 802 <<
" because a decomposition tetFace and tetPt " 803 <<
"could not be found." 810 if (debug &&
cellI_ != oldCellI)
812 WarningIn(
"void Foam::particle::initCellFacePt()")
814 <<
" searched for a cell, tetFace and tetPt." <<
nl 819 <<
" This is a different cell to that which was supplied" 820 <<
" (" << oldCellI <<
")." <<
nl dimensionedScalar sqrt(const dimensionedScalar &ds)
scalar currentTime() const
Return the particle current time.
label & cell()
Return current cell particle is in.
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
label & face()
Return current face particle is on otherwise -1.
dimensioned< scalar > mag(const dimensioned< Type > &)
A cell is defined as a list of faces with extra functionality.
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
triPointRef oldFaceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
label getNewParticleID() const
Get unique particle creation id.
label origId() const
Return const access to the particle id on originating processor.
void clear()
Clear the addressed list, i.e. set the size to zero.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void size(const label)
Override size to be inconsistent with allocated storage.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
const cellList & cells() const
dimensionedScalar sign(const dimensionedScalar &ds)
static label particleCount_
Cumulative particle counter - used to provode unique ID.
const vectorField & cellCentres() const
A 1D vector of objects of type <T> with a fixed size <Size>.
label patchFace(const label patchI, const label faceI) const
Which face of this patch is this particle on.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Point centre() const
Return centre (centroid)
label origId_
Local particle id on originating processor.
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
void findTris(const vector &position, DynamicList< label > &faceList, const tetPointRef &tet, const FixedList< vector, 4 > &tetAreas, const FixedList< label, 4 > &tetPlaneBasePtIs, const scalar tol) const
Find the tet tri faces between position and tet centre.
label tetFaceI_
Index of the face that owns the decomposed tet that the.
A face is a list of labels corresponding to mesh vertices.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static const scalar trackingCorrectionTol
Fraction of distance to tet centre to move a particle to.
#define WarningIn(functionName)
Report a warning using Foam::Warning.
scalar deltaTValue() const
Return time step value.
tetPointRef oldTet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
vector normal() const
Return the normal of the tri on tetFaceI_ for the.
vector normal() const
Return vector normal.
const vector & position() const
Return current particle position.
tetPointRef currentTet() const
Return the geometry of the current tet that the.
vector oldNormal() const
Return the normal of the tri on tetFaceI_ for the.
label cellI_
Index of the cell it is in.
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
const Point & a() const
Return vertices.
scalar tetLambda(const vector &from, const vector &to, const label triI, const vector &tetArea, const label tetPlaneBasePtI, const label cellI, const label tetFaceI, const label tetPtI, const scalar tol) const
Find the lambda value for the line to-from across the.
label faceInterpolation() const
Return the index of the face to be used in the interpolation.
const polyMesh & mesh_
Reference to the polyMesh database.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
scalar stepFraction_
Fraction of time-step completed.
const polyMesh & mesh() const
Return the mesh database.
label & tetFace()
Return current tet face particle is in.
label origProc() const
Return const access to the originating processor id.
errorManip< error > abort(error &err)
virtual const labelList & faceOwner() const
Return face owner.
bool softImpact() const
Return the impact model to be used, soft or hard (default).
bool onBoundary() const
Is the particle on the boundary/(or outside the domain)?
label otherFace(const primitiveMesh &, const label cellI, const label faceI, const label edgeI)
Return face on cell using edgeI but not faceI. Throws error.
label & tetPt()
Return current tet face particle is in.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
scalar & stepFraction()
Return the fraction of time-step completed.
label nInternalFaces() const
Mesh consisting of general polyhedral cells.
int edgeDirection(const edge &) const
Return the edge direction on the face.
label patch(const label faceI) const
Which patch is particle on.
label end() const
Return end vertex label.
bool pointInCellBB(const point &p, label celli, scalar inflationFraction=0) const
Return true if the point in the cell bounding box.
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
bool boundaryFace(const label faceI) const
Is this global face a boundary face?
const dimensionedScalar c
Speed of light in a vacuum.
void crossEdgeConnectedFace(const label &cellI, label &tetFaceI, label &tetPtI, const edge &e)
Cross the from the given face across the given edge of the.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
const Time & time() const
Return time.
label start() const
Return start vertex label.
bool internalFace(const label faceI) const
Is this global face an internal face?
virtual const faceList & faces() const
Return raw faces.
scalar movingTetLambda(const vector &from, const vector &to, const label triI, const vector &tetArea, const label tetPlaneBasePtI, const label cellI, const label tetFaceI, const label tetPtI, const scalar tol) const
Find the lambda value for a moving tri face.
label tetPtI_
Index of the point on the face that defines the decomposed.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("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]=cellShape(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]
void initCellFacePt()
Check the stored cell value (setting if necessary) and.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
label origProc_
Originating processor id.
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
void tetNeighbour(label triI)
Modify the tet owner data by crossing triI.
tetIndices currentTetIndices() const
Return the indices of the current tet that the.
label faceI_
Face index if the particle is on a face otherwise -1.
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
vector position_
Position of particle.
const Type & value() const
Return const reference to value.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
static const label labelMax