zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: Particle tracking. Resetting tracking back to existing algorithm

Browse Source 
from experimental inside/outside cell method.
This commit is contained in:
graham
2010-04-29 16:20:09 +01:00
parent b7f8f37d8d
commit 55f06f718d
6 changed files with 169 additions and 805 deletions
Download Patch File Download Diff File Expand all files Collapse all files

212
src/lagrangian/basic/Cloud/Cloud.C
View File

@ -30,189 +30,6 @@ License
#include "mapPolyMesh.H"
#include "Time.H"
#include "OFstream.H"
#include "triPointRef.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-6;
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-4;
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 0;
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
// Calculate using face properties only
template<class ParticleType>
bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
{
const cellList& cells = pMesh().cells();
const labelList& fOwner = pMesh().faceOwner();
const vectorField& fAreas = pMesh().faceAreas();
const pointField& fCentres = pMesh().faceCentres();
const pointField& cCentres = pMesh().cellCentres();
const pointField& points = pMesh().points();
const cell& cFaces = cells[cellI];
bool concave = false;
forAll(cFaces, i)
{
label f0I = cFaces[i];
// Get f0 centre
const point& fc0 = fCentres[f0I];
const face& f0 = pMesh().faces()[f0I];
vector fn0 = fAreas[f0I];
fn0 /= (mag(fn0) + VSMALL);
// Flip normal if required so that it is always pointing out of
// the cell
if (fOwner[f0I] != cellI)
{
fn0 *= -1;
}
// Check if the cell centre is visible from the plane of the face
if (((fc0 - cCentres[cellI]) & fn0) < 0)
{
Pout<< "Face " << f0I
<< " is not visible from the centre of cell " << cellI
<< endl;
}
forAll(cFaces, j)
{
if (j != i)
{
label f1I = cFaces[j];
const face& f1 = pMesh().faces()[f1I];
// Is any vertex of f1 on wrong side of the plane of f0?
forAll(f1, f1pI)
{
label ptI = f1[f1pI];
// Skip points that are shared between f1 and f0
if (findIndex(f0, ptI) > -1)
{
continue;
}
const point& pt = points[ptI];
// If the cell is concave, the point on f1 will be
// on the outside of the plane of f0, defined by
// its centre and normal, and the angle between
// (fc0 -pt) and fn0 will be greater than 90
// degrees, so the dot product will be negative.
scalar d = ((fc0 - pt) & fn0);
if (d < 0)
{
// Concave face
concave = true;
}
}
// Check for co-planar faces, which are also treated
// as concave, as they are not strictly convex.
vector fn1 = fAreas[f1I];
fn1 /= (mag(fn1) + VSMALL);
// Flip normal if required so that it is always pointing out of
// the cell
if (fOwner[f1I] != cellI)
{
fn1 *= -1;
}
if ((fn0 & fn1) > planarCosAngle)
{
// Planar face
concave = true;
}
}
}
}
return concave;
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::calcConcaveCells() const
{
const cellList& cells = pMesh().cells();
concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
PackedBoolList& concaveCell = concaveCellPtr_();
forAll(cells, cellI)
{
// if (isConcaveCell(cellI))
// {
// concaveCell[cellI] = 1;
// // TODO: extend selection to include any point connected
// // cell if there are conflicts between cells using the two
// // different methods
// }
// Force all cells to be treated exactly
concaveCell[cellI] = 1;
// Force all cells to be treated by planes
// concaveCell[cellI] = 0;
}
// {
// // Write cells that are concave to file
// DynamicList<label> tmpConcaveCells;
// forAll(cells, cellI)
// {
// if (concaveCell[cellI])
// {
// tmpConcaveCells.append(cellI);
// }
// }
// Pout<< "Cloud<ParticleType>::calcConcaveCells() :"
// << " overall cells in mesh : " << pMesh().nCells() << nl
// << " of these concave : " << tmpConcaveCells.size() << endl;
// fileName fName = pMesh().time().path()/"concaveCells";
// Pout<< " Writing " << fName.name() << endl;
// OFstream file(fName);
// file << tmpConcaveCells;
// file.flush();
// }
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -220,14 +37,12 @@ template<class ParticleType>
Foam::Cloud<ParticleType>::Cloud
(
const polyMesh& pMesh,
const IDLList<ParticleType>& particles,
const bool concaveCheck
const IDLList<ParticleType>& particles
)
:
cloud(pMesh),
IDLList<ParticleType>(),
polyMesh_(pMesh),
concaveCheck_(concaveCheck),
particleCount_(0)
{
IDLList<ParticleType>::operator=(particles);
@ -239,43 +54,20 @@ Foam::Cloud<ParticleType>::Cloud
(
const polyMesh& pMesh,
const word& cloudName,
const IDLList<ParticleType>& particles,
const bool concaveCheck
const IDLList<ParticleType>& particles
)
:
cloud(pMesh, cloudName),
IDLList<ParticleType>(),
polyMesh_(pMesh),
concaveCheck_(concaveCheck),
particleCount_(0)
{
IDLList<ParticleType>::operator=(particles);
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::clearOut()
{
concaveCellPtr_.clear();
labels_.clearStorage();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParticleType>
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::concaveCell()
const
{
if (!concaveCellPtr_.valid())
{
calcConcaveCells();
}
return concaveCellPtr_();
}
template<class ParticleType>
Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
{

47
src/lagrangian/basic/Cloud/Cloud.H
View File

@ -39,7 +39,6 @@ SourceFiles
#include "IDLList.H"
#include "IOField.H"
#include "polyMesh.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -75,11 +74,6 @@ class Cloud
const polyMesh& polyMesh_;
const bool concaveCheck_;
//- Is the cell concave
mutable autoPtr<PackedBoolList> concaveCellPtr_;
//- Overall count of particles ever created. Never decreases.
mutable label particleCount_;
@ -89,12 +83,6 @@ class Cloud
// Private Member Functions
//- Determine if a particular cell is concave
bool isConcaveCell(const label cellI) const;
//- Analyse mesh for faces not being convex.
void calcConcaveCells() const;
//- Initialise cloud on IO constructor
void initCloud(const bool checkClass);
@ -123,20 +111,6 @@ public:
// Static data
//- The smallest allowable trackFraction in trackToFace before
// an intervention is required
static const scalar minValidTrackFraction;
//- Fraction of distance to cell centre to move a particle to
// 'rescue' it from a tracking problem
static const scalar trackingRescueTolerance;
//- Overlap tolerance for face intersections.
static const scalar intersectionTolerance;
//- cosine of angle for faces to be planar.
static const scalar planarCosAngle;
//- Name of cloud properties dictionary
static word cloudPropertiesName;
@ -147,8 +121,7 @@ public:
Cloud
(
const polyMesh& mesh,
const IDLList<ParticleType>& particles,
const bool concaveCheck = true
const IDLList<ParticleType>& particles
);
//- Construct from mesh, cloud name, and a list of particles
@ -156,8 +129,7 @@ public:
(
const polyMesh& mesh,
const word& cloudName,
const IDLList<ParticleType>& particles,
const bool concaveCheck = true
const IDLList<ParticleType>& particles
);
//- Construct from mesh by reading from file
@ -165,8 +137,7 @@ public:
Cloud
(
const polyMesh& mesh,
const bool checkClass = true,
const bool concaveCheck = true
const bool checkClass = true
);
@ -176,17 +147,10 @@ public:
(
const polyMesh& pMesh,
const word& cloudName,
const bool checkClass = true,
const bool concaveCheck = true
const bool checkClass = true
);
// Destructor
//- Destroy any demand-driven data
void clearOut();
// Member Functions
// Access
@ -226,9 +190,6 @@ public:
return IDLList<ParticleType>::size();
};
//- Whether each cell is concave (demand driven data)
const PackedBoolList& concaveCell() const;
// Iterators

8
src/lagrangian/basic/Cloud/CloudIO.C
View File

@ -135,13 +135,11 @@ template<class ParticleType>
Foam::Cloud<ParticleType>::Cloud
(
const polyMesh& pMesh,
const bool checkClass,
const bool concaveCheck
const bool checkClass
)
:
cloud(pMesh),
polyMesh_(pMesh),
concaveCheck_(concaveCheck),
particleCount_(0)
{
initCloud(checkClass);
@ -153,13 +151,11 @@ Foam::Cloud<ParticleType>::Cloud
(
const polyMesh& pMesh,
const word& cloudName,
const bool checkClass,
const bool concaveCheck
const bool checkClass
)
:
cloud(pMesh, cloudName),
polyMesh_(pMesh),
concaveCheck_(concaveCheck),
particleCount_(0)
{
initCloud(checkClass);

654
src/lagrangian/basic/Particle/Particle.C
View File

@ -37,7 +37,7 @@ License
template<class ParticleType>
void Foam::Particle<ParticleType>::findFaces
(
const vector& endPosition,
const vector& position,
DynamicList<label>& faceList
) const
{
@ -46,11 +46,10 @@ void Foam::Particle<ParticleType>::findFaces
const vector& C = mesh.cellCentres()[celli_];
faceList.clear();
forAll(faces, i)
{
label facei = faces[i];
scalar lam = lambda(C, endPosition, facei);
scalar lam = lambda(C, position, facei);
if ((lam > 0) && (lam < 1.0))
{
@ -63,7 +62,7 @@ void Foam::Particle<ParticleType>::findFaces
template<class ParticleType>
void Foam::Particle<ParticleType>::findFaces
(
const vector& endPosition,
const vector& position,
const label celli,
const scalar stepFraction,
DynamicList<label>& faceList
@ -77,7 +76,7 @@ void Foam::Particle<ParticleType>::findFaces
forAll(faces, i)
{
label facei = faces[i];
scalar lam = lambda(C, endPosition, facei, stepFraction);
scalar lam = lambda(C, position, facei, stepFraction);
if ((lam > 0) && (lam < 1.0))
{
@ -87,490 +86,6 @@ void Foam::Particle<ParticleType>::findFaces
}
template<class ParticleType>
bool Foam::Particle<ParticleType>::insideCellExact
(
const vector& testPt,
const label celli,
bool beingOnAFaceMeansOutside
) const
{
const polyMesh& mesh = cloud_.pMesh();
const labelList& faces = mesh.cells()[celli];
const vector& C = mesh.cellCentres()[celli];
label nFaceCrossings = 0;
// The vector from the cell centre to the end point
vector delta = testPt - C;
forAll (faces, i)
{
label facei = faces[i];
pointHit inter = mesh.faces()[facei].intersection
(
C,
delta,
mesh.faceCentres()[facei],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
if (inter.hit())
{
if (beingOnAFaceMeansOutside)
{
if (inter.distance() <= 1.0)
{
// This face was actually crossed.
nFaceCrossings++;
}
}
else
{
if (inter.distance() < 1.0)
{
// This face was actually crossed.
nFaceCrossings++;
}
}
}
}
// if (nFaceCrossings > 1)
// {
// Pout<< "In cell " << celli_ << " there were " << nFaceCrossings
// << " face crossings detected tracking from concave cell "
// << " centre to endPosition"
// << endl;
// }
if (nFaceCrossings % 2 == 0)
{
// Even number of face crossings, so the testPt must be in the
// cell.
return true;
}
return false;
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::trackToFaceExact
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
)
{
facei_ = -1;
const polyMesh& mesh = cloud_.pMesh();
const labelList& faces = mesh.cells()[celli_];
// Check all possible face crossings to see if they are actually
// crossed, determining if endPosition is outside the current
// cell. This allows situations where the cell is outside the
// cell to start with and enters the cell at the end of the track
// to be identified.
if (insideCellExact(endPosition, celli_, false))
{
// Even number of face crossings, so the particle must end up
// still in the cell.
position_ = endPosition;
trackFraction = 1.0;
return;
}
// The particle *must* have left the cell.
// a) It may have crossed a face not yet identified by testing
// faces using the cell centre to endPosition line, so the
// potentially crossed faces of the position to endPosition
// line must be assessed.
// b) It may have been outside the cell in the first place, and, despite
// trying to pick up more faces using a) the correct face to be crossed
// is not knowable. A best guess will be used, with the expectation that
// the tracking in the destination cell will be able to recover form a
// bad guess.
// For all face assessments, a full intersection test is required,
// as nothing can be assumed about the order of crossing the
// planes of faces.
const vector deltaPosition = endPosition - position_;
if (insideCellExact(position_, celli_, false))
{
// The particle started inside the cell and finished outside
// of it, find which face to cross
scalar tmpLambda = GREAT;
scalar correctLambda = GREAT;
forAll(faces, i)
{
label facei = faces[i];
// Use exact intersection.
// TODO: A correction is required for moving meshes to
// calculate the correct lambda value.
pointHit inter = mesh.faces()[facei].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[facei],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
if (inter.hit())
{
tmpLambda = inter.distance();
if
(
tmpLambda <= 1.0
&& tmpLambda < correctLambda
)
{
// This face is crossed before any other that has
// been found so far
correctLambda = tmpLambda;
facei_ = facei;
}
}
}
if (facei_ > -1)
{
if (cloud_.boundaryFace(facei_))
{
label patchi = patch(facei_);
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (isA<wallPolyPatch>(patch))
{
if ((mesh.faceAreas()[facei_] & deltaPosition) <= 0)
{
// The particle has hit a wall face but it is
// heading in the wrong direction with respect to
// the face normal
// Do not trigger a face hit and move the position
// towards the cell centre
const point& cc = mesh.cellCentres()[celli_];
position_ +=
Cloud<ParticleType>::trackingRescueTolerance
*(cc - position_);
facei_ = -1;
}
}
}
else
{
if (correctLambda < Cloud<ParticleType>::minValidTrackFraction)
{
// The particle is not far enough away from the
// face to decide if it is a valid crossing. Let
// it move a little without crossing the face to
// resolve the ambiguity.
facei_ = -1;
}
// If the face hit was not on a patch, add a small
// amount to the track to move it off the face, If it
// was not an ambiguous face crossing, this makes sure
// the face is not ambiguous next tracking step. If
// it was ambiguous, this should resolve it.
correctLambda += Cloud<ParticleType>::minValidTrackFraction;
}
trackFraction = correctLambda;
position_ += trackFraction*(endPosition - position_);
}
else
{
// The particle started inside of the cell and finished
// outside of it, but did not find a face to cross.
// Applying a rescuing correction.
const point& cc = mesh.cellCentres()[celli_];
position_ +=
Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
}
}
else
{
// The particle started outside of the cell. Find which cell
// it should be in.
const labelList& cPts = mesh.cellPoints(celli_);
DynamicList<label> checkedCells;
bool found = false;
forAll(cPts, cPtI)
{
label ptI = cPts[cPtI];
const labelList& pCs = mesh.pointCells(ptI);
forAll(pCs, pCI)
{
label cellI = pCs[pCI];
if (findIndex(checkedCells, cellI) == -1)
{
checkedCells.append(cellI);
if (insideCellExact(position_, cellI, false))
{
found = true;
celli_ = cellI;
break;
}
}
}
if (found)
{
break;
}
}
if (!found)
{
// Didn't find a new cell after searching point connected
// cells. Applying a rescuing correction.
const point& cc = mesh.cellCentres()[celli_];
position_ +=
Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
}
}
if (facei_ > -1)
{
faceAction(trackFraction, endPosition, td);
}
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::trackToFacePlanes
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
)
{
facei_ = -1;
DynamicList<label>& faces = cloud_.labels_;
findFaces(endPosition, faces);
if (faces.empty())
{
// endPosition is inside the cell
position_ = endPosition;
trackFraction = 1.0;
return;
}
// A face has been hit
scalar lambdaMin = GREAT;
if (faces.size() == 1)
{
lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
facei_ = faces[0];
}
else
{
forAll(faces, i)
{
scalar lam =
lambda(position_, endPosition, faces[i], stepFraction_);
if (lam < lambdaMin)
{
lambdaMin = lam;
facei_ = faces[i];
}
}
}
if (static_cast<ParticleType&>(*this).softImpact())
{
// Soft-sphere particles can travel outside the domain
// but we don't use lambda since this the particle
// is going away from face
trackFraction = 1.0;
position_ = endPosition;
}
else if (lambdaMin <= 0.0 && cloud_.internalFace(facei_))
{
// For warped faces the particle can be 'outside' the cell.
// This will yield a lambda larger than 1, or smaller than 0.
// For values < 0, the particle travels away from the cell and
// we don't move the particle (except by a small value to move
// it off the face if it is an internal face), only change
// cell.
trackFraction = Cloud<ParticleType>::minValidTrackFraction;
position_ += trackFraction*(endPosition - position_);
}
else
{
if (lambdaMin <= 1.0)
{
trackFraction = lambdaMin;
position_ += trackFraction*(endPosition - position_);
}
else
{
// For values larger than 1, we move the particle to endPosition
// only.
trackFraction = 1.0;
position_ = endPosition;
}
}
faceAction(trackFraction, endPosition, td);
// If the trackFraction = 0 something went wrong.
// Either the particle is flipping back and forth across a face perhaps
// due to velocity interpolation errors or it is in a "hole" in the mesh
// caused by face warpage.
// In both cases resolve the positional ambiguity by moving the particle
// slightly towards the cell-centre.
if (trackFraction < SMALL)
{
const polyMesh& mesh = cloud_.pMesh();
const point& cc = mesh.cellCentres()[celli_];
position_ +=
Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
}
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::faceAction
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
)
{
const polyMesh& mesh = cloud_.pMesh();
if (cloud_.internalFace(facei_))
{
// Internal face, change cell
if (celli_ == mesh.faceOwner()[facei_])
{
celli_ = mesh.faceNeighbour()[facei_];
}
else if (celli_ == mesh.faceNeighbour()[facei_])
{
celli_ = mesh.faceOwner()[facei_];
}
else
{
FatalErrorIn("Particle::faceAction")
<< "face-cell addressing failure" << nl
<< abort(FatalError);
}
}
else
{
ParticleType& p = static_cast<ParticleType&>(*this);
// Soft-sphere algorithm ignores the boundary
if (p.softImpact())
{
trackFraction = 1.0;
position_ = endPosition;
}
label patchi = patch(facei_);
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (!p.hitPatch(patch, td, patchi))
{
if (isA<wedgePolyPatch>(patch))
{
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td
);
}
else
{
p.hitPatch(patch, td);
}
}
}
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::prepareForParallelTransfer
@ -701,6 +216,7 @@ Foam::label Foam::Particle<ParticleType>::track
}
template<class ParticleType>
Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
{
@ -708,7 +224,6 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
return track(endPosition, dummyTd);
}
template<class ParticleType>
template<class TrackData>
Foam::scalar Foam::Particle<ParticleType>::trackToFace
@ -717,30 +232,170 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
TrackData& td
)
{
const polyMesh& mesh = cloud_.polyMesh_;
DynamicList<label>& faces = cloud_.labels_;
findFaces(endPosition, faces);
facei_ = -1;
scalar trackFraction = 0.0;
if (cloud_.concaveCheck_)
if (faces.empty()) // inside cell
{
if (cloud_.concaveCell()[celli_])
trackFraction = 1.0;
position_ = endPosition;
}
else // hit face
{
scalar lambdaMin = GREAT;
if (faces.size() == 1)
{
// Use a more careful tracking algorithm if the cell is concave
trackToFaceExact(trackFraction, endPosition, td);
lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
facei_ = faces[0];
}
else
{
// Use the original tracking algorithm if the cell is convex
trackToFacePlanes(trackFraction, endPosition, td);
// If the particle has to cross more than one cell to reach the
// endPosition, we check which way to go.
// If one of the faces is a boundary face and the particle is
// outside, we choose the boundary face.
// The particle is outside if one of the lambda's is > 1 or < 0
forAll(faces, i)
{
scalar lam =
lambda(position_, endPosition, faces[i], stepFraction_);
if (lam < lambdaMin)
{
lambdaMin = lam;
facei_ = faces[i];
}
}
}
bool internalFace = cloud_.internalFace(facei_);
// For warped faces the particle can be 'outside' the cell.
// This will yield a lambda larger than 1, or smaller than 0
// For values < 0, the particle travels away from the cell
// and we don't move the particle, only change cell.
// For values larger than 1, we move the particle to endPosition only.
if (lambdaMin > 0.0)
{
if (lambdaMin <= 1.0)
{
trackFraction = lambdaMin;
position_ += trackFraction*(endPosition - position_);
}
else
{
trackFraction = 1.0;
position_ = endPosition;
}
}
else if (static_cast<ParticleType&>(*this).softImpact())
{
// Soft-sphere particles can travel outside the domain
// but we don't use lambda since this the particle
// is going away from face
trackFraction = 1.0;
position_ = endPosition;
}
// change cell
if (internalFace) // Internal face
{
if (celli_ == mesh.faceOwner()[facei_])
{
celli_ = mesh.faceNeighbour()[facei_];
}
else if (celli_ == mesh.faceNeighbour()[facei_])
{
celli_ = mesh.faceOwner()[facei_];
}
else
{
FatalErrorIn
(
"Particle::trackToFace(const vector&, TrackData&)"
)<< "addressing failure" << nl
<< abort(FatalError);
}
}
else
{
ParticleType& p = static_cast<ParticleType&>(*this);
// Soft-sphere algorithm ignores the boundary
if (p.softImpact())
{
trackFraction = 1.0;
position_ = endPosition;
}
label patchi = patch(facei_);
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (!p.hitPatch(patch, td, patchi))
{
if (isA<wedgePolyPatch>(patch))
{
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td
);
}
else
{
p.hitPatch(patch, td);
}
}
}
}
else
// If the trackFraction = 0 something went wrong.
// Either the particle is flipping back and forth across a face perhaps
// due to velocity interpolation errors or it is in a "hole" in the mesh
// caused by face warpage.
// In both cases resolve the positional ambiguity by moving the particle
// slightly towards the cell-centre.
if (trackFraction < SMALL)
{
trackToFacePlanes(trackFraction, endPosition, td);
position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
}
return trackFraction;
}
template<class ParticleType>
Foam::scalar Foam::Particle<ParticleType>::trackToFace
(
@ -751,7 +406,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
return trackToFace(endPosition, dummyTd);
}
template<class ParticleType>
void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
{

46
src/lagrangian/basic/Particle/Particle.H
View File

@ -163,60 +163,22 @@ protected:
const label facei
) const;
//- Find the faces between endPosition and cell centre
//- Find the faces between position and cell centre
void findFaces
(
const vector& endPosition,
const vector& position,
DynamicList<label>& faceList
) const;
//- Find the faces between endPosition and cell centre
//- Find the faces between position and cell centre
void findFaces
(
const vector& endPosition,
const vector& position,
const label celli,
const scalar stepFraction,
DynamicList<label>& faceList
) const;
//- Determine if the testPt is inside the cell by an exact
// intersection test
bool insideCellExact
(
const vector& testPt,
const label celli,
bool beingOnAFaceMeansOutside
) const;
//- Track to cell face using face decomposition, used for
// concave cells.
template<class TrackData>
void trackToFaceExact
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
);
//- Track to cell face using the infinite planes of the faces.
// The cell *must* be convex.
template<class TrackData>
void trackToFacePlanes
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
);
//- Change cell or interact with a patch as required
template<class TrackData>
void faceAction
(
scalar& trackFraction,
const vector& endPosition,
TrackData& td
);
// Patch interactions

7
src/lagrangian/basic/Particle/ParticleI.H
View File

@ -47,7 +47,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
vector Cf = mesh.faceCentres()[facei];
// patch interaction
if (cloud_.boundaryFace(facei))
if (!cloud_.internalFace(facei))
{
label patchi = patch(facei);
const polyPatch& patch = mesh.boundaryMesh()[patchi];
@ -195,7 +195,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
vector Cf = mesh.faceCentres()[facei];
// patch interaction
if (cloud_.boundaryFace(facei))
if (!cloud_.internalFace(facei))
{
label patchi = patch(facei);
const polyPatch& patch = mesh.boundaryMesh()[patchi];
@ -207,11 +207,10 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
scalar CCf = mag((C - Cf) & Sf);
// check if distance between cell centre and face centre
// is larger than wallImpactDistance
const ParticleType& p = static_cast<const ParticleType&>(*this);
if (CCf > p.wallImpactDistance(Sf))
{
Cf -= p.wallImpactDistance(Sf)*Sf;
Cf -=p.wallImpactDistance(Sf)*Sf;
}
}
}