streamlines: Corrected handling of cyclic boundaries

By default a streamline now stops at the cyclic and starts again at the
coupled location on the opposite cyclic.

There is also now an "outside" option that can be passed to the
streamlines function. This changes the default behaviour so that the
streamline continues outside of the mesh when it encounters a cyclic
patch. The following postProcess command uses the "outside" option in
this way:

    postProcess -latestTime -func "
        streamlinesPatch
        (
            patch=inlet,
            nPoints=50,
            outside=true,
            fields=(p U)
        )"

src/functionObjects/field/streamlines/streamlines.C

@@ -270,6 +270,7 @@ void Foam::functionObjects::streamlines::track()
         UIndex,         // index of U in vvInterp
 
         trackDirection_ == trackDirection::forward,
+        trackOutside_,
 
         nSubCycle_,     // automatic track control:step through cells in steps?
         trackLength_,   // fixed track length
@@ -361,6 +362,8 @@ bool Foam::functionObjects::streamlines::read(const dictionary& dict)
         trackDirection_ = trackDirectionNames_[word(dict.lookup("direction"))];
     }
 
+    trackOutside_ = dict.lookupOrDefault<Switch>("outside", false);
+
     dict.lookup("lifeTime") >> lifeTime_;
     if (lifeTime_ < 1)
     {

src/functionObjects/field/streamlines/streamlines.H

@@ -70,6 +70,7 @@ Usage
         setFormat    | Output data type        | yes         |
         U            | Tracking velocity field name | no     | U
         direction    | Direction in which to track | yes     |
+        outside      | Track outside of periodic meshes | no | no
         fields       | Fields to sample        | yes         |
         writeTime    | Write the flow time along the streamlines | no | no
         lifetime     | Maximum number of particle tracking steps | yes |
@@ -167,6 +168,9 @@ private:
         //- The direction in which to track
         trackDirection trackDirection_;
 
+        //- Whether or not to track outside of the mesh in periodic geometries
+        Switch trackOutside_;
+
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
 

src/functionObjects/field/streamlines/streamlinesParticle.C

@@ -47,34 +47,65 @@ Foam::vector Foam::streamlinesParticle::interpolateFields
     sampledScalars_.setSize(td.vsInterp_.size());
     forAll(td.vsInterp_, scalarI)
     {
+        const scalar s =
+            td.vsInterp_[scalarI].interpolate(position, celli, facei);
+
         sampledScalars_[scalarI].append
         (
-            td.vsInterp_[scalarI].interpolate
+            td.trackOutside_ ? transform_.invTransform(s) : s
-            (
-                position,
-                celli,
-                facei
-            )
         );
     }
 
+    vector U = vector::uniform(NaN);
+
     sampledVectors_.setSize(td.vvInterp_.size());
     forAll(td.vvInterp_, vectorI)
     {
+        const vector v =
+            td.vvInterp_[vectorI].interpolate(position, celli, facei);
+
+        if (vectorI == td.UIndex_)
+        {
+            U = v;
+        }
+
         sampledVectors_[vectorI].append
         (
-            td.vvInterp_[vectorI].interpolate
+            td.trackOutside_ ? transform_.invTransform(v) : v
-            (
-                position,
-                celli,
-                facei
-            )
         );
     }
 
-    const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
+    return U;
+}
 
-    return U.last();
+
+void Foam::streamlinesParticle::endTrack(trackingData& td)
+{
+    {
+        td.allPositions_.append(vectorList());
+        vectorList& top = td.allPositions_.last();
+        top.transfer(sampledPositions_);
+    }
+
+    {
+        td.allTimes_.append(scalarList());
+        scalarList& top = td.allTimes_.last();
+        top.transfer(sampledTimes_);
+    }
+
+    forAll(sampledScalars_, i)
+    {
+        td.allScalars_[i].append(scalarList());
+        scalarList& top = td.allScalars_[i].last();
+        top.transfer(sampledScalars_[i]);
+    }
+
+    forAll(sampledVectors_, i)
+    {
+        td.allVectors_[i].append(vectorList());
+        vectorList& top = td.allVectors_[i].last();
+        top.transfer(sampledVectors_[i]);
+    }
 }
 
 
@@ -90,6 +121,7 @@ Foam::streamlinesParticle::streamlinesParticle
 :
     particle(mesh, position, celli),
     lifeTime_(lifeTime),
+    transform_(transformer::I),
     age_(0)
 {}
 
@@ -108,8 +140,8 @@ Foam::streamlinesParticle::streamlinesParticle
         List<scalarList> sampledScalars;
         List<vectorList> sampledVectors;
 
-        is  >> lifeTime_ >> age_ >> sampledPositions_ >> sampledTimes_
+        is  >> lifeTime_ >> transform_ >> age_ >> sampledPositions_
-            >> sampledScalars >> sampledVectors;
+            >> sampledTimes_ >> sampledScalars >> sampledVectors;
 
         sampledScalars_.setSize(sampledScalars.size());
         forAll(sampledScalars, i)
@@ -139,6 +171,7 @@ Foam::streamlinesParticle::streamlinesParticle
 :
     particle(p),
     lifeTime_(p.lifeTime_),
+    transform_(p.transform_),
     age_(p.age_),
     sampledPositions_(p.sampledPositions_),
     sampledTimes_(p.sampledTimes_),
@@ -173,7 +206,12 @@ bool Foam::streamlinesParticle::move
             --lifeTime_;
 
             // Store current position and sampled velocity.
-            sampledPositions_.append(position());
+            sampledPositions_.append
+            (
+                td.trackOutside_
+              ? transform_.invTransformPosition(position())
+              : position()
+            );
             sampledTimes_.append(age_);
             vector U = interpolateFields(td, position(), cell(), face());
 
@@ -216,13 +254,7 @@ bool Foam::streamlinesParticle::move
                *dt
                *(1 - trackToAndHitFace(dt*U, 0, cloud, td));
 
-            if
+            if (!td.keepParticle || td.switchProcessor || lifeTime_ == 0)
-            (
-                onFace()
-            || !td.keepParticle
-            ||  td.switchProcessor
-            ||  lifeTime_ == 0
-            )
             {
                 break;
             }
@@ -245,7 +277,12 @@ bool Foam::streamlinesParticle::move
         else
         {
             // Normal exit. Store last position and fields
-            sampledPositions_.append(position());
+            sampledPositions_.append
+            (
+                td.trackOutside_
+              ? transform_.invTransformPosition(position())
+              : position()
+            );
             sampledTimes_.append(age_);
             interpolateFields(td, position(), cell(), face());
 
@@ -257,29 +294,8 @@ bool Foam::streamlinesParticle::move
             }
         }
 
-        // Transfer particle data into trackingData.
+        // End this track
-        {
+        endTrack(td);
-            td.allPositions_.append(vectorList());
-            vectorList& top = td.allPositions_.last();
-            top.transfer(sampledPositions_);
-        }
-        {
-            td.allTimes_.append(scalarList());
-            scalarList& top = td.allTimes_.last();
-            top.transfer(sampledTimes_);
-        }
-        forAll(sampledScalars_, i)
-        {
-            td.allScalars_[i].append(scalarList());
-            scalarList& top = td.allScalars_[i].last();
-            top.transfer(sampledScalars_[i]);
-        }
-        forAll(sampledVectors_, i)
-        {
-            td.allVectors_[i].append(vectorList());
-            vectorList& top = td.allVectors_[i].last();
-            top.transfer(sampledVectors_[i]);
-        }
     }
 
     return td.keepParticle;
@@ -328,51 +344,69 @@ void Foam::streamlinesParticle::hitSymmetryPatch
 
 void Foam::streamlinesParticle::hitCyclicPatch
 (
-    streamlinesCloud&,
+    streamlinesCloud& cloud,
     trackingData& td
 )
 {
-    // Remove particle
+    const cyclicPolyPatch& cpp =
-    td.keepParticle = false;
+        static_cast<const cyclicPolyPatch&>(mesh().boundaryMesh()[patch()]);
+
+    // End this track
+    if (!td.trackOutside_ && cpp.transform().transformsPosition())
+    {
+        endTrack(td);
+    }
+
+    // Move across the cyclic
+    particle::hitCyclicPatch(cloud, td);
 }
 
 
 void Foam::streamlinesParticle::hitCyclicAMIPatch
 (
-    const vector&,
+    const vector& displacement,
-    const scalar,
+    const scalar fraction,
-    streamlinesCloud&,
+    streamlinesCloud& cloud,
     trackingData& td
 )
 {
-    // Remove particle
+    // End this track
-    td.keepParticle = false;
+    endTrack(td);
+
+    // Move across the cyclic
+    particle::hitCyclicAMIPatch(displacement, fraction, cloud, td);
 }
 
 
 void Foam::streamlinesParticle::hitCyclicACMIPatch
 (
-    const vector&,
+    const vector& displacement,
-    const scalar,
+    const scalar fraction,
-    streamlinesCloud&,
+    streamlinesCloud& cloud,
     trackingData& td
 )
 {
-    // Remove particle
+    // End this track
-    td.keepParticle = false;
+    endTrack(td);
+
+    // Move across the cyclic
+    particle::hitCyclicACMIPatch(displacement, fraction, cloud, td);
 }
 
 
 void Foam::streamlinesParticle::hitCyclicRepeatAMIPatch
 (
-    const vector&,
+    const vector& displacement,
-    const scalar,
+    const scalar fraction,
-    streamlinesCloud&,
+    streamlinesCloud& cloud,
     trackingData& td
 )
 {
-    // Remove particle
+    // End this track
-    td.keepParticle = false;
+    endTrack(td);
+
+    // Move across the cyclic
+    particle::hitCyclicRepeatAMIPatch(displacement, fraction, cloud, td);
 }
 
 
@@ -382,7 +416,16 @@ void Foam::streamlinesParticle::hitProcessorPatch
     trackingData& td
 )
 {
-    // Switch particle
+    const processorPolyPatch& ppp =
+        static_cast<const processorPolyPatch&>(mesh().boundaryMesh()[patch()]);
+
+    // End this track
+    if (!td.trackOutside_ && ppp.transform().transformsPosition())
+    {
+        endTrack(td);
+    }
+
+    // Switch processor
     td.switchProcessor = true;
 }
 
@@ -398,12 +441,17 @@ void Foam::streamlinesParticle::hitWallPatch
 }
 
 
+void Foam::streamlinesParticle::transformProperties
+(
+    const transformer& transform
+)
+{
+    transform_ = transform & transform_;
+}
+
+
 void Foam::streamlinesParticle::readFields(Cloud<streamlinesParticle>& c)
 {
-//    if (!c.size())
-//    {
-//        return;
-//    }
     bool valid = c.size();
 
     particle::readFields(c);
@@ -415,6 +463,12 @@ void Foam::streamlinesParticle::readFields(Cloud<streamlinesParticle>& c)
     );
     c.checkFieldIOobject(c, lifeTime);
 
+    IOList<transformer> transform
+    (
+        c.fieldIOobject("transform", IOobject::MUST_READ),
+        valid
+    );
+
     vectorFieldIOField sampledPositions
     (
         c.fieldIOobject("sampledPositions", IOobject::MUST_READ),
@@ -433,6 +487,7 @@ void Foam::streamlinesParticle::readFields(Cloud<streamlinesParticle>& c)
     forAllIter(Cloud<streamlinesParticle>, c, iter)
     {
         iter().lifeTime_ = lifeTime[i];
+        iter().transform_ = transform[i];
         iter().sampledPositions_.transfer(sampledPositions[i]);
         iter().sampledTimes_.transfer(sampledTimes[i]);
         i++;
@@ -451,11 +506,19 @@ void Foam::streamlinesParticle::writeFields(const Cloud<streamlinesParticle>& c)
         c.fieldIOobject("lifeTime", IOobject::NO_READ),
         np
     );
+
+    IOList<transformer> transform
+    (
+        c.fieldIOobject("transform", IOobject::NO_READ),
+        np
+    );
+
     vectorFieldIOField sampledPositions
     (
         c.fieldIOobject("sampledPositions", IOobject::NO_READ),
         np
     );
+
     scalarFieldIOField sampledTimes
     (
         c.fieldIOobject("sampledTimes", IOobject::NO_READ),
@@ -466,6 +529,7 @@ void Foam::streamlinesParticle::writeFields(const Cloud<streamlinesParticle>& c)
     forAllConstIter(Cloud<streamlinesParticle>, c, iter)
     {
         lifeTime[i] = iter().lifeTime_;
+        transform[i] = iter().transform_;
         sampledPositions[i] = iter().sampledPositions_;
         sampledTimes[i] = iter().sampledTimes_;
         i++;
@@ -483,6 +547,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const streamlinesParticle& p)
 {
     os  << static_cast<const particle&>(p)
         << token::SPACE << p.lifeTime_
+        << token::SPACE << p.transform_
         << token::SPACE << p.age_
         << token::SPACE << p.sampledPositions_
         << token::SPACE << p.sampledTimes_

src/functionObjects/field/streamlines/streamlinesParticle.H

@@ -80,6 +80,8 @@ public:
 
             bool trackForward_;
 
+            bool trackOutside_;
+
             const label nSubCycle_;
 
             const scalar trackLength_;
@@ -103,6 +105,7 @@ public:
                 const PtrList<interpolation<vector>>& vvInterp,
                 const label UIndex,
                 const bool trackForward,
+                const bool trackOutside,
                 const label nSubCycle,
                 const scalar trackLength,
                 DynamicList<List<point>>& allPositions,
@@ -116,6 +119,7 @@ public:
                 vvInterp_(vvInterp),
                 UIndex_(UIndex),
                 trackForward_(trackForward),
+                trackOutside_(trackOutside),
                 nSubCycle_(nSubCycle),
                 trackLength_(trackLength),
                 allPositions_(allPositions),
@@ -133,6 +137,9 @@ private:
         //- Lifetime of particle. Particle dies when reaches 0.
         label lifeTime_;
 
+        //- Current compound transform
+        transformer transform_;
+
         //- Age of the particle
         scalar age_;
 
@@ -160,6 +167,9 @@ private:
             const label facei
         );
 
+        //- End the current track
+        void endTrack(trackingData&);
+
 
 public:
 
@@ -276,6 +286,13 @@ public:
             void hitWallPatch(streamlinesCloud&, trackingData&);
 
 
+        // Transformations
+
+            //- Transform the physical properties of the particle
+            //  according to the given transformation tensor
+            virtual void transformProperties(const transformer&);
+
+
         // I-O
 
             //- Read

src/lagrangian/basic/particle/particleTemplates.C

@@ -412,6 +412,8 @@ void Foam::particle::hitCyclicACMIPatch
 {
     typename TrackCloudType::particleType& p =
         static_cast<typename TrackCloudType::particleType&>(*this);
+    typename TrackCloudType::particleType::trackingData& ttd =
+        static_cast<typename TrackCloudType::particleType::trackingData&>(td);
 
     const cyclicACMIPolyPatch& cpp =
         static_cast<const cyclicACMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
@@ -445,7 +447,7 @@ void Foam::particle::hitCyclicACMIPatch
 
     if (couple)
     {
-        p.hitCyclicAMIPatch(displacement, fraction, cloud, td);
+        p.hitCyclicAMIPatch(displacement, fraction, cloud, ttd);
     }
     else
     {
@@ -468,8 +470,10 @@ void Foam::particle::hitCyclicRepeatAMIPatch
 {
     typename TrackCloudType::particleType& p =
         static_cast<typename TrackCloudType::particleType&>(*this);
+    typename TrackCloudType::particleType::trackingData& ttd =
+        static_cast<typename TrackCloudType::particleType::trackingData&>(td);
 
-    p.hitCyclicAMIPatch(displacement, fraction, cloud, td);
+    p.hitCyclicAMIPatch(displacement, fraction, cloud, ttd);
 }