43 const Foam::label Foam::particle::maxNBehind_ = 10;
56 "writeLagrangianPositions",
64 void Foam::particle::stationaryTetReverseTransform
81 detA = ab & (ac ^ ad);
93 void Foam::particle::movingTetReverseTransform
95 const scalar fraction,
97 FixedList<scalar, 4>& detA,
98 FixedList<barycentricTensor, 3>&
T 101 Pair<barycentricTensor>
A = movingTetTransform(fraction);
103 const Pair<vector> ab(
A[0].
b() -
A[0].a(),
A[1].
b() -
A[1].a());
104 const Pair<vector> ac(
A[0].
c() -
A[0].a(),
A[1].
c() -
A[1].a());
105 const Pair<vector> ad(
A[0].d() -
A[0].a(),
A[1].d() -
A[1].a());
106 const Pair<vector> bc(
A[0].
c() -
A[0].
b(),
A[1].
c() -
A[1].
b());
107 const Pair<vector> bd(
A[0].d() -
A[0].
b(),
A[1].d() -
A[1].
b());
109 centre[0] =
A[0].a();
110 centre[1] =
A[1].a();
112 detA[0] = ab[0] & (ac[0] ^ ad[0]);
114 (ab[1] & (ac[0] ^ ad[0]))
115 + (ab[0] & (ac[1] ^ ad[0]))
116 + (ab[0] & (ac[0] ^ ad[1]));
118 (ab[0] & (ac[1] ^ ad[1]))
119 + (ab[1] & (ac[0] ^ ad[1]))
120 + (ab[1] & (ac[1] ^ ad[0]));
121 detA[3] = ab[1] & (ac[1] ^ ad[1]);
132 (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
133 (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
134 (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
135 (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
147 void Foam::particle::reflect()
149 std::swap(coordinates_.c(), coordinates_.d());
153 void Foam::particle::rotate(
const bool reverse)
157 scalar temp = coordinates_.b();
158 coordinates_.b() = coordinates_.c();
159 coordinates_.c() = coordinates_.d();
160 coordinates_.d() = temp;
164 scalar temp = coordinates_.d();
165 coordinates_.d() = coordinates_.c();
166 coordinates_.c() = coordinates_.b();
167 coordinates_.b() = temp;
172 void Foam::particle::changeTet(
const label tetTriI)
174 const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
176 const label firstTetPtI = 1;
177 const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
183 else if (tetTriI == 2)
187 if (tetPti_ == lastTetPtI)
199 if (tetPti_ == firstTetPtI)
210 else if (tetTriI == 3)
214 if (tetPti_ == firstTetPtI)
226 if (tetPti_ == lastTetPtI)
240 <<
"Changing tet without changing cell should only happen when the " 241 <<
"track is on triangle 1, 2 or 3." 247 void Foam::particle::changeFace(
const label tetTriI)
250 const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
256 sharedEdge = edge(triOldIs[1], triOldIs[2]);
258 else if (tetTriI == 2)
260 sharedEdge = edge(triOldIs[2], triOldIs[0]);
262 else if (tetTriI == 3)
264 sharedEdge = edge(triOldIs[0], triOldIs[1]);
269 <<
"Changing face without changing cell should only happen when the" 270 <<
" track is on triangle 1, 2 or 3." 273 sharedEdge = edge(-1, -1);
279 for (
const label newFaceI : mesh_.cells()[celli_])
281 const Foam::face& newFace = mesh_.faces()[newFaceI];
282 const label newOwner = mesh_.faceOwner()[newFaceI];
285 if (tetFacei_ == newFaceI)
293 const label edgeComp = newOwner == celli_ ? -1 : +1;
298 edgeI < newFace.
size()
304 if (edgeI >= newFace.
size())
310 const label newBaseI =
max(0, mesh_.tetBasePtIs()[newFaceI]);
311 edgeI = (edgeI - newBaseI + newFace.
size()) % newFace.
size();
316 edgeI =
min(
max(1, edgeI), newFace.
size() - 2);
319 tetFacei_ = newFaceI;
329 <<
"The search for an edge-connected face and tet-point failed." 334 if (sharedEdge.otherVertex(triOldIs[1]) == -1)
338 else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
344 const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
350 if (sharedEdge.otherVertex(triNewIs[1]) == -1)
354 else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
361 void Foam::particle::changeCell()
364 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
365 const bool isOwner = celli_ == ownerCellI;
366 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
372 void Foam::particle::changeToMasterPatch()
374 label thisPatch =
patch();
376 for (
const label otherFaceI : mesh_.cells()[celli_])
379 if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
387 const Foam::face& thisFace = mesh_.faces()[facei_];
391 const label otherPatch =
392 mesh_.boundaryMesh().whichPatch(otherFaceI);
393 if (thisPatch > otherPatch)
396 thisPatch = otherPatch;
405 void Foam::particle::locate
410 const bool boundaryFail,
411 const string& boundaryMsg
424 celli_ = mesh_.cellTree().findInside(position);
428 <<
"Cell not found for particle position " << position <<
"." 435 const vector displacement =
437 - mesh_.cellCentres()[celli_]
438 +
vector(VSMALL, VSMALL, VSMALL);
444 scalar minF = VGREAT;
445 label minTetFacei = -1, minTetPti = -1;
449 for (label tetPti = 1; tetPti <
f.
size() - 1; ++tetPti)
452 tetFacei_ =
c[cellTetFacei];
457 const scalar
f = trackToTri(displacement, 0, tetTriI);
467 minTetFacei = tetFacei_;
476 tetFacei_ = minTetFacei;
479 track(displacement, 0);
492 static label nWarnings = 0;
493 static const label maxNWarnings = 100;
494 if ((nWarnings < maxNWarnings) && boundaryFail)
499 if (nWarnings == maxNWarnings)
502 <<
"Suppressing any further warnings about particles being " 503 <<
"located outside of the mesh." <<
endl;
518 const label tetFacei,
531 origProc_(
Pstream::myProcNo()),
532 origId_(getNewParticleID())
544 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
552 origProc_(
Pstream::myProcNo()),
553 origId_(getNewParticleID())
561 "Particle initialised with a location outside of the mesh" 571 const label tetFacei,
577 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
585 origProc_(
Pstream::myProcNo()),
586 origId_(getNewParticleID())
596 "Particle initialised with a location outside of the mesh" 605 coordinates_(
p.coordinates_),
607 tetFacei_(
p.tetFacei_),
610 stepFraction_(
p.stepFraction_),
612 nBehind_(
p.nBehind_),
613 origProc_(
p.origProc_),
628 const vector& displacement,
629 const scalar fraction
632 scalar
f = trackToFace(displacement, fraction);
634 while (onInternalFace())
638 f *= trackToFace(
f*displacement,
f*fraction);
647 const vector& displacement,
648 const scalar fraction
653 label tetTriI = onFace() ? 0 : -1;
658 while (nBehind_ < maxNBehind_)
660 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
667 else if (tetTriI == 0)
681 static label stuckID = -1, stuckProc = -1;
682 if (origId_ != stuckID && origProc_ != stuckProc)
685 <<
"Particle #" << origId_ <<
" got stuck at " << position()
690 stuckProc = origProc_;
692 stepFraction_ +=
f*fraction;
703 const vector& displacement,
704 const scalar fraction,
708 const vector x0 = position();
709 const vector x1 = displacement;
714 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
715 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
722 stationaryTetReverseTransform(centre, detA,
T);
726 Pout<<
"Tet points: " << currentTetIndices().tet(mesh_) <<
endl 727 <<
"Tet determinant = " << detA <<
endl 728 <<
"Start local coordinates = " <<
y0 <<
endl;
736 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
741 scalar muH = detA <= 0 ? VGREAT : 1/detA;
742 for (label i = 0; i < 4; ++ i)
744 if (Tx1[i] < - detA*SMALL)
746 scalar
mu = -
y0[i]/Tx1[i];
750 Pout<<
"Hit on tet face " << i <<
" at local coordinate " 751 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the " 752 <<
"way along the track" <<
endl;
755 if (0 <=
mu &&
mu < muH)
767 for (label i = 0; i < 4; ++ i)
786 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
790 Pout<<
"Track hit no tet faces" <<
endl;
792 Pout<<
"End local coordinates = " << yH <<
endl 793 <<
"End global coordinates = " << position() <<
endl 794 <<
"Tracking displacement = " << position() - x0 <<
endl 795 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to " 796 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
800 stepFraction_ += fraction*muH*detA;
804 Pout <<
"Step Fraction : " << stepFraction_ <<
endl;
805 Pout <<
"fraction*muH*detA : " << fraction*muH*detA <<
endl;
806 Pout <<
"static muH : " << muH <<
endl;
807 Pout <<
"origId() : " << origId() <<
endl;
811 if (detA <= 0 || nBehind_ > 0)
813 behind_ += muH*detA*
mag(displacement);
826 return iH != -1 ? 1 - muH*detA : 0;
832 const vector& displacement,
833 const scalar fraction,
837 const vector x0 = position();
838 const vector x1 = displacement;
843 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
844 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
849 FixedList<scalar, 4> detA;
850 FixedList<barycentricTensor, 3>
T;
851 movingTetReverseTransform(fraction, centre, detA,
T);
855 Pair<vector> o,
b, v1, v2;
856 movingTetGeometry(fraction, o,
b, v1, v2);
857 Pout<<
"Tet points o=" << o[0] <<
", b=" <<
b[0]
858 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl 859 <<
"Tet determinant = " << detA[0] <<
endl 860 <<
"Start local coordinates = " <<
y0[0] <<
endl;
865 const vector x0Rel = x0 - centre[0];
866 const vector x1Rel = x1 - centre[1];
869 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
872 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
874 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
876 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
878 FixedList<cubicEqn, 4> hitEqn;
881 hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
886 for (label i = 0; i < 4; ++ i)
888 Pout<< (i ?
" " :
"Hit equation ") << i <<
" = " 889 << hitEqn[i] <<
endl;
891 Pout<<
" DetA equation = " << detA <<
endl;
896 scalar muH = detA[0] <= 0 ? VGREAT : 1/detA[0];
897 for (label i = 0; i < 4; ++i)
899 const Roots<3>
mu = hitEqn[i].roots();
901 for (label j = 0; j < 3; ++j)
906 && hitEqn[i].derivative(
mu[j]) < - detA[0]*SMALL
913 hitEqn[0].value(
mu[j]),
914 hitEqn[1].value(
mu[j]),
915 hitEqn[2].value(
mu[j]),
916 hitEqn[3].value(
mu[j])
918 const scalar detAH = detAEqn.value(
mu[j]);
920 Pout<<
"Hit on tet face " << i <<
" at local coordinate " 921 << (std::isnormal(detAH) ?
name(yH/detAH) :
"???")
922 <<
", " <<
mu[j]*detA[0]*100 <<
"% of the " 923 <<
"way along the track" <<
endl;
925 Pout<<
"derivative : " << hitEqn[i].derivative(
mu[j]) <<
nl 926 <<
" coord " << j <<
" mu[j]: " <<
mu[j] <<
nl 927 <<
" hitEq " << i <<
endl;
930 if (0 <=
mu[j] &&
mu[j] < muH)
942 hitEqn[0].value(muH),
943 hitEqn[1].value(muH),
944 hitEqn[2].value(muH),
954 const scalar detAH = detAEqn.value(muH);
955 if (!std::isnormal(detAH))
958 <<
"A moving tet collapsed onto a particle. This is not supported. " 959 <<
"The mesh is too poor, or the motion too severe, for particle " 965 for (label i = 0; i < 4; ++ i)
967 yH.replace(i, i == iH ? 0 :
max(0, yH[i]));
980 scalar advance = muH*detA[0];
983 stepFraction_ += fraction*advance;
986 if (detA[0] <= 0 || nBehind_ > 0)
988 behind_ += muH*detA[0]*
mag(displacement);
1005 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
1009 Pout<<
"Track hit no tet faces" <<
endl;
1020 return iH != -1 ? 1 - muH*detA[0] : 0;
1026 const vector& displacement,
1027 const scalar fraction,
1031 if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
1033 return trackToMovingTri(displacement, fraction, tetTriI);
1037 return trackToStationaryTri(displacement, fraction, tetTriI);
1044 if (
cmptMin(mesh_.geometricD()) == -1)
1052 return vector::zero;
1068 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
1079 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
1089 transformProperties(
T);
1099 transformProperties(-
s);
1104 facei_ += ppp.
start();
1109 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1150 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1154 tetFacei_ = mesh_.cells()[celli_][0];
1165 if (mesh_.moving() && stepFraction_ != 1)
1167 Pair<vector> centre;
1168 FixedList<scalar, 4> detA;
1169 FixedList<barycentricTensor, 3>
T;
1170 movingTetReverseTransform(0, centre, detA,
T);
1171 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1178 stationaryTetReverseTransform(centre, detA,
T);
1179 coordinates_ += (
pos - centre) &
T/detA;
1186 const polyMesh& procMesh,
1187 const label procCell,
1188 const label procTetFace
1197 (mesh_.faceOwner()[tetFacei_] == celli_)
1198 == (procMesh.faceOwner()[procTetFace] == procCell)
1205 return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
1213 const mapPolyMesh& mapper
1220 mapper.reverseCellMap()[celli_],
1222 "Particle mapped to a location outside of the mesh" 1235 "Particle mapped to a location outside of the mesh"
void correctAfterParallelTransfer(const label patchi, trackingData &td)
Convert processor patch addressing to the global equivalents.
void size(const label n)
Older name for setAddressableSize.
vector deviationFromMeshCentre() const
Get the displacement from the mesh centre. Used to correct the.
virtual bool separated() const
Are the planes separated.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
label start() const noexcept
Return start label of this patch in the polyMesh face list.
#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.
label origId() const noexcept
Return the particle ID on the originating processor.
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but also stops on internal faces.
scalar trackToTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToFace, but also stops on tet triangles. On.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool writeLagrangianCoordinates
Write particle coordinates file (v1712 and later) Default is true.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
dimensionedScalar y0(const dimensionedScalar &ds)
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
int infoSwitch(const char *name, const int deflt=0)
Lookup info switch or add default value.
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.
virtual const tensorField & forwardT() const
Return face transformation tensor.
#define forAll(list, i)
Loop across all elements in list.
virtual bool parallel() const
Are the cyclic planes parallel.
dimensionedScalar pos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
const labelUList & faceCells() const
Return face-cell addressing.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
Inter-processor communications stream.
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. The particle may still be.
void prepareForParallelTransfer()
Convert global addressing to the processor patch local equivalents.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void autoMap(const vector &position, const mapPolyMesh &mapper)
Map after a topology change.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Foam::edge edge(const label edgei) const
Return i-th face edge (forward walk order).
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
const dimensionedScalar mu
Atomic mass unit.
label origProc() const noexcept
Return the originating processor ID.
scalar trackToMovingTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for moving meshes.
PtrList< coordinateSystem > coordinates(solidRegions.size())
#define WarningInFunction
Report a warning using Foam::Warning.
static bool writeLagrangianPositions
Write particle positions file (v1706 format and earlier) Default is true (disable in etc/controlDict)...
A cell is defined as a list of faces with extra functionality.
const dimensionedScalar c
Speed of light in a vacuum.
const std::string patch
OpenFOAM patch number as a std::string.
void relocate(const point &position, const label celli=-1)
Set the addressing based on the provided position.
registerInfoSwitch("writeLagrangianPositions", bool, Foam::particle::writeLagrangianPositions)
static label particleCount_
Cumulative particle counter - used to provide unique ID.
bool operator!=(const eddy &a, const eddy &b)
Mesh consisting of general polyhedral cells.
static int compare(const face &a, const face &b)
Compare faces.
void replace(const direction d, const dimensioned< cmptType > &dc)
Return a component with a dimensioned<cmptType>
scalar trackToStationaryTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for stationary meshes.
tmp< faMatrix< Type > > operator==(const faMatrix< Type > &, const faMatrix< Type > &)
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
scalar track(const vector &displacement, const scalar fraction)
Track along the displacement for a given fraction of the overall.
Tensor of scalars, i.e. Tensor<scalar>.
vector position() const
Return current particle position.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static int compare(const edge &a, const edge &b)
Compare edges.