43const Foam::label Foam::particle::maxNBehind_ = 10;
56 "writeLagrangianPositions",
64void Foam::particle::stationaryTetReverseTransform
81 detA = ab & (ac ^ ad);
93void 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])
147void Foam::particle::reflect()
149 std::swap(coordinates_.c(), coordinates_.d());
153void 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;
172void 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."
247void 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 forAll(mesh_.cells()[celli_], cellFaceI)
281 const label newFaceI = mesh_.cells()[celli_][cellFaceI];
282 const class face& newFace = mesh_.faces()[newFaceI];
283 const label newOwner = mesh_.faceOwner()[newFaceI];
286 if (tetFacei_ == newFaceI)
294 const label edgeComp = newOwner == celli_ ? -1 : +1;
299 edgeI < newFace.size()
300 &&
edge::compare(sharedEdge, newFace.edge(edgeI)) != edgeComp;
305 if (edgeI >= newFace.size())
311 const label newBaseI =
max(0, mesh_.tetBasePtIs()[newFaceI]);
312 edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
317 edgeI =
min(
max(1, edgeI), newFace.size() - 2);
320 tetFacei_ = newFaceI;
330 <<
"The search for an edge-connected face and tet-point failed."
335 if (sharedEdge.otherVertex(triOldIs[1]) == -1)
339 else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
345 const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
351 if (sharedEdge.otherVertex(triNewIs[1]) == -1)
355 else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
362void Foam::particle::changeCell()
365 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
366 const bool isOwner = celli_ == ownerCellI;
367 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
373void Foam::particle::changeToMasterPatch()
375 label thisPatch =
patch();
377 forAll(mesh_.cells()[celli_], cellFaceI)
380 const label otherFaceI = mesh_.cells()[celli_][cellFaceI];
381 if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
389 const class face& thisFace = mesh_.faces()[facei_];
390 const class face&
otherFace = mesh_.faces()[otherFaceI];
393 const label otherPatch =
394 mesh_.boundaryMesh().whichPatch(otherFaceI);
395 if (thisPatch > otherPatch)
398 thisPatch = otherPatch;
407void Foam::particle::locate
412 const bool boundaryFail,
413 const string& boundaryMsg
426 celli_ = mesh_.cellTree().findInside(position);
430 <<
"Cell not found for particle position " << position <<
"."
437 const vector displacement =
439 - mesh_.cellCentres()[celli_]
440 +
vector(VSMALL, VSMALL, VSMALL);
445 const class cell&
c = mesh_.cells()[celli_];
446 scalar minF = VGREAT;
447 label minTetFacei = -1, minTetPti = -1;
450 const class face&
f = mesh_.faces()[
c[cellTetFacei]];
451 for (label tetPti = 1; tetPti <
f.
size() - 1; ++tetPti)
454 tetFacei_ =
c[cellTetFacei];
459 const scalar
f = trackToTri(displacement, 0, tetTriI);
469 minTetFacei = tetFacei_;
478 tetFacei_ = minTetFacei;
481 track(displacement, 0);
494 static label nWarnings = 0;
495 static const label maxNWarnings = 100;
496 if ((nWarnings < maxNWarnings) && boundaryFail)
501 if (nWarnings == maxNWarnings)
504 <<
"Suppressing any further warnings about particles being "
505 <<
"located outside of the mesh." <<
endl;
520 const label tetFacei,
533 origProc_(
Pstream::myProcNo()),
534 origId_(getNewParticleID())
546 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
554 origProc_(
Pstream::myProcNo()),
555 origId_(getNewParticleID())
563 "Particle initialised with a location outside of the mesh"
573 const label tetFacei,
579 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
587 origProc_(
Pstream::myProcNo()),
588 origId_(getNewParticleID())
598 "Particle initialised with a location outside of the mesh"
607 coordinates_(
p.coordinates_),
609 tetFacei_(
p.tetFacei_),
612 stepFraction_(
p.stepFraction_),
614 nBehind_(
p.nBehind_),
615 origProc_(
p.origProc_),
623 coordinates_(
p.coordinates_),
625 tetFacei_(
p.tetFacei_),
628 stepFraction_(
p.stepFraction_),
630 nBehind_(
p.nBehind_),
631 origProc_(
p.origProc_),
640 const vector& displacement,
641 const scalar fraction
644 scalar
f = trackToFace(displacement, fraction);
646 while (onInternalFace())
650 f *= trackToFace(
f*displacement,
f*fraction);
659 const vector& displacement,
660 const scalar fraction
665 label tetTriI = onFace() ? 0 : -1;
670 while (nBehind_ < maxNBehind_)
672 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
679 else if (tetTriI == 0)
693 static label stuckID = -1, stuckProc = -1;
694 if (origId_ != stuckID && origProc_ != stuckProc)
697 <<
"Particle #" << origId_ <<
" got stuck at " << position()
702 stuckProc = origProc_;
704 stepFraction_ +=
f*fraction;
715 const vector& displacement,
716 const scalar fraction,
720 const vector x0 = position();
721 const vector x1 = displacement;
726 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
727 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
734 stationaryTetReverseTransform(centre, detA,
T);
739 stationaryTetGeometry(o,
b, v1, v2);
740 Pout<<
"Tet points o=" << o <<
", b=" <<
b
741 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl
742 <<
"Tet determinant = " << detA <<
endl
743 <<
"Start local coordinates = " <<
y0 <<
endl;
751 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
756 scalar muH = detA <= 0 ? VGREAT : 1/detA;
757 for (label i = 0; i < 4; ++ i)
759 if (Tx1[i] < - detA*SMALL)
761 scalar
mu = -
y0[i]/Tx1[i];
765 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
766 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the "
767 <<
"way along the track" <<
endl;
770 if (0 <=
mu &&
mu < muH)
782 for (label i = 0; i < 4; ++ i)
801 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
805 Pout<<
"Track hit no tet faces" <<
endl;
807 Pout<<
"End local coordinates = " << yH <<
endl
808 <<
"End global coordinates = " << position() <<
endl
809 <<
"Tracking displacement = " << position() - x0 <<
endl
810 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to "
811 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
815 stepFraction_ += fraction*muH*detA;
819 Pout <<
"Step Fraction : " << stepFraction_ <<
endl;
820 Pout <<
"fraction*muH*detA : " << fraction*muH*detA <<
endl;
821 Pout <<
"static muH : " << muH <<
endl;
822 Pout <<
"origId() : " << origId() <<
endl;
826 if (detA <= 0 || nBehind_ > 0)
828 behind_ += muH*detA*
mag(displacement);
841 return iH != -1 ? 1 - muH*detA : 0;
847 const vector& displacement,
848 const scalar fraction,
852 const vector x0 = position();
853 const vector x1 = displacement;
858 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
859 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
866 movingTetReverseTransform(fraction, centre, detA,
T);
871 movingTetGeometry(fraction, o,
b, v1, v2);
872 Pout<<
"Tet points o=" << o[0] <<
", b=" <<
b[0]
873 <<
", v1=" << v1[0] <<
", v2=" << v2[0] <<
endl
874 <<
"Tet determinant = " << detA[0] <<
endl
875 <<
"Start local coordinates = " <<
y0[0] <<
endl;
880 const vector x0Rel = x0 - centre[0];
881 const vector x1Rel = x1 - centre[1];
884 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
887 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
889 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
891 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
896 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
901 for (label i = 0; i < 4; ++ i)
903 Pout<< (i ?
" " :
"Hit equation ") << i <<
" = "
904 << hitEqn[i] <<
endl;
906 Pout<<
" DetA equation = " << detA <<
endl;
911 scalar muH = detA[0] <= 0 ? VGREAT : 1/detA[0];
912 for (label i = 0; i < 4; ++i)
916 for (label j = 0; j < 3; ++j)
921 && hitEqn[i].derivative(
mu[j]) < - detA[0]*SMALL
928 hitEqn[0].value(
mu[j]),
929 hitEqn[1].value(
mu[j]),
930 hitEqn[2].value(
mu[j]),
931 hitEqn[3].value(
mu[j])
933 const scalar detAH = detAEqn.
value(
mu[j]);
935 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
936 << (std::isnormal(detAH) ?
name(yH/detAH) :
"???")
937 <<
", " <<
mu[j]*detA[0]*100 <<
"% of the "
938 <<
"way along the track" <<
endl;
940 Pout<<
"derivative : " << hitEqn[i].derivative(
mu[j]) <<
nl
941 <<
" coord " << j <<
" mu[j]: " <<
mu[j] <<
nl
942 <<
" hitEq " << i <<
endl;
945 if (0 <=
mu[j] &&
mu[j] < muH)
957 hitEqn[0].value(muH),
958 hitEqn[1].value(muH),
959 hitEqn[2].value(muH),
969 const scalar detAH = detAEqn.
value(muH);
970 if (!std::isnormal(detAH))
973 <<
"A moving tet collapsed onto a particle. This is not supported. "
974 <<
"The mesh is too poor, or the motion too severe, for particle "
980 for (label i = 0; i < 4; ++ i)
995 scalar advance = muH*detA[0];
998 stepFraction_ += fraction*advance;
1001 if (detA[0] <= 0 || nBehind_ > 0)
1003 behind_ += muH*detA[0]*
mag(displacement);
1020 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
1024 Pout<<
"Track hit no tet faces" <<
endl;
1035 return iH != -1 ? 1 - muH*detA[0] : 0;
1041 const vector& displacement,
1042 const scalar fraction,
1046 if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
1048 return trackToMovingTri(displacement, fraction, tetTriI);
1052 return trackToStationaryTri(displacement, fraction, tetTriI);
1059 if (
cmptMin(mesh_.geometricD()) == -1)
1083 facei_ = mesh_.boundaryMesh()[patch()].whichFace(facei_);
1094 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
1104 transformProperties(
T);
1114 transformProperties(-
s);
1119 facei_ += ppp.
start();
1124 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1165 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1169 tetFacei_ = mesh_.cells()[celli_][0];
1180 if (mesh_.moving() && stepFraction_ != 1)
1185 movingTetReverseTransform(0, centre, detA,
T);
1186 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1193 stationaryTetReverseTransform(centre, detA,
T);
1194 coordinates_ += (
pos - centre) &
T/detA;
1202 const label procCell,
1203 const label procTetFace
1212 (mesh_.faceOwner()[tetFacei_] == celli_)
1213 == (procMesh.
faceOwner()[procTetFace] == procCell)
1220 return procMesh.
faces()[procTetFace].
size() - 1 - tetPti_;
1237 "Particle mapped to a location outside of the mesh"
1250 "Particle mapped to a location outside of the mesh"
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
A 1D vector of objects of type <T> with a fixed length <N>.
An ordered pair of two objects of type <T> with first() and second() elements.
Inter-processor communications stream.
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
void size(const label n)
Older name for setAddressableSize.
void replace(const direction, const Cmpt &)
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
virtual bool separated() const
Are the planes separated.
virtual bool parallel() const
Are the cyclic planes parallel.
virtual const vectorField & separation() const
If the planes are separated the separation vector.
virtual const tensorField & forwardT() const
Return face transformation tensor.
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
scalar value(const scalar x) const
Evaluate the cubic equation at x.
static int compare(const edge &a, const edge &b)
Compare edges.
static int compare(const face &a, const face &b)
Compare faces.
virtual void track()
Do the actual tracking to fill the track data.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const labelList & reverseCellMap() const
Reverse cell map.
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. The particle may still be.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
label procTetPt(const polyMesh &procMesh, const label procCell, const label procTetFace) const
Return the tet point appropriate for decomposition or reconstruction.
vector deviationFromMeshCentre() const
Get the displacement from the mesh centre. Used to correct the.
scalar trackToStationaryTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for stationary meshes.
void correctAfterParallelTransfer(const label patchi, trackingData &td)
Convert processor patch addressing to the global equivalents.
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but also stops on internal faces.
vector position() const
Return current particle position.
void prepareForInteractionListReferral(const vectorTensorTransform &transform)
Break the topology and store the particle position so that the.
void autoMap(const vector &position, const mapPolyMesh &mapper)
Map after a topology change.
scalar trackToMovingTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for moving meshes.
void prepareForParallelTransfer()
Convert global addressing to the processor patch local equivalents.
static label particleCount_
Cumulative particle counter - used to provide unique ID.
static bool writeLagrangianPositions
label origId() const noexcept
Return the particle ID on the originating processor.
void relocate(const point &position, const label celli=-1)
Set the addressing based on the provided position.
scalar trackToTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToFace, but also stops on tet triangles. On.
static bool writeLagrangianCoordinates
label origProc() const noexcept
Return the originating processor ID.
Mesh consisting of general polyhedral cells.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
label start() const
Return start label of this patch in the polyMesh face list.
const labelUList & faceCells() const
Return face-cell addressing.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & mu
PtrList< coordinateSystem > coordinates(solidRegions.size())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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))
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
int infoSwitch(const char *name, const int deflt=0)
Lookup info switch or add default value.
const std::string patch
OpenFOAM patch number as a std::string.
bool operator!=(const eddy &a, const eddy &b)
dimensionedScalar pos(const dimensionedScalar &ds)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar y0(const dimensionedScalar &ds)
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< faMatrix< Type > > operator==(const faMatrix< Type > &, const faMatrix< Type > &)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
#define registerInfoSwitch(Name, Type, SwitchVar)
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.