45 void Foam::particle::stationaryTetReverseTransform
54 const vector ab = A.b() - A.a();
55 const vector ac = A.c() - A.a();
56 const vector ad = A.d() - A.a();
57 const vector bc = A.c() - A.b();
58 const vector bd = A.d() - A.b();
62 detA = ab & (ac ^ ad);
74 void Foam::particle::movingTetReverseTransform
76 const scalar fraction,
78 FixedList<scalar, 4>& detA,
79 FixedList<barycentricTensor, 3>& T
82 Pair<barycentricTensor> A = movingTetTransform(fraction);
84 const Pair<vector> ab(A[0].
b() - A[0].a(), A[1].
b() - A[1].a());
85 const Pair<vector> ac(A[0].
c() - A[0].a(), A[1].
c() - A[1].a());
86 const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
87 const Pair<vector> bc(A[0].
c() - A[0].
b(), A[1].
c() - A[1].
b());
88 const Pair<vector> bd(A[0].d() - A[0].
b(), A[1].d() - A[1].
b());
93 detA[0] = ab[0] & (ac[0] ^ ad[0]);
95 (ab[1] & (ac[0] ^ ad[0]))
96 + (ab[0] & (ac[1] ^ ad[0]))
97 + (ab[0] & (ac[0] ^ ad[1]));
99 (ab[0] & (ac[1] ^ ad[1]))
100 + (ab[1] & (ac[0] ^ ad[1]))
101 + (ab[1] & (ac[1] ^ ad[0]));
102 detA[3] = ab[1] & (ac[1] ^ ad[1]);
113 (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
114 (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
115 (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
116 (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
128 void Foam::particle::reflect()
130 Swap(coordinates_.c(), coordinates_.d());
134 void Foam::particle::rotate(
const bool reverse)
138 scalar temp = coordinates_.b();
139 coordinates_.b() = coordinates_.c();
140 coordinates_.c() = coordinates_.d();
141 coordinates_.d() = temp;
145 scalar temp = coordinates_.d();
146 coordinates_.d() = coordinates_.c();
147 coordinates_.c() = coordinates_.b();
148 coordinates_.b() = temp;
153 void Foam::particle::changeTet(
const label tetTriI)
160 const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
162 const label firstTetPtI = 1;
163 const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
169 else if (tetTriI == 2)
173 if (tetPti_ == lastTetPtI)
185 if (tetPti_ == firstTetPtI)
196 else if (tetTriI == 3)
200 if (tetPti_ == firstTetPtI)
212 if (tetPti_ == lastTetPtI)
226 <<
"Changing tet without changing cell should only happen when the " 227 <<
"track is on triangle 1, 2 or 3." 233 void Foam::particle::changeFace(
const label tetTriI)
241 const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
247 sharedEdge = edge(triOldIs[1], triOldIs[2]);
249 else if (tetTriI == 2)
251 sharedEdge = edge(triOldIs[2], triOldIs[0]);
253 else if (tetTriI == 3)
255 sharedEdge = edge(triOldIs[0], triOldIs[1]);
260 <<
"Changing face without changing cell should only happen when the" 261 <<
" track is on triangle 1, 2 or 3." 264 sharedEdge = edge(-1, -1);
270 forAll(mesh_.cells()[celli_], cellFaceI)
272 const label newFaceI = mesh_.cells()[celli_][cellFaceI];
273 const class face& newFace = mesh_.faces()[newFaceI];
274 const label newOwner = mesh_.faceOwner()[newFaceI];
277 if (tetFacei_ == newFaceI)
285 const label edgeComp = newOwner == celli_ ? -1 : +1;
290 edgeI < newFace.size()
291 &&
edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
296 if (edgeI >= newFace.size())
302 const label newBaseI =
max(0, mesh_.tetBasePtIs()[newFaceI]);
303 edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
308 edgeI =
min(
max(1, edgeI), newFace.size() - 2);
311 tetFacei_ = newFaceI;
321 <<
"The search for an edge-connected face and tet-point failed." 326 if (sharedEdge.otherVertex(triOldIs[1]) == -1)
330 else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
336 const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
342 if (sharedEdge.otherVertex(triNewIs[1]) == -1)
346 else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
353 void Foam::particle::changeCell()
361 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
362 const bool isOwner = celli_ == ownerCellI;
363 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
370 void Foam::particle::changeToMasterPatch()
377 label thisPatch = patch();
379 forAll(mesh_.cells()[celli_], cellFaceI)
382 const label otherFaceI = mesh_.cells()[celli_][cellFaceI];
383 if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
391 const class face& thisFace = mesh_.faces()[facei_];
392 const class face& otherFace = mesh_.faces()[otherFaceI];
395 const label otherPatch =
396 mesh_.boundaryMesh().whichPatch(otherFaceI);
397 if (thisPatch > otherPatch)
400 thisPatch = otherPatch;
409 void Foam::particle::locate
413 const bool boundaryFail,
414 const string boundaryMsg
425 celli = mesh_.cellTree().findInside(position);
430 <<
"Cell not found for particle position " << position <<
"." 436 const vector displacement = position - mesh_.cellCentres()[celli_];
441 const class cell& c = mesh_.cells()[celli_];
442 scalar minF = vGreat;
443 label minTetFacei = -1, minTetPti = -1;
446 const class face& f = mesh_.faces()[c[cellTetFacei]];
447 for (
label tetPti = 1; tetPti < f.size() - 1; ++ tetPti)
450 tetFacei_ = c[cellTetFacei];
455 const scalar f = trackToTri(displacement, 0, tetTriI);
465 minTetFacei = tetFacei_;
474 tetFacei_ = minTetFacei;
478 track(displacement, 0);
491 static label nWarnings = 0;
492 static const label maxNWarnings = 100;
493 if (nWarnings < maxNWarnings)
498 if (nWarnings == maxNWarnings)
501 <<
"Suppressing any further warnings about particles being " 502 <<
"located outside of the mesh." <<
endl;
516 const label tetFacei,
521 coordinates_(coordinates),
530 origId_(getNewParticleID())
542 coordinates_(- vGreat, - vGreat, - vGreat, - vGreat),
551 origId_(getNewParticleID())
558 "Particle initialised with a location outside of the mesh." 566 coordinates_(p.coordinates_),
568 tetFacei_(p.tetFacei_),
571 stepFraction_(p.stepFraction_),
573 nBehind_(p.nBehind_),
574 origProc_(p.origProc_),
582 coordinates_(p.coordinates_),
584 tetFacei_(p.tetFacei_),
587 stepFraction_(p.stepFraction_),
589 nBehind_(p.nBehind_),
590 origProc_(p.origProc_),
599 const vector& displacement,
600 const scalar fraction
623 const vector& displacement,
624 const scalar fraction
632 const scalar f =
trackToFace(displacement, fraction);
645 const vector& displacement,
646 const scalar fraction
661 while (nBehind_ < maxNBehind_)
663 f *=
trackToTri(f*displacement, f*fraction, tetTriI);
670 else if (tetTriI == 0)
684 static label stuckID = -1, stuckProc = -1;
685 if (origId_ != stuckID && origProc_ != stuckProc)
688 <<
"Particle #" << origId_ <<
" got stuck at " <<
position()
693 stuckProc = origProc_;
695 stepFraction_ += f*fraction;
706 const vector& displacement,
707 const scalar fraction,
712 const vector x1 = displacement;
718 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
725 stationaryTetReverseTransform(centre, detA, T);
730 stationaryTetGeometry(o, b, v1, v2);
731 Info<<
"Tet points o=" << o <<
", b=" << b
732 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl 733 <<
"Tet determinant = " << detA <<
endl 734 <<
"Start local coordinates = " << y0 <<
endl;
742 Info<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
747 scalar muH = std::isnormal(detA) && detA <= 0 ? vGreat : 1/detA;
748 for (
label i = 0; i < 4; ++ i)
750 if (Tx1[i] < - detA*small)
752 scalar
mu = - y0[i]/Tx1[i];
756 Info<<
"Hit on tet face " << i <<
" at local coordinate " 757 << y0 + mu*Tx1 <<
", " << mu*detA*100 <<
"% of the " 758 <<
"way along the track" <<
endl;
761 if (0 <= mu && mu < muH)
773 for (
label i = 0; i < 4; ++ i)
792 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
796 Info<<
"Track hit no tet faces" <<
endl;
798 Info<<
"End local coordinates = " << yH <<
endl 800 <<
"Tracking displacement = " <<
position() - x0 <<
endl 801 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to " 802 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
806 stepFraction_ += fraction*muH*detA;
809 if (detA <= 0 || nBehind_ > 0)
811 behind_ += muH*detA*
mag(displacement);
824 return iH != -1 ? 1 - muH*detA : 0;
830 const vector& displacement,
831 const scalar fraction,
836 const vector x1 = displacement;
842 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
849 movingTetReverseTransform(fraction, centre, detA, T);
854 movingTetGeometry(fraction, o, b, v1, v2);
855 Info<<
"Tet points o=" << o[0] <<
", b=" << b[0]
856 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl 857 <<
"Tet determinant = " << detA[0] <<
endl 858 <<
"Start local coordinates = " << y0[0] <<
endl;
862 const vector x0Rel = x0 - centre[0];
863 const vector x1Rel = x1 - centre[1];
866 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
869 ((x1Rel & T[2]) + detA[3]*yC)*
sqr(detA[0]);
871 ((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
873 ((x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC);
878 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
883 for (
label i = 0; i < 4; ++ i)
885 Info<< (i ?
" " :
"Hit equation ") << i <<
" = " 886 << hitEqn[i] <<
endl;
888 Info<<
" DetA equation = " << detA <<
endl;
893 scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? vGreat : 1/detA[0];
894 for (
label i = 0; i < 4; ++ i)
898 for (
label j = 0; j < 3; ++ j)
903 && hitEqn[i].derivative(mu[j]) < - detA[0]*small
910 hitEqn[0].value(mu[j]),
911 hitEqn[1].value(mu[j]),
912 hitEqn[2].value(mu[j]),
913 hitEqn[3].value(mu[j])
915 const scalar detAH = detAEqn.
value(mu[j]);
917 Info<<
"Hit on tet face " << i <<
" at local coordinate " 918 << yH/detAH <<
", " << mu[j]*detA[0]*100 <<
"% of the " 919 <<
"way along the track" <<
endl;
922 if (0 <= mu[j] && mu[j] < muH)
934 hitEqn[0].value(muH),
935 hitEqn[1].value(muH),
936 hitEqn[2].value(muH),
946 const scalar detAH = detAEqn.
value(muH);
947 if (!std::isnormal(detAH))
950 <<
"A moving tet collapsed onto a particle. This is not supported. " 951 <<
"The mesh is too poor, or the motion too severe, for particle " 957 for (
label i = 0; i < 4; ++ i)
973 stepFraction_ += fraction*muH*detA[0];
976 if (detA[0] <= 0 || nBehind_ > 0)
978 behind_ += muH*detA[0]*
mag(displacement);
995 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
999 Info<<
"Track hit no tet faces" <<
endl;
1001 Info<<
"End local coordinates = " << yH <<
endl 1003 <<
"Tracking displacement = " <<
position() - x0 <<
endl 1004 << muH*detA[0]*100 <<
"% of the step from " << stepFraction_
1005 <<
" to " << stepFraction_ + fraction <<
" completed" <<
endl 1009 return iH != -1 ? 1 - muH*detA[0] : 0;
1015 const vector& displacement,
1016 const scalar fraction,
1020 if (mesh_.
moving() && (stepFraction_ != 1 || fraction != 0))
1093 facei_ += ppp.
start();
1098 tetPti_ = mesh_.
faces()[tetFacei_].
size() - 1 - tetPti_;
1129 if (transform.
hasR())
1139 const vector pos(coordinates_.
b(), coordinates_.
c(), coordinates_.
d());
1143 tetFacei_ = mesh_.
cells()[celli_][0];
1154 if (mesh_.
moving() && stepFraction_ != 1)
1159 movingTetReverseTransform(0, centre, detA, T);
1160 coordinates_ += (
pos - centre[0]) & T[0]/detA[0];
1167 stationaryTetReverseTransform(centre, detA, T);
1168 coordinates_ += (
pos - centre) & T/detA;
1176 const label procCell,
1177 const label procTetFace
1186 (mesh_.
faceOwner()[tetFacei_] == celli_)
1187 == (procMesh.
faceOwner()[procTetFace] == procCell)
1194 return procMesh.
faces()[procTetFace].
size() - 1 - tetPti_;
1210 "Particle mapped to a location outside of the mesh."
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
label origProc() const
Return the originating processor ID.
void correctAfterParallelTransfer(const label patchi, trackingData &td)
Convert processor patch addressing to the global equivalents.
#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.
vector deviationFromMeshCentre() const
Get the displacement from the mesh centre. Used to correct the.
virtual bool separated() const
Are the planes separated.
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>.
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 appropriate for decomposition or reconstruction.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but stops when a face is hit.
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 stops when a tet triangle is hit. 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.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
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 &)
scalar trackToCell(const vector &displacement, const scalar fraction)
As particle::track, but stops when a new cell is reached.
virtual const vectorField & separation() const
If the planes are separated the separation vector.
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.
virtual const tensorField & forwardT() const
Return face transformation tensor.
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
virtual bool parallel() const
Are the cyclic planes parallel.
dimensionedScalar pos(const dimensionedScalar &ds)
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
An ordered pair of two objects of type <T> with first() and second() elements.
const labelUList & faceCells() const
Return face-cell addressing.
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.
static int compare(const face &, const face &)
Compare faces.
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. The particle may still be.
virtual const labelList & faceOwner() const
Return face owner.
void prepareForParallelTransfer()
Convert global addressing to the processor patch local equivalents.
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.
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.
const dimensionedScalar c
Speed of light in a vacuum.
label start() const
Return start label of this patch in the polyMesh face list.
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 provide unique ID.
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.
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
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 &)
scalar value(const scalar x) const
Evaluate at x.
vector position() const
Return current particle position.
static const Vector< scalar > zero
label patch() const
Return the index of patch that the particle is on.
void type(const direction i, const rootType t)
Set the type of the i-th root.