diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
new file mode 100644
index 0000000000000000000000000000000000000000..e73f4e72f4884857614ab45c4aca569ce215a760
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
@@ -0,0 +1,659 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*----------------------------------------------------------------------------*/
+
+#include "wallBoundedParticle.H"
+#include "vectorFieldIOField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+// defineParticleTypeNameAndDebug(wallBoundedParticle, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+//// Check position is inside tet
+//void Foam::wallBoundedParticle::checkInside() const
+//{
+// const tetIndices ti(currentTetIndices());
+// const tetPointRef tpr(ti.tet(mesh_));
+// if (!tpr.inside(position()))
+// {
+// FatalErrorIn("wallBoundedParticle::checkInside(..)")
+// << "Particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "is not inside " << tpr
+// << abort(FatalError);
+// }
+//}
+//
+//
+//void Foam::wallBoundedParticle::checkOnEdge() const
+//{
+// // Check that edge (as indicated by meshEdgeStart_, diagEdge_) is
+// // indeed one that contains the position.
+// const edge e = currentEdge();
+//
+// linePointRef ln(e.line(mesh_.points()));
+//
+// pointHit ph(ln.nearestDist(position()));
+//
+// if (ph.distance() > 1E-6)
+// {
+// FatalErrorIn
+// (
+// "wallBoundedParticle::checkOnEdge()"
+// ) << "Problem :"
+// << " particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "edge:" << e
+// << " at:" << ln
+// << " distance:" << ph.distance()
+// << abort(FatalError);
+// }
+//}
+//
+//
+//void Foam::wallBoundedParticle::checkOnTriangle(const point& p)
+//const
+//{
+// const triFace tri(currentTetIndices().faceTriIs(mesh_));
+// pointHit ph = tri.nearestPoint(p, mesh_.points());
+// if (ph.distance() > 1e-9)
+// {
+// FatalErrorIn
+// (
+// "wallBoundedParticle::checkOnTriangle(const point&)"
+// ) << "Problem :"
+// << " particle:" //<< static_cast<const particle&>(*this)
+// << info()
+// << "point:" << p
+// << " distance:" << ph.distance()
+// << abort(FatalError);
+// }
+//}
+
+
+// Construct the edge the particle is on (according to meshEdgeStart_,
+// diagEdge_)
+Foam::edge Foam::wallBoundedParticle::currentEdge() const
+{
+ if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
+ {
+ FatalErrorIn("wallBoundedParticle::currentEdge() const")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "cannot both be on a mesh edge and a face-diagonal edge."
+ << " meshEdgeStart_:" << meshEdgeStart_
+ << " diagEdge_:" << diagEdge_
+ << abort(FatalError);
+ }
+
+ const Foam::face& f = mesh_.faces()[tetFace()];
+
+ if (meshEdgeStart_ != -1)
+ {
+ return edge(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ }
+ else
+ {
+ label faceBasePtI = mesh_.tetBasePtIs()[tetFace()];
+ label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+ return edge(f[faceBasePtI], f[diagPtI]);
+ }
+}
+
+
+void Foam::wallBoundedParticle::crossEdgeConnectedFace
+(
+ const edge& meshEdge
+)
+{
+ //label oldFaceI = tetFace();
+
+ // Update tetFace, tetPt
+ particle::crossEdgeConnectedFace(cell(), tetFace(), tetPt(), meshEdge);
+
+ // Update face to be same as tracking one
+ face() = tetFace();
+
+
+ // And adapt meshEdgeStart_.
+ const Foam::face& f = mesh_.faces()[tetFace()];
+ label fp = findIndex(f, meshEdge[0]);
+
+ if (f.nextLabel(fp) == meshEdge[1])
+ {
+ meshEdgeStart_ = fp;
+ }
+ else
+ {
+ label fpMin1 = f.rcIndex(fp);
+
+ if (f[fpMin1] == meshEdge[1])
+ {
+ meshEdgeStart_ = fpMin1;
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "wallBoundedParticle::crossEdgeConnectedFace"
+ "(const edge&)"
+ ) << "Problem :"
+ << " particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "face:" << tetFace()
+ << " verts:" << f
+ << " meshEdge:" << meshEdge
+ << abort(FatalError);
+ }
+ }
+
+ diagEdge_ = -1;
+
+ //Pout<< " crossed meshEdge "
+ // << meshEdge.line(mesh().points())
+ // << " from face:" << oldFaceI
+ // << " to face:" << tetFace() << endl;
+
+
+ // Check that still on same mesh edge
+
+ const edge eNew(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ if (eNew != meshEdge)
+ {
+ FatalErrorIn
+ (
+ "wallBoundedParticle::crossEdgeConnectedFace"
+ "(const edge&)"
+ ) << "Problem" << abort(FatalError);
+ }
+
+
+ // Check that edge (as indicated by meshEdgeStart_) is indeed one that
+ // contains the position.
+ //checkOnEdge();
+}
+
+
+void Foam::wallBoundedParticle::crossDiagonalEdge()
+{
+ if (diagEdge_ == -1)
+ {
+ FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "not on a diagonal edge" << abort(FatalError);
+ }
+ if (meshEdgeStart_ != -1)
+ {
+ FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
+ }
+
+ //label oldTetPt = tetPt();
+
+ const Foam::face& f = mesh_.faces()[tetFace()];
+
+ // tetPtI starts from 1, goes up to f.size()-2
+
+ if (tetPt() == diagEdge_)
+ {
+ tetPt() = f.rcIndex(tetPt());
+ }
+ else
+ {
+ label nextTetPt = f.fcIndex(tetPt());
+ if (diagEdge_ == nextTetPt)
+ {
+ tetPt() = nextTetPt;
+ }
+ else
+ {
+ FatalErrorIn("wallBoundedParticle::crossDiagonalEdge()")
+ << "Particle:" //<< static_cast<const particle&>(*this)
+ << info()
+ << "tetPt:" << tetPt()
+ << " diagEdge:" << diagEdge_ << abort(FatalError);
+ }
+ }
+
+ meshEdgeStart_ = -1;
+
+ //Pout<< " crossed diagonalEdge "
+ // << currentEdge().line(mesh().points())
+ // << " from tetPt:" << oldTetPt
+ // << " to tetPt:" << tetPt() << endl;
+}
+
+
+//- Track through a single triangle.
+// Gets passed tet+triangle the particle is in. Updates position() but nothing
+// else. Returns the triangle edge the particle is now on.
+Foam::scalar Foam::wallBoundedParticle::trackFaceTri
+(
+ const vector& endPosition,
+ label& minEdgeI
+)
+{
+ // Track p from position to endPosition
+ const triFace tri(currentTetIndices().faceTriIs(mesh_));
+ vector n = tri.normal(mesh_.points());
+ //if (mag(n) < sqr(SMALL))
+ //{
+ // FatalErrorIn("wallBoundedParticle::trackFaceTri(..)")
+ // << "Small triangle." //<< static_cast<const particle&>(*this)
+ // << info()
+ // << "n:" << n
+ // << abort(FatalError);
+ //}
+ n /= mag(n)+VSMALL;
+
+ // Check which edge intersects the trajectory.
+ // Project trajectory onto triangle
+ minEdgeI = -1;
+ scalar minS = 1; // end position
+
+ //const point oldPosition(position());
+
+
+ edge currentE(-1, -1);
+ if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ {
+ currentE = currentEdge();
+ }
+
+ // Determine path along line position+s*d to see where intersections
+ // are.
+
+ forAll(tri, i)
+ {
+ label j = tri.fcIndex(i);
+
+ const point& pt0 = mesh_.points()[tri[i]];
+ const point& pt1 = mesh_.points()[tri[j]];
+
+ if (edge(tri[i], tri[j]) == currentE)
+ {
+ // Do not check particle is on
+ continue;
+ }
+
+ // Outwards pointing normal
+ vector edgeNormal = (pt1-pt0)^n;
+
+ //if (mag(edgeNormal) < SMALL)
+ //{
+ // FatalErrorIn("wallBoundedParticle::trackFaceTri(..)")
+ // << "Edge not perpendicular to triangle."
+ // //<< static_cast<const particle&>(*this)
+ // << info()
+ // << "triangle n:" << n
+ // << " edgeNormal:" << edgeNormal
+ // << " on tri:" << tri
+ // << " at:" << pt0
+ // << " at:" << pt1
+ // << abort(FatalError);
+ //}
+
+
+ edgeNormal /= mag(edgeNormal)+VSMALL;
+
+ // Determine whether position and end point on either side of edge.
+ scalar sEnd = (endPosition-pt0)&edgeNormal;
+ if (sEnd >= 0)
+ {
+ // endPos is outside triangle. position() should always be
+ // inside.
+ scalar sStart = (position()-pt0)&edgeNormal;
+ if (mag(sEnd-sStart) > VSMALL)
+ {
+ scalar s = sStart/(sStart-sEnd);
+
+ if (s >= 0 && s < minS)
+ {
+ minS = s;
+ minEdgeI = i;
+ }
+ }
+ }
+ }
+
+ if (minEdgeI != -1)
+ {
+ position() += minS*(endPosition-position());
+ }
+ else
+ {
+ // Did not hit any edge so tracked to the end position. Set position
+ // without any calculation to avoid truncation errors.
+ position() = endPosition;
+ minS = 1.0;
+ }
+
+ // Project position() back onto plane of triangle
+ const point& triPt = mesh_.points()[tri[0]];
+ position() -= ((position()-triPt)&n)*n;
+
+
+ //Pout<< " tracked from:" << oldPosition << " to:" << position()
+ // << " projectedEnd:" << endPosition
+ // << " at s:" << minS << endl;
+ //if (minEdgeI != -1)
+ //{
+ // Pout<< " on edge:" << minEdgeI
+ // << " on edge:"
+ // << mesh_.points()[tri[minEdgeI]]
+ // << mesh_.points()[tri[tri.fcIndex(minEdgeI)]]
+ // << endl;
+ //}
+
+ return minS;
+}
+
+
+// See if the current triangle has got a point on the
+// correct side of the edge.
+bool Foam::wallBoundedParticle::isTriAlongTrack
+(
+ const point& endPosition
+) const
+{
+ const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
+ const edge currentE = currentEdge();
+
+ //if (debug)
+ {
+ if
+ (
+ currentE[0] == currentE[1]
+ || findIndex(triVerts, currentE[0]) == -1
+ || findIndex(triVerts, currentE[1]) == -1
+ )
+ {
+ FatalErrorIn
+ (
+ "wallBoundedParticle::isTriAlongTrack"
+ "(const point&)"
+ ) << "Edge " << currentE << " not on triangle " << triVerts
+ << info()
+ << abort(FatalError);
+ }
+ }
+
+
+ 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);
+ const point& pt0 = mesh_.points()[triVerts[i]];
+ const point& pt1 = mesh_.points()[triVerts[j]];
+
+ if (edge(triVerts[i], triVerts[j]) == currentE)
+ {
+ vector edgeNormal = (pt1-pt0)^n;
+ return (dir&edgeNormal) < 0;
+ }
+ }
+
+ FatalErrorIn
+ (
+ "wallBoundedParticle::isTriAlongTrack"
+ "(const point&)"
+ ) << "Problem" << abort(FatalError);
+
+ return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+ const polyMesh& mesh,
+ const vector& position,
+ const label cellI,
+ const label tetFaceI,
+ const label tetPtI,
+ const label meshEdgeStart,
+ const label diagEdge
+)
+:
+ particle(mesh, position, cellI, tetFaceI, tetPtI),
+ meshEdgeStart_(meshEdgeStart),
+ diagEdge_(diagEdge)
+{
+ //checkInside();
+
+ //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ //{
+ // checkOnEdge();
+ //}
+
+ // Unfortunately have no access to trackdata so cannot check if particle
+ // is on a wallPatch or has an mesh edge set (either of which is
+ // a requirement).
+}
+
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+ const polyMesh& mesh,
+ Istream& is,
+ bool readFields
+)
+:
+ particle(mesh, is, readFields)
+{
+ if (readFields)
+ {
+ is >> meshEdgeStart_ >> diagEdge_;
+ }
+
+ // Check state of Istream
+ is.check
+ (
+ "wallBoundedParticle::wallBoundedParticle"
+ "(const Cloud<wallBoundedParticle>&, Istream&, bool)"
+ );
+}
+
+
+Foam::wallBoundedParticle::wallBoundedParticle
+(
+ const wallBoundedParticle& p
+)
+:
+ particle(p),
+ meshEdgeStart_(p.meshEdgeStart_),
+ diagEdge_(p.diagEdge_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+//void Foam::wallBoundedParticle::readFields
+//(
+// Cloud<wallBoundedParticle>& c
+//)
+//{
+// if (!c.size())
+// {
+// return;
+// }
+//
+// particle::readFields(c);
+//
+// IOField<label> meshEdgeStart
+// (
+// c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
+// );
+//
+// IOField<label> diagEdge
+// (
+// c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
+// );
+// c.checkFieldIOobject(c, diagEdge);
+//
+// label i = 0;
+// forAllIter(Cloud<wallBoundedParticle>, c, iter)
+// {
+// iter().meshEdgeStart_ = meshEdgeStart[i];
+// iter().diagEdge_ = diagEdge[i];
+// i++;
+// }
+//}
+//
+//
+//void Foam::wallBoundedParticle::writeFields
+//(
+// const Cloud<wallBoundedParticle>& c
+//)
+//{
+// particle::writeFields(c);
+//
+// label np = c.size();
+//
+// IOField<label> meshEdgeStart
+// (
+// c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
+// np
+// );
+// IOField<label> diagEdge
+// (
+// c.fieldIOobject("diagEdge", IOobject::NO_READ),
+// np
+// );
+//
+// label i = 0;
+// forAllConstIter(Cloud<wallBoundedParticle>, c, iter)
+// {
+// meshEdgeStart[i] = iter().meshEdgeStart_;
+// diagEdge[i] = iter().diagEdge_;
+// i++;
+// }
+//
+// meshEdgeStart.write();
+// diagEdge.write();
+//}
+
+void Foam::wallBoundedParticle::write(Ostream& os, bool writeFields) const
+{
+ const particle& p = static_cast<const particle&>(*this);
+
+ if (os.format() == IOstream::ASCII)
+ {
+ // Write base particle
+ p.write(os, writeFields);
+
+ if (writeFields)
+ {
+ // Write the additional entries
+ os << token::SPACE << meshEdgeStart_
+ << token::SPACE << diagEdge_;
+ }
+ }
+ else
+ {
+ // Write base particle
+ p.write(os, writeFields);
+
+ // Write additional entries
+ if (writeFields)
+ {
+ os.write
+ (
+ reinterpret_cast<const char*>(&meshEdgeStart_),
+ sizeof(meshEdgeStart_)
+ + sizeof(diagEdge_)
+ );
+ }
+ }
+
+ // Check state of Ostream
+ os.check("wallBoundedParticle::write(Ostream& os, bool) const");
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const wallBoundedParticle& p
+)
+{
+ // Write all data
+ p.write(os, true);
+
+ return os;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+ Ostream& os,
+ const InfoProxy<wallBoundedParticle>& ip
+)
+{
+ const wallBoundedParticle& p = ip.t_;
+
+ tetPointRef tpr(p.currentTet());
+
+ os << " " << static_cast<const particle&>(p) << nl
+ << " tet:" << nl;
+ os << " ";
+ meshTools::writeOBJ(os, tpr.a());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.b());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.c());
+ os << " ";
+ meshTools::writeOBJ(os, tpr.d());
+ os << " l 1 2" << nl
+ << " l 1 3" << nl
+ << " l 1 4" << nl
+ << " l 2 3" << nl
+ << " l 2 4" << nl
+ << " l 3 4" << nl;
+ os << " ";
+ meshTools::writeOBJ(os, p.position());
+
+ return os;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
new file mode 100644
index 0000000000000000000000000000000000000000..df93215be1a6d53747591f5c8e162239be34a405
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
@@ -0,0 +1,345 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::wallBoundedParticle
+
+Description
+ Particle class that tracks on triangles of boundary faces. Use
+ trackToEdge similar to trackToFace on particle.
+
+SourceFiles
+ wallBoundedParticle.C
+ wallBoundedParticleTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallBoundedParticle_H
+#define wallBoundedParticle_H
+
+#include "particle.H"
+#include "autoPtr.H"
+#include "InfoProxy.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class wallBoundedParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallBoundedParticle
+:
+ public particle
+{
+
+public:
+
+ //- Class used to pass tracking data to the trackToFace function
+ template<class CloudType>
+ class TrackingData
+ :
+ public particle::TrackingData<CloudType>
+ {
+
+ public:
+
+ const PackedBoolList& isWallPatch_;
+
+ // Constructors
+
+ inline TrackingData
+ (
+ CloudType& cloud,
+ const PackedBoolList& isWallPatch
+ )
+ :
+ particle::TrackingData<CloudType>
+ (
+ cloud
+ ),
+ isWallPatch_(isWallPatch)
+ {}
+ };
+
+
+protected:
+
+ // Protected data
+
+ //- Particle is on mesh edge:
+ // const face& f = mesh.faces()[tetFace()]
+ // const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
+ // Note that this real edge
+ // is also one of the edges of the face-triangle (from
+ // tetFace()+tetPt()).
+ label meshEdgeStart_;
+
+ //- Particle is on diagonal edge:
+ // const face& f = mesh.faces()[tetFace()]
+ // label faceBasePtI = mesh.tetBasePtIs()[faceI];
+ // label diagPtI = (faceBasePtI+diagEdge_)%f.size();
+ // const edge e(f[faceBasePtI], f[diagPtI]);
+ label diagEdge_;
+
+
+ // Protected Member Functions
+
+ //- Construct current edge
+ edge currentEdge() const;
+
+ //- Check if inside current tet
+ //void checkInside() const;
+
+ //- Check if on current edge
+ //void checkOnEdge() const;
+
+ //- Check if point on triangle
+ //void checkOnTriangle(const point&) const;
+
+ //- Cross mesh edge into different face on same cell
+ void crossEdgeConnectedFace(const edge& meshEdge);
+
+ //- Cross diagonal edge into different triangle on same face,cell
+ void crossDiagonalEdge();
+
+ //- Track through single triangle
+ scalar trackFaceTri(const vector& endPosition, label& minEdgeI);
+
+ //- Is current triangle in the track direction
+ bool isTriAlongTrack(const point& endPosition) const;
+
+
+ // Patch interactions
+
+ //- Do all patch interaction
+ template<class TrackData>
+ void patchInteraction(TrackData& td, const scalar trackFraction);
+
+ //- Overridable function to handle the particle hitting a patch
+ // Executed before other patch-hitting functions
+ template<class TrackData>
+ bool hitPatch
+ (
+ const polyPatch&,
+ TrackData& td,
+ const label patchI,
+ const scalar trackFraction,
+ const tetIndices& tetIs
+ );
+
+ //- Overridable function to handle the particle hitting a wedge
+ template<class TrackData>
+ void hitWedgePatch
+ (
+ const wedgePolyPatch&,
+ TrackData& td
+ );
+
+ //- Overridable function to handle the particle hitting a
+ // symmetryPlane
+ template<class TrackData>
+ void hitSymmetryPatch
+ (
+ const symmetryPolyPatch&,
+ TrackData& td
+ );
+
+ //- Overridable function to handle the particle hitting a cyclic
+ template<class TrackData>
+ void hitCyclicPatch
+ (
+ const cyclicPolyPatch&,
+ TrackData& td
+ );
+
+ //- Overridable function to handle the particle hitting a
+ //- processorPatch
+ template<class TrackData>
+ void hitProcessorPatch
+ (
+ const processorPolyPatch&,
+ TrackData& td
+ );
+
+ //- Overridable function to handle the particle hitting a wallPatch
+ template<class TrackData>
+ void hitWallPatch
+ (
+ const wallPolyPatch&,
+ TrackData& td,
+ const tetIndices&
+ );
+
+ //- Overridable function to handle the particle hitting a polyPatch
+ template<class TrackData>
+ void hitPatch
+ (
+ const polyPatch&,
+ TrackData& td
+ );
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ wallBoundedParticle
+ (
+ const polyMesh& c,
+ const vector& position,
+ const label cellI,
+ const label tetFaceI,
+ const label tetPtI,
+ const label meshEdgeStart,
+ const label diagEdge
+ );
+
+ //- Construct from Istream
+ wallBoundedParticle
+ (
+ const polyMesh& c,
+ Istream& is,
+ bool readFields = true
+ );
+
+ //- Construct copy
+ wallBoundedParticle(const wallBoundedParticle& p);
+
+ //- Construct and return a clone
+ autoPtr<particle> clone() const
+ {
+ return autoPtr<particle>(new wallBoundedParticle(*this));
+ }
+
+ //- Factory class to read-construct particles used for
+ // parallel transfer
+ class iNew
+ {
+ const polyMesh& mesh_;
+
+ public:
+
+ iNew(const polyMesh& mesh)
+ :
+ mesh_(mesh)
+ {}
+
+ autoPtr<wallBoundedParticle> operator()
+ (
+ Istream& is
+ ) const
+ {
+ return autoPtr<wallBoundedParticle>
+ (
+ new wallBoundedParticle(mesh_, is, true)
+ );
+ }
+ };
+
+
+ // Member Functions
+
+ // Access
+
+ //- -1 or label of mesh edge
+ inline label meshEdgeStart() const
+ {
+ return meshEdgeStart_;
+ }
+
+ //- -1 or diagonal edge
+ inline label diagEdge() const
+ {
+ return diagEdge_;
+ }
+
+
+ // Track
+
+ //- Equivalent of trackToFace
+ template<class TrackData>
+ scalar trackToEdge
+ (
+ TrackData& td,
+ const vector& endPosition
+ );
+
+
+ // Info
+
+ //- Return info proxy.
+ // Used to print particle information to a stream
+ inline InfoProxy<wallBoundedParticle> info() const
+ {
+ return *this;
+ }
+
+
+ // I-O
+
+ //- Read
+ template<class CloudType>
+ static void readFields(CloudType&);
+
+ //- Write
+ template<class CloudType>
+ static void writeFields(const CloudType&);
+
+ //- Write the particle data
+ void write(Ostream& os, bool writeFields) const;
+
+
+ // Ostream Operator
+
+ friend Ostream& operator<<
+ (
+ Ostream&,
+ const wallBoundedParticle&
+ );
+
+ friend Ostream& operator<<
+ (
+ Ostream&,
+ const InfoProxy<wallBoundedParticle>&
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "wallBoundedParticleTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..f79d964c7be03dd594752ea04c613a6c099fddd4
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
@@ -0,0 +1,543 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*----------------------------------------------------------------------------*/
+
+#include "wallBoundedParticle.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template<class TrackData>
+void Foam::wallBoundedParticle::patchInteraction
+(
+ TrackData& td,
+ const scalar trackFraction
+)
+{
+// typedef TrackData::CloudType cloudType;
+ typedef typename TrackData::cloudType::particleType particleType;
+
+ particleType& p = static_cast<particleType&>(*this);
+ p.hitFace(td);
+
+ if (!internalFace(faceI_))
+ {
+ label origFaceI = faceI_;
+ label patchI = patch(faceI_);
+
+ // No action taken for tetPtI_ for tetFaceI_ here, handled by
+ // patch interaction call or later during processor transfer.
+
+
+ // Dummy tet indices. What to do here?
+ tetIndices faceHitTetIs;
+
+ if
+ (
+ !p.hitPatch
+ (
+ mesh_.boundaryMesh()[patchI],
+ td,
+ patchI,
+ trackFraction,
+ faceHitTetIs
+ )
+ )
+ {
+ // Did patch interaction model switch patches?
+ // Note: recalculate meshEdgeStart_, diagEdge_!
+ if (faceI_ != origFaceI)
+ {
+ patchI = patch(faceI_);
+ }
+
+ const polyPatch& patch = mesh_.boundaryMesh()[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, faceHitTetIs
+ );
+ }
+ else
+ {
+ p.hitPatch(patch, td);
+ }
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+//- Track particle to a given position and returns 1.0 if the
+// trajectory is completed without hitting a face otherwise
+// stops at the face and returns the fraction of the trajectory
+// completed.
+// on entry 'stepFraction()' should be set to the fraction of the
+// time-step at which the tracking starts.
+template<class TrackData>
+Foam::scalar Foam::wallBoundedParticle::trackToEdge
+(
+ TrackData& td,
+ const vector& endPosition
+)
+{
+ // Are we on a track face? If not we do a topological walk.
+
+ // Particle:
+ // - cell_ always set
+ // - tetFace_, tetPt_ always set (these identify tet particle is in)
+ // - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
+
+ //checkInside();
+ //checkOnTriangle(position());
+ //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
+ //{
+ // checkOnEdge();
+ //}
+
+ scalar trackFraction = 0.0;
+
+ if (!td.isWallPatch_[tetFace()])
+ {
+ // Don't track across face. Just walk in cell. Particle is on
+ // mesh edge (as indicated by meshEdgeStart_).
+ const edge meshEdge(currentEdge());
+
+ // If internal face check whether to go to neighbour cell or just
+ // check to the other internal tet on the edge.
+ if (mesh_.isInternalFace(tetFace()))
+ {
+ label nbrCellI =
+ (
+ cellI_ == mesh_.faceOwner()[faceI_]
+ ? mesh_.faceNeighbour()[faceI_]
+ : mesh_.faceOwner()[faceI_]
+ );
+ // Check angle to nbrCell tet. Is it in the direction of the
+ // 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)
+ {
+ // Change into nbrCell. No need to change tetFace, tetPt.
+ //Pout<< " crossed from cell:" << cellI_
+ // << " into " << nbrCellI << endl;
+ cellI_ = nbrCellI;
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Walk to other face on edge. Changes tetFace, tetPt but not
+ // cell.
+ crossEdgeConnectedFace(meshEdge);
+ patchInteraction(td, trackFraction);
+ }
+ }
+ else
+ {
+ // Walk to other face on edge. This might give loop since
+ // particle should have been removed?
+ crossEdgeConnectedFace(meshEdge);
+ patchInteraction(td, trackFraction);
+ }
+ }
+ else
+ {
+ // We're inside a tet on the wall. Check if the current tet is
+ // the one to cross. If not we cross into the neighbouring triangle.
+
+ if (mesh_.isInternalFace(tetFace()))
+ {
+ FatalErrorIn
+ (
+ "wallBoundedParticle::trackToEdge"
+ "(TrackData&, const vector&)"
+ ) << "Can only track on boundary faces."
+ << " Face:" << tetFace()
+ << " at:" << mesh_.faceCentres()[tetFace()]
+ << abort(FatalError);
+ }
+
+ point projectedEndPosition = endPosition;
+ // 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;
+ }
+
+
+ bool doTrack = false;
+ if (meshEdgeStart_ == -1 && diagEdge_ == -1)
+ {
+ // We're starting and not yet on an edge.
+ doTrack = true;
+ }
+ else
+ {
+ // See if the current triangle has got a point on the
+ // correct side of the edge.
+ doTrack = isTriAlongTrack(projectedEndPosition);
+ }
+
+
+ if (doTrack)
+ {
+ // Track across triangle. Return triangle edge crossed.
+ label triEdgeI = -1;
+ trackFraction = trackFaceTri(projectedEndPosition, triEdgeI);
+
+ if (triEdgeI == -1)
+ {
+ // Reached endpoint
+ //checkInside();
+ diagEdge_ = -1;
+ meshEdgeStart_ = -1;
+ return trackFraction;
+ }
+
+ const tetIndices ti(currentTetIndices());
+
+ // Triangle (faceTriIs) gets constructed from
+ // f[faceBasePtI_],
+ // f[facePtAI_],
+ // f[facePtBI_]
+ //
+ // So edge indices are:
+ // 0 : edge between faceBasePtI_ and facePtAI_
+ // 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
+ // 2 : edge between facePtBI_ and faceBasePtI_
+
+ const Foam::face& f = mesh_.faces()[ti.face()];
+ const label fp0 = ti.faceBasePt();
+
+ if (triEdgeI == 0)
+ {
+ if (ti.facePtA() == f.fcIndex(fp0))
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = fp0;
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else if (ti.facePtA() == f.rcIndex(fp0))
+ {
+ //Note: should not happen since boundary face so owner
+ //Pout<< "Real edge." << endl;
+ FatalErrorIn("shold not happend") << info()
+ << abort(FatalError);
+
+ diagEdge_ = -1;
+ meshEdgeStart_ = f.rcIndex(fp0);
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Get index of triangle on other side of edge.
+ diagEdge_ = ti.facePtA()-fp0;
+ if (diagEdge_ < 0)
+ {
+ diagEdge_ += f.size();
+ }
+ meshEdgeStart_ = -1;
+ //checkOnEdge();
+ crossDiagonalEdge();
+ }
+ }
+ else if (triEdgeI == 1)
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = ti.facePtA();
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else // if (triEdgeI == 2)
+ {
+ if (ti.facePtB() == f.rcIndex(fp0))
+ {
+ //Pout<< "Real edge." << endl;
+ diagEdge_ = -1;
+ meshEdgeStart_ = ti.facePtB();
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else if (ti.facePtB() == f.fcIndex(fp0))
+ {
+ //Note: should not happen since boundary face so owner
+ //Pout<< "Real edge." << endl;
+ FatalErrorIn("shold not happend") << info()
+ << abort(FatalError);
+
+ diagEdge_ = -1;
+ meshEdgeStart_ = fp0;
+ //checkOnEdge();
+ crossEdgeConnectedFace(currentEdge());
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ //Pout<< "Triangle edge." << endl;
+ // Get index of triangle on other side of edge.
+ diagEdge_ = ti.facePtB()-fp0;
+ if (diagEdge_ < 0)
+ {
+ diagEdge_ += f.size();
+ }
+ meshEdgeStart_ = -1;
+ //checkOnEdge();
+ crossDiagonalEdge();
+ }
+ }
+ }
+ else
+ {
+ // Current tet is not the right one. Check the neighbour tet.
+
+ if (meshEdgeStart_ != -1)
+ {
+ // Particle is on mesh edge so change into other face on cell
+ crossEdgeConnectedFace(currentEdge());
+ //checkOnEdge();
+ patchInteraction(td, trackFraction);
+ }
+ else
+ {
+ // Particle is on diagonal edge so change into the other
+ // triangle.
+ crossDiagonalEdge();
+ //checkOnEdge();
+ }
+ }
+ }
+
+ //checkInside();
+
+ return trackFraction;
+}
+
+
+template<class TrackData>
+bool Foam::wallBoundedParticle::hitPatch
+(
+ const polyPatch&,
+ TrackData& td,
+ const label patchI,
+ const scalar trackFraction,
+ const tetIndices& tetIs
+)
+{
+ // Disable generic patch interaction
+ return false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitWedgePatch
+(
+ const wedgePolyPatch& pp,
+ TrackData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitSymmetryPatch
+(
+ const symmetryPolyPatch& pp,
+ TrackData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitCyclicPatch
+(
+ const cyclicPolyPatch& pp,
+ TrackData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitProcessorPatch
+(
+ const processorPolyPatch& pp,
+ TrackData& td
+)
+{
+ // Switch particle
+ td.switchProcessor = true;
+
+ // Adapt edgeStart_ for other side.
+ // E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
+ // on the other side between 2 and 3 so edgeStart_ should be
+ // f.size()-edgeStart_-1.
+
+ const Foam::face& f = mesh_.faces()[face()];
+
+ if (meshEdgeStart_ != -1)
+ {
+ meshEdgeStart_ = f.size()-meshEdgeStart_-1;
+ }
+ else
+ {
+ // diagEdge_ is relative to faceBasePt
+ diagEdge_ = f.size()-diagEdge_;
+ }
+}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitWallPatch
+(
+ const wallPolyPatch& wpp,
+ TrackData& td,
+ const tetIndices&
+)
+{}
+
+
+template<class TrackData>
+void Foam::wallBoundedParticle::hitPatch
+(
+ const polyPatch& wpp,
+ TrackData& td
+)
+{
+ // Remove particle
+ td.keepParticle = false;
+}
+
+
+template<class CloudType>
+void Foam::wallBoundedParticle::readFields(CloudType& c)
+{
+ if (!c.size())
+ {
+ return;
+ }
+
+ particle::readFields(c);
+
+ IOField<label> meshEdgeStart
+ (
+ c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
+ );
+
+ IOField<label> diagEdge
+ (
+ c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
+ );
+ c.checkFieldIOobject(c, diagEdge);
+
+ label i = 0;
+ forAllIter(typename CloudType, c, iter)
+ {
+ iter().meshEdgeStart_ = meshEdgeStart[i];
+ iter().diagEdge_ = diagEdge[i];
+ i++;
+ }
+}
+
+
+template<class CloudType>
+void Foam::wallBoundedParticle::writeFields(const CloudType& c)
+{
+ particle::writeFields(c);
+
+ label np = c.size();
+
+ IOField<label> meshEdgeStart
+ (
+ c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
+ np
+ );
+ IOField<label> diagEdge
+ (
+ c.fieldIOobject("diagEdge", IOobject::NO_READ),
+ np
+ );
+
+ label i = 0;
+ forAllConstIter(typename CloudType, c, iter)
+ {
+ meshEdgeStart[i] = iter().meshEdgeStart_;
+ diagEdge[i] = iter().diagEdge_;
+ i++;
+ }
+
+ meshEdgeStart.write();
+ diagEdge.write();
+}
+
+
+// ************************************************************************* //