33 const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
45 void Foam::particle::stationaryTetGeometry
53 const triFace triIs(currentTetIndices().faceTriIs(mesh_));
59 vertex1 = pts[triIs[1]];
60 vertex2 = pts[triIs[2]];
66 vector centre, base, vertex1, vertex2;
67 stationaryTetGeometry(centre, base, vertex1, vertex2);
73 void Foam::particle::stationaryTetReverseTransform
82 const vector ab = A.b() - A.a();
83 const vector ac = A.c() - A.a();
84 const vector ad = A.d() - A.a();
85 const vector bc = A.c() - A.b();
86 const vector bd = A.d() - A.b();
90 detA = ab & (ac ^ ad);
102 void Foam::particle::movingTetGeometry
104 const scalar fraction,
105 Pair<vector>& centre,
107 Pair<vector>& vertex1,
108 Pair<vector>& vertex2
111 const triFace triIs(currentTetIndices().faceTriIs(mesh_));
119 const vector ccOld = mesh_.cells()[celli_].
centre(ptsOld, mesh_.faces());
120 const vector ccNew = mesh_.cells()[celli_].
centre(ptsNew, mesh_.faces());
125 const Pair<scalar> s = stepFractionSpan();
126 const scalar f0 = s[0] + stepFraction_*s[1],
f1 = fraction*s[1];
128 centre[0] = ccOld + f0*(ccNew - ccOld);
129 base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
130 vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
131 vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
133 centre[1] = f1*(ccNew - ccOld);
134 base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
135 vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
136 vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
142 const scalar fraction
145 Pair<vector> centre, base, vertex1, vertex2;
146 movingTetGeometry(fraction, centre, base, vertex1, vertex2);
149 Pair<barycentricTensor>
157 void Foam::particle::movingTetReverseTransform
159 const scalar fraction,
160 Pair<vector>& centre,
161 FixedList<scalar, 4>& detA,
162 FixedList<barycentricTensor, 3>& T
165 Pair<barycentricTensor> A = movingTetTransform(fraction);
167 const Pair<vector> ab(A[0].
b() - A[0].a(), A[1].
b() - A[1].a());
168 const Pair<vector> ac(A[0].
c() - A[0].a(), A[1].
c() - A[1].a());
169 const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
170 const Pair<vector> bc(A[0].
c() - A[0].
b(), A[1].
c() - A[1].
b());
171 const Pair<vector> bd(A[0].d() - A[0].
b(), A[1].d() - A[1].
b());
173 centre[0] = A[0].a();
174 centre[1] = A[1].a();
176 detA[0] = ab[0] & (ac[0] ^ ad[0]);
178 (ab[1] & (ac[0] ^ ad[0]))
179 + (ab[0] & (ac[1] ^ ad[0]))
180 + (ab[0] & (ac[0] ^ ad[1]));
182 (ab[0] & (ac[1] ^ ad[1]))
183 + (ab[1] & (ac[0] ^ ad[1]))
184 + (ab[1] & (ac[1] ^ ad[0]));
185 detA[3] = ab[1] & (ac[1] ^ ad[1]);
196 (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
197 (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
198 (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
199 (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
211 void Foam::particle::reflect()
213 Swap(coordinates_.c(), coordinates_.d());
217 void Foam::particle::rotate(
const bool reverse)
221 scalar temp = coordinates_.b();
222 coordinates_.b() = coordinates_.c();
223 coordinates_.c() = coordinates_.d();
224 coordinates_.d() = temp;
228 scalar temp = coordinates_.d();
229 coordinates_.d() = coordinates_.c();
230 coordinates_.c() = coordinates_.b();
231 coordinates_.b() = temp;
236 void Foam::particle::changeTet(
const label tetTriI)
238 const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
240 const label firstTetPtI = 1;
241 const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
247 else if (tetTriI == 2)
251 if (tetPti_ == lastTetPtI)
263 if (tetPti_ == firstTetPtI)
274 else if (tetTriI == 3)
278 if (tetPti_ == firstTetPtI)
290 if (tetPti_ == lastTetPtI)
304 <<
"Changing tet without changing cell should only happen when the " 305 <<
"track is on triangle 1, 2 or 3." 311 void Foam::particle::changeFace(
const label tetTriI)
314 const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
320 sharedEdge = edge(triOldIs[1], triOldIs[2]);
322 else if (tetTriI == 2)
324 sharedEdge = edge(triOldIs[2], triOldIs[0]);
326 else if (tetTriI == 3)
328 sharedEdge = edge(triOldIs[0], triOldIs[1]);
333 <<
"Changing face without changing cell should only happen when the" 334 <<
" track is on triangle 1, 2 or 3." 337 sharedEdge = edge(-1, -1);
343 forAll(mesh_.cells()[celli_], cellFaceI)
345 const label newFaceI = mesh_.cells()[celli_][cellFaceI];
346 const class face& newFace = mesh_.faces()[newFaceI];
349 if (tetFacei_ == newFaceI)
357 for (; edgeI < newFace.size(); ++ edgeI)
359 edgeComp =
edge::compare(sharedEdge, newFace.faceEdge(edgeI));
374 const label newBaseI =
max(0, mesh_.tetBasePtIs()[newFaceI]);
375 edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
380 edgeI =
min(
max(1, edgeI), newFace.size() - 2);
383 tetFacei_ = newFaceI;
393 <<
"The search for an edge-connected face and tet-point failed." 398 if (sharedEdge.otherVertex(triOldIs[1]) == -1)
402 else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
408 const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
414 if (sharedEdge.otherVertex(triNewIs[1]) == -1)
418 else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
425 void Foam::particle::changeCell()
428 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
429 const bool isOwner = celli_ == ownerCellI;
430 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
437 void Foam::particle::locate
442 const bool boundaryFail,
443 const string boundaryMsg
451 celli_ = mesh_.cellTree().findInside(position);
456 <<
"Cell not found for particle position " << position <<
"." 462 tetFacei_ = mesh_.cells()[celli_][0];
467 track(position - mesh_.cellCentres()[celli_], 0);
484 if (direction ==
nullptr)
486 const polyPatch& p = mesh_.boundaryMesh()[patch()];
487 direction = &p.faceNormals()[p.whichFace(facei_)];
490 const vector n = *direction/
mag(*direction);
491 const vector s = position - mesh_.cellCentres()[celli_];
496 tetFacei_ = mesh_.cells()[celli_][0];
503 static label nWarnings = 0;
504 static const label maxNWarnings = 100;
505 if (nWarnings < maxNWarnings)
510 if (nWarnings == maxNWarnings)
513 <<
"Suppressing any further warnings about particles being " 514 <<
"located outside of the mesh." <<
endl;
528 const label tetFacei,
533 coordinates_(coordinates),
540 origId_(getNewParticleID())
552 coordinates_(- VGREAT, - VGREAT, - VGREAT, - VGREAT),
559 origId_(getNewParticleID())
567 "Particle initialised with a location outside of the mesh." 575 coordinates_(p.coordinates_),
577 tetFacei_(p.tetFacei_),
580 stepFraction_(p.stepFraction_),
581 origProc_(p.origProc_),
589 coordinates_(p.coordinates_),
591 tetFacei_(p.tetFacei_),
594 stepFraction_(p.stepFraction_),
595 origProc_(p.origProc_),
604 const vector& displacement,
605 const scalar fraction
623 const vector& displacement,
624 const scalar fraction
635 f *=
trackToTri(f*displacement, f*fraction, tetTriI);
642 else if (tetTriI == 0)
659 const vector& displacement,
660 const scalar fraction,
665 const vector x1 = displacement;
671 <<
"Tracking from " << x0 <<
" to " << x0 + x1 <<
endl;
678 stationaryTetReverseTransform(centre, detA, T);
683 stationaryTetGeometry(o, b, v1, v2);
684 Info<<
"Tet points o=" << o <<
", b=" << b
685 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl 686 <<
"Tet determinant = " << detA <<
endl 687 <<
"Start local coordinates = " << y0 <<
endl;
691 const scalar
f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
698 Info<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
703 scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
704 for (
label i = 0; i < 4; ++ i)
706 if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
708 scalar
mu = - y0[i]/Tx1[i];
712 Info<<
"Hit on tet face " << i <<
" at local coordinate " 713 << y0 + mu*Tx1 <<
", " << mu*detA*100 <<
"\% of the " 714 <<
"way along the track" <<
endl;
717 if (0 <= mu && mu < muH)
729 for (
label i = 0; i < 4; ++ i)
748 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
752 Info<<
"Track hit no tet faces" <<
endl;
754 Info<<
"End local coordinates = " << yH <<
endl 756 <<
"Tracking displacement = " <<
position() - x0 <<
endl 757 << muH*detA*100 <<
"\% of the step from " << stepFraction_ <<
" to " 758 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
762 stepFraction_ += fraction*muH*detA;
764 return iH != -1 ? 1 - muH*detA : 0;
770 const vector& displacement,
771 const scalar fraction,
776 const vector x1 = displacement;
782 <<
"Tracking from " << x0 <<
" to " << x0 + x1 <<
endl;
789 movingTetReverseTransform(fraction, centre, detA, T);
792 const scalar
f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
795 const vector x0Rel = x0 - centre[0];
796 const vector x1Rel = f*x1 - centre[1];
799 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
802 ((x1Rel & T[2]) + detA[3]*yC)*
sqr(detA[0]);
804 ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
806 ((x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC);
811 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
816 scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
817 for (
label i = 0; i < 4; ++ i)
821 for (
label j = 0; j < 3; ++ j)
825 if (0 <= mu[j] && mu[j] < muH)
837 hitEqn[0].value(muH),
838 hitEqn[1].value(muH),
839 hitEqn[2].value(muH),
849 const scalar detAH = detAEqn.
value(muH);
850 if (!std::isnormal(detAH))
853 <<
"A moving tet collapsed onto a particle. This is not supported. " 854 <<
"The mesh is too poor, or the motion too severe, for particle " 860 for (
label i = 0; i < 4; ++ i)
876 stepFraction_ += fraction*muH*detA[0];
882 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
886 Info<<
"Track hit no tet faces" <<
endl;
888 Info<<
"End local coordinates = " << yH <<
endl 892 return iH != -1 ? 1 - muH*detA[0] : 0;
898 const vector& displacement,
899 const scalar fraction,
918 bool isConstrained =
false;
921 if (dirs[dirI] == -1)
923 isConstrained =
true;
942 <<
"Patch data was requested for a particle that isn't on a patch" 949 movingTetGeometry(1, centre, base, vertex1, vertex2);
951 n =
triPointRef(base[0], vertex1[0], vertex2[0]).normal();
957 coordinates_.
b()*base[1]
958 + coordinates_.
c()*vertex1[1]
959 + coordinates_.
d()*vertex2[1];
967 vector centre, base, vertex1, vertex2;
968 stationaryTetGeometry(centre, base, vertex1, vertex2);
1011 if (transform.
hasR())
1021 const vector pos(coordinates_.
b(), coordinates_.
c(), coordinates_.
d());
1025 tetFacei_ = mesh_.
cells()[celli_][0];
1041 movingTetReverseTransform(0, centre, detA, T);
1042 coordinates_ += (
pos - centre[0]) & T[0]/detA[0];
1049 stationaryTetReverseTransform(centre, detA, T);
1050 coordinates_ += (
pos - centre) & T/detA;
1058 const label procCell,
1059 const label procTetFace
1068 (mesh_.
faceOwner()[tetFacei_] == celli_)
1069 == (procMesh.
faceOwner()[procTetFace] == procCell)
1076 return procMesh.
faces()[procTetFace].
size() - 1 - tetPti_;
1093 "Particle mapped to a location outside of the mesh."
label origProc() const
Return the originating processor ID.
#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.
bool moving() const
Is mesh moving.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A 1D vector of objects of type <T> with a fixed size <Size>.
bool onBoundaryFace() const
Is the particle on a boundary face?
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label procTetPt(const polyMesh &procMesh, const label procCell, const label procTetFace) const
Return the tet point approproate for decomposition or reconstruction.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but also stops on internal faces.
void size(const label)
Override size to be inconsistent with allocated storage.
Tensor< Cmpt > T() const
Return transpose.
static int compare(const edge &, const edge &)
Compare edges.
scalar trackToTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToFace, but also stops on tet triangles. On.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar y0(const dimensionedScalar &ds)
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
const cellList & cells() const
Vector< scalar > vector
A scalar version of the templated Vector.
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
void replace(const direction, const Cmpt &)
void prepareForInteractionListReferral(const vectorTensorTransform &transform)
Break the topology and store the particle position so that the.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
dimensionedScalar pos(const dimensionedScalar &ds)
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
vectorField pointField
pointField is a vectorField.
An ordered pair of two objects of type <T> with first() and second() elements.
void patchData(vector &n, vector &U) const
Get the normal and velocity of the current patch location.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
scalar deltaTValue() const
Return time step value.
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. The particle may still be.
virtual const labelList & faceOwner() const
Return face owner.
virtual scalar wallImpactDistance(const vector &n) const
The nearest distance to a wall that.
virtual const faceList & faces() const
Return raw faces.
label origId() const
Return the particle ID on the originating processor.
void autoMap(const vector &position, const mapPolyMesh &mapper)
Map after a topology change.
bool onFace() const
Is the particle on a face?
const labelList & reverseCellMap() const
Reverse cell map.
const Time & time() const
Return time.
defineTypeNameAndDebug(combustionModel, 0)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
const dimensionedScalar mu
Atomic mass unit.
scalar trackToMovingTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for moving meshes.
#define WarningInFunction
Report a warning using Foam::Warning.
void type(const direction i, const roots::type t)
Set the type of the i-th root.
const dimensionedScalar c
Speed of light in a vacuum.
triangle< point, const point & > triPointRef
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
dimensioned< scalar > mag(const dimensioned< Type > &)
static label particleCount_
Cumulative particle counter - used to provode unique ID.
Field< vector > vectorField
Specialisation of Field<T> for vector.
bool onInternalFace() const
Is the particle on an internal face?
Mesh consisting of general polyhedral cells.
scalar trackToStationaryTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for stationary meshes.
scalar track(const vector &displacement, const scalar fraction)
Track along the displacement for a given fraction of the overall.
bool operator!=(const particle &, const particle &)
void constrainToMeshCentre()
Set the constrained components of the particle position to the.
scalar value(const scalar x) const
Evaluate at x.
vector position() const
Return current particle position.