560 lines
16 KiB
C
560 lines
16 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
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/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "wallBoundedParticle.H"
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
template<class TrackData>
|
|
void Foam::wallBoundedParticle::patchInteraction
|
|
(
|
|
TrackData& td,
|
|
const scalar trackFraction
|
|
)
|
|
{
|
|
typedef typename TrackData::cloudType::particleType particleType;
|
|
|
|
particleType& p = static_cast<particleType&>(*this);
|
|
p.hitFace(td);
|
|
|
|
if (!internalFace(facei_))
|
|
{
|
|
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.
|
|
|
|
|
|
// 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 (facei_ != origFacei)
|
|
{
|
|
patchi = patch(facei_);
|
|
}
|
|
|
|
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);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
template<class TrackData>
|
|
Foam::scalar Foam::wallBoundedParticle::trackToEdge
|
|
(
|
|
TrackData& td,
|
|
const vector& endPosition
|
|
)
|
|
{
|
|
// Track particle to a given position and returns 1.0 if the
|
|
// trajectory is completed without hitting a face otherwise
|
|
// stops at the face and returns the fraction of the trajectory
|
|
// completed.
|
|
// on entry 'stepFraction()' should be set to the fraction of the
|
|
// time-step at which the tracking starts.
|
|
|
|
// Are we on a track face? If not we do a topological walk.
|
|
|
|
// Particle:
|
|
// - cell_ always set
|
|
// - tetFace_, tetPt_ always set (these identify tet particle is in)
|
|
// - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
|
|
|
|
//checkInside();
|
|
//checkOnTriangle(position());
|
|
//if (meshEdgeStart_ != -1 || diagEdge_ != -1)
|
|
//{
|
|
// checkOnEdge();
|
|
//}
|
|
|
|
scalar trackFraction = 0.0;
|
|
|
|
if (!td.isWallPatch_[tetFace()])
|
|
{
|
|
// Don't track across face. Just walk in cell. Particle is on
|
|
// mesh edge (as indicated by meshEdgeStart_).
|
|
const edge meshEdge(currentEdge());
|
|
|
|
// 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()))
|
|
{
|
|
label nbrCelli =
|
|
(
|
|
celli_ == mesh_.faceOwner()[facei_]
|
|
? mesh_.faceNeighbour()[facei_]
|
|
: mesh_.faceOwner()[facei_]
|
|
);
|
|
// 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_);
|
|
if ((nbrTi.faceTri(mesh_).normal() & (endPosition-position())) < 0)
|
|
{
|
|
// Change into nbrCell. No need to change tetFace, tetPt.
|
|
//Pout<< " crossed from cell:" << celli_
|
|
// << " into " << nbrCelli << endl;
|
|
celli_ = nbrCelli;
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else
|
|
{
|
|
// Walk to other face on edge. Changes tetFace, tetPt but not
|
|
// cell.
|
|
crossEdgeConnectedFace(meshEdge);
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Walk to other face on edge. This might give loop since
|
|
// particle should have been removed?
|
|
crossEdgeConnectedFace(meshEdge);
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// 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()))
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Can only track on boundary faces."
|
|
<< " Face:" << tetFace()
|
|
<< " at:" << mesh_.faceCentres()[tetFace()]
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
point projectedEndPosition = endPosition;
|
|
// Remove normal component
|
|
{
|
|
const triFace tri(currentTetIndices().faceTriIs(mesh_));
|
|
vector n = tri.normal(mesh_.points());
|
|
n /= mag(n);
|
|
const point& basePt = mesh_.points()[tri[0]];
|
|
projectedEndPosition -= ((projectedEndPosition-basePt)&n)*n;
|
|
}
|
|
|
|
|
|
bool doTrack = false;
|
|
if (meshEdgeStart_ == -1 && diagEdge_ == -1)
|
|
{
|
|
// We're starting and not yet on an edge.
|
|
doTrack = true;
|
|
}
|
|
else
|
|
{
|
|
// See if the current triangle has got a point on the
|
|
// correct side of the edge.
|
|
doTrack = isTriAlongTrack(projectedEndPosition);
|
|
}
|
|
|
|
|
|
if (doTrack)
|
|
{
|
|
// Track across triangle. Return triangle edge crossed.
|
|
label triEdgeI = -1;
|
|
trackFraction = trackFaceTri(projectedEndPosition, triEdgeI);
|
|
|
|
if (triEdgeI == -1)
|
|
{
|
|
// Reached endpoint
|
|
//checkInside();
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = -1;
|
|
return trackFraction;
|
|
}
|
|
|
|
const tetIndices ti(currentTetIndices());
|
|
|
|
// Triangle (faceTriIs) gets constructed from
|
|
// f[faceBasePtI_],
|
|
// f[facePtAI_],
|
|
// f[facePtBI_]
|
|
//
|
|
// So edge indices are:
|
|
// 0 : edge between faceBasePtI_ and facePtAI_
|
|
// 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();
|
|
|
|
if (triEdgeI == 0)
|
|
{
|
|
if (ti.facePtA() == f.fcIndex(fp0))
|
|
{
|
|
//Pout<< "Real edge." << endl;
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = fp0;
|
|
//checkOnEdge();
|
|
crossEdgeConnectedFace(currentEdge());
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else if (ti.facePtA() == f.rcIndex(fp0))
|
|
{
|
|
//Note: should not happen since boundary face so owner
|
|
//Pout<< "Real edge." << endl;
|
|
FatalErrorInFunction
|
|
<< abort(FatalError);
|
|
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = f.rcIndex(fp0);
|
|
//checkOnEdge();
|
|
crossEdgeConnectedFace(currentEdge());
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else
|
|
{
|
|
// Get index of triangle on other side of edge.
|
|
diagEdge_ = ti.facePtA()-fp0;
|
|
if (diagEdge_ < 0)
|
|
{
|
|
diagEdge_ += f.size();
|
|
}
|
|
meshEdgeStart_ = -1;
|
|
//checkOnEdge();
|
|
crossDiagonalEdge();
|
|
}
|
|
}
|
|
else if (triEdgeI == 1)
|
|
{
|
|
//Pout<< "Real edge." << endl;
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = ti.facePtA();
|
|
//checkOnEdge();
|
|
crossEdgeConnectedFace(currentEdge());
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else // if (triEdgeI == 2)
|
|
{
|
|
if (ti.facePtB() == f.rcIndex(fp0))
|
|
{
|
|
//Pout<< "Real edge." << endl;
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = ti.facePtB();
|
|
//checkOnEdge();
|
|
crossEdgeConnectedFace(currentEdge());
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else if (ti.facePtB() == f.fcIndex(fp0))
|
|
{
|
|
//Note: should not happen since boundary face so owner
|
|
//Pout<< "Real edge." << endl;
|
|
FatalErrorInFunction
|
|
<< abort(FatalError);
|
|
|
|
diagEdge_ = -1;
|
|
meshEdgeStart_ = fp0;
|
|
//checkOnEdge();
|
|
crossEdgeConnectedFace(currentEdge());
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else
|
|
{
|
|
//Pout<< "Triangle edge." << endl;
|
|
// Get index of triangle on other side of edge.
|
|
diagEdge_ = ti.facePtB()-fp0;
|
|
if (diagEdge_ < 0)
|
|
{
|
|
diagEdge_ += f.size();
|
|
}
|
|
meshEdgeStart_ = -1;
|
|
//checkOnEdge();
|
|
crossDiagonalEdge();
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Current tet is not the right one. Check the neighbour tet.
|
|
|
|
if (meshEdgeStart_ != -1)
|
|
{
|
|
// Particle is on mesh edge so change into other face on cell
|
|
crossEdgeConnectedFace(currentEdge());
|
|
//checkOnEdge();
|
|
patchInteraction(td, trackFraction);
|
|
}
|
|
else
|
|
{
|
|
// Particle is on diagonal edge so change into the other
|
|
// triangle.
|
|
crossDiagonalEdge();
|
|
//checkOnEdge();
|
|
}
|
|
}
|
|
}
|
|
|
|
//checkInside();
|
|
|
|
return trackFraction;
|
|
}
|
|
|
|
|
|
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>
|
|
void Foam::wallBoundedParticle::hitProcessorPatch
|
|
(
|
|
const processorPolyPatch& pp,
|
|
TrackData& td
|
|
)
|
|
{
|
|
// Switch particle
|
|
td.switchProcessor = true;
|
|
|
|
// Adapt edgeStart_ for other side.
|
|
// E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
|
|
// on the other side between 2 and 3 so edgeStart_ should be
|
|
// f.size()-edgeStart_-1.
|
|
|
|
const Foam::face& f = mesh_.faces()[face()];
|
|
|
|
if (meshEdgeStart_ != -1)
|
|
{
|
|
meshEdgeStart_ = f.size()-meshEdgeStart_-1;
|
|
}
|
|
else
|
|
{
|
|
// diagEdge_ is relative to faceBasePt
|
|
diagEdge_ = f.size()-diagEdge_;
|
|
}
|
|
}
|
|
|
|
|
|
template<class TrackData>
|
|
void Foam::wallBoundedParticle::hitWallPatch
|
|
(
|
|
const wallPolyPatch& wpp,
|
|
TrackData& td,
|
|
const tetIndices&
|
|
)
|
|
{}
|
|
|
|
|
|
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)
|
|
{
|
|
if (!c.size())
|
|
{
|
|
return;
|
|
}
|
|
|
|
particle::readFields(c);
|
|
|
|
IOField<label> meshEdgeStart
|
|
(
|
|
c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
|
|
);
|
|
|
|
IOField<label> diagEdge
|
|
(
|
|
c.fieldIOobject("diagEdge_", IOobject::MUST_READ)
|
|
);
|
|
c.checkFieldIOobject(c, diagEdge);
|
|
|
|
label i = 0;
|
|
forAllIter(typename CloudType, c, iter)
|
|
{
|
|
iter().meshEdgeStart_ = meshEdgeStart[i];
|
|
iter().diagEdge_ = diagEdge[i];
|
|
i++;
|
|
}
|
|
}
|
|
|
|
|
|
template<class CloudType>
|
|
void Foam::wallBoundedParticle::writeFields(const CloudType& c)
|
|
{
|
|
particle::writeFields(c);
|
|
|
|
label np = c.size();
|
|
|
|
IOField<label> meshEdgeStart
|
|
(
|
|
c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
|
|
np
|
|
);
|
|
IOField<label> diagEdge
|
|
(
|
|
c.fieldIOobject("diagEdge", IOobject::NO_READ),
|
|
np
|
|
);
|
|
|
|
label i = 0;
|
|
forAllConstIter(typename CloudType, c, iter)
|
|
{
|
|
meshEdgeStart[i] = iter().meshEdgeStart_;
|
|
diagEdge[i] = iter().diagEdge_;
|
|
i++;
|
|
}
|
|
|
|
meshEdgeStart.write();
|
|
diagEdge.write();
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|