From c7edc26dbfa77bb27b02956920c7b3cf3333def4 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 14 Sep 2017 12:02:03 +0100
Subject: [PATCH] ENH: Reinstated the wallBoundedStreamline function object
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.
---
.../polyMeshTetDecomposition/tetIndices.H | 20 +-
.../polyMeshTetDecomposition/tetIndicesI.H | 76 +++-
src/functionObjects/field/Make/files | 2 -
.../field/nearWallFields/nearWallFields.C | 9 +-
.../wallBoundedParticle.C | 137 +++-----
.../wallBoundedParticle.H | 137 +++-----
.../wallBoundedParticleTemplates.C | 332 ++++++------------
.../wallBoundedStreamLine.C | 27 +-
.../wallBoundedStreamLineParticle.C | 153 +-------
.../wallBoundedStreamLineParticle.H | 33 +-
.../wallBoundedStreamLineParticleTemplates.C | 150 ++++++++
src/lagrangian/basic/particle/particle.C | 36 +-
src/lagrangian/basic/particle/particle.H | 21 ++
src/lagrangian/basic/particle/particleI.H | 18 +
.../simpleFoam/motorBike/system/controlDict | 1 +
15 files changed, 544 insertions(+), 608 deletions(-)
create mode 100644 src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleTemplates.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
index 8ab5d1e2e8..0ddf19294c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
@@ -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;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
index 70edef09b1..d09eeeefea 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H
@@ -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,70 @@ 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)
+ {
+ 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)
+ }
+
+ label facePtI = (tetPti_ + faceBasePtI) % f.size();
+ label faceOtherPtI = f.fcIndex(facePtI);
+
+ if (mesh.faceOwner()[face()] != cell())
+ {
+ Swap(facePtI, faceOtherPtI);
+ }
+
+ return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]);
+}
+
+
+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)
{
- 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())
@@ -93,7 +145,7 @@ inline Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
Swap(facePtI, faceOtherPtI);
}
- return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]);
+ return triFace(faceBasePtI, facePtI, faceOtherPtI);
}
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 3362b2dbca..b44dc3819f 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -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
diff --git a/src/functionObjects/field/nearWallFields/nearWallFields.C b/src/functionObjects/field/nearWallFields/nearWallFields.C
index 681956f30f..03c1aea82b 100644
--- a/src/functionObjects/field/nearWallFields/nearWallFields.C
+++ b/src/functionObjects/field/nearWallFields/nearWallFields.C
@@ -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
{
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
index 6996e1500c..8a21bba822 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
@@ -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;
}
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
index c648fc5c3c..52f006ca1c 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
@@ -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.
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
index b2febe5203..23f0a6044c 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
@@ -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_);
-
- // No action taken for tetPti_ for tetFacei_ here, handled by
- // patch interaction call or later during processor transfer.
-
+ label origFacei = face();
+ label patchi = patch();
- // 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();
}
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
index 66ff518d12..ce96d78714 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLine.C
@@ -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);
}
}
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
index 250e4a21b1..c9ab87bf5f 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -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
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
index 09b0b21229..406468e1cb 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
@@ -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
// ************************************************************************* //
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleTemplates.C
new file mode 100644
index 0000000000..0450866f9f
--- /dev/null
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticleTemplates.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 5561ac9f15..e52ab0ba19 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -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_),
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 9ddea7bbd7..21a8536b05 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -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;
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index d8116933ca..8c223fd876 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -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_;
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
index c767564e02..a626f4d688 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
@@ -47,6 +47,7 @@ runTimeModifiable true;
functions
{
#include "streamLines"
+ #include "wallBoundedStreamLines"
#include "cuttingPlane"
#include "forceCoeffs"
#include "ensightWrite"
--
GitLab