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_),
621 coordinates_(
p.coordinates_),
623 tetFacei_(
p.tetFacei_),
626 stepFraction_(
p.stepFraction_),
628 nBehind_(
p.nBehind_),
629 origProc_(
p.origProc_),
638 const vector& displacement,
639 const scalar fraction
642 scalar
f = trackToFace(displacement, fraction);
644 while (onInternalFace())
648 f *= trackToFace(
f*displacement,
f*fraction);
657 const vector& displacement,
658 const scalar fraction
663 label tetTriI = onFace() ? 0 : -1;
668 while (nBehind_ < maxNBehind_)
670 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
677 else if (tetTriI == 0)
691 static label stuckID = -1, stuckProc = -1;
692 if (origId_ != stuckID && origProc_ != stuckProc)
695 <<
"Particle #" << origId_ <<
" got stuck at " << position()
700 stuckProc = origProc_;
702 stepFraction_ +=
f*fraction;
713 const vector& displacement,
714 const scalar fraction,
718 const vector x0 = position();
719 const vector x1 = displacement;
724 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
725 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
732 stationaryTetReverseTransform(centre, detA,
T);
736 Pout<<
"Tet points: " << currentTetIndices().tet(mesh_) <<
endl 737 <<
"Tet determinant = " << detA <<
endl 738 <<
"Start local coordinates = " <<
y0 <<
endl;
746 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
751 scalar muH = detA <= 0 ? VGREAT : 1/detA;
752 for (label i = 0; i < 4; ++ i)
754 if (Tx1[i] < - detA*SMALL)
756 scalar
mu = -
y0[i]/Tx1[i];
760 Pout<<
"Hit on tet face " << i <<
" at local coordinate " 761 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the " 762 <<
"way along the track" <<
endl;
765 if (0 <=
mu &&
mu < muH)
777 for (label i = 0; i < 4; ++ i)
796 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
800 Pout<<
"Track hit no tet faces" <<
endl;
802 Pout<<
"End local coordinates = " << yH <<
endl 803 <<
"End global coordinates = " << position() <<
endl 804 <<
"Tracking displacement = " << position() - x0 <<
endl 805 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to " 806 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
810 stepFraction_ += fraction*muH*detA;
814 Pout <<
"Step Fraction : " << stepFraction_ <<
endl;
815 Pout <<
"fraction*muH*detA : " << fraction*muH*detA <<
endl;
816 Pout <<
"static muH : " << muH <<
endl;
817 Pout <<
"origId() : " << origId() <<
endl;
821 if (detA <= 0 || nBehind_ > 0)
823 behind_ += muH*detA*
mag(displacement);
836 return iH != -1 ? 1 - muH*detA : 0;
842 const vector& displacement,
843 const scalar fraction,
847 const vector x0 = position();
848 const vector x1 = displacement;
853 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
854 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
859 FixedList<scalar, 4> detA;
860 FixedList<barycentricTensor, 3>
T;
861 movingTetReverseTransform(fraction, centre, detA,
T);
865 Pair<vector> o,
b, v1, v2;
866 movingTetGeometry(fraction, o,
b, v1, v2);
867 Pout<<
"Tet points o=" << o[0] <<
", b=" <<
b[0]
868 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl 869 <<
"Tet determinant = " << detA[0] <<
endl 870 <<
"Start local coordinates = " <<
y0[0] <<
endl;
875 const vector x0Rel = x0 - centre[0];
876 const vector x1Rel = x1 - centre[1];
879 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
882 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
884 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
886 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
888 FixedList<cubicEqn, 4> hitEqn;
891 hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
896 for (label i = 0; i < 4; ++ i)
898 Pout<< (i ?
" " :
"Hit equation ") << i <<
" = " 899 << hitEqn[i] <<
endl;
901 Pout<<
" DetA equation = " << detA <<
endl;
906 scalar muH = detA[0] <= 0 ? VGREAT : 1/detA[0];
907 for (label i = 0; i < 4; ++i)
909 const Roots<3>
mu = hitEqn[i].roots();
911 for (label j = 0; j < 3; ++j)
916 && hitEqn[i].derivative(
mu[j]) < - detA[0]*SMALL
923 hitEqn[0].value(
mu[j]),
924 hitEqn[1].value(
mu[j]),
925 hitEqn[2].value(
mu[j]),
926 hitEqn[3].value(
mu[j])
928 const scalar detAH = detAEqn.value(
mu[j]);
930 Pout<<
"Hit on tet face " << i <<
" at local coordinate " 931 << (std::isnormal(detAH) ?
name(yH/detAH) :
"???")
932 <<
", " <<
mu[j]*detA[0]*100 <<
"% of the " 933 <<
"way along the track" <<
endl;
935 Pout<<
"derivative : " << hitEqn[i].derivative(
mu[j]) <<
nl 936 <<
" coord " << j <<
" mu[j]: " <<
mu[j] <<
nl 937 <<
" hitEq " << i <<
endl;
940 if (0 <=
mu[j] &&
mu[j] < muH)
952 hitEqn[0].value(muH),
953 hitEqn[1].value(muH),
954 hitEqn[2].value(muH),
964 const scalar detAH = detAEqn.value(muH);
965 if (!std::isnormal(detAH))
968 <<
"A moving tet collapsed onto a particle. This is not supported. " 969 <<
"The mesh is too poor, or the motion too severe, for particle " 975 for (label i = 0; i < 4; ++ i)
977 yH.replace(i, i == iH ? 0 :
max(0, yH[i]));
990 scalar advance = muH*detA[0];
993 stepFraction_ += fraction*advance;
996 if (detA[0] <= 0 || nBehind_ > 0)
998 behind_ += muH*detA[0]*
mag(displacement);
1015 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
1019 Pout<<
"Track hit no tet faces" <<
endl;
1030 return iH != -1 ? 1 - muH*detA[0] : 0;
1036 const vector& displacement,
1037 const scalar fraction,
1041 if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
1043 return trackToMovingTri(displacement, fraction, tetTriI);
1047 return trackToStationaryTri(displacement, fraction, tetTriI);
1054 if (
cmptMin(mesh_.geometricD()) == -1)
1062 return vector::zero;
1078 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
1089 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
1099 transformProperties(
T);
1109 transformProperties(-
s);
1114 facei_ += ppp.
start();
1119 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1160 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1164 tetFacei_ = mesh_.cells()[celli_][0];
1175 if (mesh_.moving() && stepFraction_ != 1)
1177 Pair<vector> centre;
1178 FixedList<scalar, 4> detA;
1179 FixedList<barycentricTensor, 3>
T;
1180 movingTetReverseTransform(0, centre, detA,
T);
1181 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1188 stationaryTetReverseTransform(centre, detA,
T);
1189 coordinates_ += (
pos - centre) &
T/detA;
1196 const polyMesh& procMesh,
1197 const label procCell,
1198 const label procTetFace
1207 (mesh_.faceOwner()[tetFacei_] == celli_)
1208 == (procMesh.faceOwner()[procTetFace] == procCell)
1215 return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
1223 const mapPolyMesh& mapper
1230 mapper.reverseCellMap()[celli_],
1232 "Particle mapped to a location outside of the mesh" 1245 "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.
#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.
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 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.
label start() const
Return start label of this patch in the polyMesh face list.
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.