Go to the documentation of this file.
42 const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
55 "writeLagrangianPositions",
63 void Foam::particle::stationaryTetReverseTransform
80 detA = ab & (ac ^ ad);
92 void Foam::particle::movingTetReverseTransform
94 const scalar fraction,
96 FixedList<scalar, 4>& detA,
97 FixedList<barycentricTensor, 3>&
T
100 Pair<barycentricTensor>
A = movingTetTransform(fraction);
102 const Pair<vector> ab(
A[0].
b() -
A[0].a(),
A[1].
b() -
A[1].a());
103 const Pair<vector> ac(
A[0].
c() -
A[0].a(),
A[1].
c() -
A[1].a());
104 const Pair<vector> ad(
A[0].d() -
A[0].a(),
A[1].d() -
A[1].a());
105 const Pair<vector> bc(
A[0].
c() -
A[0].
b(),
A[1].
c() -
A[1].
b());
106 const Pair<vector> bd(
A[0].d() -
A[0].
b(),
A[1].d() -
A[1].
b());
108 centre[0] =
A[0].a();
109 centre[1] =
A[1].a();
111 detA[0] = ab[0] & (ac[0] ^ ad[0]);
113 (ab[1] & (ac[0] ^ ad[0]))
114 + (ab[0] & (ac[1] ^ ad[0]))
115 + (ab[0] & (ac[0] ^ ad[1]));
117 (ab[0] & (ac[1] ^ ad[1]))
118 + (ab[1] & (ac[0] ^ ad[1]))
119 + (ab[1] & (ac[1] ^ ad[0]));
120 detA[3] = ab[1] & (ac[1] ^ ad[1]);
131 (bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
132 (ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
133 (ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
134 (ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
146 void Foam::particle::reflect()
148 Swap(coordinates_.c(), coordinates_.d());
152 void Foam::particle::rotate(
const bool reverse)
156 scalar temp = coordinates_.b();
157 coordinates_.b() = coordinates_.c();
158 coordinates_.c() = coordinates_.d();
159 coordinates_.d() = temp;
163 scalar temp = coordinates_.d();
164 coordinates_.d() = coordinates_.c();
165 coordinates_.c() = coordinates_.b();
166 coordinates_.b() = temp;
171 void Foam::particle::changeTet(
const label tetTriI)
173 const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
175 const label firstTetPtI = 1;
176 const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
182 else if (tetTriI == 2)
186 if (tetPti_ == lastTetPtI)
198 if (tetPti_ == firstTetPtI)
209 else if (tetTriI == 3)
213 if (tetPti_ == firstTetPtI)
225 if (tetPti_ == lastTetPtI)
239 <<
"Changing tet without changing cell should only happen when the "
240 <<
"track is on triangle 1, 2 or 3."
246 void Foam::particle::changeFace(
const label tetTriI)
249 const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
255 sharedEdge = edge(triOldIs[1], triOldIs[2]);
257 else if (tetTriI == 2)
259 sharedEdge = edge(triOldIs[2], triOldIs[0]);
261 else if (tetTriI == 3)
263 sharedEdge = edge(triOldIs[0], triOldIs[1]);
268 <<
"Changing face without changing cell should only happen when the"
269 <<
" track is on triangle 1, 2 or 3."
272 sharedEdge = edge(-1, -1);
278 forAll(mesh_.cells()[celli_], cellFaceI)
280 const label newFaceI = mesh_.cells()[celli_][cellFaceI];
281 const class 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()
299 &&
edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
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;
373 void Foam::particle::changeToMasterPatch()
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;
407 void Foam::particle::locate
412 const bool boundaryFail,
413 const string& boundaryMsg
421 celli_ = mesh_.cellTree().findInside(position);
426 <<
"Cell not found for particle position " << position <<
"."
434 coordinates_ =
barycentric(1-3*SMALL, SMALL, SMALL, SMALL);
435 tetFacei_ = mesh_.cells()[celli_][0];
440 track(position - this->position(), 0);
451 <<
" when tracking from centre " << mesh_.cellCentres()[celli_]
452 <<
" of cell " << celli_ <<
" to position " << position
463 const polyPatch&
p = mesh_.boundaryMesh()[
patch()];
468 const vector s = position - mesh_.cellCentres()[celli_];
473 tetFacei_ = mesh_.cells()[celli_][0];
480 static label nWarnings = 0;
481 static const label maxNWarnings = 100;
482 if (nWarnings < maxNWarnings)
485 <<
" when tracking from centre " << mesh_.cellCentres()[celli_]
486 <<
" of cell " << celli_ <<
" to position " << position
490 if (nWarnings == maxNWarnings)
493 <<
"Suppressing any further warnings about particles being "
494 <<
"located outside of the mesh." <<
endl;
508 const label tetFacei,
520 origId_(getNewParticleID())
532 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
539 origId_(getNewParticleID())
547 "Particle initialised with a location outside of the mesh"
557 const label tetFacei,
563 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
570 origId_(getNewParticleID())
580 "Particle initialised with a location outside of the mesh"
589 coordinates_(
p.coordinates_),
591 tetFacei_(
p.tetFacei_),
594 stepFraction_(
p.stepFraction_),
595 origProc_(
p.origProc_),
603 coordinates_(
p.coordinates_),
605 tetFacei_(
p.tetFacei_),
608 stepFraction_(
p.stepFraction_),
609 origProc_(
p.origProc_),
618 const vector& displacement,
619 const scalar fraction
622 scalar
f = trackToFace(displacement, fraction);
624 while (onInternalFace())
628 f *= trackToFace(
f*displacement,
f*fraction);
637 const vector& displacement,
638 const scalar fraction
643 label tetTriI = onFace() ? 0 : -1;
649 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
656 else if (tetTriI == 0)
673 const vector& displacement,
674 const scalar fraction,
678 const vector x0 = position();
679 const vector x1 = displacement;
684 Info<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
685 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
692 stationaryTetReverseTransform(centre, detA,
T);
697 stationaryTetGeometry(o,
b, v1, v2);
698 Info<<
"Tet points o=" << o <<
", b=" <<
b
699 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl
700 <<
"Tet determinant = " << detA <<
endl
701 <<
"Start local coordinates = " <<
y0 <<
endl;
705 const scalar
f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
712 Info<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
717 scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
718 for (
label i = 0; i < 4; ++ i)
720 if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
722 scalar
mu = -
y0[i]/Tx1[i];
726 Info<<
"Hit on tet face " << i <<
" at local coordinate "
727 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the "
728 <<
"way along the track" <<
endl;
731 if (0 <=
mu &&
mu < muH)
743 for (
label i = 0; i < 4; ++ i)
762 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
766 Info<<
"Track hit no tet faces" <<
endl;
768 Info<<
"End local coordinates = " << yH <<
endl
769 <<
"End global coordinates = " << position() <<
endl
770 <<
"Tracking displacement = " << position() - x0 <<
endl
771 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to "
772 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
776 stepFraction_ += fraction*muH*detA;
778 return iH != -1 ? 1 - muH*detA : 0;
784 const vector& displacement,
785 const scalar fraction,
789 const vector x0 = position();
790 const vector x1 = displacement;
795 Info<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
796 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
803 movingTetReverseTransform(fraction, centre, detA,
T);
806 const scalar
f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
809 const vector x0Rel = x0 - centre[0];
810 const vector x1Rel =
f*x1 - centre[1];
813 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
816 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
818 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
820 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
825 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
830 scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
831 for (
label i = 0; i < 4; ++i)
835 for (
label j = 0; j < 3; ++j)
839 if (0 <=
mu[j] &&
mu[j] < muH)
851 hitEqn[0].value(muH),
852 hitEqn[1].value(muH),
853 hitEqn[2].value(muH),
863 const scalar detAH = detAEqn.
value(muH);
864 if (!std::isnormal(detAH))
867 <<
"A moving tet collapsed onto a particle. This is not supported. "
868 <<
"The mesh is too poor, or the motion too severe, for particle "
874 for (
label i = 0; i < 4; ++ i)
890 stepFraction_ += fraction*muH*detA[0];
896 Info<<
"Track hit tet face " << iH <<
" first" <<
endl;
900 Info<<
"Track hit no tet faces" <<
endl;
902 Info<<
"End local coordinates = " << yH <<
endl
903 <<
"End global coordinates = " << position() <<
endl;
906 return iH != -1 ? 1 - muH*detA[0] : 0;
912 const vector& displacement,
913 const scalar fraction,
919 return trackToMovingTri(displacement, fraction, tetTriI);
923 return trackToStationaryTri(displacement, fraction, tetTriI);
930 if (
cmptMin(mesh_.geometricD()) == -1)
954 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
965 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
975 transformProperties(
T);
985 transformProperties(-
s);
990 facei_ += ppp.
start();
995 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1036 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1040 tetFacei_ = mesh_.cells()[celli_][0];
1056 movingTetReverseTransform(0, centre, detA,
T);
1057 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1064 stationaryTetReverseTransform(centre, detA,
T);
1065 coordinates_ += (
pos - centre) &
T/detA;
1073 const label procCell,
1074 const label procTetFace
1083 (mesh_.faceOwner()[tetFacei_] == celli_)
1084 == (procMesh.
faceOwner()[procTetFace] == procCell)
1091 return procMesh.
faces()[procTetFace].size() - 1 - tetPti_;
1108 "Particle mapped to a location outside of the mesh"
1121 "Particle mapped to a location outside of the mesh"
int debug
Static debugging option.
void prepareForParallelTransfer()
Convert global addressing to the processor patch local equivalents.
void reverse(UList< T > &list, const label n)
virtual const vectorField & separation() const
If the planes are separated the separation vector.
void correctAfterInteractionListReferral(const label celli)
Correct the topology after referral. The particle may still be.
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
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))
const dimensionedScalar mu
Atomic mass unit.
label origId() const
Return the particle ID on the originating processor.
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.
virtual const tensorField & forwardT() const
Return face transformation tensor.
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
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...
void relocate(const point &position, const label celli=-1)
Set the addressing based on the provided position.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
void Swap(DynamicList< T, SizeMin1 > &a, DynamicList< T, SizeMin2 > &b)
static label particleCount_
Cumulative particle counter - used to provide unique ID.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Mesh consisting of general polyhedral cells.
static bool writeLagrangianPositions
#define forAll(list, i)
Loop across all elements in list.
static int compare(const face &a, const face &b)
Compare faces.
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
virtual bool parallel() const
Are the cyclic planes parallel.
int infoSwitch(const char *name, const int deflt=0)
Lookup info switch or add default value.
scalar track(const vector &displacement, const scalar fraction)
Track along the displacement for a given fraction of the overall.
bool operator!=(const eddy &a, const eddy &b)
virtual bool separated() const
Are the planes separated.
scalar trackToStationaryTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToTri, but for stationary meshes.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
tmp< faMatrix< Type > > operator==(const faMatrix< Type > &, const faMatrix< Type > &)
messageStream Info
Information stream (uses stdout - output is on the master only)
scalar value(const scalar x) const
Evaluate at x.
PtrList< coordinateSystem > coordinates(solidRegions.size())
virtual const labelList & faceOwner() const
Return face owner.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensionedScalar y0(const dimensionedScalar &ds)
vector deviationFromMeshCentre() const
Get the displacement from the mesh centre. Used to correct the.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
void replace(const direction, const Cmpt &)
label origProc() const
Return the originating processor ID.
void correctAfterParallelTransfer(const label patchi, trackingData &td)
Convert processor patch addressing to the global equivalents.
label procTetPt(const polyMesh &procMesh, const label procCell, const label procTetFace) const
Return the tet point appropriate for decomposition or reconstruction.
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
const labelUList & faceCells() const
Return face-cell addressing.
Vector< scalar > vector
A scalar version of the templated Vector.
label start() const
Return start label of this patch in the polyMesh face list.
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but also stops on internal faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
static bool writeLagrangianCoordinates
const labelList & reverseCellMap() const
Reverse cell map.
virtual const faceList & faces() const
Return raw faces.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static int compare(const edge &a, const edge &b)
Compare edges.
scalar trackToTri(const vector &displacement, const scalar fraction, label &tetTriI)
As particle::trackToFace, but also stops on tet triangles. On.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
An ordered pair of two objects of type <T> with first() and second() elements.
const std::string patch
OpenFOAM patch number as a std::string.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A 1D vector of objects of type <T> with a fixed length <N>.
registerInfoSwitch("writeLagrangianPositions", bool, Foam::particle::writeLagrangianPositions)
void prepareForInteractionListReferral(const vectorTensorTransform &transform)
Break the topology and store the particle position so that the.
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
const dimensionedScalar c
Speed of light in a vacuum.
static const Vector< scalar > zero
Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0.
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
defineTypeNameAndDebug(combustionModel, 0)
#define WarningInFunction
Report a warning using Foam::Warning.
particle(const polyMesh &mesh, const barycentric &coordinates, const label celli, const label tetFacei, const label tetPti)
Construct from components.
dimensionedScalar pos(const dimensionedScalar &ds)