Skip to content
GitLab
Commit fb7698f6 authored by mattijs's avatar mattijs Committed by Andrew Heather
Browse files

ENH: streamLine: add bidirectional tracking. Fixes #1466.

parent 88ebf110
Hide whitespace changes
Inline Side-by-side
src/functionObjects/field/streamLine/streamLine.C
View file @ fb7698f6
......@@ -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
......
src/functionObjects/field/streamLine/streamLine.H
View file @ fb7698f6
......@@ -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|
......
src/functionObjects/field/streamLine/streamLineBase.C
View file @ fb7698f6
......@@ -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)
{
......
src/functionObjects/field/streamLine/streamLineBase.H
View file @ fb7698f6
......@@ -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_;
......
src/functionObjects/field/streamLine/streamLineParticle.C
View file @ fb7698f6
......@@ -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_
......
src/functionObjects/field/streamLine/streamLineParticle.H
View file @ fb7698f6
......@@ -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
......
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
View file @ fb7698f6
......@@ -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
......
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.H
View file @ fb7698f6
......@@ -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:
......
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
View file @ fb7698f6
......@@ -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_
......
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
View file @ fb7698f6
......@@ -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
);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment