mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Modified version of Mattijs' concave tracking modifications. Changes
to Mattijs' work are: + Correct use of normals and cosine (dot product) to identify planar faces. + Correct use for lambda in tracking - it is a fraction, not a distance. + Not doing a reduce on demand driven construction data - not all processors call it at the same time, so crashes. This implementation contains an attempt at making the calculation of lambdaC (from the cell centre) use decomposed triangles for faces, but this is a bad approach, the concept of using lambdaC relies on the definition of a convex polyhedron, i.e. none of the planes of the faces of the polyhedron are inside the volume. Can't use this method and will need to treat a convex cell completely differently, not just some of its faces.
This commit is contained in:
graham
@ -31,6 +31,250 @@ 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-10;
|
||||
|
||||
template<class ParticleType>
|
||||
const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
|
||||
|
||||
template<class ParticleType>
|
||||
const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-6;
|
||||
|
||||
template<class ParticleType>
|
||||
const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
|
||||
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
// Calculate using face properties only
|
||||
template<class ParticleType>
|
||||
Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
|
||||
(
|
||||
const label cellI,
|
||||
const label f0I,
|
||||
const label f1I
|
||||
) const
|
||||
{
|
||||
const vectorField& fAreas = pMesh().faceAreas();
|
||||
const pointField& fCentres = pMesh().faceCentres();
|
||||
const labelList& fOwner = pMesh().faceOwner();
|
||||
const pointField& points = pMesh().points();
|
||||
|
||||
// Get f0 centre
|
||||
const point& fc0 = fCentres[f0I];
|
||||
|
||||
const face& f0 = pMesh().faces()[f0I];
|
||||
|
||||
const face& f1 = pMesh().faces()[f1I];
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
// 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];
|
||||
|
||||
// Normal distance from fc0 to pt
|
||||
scalar d = ((fc0 - pt) & fn0);
|
||||
|
||||
if (d < 0)
|
||||
{
|
||||
// Proper concave: f1 on wrong side of f0
|
||||
return 2;
|
||||
}
|
||||
}
|
||||
|
||||
// Could be a convex cell, but check angle for planar faces.
|
||||
|
||||
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
|
||||
return 1;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
void Foam::Cloud<ParticleType>::calcIntrudingFaces() const
|
||||
{
|
||||
const labelList& fOwner = pMesh().faceOwner();
|
||||
const cellList& cells = pMesh().cells();
|
||||
|
||||
intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
|
||||
PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
|
||||
intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
|
||||
PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
|
||||
|
||||
PackedBoolList cellsWithProblemFaces(cells.size());
|
||||
DynamicList<label> planarFaces;
|
||||
DynamicList<label> intrudingFaces;
|
||||
|
||||
|
||||
label nIntruded = 0;
|
||||
label nPlanar = 0;
|
||||
|
||||
forAll(cells, cellI)
|
||||
{
|
||||
const cell& cFaces = cells[cellI];
|
||||
|
||||
forAll(cFaces, i)
|
||||
{
|
||||
label f0 = cFaces[i];
|
||||
|
||||
bool own0 = (fOwner[f0] == cellI);
|
||||
|
||||
// Now check that all other faces of the cell are on the 'inside'
|
||||
// of the face.
|
||||
|
||||
bool intrudes = false;
|
||||
|
||||
forAll(cFaces, j)
|
||||
{
|
||||
if (j != i)
|
||||
{
|
||||
label state = isIntrudingFace
|
||||
(
|
||||
cellI,
|
||||
f0,
|
||||
cFaces[j]
|
||||
);
|
||||
|
||||
if (state == 1)
|
||||
{
|
||||
// planar
|
||||
planarFaces.append(f0);
|
||||
|
||||
intrudes = true;
|
||||
nPlanar++;
|
||||
break;
|
||||
}
|
||||
else if (state == 2)
|
||||
{
|
||||
// concave
|
||||
intrudingFaces.append(f0);
|
||||
|
||||
intrudes = true;
|
||||
nIntruded++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (intrudes)
|
||||
{
|
||||
cellsWithProblemFaces[cellI] = 1;
|
||||
|
||||
if (own0)
|
||||
{
|
||||
intrudesIntoOwner[f0] = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
intrudesIntoNeighbour[f0] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// intrudesIntoOwner[f0] = 1;
|
||||
// intrudesIntoNeighbour[f0] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// if (debug)
|
||||
{
|
||||
Pout<< "Cloud<ParticleType>::calcIntrudingFaces() :"
|
||||
<< " overall faces in mesh : "
|
||||
<< pMesh().nFaces() << nl
|
||||
<< " of these planar : "
|
||||
<< nPlanar << nl
|
||||
<< " of these intruding into their cell (concave) : "
|
||||
<< nIntruded << endl;
|
||||
|
||||
// Info<< "Hard coded picking up every face." << endl;
|
||||
|
||||
// Write the faces and cells that are a problem to file to be
|
||||
// compared and made into sets
|
||||
|
||||
{
|
||||
DynamicList<label> cellsWithIntrudingFaces;
|
||||
|
||||
forAll(cells, cellI)
|
||||
{
|
||||
if (cellsWithProblemFaces[cellI])
|
||||
{
|
||||
cellsWithIntrudingFaces.append(cellI);
|
||||
}
|
||||
}
|
||||
|
||||
fileName fName = pMesh().time().path()/"cellsWithProblemFaces";
|
||||
|
||||
Pout<< " Writing " << fName.name() << endl;
|
||||
|
||||
OFstream file(fName);
|
||||
|
||||
file << cellsWithIntrudingFaces;
|
||||
|
||||
file.flush();
|
||||
}
|
||||
|
||||
{
|
||||
fileName fName = pMesh().time().path()/"planarFaces";
|
||||
|
||||
Pout<< " Writing " << fName.name() << endl;
|
||||
|
||||
OFstream file(fName);
|
||||
|
||||
file << planarFaces;
|
||||
|
||||
file.flush();
|
||||
}
|
||||
|
||||
{
|
||||
fileName fName = pMesh().time().path()/"intrudingFaces";
|
||||
|
||||
Pout<< " Writing " << fName.name() << endl;
|
||||
|
||||
OFstream file(fName);
|
||||
|
||||
file << intrudingFaces;
|
||||
|
||||
file.flush();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
@ -38,12 +282,14 @@ template<class ParticleType>
|
||||
Foam::Cloud<ParticleType>::Cloud
|
||||
(
|
||||
const polyMesh& pMesh,
|
||||
const IDLList<ParticleType>& particles
|
||||
const IDLList<ParticleType>& particles,
|
||||
const bool concaveCheck
|
||||
)
|
||||
:
|
||||
cloud(pMesh),
|
||||
IDLList<ParticleType>(),
|
||||
polyMesh_(pMesh),
|
||||
concaveCheck_(concaveCheck),
|
||||
particleCount_(0)
|
||||
{
|
||||
IDLList<ParticleType>::operator=(particles);
|
||||
@ -55,20 +301,55 @@ Foam::Cloud<ParticleType>::Cloud
|
||||
(
|
||||
const polyMesh& pMesh,
|
||||
const word& cloudName,
|
||||
const IDLList<ParticleType>& particles
|
||||
const IDLList<ParticleType>& particles,
|
||||
const bool concaveCheck
|
||||
)
|
||||
:
|
||||
cloud(pMesh, cloudName),
|
||||
IDLList<ParticleType>(),
|
||||
polyMesh_(pMesh),
|
||||
concaveCheck_(concaveCheck),
|
||||
particleCount_(0)
|
||||
{
|
||||
IDLList<ParticleType>::operator=(particles);
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
void Foam::Cloud<ParticleType>::clearOut()
|
||||
{
|
||||
intrudesIntoOwnerPtr_.clear();
|
||||
intrudesIntoNeighbourPtr_.clear();
|
||||
labels_.clearStorage();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class ParticleType>
|
||||
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
|
||||
const
|
||||
{
|
||||
if (!intrudesIntoOwnerPtr_.valid())
|
||||
{
|
||||
calcIntrudingFaces();
|
||||
}
|
||||
return intrudesIntoOwnerPtr_();
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
|
||||
const
|
||||
{
|
||||
if (!intrudesIntoNeighbourPtr_.valid())
|
||||
{
|
||||
calcIntrudingFaces();
|
||||
}
|
||||
return intrudesIntoNeighbourPtr_();
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
|
||||
{
|
||||
|
||||
@ -40,6 +40,7 @@ SourceFiles
|
||||
#include "IDLList.H"
|
||||
#include "IOField.H"
|
||||
#include "polyMesh.H"
|
||||
#include "PackedBoolList.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
@ -75,6 +76,13 @@ class Cloud
|
||||
|
||||
const polyMesh& polyMesh_;
|
||||
|
||||
const bool concaveCheck_;
|
||||
|
||||
//- Does face intrude into owner cell
|
||||
mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
|
||||
//- Does face intrude into neighbour cell
|
||||
mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
|
||||
|
||||
//- Overall count of particles ever created. Never decreases.
|
||||
mutable label particleCount_;
|
||||
|
||||
@ -84,6 +92,21 @@ class Cloud
|
||||
|
||||
// Private member functions
|
||||
|
||||
//- Helper: is face f1 on wrong side of f0 (both faces are on same cell)
|
||||
// 0 : convex
|
||||
// 1 : planar to within planarCosAngle
|
||||
// 2 : concave
|
||||
label isIntrudingFace
|
||||
(
|
||||
const label cellI,
|
||||
const label f0,
|
||||
const label f1
|
||||
) const;
|
||||
|
||||
//- Analyse mesh for faces not being convex.
|
||||
// Sets intrudesIntoOwner_,intrudesIntoNeighbour_.
|
||||
void calcIntrudingFaces() const;
|
||||
|
||||
//- Initialise cloud on IO constructor
|
||||
void initCloud(const bool checkClass);
|
||||
|
||||
@ -104,13 +127,31 @@ public:
|
||||
TypeName("Cloud");
|
||||
|
||||
|
||||
// Static data members
|
||||
|
||||
//- 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;
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from mesh and a list of particles
|
||||
Cloud
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const IDLList<ParticleType>& particles
|
||||
const IDLList<ParticleType>& particles,
|
||||
const bool concaveCheck = true
|
||||
);
|
||||
|
||||
//- Construct from mesh, cloud name, and a list of particles
|
||||
@ -118,7 +159,8 @@ public:
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const word& cloudName,
|
||||
const IDLList<ParticleType>& particles
|
||||
const IDLList<ParticleType>& particles,
|
||||
const bool concaveCheck = true
|
||||
);
|
||||
|
||||
//- Construct from mesh by reading from file
|
||||
@ -126,7 +168,8 @@ public:
|
||||
Cloud
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const bool checkClass = true
|
||||
const bool checkClass = true,
|
||||
const bool concaveCheck = true
|
||||
);
|
||||
|
||||
|
||||
@ -136,10 +179,17 @@ public:
|
||||
(
|
||||
const polyMesh& pMesh,
|
||||
const word& cloudName,
|
||||
const bool checkClass = true
|
||||
const bool checkClass = true,
|
||||
const bool concaveCheck = true
|
||||
);
|
||||
|
||||
|
||||
// Destructor
|
||||
|
||||
//- Destroy any demand-driven data
|
||||
void clearOut();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Access
|
||||
@ -179,6 +229,13 @@ public:
|
||||
return IDLList<ParticleType>::size();
|
||||
};
|
||||
|
||||
//- Per face whether it intrudes into the owner cell.(demand driven)
|
||||
const PackedBoolList& intrudesIntoOwner() const;
|
||||
|
||||
//- Per face whether it intrudes into the neighbour cell. (demand
|
||||
// driven.
|
||||
const PackedBoolList& intrudesIntoNeighbour() const;
|
||||
|
||||
|
||||
// Iterators
|
||||
|
||||
|
||||
@ -62,11 +62,13 @@ template<class ParticleType>
|
||||
Foam::Cloud<ParticleType>::Cloud
|
||||
(
|
||||
const polyMesh& pMesh,
|
||||
const bool checkClass
|
||||
const bool checkClass,
|
||||
const bool concaveCheck
|
||||
)
|
||||
:
|
||||
cloud(pMesh),
|
||||
polyMesh_(pMesh),
|
||||
concaveCheck_(concaveCheck),
|
||||
particleCount_(0)
|
||||
{
|
||||
initCloud(checkClass);
|
||||
@ -78,14 +80,21 @@ Foam::Cloud<ParticleType>::Cloud
|
||||
(
|
||||
const polyMesh& pMesh,
|
||||
const word& cloudName,
|
||||
const bool checkClass
|
||||
const bool checkClass,
|
||||
const bool concaveCheck
|
||||
)
|
||||
:
|
||||
cloud(pMesh, cloudName),
|
||||
polyMesh_(pMesh),
|
||||
concaveCheck_(concaveCheck),
|
||||
particleCount_(0)
|
||||
{
|
||||
initCloud(checkClass);
|
||||
|
||||
// if (concaveCheck_)
|
||||
// {
|
||||
// intrudesIntoOwner();
|
||||
// }
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -47,10 +47,68 @@ void Foam::Particle<ParticleType>::findFaces
|
||||
const vector& C = mesh.cellCentres()[celli_];
|
||||
|
||||
faceList.clear();
|
||||
|
||||
// forAll(faces, i)
|
||||
// {
|
||||
// label facei = faces[i];
|
||||
// scalar lam = lambda(C, position, facei);
|
||||
|
||||
// if ((lam > 0) && (lam < 1.0))
|
||||
// {
|
||||
// faceList.append(facei);
|
||||
// }
|
||||
// }
|
||||
|
||||
vector deltaPosition = position - C;
|
||||
scalar deltaPositionMag = mag(deltaPosition);
|
||||
scalar lam;
|
||||
|
||||
forAll(faces, i)
|
||||
{
|
||||
label facei = faces[i];
|
||||
scalar lam = lambda(C, position, facei);
|
||||
|
||||
lam = GREAT;
|
||||
|
||||
if (cloud_.concaveCheck_)
|
||||
{
|
||||
// Check for intruding face
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
label isOwner = (celli_ == mesh.faceOwner()[facei]);
|
||||
|
||||
if
|
||||
(
|
||||
(isOwner && cloud_.intrudesIntoOwner()[facei])
|
||||
|| (!isOwner && cloud_.intrudesIntoNeighbour()[facei])
|
||||
)
|
||||
{
|
||||
// Concave face. Use exact intersection.
|
||||
|
||||
pointHit inter = mesh.faces()[facei].intersection
|
||||
(
|
||||
position_,
|
||||
deltaPosition,
|
||||
mesh.faceCentres()[facei],
|
||||
mesh.points(),
|
||||
intersection::HALF_RAY,
|
||||
Cloud<ParticleType>::intersectionTolerance
|
||||
);
|
||||
|
||||
if (inter.hit())
|
||||
{
|
||||
lam = inter.distance()/deltaPositionMag;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Convex face. Can use lambda.
|
||||
lam = lambda(C, position, facei);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
lam = lambda(C, position, facei);
|
||||
}
|
||||
|
||||
if ((lam > 0) && (lam < 1.0))
|
||||
{
|
||||
@ -241,6 +299,17 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
facei_ = -1;
|
||||
scalar trackFraction = 0.0;
|
||||
|
||||
scalar lambdaMin = GREAT;
|
||||
|
||||
// Temporary addition to report on infinite loops ->
|
||||
List<scalar> lambdas(faces.size());
|
||||
|
||||
forAll(faces, i)
|
||||
{
|
||||
lambdas[i] = lambda(position_, endPosition, faces[i], stepFraction_);
|
||||
}
|
||||
// <- Temporary addition to report on infinite loops
|
||||
|
||||
if (faces.empty()) // inside cell
|
||||
{
|
||||
trackFraction = 1.0;
|
||||
@ -248,7 +317,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
}
|
||||
else // hit face
|
||||
{
|
||||
scalar lambdaMin = GREAT;
|
||||
// scalar lambdaMin = GREAT;
|
||||
|
||||
if (faces.size() == 1)
|
||||
{
|
||||
@ -257,21 +326,86 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
}
|
||||
else
|
||||
{
|
||||
// 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.
|
||||
// // 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
|
||||
|
||||
scalar fallBackLambdaMin = GREAT;
|
||||
label fallBackFacei = -1;
|
||||
|
||||
vector deltaPosition = endPosition - position_;
|
||||
scalar deltaPositionMag = mag(deltaPosition);
|
||||
scalar hitLambda = GREAT;
|
||||
|
||||
forAll(faces, i)
|
||||
{
|
||||
scalar lam =
|
||||
lambda(position_, endPosition, faces[i], stepFraction_);
|
||||
|
||||
if (lam < lambdaMin)
|
||||
// Calculate fallback
|
||||
if (lam < fallBackLambdaMin)
|
||||
{
|
||||
lambdaMin = lam;
|
||||
facei_ = faces[i];
|
||||
fallBackLambdaMin = lam;
|
||||
fallBackFacei = faces[i];
|
||||
}
|
||||
|
||||
if (cloud_.concaveCheck_)
|
||||
{
|
||||
// Check for intruding face
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
label isOwner = (celli_ == mesh.faceOwner()[faces[i]]);
|
||||
|
||||
if
|
||||
(
|
||||
(isOwner && cloud_.intrudesIntoOwner()[faces[i]])
|
||||
|| (!isOwner && cloud_.intrudesIntoNeighbour()[faces[i]])
|
||||
)
|
||||
{
|
||||
// Concave face. Use exact intersection.
|
||||
|
||||
pointHit inter = mesh.faces()[faces[i]].intersection
|
||||
(
|
||||
position_,
|
||||
deltaPosition,
|
||||
mesh.faceCentres()[faces[i]],
|
||||
mesh.points(),
|
||||
intersection::HALF_RAY,
|
||||
Cloud<ParticleType>::intersectionTolerance
|
||||
);
|
||||
|
||||
if (inter.hit())
|
||||
{
|
||||
hitLambda = inter.distance()/deltaPositionMag;
|
||||
|
||||
if (hitLambda < lambdaMin)
|
||||
{
|
||||
lambdaMin = hitLambda;
|
||||
facei_ = faces[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Convex face. Can use lambda.
|
||||
|
||||
if (lam < lambdaMin)
|
||||
{
|
||||
lambdaMin = lam;
|
||||
facei_ = faces[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Fall back to the original algorithm if no
|
||||
// intersection was found
|
||||
if (facei_ == -1)
|
||||
{
|
||||
lambdaMin = fallBackLambdaMin;
|
||||
facei_ = fallBackFacei;
|
||||
}
|
||||
}
|
||||
|
||||
@ -320,9 +454,8 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
FatalErrorIn
|
||||
(
|
||||
"Particle::trackToFace(const vector&, TrackData&)"
|
||||
)
|
||||
<< "addressing failure" << nl
|
||||
<< abort(FatalError);
|
||||
)<< "addressing failure" << nl
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
else
|
||||
@ -381,7 +514,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
p.hitPatch(patch, td);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@ -391,14 +523,44 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
// caused by face warpage.
|
||||
// In both cases resolve the positional ambiguity by moving the particle
|
||||
// slightly towards the cell-centre.
|
||||
if (trackFraction < SMALL)
|
||||
|
||||
if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
|
||||
{
|
||||
position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
|
||||
Pout<< "tracking error "
|
||||
<< origId_ << " "
|
||||
<< origProc_ << " "
|
||||
<< position_ << " "
|
||||
<< endPosition << " "
|
||||
<< position_ + 1000*(endPosition - position_) << " "
|
||||
<< position_ - 1000*(endPosition - position_) << " "
|
||||
<< lambdaMin << " "
|
||||
<< trackFraction << " "
|
||||
<< stepFraction_ << " "
|
||||
<< celli_ << " "
|
||||
<< facei_ << " "
|
||||
<< cloud_.intrudesIntoOwner()[facei_] << " "
|
||||
<< cloud_.intrudesIntoNeighbour()[facei_] << " "
|
||||
<< faces << " "
|
||||
<< lambdas
|
||||
<< endl;
|
||||
|
||||
// if (stepFraction_ < SMALL)
|
||||
// {
|
||||
// Pout<< "Intervene to mark face" << endl;
|
||||
|
||||
// cloud_.intrudesIntoOwnerPtr_()[facei_] = 1;
|
||||
// cloud_.intrudesIntoNeighbourPtr_()[facei_] = 1;
|
||||
// }
|
||||
|
||||
const point& cc = mesh.cellCentres()[celli_];
|
||||
position_ +=
|
||||
Cloud<ParticleType>::trackingRescueTolerance*(cc - position_);
|
||||
}
|
||||
|
||||
return trackFraction;
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
(
|
||||
@ -409,6 +571,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
|
||||
return trackToFace(endPosition, dummyTd);
|
||||
}
|
||||
|
||||
|
||||
template<class ParticleType>
|
||||
void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
|
||||
{
|
||||
|
||||
Reference in New Issue
Block a user