zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

particle: Improvements to tracking across AMIs

Browse Source 
The projection direction has been corrected to include the effect of
mesh motion. The case where the source and receiving faces are of
differing orientations and the particle displacement points into both is
now detected and handled.
This commit is contained in:
Will Bainbridge
2018-06-26 16:48:53 +01:00
parent fb4cb6e57a
commit 35dd88e553
15 changed files with 293 additions and 128 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
src/functionObjects/field/nearWallFields/findCellParticle.C
View File

@ -176,9 +176,10 @@ void Foam::findCellParticle::hitCyclicPatch
void Foam::findCellParticle::hitCyclicAMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData& td,
const vector&
trackingData& td
)
{
// Remove particle
@ -188,9 +189,23 @@ void Foam::findCellParticle::hitCyclicAMIPatch
void Foam::findCellParticle::hitCyclicACMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData& td,
const vector&
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::findCellParticle::hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData& td
)
{
// Remove particle

20
src/functionObjects/field/nearWallFields/findCellParticle.H
View File

@ -228,17 +228,29 @@ public:
//- Overridable function to handle the particle hitting a cyclicAMI
void hitCyclicAMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData&,
const vector&
trackingData&
);
//- Overridable function to handle the particle hitting a cyclicACMI
void hitCyclicACMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData&,
const vector&
trackingData&
);
//- Overridable function to handle the particle hitting a
//cyclicRepeatAMI
void hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
Cloud<findCellParticle>&,
trackingData&
);
//- Overridable function to handle the particle hitting a

23
src/functionObjects/field/streamLine/streamLineParticle.C
View File

@ -323,9 +323,10 @@ void Foam::streamLineParticle::hitCyclicPatch
void Foam::streamLineParticle::hitCyclicAMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData& td,
const vector&
trackingData& td
)
{
// Remove particle
@ -335,9 +336,23 @@ void Foam::streamLineParticle::hitCyclicAMIPatch
void Foam::streamLineParticle::hitCyclicACMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData& td,
const vector&
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::streamLineParticle::hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData& td
)
{
// Remove particle

20
src/functionObjects/field/streamLine/streamLineParticle.H
View File

@ -231,18 +231,30 @@ public:
// cyclicAMIPatch
void hitCyclicAMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData&,
const vector& direction
trackingData&
);
//- Overridable function to handle the particle hitting a
// cyclicACMIPatch
void hitCyclicACMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData&,
const vector& direction
trackingData&
);
//- Overridable function to handle the particle hitting a
// cyclicRepeatAMIPatch
void hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
streamLineParticleCloud&,
trackingData&
);
//- Overridable function to handle the particle hitting a

28
src/lagrangian/basic/particle/particle.C
View File

@ -409,7 +409,6 @@ void Foam::particle::changeToMasterPatch()
void Foam::particle::locate
(
const vector& position,
const vector* direction,
label celli,
const bool boundaryFail,
const string boundaryMsg
@ -488,34 +487,11 @@ void Foam::particle::locate
}
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 sN = (displacement & n)*n;
const vector sT = displacement - sN;
coordinates_ = barycentric(1, 0, 0, 0);
celli_ = celli;
tetFacei_ = minTetFacei;
tetPti_ = minTetPti;
facei_ = -1;
track(sT, 0);
track(sN, 0);
static label nWarnings = 0;
static const label maxNWarnings = 100;
if (nWarnings < maxNWarnings)
{
WarningInFunction << boundaryMsg << endl;
WarningInFunction << boundaryMsg.c_str() << endl;
++ nWarnings;
}
if (nWarnings == maxNWarnings)
@ -572,7 +548,6 @@ Foam::particle::particle
locate
(
position,
nullptr,
celli,
false,
"Particle initialised with a location outside of the mesh."
@ -1171,7 +1146,6 @@ void Foam::particle::autoMap
locate
(
position,
nullptr,
mapper.reverseCellMap()[celli_],
true,
"Particle mapped to a location outside of the mesh."

35
src/lagrangian/basic/particle/particle.H
View File

@ -274,7 +274,6 @@ private:
void locate
(
const vector& position,
const vector* direction,
label celli,
const bool boundaryFail,
const string boundaryMsg
@ -309,21 +308,34 @@ protected:
//- Overridable function to handle the particle hitting a cyclicAMIPatch
template<class TrackCloudType>
void hitCyclicAMIPatch(TrackCloudType&, trackingData&, const vector&);
void hitCyclicAMIPatch
(
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
);
//- Overridable function to handle the particle hitting a
// cyclicACMIPatch
template<class TrackCloudType>
void hitCyclicACMIPatch(TrackCloudType&, trackingData&, const vector&);
void hitCyclicACMIPatch
(
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
);
//- Overridable function to handle the particle hitting an
// cyclicRepeatAMIPolyPatch
template<class TrackCloudType>
void hitCyclicRepeatAMIPatch
(
TrackCloudType&,
trackingData&,
const vector&
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
);
//- Overridable function to handle the particle hitting a processorPatch
@ -559,16 +571,17 @@ public:
template<class TrackCloudType>
void hitFace
(
const vector& direction,
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
);
//- Convenience function. Cobines trackToFace and hitFace
//- Convenience function. Combines trackToFace and hitFace
template<class TrackCloudType>
scalar trackToAndHitFace
(
const vector& direction,
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
@ -582,8 +595,8 @@ public:
// Patch data
//- Get the normal and velocity of the current patch location
inline void patchData(vector& n, vector& U) const;
//- Get the normal and displacement of the current patch location
inline void patchData(vector& normal, vector& displacement) const;
// Transformations

14
src/lagrangian/basic/particle/particleI.H
View File

@ -286,7 +286,7 @@ inline Foam::vector Foam::particle::position() const
}
void Foam::particle::patchData(vector& n, vector& U) const
void Foam::particle::patchData(vector& normal, vector& displacement) const
{
if (!onBoundaryFace())
{
@ -300,27 +300,23 @@ void Foam::particle::patchData(vector& n, vector& U) const
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);
n = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
normal = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
// Interpolate the motion of the three face vertices to the current
// coordinates
U =
displacement =
coordinates_.b()*base[1]
+ coordinates_.c()*vertex1[1]
+ coordinates_.d()*vertex2[1];
// The movingTetGeometry method gives the motion as a displacement
// across the time-step, so we divide by the time-step to get velocity
U /= mesh_.time().deltaTValue();
}
else
{
vector centre, base, vertex1, vertex2;
stationaryTetGeometry(centre, base, vertex1, vertex2);
n = triPointRef(base, vertex1, vertex2).normal();
normal = triPointRef(base, vertex1, vertex2).normal();
U = Zero;
displacement = Zero;
}
}

155
src/lagrangian/basic/particle/particleTemplates.C
View File

@ -103,7 +103,8 @@ void Foam::particle::writeFields(const TrackCloudType& c)
template<class TrackCloudType>
void Foam::particle::hitFace
(
const vector& direction,
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
)
@ -152,15 +153,15 @@ void Foam::particle::hitFace
}
else if (isA<cyclicACMIPolyPatch>(patch))
{
p.hitCyclicACMIPatch(cloud, ttd, direction);
p.hitCyclicACMIPatch(displacement, fraction, cloud, ttd);
}
else if (isA<cyclicAMIPolyPatch>(patch))
{
p.hitCyclicAMIPatch(cloud, ttd, direction);
p.hitCyclicAMIPatch(displacement, fraction, cloud, ttd);
}
else if (isA<cyclicRepeatAMIPolyPatch>(patch))
{
p.hitCyclicRepeatAMIPatch(cloud, ttd, direction);
p.hitCyclicRepeatAMIPatch(displacement, fraction, cloud, ttd);
}
else if (isA<processorPolyPatch>(patch))
{
@ -182,7 +183,7 @@ void Foam::particle::hitFace
template<class TrackCloudType>
Foam::scalar Foam::particle::trackToAndHitFace
(
const vector& direction,
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
@ -193,9 +194,9 @@ Foam::scalar Foam::particle::trackToAndHitFace
Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
}
const scalar f = trackToFace(direction, fraction);
const scalar f = trackToFace(displacement, fraction);
hitFace(direction, cloud, td);
hitFace(displacement, fraction, cloud, td);
return f;
}
@ -282,30 +283,50 @@ void Foam::particle::hitCyclicPatch(TrackCloudType&, trackingData&)
template<class TrackCloudType>
void Foam::particle::hitCyclicAMIPatch
(
TrackCloudType&,
trackingData& td,
const vector& direction
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td
)
{
vector pos = position();
const cyclicAMIPolyPatch& cpp =
static_cast<const cyclicAMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
const label sendFacei = cpp.whichFace(facei_);
const labelPair receiveIs = cpp.pointAMIAndFace(sendFacei, direction, pos);
if (debug)
{
Info<< "Particle " << origId() << " crossing AMI from " << cpp.name()
<< " to " << receiveCpp.name() << endl << endl;
}
// Get the send patch data
vector sendNormal, sendDisplacement;
patchData(sendNormal, sendDisplacement);
vector pos = position();
const labelPair receiveIs =
cpp.pointAMIAndFace
(
cpp.whichFace(facei_),
displacement - fraction*sendDisplacement,
pos
);
const label receiveAMIi = receiveIs.first();
const label receiveFacei = receiveIs.second();
// If the receiving face could not be found then issue a warning and remove
// the particle
if (receiveFacei < 0)
{
// If the patch face of the particle is not known assume that the
// particle is lost and mark it to be deleted.
td.keepParticle = false;
WarningInFunction
<< "Particle lost across " << cyclicAMIPolyPatch::typeName
<< " patches " << cpp.name() << " and " << receiveCpp.name()
<< " at position " << pos << endl;
<< "Particle transfer from " << cyclicAMIPolyPatch::typeName
<< " patches " << cpp.name() << " to " << receiveCpp.name()
<< " failed at position " << pos << " and with displacement "
<< (displacement - fraction*sendDisplacement) << nl
<< " A receiving face could not be found" << nl
<< " The particle has been removed" << nl << endl;
return;
}
@ -313,12 +334,9 @@ void Foam::particle::hitCyclicAMIPatch
facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
// Locate the particle on the receiving side
vector directionT = direction;
cpp.reverseTransformDirection(directionT, sendFacei);
locate
(
pos,
&directionT,
mesh_.faceOwner()[facei_],
false,
"Particle crossed between " + cyclicAMIPolyPatch::typeName +
@ -330,14 +348,58 @@ void Foam::particle::hitCyclicAMIPatch
// register as incomplete
facei_ = tetFacei_;
// Get the receive patch data
vector receiveNormal, receiveDisplacement;
if (onBoundaryFace())
{
patchData(receiveNormal, receiveDisplacement);
}
else
{
receiveNormal = receiveCpp.faceNormals()[receiveFacei];
}
// If the displacement points into the receiving face then issue a warning
// and remove the particle
if
(
onBoundaryFace()
&& ((displacement - 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 "
<< (displacement - 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;
}
// Transform the properties
const vectorTensorTransform AMITransform =
receiveCpp.owner()
? receiveCpp.AMITransforms()[receiveAMIi]
: inv(cpp.AMITransforms()[receiveAMIi]);
if (AMITransform.hasR())
{
transformProperties(AMITransform.R());
}
else if (AMITransform.t() != vector::zero)
{
transformProperties(AMITransform.t());
}
if (!receiveCpp.parallel())
{
const tensor& T =
(
receiveCpp.forwardT().size() == 1
? receiveCpp.forwardT()[0]
: receiveCpp.forwardT()[receiveFacei]
: receiveCpp.forwardT()[facei_]
);
transformProperties(T);
}
@ -347,36 +409,31 @@ void Foam::particle::hitCyclicAMIPatch
(
(receiveCpp.separation().size() == 1)
? receiveCpp.separation()[0]
: receiveCpp.separation()[receiveFacei]
: receiveCpp.separation()[facei_]
);
transformProperties(-s);
}
const vectorTensorTransform& T =
cpp.owner()
? cpp.AMITransforms()[receiveAMIi]
: cpp.neighbPatch().AMITransforms()[receiveAMIi];
if (T.hasR())
{
transformProperties(T.R());
}
else if (T.t() != vector::zero)
{
transformProperties(T.t());
}
}
template<class TrackCloudType>
void Foam::particle::hitCyclicACMIPatch
(
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td,
const vector& direction
trackingData& td
)
{
typename TrackCloudType::particleType& p =
static_cast<typename TrackCloudType::particleType&>(*this);
const cyclicACMIPolyPatch& cpp =
static_cast<const cyclicACMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
vector patchNormal, patchDisplacement;
patchData(patchNormal, patchDisplacement);
const label localFacei = cpp.whichFace(facei_);
// If the mask is within the patch tolerance at either end, then we can
@ -391,20 +448,26 @@ void Foam::particle::hitCyclicACMIPatch
if (!couple && !nonOverlap)
{
vector pos = position();
couple = cpp.pointAMIAndFace(localFacei, direction, pos).first() >= 0;
couple =
cpp.pointAMIAndFace
(
localFacei,
displacement - fraction*patchDisplacement,
pos
).first() >= 0;
nonOverlap = !couple;
}
if (couple)
{
hitCyclicAMIPatch(cloud, td, direction);
p.hitCyclicAMIPatch(displacement, fraction, cloud, td);
}
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);
p.hitFace(displacement, fraction, cloud, td);
}
}
@ -412,12 +475,16 @@ void Foam::particle::hitCyclicACMIPatch
template<class TrackCloudType>
void Foam::particle::hitCyclicRepeatAMIPatch
(
const vector& displacement,
const scalar fraction,
TrackCloudType& cloud,
trackingData& td,
const vector& direction
trackingData& td
)
{
hitCyclicAMIPatch(cloud, td, direction);
typename TrackCloudType::particleType& p =
static_cast<typename TrackCloudType::particleType&>(*this);
p.hitCyclicAMIPatch(displacement, fraction, cloud, td);
}

1
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -706,6 +706,7 @@ void Foam::KinematicCloud<CloudType>::patchData
) const
{
p.patchData(nw, Up);
Up /= p.mesh().time().deltaTValue();
// If this is a wall patch, then there may be a non-zero tangential velocity
// component; the lid velocity in a lid-driven cavity case, for example. We

2
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -358,7 +358,7 @@ bool Foam::KinematicParcel<ParcelType>::move
if (p.onFace() && ttd.keepParticle)
{
p.hitFace(s, cloud, ttd);
p.hitFace(f*s - d, f, cloud, ttd);
}
}

23
src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
View File

@ -206,9 +206,10 @@ void Foam::trackedParticle::hitCyclicPatch
void Foam::trackedParticle::hitCyclicAMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData& td,
const vector& direction
trackingData& td
)
{
// Remove particle
@ -218,9 +219,23 @@ void Foam::trackedParticle::hitCyclicAMIPatch
void Foam::trackedParticle::hitCyclicACMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData& td,
const vector&
trackingData& td
)
{
// Remove particle
td.keepParticle = false;
}
void Foam::trackedParticle::hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData& td
)
{
// Remove particle

20
src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
View File

@ -263,17 +263,29 @@ public:
//- Overridable function to handle the particle hitting a cyclicAMI
void hitCyclicAMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData&,
const vector&
trackingData&
);
//- Overridable function to handle the particle hitting a cyclicACMI
void hitCyclicACMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData&,
const vector&
trackingData&
);
//- Overridable function to handle the particle hitting a
// cyclicRepeatAMI
void hitCyclicRepeatAMIPatch
(
const vector&,
const scalar,
Cloud<trackedParticle>&,
trackingData&
);
//- Overridable function to handle the particle hitting a

30
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
View File

@ -1232,7 +1232,6 @@ const
const labelList& addr = tgtAddress_[tgtFacei];
pointHit nearest;
nearest.setDistance(great);
label nearestFacei = -1;
forAll(addr, i)
@ -1240,23 +1239,27 @@ const
const label srcFacei = addr[i];
const face& f = srcPatch[srcFacei];
pointHit ray = f.ray(tgtPoint, n, srcPoints);
const pointHit ray =
f.ray(tgtPoint, n, srcPoints, intersection::VISIBLE);
if (ray.hit())
{
// tgtPoint = ray.rawPoint();
tgtPoint = ray.rawPoint();
return srcFacei;
}
else if (ray.distance() < nearest.distance())
const pointHit near = f.nearestPoint(tgtPoint, srcPoints);
if (near.distance() < nearest.distance())
{
nearest = ray;
nearest = near;
nearestFacei = srcFacei;
}
}
if (nearest.hit() || nearest.eligibleMiss())
{
// tgtPoint = nearest.rawPoint();
tgtPoint = nearest.rawPoint();
return nearestFacei;
}
@ -1277,7 +1280,6 @@ const
const pointField& tgtPoints = tgtPatch.points();
pointHit nearest;
nearest.setDistance(great);
label nearestFacei = -1;
// Target face addresses that intersect source face srcFacei
@ -1288,23 +1290,27 @@ const
const label tgtFacei = addr[i];
const face& f = tgtPatch[tgtFacei];
pointHit ray = f.ray(srcPoint, n, tgtPoints);
const pointHit ray =
f.ray(srcPoint, n, tgtPoints, intersection::VISIBLE);
if (ray.hit())
{
// srcPoint = ray.rawPoint();
srcPoint = ray.rawPoint();
return tgtFacei;
}
else if (ray.distance() < nearest.distance())
const pointHit near = f.nearestPoint(srcPoint, tgtPoints);
if (near.distance() < nearest.distance())
{
nearest = ray;
nearest = near;
nearestFacei = tgtFacei;
}
}
if (nearest.hit() || nearest.eligibleMiss())
{
// srcPoint = nearest.rawPoint();
srcPoint = nearest.rawPoint();
return nearestFacei;
}

20
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
View File

@ -952,6 +952,26 @@ void Foam::cyclicAMIPolyPatch::transformPosition
}
void Foam::cyclicAMIPolyPatch::transformDirection
(
vector& d,
const label facei
) const
{
if (!parallel())
{
const tensor& T =
(
forwardT().size() == 1
? forwardT()[0]
: forwardT()[facei]
);
d = Foam::transform(T, d);
}
}
void Foam::cyclicAMIPolyPatch::reverseTransformPosition
(
point& l,

7
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
View File

@ -344,6 +344,13 @@ public:
const label facei
) const;
//- Transform a patch-based direction from nbr side to this side
virtual void transformDirection
(
vector& d,
const label facei
) const;
//- Transform a patch-based position from this side to nbr side
virtual void reverseTransformPosition
(