From fb7698f6a1578e5ef654e0978710eac59f4a56ec Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 23 Oct 2019 16:52:27 +0100
Subject: [PATCH] ENH: streamLine: add bidirectional tracking. Fixes #1466.

---
 .../field/streamLine/streamLine.C             | 24 +++++-
 .../field/streamLine/streamLine.H             |  4 +-
 .../field/streamLine/streamLineBase.C         | 36 +++++++-
 .../field/streamLine/streamLineBase.H         | 21 ++++-
 .../field/streamLine/streamLineParticle.C     | 23 +++---
 .../field/streamLine/streamLineParticle.H     | 11 +--
 .../wallBoundedStreamLine.C                   | 82 +++++++++++++++----
 .../wallBoundedStreamLine.H                   | 14 +++-
 .../wallBoundedStreamLineParticle.C           |  8 +-
 .../wallBoundedStreamLineParticle.H           |  7 +-
 10 files changed, 184 insertions(+), 46 deletions(-)

diff --git a/src/functionObjects/field/streamLine/streamLine.C b/src/functionObjects/field/streamLine/streamLine.C
index 881330886d..f2bbacfdf0 100644
--- a/src/functionObjects/field/streamLine/streamLine.C
+++ b/src/functionObjects/field/streamLine/streamLine.C
@@ -56,18 +56,35 @@ void Foam::functionObjects::streamLine::track()
 
     const sampledSet& seedPoints = sampledSetPoints();
 
-    forAll(seedPoints, i)
+    forAll(seedPoints, seedi)
     {
         particles.addParticle
         (
             new streamLineParticle
             (
                 mesh_,
-                seedPoints[i],
-                seedPoints.cells()[i],
+                seedPoints[seedi],
+                seedPoints.cells()[seedi],
+                (trackDir_ == trackDirType::FORWARD),
                 lifeTime_
             )
         );
+
+        if (trackDir_ == trackDirType::BIDIRECTIONAL)
+        {
+            // Add additional particle for the forward bit of the track
+            particles.addParticle
+            (
+                new streamLineParticle
+                (
+                    mesh_,
+                    seedPoints[seedi],
+                    seedPoints.cells()[seedi],
+                    true,
+                    lifeTime_
+                )
+            );
+        }
     }
 
     label nSeeds = returnReduce(particles.size(), sumOp<label>());
@@ -99,7 +116,6 @@ void Foam::functionObjects::streamLine::track()
         vsInterp,
         vvInterp,
         UIndex,         // index of U in vvInterp
-        trackForward_,  // track in +u direction?
         nSubCycle_,     // automatic track control:step through cells in steps?
         trackLength_,   // fixed track length
 
diff --git a/src/functionObjects/field/streamLine/streamLine.H b/src/functionObjects/field/streamLine/streamLine.H
index 0f4885a7c7..9457abcf66 100644
--- a/src/functionObjects/field/streamLine/streamLine.H
+++ b/src/functionObjects/field/streamLine/streamLine.H
@@ -45,7 +45,8 @@ Usage
 
         setFormat       vtk;
         U               U;
-        trackForward    yes;
+        //Deprecated: trackForward    yes;
+        direction       bidirectional;    // or forward/backward
 
         fields
         (
@@ -77,6 +78,7 @@ Usage
         setFormat    | Output data type        | yes         |
         U            | Tracking velocity field name | no     | U
         fields       | Fields to sample        | yes         |
+        direction    | Direction to track    | yes         |
         lifetime     | Maximum number of particle tracking steps | yes |
         trackLength  | Tracking segment length | no          |
         nSubCycle    | Number of tracking steps per cell | no|
diff --git a/src/functionObjects/field/streamLine/streamLineBase.C b/src/functionObjects/field/streamLine/streamLineBase.C
index c1f31c200b..055ec83bfc 100644
--- a/src/functionObjects/field/streamLine/streamLineBase.C
+++ b/src/functionObjects/field/streamLine/streamLineBase.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2015 OpenFOAM Foundation
@@ -47,6 +47,18 @@ namespace functionObjects
 }
 
 
+const Foam::Enum
+<
+    Foam::functionObjects::streamLineBase::trackDirType
+>
+Foam::functionObjects::streamLineBase::trackDirTypeNames
+({
+    { trackDirType::FORWARD, "forward" },
+    { trackDirType::BACKWARD, "backward" },
+    { trackDirType::BIDIRECTIONAL, "bidirectional" }
+});
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 const Foam::word&
@@ -886,7 +898,27 @@ bool Foam::functionObjects::streamLineBase::read(const dictionary& dict)
 
     Info<< "    Employing velocity field " << UName_ << endl;
 
-    dict.readEntry("trackForward", trackForward_);
+    bool trackForward;
+    if (dict.readIfPresent("trackForward", trackForward))
+    {
+        trackDir_ =
+        (
+            trackForward
+          ? trackDirType::FORWARD
+          : trackDirType::BACKWARD
+        );
+
+        if (dict.found("direction"))
+        {
+            FatalIOErrorInFunction(dict) << "Cannot specify both "
+                << "\"trackForward\" and \"direction\""
+                << exit(FatalIOError);
+        }
+    }
+    else
+    {
+        trackDir_ = trackDirTypeNames.get("direction", dict);
+    }
     dict.readEntry("lifeTime", lifeTime_);
     if (lifeTime_ < 1)
     {
diff --git a/src/functionObjects/field/streamLine/streamLineBase.H b/src/functionObjects/field/streamLine/streamLineBase.H
index 776696f768..9ab8ec4998 100644
--- a/src/functionObjects/field/streamLine/streamLineBase.H
+++ b/src/functionObjects/field/streamLine/streamLineBase.H
@@ -45,6 +45,7 @@ SourceFiles
 #include "writer.H"
 #include "indirectPrimitivePatch.H"
 #include "interpolation.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,6 +66,22 @@ class streamLineBase
 :
     public fvMeshFunctionObject
 {
+public:
+
+    // Public data types
+
+        //- Enumeration defining the track direction
+        enum trackDirType : char
+        {
+            FORWARD,
+            BACKWARD,
+            BIDIRECTIONAL
+        };
+
+        //- Names for the trackDir
+        static const Enum<trackDirType> trackDirTypeNames;
+
+
 protected:
 
         //- Seed set engine
@@ -85,8 +102,8 @@ protected:
         //- Interpolation scheme to use
         word interpolationScheme_;
 
-        //- Whether to use +u or -u
-        bool trackForward_;
+        //- Whether to use +u or -u or both
+        trackDirType trackDir_;
 
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.C b/src/functionObjects/field/streamLine/streamLineParticle.C
index 385747ca1a..d7c1a7db26 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/functionObjects/field/streamLine/streamLineParticle.C
@@ -86,10 +86,12 @@ Foam::streamLineParticle::streamLineParticle
     const polyMesh& mesh,
     const vector& position,
     const label celli,
+    const bool trackForward,
     const label lifeTime
 )
 :
     particle(mesh, position, celli),
+    trackForward_(trackForward),
     lifeTime_(lifeTime)
 {}
 
@@ -98,17 +100,19 @@ Foam::streamLineParticle::streamLineParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields)
+    particle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
         List<scalarList> sampledScalars;
         List<vectorList> sampledVectors;
 
-        is  >> lifeTime_ >> sampledPositions_ >> sampledScalars
+        is  >> trackForward_ >> lifeTime_
+            >> sampledPositions_ >> sampledScalars
             >> sampledVectors;
 
         sampledScalars_.setSize(sampledScalars.size());
@@ -133,9 +137,11 @@ Foam::streamLineParticle::streamLineParticle
 )
 :
     particle(p),
+    trackForward_(p.trackForward_),
     lifeTime_(p.lifeTime_),
     sampledPositions_(p.sampledPositions_),
-    sampledScalars_(p.sampledScalars_)
+    sampledScalars_(p.sampledScalars_),
+    sampledVectors_(p.sampledVectors_)
 {}
 
 
@@ -168,7 +174,7 @@ bool Foam::streamLineParticle::move
             sampledPositions_.append(position());
             vector U = interpolateFields(td, position(), cell(), face());
 
-            if (!td.trackForward_)
+            if (!trackForward_)
             {
                 U = -U;
             }
@@ -366,11 +372,7 @@ void Foam::streamLineParticle::hitWallPatch
 
 void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
 {
-//    if (!c.size())
-//    {
-//        return;
-//    }
-    bool valid = c.size();
+    const bool valid = c.size();
 
     particle::readFields(c);
 
@@ -433,6 +435,7 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
 Foam::Ostream& Foam::operator<<(Ostream& os, const streamLineParticle& p)
 {
     os  << static_cast<const particle&>(p)
+        << token::SPACE << p.trackForward_
         << token::SPACE << p.lifeTime_
         << token::SPACE << p.sampledPositions_
         << token::SPACE << p.sampledScalars_
diff --git a/src/functionObjects/field/streamLine/streamLineParticle.H b/src/functionObjects/field/streamLine/streamLineParticle.H
index 3763472bd1..bc6acb6bb7 100644
--- a/src/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/functionObjects/field/streamLine/streamLineParticle.H
@@ -79,8 +79,6 @@ public:
 
             const label UIndex_;
 
-            const bool trackForward_;
-
             const label nSubCycle_;
 
             const scalar trackLength_;
@@ -101,7 +99,6 @@ public:
                 const PtrList<interpolation<scalar>>& vsInterp,
                 const PtrList<interpolation<vector>>& vvInterp,
                 const label UIndex,
-                const bool trackForward,
                 const label nSubCycle,
                 const scalar trackLength,
                 DynamicList<List<point>>& allPositions,
@@ -113,7 +110,6 @@ public:
                 vsInterp_(vsInterp),
                 vvInterp_(vvInterp),
                 UIndex_(UIndex),
-                trackForward_(trackForward),
                 nSubCycle_(nSubCycle),
                 trackLength_(trackLength),
                 allPositions_(allPositions),
@@ -127,6 +123,9 @@ private:
 
     // Private data
 
+        //- Whether particle transports with +U or -U
+        bool trackForward_;
+
         //- Lifetime of particle. Particle dies when reaches 0.
         label lifeTime_;
 
@@ -162,6 +161,7 @@ public:
             const polyMesh& c,
             const vector& position,
             const label celli,
+            const bool trackForward,
             const label lifeTime
         );
 
@@ -170,7 +170,8 @@ public:
         (
             const polyMesh& c,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct copy
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
index f4e61cbdb4..3d49eb97a9 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
@@ -51,7 +51,8 @@ namespace functionObjects
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
+Foam::Tuple2<Foam::tetIndices, Foam::point>
+Foam::functionObjects::wallBoundedStreamLine::findNearestTet
 (
     const bitSet& isWallPatch,
     const point& seedPt,
@@ -63,6 +64,7 @@ Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
     label minFacei = -1;
     label minTetPti = -1;
     scalar minDistSqr = sqr(GREAT);
+    point nearestPt(GREAT, GREAT, GREAT);
 
     for (label facei : cFaces)
     {
@@ -76,14 +78,17 @@ Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
             for (label i = 2; i < f.size(); i++)
             {
                 const point& thisPoint = mesh_.points()[f[fp]];
-                label nextFp = f.fcIndex(fp);
+                const label nextFp = f.fcIndex(fp);
                 const point& nextPoint = mesh_.points()[f[nextFp]];
 
                 const triPointRef tri(basePoint, thisPoint, nextPoint);
 
-                scalar d2 = magSqr(tri.centre() - seedPt);
+                //const scalar d2 = magSqr(tri.centre() - seedPt);
+                const pointHit nearInfo(tri.nearestPoint(seedPt));
+                const scalar d2 = nearInfo.distance();
                 if (d2 < minDistSqr)
                 {
+                    nearestPt = nearInfo.rawPoint();
                     minDistSqr = d2;
                     minFacei = facei;
                     minTetPti = i-1;
@@ -93,16 +98,37 @@ Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
         }
     }
 
-    // Put particle in tet
-    return tetIndices
+    // Return tet and nearest point on wall triangle
+    return Tuple2<Foam::tetIndices, Foam::point>
     (
-        celli,
-        minFacei,
-        minTetPti
+        tetIndices
+        (
+            celli,
+            minFacei,
+            minTetPti
+        ),
+        nearestPt
     );
 }
 
 
+Foam::point Foam::functionObjects::wallBoundedStreamLine::pushIn
+(
+    const triPointRef& tri,
+    const point& pt
+) const
+{
+    //pointHit nearPt(tri.nearestPoint(pt));
+    //if (nearPt.distance() > ROOTSMALL)
+    //{
+    //    FatalErrorInFunction << "tri:" << tri
+    //        << " seed:" << pt << exit(FatalError);
+    //}
+
+    return (1.0-ROOTSMALL)*pt+ROOTSMALL*tri.centre();
+}
+
+
 void Foam::functionObjects::wallBoundedStreamLine::track()
 {
     // Determine the 'wall' patches
@@ -130,14 +156,18 @@ void Foam::functionObjects::wallBoundedStreamLine::track()
 
         const sampledSet& seedPoints = sampledSetPoints();
 
-        forAll(seedPoints, i)
+        forAll(seedPoints, seedi)
         {
-            const label celli = seedPoints.cells()[i];
+            const label celli = seedPoints.cells()[seedi];
 
             if (celli != -1)
             {
-                const point& seedPt = seedPoints[i];
-                tetIndices ids(findNearestTet(isWallPatch, seedPt, celli));
+                const point& seedPt = seedPoints[seedi];
+                Tuple2<tetIndices, point> nearestId
+                (
+                    findNearestTet(isWallPatch, seedPt, celli)
+                );
+                const tetIndices& ids = nearestId.first();
 
                 if (ids.face() != -1 && isWallPatch[ids.face()])
                 {
@@ -154,15 +184,38 @@ void Foam::functionObjects::wallBoundedStreamLine::track()
                         new wallBoundedStreamLineParticle
                         (
                             mesh_,
-                            ids.faceTri(mesh_).centre(),
+                            // Perturb seed point to be inside triangle
+                            pushIn(ids.faceTri(mesh_), nearestId.second()),
                             ids.cell(),
                             ids.face(),     // tetFace
                             ids.tetPt(),
                             -1,             // not on a mesh edge
                             -1,             // not on a diagonal edge
+                            (trackDir_ == trackDirType::FORWARD), // forward?
                             lifeTime_       // lifetime
                         )
                     );
+
+                    if (trackDir_ == trackDirType::BIDIRECTIONAL)
+                    {
+                        // Additional particle for other half of track
+                        particles.addParticle
+                        (
+                            new wallBoundedStreamLineParticle
+                            (
+                                mesh_,
+                                // Perturb seed point to be inside triangle
+                                pushIn(ids.faceTri(mesh_), nearestId.second()),
+                                ids.cell(),
+                                ids.face(),     // tetFace
+                                ids.tetPt(),
+                                -1,             // not on a mesh edge
+                                -1,             // not on a diagonal edge
+                                true,           // forward
+                                lifeTime_       // lifetime
+                            )
+                        );
+                    }
                 }
                 else
                 {
@@ -173,7 +226,7 @@ void Foam::functionObjects::wallBoundedStreamLine::track()
         }
     }
 
-    label nSeeds = returnReduce(particles.size(), sumOp<label>());
+    const label nSeeds = returnReduce(particles.size(), sumOp<label>());
 
     Log << type() << " : seeded " << nSeeds << " particles." << endl;
 
@@ -204,7 +257,6 @@ void Foam::functionObjects::wallBoundedStreamLine::track()
         vsInterp,
         vvInterp,
         UIndex,         // index of U in vvInterp
-        trackForward_,  // track in +u direction?
         trackLength_,   // fixed track length
         isWallPatch,    // which faces are to follow
 
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
index 1b64c51722..078ed93b48 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
@@ -46,7 +46,8 @@ Usage
 
         setFormat       vtk;
         U               UNear;
-        trackForward    yes;
+        //Deprecated: trackForward    yes;
+        direction       bidirectional;    // or forward/backward
 
         fields
         (
@@ -76,6 +77,7 @@ Usage
         setFormat    | Output data type        | yes         |
         U            | Tracking velocity field name | yes    |
         fields       | Fields to sample        | yes         |
+        direction    | Direction to track    | yes         |
         lifetime     | Maximum number of particle tracking steps | yes |
         trackLength  | Tracking segment length | no          |
         nSubCycle    | Number of tracking steps per cell | no|
@@ -139,13 +141,21 @@ protected:
     // Protected Member Functions
 
         //- Find wall tet on cell
-        tetIndices findNearestTet
+        Tuple2<tetIndices, point> findNearestTet
         (
             const bitSet& isWallPatch,
             const point& seedPt,
             const label celli
         ) const;
 
+        //- Push a point a tiny bit towards the centre of the triangle it is in
+        //  to avoid tracking problems
+        point pushIn
+        (
+            const triPointRef& tri,
+            const point& pt
+        ) const;
+
 
 public:
 
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
index 7e5a925380..1e17dff017 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -96,7 +96,7 @@ Foam::vector Foam::wallBoundedStreamLineParticle::sample
 {
     vector U = interpolateFields(td, localPosition_, cell(), face());
 
-    if (!td.trackForward_)
+    if (!trackForward_)
     {
         U = -U;
     }
@@ -127,6 +127,7 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
     const label tetPti,
     const label meshEdgeStart,
     const label diagEdge,
+    const bool trackForward,
     const label lifeTime
 )
 :
@@ -140,6 +141,7 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
         meshEdgeStart,
         diagEdge
     ),
+    trackForward_(trackForward),
     lifeTime_(lifeTime)
 {}
 
@@ -159,7 +161,7 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
         List<scalarList> sampledScalars;
         List<vectorList> sampledVectors;
 
-        is  >> lifeTime_
+        is  >> trackForward_ >> lifeTime_
             >> sampledPositions_ >> sampledScalars >> sampledVectors;
 
         sampledScalars_.setSize(sampledScalars.size());
@@ -184,6 +186,7 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
 )
 :
     wallBoundedParticle(p),
+    trackForward_(p.trackForward_),
     lifeTime_(p.lifeTime_),
     sampledPositions_(p.sampledPositions_),
     sampledScalars_(p.sampledScalars_),
@@ -269,6 +272,7 @@ Foam::Ostream& Foam::operator<<
 )
 {
     os  << static_cast<const wallBoundedParticle&>(p)
+        << token::SPACE << p.trackForward_
         << token::SPACE << p.lifeTime_
         << token::SPACE << p.sampledPositions_
         << token::SPACE << p.sampledScalars_
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
index 4a0f74b30a..fc2cef0c95 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
@@ -82,7 +82,6 @@ public:
         const PtrList<interpolation<scalar>>& vsInterp_;
         const PtrList<interpolation<vector>>& vvInterp_;
         const label UIndex_;
-        const bool trackForward_;
         const scalar trackLength_;
 
         DynamicList<vectorList>& allPositions_;
@@ -99,7 +98,6 @@ public:
                 const PtrList<interpolation<scalar>>& vsInterp,
                 const PtrList<interpolation<vector>>& vvInterp,
                 const label UIndex,
-                const bool trackForward,
                 const scalar trackLength,
                 const bitSet& isWallPatch,
 
@@ -116,7 +114,6 @@ public:
                 vsInterp_(vsInterp),
                 vvInterp_(vvInterp),
                 UIndex_(UIndex),
-                trackForward_(trackForward),
                 trackLength_(trackLength),
 
                 allPositions_(allPositions),
@@ -132,6 +129,9 @@ protected:
 
     // Protected data
 
+        //- Track with +U or -U
+        bool trackForward_;
+
         //- Lifetime of particle. Particle dies when reaches 0.
         label lifeTime_;
 
@@ -172,6 +172,7 @@ public:
             const label tetPti,
             const label meshEdgeStart,
             const label diagEdge,
+            const bool trackForward,
             const label lifeTime
         );