Adding instrumentation Pout statements.

Using an additional test for best guess face, which is the nearest
face, although this is not fully thought out, and can be a problem if
the nearest face is on a patch and this is not appropriate.
This commit is contained in:
graham
2009-10-08 18:50:46 +01:00
parent 2d8e739837
commit 2d6c691e75
2 changed files with 127 additions and 38 deletions

View File

@ -39,10 +39,10 @@ template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-10; const Foam::scalar Foam::Cloud<ParticleType>::minValidTrackFraction = 1e-10;
template<class ParticleType> template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3; const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-6;
template<class ParticleType> template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = SMALL; const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-3;
template<class ParticleType> template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6); const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
@ -189,14 +189,15 @@ void Foam::Cloud<ParticleType>::calcConcaveCells() const
intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces())); intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces())); intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
forAll(cells, cellI) forAll(cells, cellI)
{ {
if (isConcaveCell(cellI)) // if (isConcaveCell(cellI))
{ // {
// concaveCell[cellI] = 1;
// }
concaveCell[cellI] = 1; concaveCell[cellI] = 1;
} }
}
{ {
// Write cells that are a problem to file // Write cells that are a problem to file

View File

@ -118,6 +118,14 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
return; return;
} }
Pout<< nl << origProc_ << " "
<< origId_ << " "
<< position_ << " "
<< endPosition << " "
<< celli_ << nl
<< potentiallyCrossedFaces
<< endl;
label nFaceCrossings = 0; label nFaceCrossings = 0;
forAll (potentiallyCrossedFaces, pCFI) forAll (potentiallyCrossedFaces, pCFI)
@ -156,11 +164,13 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
if (nFaceCrossings > 1) if (nFaceCrossings > 1)
{ {
Pout<< "In cell " << celli_ << " there were " << nFaceCrossings Pout<< "In cell " << celli_ << " there were " << nFaceCrossings
<< " face crossings detected tracking from cell centre to " << " face crossings detected tracking from concave cell centre to "
<< " endPosition" << " endPosition"
<< endl; << endl;
} }
Pout<< nFaceCrossings << " " << (nFaceCrossings % 2) << endl;
if (nFaceCrossings % 2 == 0) if (nFaceCrossings % 2 == 0)
{ {
// Even number of face crossings, so the particle must end up // Even number of face crossings, so the particle must end up
@ -201,7 +211,7 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
// has been found // has been found
Pout<< "New face crossing " << facei Pout<< "New face crossing " << facei
<< " of cell " << celli_ << " found" << " of concave cell " << celli_ << " found"
<< endl; << endl;
potentiallyCrossedFaces.append(facei); potentiallyCrossedFaces.append(facei);
@ -209,6 +219,8 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
} }
} }
Pout<< potentiallyCrossedFaces << endl;
vector deltaPosition = endPosition - position_; vector deltaPosition = endPosition - position_;
scalar tmpLambda = GREAT; scalar tmpLambda = GREAT;
@ -261,9 +273,19 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
{ {
// No face has been identified to be crossed yet, but the cell // No face has been identified to be crossed yet, but the cell
// must have been left, so the best guess of which face to // must have been left, so the best guess of which face to
// cross is required. This face will not be one of those // cross is required.
// faces already tested, and could be any other face of the
// cell. vector delta =
mag(mesh.cellCentres()[celli_] - position_)
*deltaPosition/mag(deltaPosition);
Pout<< "Need a best guess " << nl
<< position_ << nl
<< position_ + delta << nl
<< position_ - delta
<< endl;
// ---- Best guess by projection --
DynamicList<scalar> tmpLambdas; DynamicList<scalar> tmpLambdas;
DynamicList<label> tmpLambdaFaceIs; DynamicList<label> tmpLambdaFaceIs;
@ -272,10 +294,6 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
{ {
label facei = faces[i]; label facei = faces[i];
if (findIndex(potentiallyCrossedFaces, facei) == -1)
{
// This face has not been assessed yet.
// Use exact FULL_RAY intersection to allow negative // Use exact FULL_RAY intersection to allow negative
// values // values
@ -297,8 +315,37 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
tmpLambdas.append(inter.distance()); tmpLambdas.append(inter.distance());
tmpLambdaFaceIs.append(facei); tmpLambdaFaceIs.append(facei);
} }
inter = mesh.faces()[facei].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[facei],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
Pout<< "Test repeat HALF "
<< facei << " " << inter.distance() << endl;
inter = mesh.faces()[facei].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[facei],
mesh.points(),
intersection::FULL_RAY,
Cloud<ParticleType>::intersectionTolerance
);
Pout<< "Test repeat FULL "
<< facei << " " << inter.distance() << nl
<< (inter.distance()*deltaPosition) + position_
<< endl;
} }
}
Pout<< tmpLambdaFaceIs << " " << tmpLambdas << endl;
// The closest lambda value (less than zero) or (greater than or // The closest lambda value (less than zero) or (greater than or
// equal to one), and corresponding face. // equal to one), and corresponding face.
@ -343,12 +390,48 @@ void Foam::Particle<ParticleType>::trackToFaceConcave
} }
else else
{ {
Pout<< "Tracking failure in concave cell - didn't find a " Pout<< "Tracking in concave cell - didn't find a "
<< "best guess face" << endl; << "best guess face" << nl
<< origId_ << " "
<< origProc_ << " "
<< position_ << " "
<< endPosition << " "
<< stepFraction_ << " "
<< celli_ << " "
<< endl;
// ---- Best guess by nearest ----
scalar minDistance = GREAT;
forAll(faces, i)
{
label facei = faces[i];
if (cloud_.internalFace(facei))
{
pointHit nearest = mesh.faces()[facei].nearestPoint
(
position_,
mesh.points()
);
if (nearest.distance() < minDistance)
{
minDistance = nearest.distance();
facei_ = facei;
}
} }
} }
// ---- End of best guess by nearest ----
}
}
Pout<< facei_ << " " << celli_ << endl;
faceAction(trackFraction, endPosition, td); faceAction(trackFraction, endPosition, td);
Pout<< facei_ << " " << celli_ << endl;
} }
@ -722,6 +805,11 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
// Use a more careful tracking algorithm if the cell is concave // Use a more careful tracking algorithm if the cell is concave
trackToFaceConcave(trackFraction, endPosition, td); trackToFaceConcave(trackFraction, endPosition, td);
} }
else
{
// Use a more careful tracking algorithm if the cell is concave
trackToFaceConvex(trackFraction, endPosition, td);
}
} }
else else
{ {