zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: streamLine: added preliminary support for surface following

Browse Source
This commit is contained in:
mattijs
2011-12-22 11:46:14 +00:00
parent 06699987a8
commit b752c60cd1
4 changed files with 336 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

235
src/postProcessing/functionObjects/field/streamLine/streamLine.C
View File

@ -34,6 +34,7 @@ License
#include "globalIndex.H"
#include "mapDistribute.H"
#include "interpolationCellPoint.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,6 +42,57 @@ defineTypeNameAndDebug(Foam::streamLine, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::streamLine::wallPatch() const
{
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
nFaces += patches[patchI].size();
}
}
labelList addressing(nFaces);
nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
addressing[nFaces++] = pp.start()+i;
}
}
}
return autoPtr<indirectPrimitivePatch>
(
new indirectPrimitivePatch
(
IndirectList<face>
(
mesh.faces(),
addressing
),
mesh.points()
)
);
}
void Foam::streamLine::track()
{
const Time& runTime = obr_.time();
@ -72,7 +124,7 @@ void Foam::streamLine::track()
label nSeeds = returnReduce(particles.size(), sumOp<label>());
Info<< "Seeded " << nSeeds << " particles." << endl;
Info<< type() << " : seeded " << nSeeds << " particles." << endl;
// Read or lookup fields
PtrList<volScalarField> vsFlds;
@ -163,7 +215,6 @@ void Foam::streamLine::track()
vsInterp.set
(
nScalar++,
//new interpolationCellPoint<scalar>(f)
interpolation<scalar>::New
(
interpolationScheme_,
@ -186,7 +237,6 @@ void Foam::streamLine::track()
vvInterp.set
(
nVector++,
//new interpolationCellPoint<vector>(f)
interpolation<vector>::New
(
interpolationScheme_,
@ -241,6 +291,120 @@ void Foam::streamLine::track()
allVectors_[i].setCapacity(nSeeds);
}
//- For surface following: per face the normal to constrain the velocity
// with.
pointField patchNormals;
if (followSurface_)
{
// Determine geometric data on the non-constraint patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - point normals
// - neighbouring cells
autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
// Calculate parallel consistent normal (on boundaryPatch points only)
pointField boundaryNormals
(
PatchTools::pointNormals
(
mesh,
boundaryPatch(),
boundaryPatch().addressing()
)
);
// Make sure all points know about point normal. There might be points
// which are not on boundaryPatch() but coupled to points that are.
pointField pointNormals(mesh.nPoints(), vector::zero);
UIndirectList<point>
(
pointNormals,
boundaryPatch().meshPoints()
) = boundaryNormals;
boundaryNormals.clear();
mesh.globalData().syncPointData
(
pointNormals,
maxMagSqrEqOp<point>(),
mapDistribute::transform()
);
// Calculate per-face a single constraint normal:
// - faces on patch: face-normal
// - faces on cell with point on patch: point-normal
// - rest: vector::zero
//
// Note: the disadvantage is that we can only handle faces with one
// point on the patch - otherwise we might pick up the wrong
// velocity.
patchNormals.setSize(mesh.nFaces(), vector::zero);
label nPointNormals = 0;
// 1. faces on patch
UIndirectList<point>(patchNormals, boundaryPatch().addressing()) =
boundaryPatch().faceNormals();
// 2. faces with a point on the patch
forAll(mesh.faces(), faceI)
{
if (patchNormals[faceI] == vector::zero)
{
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
label pointI = f[fp];
if (pointNormals[pointI] != vector::zero)
{
patchNormals[faceI] = pointNormals[pointI];
nPointNormals++;
}
}
}
}
// 3. cells on faces with a point on the patch
forAll(mesh.points(), pointI)
{
if (pointNormals[pointI] != vector::zero)
{
const labelList& pCells = mesh.pointCells(pointI);
forAll(pCells, i)
{
const cell& cFaces = mesh.cells()[pCells[i]];
forAll(cFaces, j)
{
label faceI = cFaces[j];
if (patchNormals[faceI] == vector::zero)
{
patchNormals[faceI] = pointNormals[pointI];
nPointNormals++;
}
}
}
}
}
Info<< type() << " : calculated constrained-tracking normals" << nl
<< "Number of faces in mesh : "
<< returnReduce(mesh.nFaces(), sumOp<label>()) << nl
<< " of which on walls : "
<< returnReduce
(
boundaryPatch().addressing().size(),
sumOp<label>()
) << nl
<< " of which on cell with point on boundary : "
<< returnReduce(nPointNormals, sumOp<label>()) << nl
<< endl;
}
// additional particle info
streamLineParticle::trackingData td
(
@ -248,8 +412,16 @@ void Foam::streamLine::track()
vsInterp,
vvInterp,
UIndex, // index of U in vvInterp
trackForward_, // track in +u direction
nSubCycle_, // step through cells in steps?
trackForward_, // track in +u direction?
nSubCycle_, // automatic track control:step through cells in steps?
trackLength_, // fixed track length
//- For surface-constrained streamlines
followSurface_, // stay on surface?
patchNormals, // per face normal to constrain tracking to
//pointNormals, // per point ,,
trackHeight_, // track height
allTracks_,
allScalars_,
allVectors_
@ -279,7 +451,8 @@ Foam::streamLine::streamLine
name_(name),
obr_(obr),
loadFromFiles_(loadFromFiles),
active_(true)
active_(true),
nSubCycle_(0)
{
// Only active if a fvMesh is available
if (isA<fvMesh>(obr_))
@ -334,6 +507,16 @@ void Foam::streamLine::read(const dictionary& dict)
dict.lookup("U") >> UName_;
}
}
if (findIndex(fields_, UName_) == -1)
{
FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
<< "Velocity field for tracking " << UName_
<< " should be present in the list of fields " << fields_
<< exit(FatalIOError);
}
dict.lookup("trackForward") >> trackForward_;
dict.lookup("lifeTime") >> lifeTime_;
if (lifeTime_ < 1)
@ -343,11 +526,49 @@ void Foam::streamLine::read(const dictionary& dict)
<< exit(FatalError);
}
dict.lookup("nSubCycle") >> nSubCycle_;
bool subCycling = dict.found("nSubCycle");
bool fixedLength = dict.found("trackLength");
if (subCycling && fixedLength)
{
FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
<< "Cannot both specify automatic time stepping (through '"
<< "nSubCycle' specification) and fixed track length (through '"
<< "trackLength')"
<< exit(FatalIOError);
}
nSubCycle_ = 1;
if (dict.readIfPresent("nSubCycle", nSubCycle_))
{
trackLength_ = VGREAT;
if (nSubCycle_ < 1)
{
nSubCycle_ = 1;
}
Info<< type() << " : automatic track length specified through"
<< " number of sub cycles : " << nSubCycle_ << nl << endl;
}
else
{
dict.lookup("trackLength") >> trackLength_;
Info<< type() << " : fixed track length specified : "
<< trackLength_ << nl << endl;
}
followSurface_ = false;
trackHeight_ = VGREAT;
if (dict.readIfPresent("followSurface", followSurface_))
{
dict.lookup("trackHeight") >> trackHeight_;
Info<< type() << " : surface-constrained streamline at "
<< trackHeight_ << " distance above the surface" << nl << endl;
}
interpolationScheme_ = dict.lookupOrDefault
(

11
src/postProcessing/functionObjects/field/streamLine/streamLine.H
View File

@ -43,6 +43,7 @@ SourceFiles
#include "vectorList.H"
#include "polyMesh.H"
#include "writer.H"
#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -95,9 +96,17 @@ class streamLine
//- Maximum lifetime (= number of cells) of particle
label lifeTime_;
//- Surface-constrained?
bool followSurface_;
scalar trackHeight_;
//- Number of subcycling steps
label nSubCycle_;
//- Track length
scalar trackLength_;
//- Optional specified name of particles
word cloudName_;
@ -141,6 +150,8 @@ class streamLine
List<DynamicList<vectorList> > allVectors_;
//- Construct patch out of all wall patch faces
autoPtr<indirectPrimitivePatch> wallPatch() const;
//- Do all seeding and tracking
void track();

74
src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
View File

@ -43,7 +43,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
const vector& U
) const
{
streamLineParticle testParticle(*this);
particle testParticle(*this);
bool oldKeepParticle = td.keepParticle;
bool oldSwitchProcessor = td.switchProcessor;
@ -59,7 +59,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
(
const trackingData& td,
const point& position,
const label cellI
const label cellI,
const label faceI
)
{
if (cellI == -1)
@ -76,7 +77,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vsInterp_[scalarI].interpolate
(
position,
cellI
cellI,
faceI
)
);
}
@ -89,7 +91,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
td.vvInterp_[vectorI].interpolate
(
position,
cellI
cellI,
faceI
)
);
}
@ -201,13 +204,24 @@ bool Foam::streamLineParticle::move
// Store current position and sampled velocity.
sampledPositions_.append(position());
vector U = interpolateFields(td, position(), cell());
vector U = interpolateFields(td, position(), cell(), face());
if (!td.trackForward_)
{
U = -U;
}
if (td.followSurface_)
{
//- Simple: remove wall-normal component. Should keep particle
// at same height. Does not work well on warped faces.
const vector& n = td.patchNormals_[tetFaceI_];
if (n != vector::zero)
{
U -= (U&n)*n;
}
}
scalar magU = mag(U);
if (magU < SMALL)
@ -219,7 +233,13 @@ bool Foam::streamLineParticle::move
U /= magU;
if (subIter == 0 && td.nSubCycle_ > 1)
if (td.trackLength_ < GREAT)
{
dt = td.trackLength_;
//Pout<< " subiteration " << subIter
// << " : fixed length: updated dt:" << dt << endl;
}
else if (subIter == 0 && td.nSubCycle_ > 1)
{
// Adapt dt to cross cell in a few steps
dt = calcSubCycleDeltaT(td, dt, U);
@ -264,7 +284,8 @@ bool Foam::streamLineParticle::move
if (debug)
{
Pout<< "streamLineParticle : Removing stagnant particle:"
<< p << " sampled positions:" << sampledPositions_.size()
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
td.keepParticle = false;
@ -273,12 +294,13 @@ bool Foam::streamLineParticle::move
{
// Normal exit. Store last position and fields
sampledPositions_.append(position());
interpolateFields(td, position(), cell());
interpolateFields(td, position(), cell(), face());
if (debug)
{
Pout<< "streamLineParticle : Removing particle:"
<< p << " sampled positions:" << sampledPositions_.size()
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
}
@ -373,10 +395,44 @@ void Foam::streamLineParticle::hitWallPatch
trackingData& td,
const tetIndices&
)
{
if (td.followSurface_)
{
// Push slightly in. Since we're pushing to tet centre we should still
// be in the same tet.
tetPointRef tet = currentTet();
const point oldPosition(position());
position() += trackingCorrectionTol*(tet.centre() - position());
if (debug)
{
label patchID = patch(face());
Pout<< "hitWallPatch : on wall "
<< mesh().boundaryMesh()[patchID].name()
<< " moved from " << oldPosition
<< " to " << position()
<< " due to centre " << tet.centre() << endl;
}
// Check that still in tet (should always be the case
if (debug && !tet.inside(position()))
{
FatalErrorIn("streamLineParticle::hitWallPatch()")
<< "Pushing particle " << *this
<< " away from wall " << wpp.name()
<< " moved it outside of tet."
<< exit(FatalError);
}
}
else
{
// Remove particle
td.keepParticle = false;
}
}
void Foam::streamLineParticle::hitPatch

29
src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
View File

@ -73,6 +73,13 @@ public:
const label UIndex_;
const bool trackForward_;
const label nSubCycle_;
const scalar trackLength_;
//- For surface-constrained tracking
const bool followSurface_;
const pointField& patchNormals_;
const scalar trackHeight_;
DynamicList<vectorList>& allPositions_;
List<DynamicList<scalarList> >& allScalars_;
@ -89,6 +96,12 @@ public:
const label UIndex,
const bool trackForward,
const label nSubCycle,
const scalar trackLength,
const bool followSurface,
const pointField& patchNormals,
const scalar trackHeight,
DynamicList<List<point> >& allPositions,
List<DynamicList<scalarList> >& allScalars,
List<DynamicList<vectorList> >& allVectors
@ -100,6 +113,12 @@ public:
UIndex_(UIndex),
trackForward_(trackForward),
nSubCycle_(nSubCycle),
trackLength_(trackLength),
followSurface_(followSurface),
patchNormals_(patchNormals),
trackHeight_(trackHeight),
allPositions_(allPositions),
allScalars_(allScalars),
allVectors_(allVectors)
@ -135,12 +154,20 @@ private:
const vector& U
) const;
void constrainVelocity
(
trackingData& td,
const scalar dt,
vector& U
);
//- Interpolate all quantities; return interpolated velocity.
vector interpolateFields
(
const trackingData&,
const point&,
const label cellI
const label cellI,
const label faceI
);