From 701fc2b6a58dffe14c64b068f6cc13ab8c360581 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 23 Nov 2021 17:02:43 +0000 Subject: [PATCH] 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) )" --- .../field/streamlines/streamlines.C | 3 + .../field/streamlines/streamlines.H | 4 + .../field/streamlines/streamlinesParticle.C | 207 ++++++++++++------ .../field/streamlines/streamlinesParticle.H | 17 ++ .../basic/particle/particleTemplates.C | 8 +- 5 files changed, 166 insertions(+), 73 deletions(-) diff --git a/src/functionObjects/field/streamlines/streamlines.C b/src/functionObjects/field/streamlines/streamlines.C index e7fc94ed52..8b45e6f4f2 100644 --- a/src/functionObjects/field/streamlines/streamlines.C +++ b/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("outside", false); + dict.lookup("lifeTime") >> lifeTime_; if (lifeTime_ < 1) { diff --git a/src/functionObjects/field/streamlines/streamlines.H b/src/functionObjects/field/streamlines/streamlines.H index 53861c5351..7d9e408dc4 100644 --- a/src/functionObjects/field/streamlines/streamlines.H +++ b/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_; diff --git a/src/functionObjects/field/streamlines/streamlinesParticle.C b/src/functionObjects/field/streamlines/streamlinesParticle.C index d42d5f22e9..6dc7785c49 100644 --- a/src/functionObjects/field/streamlines/streamlinesParticle.C +++ b/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 - ( - position, - celli, - facei - ) + td.trackOutside_ ? transform_.invTransform(s) : s ); } + 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 - ( - position, - celli, - facei - ) + td.trackOutside_ ? transform_.invTransform(v) : v ); } - const DynamicList& 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 sampledScalars; List sampledVectors; - is >> lifeTime_ >> age_ >> sampledPositions_ >> sampledTimes_ - >> sampledScalars >> sampledVectors; + is >> lifeTime_ >> transform_ >> age_ >> sampledPositions_ + >> 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 - ( - onFace() - || !td.keepParticle - || td.switchProcessor - || lifeTime_ == 0 - ) + if (!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. - { - 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]); - } + // End this track + endTrack(td); } return td.keepParticle; @@ -328,51 +344,69 @@ void Foam::streamlinesParticle::hitSymmetryPatch void Foam::streamlinesParticle::hitCyclicPatch ( - streamlinesCloud&, + streamlinesCloud& cloud, trackingData& td ) { - // Remove particle - td.keepParticle = false; + const cyclicPolyPatch& cpp = + static_cast(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 scalar, - streamlinesCloud&, + const vector& displacement, + const scalar fraction, + streamlinesCloud& cloud, trackingData& td ) { - // Remove particle - td.keepParticle = false; + // End this track + endTrack(td); + + // Move across the cyclic + particle::hitCyclicAMIPatch(displacement, fraction, cloud, td); } void Foam::streamlinesParticle::hitCyclicACMIPatch ( - const vector&, - const scalar, - streamlinesCloud&, + const vector& displacement, + const scalar fraction, + streamlinesCloud& cloud, trackingData& td ) { - // Remove particle - td.keepParticle = false; + // End this track + endTrack(td); + + // Move across the cyclic + particle::hitCyclicACMIPatch(displacement, fraction, cloud, td); } void Foam::streamlinesParticle::hitCyclicRepeatAMIPatch ( - const vector&, - const scalar, - streamlinesCloud&, + const vector& displacement, + const scalar fraction, + streamlinesCloud& cloud, trackingData& td ) { - // Remove particle - td.keepParticle = false; + // End this track + 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(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& c) { -// if (!c.size()) -// { -// return; -// } bool valid = c.size(); particle::readFields(c); @@ -415,6 +463,12 @@ void Foam::streamlinesParticle::readFields(Cloud& c) ); c.checkFieldIOobject(c, lifeTime); + IOList 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& c) forAllIter(Cloud, 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& c) c.fieldIOobject("lifeTime", IOobject::NO_READ), np ); + + IOList 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& c) forAllConstIter(Cloud, 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(p) << token::SPACE << p.lifeTime_ + << token::SPACE << p.transform_ << token::SPACE << p.age_ << token::SPACE << p.sampledPositions_ << token::SPACE << p.sampledTimes_ diff --git a/src/functionObjects/field/streamlines/streamlinesParticle.H b/src/functionObjects/field/streamlines/streamlinesParticle.H index c4fc5f315e..c4055ab305 100644 --- a/src/functionObjects/field/streamlines/streamlinesParticle.H +++ b/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>& vvInterp, const label UIndex, const bool trackForward, + const bool trackOutside, const label nSubCycle, const scalar trackLength, DynamicList>& 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 diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index 0578a38ec5..a2367a1275 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -412,6 +412,8 @@ void Foam::particle::hitCyclicACMIPatch { typename TrackCloudType::particleType& p = static_cast(*this); + typename TrackCloudType::particleType::trackingData& ttd = + static_cast(td); const cyclicACMIPolyPatch& cpp = static_cast(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(*this); + typename TrackCloudType::particleType::trackingData& ttd = + static_cast(td); - p.hitCyclicAMIPatch(displacement, fraction, cloud, td); + p.hitCyclicAMIPatch(displacement, fraction, cloud, ttd); }