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: Reinstated the wallBoundedStreamline function object

Browse Source 
Note: performs its own tracking and does not rely on the base
particle::trackXXX functions, and uses a local particle position.

Look to update to barycentric tracking in the future.
This commit is contained in:
Andrew Heather
2017-09-14 12:02:03 +01:00
parent 2defba00a9
commit d7fd550e61
15 changed files with 545 additions and 609 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View File

@ -142,18 +142,30 @@ public:
inline label& tetPt();
//- Return the indices corresponding to the tri on the face for
// this tet. The normal of the tri points out of the cell.
inline triFace faceTriIs(const polyMesh& mesh) const;
// this tet. The normal of the tri points out of the cell
inline triFace faceTriIs
(
const polyMesh& mesh,
const bool warn = true
) const;
//- Return the local indices corresponding to the tri on the face
// for this tet. The normal of the tri points out of the cell
inline triFace triIs
(
const polyMesh& mesh,
const bool warn = true
) const;
//- Return the geometry corresponding to this tet
inline tetPointRef tet(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the face for
// this tet. The normal of the tri points out of the cell.
// this tet. The normal of the tri points out of the cell
inline triPointRef faceTri(const polyMesh& mesh) const;
//- Return the geometry corresponding to the tri on the face for
// this tet using the old positions.
// this tet using the old positions
inline triPointRef oldFaceTri(const polyMesh& mesh) const;

78
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
View File

@ -61,7 +61,11 @@ inline Foam::label& Foam::tetIndices::tetPt()
}
inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
inline Foam::triFace Foam::tetIndices::faceTriIs
(
const polyMesh& mesh,
const bool warn
) const
{
const Foam::face& f = mesh.faces()[face()];
@ -70,22 +74,26 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
if (faceBasePtI < 0)
{
faceBasePtI = 0;
if (nWarnings < maxNWarnings)
if (warn)
{
WarningInFunction
<< "No base point for face " << face() << ", " << f
<< ", produces a valid tet decomposition." << endl;
++ nWarnings;
}
if (nWarnings == maxNWarnings)
{
Warning
<< "Suppressing any further warnings." << endl;
++ nWarnings;
if (nWarnings < maxNWarnings)
{
WarningInFunction
<< "No base point for face " << face() << ", " << f
<< ", produces a valid tet decomposition." << endl;
++nWarnings;
}
if (nWarnings == maxNWarnings)
{
Warning
<< "Suppressing any further warnings." << endl;
++nWarnings;
}
}
}
label facePtI = (tetPt() + faceBasePtI) % f.size();
label facePtI = (tetPti_ + faceBasePtI) % f.size();
label faceOtherPtI = f.fcIndex(facePtI);
if (mesh.faceOwner()[face()] != cell())
@ -97,6 +105,50 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
}
inline Foam::triFace Foam::tetIndices::triIs
(
const polyMesh& mesh,
const bool warn
) const
{
const Foam::face& f = mesh.faces()[face()];
label faceBasePtI = mesh.tetBasePtIs()[face()];
if (faceBasePtI < 0)
{
faceBasePtI = 0;
if (warn)
{
if (nWarnings < maxNWarnings)
{
WarningInFunction
<< "No base point for face " << face() << ", " << f
<< ", produces a valid tet decomposition." << endl;
++nWarnings;
}
if (nWarnings == maxNWarnings)
{
Warning
<< "Suppressing any further warnings." << endl;
++nWarnings;
}
}
}
label facePtI = (tetPti_ + faceBasePtI) % f.size();
label faceOtherPtI = f.fcIndex(facePtI);
if (mesh.faceOwner()[face()] != cell())
{
Swap(facePtI, faceOtherPtI);
}
return triFace(faceBasePtI, facePtI, faceOtherPtI);
}
inline Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const
{
const pointField& meshPoints = mesh.points();

2
src/functionObjects/field/Make/files
View File

@ -25,12 +25,10 @@ streamLine/streamLineBase.C
streamLine/streamLineParticle.C
streamLine/streamLineParticleCloud.C
/*
wallBoundedStreamLine/wallBoundedStreamLine.C
wallBoundedStreamLine/wallBoundedStreamLineParticle.C
wallBoundedStreamLine/wallBoundedStreamLineParticleCloud.C
wallBoundedStreamLine/wallBoundedParticle.C
*/
surfaceInterpolate/surfaceInterpolate.C

9
src/functionObjects/field/nearWallFields/nearWallFields.C
View File

@ -111,11 +111,12 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
tetPti = (startInfo.index()+1) % f.size();
start = startInfo.hitPoint();
//// Uncomment below to shift slightly in:
//tetIndices tet(celli, meshFacei, tetPti, mesh_);
//start =
// (1.0-1e-6)*startInfo.hitPoint()
// + 1e-6*tet.tet(mesh_).centre();
tetIndices tet(celli, meshFacei, tetPti);
start =
(1.0 - 1e-6)*startInfo.hitPoint()
+ 1e-6*tet.tet(mesh_).centre();
}
else
{

137
src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
View File

@ -35,57 +35,6 @@ const std::size_t Foam::wallBoundedParticle::sizeofFields_
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::tetIndices Foam::wallBoundedParticle::currentTetIndices() const
{
// Replacement for particle::currentTetIndices that avoids error
// upon invalid tetBasePtIs
const faceList& pFaces = mesh_.faces();
const labelList& pOwner = mesh_.faceOwner();
const Foam::face& f = pFaces[tetFacei_];
bool own = (pOwner[tetFacei_] == celli_);
label faceBasePtI = mesh_.tetBasePtIs()[tetFacei_];
if (faceBasePtI == -1)
{
//WarningInFunction
// << "No base point for face " << tetFacei_ << ", " << f
// << ", produces a decomposition that has a minimum "
// << "volume greater than tolerance."
// << endl;
faceBasePtI = 0;
}
label facePtI = (tetPti_ + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label facePtAI;
label facePtBI;
if (own)
{
facePtAI = facePtI;
facePtBI = otherFacePtI;
}
else
{
facePtAI = otherFacePtI;
facePtBI = facePtI;
}
return tetIndices
(
celli_,
tetFacei_,
faceBasePtI,
facePtAI,
facePtBI,
tetPti_
);
}
Foam::edge Foam::wallBoundedParticle::currentEdge() const
{
if ((meshEdgeStart_ != -1) == (diagEdge_ != -1))
@ -99,7 +48,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
<< abort(FatalError);
}
const Foam::face& f = mesh_.faces()[tetFace()];
const Foam::face& f = mesh().faces()[tetFace()];
if (meshEdgeStart_ != -1)
{
@ -107,7 +56,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
}
else
{
label faceBasePti = mesh_.tetBasePtIs()[tetFace()];
label faceBasePti = mesh().tetBasePtIs()[tetFace()];
if (faceBasePti == -1)
{
//FatalErrorInFunction
@ -120,7 +69,7 @@ Foam::edge Foam::wallBoundedParticle::currentEdge() const
faceBasePti = 0;
}
label diagPti = (faceBasePti+diagEdge_)%f.size();
label diagPti = (faceBasePti + diagEdge_)%f.size();
return edge(f[faceBasePti], f[diagPti]);
}
@ -135,26 +84,24 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace
const edge& e
)
{
const faceList& pFaces = mesh_.faces();
const cellList& pCells = mesh_.cells();
const faceList& pFaces = mesh().faces();
const cellList& pCells = mesh().cells();
const Foam::face& f = pFaces[tetFacei];
const Foam::cell& thisCell = pCells[celli];
forAll(thisCell, cFI)
for (const label facei : thisCell)
{
// Loop over all other faces of this cell and
// find the one that shares this edge
label fI = thisCell[cFI];
if (tetFacei == fI)
if (tetFacei == facei)
{
continue;
}
const Foam::face& otherFace = pFaces[fI];
const Foam::face& otherFace = pFaces[facei];
label edDir = otherFace.edgeDirection(e);
@ -162,7 +109,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace
{
continue;
}
else if (f == pFaces[fI])
else if (f == pFaces[facei])
{
// This is a necessary condition if using duplicate baffles
// (so coincident faces). We need to make sure we don't cross into
@ -174,7 +121,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace
else
{
//Found edge on other face
tetFacei = fI;
tetFacei = facei;
label eIndex = -1;
@ -192,7 +139,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace
eIndex = findIndex(otherFace, e.end());
}
label tetBasePtI = mesh_.tetBasePtIs()[fI];
label tetBasePtI = mesh().tetBasePtIs()[facei];
if (tetBasePtI == -1)
{
@ -246,7 +193,7 @@ void Foam::wallBoundedParticle::crossEdgeConnectedFace(const edge& meshEdge)
face() = tetFace();
// And adapt meshEdgeStart_.
const Foam::face& f = mesh_.faces()[tetFace()];
const Foam::face& f = mesh().faces()[tetFace()];
label fp = findIndex(f, meshEdge[0]);
if (f.nextLabel(fp) == meshEdge[1])
@ -303,7 +250,7 @@ void Foam::wallBoundedParticle::crossDiagonalEdge()
<< "meshEdgeStart_:" << meshEdgeStart_ << abort(FatalError);
}
const Foam::face& f = mesh_.faces()[tetFace()];
const Foam::face& f = mesh().faces()[tetFace()];
// tetPti starts from 1, goes up to f.size()-2
@ -340,7 +287,7 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
)
{
// Track p from position to endPosition
const triFace tri(currentTetIndices().faceTriIs(mesh_));
const triFace tri(currentTetIndices().faceTriIs(mesh(), false));
// Check which edge intersects the trajectory.
// Project trajectory onto triangle
@ -358,8 +305,8 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
{
label j = tri.fcIndex(i);
const point& pt0 = mesh_.points()[tri[i]];
const point& pt1 = mesh_.points()[tri[j]];
const point& pt0 = mesh().points()[tri[i]];
const point& pt1 = mesh().points()[tri[j]];
if (edge(tri[i], tri[j]) == currentE)
{
@ -368,20 +315,20 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
}
// Outwards pointing normal
vector edgeNormal = (pt1-pt0)^n;
vector edgeNormal = (pt1 - pt0)^n;
edgeNormal /= mag(edgeNormal)+VSMALL;
edgeNormal /= mag(edgeNormal) + VSMALL;
// Determine whether position and end point on either side of edge.
scalar sEnd = (endPosition-pt0)&edgeNormal;
scalar sEnd = (endPosition - pt0)&edgeNormal;
if (sEnd >= 0)
{
// endPos is outside triangle. position() should always be
// endPos is outside triangle. localPosition_ should always be
// inside.
scalar sStart = (position()-pt0)&edgeNormal;
if (mag(sEnd-sStart) > VSMALL)
scalar sStart = (localPosition_ - pt0)&edgeNormal;
if (mag(sEnd - sStart) > VSMALL)
{
scalar s = sStart/(sStart-sEnd);
scalar s = sStart/(sStart - sEnd);
if (s >= 0 && s < minS)
{
@ -394,19 +341,19 @@ Foam::scalar Foam::wallBoundedParticle::trackFaceTri
if (minEdgei != -1)
{
position() += minS*(endPosition-position());
localPosition_ += minS*(endPosition - localPosition_);
}
else
{
// Did not hit any edge so tracked to the end position. Set position
// without any calculation to avoid truncation errors.
position() = endPosition;
localPosition_ = endPosition;
minS = 1.0;
}
// Project position() back onto plane of triangle
const point& triPt = mesh_.points()[tri[0]];
position() -= ((position()-triPt)&n)*n;
// Project localPosition_ back onto plane of triangle
const point& triPt = mesh().points()[tri[0]];
localPosition_ -= ((localPosition_ - triPt)&n)*n;
return minS;
}
@ -418,7 +365,7 @@ bool Foam::wallBoundedParticle::isTriAlongTrack
const point& endPosition
) const
{
const triFace triVerts(currentTetIndices().faceTriIs(mesh_));
const triFace triVerts(currentTetIndices().faceTriIs(mesh(), false));
const edge currentE = currentEdge();
if
@ -435,13 +382,13 @@ bool Foam::wallBoundedParticle::isTriAlongTrack
}
const vector dir = endPosition-position();
const vector dir = endPosition - localPosition_;
forAll(triVerts, i)
{
label j = triVerts.fcIndex(i);
const point& pt0 = mesh_.points()[triVerts[i]];
const point& pt1 = mesh_.points()[triVerts[j]];
const point& pt0 = mesh().points()[triVerts[i]];
const point& pt1 = mesh().points()[triVerts[j]];
if (edge(triVerts[i], triVerts[j]) == currentE)
{
@ -462,7 +409,7 @@ bool Foam::wallBoundedParticle::isTriAlongTrack
Foam::wallBoundedParticle::wallBoundedParticle
(
const polyMesh& mesh,
const vector& position,
const point& position,
const label celli,
const label tetFacei,
const label tetPti,
@ -470,7 +417,9 @@ Foam::wallBoundedParticle::wallBoundedParticle
const label diagEdge
)
:
particle(mesh, position, celli, tetFacei, tetPti),
// particle(mesh, barycentric(1, 0, 0, 0), celli, tetFacei, tetPti),
particle(mesh, position, celli, tetFacei, tetPti, false),
localPosition_(position),
meshEdgeStart_(meshEdgeStart),
diagEdge_(diagEdge)
{}
@ -490,11 +439,11 @@ Foam::wallBoundedParticle::wallBoundedParticle
{
if (is.format() == IOstream::ASCII)
{
is >> meshEdgeStart_ >> diagEdge_;
is >> localPosition_ >> meshEdgeStart_ >> diagEdge_;
}
else
{
is.read(reinterpret_cast<char*>(&meshEdgeStart_), sizeofFields_);
is.read(reinterpret_cast<char*>(&localPosition_), sizeofFields_);
}
}
@ -508,6 +457,7 @@ Foam::wallBoundedParticle::wallBoundedParticle
)
:
particle(p),
localPosition_(p.localPosition_),
meshEdgeStart_(p.meshEdgeStart_),
diagEdge_(p.diagEdge_)
{}
@ -524,15 +474,16 @@ Foam::Ostream& Foam::operator<<
if (os.format() == IOstream::ASCII)
{
os << static_cast<const particle&>(p)
<< token::SPACE << p.meshEdgeStart_
<< token::SPACE << p.diagEdge_;
<< token::SPACE << p.localPosition_
<< token::SPACE << p.meshEdgeStart_
<< token::SPACE << p.diagEdge_;
}
else
{
os << static_cast<const particle&>(p);
os.write
(
reinterpret_cast<const char*>(&p.meshEdgeStart_),
reinterpret_cast<const char*>(&p.localPosition_),
wallBoundedParticle::sizeofFields_
);
}
@ -568,7 +519,7 @@ Foam::Ostream& Foam::operator<<
<< " l 2 4" << nl
<< " l 3 4" << nl;
os << " ";
meshTools::writeOBJ(os, p.position());
meshTools::writeOBJ(os, p.localPosition_);
return os;
}

137
src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
View File

@ -71,10 +71,9 @@ class wallBoundedParticle
public:
//- Class used to pass tracking data to the trackToFace function
template<class CloudType>
class TrackingData
class trackingData
:
public particle::TrackingData<CloudType>
public particle::trackingData
{
public:
@ -83,16 +82,14 @@ public:
// Constructors
inline TrackingData
template <class TrackCloudType>
inline trackingData
(
CloudType& cloud,
const TrackCloudType& cloud,
const PackedBoolList& isWallPatch
)
:
particle::TrackingData<CloudType>
(
cloud
),
particle::trackingData(cloud),
isWallPatch_(isWallPatch)
{}
};
@ -102,6 +99,10 @@ protected:
// Protected data
//- Particle position is updated locally as opposed to via track
// functions of the base Foam::particle class
point localPosition_;
//- Particle is on mesh edge:
// const face& f = mesh.faces()[tetFace()]
// const edge e(f[meshEdgeStart_], f.nextLabel(meshEdgeStart_));
@ -123,10 +124,6 @@ protected:
//- Construct current edge
edge currentEdge() const;
//- Replacement for particle::currentTetIndices() that avoids bombing
// out on invalid tet decomposition (tetBasePtIs = -1)
tetIndices currentTetIndices() const;
//- Replacement for particle::crossEdgeConnectedFace that avoids bombing
// out on invalid tet decomposition (tetBasePtIs = -1)
void crossEdgeConnectedFace
@ -150,84 +147,6 @@ protected:
bool isTriAlongTrack(const vector& n, const point& endPosition) const;
// Patch interactions
//- Do all patch interaction
template<class TrackData>
void patchInteraction(TrackData& td, const scalar trackFraction);
//- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions
template<class TrackData>
bool hitPatch
(
const polyPatch&,
TrackData& td,
const label patchi,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a wedge
template<class TrackData>
void hitWedgePatch
(
const wedgePolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
// symmetry plane
template<class TrackData>
void hitSymmetryPlanePatch
(
const symmetryPlanePolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
// symmetry patch
template<class TrackData>
void hitSymmetryPatch
(
const symmetryPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a cyclic
template<class TrackData>
void hitCyclicPatch
(
const cyclicPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a
//- processorPatch
template<class TrackData>
void hitProcessorPatch
(
const processorPolyPatch&,
TrackData& td
);
//- Overridable function to handle the particle hitting a wallPatch
template<class TrackData>
void hitWallPatch
(
const wallPolyPatch&,
TrackData& td,
const tetIndices&
);
//- Overridable function to handle the particle hitting a polyPatch
template<class TrackData>
void hitPatch
(
const polyPatch&,
TrackData& td
);
public:
// Constructors
@ -236,7 +155,7 @@ public:
wallBoundedParticle
(
const polyMesh& c,
const vector& position,
const point& position,
const label celli,
const label tetFacei,
const label tetPti,
@ -308,14 +227,44 @@ public:
// Track
//- Equivalent of trackToFace
template<class TrackData>
template<class TrackCloudType>
scalar trackToEdge
(
TrackData& td,
TrackCloudType& cloud,
trackingData& td,
const vector& endPosition
);
// Patch interactions
//- Do all patch interaction
template<class TrackCloudType>
void patchInteraction
(
TrackCloudType& cloud,
trackingData& td,
const scalar trackFraction
);
//- Overridable function to handle the particle hitting a
//- processorPatch
template<class TrackCloudType>
void hitProcessorPatch
(
TrackCloudType& cloud,
trackingData& td
);
//- Overridable function to handle the particle hitting a wallPatch
template<class TrackCloudType>
void hitWallPatch
(
TrackCloudType& cloud,
trackingData& td
);
// Info
//- Return info proxy.

332
src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
View File

@ -27,97 +27,44 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class TrackData>
template<class TrackCloudType>
void Foam::wallBoundedParticle::patchInteraction
(
TrackData& td,
TrackCloudType& cloud,
trackingData& td,
const scalar trackFraction
)
{
typedef typename TrackData::cloudType::particleType particleType;
typename TrackCloudType::particleType& p =
static_cast<typename TrackCloudType::particleType&>(*this);
typename TrackCloudType::particleType::trackingData& ttd =
static_cast<typename TrackCloudType::particleType::trackingData&>(td);
particleType& p = static_cast<particleType&>(*this);
p.hitFace(td);
if (!internalFace(facei_))
if (!mesh().isInternalFace(face()))
{
label origFacei = facei_;
label patchi = patch(facei_);
label origFacei = face();
label patchi = patch();
// No action taken for tetPti_ for tetFacei_ here, handled by
// patch interaction call or later during processor transfer.
// Dummy tet indices. What to do here?
tetIndices faceHitTetIs;
if
(
!p.hitPatch
(
mesh_.boundaryMesh()[patchi],
td,
patchi,
trackFraction,
faceHitTetIs
)
)
// Did patch interaction model switch patches?
// Note: recalculate meshEdgeStart_, diagEdge_!
if (face() != origFacei)
{
// Did patch interaction model switch patches?
// Note: recalculate meshEdgeStart_, diagEdge_!
if (facei_ != origFacei)
{
patchi = patch(facei_);
}
patchi = patch();
}
const polyPatch& patch = mesh_.boundaryMesh()[patchi];
const polyPatch& patch = mesh().boundaryMesh()[patchi];
if (isA<wedgePolyPatch>(patch))
{
p.hitWedgePatch
(
static_cast<const wedgePolyPatch&>(patch), td
);
}
else if (isA<symmetryPlanePolyPatch>(patch))
{
p.hitSymmetryPlanePatch
(
static_cast<const symmetryPlanePolyPatch&>(patch), td
);
}
else if (isA<symmetryPolyPatch>(patch))
{
p.hitSymmetryPatch
(
static_cast<const symmetryPolyPatch&>(patch), td
);
}
else if (isA<cyclicPolyPatch>(patch))
{
p.hitCyclicPatch
(
static_cast<const cyclicPolyPatch&>(patch), td
);
}
else if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch
(
static_cast<const processorPolyPatch&>(patch), td
);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch
(
static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
);
}
else
{
p.hitPatch(patch, td);
}
if (isA<processorPolyPatch>(patch))
{
p.hitProcessorPatch(cloud, ttd);
}
else if (isA<wallPolyPatch>(patch))
{
p.hitWallPatch(cloud, ttd);
}
else
{
td.keepParticle = false;
}
}
}
@ -125,10 +72,11 @@ void Foam::wallBoundedParticle::patchInteraction
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class TrackData>
template<class TrackCloudType>
Foam::scalar Foam::wallBoundedParticle::trackToEdge
(
TrackData& td,
TrackCloudType& cloud,
trackingData& td,
const vector& endPosition
)
{
@ -147,7 +95,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
//checkInside();
//checkOnTriangle(position());
//checkOnTriangle(localPosition_);
//if (meshEdgeStart_ != -1 || diagEdge_ != -1)
//{
// checkOnEdge();
@ -163,39 +111,36 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// If internal face check whether to go to neighbour cell or just
// check to the other internal tet on the edge.
if (mesh_.isInternalFace(tetFace()))
if (mesh().isInternalFace(tetFace()))
{
label nbrCelli =
(
celli_ == mesh_.faceOwner()[facei_]
? mesh_.faceNeighbour()[facei_]
: mesh_.faceOwner()[facei_]
cell() == mesh().faceOwner()[face()]
? mesh().faceNeighbour()[face()]
: mesh().faceOwner()[face()]
);
// Check angle to nbrCell tet. Is it in the direction of the
// endposition? i.e. since volume of nbr tet is positive the
// tracking direction should be into the tet.
tetIndices nbrTi(nbrCelli, tetFacei_, tetPti_, mesh_);
tetIndices nbrTi(nbrCelli, tetFace(), tetPt());
bool posVol = (nbrTi.tet(mesh_).mag() > 0);
const bool posVol = (nbrTi.tet(mesh()).mag() > 0);
const vector path(endPosition - localPosition_);
if
(
posVol
== ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
)
if (posVol == ((nbrTi.faceTri(mesh()).normal() & path) < 0))
{
// Change into nbrCell. No need to change tetFace, tetPt.
//Pout<< " crossed from cell:" << celli_
// << " into " << nbrCelli << endl;
celli_ = nbrCelli;
patchInteraction(td, trackFraction);
cell() = nbrCelli;
patchInteraction(cloud, td, trackFraction);
}
else
{
// Walk to other face on edge. Changes tetFace, tetPt but not
// cell.
crossEdgeConnectedFace(meshEdge);
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
}
else
@ -203,7 +148,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// Walk to other face on edge. This might give loop since
// particle should have been removed?
crossEdgeConnectedFace(meshEdge);
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
}
else
@ -211,41 +156,42 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// We're inside a tet on the wall. Check if the current tet is
// the one to cross. If not we cross into the neighbouring triangle.
if (mesh_.isInternalFace(tetFace()))
if (mesh().isInternalFace(tetFace()))
{
FatalErrorInFunction
<< "Can only track on boundary faces."
<< " Face:" << tetFace()
<< " at:" << mesh_.faceCentres()[tetFace()]
<< " at:" << mesh().faceCentres()[tetFace()]
<< abort(FatalError);
}
const triFace tri(currentTetIndices().faceTriIs(mesh_));
vector n = tri.normal(mesh_.points());
const triFace tri(currentTetIndices().faceTriIs(mesh(), false));
vector n = tri.normal(mesh().points());
n /= mag(n);
point projectedEndPosition = endPosition;
const bool posVol = (currentTetIndices().tet(mesh_).mag() > 0);
const bool posVol = (currentTetIndices().tet(mesh()).mag() > 0);
if (!posVol)
{
// Negative tet volume. Track back by setting the end point
projectedEndPosition = position() - (endPosition-position());
projectedEndPosition =
localPosition_ - (endPosition - localPosition_);
// Make sure to use a large enough vector to cross the negative
// face. Bit overkill.
const vector d(endPosition-position());
const vector d(endPosition - localPosition_);
const scalar magD(mag(d));
if (magD > ROOTVSMALL)
{
// Get overall mesh bounding box
treeBoundBox meshBb(mesh_.bounds());
treeBoundBox meshBb(mesh().bounds());
// Extend to make 3D
meshBb.inflate(ROOTSMALL);
// Create vector guaranteed to cross mesh bounds
projectedEndPosition = position()-meshBb.mag()*d/magD;
projectedEndPosition = localPosition_ - meshBb.mag()*d/magD;
// Clip to mesh bounds
point intPt;
@ -253,9 +199,9 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
bool ok = meshBb.intersects
(
projectedEndPosition,
position()-projectedEndPosition,
localPosition_ - projectedEndPosition,
projectedEndPosition,
position(),
localPosition_,
intPt,
intPtBits
);
@ -269,8 +215,8 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// Remove normal component
{
const point& basePt = mesh_.points()[tri[0]];
projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
const point& basePt = mesh().points()[tri[0]];
projectedEndPosition -= ((projectedEndPosition - basePt)&n)*n;
}
@ -304,7 +250,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
}
const tetIndices ti(currentTetIndices());
const triFace trif(ti.triIs(mesh(), false));
// Triangle (faceTriIs) gets constructed from
// f[faceBasePtI_],
// f[facePtAI_],
@ -315,21 +261,21 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
// 2 : edge between facePtBI_ and faceBasePtI_
const Foam::face& f = mesh_.faces()[ti.face()];
const label fp0 = ti.faceBasePt();
const Foam::face& f = mesh().faces()[ti.face()];
const label fp0 = trif[0];
if (triEdgei == 0)
{
if (ti.facePtA() == f.fcIndex(fp0))
if (trif[1] == f.fcIndex(fp0))
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = fp0;
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else if (ti.facePtA() == f.rcIndex(fp0))
else if (trif[1] == f.rcIndex(fp0))
{
//Note: should not happen since boundary face so owner
//Pout<< "Real edge." << endl;
@ -340,12 +286,12 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
meshEdgeStart_ = f.rcIndex(fp0);
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else
{
// Get index of triangle on other side of edge.
diagEdge_ = ti.facePtA()-fp0;
diagEdge_ = trif[1] - fp0;
if (diagEdge_ < 0)
{
diagEdge_ += f.size();
@ -359,40 +305,39 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = ti.facePtA();
meshEdgeStart_ = trif[1];
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else // if (triEdgei == 2)
{
if (ti.facePtB() == f.rcIndex(fp0))
if (trif[2] == f.rcIndex(fp0))
{
//Pout<< "Real edge." << endl;
diagEdge_ = -1;
meshEdgeStart_ = ti.facePtB();
meshEdgeStart_ = trif[2];
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else if (ti.facePtB() == f.fcIndex(fp0))
else if (trif[2] == f.fcIndex(fp0))
{
//Note: should not happen since boundary face so owner
//Pout<< "Real edge." << endl;
FatalErrorInFunction
<< abort(FatalError);
FatalErrorInFunction << abort(FatalError);
diagEdge_ = -1;
meshEdgeStart_ = fp0;
//checkOnEdge();
crossEdgeConnectedFace(currentEdge());
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else
{
//Pout<< "Triangle edge." << endl;
// Get index of triangle on other side of edge.
diagEdge_ = ti.facePtB()-fp0;
diagEdge_ = trif[2] - fp0;
if (diagEdge_ < 0)
{
diagEdge_ += f.size();
@ -412,7 +357,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
// Particle is on mesh edge so change into other face on cell
crossEdgeConnectedFace(currentEdge());
//checkOnEdge();
patchInteraction(td, trackFraction);
patchInteraction(cloud, td, trackFraction);
}
else
{
@ -430,74 +375,11 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
}
template<class TrackData>
bool Foam::wallBoundedParticle::hitPatch
(
const polyPatch&,
TrackData& td,
const label patchi,
const scalar trackFraction,
const tetIndices& tetIs
)
{
// Disable generic patch interaction
return false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitWedgePatch
(
const wedgePolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitSymmetryPlanePatch
(
const symmetryPlanePolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitSymmetryPatch
(
const symmetryPolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
void Foam::wallBoundedParticle::hitCyclicPatch
(
const cyclicPolyPatch& pp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class TrackData>
template<class TrackCloudType>
void Foam::wallBoundedParticle::hitProcessorPatch
(
const processorPolyPatch& pp,
TrackData& td
TrackCloudType& cloud,
trackingData& td
)
{
// Switch particle
@ -508,44 +390,31 @@ void Foam::wallBoundedParticle::hitProcessorPatch
// on the other side between 2 and 3 so edgeStart_ should be
// f.size()-edgeStart_-1.
const Foam::face& f = mesh_.faces()[face()];
const Foam::face& f = mesh().faces()[face()];
if (meshEdgeStart_ != -1)
{
meshEdgeStart_ = f.size()-meshEdgeStart_-1;
meshEdgeStart_ = f.size() - meshEdgeStart_-1;
}
else
{
// diagEdge_ is relative to faceBasePt
diagEdge_ = f.size()-diagEdge_;
diagEdge_ = f.size() - diagEdge_;
}
}
template<class TrackData>
template<class TrackCloudType>
void Foam::wallBoundedParticle::hitWallPatch
(
const wallPolyPatch& wpp,
TrackData& td,
const tetIndices&
TrackCloudType& cloud,
trackingData& td
)
{}
template<class TrackData>
void Foam::wallBoundedParticle::hitPatch
(
const polyPatch& wpp,
TrackData& td
)
{
// Remove particle
td.keepParticle = false;
}
template<class CloudType>
void Foam::wallBoundedParticle::readFields(CloudType& c)
template<class TrackCloudType>
void Foam::wallBoundedParticle::readFields(TrackCloudType& c)
{
if (!c.size())
{
@ -554,34 +423,47 @@ void Foam::wallBoundedParticle::readFields(CloudType& c)
particle::readFields(c);
IOField<point> localPosition
(
c.fieldIOobject("position", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, localPosition);
IOField<label> meshEdgeStart
(
c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, meshEdgeStart);
IOField<label> diagEdge
(
c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
c.fieldIOobject("diagEdge", IOobject::MUST_READ)
);
c.checkFieldIOobject(c, diagEdge);
label i = 0;
forAllIter(typename CloudType, c, iter)
forAllIters(c, iter)
{
iter().localPosition_ = localPosition[i];
iter().meshEdgeStart_ = meshEdgeStart[i];
iter().diagEdge_ = diagEdge[i];
i++;
++i;
}
}
template<class CloudType>
void Foam::wallBoundedParticle::writeFields(const CloudType& c)
template<class TrackCloudType>
void Foam::wallBoundedParticle::writeFields(const TrackCloudType& c)
{
particle::writeFields(c);
label np = c.size();
IOField<point> localPosition
(
c.fieldIOobject("position", IOobject::NO_READ),
np
);
IOField<label> meshEdgeStart
(
c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
@ -594,13 +476,15 @@ void Foam::wallBoundedParticle::writeFields(const CloudType& c)
);
label i = 0;
forAllConstIter(typename CloudType, c, iter)
forAllConstIters(c, iter)
{
localPosition[i] = iter().localPosition_;
meshEdgeStart[i] = iter().meshEdgeStart_;
diagEdge[i] = iter().diagEdge_;
i++;
++i;
}
localPosition.write();
meshEdgeStart.write();
diagEdge.write();
}

27
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -62,10 +62,8 @@ Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
label minTetPti = -1;
scalar minDistSqr = sqr(GREAT);
forAll(cFaces, cFacei)
for (label facei : cFaces)
{
label facei = cFaces[cFacei];
if (isWallPatch[facei])
{
const face& f = mesh_.faces()[facei];
@ -98,8 +96,7 @@ Foam::tetIndices Foam::functionObjects::wallBoundedStreamLine::findNearestTet
(
celli,
minFacei,
minTetPti,
mesh_
minTetPti
);
}
@ -226,7 +223,7 @@ void Foam::functionObjects::wallBoundedStreamLine::track()
const scalar trackTime = Foam::sqrt(GREAT);
// Track
particles.move(td, trackTime);
particles.move(particles, td, trackTime);
}
@ -265,6 +262,7 @@ bool Foam::functionObjects::wallBoundedStreamLine::read(const dictionary& dict)
{
// 1. Positive volume decomposition tets
faceSet faces(mesh_, "lowQualityTetFaces", mesh_.nFaces()/100+1);
if
(
polyMeshTetDecomposition::checkFaceTets
@ -287,21 +285,18 @@ bool Foam::functionObjects::wallBoundedStreamLine::read(const dictionary& dict)
// 2. All edges on a cell having two faces
EdgeMap<label> numFacesPerEdge;
forAll(mesh_.cells(), celli)
for (const cell& cFaces : mesh_.cells())
{
const cell& cFaces = mesh_.cells()[celli];
numFacesPerEdge.clear();
forAll(cFaces, cFacei)
for (const label facei : cFaces)
{
label facei = cFaces[cFacei];
const face& f = mesh_.faces()[facei];
forAll(f, fp)
{
const edge e(f[fp], f.nextLabel(fp));
EdgeMap<label>::iterator eFnd =
numFacesPerEdge.find(e);
EdgeMap<label>::iterator eFnd = numFacesPerEdge.find(e);
if (eFnd != numFacesPerEdge.end())
{
eFnd()++;
@ -313,12 +308,12 @@ bool Foam::functionObjects::wallBoundedStreamLine::read(const dictionary& dict)
}
}
forAllConstIter(EdgeMap<label>, numFacesPerEdge, iter)
forAllConstIters(numFacesPerEdge, iter)
{
if (iter() != 2)
{
FatalErrorInFunction
<< "problem cell:" << celli
<< "problem cell:" << cFaces
<< abort(FatalError);
}
}

153
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
View File

@ -42,20 +42,14 @@ Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
<< "Cell:" << celli << abort(FatalError);
}
const tetIndices ti = currentTetIndices();
const vector U = td.vvInterp_[td.UIndex_].interpolate
(
position,
ti, //celli,
facei
);
const vector U =
td.vvInterp_[td.UIndex_].interpolate(position, celli, facei);
// Check if at different position
if
(
!sampledPositions_.size()
|| magSqr(sampledPositions_.last()-position) > Foam::sqr(SMALL)
|| magSqr(sampledPositions_.last() - position) > Foam::sqr(SMALL)
)
{
// Store the location
@ -67,12 +61,7 @@ Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
{
sampledScalars_[scalari].append
(
td.vsInterp_[scalari].interpolate
(
position,
ti, //celli,
facei
)
td.vsInterp_[scalari].interpolate(position, celli, facei)
);
}
@ -87,12 +76,8 @@ Foam::vector Foam::wallBoundedStreamLineParticle::interpolateFields
}
else
{
positionU = td.vvInterp_[vectori].interpolate
(
position,
ti, //celli,
facei
);
positionU =
td.vvInterp_[vectori].interpolate(position, celli, facei);
}
sampledVectors_[vectori].append(positionU);
}
@ -107,7 +92,7 @@ Foam::vector Foam::wallBoundedStreamLineParticle::sample
trackingData& td
)
{
vector U = interpolateFields(td, position(), cell(), tetFace());
vector U = interpolateFields(td, localPosition_, cell(), face());
if (!td.trackForward_)
{
@ -134,7 +119,7 @@ Foam::vector Foam::wallBoundedStreamLineParticle::sample
Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
(
const polyMesh& mesh,
const vector& position,
const point& position,
const label celli,
const label tetFacei,
const label tetPti,
@ -206,128 +191,6 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::wallBoundedStreamLineParticle::move
(
trackingData& td,
const scalar trackTime
)
{
wallBoundedStreamLineParticle& p = static_cast
<
wallBoundedStreamLineParticle&
>(*this);
// Check position is inside tet
//checkInside();
td.switchProcessor = false;
td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar maxDt = mesh_.bounds().mag();
while
(
td.keepParticle
&& !td.switchProcessor
&& lifeTime_ > 0
)
{
// set the lagrangian time-step
scalar dt = maxDt;
--lifeTime_;
// Get sampled velocity and fields. Store if position changed.
vector U = sample(td);
// !user parameter!
if (dt < SMALL)
{
// Force removal
lifeTime_ = 0;
break;
}
if (td.trackLength_ < GREAT)
{
dt = td.trackLength_;
}
scalar fraction = trackToEdge(td, position() + dt*U);
dt *= fraction;
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
if (tEnd <= ROOTVSMALL)
{
// Force removal
lifeTime_ = 0;
}
}
if (!td.keepParticle || lifeTime_ == 0)
{
if (lifeTime_ == 0)
{
if (debug)
{
Pout<< "wallBoundedStreamLineParticle :"
<< " Removing stagnant particle:"
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
td.keepParticle = false;
}
else
{
// Normal exit. Store last position and fields
sample(td);
if (debug)
{
Pout<< "wallBoundedStreamLineParticle : Removing particle:"
<< p.position()
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
}
// Transfer particle data into trackingData.
{
//td.allPositions_.append(sampledPositions_);
td.allPositions_.append(vectorList());
vectorList& top = td.allPositions_.last();
top.transfer(sampledPositions_);
}
forAll(sampledScalars_, i)
{
//td.allScalars_[i].append(sampledScalars_[i]);
td.allScalars_[i].append(scalarList());
scalarList& top = td.allScalars_[i].last();
top.transfer(sampledScalars_[i]);
}
forAll(sampledVectors_, i)
{
//td.allVectors_[i].append(sampledVectors_[i]);
td.allVectors_[i].append(vectorList());
vectorList& top = td.allVectors_[i].last();
top.transfer(sampledVectors_[i]);
}
}
return td.keepParticle;
}
void Foam::wallBoundedStreamLineParticle::readFields
(
Cloud<wallBoundedStreamLineParticle>& c

33
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -58,7 +58,7 @@ Ostream& operator<<(Ostream&, const wallBoundedStreamLineParticle&);
/*---------------------------------------------------------------------------*\
Class wallBoundedStreamLineParticle Declaration
Class wallBoundedStreamLineParticle Declaration
\*---------------------------------------------------------------------------*/
class wallBoundedStreamLineParticle
@ -71,10 +71,7 @@ public:
//- Class used to pass tracking data to the trackToEdge function
class trackingData
:
public wallBoundedParticle::TrackingData
<
Cloud<wallBoundedStreamLineParticle>
>
public wallBoundedParticle::trackingData
{
public:
@ -93,9 +90,10 @@ public:
// Constructors
template <class TrackCloudType>
trackingData
(
Cloud<wallBoundedStreamLineParticle>& cloud,
TrackCloudType& cloud,
const PtrList<interpolation<scalar>>& vsInterp,
const PtrList<interpolation<vector>>& vvInterp,
const label UIndex,
@ -108,10 +106,7 @@ public:
List<DynamicList<vectorList>>& allVectors
)
:
wallBoundedParticle::TrackingData
<
Cloud<wallBoundedStreamLineParticle>
>
wallBoundedParticle::trackingData
(
cloud,
isWallPatch
@ -167,7 +162,7 @@ public:
wallBoundedStreamLineParticle
(
const polyMesh& c,
const vector& position,
const point& position,
const label celli,
const label tetFacei,
const label tetPti,
@ -225,7 +220,13 @@ public:
// Tracking
//- Track all particles to their end point
bool move(trackingData&, const scalar trackTime);
template<class TrackCloudType>
bool move
(
TrackCloudType& cloud,
trackingData& td,
const scalar trackTime
);
// I-O
@ -256,6 +257,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "wallBoundedStreamLineParticleTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

150
src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleTemplates.C Normal file
View File

@ -0,0 +1,150 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class TrackCloudType>
bool Foam::wallBoundedStreamLineParticle::move
(
TrackCloudType& cloud,
trackingData& td,
const scalar trackTime
)
{
typename TrackCloudType::particleType& p =
static_cast<typename TrackCloudType::particleType&>(*this);
// Check position is inside tet
//checkInside();
td.switchProcessor = false;
td.keepParticle = true;
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar maxDt = mesh().bounds().mag();
while
(
td.keepParticle
&& !td.switchProcessor
&& lifeTime_ > 0
)
{
// set the lagrangian time-step
scalar dt = maxDt;
--lifeTime_;
// Get sampled velocity and fields. Store if position changed.
vector U = p.sample(td);
// !user parameter!
if (dt < SMALL)
{
// Force removal
lifeTime_ = 0;
break;
}
if (td.trackLength_ < GREAT)
{
dt = td.trackLength_;
}
scalar fraction = trackToEdge(cloud, td, localPosition_ + dt*U);
dt *= fraction;
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
if (tEnd <= ROOTVSMALL)
{
// Force removal
lifeTime_ = 0;
}
}
if (!td.keepParticle || lifeTime_ == 0)
{
if (lifeTime_ == 0)
{
if (debug)
{
Pout<< "wallBoundedStreamLineParticle :"
<< " Removing stagnant particle:"
<< localPosition_
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
td.keepParticle = false;
}
else
{
// Normal exit. Store last position and fields
sample(td);
if (debug)
{
Pout<< "wallBoundedStreamLineParticle : Removing particle:"
<< localPosition_
<< " sampled positions:" << sampledPositions_.size()
<< endl;
}
}
// Transfer particle data into trackingData.
{
//td.allPositions_.append(sampledPositions_);
td.allPositions_.append(vectorList());
vectorList& top = td.allPositions_.last();
top.transfer(sampledPositions_);
}
forAll(sampledScalars_, i)
{
//td.allScalars_[i].append(sampledScalars_[i]);
td.allScalars_[i].append(scalarList());
scalarList& top = td.allScalars_[i].last();
top.transfer(sampledScalars_[i]);
}
forAll(sampledVectors_, i)
{
//td.allVectors_[i].append(sampledVectors_[i]);
td.allVectors_[i].append(vectorList());
vectorList& top = td.allVectors_[i].last();
top.transfer(sampledVectors_[i]);
}
}
return td.keepParticle;
}
// ************************************************************************* //

36
src/lagrangian/basic/particle/particle.C
View File

@ -588,7 +588,7 @@ Foam::particle::particle
)
:
mesh_(mesh),
coordinates_(- VGREAT, - VGREAT, - VGREAT, - VGREAT),
coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
celli_(celli),
tetFacei_(-1),
tetPti_(-1),
@ -608,6 +608,40 @@ Foam::particle::particle
}
Foam::particle::particle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const label tetFacei,
const label tetPti,
bool doLocate
)
:
mesh_(mesh),
coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
celli_(celli),
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
stepFraction_(0.0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
if (doLocate)
{
locate
(
position,
nullptr,
celli,
false,
"Particle initialised with a location outside of the mesh."
);
}
}
Foam::particle::particle(const particle& p)
:
mesh_(p.mesh_),

21
src/lagrangian/basic/particle/particle.H
View File

@ -370,6 +370,18 @@ public:
const label celli
);
//- Construct from components
particle
(
const polyMesh& mesh,
const vector& position,
const label celli,
const label tetFacei,
const label tetPti,
const bool doLocate = true
);
//- Construct from Istream
particle
(
@ -438,12 +450,21 @@ public:
//- Return current tet face particle is in
inline label tetFace() const;
//- Return current tet face particle is in for manipulation
inline label& tetFace();
//- Return current tet face particle is in
inline label tetPt() const;
//- Return current tet face particle is in for manipulation
inline label& tetPt();
//- Return current face particle is on otherwise -1
inline label face() const;
//- Return current face particle is on for manipulation
inline label& face();
//- Return the fraction of time-step completed
inline scalar stepFraction() const;

18
src/lagrangian/basic/particle/particleI.H
View File

@ -72,18 +72,36 @@ inline Foam::label Foam::particle::tetFace() const
}
inline Foam::label& Foam::particle::tetFace()
{
return tetFacei_;
}
inline Foam::label Foam::particle::tetPt() const
{
return tetPti_;
}
inline Foam::label& Foam::particle::tetPt()
{
return tetPti_;
}
inline Foam::label Foam::particle::face() const
{
return facei_;
}
inline Foam::label& Foam::particle::face()
{
return facei_;
}
inline Foam::scalar Foam::particle::stepFraction() const
{
return stepFraction_;

1
tutorials/incompressible/simpleFoam/motorBike/system/controlDict
View File

@ -47,6 +47,7 @@ runTimeModifiable true;
functions
{
#include "streamLines"
#include "wallBoundedStreamLines"
#include "cuttingPlane"
#include "forceCoeffs"
#include "ensightWrite"