diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
index 8166029dbd4335426cbf0b4605a6180c383e174f..02d87225606b2b92ebbe4b19bdbc0c5f047ebeae 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -35,6 +35,57 @@ const std::size_t Foam::wallBoundedParticle::sizeofFields_
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+Foam::tetIndices Foam::wallBoundedParticle::currentTetIndices() const
+{
+ // Replacement for particle::currentTetIndices that avoids error
+ // upon invalid tetBasePtIs
+
+ const faceList& pFaces = mesh_.faces();
+ const labelList& pOwner = mesh_.faceOwner();
+
+ const Foam::face& f = pFaces[tetFacei_];
+
+ bool own = (pOwner[tetFacei_] == celli_);
+
+ label faceBasePtI = mesh_.tetBasePtIs()[tetFacei_];
+ if (faceBasePtI == -1)
+ {
+ //WarningInFunction
+ // << "No base point for face " << tetFacei_ << ", " << f
+ // << ", produces a decomposition that has a minimum "
+ // << "volume greater than tolerance."
+ // << endl;
+ faceBasePtI = 0;
+ }
+
+ label facePtI = (tetPti_ + faceBasePtI) % f.size();
+ label otherFacePtI = f.fcIndex(facePtI);
+
+ label facePtAI;
+ label facePtBI;
+ if (own)
+ {
+ facePtAI = facePtI;
+ facePtBI = otherFacePtI;
+ }
+ else
+ {
+ facePtAI = otherFacePtI;
+ facePtBI = facePtI;
+ }
+
+ return tetIndices
+ (
+ celli_,
+ tetFacei_,
+ faceBasePtI,
+ facePtAI,
+ facePtBI,
+ tetPti_
+ );
+}
+
+
Foam::edge Foam::wallBoundedParticle::currentEdge() const
{
if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
@@ -57,6 +108,18 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
else
{
label faceBasePti = mesh_.tetBasePtIs()[tetFace()];
+ if (faceBasePti == -1)
+ {
+ //FatalErrorInFunction
+ //WarningInFunction
+ // << "No base point for face " << tetFace() << ", " << f
+ // << ", produces a decomposition that has a minimum "
+ // << "volume greater than tolerance."
+ // //<< abort(FatalError);
+ // << endl;
+ faceBasePti = 0;
+ }
+
label diagPti = (faceBasePti+diagEdge_)%f.size();
return edge(f[faceBasePti], f[diagPti]);
@@ -66,11 +129,118 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
void Foam::wallBoundedParticle::crossEdgeConnectedFace
(
- const edge& meshEdge
+ const label& celli,
+ label& tetFacei,
+ label& tetPti,
+ const edge& e
)
+{
+ const faceList& pFaces = mesh_.faces();
+ const cellList& pCells = mesh_.cells();
+
+ const Foam::face& f = pFaces[tetFacei];
+
+ const Foam::cell& thisCell = pCells[celli];
+
+ forAll(thisCell, cFI)
+ {
+ // Loop over all other faces of this cell and
+ // find the one that shares this edge
+
+ label fI = thisCell[cFI];
+
+ if (tetFacei == fI)
+ {
+ continue;
+ }
+
+ const Foam::face& otherFace = pFaces[fI];
+
+ label edDir = otherFace.edgeDirection(e);
+
+ if (edDir == 0)
+ {
+ continue;
+ }
+ else if (f == pFaces[fI])
+ {
+ // This is a necessary condition if using duplicate baffles
+ // (so coincident faces). We need to make sure we don't cross into
+ // the face with the same vertices since we might enter a tracking
+ // loop where it never exits. This test should be cheap
+ // for most meshes so can be left in for 'normal' meshes.
+ continue;
+ }
+ else
+ {
+ //Found edge on other face
+ tetFacei = fI;
+
+ label eIndex = -1;
+
+ if (edDir == 1)
+ {
+ // Edge is in the forward circulation of this face, so
+ // work with the start point of the edge
+ eIndex = findIndex(otherFace, e.start());
+ }
+ else
+ {
+ // edDir == -1, so the edge is in the reverse
+ // circulation of this face, so work with the end
+ // point of the edge
+ eIndex = findIndex(otherFace, e.end());
+ }
+
+ label tetBasePtI = mesh_.tetBasePtIs()[fI];
+
+ if (tetBasePtI == -1)
+ {
+ //FatalErrorInFunction
+ //WarningInFunction
+ // << "No base point for face " << fI << ", " << f
+ // << ", produces a decomposition that has a minimum "
+ // << "volume greater than tolerance."
+ // //<< abort(FatalError);
+ // << endl;
+ tetBasePtI = 0;
+ }
+
+ // Find eIndex relative to the base point on new face
+ eIndex -= tetBasePtI;
+
+ if (neg(eIndex))
+ {
+ eIndex = (eIndex + otherFace.size()) % otherFace.size();
+ }
+
+ if (eIndex == 0)
+ {
+ // The point is the base point, so this is first tet
+ // in the face circulation
+ tetPti = 1;
+ }
+ else if (eIndex == otherFace.size() - 1)
+ {
+ // The point is the last before the base point, so
+ // this is the last tet in the face circulation
+ tetPti = otherFace.size() - 2;
+ }
+ else
+ {
+ tetPti = eIndex;
+ }
+
+ break;
+ }
+ }
+}
+
+
+void Foam::wallBoundedParticle::crossEdgeConnectedFace(const edge& meshEdge)
{
// Update tetFace, tetPt
- particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
+ crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
// Update face to be same as tracking one
face() = tetFace();
@@ -164,14 +334,13 @@ void Foam::wallBoundedParticle::crossDiagonalEdge()
Foam::scalar Foam::wallBoundedParticle::trackFaceTri
(
+ const vector& n,
const vector& endPosition,
label& minEdgei
)
{
// Track p from position to endPosition
const triFace tri(currentTetIndices().faceTriIs(mesh_));
- vector n = tri.normal(mesh_.points());
- n /= mag(n)+VSMALL;
// Check which edge intersects the trajectory.
// Project trajectory onto triangle
@@ -245,6 +414,7 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
bool Foam::wallBoundedParticle::isTriAlongTrack
(
+ const vector& n,
const point& endPosition
) const
{
@@ -267,10 +437,6 @@ bool Foam::wallBoundedParticle::isTriAlongTrack
const vector dir = endPosition-position();
- // Get normal of currentE
- vector n = triVerts.normal(mesh_.points());
- n /= mag(n);
-
forAll(triVerts, i)
{
label j = triVerts.fcIndex(i);
@@ -382,7 +548,7 @@ Foam::Ostream& Foam::operator<<
{
const wallBoundedParticle& p = ip.t_;
- tetPointRef tpr(p.currentTet());
+ tetPointRef tpr(p.currentTetIndices().tet(p.mesh()));
os << " " << static_cast<const particle&>(p) << nl
<< " tet:" << nl;
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
index a82b3f5b64cd0f50e10dd4a3e7849075ae831a9b..654ee8ae5755e4d2400bc8592b6a80ada4178ca7 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -123,6 +123,20 @@ protected:
//- Construct current edge
edge currentEdge() const;
+ //- Replacement for particle::currentTetIndices() that avoids bombing
+ // out on invalid tet decomposition (tetBasePtIs = -1)
+ tetIndices currentTetIndices() const;
+
+ //- Replacement for particle::crossEdgeConnectedFace that avoids bombing
+ // out on invalid tet decomposition (tetBasePtIs = -1)
+ void crossEdgeConnectedFace
+ (
+ const label& celli,
+ label& tetFacei,
+ label& tetPti,
+ const edge& e
+ );
+
//- Cross mesh edge into different face on same cell
void crossEdgeConnectedFace(const edge& meshEdge);
@@ -130,10 +144,10 @@ protected:
void crossDiagonalEdge();
//- Track through single triangle
- scalar trackFaceTri(const vector& endPosition, label& minEdgei);
+ scalar trackFaceTri(const vector& n, const vector& endPosition, label&);
//- Is current triangle in the track direction
- bool isTriAlongTrack(const point& endPosition) const;
+ bool isTriAlongTrack(const vector& n, const point& endPosition) const;
// Patch interactions
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
index e4f151fe26d9169aa57e50318231a4ebe7b26534..731c4ad0b3971a58eabda3ba0e0c943e6d64d26a 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -175,7 +175,14 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// endposition? i.e. since volume of nbr tet is positive the
// tracking direction should be into the tet.
tetIndices nbrTi(nbrCelli, tetFacei_, tetPti_, mesh_);
- if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
+
+ bool posVol = (nbrTi.tet(mesh_).mag() > 0);
+
+ if
+ (
+ posVol
+ == ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
+ )
{
// Change into nbrCell. No need to change tetFace, tetPt.
//Pout<< " crossed from cell:" << celli_
@@ -213,12 +220,22 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
<< abort(FatalError);
}
+ const triFace tri(currentTetIndices().faceTriIs(mesh_));
+ vector n = tri.normal(mesh_.points());
+ n /= mag(n);
+
point projectedEndPosition = endPosition;
+
+ const bool posVol = (currentTetIndices().tet(mesh_).mag() > 0);
+
+ if (!posVol)
+ {
+ // Negative tet volume. Track back by setting the end point
+ projectedEndPosition = position() - (endPosition-position());
+ }
+
// Remove normal component
{
- const triFace tri(currentTetIndices().faceTriIs(mesh_));
- vector n = tri.normal(mesh_.points());
- n /= mag(n);
const point& basePt = mesh_.points()[tri[0]];
projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
}
@@ -234,7 +251,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
{
// See if the current triangle has got a point on the
// correct side of the edge.
- doTrack = isTriAlongTrack(projectedEndPosition);
+ doTrack = isTriAlongTrack(n, projectedEndPosition);
}
@@ -242,7 +259,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
{
// Track across triangle. Return triangle edge crossed.
label triEdgei = -1;
- trackFraction = trackFaceTri(projectedEndPosition, triEdgei);
+ trackFraction = trackFaceTri(n, projectedEndPosition, triEdgei);
if (triEdgei == -1)
{
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
index d626c63fcb9c6d94594e597456b7ee011cef6e7a..ef914b5b3209d4d6b8cf89be194746db2e8bc54e 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -120,11 +120,12 @@ Foam::vector Foam::wallBoundedStreamLineParticle::sample
{
// Stagnant particle. Might as well stop
lifeTime_ = 0;
+ return vector::zero;
+ }
+ else
+ {
+ return U/magU;
}
-
- U /= magU;
-
- return U;
}
@@ -267,16 +268,6 @@ bool Foam::wallBoundedStreamLineParticle::move
// Force removal
lifeTime_ = 0;
}
-
- if
- (
- !td.keepParticle
- || td.switchProcessor
- || lifeTime_ == 0
- )
- {
- break;
- }
}