Go to the documentation of this file.
43 const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
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 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 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.faceEdge(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)
362 void Foam::particle::changeCell()
365 const label ownerCellI = mesh_.faceOwner()[tetFacei_];
366 const bool isOwner = celli_ == ownerCellI;
367 celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
374 void Foam::particle::changeToMasterPatch()
376 label thisPatch =
patch();
378 forAll(mesh_.cells()[celli_], cellFaceI)
381 const label otherFaceI = mesh_.cells()[celli_][cellFaceI];
382 if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
390 const class face& thisFace = mesh_.faces()[facei_];
391 const class face&
otherFace = mesh_.faces()[otherFaceI];
394 const label otherPatch =
395 mesh_.boundaryMesh().whichPatch(otherFaceI);
396 if (thisPatch > otherPatch)
399 thisPatch = otherPatch;
408 void Foam::particle::locate
413 const bool boundaryFail,
414 const string& boundaryMsg
422 celli_ = mesh_.cellTree().findInside(position);
427 <<
"Cell not found for particle position " << position <<
"."
435 coordinates_ =
barycentric(1-3*SMALL, SMALL, SMALL, SMALL);
436 tetFacei_ = mesh_.cells()[celli_][0];
441 track(position - this->position(), 0);
452 <<
" when tracking from centre " << mesh_.cellCentres()[celli_]
453 <<
" of cell " << celli_ <<
" to position " << position
464 const polyPatch&
p = mesh_.boundaryMesh()[
patch()];
469 const vector s = position - mesh_.cellCentres()[celli_];
474 tetFacei_ = mesh_.cells()[celli_][0];
481 static label nWarnings = 0;
482 static const label maxNWarnings = 100;
483 if (nWarnings < maxNWarnings)
486 <<
" when tracking from centre " << mesh_.cellCentres()[celli_]
487 <<
" of cell " << celli_ <<
" to position " << position
491 if (nWarnings == maxNWarnings)
494 <<
"Suppressing any further warnings about particles being "
495 <<
"located outside of the mesh." <<
endl;
509 const label tetFacei,
521 origId_(getNewParticleID())
533 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
540 origId_(getNewParticleID())
548 "Particle initialised with a location outside of the mesh"
558 const label tetFacei,
564 coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
571 origId_(getNewParticleID())
581 "Particle initialised with a location outside of the mesh"
590 coordinates_(
p.coordinates_),
592 tetFacei_(
p.tetFacei_),
595 stepFraction_(
p.stepFraction_),
596 origProc_(
p.origProc_),
604 coordinates_(
p.coordinates_),
606 tetFacei_(
p.tetFacei_),
609 stepFraction_(
p.stepFraction_),
610 origProc_(
p.origProc_),
619 const vector& displacement,
620 const scalar fraction
623 scalar
f = trackToFace(displacement, fraction);
625 while (onInternalFace())
629 f *= trackToFace(
f*displacement,
f*fraction);
638 const vector& displacement,
639 const scalar fraction
644 label tetTriI = onFace() ? 0 : -1;
650 f *= trackToTri(
f*displacement,
f*fraction, tetTriI);
657 else if (tetTriI == 0)
674 const vector& displacement,
675 const scalar fraction,
679 const vector x0 = position();
680 const vector x1 = displacement;
685 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
686 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
693 stationaryTetReverseTransform(centre, detA,
T);
698 stationaryTetGeometry(o,
b, v1, v2);
699 Pout<<
"Tet points o=" << o <<
", b=" <<
b
700 <<
", v1=" << v1 <<
", v2=" << v2 <<
endl
701 <<
"Tet determinant = " << detA <<
endl
702 <<
"Start local coordinates = " <<
y0 <<
endl;
706 const scalar
f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
713 Pout<<
"Local displacement = " << Tx1 <<
"/" << detA <<
endl;
718 scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
719 for (label i = 0; i < 4; ++ i)
721 if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
723 scalar
mu = -
y0[i]/Tx1[i];
727 Pout<<
"Hit on tet face " << i <<
" at local coordinate "
728 <<
y0 +
mu*Tx1 <<
", " <<
mu*detA*100 <<
"% of the "
729 <<
"way along the track" <<
endl;
732 if (0 <=
mu &&
mu < muH)
744 for (label i = 0; i < 4; ++ i)
763 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
767 Pout<<
"Track hit no tet faces" <<
endl;
769 Pout<<
"End local coordinates = " << yH <<
endl
770 <<
"End global coordinates = " << position() <<
endl
771 <<
"Tracking displacement = " << position() - x0 <<
endl
772 << muH*detA*100 <<
"% of the step from " << stepFraction_ <<
" to "
773 << stepFraction_ + fraction <<
" completed" <<
endl <<
endl;
777 stepFraction_ += fraction*muH*detA;
779 return iH != -1 ? 1 - muH*detA : 0;
785 const vector& displacement,
786 const scalar fraction,
790 const vector x0 = position();
791 const vector x1 = displacement;
796 Pout<<
"Particle " << origId() <<
endl <<
"Tracking from " << x0
797 <<
" along " << x1 <<
" to " << x0 + x1 <<
endl;
804 movingTetReverseTransform(fraction, centre, detA,
T);
807 const scalar
f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
810 const vector x0Rel = x0 - centre[0];
811 const vector x1Rel =
f*x1 - centre[1];
814 cubicEqn detAEqn(
sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
817 ((x1Rel &
T[2]) + detA[3]*yC)*
sqr(detA[0]);
819 ((x1Rel &
T[1]) + (x0Rel &
T[2]) + detA[2]*yC)*detA[0];
821 ((x1Rel &
T[0]) + (x0Rel &
T[1]) + detA[1]*yC);
826 hitEqn[i] =
cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
831 scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
832 for (label i = 0; i < 4; ++i)
836 for (label j = 0; j < 3; ++j)
840 if (0 <=
mu[j] &&
mu[j] < muH)
852 hitEqn[0].value(muH),
853 hitEqn[1].value(muH),
854 hitEqn[2].value(muH),
864 const scalar detAH = detAEqn.
value(muH);
865 if (!std::isnormal(detAH))
868 <<
"A moving tet collapsed onto a particle. This is not supported. "
869 <<
"The mesh is too poor, or the motion too severe, for particle "
875 for (label i = 0; i < 4; ++ i)
891 stepFraction_ += fraction*muH*detA[0];
897 Pout<<
"Track hit tet face " << iH <<
" first" <<
endl;
901 Pout<<
"Track hit no tet faces" <<
endl;
903 Pout<<
"End local coordinates = " << yH <<
endl
904 <<
"End global coordinates = " << position() <<
endl;
907 return iH != -1 ? 1 - muH*detA[0] : 0;
913 const vector& displacement,
914 const scalar fraction,
920 return trackToMovingTri(displacement, fraction, tetTriI);
924 return trackToStationaryTri(displacement, fraction, tetTriI);
931 if (
cmptMin(mesh_.geometricD()) == -1)
955 facei_ = mesh_.boundaryMesh()[
patch()].whichFace(facei_);
966 refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
976 transformProperties(
T);
986 transformProperties(-
s);
991 facei_ += ppp.
start();
996 tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
1037 const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
1041 tetFacei_ = mesh_.cells()[celli_][0];
1057 movingTetReverseTransform(0, centre, detA,
T);
1058 coordinates_ += (
pos - centre[0]) &
T[0]/detA[0];
1065 stationaryTetReverseTransform(centre, detA,
T);
1066 coordinates_ += (
pos - centre) &
T/detA;
1074 const label procCell,
1075 const label procTetFace
1084 (mesh_.faceOwner()[tetFacei_] == celli_)
1085 == (procMesh.
faceOwner()[procTetFace] == procCell)
1092 return procMesh.
faces()[procTetFace].size() - 1 - tetPti_;
1109 "Particle mapped to a location outside of the mesh"
1122 "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.
prefixOSstream Pout
An Ostream wrapper for parallel output to std::cout.
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.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
tmp< faMatrix< Type > > operator==(const faMatrix< Type > &, const faMatrix< Type > &)
scalar value(const scalar x) const
Evaluate the cubic equation 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)
Return the sum of components of a SphericalTensor.
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
Container to encapsulate various operations for cubic equation of the forms with real coefficients:
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)