diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 6588607fbb283331e8d2bac294543a6a665ecb7b..511642690fd095b1ad756e587dedc6ff2d9a3893 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2018 OpenCFD Ltd.
+ Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -40,7 +40,7 @@ namespace Foam
defineTypeNameAndDebug(particle, 0);
}
-const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
+const Foam::label Foam::particle::maxNBehind_ = 10;
Foam::label Foam::particle::particleCount_ = 0;
@@ -365,7 +365,6 @@ void Foam::particle::changeCell()
const label ownerCellI = mesh_.faceOwner()[tetFacei_];
const bool isOwner = celli_ == ownerCellI;
celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
-
// Reflect to account for the change of triangle orientation in the new cell
reflect();
}
@@ -414,88 +413,95 @@ void Foam::particle::locate
const string& boundaryMsg
)
{
+ if (debug)
+ {
+ Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
+ }
+
celli_ = celli;
// Find the cell, if it has not been given
- if (celli_ < 0)
+ if (celli < 0)
{
celli_ = mesh_.cellTree().findInside(position);
}
- if (celli_ < 0)
+ if (celli < 0)
{
FatalErrorInFunction
<< "Cell not found for particle position " << position << "."
<< exit(FatalError);
}
- // Put the particle at (almost) the cell centre and in a random tet.
- // Note perturbing the cell centre to make sure we find at least one
- // tet containing it. With start point exactly at the cell centre very
- // occasionally it would not get found in any of the tets
- coordinates_ = barycentric(1-3*SMALL, SMALL, SMALL, SMALL);
- tetFacei_ = mesh_.cells()[celli_][0];
- tetPti_ = 1;
- facei_ = -1;
+ const vector displacement = position - mesh_.cellCentres()[celli_];
- // Track to the injection point
- track(position - this->position(), 0);
+ // Loop all cell tets to find the one containing the position. Track
+ // through each tet from the cell centre. If a tet contains the position
+ // then the track will end with a single trackToTri.
+ const class cell& c = mesh_.cells()[celli_];
+ scalar minF = VGREAT;
+ label minTetFacei = -1, minTetPti = -1;
+ forAll(c, cellTetFacei)
+ {
+ const class face& f = mesh_.faces()[c[cellTetFacei]];
+ for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
+ {
+ coordinates_ = barycentric(1, 0, 0, 0);
+ tetFacei_ = c[cellTetFacei];
+ tetPti_ = tetPti;
+ facei_ = -1;
+
+ label tetTriI = -1;
+ const scalar f = trackToTri(displacement, 0, tetTriI);
+
+ if (tetTriI == -1)
+ {
+ return;
+ }
+
+ if (f < minF)
+ {
+ minF = f;
+ minTetFacei = tetFacei_;
+ minTetPti = tetPti_;
+ }
+ }
+ }
+ // The particle must be (hopefully only slightly) outside the cell. Track
+ // into the tet which got the furthest.
+ coordinates_ = barycentric(1, 0, 0, 0);
+ tetFacei_ = minTetFacei;
+ tetPti_ = minTetPti;
+ facei_ = -1;
+ track(displacement, 0);
if (!onFace())
{
return;
}
- // We hit a boundary ...
+ // If we are here then we hit a boundary
if (boundaryFail)
{
- FatalErrorInFunction << boundaryMsg
- << " when tracking from centre " << mesh_.cellCentres()[celli_]
- << " of cell " << celli_ << " to position " << position
- << exit(FatalError);
+ FatalErrorInFunction << boundaryMsg << exit(FatalError);
}
else
{
- // Re-do the track, but this time do the bit tangential to the
- // direction/patch first. This gets us as close as possible to the
- // original path/position.
-
- if (direction == nullptr)
- {
- const polyPatch& p = mesh_.boundaryMesh()[patch()];
- direction = &p.faceNormals()[p.whichFace(facei_)];
- }
-
- const vector n = *direction/mag(*direction);
- const vector s = position - mesh_.cellCentres()[celli_];
- const vector sN = (s & n)*n;
- const vector sT = s - sN;
-
- coordinates_ = barycentric(1, 0, 0, 0);
- tetFacei_ = mesh_.cells()[celli_][0];
- tetPti_ = 1;
- facei_ = -1;
-
- track(sT, 0);
- track(sN, 0);
-
static label nWarnings = 0;
static const label maxNWarnings = 100;
- if (nWarnings < maxNWarnings)
+ if ((nWarnings < maxNWarnings) && boundaryFail)
{
- WarningInFunction << boundaryMsg.c_str()
- << " when tracking from centre " << mesh_.cellCentres()[celli_]
- << " of cell " << celli_ << " to position " << position
- << endl;
- ++nWarnings;
+ WarningInFunction << boundaryMsg.c_str() << endl;
+ ++ nWarnings;
}
if (nWarnings == maxNWarnings)
{
WarningInFunction
<< "Suppressing any further warnings about particles being "
<< "located outside of the mesh." << endl;
- ++nWarnings;
+ ++ nWarnings;
}
}
+
}
@@ -516,7 +522,9 @@ Foam::particle::particle
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
- stepFraction_(0.0),
+ stepFraction_(1.0),
+ behind_(0.0),
+ nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{}
@@ -535,7 +543,9 @@ Foam::particle::particle
tetFacei_(-1),
tetPti_(-1),
facei_(-1),
- stepFraction_(0.0),
+ stepFraction_(1.0),
+ behind_(0.0),
+ nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
@@ -566,7 +576,9 @@ Foam::particle::particle
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
- stepFraction_(0.0),
+ stepFraction_(1.0),
+ behind_(0.0),
+ nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
@@ -593,6 +605,8 @@ Foam::particle::particle(const particle& p)
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
+ behind_(p.behind_),
+ nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@@ -607,6 +621,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
+ behind_(p.behind_),
+ nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
@@ -645,7 +661,8 @@ Foam::scalar Foam::particle::trackToFace
facei_ = -1;
- while (true)
+ // Loop the tets in the current cell
+ while (nBehind_ < maxNBehind_)
{
f *= trackToTri(f*displacement, f*fraction, tetTriI);
@@ -666,6 +683,25 @@ Foam::scalar Foam::particle::trackToFace
changeTet(tetTriI);
}
}
+
+ // Warn if stuck, and incorrectly advance the step fraction to completion
+ static label stuckID = -1, stuckProc = -1;
+ if (origId_ != stuckID && origProc_ != stuckProc)
+ {
+ WarningInFunction
+ << "Particle #" << origId_ << " got stuck at " << position()
+ << endl;
+ }
+
+ stuckID = origId_;
+ stuckProc = origProc_;
+
+ stepFraction_ += f*fraction;
+
+ behind_ = 0;
+ nBehind_ = 0;
+
+ return 0;
}
@@ -702,11 +738,8 @@ Foam::scalar Foam::particle::trackToStationaryTri
<< "Start local coordinates = " << y0 << endl;
}
- // Get the factor by which the displacement is increased
- const scalar f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
-
// Calculate the local tracking displacement
- barycentric Tx1(f*x1 & T);
+ barycentric Tx1(x1 & T);
if (debug)
{
@@ -718,7 +751,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
for (label i = 0; i < 4; ++ i)
{
- if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
+ if (Tx1[i] < - detA*SMALL)
{
scalar mu = - y0[i]/Tx1[i];
@@ -776,6 +809,30 @@ Foam::scalar Foam::particle::trackToStationaryTri
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA;
+ if (debug)
+ {
+ Info << "Step Fraction : " << stepFraction_ << endl;
+ Info << "fraction*muH*detA : " << fraction*muH*detA << endl;
+ Info << "static muH : " << muH << endl;
+ Info << "origId() : " << origId() << endl;
+ }
+
+ // Accumulate displacement behind
+ if (detA <= 0 || nBehind_ > 0)
+ {
+ behind_ += muH*detA*mag(displacement);
+
+ if (behind_ > 0)
+ {
+ behind_ = 0;
+ nBehind_ = 0;
+ }
+ else
+ {
+ ++ nBehind_;
+ }
+ }
+
return iH != -1 ? 1 - muH*detA : 0;
}
@@ -803,12 +860,20 @@ Foam::scalar Foam::particle::trackToMovingTri
FixedList<barycentricTensor, 3> T;
movingTetReverseTransform(fraction, centre, detA, T);
- // Get the factor by which the displacement is increased
- const scalar f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
+ if (debug)
+ {
+ Pair<vector> o, b, v1, v2;
+ movingTetGeometry(fraction, o, b, v1, v2);
+ Info<< "Tet points o=" << o[0] << ", b=" << b[0]
+ << ", v1=" << v1[0] << ", v2=" << v2[0] << endl
+ << "Tet determinant = " << detA[0] << endl
+ << "Start local coordinates = " << y0[0] << endl;
+ }
+
// Get the relative global position
const vector x0Rel = x0 - centre[0];
- const vector x1Rel = f*x1 - centre[1];
+ const vector x1Rel = x1 - centre[1];
// Form the determinant and hit equations
cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
@@ -826,6 +891,16 @@ Foam::scalar Foam::particle::trackToMovingTri
hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
}
+ if (debug)
+ {
+ for (label i = 0; i < 4; ++ i)
+ {
+ Info<< (i ? " " : "Hit equation ") << i << " = "
+ << hitEqn[i] << endl;
+ }
+ Info<< " DetA equation = " << detA << endl;
+ }
+
// Calculate the hit fraction
label iH = -1;
scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
@@ -835,8 +910,33 @@ Foam::scalar Foam::particle::trackToMovingTri
for (label j = 0; j < 3; ++j)
{
- if (mu.type(j) == roots::real && hitEqn[i].derivative(mu[j]) < 0)
+ if
+ (
+ mu.type(j) == roots::real
+ && hitEqn[i].derivative(mu[j]) < - detA[0]*SMALL
+ )
{
+ if (debug)
+ {
+ const barycentric yH
+ (
+ hitEqn[0].value(mu[j]),
+ hitEqn[1].value(mu[j]),
+ hitEqn[2].value(mu[j]),
+ hitEqn[3].value(mu[j])
+ );
+ const scalar detAH = detAEqn.value(mu[j]);
+
+ Info<< "Hit on tet face " << i << " at local coordinate "
+ << (std::isnormal(detAH) ? name(yH/detAH) : "???")
+ << ", " << mu[j]*detA[0]*100 << "% of the "
+ << "way along the track" << endl;
+
+ Info<< "derivative : " << hitEqn[i].derivative(mu[j]) << nl
+ << " coord " << j << " mu[j]: " << mu[j] << nl
+ << " hitEq " << i << endl;
+ }
+
if (0 <= mu[j] && mu[j] < muH)
{
iH = i;
@@ -887,8 +987,26 @@ Foam::scalar Foam::particle::trackToMovingTri
coordinates_ = yH;
tetTriI = iH;
+ scalar advance = muH*detA[0];
+
// Set the proportion of the track that has been completed
- stepFraction_ += fraction*muH*detA[0];
+ stepFraction_ += fraction*advance;
+
+ // Accumulate displacement behind
+ if (detA[0] <= 0 || nBehind_ > 0)
+ {
+ behind_ += muH*detA[0]*mag(displacement);
+
+ if (behind_ > 0)
+ {
+ behind_ = 0;
+ nBehind_ = 0;
+ }
+ else
+ {
+ ++ nBehind_;
+ }
+ }
if (debug)
{
@@ -900,10 +1018,15 @@ Foam::scalar Foam::particle::trackToMovingTri
{
Pout<< "Track hit no tet faces" << endl;
}
- Pout<< "End local coordinates = " << yH << endl
- << "End global coordinates = " << position() << endl;
+// Pout<< "End local coordinates = " << yH << endl
+// << "End global coordinates = " << position() << endl
+// << "Tracking displacement = " << position() - x0 << endl
+// << muH*detA[0]*100 << "% of the step from " << stepFraction_
+// << " to " << stepFraction_ + fraction << " completed" << endl
+// << endl;
}
+
return iH != -1 ? 1 - muH*detA[0] : 0;
}
@@ -915,7 +1038,7 @@ Foam::scalar Foam::particle::trackToTri
label& tetTriI
)
{
- if (mesh_.moving())
+ if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
{
return trackToMovingTri(displacement, fraction, tetTriI);
}
@@ -1049,7 +1172,7 @@ void Foam::particle::correctAfterInteractionListReferral(const label celli)
// so this approximate topology is good enough. By using the nearby cell we
// minimise the error associated with the incorrect topology.
coordinates_ = barycentric(1, 0, 0, 0);
- if (mesh_.moving())
+ if (mesh_.moving() && stepFraction_ != 1)
{
Pair<vector> centre;
FixedList<scalar, 4> detA;
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 3345be1fa15f51c46933a6db33b237de147154ab..03f938d9544c1cd084c1bf196570d36822813dda 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2017-2019 OpenCFD Ltd.
+ Copyright (C) 2017-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -78,7 +78,7 @@ class particle
:
public IDLList<particle>::link
{
- // Private member data
+ // Private Data
//- Size in bytes of the position data
static const std::size_t sizeofPosition;
@@ -86,15 +86,9 @@ class particle
//- Size in bytes of the fields
static const std::size_t sizeofFields;
- //- The factor by which the displacement is increased when passing
- // through negative space. This should be slightly bigger than one.
- // This is needed as a straight trajectory can form a closed loop
- // through regions of overlapping positive and negative space, leading
- // to a track which never ends. By increasing the rate of displacement
- // through negative regions, the change in track fraction over this
- // part of the loop no longer exactly cancels the change over the
- // positive part, and the track therefore terminates.
- static const scalar negativeSpaceDisplacementFactor;
+ //- The value of nBehind_ at which tracking is abandoned. See the
+ // description of nBehind_.
+ static const label maxNBehind_;
public:
@@ -103,7 +97,7 @@ public:
{
public:
- // Public data
+ // Public Data
//- Flag to switch processor
bool switchProcessor;
@@ -135,7 +129,7 @@ public:
private:
- // Private data
+ // Private Data
//- Reference to the polyMesh database
const polyMesh& mesh_;
@@ -147,12 +141,12 @@ private:
label celli_;
//- Index of the face that owns the decomposed tet that the
- // particle is in
+ //- particle is in
label tetFacei_;
//- Index of the point on the face that defines the decomposed
- // tet that the particle is in. Relative to the face base
- // point.
+ //- tet that the particle is in. Relative to the face base
+ //- point.
label tetPti_;
//- Face index if the particle is on a face otherwise -1
@@ -161,6 +155,18 @@ private:
//- Fraction of time-step completed
scalar stepFraction_;
+ //- The distance behind the maximum distance reached so far
+ scalar behind_;
+
+ //- The number of tracks carried out that ended in a distance behind the
+ //- maximum distance reached so far. Once this reaches maxNBehind_,
+ // tracking is abandoned for the current step. This is needed because
+ // when tetrahedra are inverted a straight trajectory can form a closed
+ // loop through regions of overlapping positive and negative space.
+ // Without this break clause, such loops can result in a valid track
+ // which never ends.
+ label nBehind_;
+
//- Originating processor id
label origProc_;
@@ -331,7 +337,7 @@ protected:
public:
- // Static data members
+ // Static Data Members
//- Runtime type information
TypeName("particle");
@@ -526,6 +532,9 @@ public:
//- Return current particle position
inline vector position() const;
+ //- Reset particle data
+ inline void reset();
+
// Track
@@ -664,6 +673,7 @@ public:
void relocate(const point& position, const label celli = -1);
+
// I-O
//- Write the name representation to stream
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 512f1bd8e30f5a06cea8c13ae0b09410936a28fc..35d9598dbbb1755c884ae903256b3cd486a86bc2 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2018 OpenFOAM Foundation
+ Copyright (C) 2011-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -68,6 +69,7 @@ inline void Foam::particle::movingTetGeometry
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
+
const pointField& ptsOld = mesh_.oldPoints();
const pointField& ptsNew = mesh_.points();
@@ -75,8 +77,11 @@ inline void Foam::particle::movingTetGeometry
// we need to put a mesh_.oldCellCentres() method in for this to work. The
// values obtained from the mesh and those obtained from the cell do not
// necessarily match. See mantis #1993.
- const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
- const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+ //const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
+ //const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+
+ const vector ccOld = mesh_.oldCellCentres()[celli_];
+ const vector ccNew = mesh_.cellCentres()[celli_];
// Old and new points and cell centres are not sub-cycled. If we are sub-
// cycling, then we have to account for the timestep change here by
@@ -265,7 +270,7 @@ inline Foam::tetIndices Foam::particle::currentTetIndices() const
inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
{
- if (mesh_.moving())
+ if (mesh_.moving() && stepFraction_ != 1)
{
return movingTetTransform(0)[0];
}
@@ -312,6 +317,14 @@ inline Foam::vector Foam::particle::position() const
}
+inline void Foam::particle::reset()
+{
+ stepFraction_ = 0;
+ nBehind_ = 0;
+ behind_ = 0;
+}
+
+
void Foam::particle::patchData(vector& n, vector& U) const
{
if (!onBoundaryFace())
@@ -321,7 +334,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
<< exit(FatalError);
}
- if (mesh_.moving())
+ if ((mesh_.moving() && stepFraction_ != 1))
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index d9f53e383eb7967ef12036839650cd9acc0e146e..5509647fd9c32a6084fec7a84cefa6132d6a3dc7 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2016-2019 OpenCFD Ltd.
+ Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -293,7 +293,7 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
template<class TrackCloudType>
void Foam::particle::hitFace
(
- const vector& direction,
+ const vector& displacement,
TrackCloudType& cloud,
trackingData& td
)
@@ -337,11 +337,11 @@ void Foam::particle::hitFace
}
else if (isA<cyclicACMIPolyPatch>(patch))
{
- p.hitCyclicACMIPatch(cloud, ttd, direction);
+ p.hitCyclicACMIPatch(cloud, ttd, displacement);
}
else if (isA<cyclicAMIPolyPatch>(patch))
{
- p.hitCyclicAMIPatch(cloud, ttd, direction);
+ p.hitCyclicAMIPatch(cloud, ttd, displacement);
}
else if (isA<processorPolyPatch>(patch))
{
@@ -459,7 +459,7 @@ void Foam::particle::hitCyclicAMIPatch
(
TrackCloudType&,
trackingData& td,
- const vector& direction
+ const vector& displacement
)
{
vector pos = position();
@@ -468,7 +468,14 @@ void Foam::particle::hitCyclicAMIPatch
static_cast<const cyclicAMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
const label sendFacei = cpp.whichFace(facei_);
- const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
+ const label receiveFacei = cpp.pointFace(sendFacei, displacement, pos);
+
+ if (false)
+ {
+ Info<< "My new pos : " << pos << endl;
+ Info<< "Particle " << origId() << " crossing AMI from " << cpp.name()
+ << " to " << receiveCpp.name() << endl << endl;
+ }
if (receiveFacei < 0)
{
@@ -485,12 +492,13 @@ void Foam::particle::hitCyclicAMIPatch
facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
// Locate the particle on the receiving side
- vector directionT = direction;
- cpp.reverseTransformDirection(directionT, sendFacei);
+ vector displacementT = displacement;
+ cpp.reverseTransformDirection(displacementT, sendFacei);
+
locate
(
pos,
- &directionT,
+ &displacementT,
mesh_.faceOwner()[facei_],
false,
"Particle crossed between " + cyclicAMIPolyPatch::typeName +
@@ -523,6 +531,27 @@ void Foam::particle::hitCyclicAMIPatch
);
transformProperties(-s);
}
+
+ //if (onBoundaryFace())
+ {
+ //DebugVar("On boudanry")
+// vector receiveNormal, receiveDisplacement;
+// patchData(receiveNormal, receiveDisplacement);
+//
+// if (((displacementT - fraction*receiveDisplacement)&receiveNormal) > 0)
+// {
+// td.keepParticle = false;
+// WarningInFunction
+// << "Particle transfer from " << cyclicAMIPolyPatch::typeName
+// << " patches " << cpp.name() << " to " << receiveCpp.name()
+// << " failed at position " << pos << " and with displacement "
+// << (displacementT - fraction*receiveDisplacement) << nl
+// << " The displacement points into both the source and "
+// << "receiving faces, so the tracking cannot proceed" << nl
+// << " The particle has been removed" << nl << endl;
+// return;
+// }
+ }
}
@@ -531,7 +560,7 @@ void Foam::particle::hitCyclicACMIPatch
(
TrackCloudType& cloud,
trackingData& td,
- const vector& direction
+ const vector& displacement
)
{
const cyclicACMIPolyPatch& cpp =
@@ -551,20 +580,20 @@ void Foam::particle::hitCyclicACMIPatch
if (!couple && !nonOverlap)
{
vector pos = position();
- couple = cpp.pointFace(localFacei, direction, pos) >= 0;
+ couple = cpp.pointFace(localFacei, displacement, pos) >= 0;
nonOverlap = !couple;
}
if (couple)
{
- hitCyclicAMIPatch(cloud, td, direction);
+ hitCyclicAMIPatch(cloud, td, displacement);
}
else
{
// Move to the face associated with the non-overlap patch and redo the
// face interaction.
tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei;
- hitFace(direction, cloud, td);
+ hitFace(displacement, cloud, td);
}
}