songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
d3a079b4da
openfoam/src/lagrangian/basic/particle/particle.C
mattijs 177176c5ff BUG: particle: use correct celli. Fixes #1992 
Was checking the old celli instead of the result of
re-finding the position. See also Foundation commit 50a965f8866683a81d79cbc7811af7333baf9d10.
2023-10-10 15:40:47 +01:00

1265 lines
32 KiB
C
Raw Blame History

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017, 2020 OpenFOAM Foundation
Copyright (C) 2018-2021 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/>.
\*---------------------------------------------------------------------------*/
#include "particle.H"
#include "transform.H"
#include "treeDataCell.H"
#include "cubicEqn.H"
#include "registerSwitch.H"
#include "indexedOctree.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(particle, 0);
}
const Foam::label Foam::particle::maxNBehind_ = 10;
Foam::label Foam::particle::particleCount_ = 0;
bool Foam::particle::writeLagrangianCoordinates = true;
bool Foam::particle::writeLagrangianPositions
(
Foam::debug::infoSwitch("writeLagrangianPositions", 1)
);
registerInfoSwitch
(
"writeLagrangianPositions",
bool,
Foam::particle::writeLagrangianPositions
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::particle::stationaryTetReverseTransform
(
vector& centre,
scalar& detA,
barycentricTensor& T
) const
{
barycentricTensor A = stationaryTetTransform();
const vector ab = A.b() - A.a();
const vector ac = A.c() - A.a();
const vector ad = A.d() - A.a();
const vector bc = A.c() - A.b();
const vector bd = A.d() - A.b();
centre = A.a();
detA = ab & (ac ^ ad);
T = barycentricTensor
(
bd ^ bc,
ac ^ ad,
ad ^ ab,
ab ^ ac
);
}
void Foam::particle::movingTetReverseTransform
(
const scalar fraction,
Pair<vector>& centre,
FixedList<scalar, 4>& detA,
FixedList<barycentricTensor, 3>& T
) const
{
Pair<barycentricTensor> A = movingTetTransform(fraction);
const Pair<vector> ab(A[0].b() - A[0].a(), A[1].b() - A[1].a());
const Pair<vector> ac(A[0].c() - A[0].a(), A[1].c() - A[1].a());
const Pair<vector> ad(A[0].d() - A[0].a(), A[1].d() - A[1].a());
const Pair<vector> bc(A[0].c() - A[0].b(), A[1].c() - A[1].b());
const Pair<vector> bd(A[0].d() - A[0].b(), A[1].d() - A[1].b());
centre[0] = A[0].a();
centre[1] = A[1].a();
detA[0] = ab[0] & (ac[0] ^ ad[0]);
detA[1] =
(ab[1] & (ac[0] ^ ad[0]))
+ (ab[0] & (ac[1] ^ ad[0]))
+ (ab[0] & (ac[0] ^ ad[1]));
detA[2] =
(ab[0] & (ac[1] ^ ad[1]))
+ (ab[1] & (ac[0] ^ ad[1]))
+ (ab[1] & (ac[1] ^ ad[0]));
detA[3] = ab[1] & (ac[1] ^ ad[1]);
T[0] = barycentricTensor
(
bd[0] ^ bc[0],
ac[0] ^ ad[0],
ad[0] ^ ab[0],
ab[0] ^ ac[0]
);
T[1] = barycentricTensor
(
(bd[0] ^ bc[1]) + (bd[1] ^ bc[0]),
(ac[0] ^ ad[1]) + (ac[1] ^ ad[0]),
(ad[0] ^ ab[1]) + (ad[1] ^ ab[0]),
(ab[0] ^ ac[1]) + (ab[1] ^ ac[0])
);
T[2] = barycentricTensor
(
bd[1] ^ bc[1],
ac[1] ^ ad[1],
ad[1] ^ ab[1],
ab[1] ^ ac[1]
);
}
void Foam::particle::reflect()
{
Swap(coordinates_.c(), coordinates_.d());
}
void Foam::particle::rotate(const bool reverse)
{
if (!reverse)
{
scalar temp = coordinates_.b();
coordinates_.b() = coordinates_.c();
coordinates_.c() = coordinates_.d();
coordinates_.d() = temp;
}
else
{
scalar temp = coordinates_.d();
coordinates_.d() = coordinates_.c();
coordinates_.c() = coordinates_.b();
coordinates_.b() = temp;
}
}
void Foam::particle::changeTet(const label tetTriI)
{
const bool isOwner = mesh_.faceOwner()[tetFacei_] == celli_;
const label firstTetPtI = 1;
const label lastTetPtI = mesh_.faces()[tetFacei_].size() - 2;
if (tetTriI == 1)
{
changeFace(tetTriI);
}
else if (tetTriI == 2)
{
if (isOwner)
{
if (tetPti_ == lastTetPtI)
{
changeFace(tetTriI);
}
else
{
reflect();
tetPti_ += 1;
}
}
else
{
if (tetPti_ == firstTetPtI)
{
changeFace(tetTriI);
}
else
{
reflect();
tetPti_ -= 1;
}
}
}
else if (tetTriI == 3)
{
if (isOwner)
{
if (tetPti_ == firstTetPtI)
{
changeFace(tetTriI);
}
else
{
reflect();
tetPti_ -= 1;
}
}
else
{
if (tetPti_ == lastTetPtI)
{
changeFace(tetTriI);
}
else
{
reflect();
tetPti_ += 1;
}
}
}
else
{
FatalErrorInFunction
<< "Changing tet without changing cell should only happen when the "
<< "track is on triangle 1, 2 or 3."
<< exit(FatalError);
}
}
void Foam::particle::changeFace(const label tetTriI)
{
// Get the old topology
const triFace triOldIs(currentTetIndices().faceTriIs(mesh_));
// Get the shared edge and the pre-rotation
edge sharedEdge;
if (tetTriI == 1)
{
sharedEdge = edge(triOldIs[1], triOldIs[2]);
}
else if (tetTriI == 2)
{
sharedEdge = edge(triOldIs[2], triOldIs[0]);
}
else if (tetTriI == 3)
{
sharedEdge = edge(triOldIs[0], triOldIs[1]);
}
else
{
FatalErrorInFunction
<< "Changing face without changing cell should only happen when the"
<< " track is on triangle 1, 2 or 3."
<< exit(FatalError);
sharedEdge = edge(-1, -1);
}
// Find the face in the same cell that shares the edge, and the
// corresponding tetrahedra point
tetPti_ = -1;
forAll(mesh_.cells()[celli_], cellFaceI)
{
const label newFaceI = mesh_.cells()[celli_][cellFaceI];
const class face& newFace = mesh_.faces()[newFaceI];
const label newOwner = mesh_.faceOwner()[newFaceI];
// Exclude the current face
if (tetFacei_ == newFaceI)
{
continue;
}
// Loop over the edges, looking for the shared one. Note that we have to
// match the direction of the edge as well as the end points in order to
// avoid false positives when dealing with coincident ACMI faces.
const label edgeComp = newOwner == celli_ ? -1 : +1;
label edgeI = 0;
for
(
;
edgeI < newFace.size()
&& edge::compare(sharedEdge, newFace.faceEdge(edgeI)) != edgeComp;
++ edgeI
);
// If the face does not contain the edge, then move on to the next face
if (edgeI >= newFace.size())
{
continue;
}
// Make the edge index relative to the base point
const label newBaseI = max(0, mesh_.tetBasePtIs()[newFaceI]);
edgeI = (edgeI - newBaseI + newFace.size()) % newFace.size();
// If the edge is next the base point (i.e., the index is 0 or n - 1),
// then we swap it for the adjacent edge. This new edge is opposite the
// base point, and defines the tet with the original edge in it.
edgeI = min(max(1, edgeI), newFace.size() - 2);
// Set the new face and tet point
tetFacei_ = newFaceI;
tetPti_ = edgeI;
// Exit the loop now that the tet point has been found
break;
}
if (tetPti_ == -1)
{
FatalErrorInFunction
<< "The search for an edge-connected face and tet-point failed."
<< exit(FatalError);
}
// Pre-rotation puts the shared edge opposite the base of the tetrahedron
if (sharedEdge.otherVertex(triOldIs[1]) == -1)
{
rotate(false);
}
else if (sharedEdge.otherVertex(triOldIs[2]) == -1)
{
rotate(true);
}
// Get the new topology
const triFace triNewIs = currentTetIndices().faceTriIs(mesh_);
// Reflect to account for the change of triangle orientation on the new face
reflect();
// Post rotation puts the shared edge back in the correct location
if (sharedEdge.otherVertex(triNewIs[1]) == -1)
{
rotate(true);
}
else if (sharedEdge.otherVertex(triNewIs[2]) == -1)
{
rotate(false);
}
}
void Foam::particle::changeCell()
{
// Set the cell to be the one on the other side of the face
const label ownerCellI = mesh_.faceOwner()[tetFacei_];
const bool isOwner = celli_ == ownerCellI;
celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
// Reflect to account for the change of triangle orientation in the new cell
reflect();
}
void Foam::particle::changeToMasterPatch()
{
label thisPatch = patch();
forAll(mesh_.cells()[celli_], cellFaceI)
{
// Skip the current face and any internal faces
const label otherFaceI = mesh_.cells()[celli_][cellFaceI];
if (facei_ == otherFaceI || mesh_.isInternalFace(otherFaceI))
{
continue;
}
// Compare the two faces. If they are the same, chose the one with the
// lower patch index. In the case of an ACMI-wall pair, this will be
// the ACMI, which is what we want.
const class face& thisFace = mesh_.faces()[facei_];
const class face& otherFace = mesh_.faces()[otherFaceI];
if (face::compare(thisFace, otherFace) != 0)
{
const label otherPatch =
mesh_.boundaryMesh().whichPatch(otherFaceI);
if (thisPatch > otherPatch)
{
facei_ = otherFaceI;
thisPatch = otherPatch;
}
}
}
tetFacei_ = facei_;
}
void Foam::particle::locate
(
const vector& position,
const vector* direction,
const label celli,
const bool boundaryFail,
const string& boundaryMsg
)
{
if (debug)
{
Pout<< "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
}
celli_ = celli;
// Find the cell, if it has not been given
if (celli_ < 0)
{
celli_ = mesh_.cellTree().findInside(position);
if (celli_ < 0)
{
FatalErrorInFunction
<< "Cell not found for particle position " << position << "."
<< exit(FatalError);
}
}
const vector displacement = position - mesh_.cellCentres()[celli_];
// Loop all cell tets to find the one containing the position. Track
// through each tet from the cell centre. If a tet contains the position
// then the track will end with a single trackToTri.
const class cell& c = mesh_.cells()[celli_];
scalar minF = VGREAT;
label minTetFacei = -1, minTetPti = -1;
forAll(c, cellTetFacei)
{
const class face& f = mesh_.faces()[c[cellTetFacei]];
for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
{
coordinates_ = barycentric(1, 0, 0, 0);
tetFacei_ = c[cellTetFacei];
tetPti_ = tetPti;
facei_ = -1;
label tetTriI = -1;
const scalar f = trackToTri(displacement, 0, tetTriI);
if (tetTriI == -1)
{
return;
}
if (f < minF)
{
minF = f;
minTetFacei = tetFacei_;
minTetPti = tetPti_;
}
}
}
// The particle must be (hopefully only slightly) outside the cell. Track
// into the tet which got the furthest.
coordinates_ = barycentric(1, 0, 0, 0);
tetFacei_ = minTetFacei;
tetPti_ = minTetPti;
facei_ = -1;
track(displacement, 0);
if (!onFace())
{
return;
}
// If we are here then we hit a boundary
if (boundaryFail)
{
FatalErrorInFunction << boundaryMsg << exit(FatalError);
}
else
{
static label nWarnings = 0;
static const label maxNWarnings = 100;
if ((nWarnings < maxNWarnings) && boundaryFail)
{
WarningInFunction << boundaryMsg.c_str() << endl;
++ nWarnings;
}
if (nWarnings == maxNWarnings)
{
WarningInFunction
<< "Suppressing any further warnings about particles being "
<< "located outside of the mesh." << endl;
++ nWarnings;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::particle::particle
(
const polyMesh& mesh,
const barycentric& coordinates,
const label celli,
const label tetFacei,
const label tetPti
)
:
mesh_(mesh),
coordinates_(coordinates),
celli_(celli),
tetFacei_(tetFacei),
tetPti_(tetPti),
facei_(-1),
stepFraction_(1.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{}
Foam::particle::particle
(
const polyMesh& mesh,
const vector& position,
const label celli
)
:
mesh_(mesh),
coordinates_(-VGREAT, -VGREAT, -VGREAT, -VGREAT),
celli_(celli),
tetFacei_(-1),
tetPti_(-1),
facei_(-1),
stepFraction_(1.0),
behind_(0.0),
nBehind_(0),
origProc_(Pstream::myProcNo()),
origId_(getNewParticleID())
{
locate
(
position,
nullptr,
celli,
false,
"Particle initialised with a location outside of the mesh"
);
}
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_(1.0),
behind_(0.0),
nBehind_(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_),
coordinates_(p.coordinates_),
celli_(p.celli_),
tetFacei_(p.tetFacei_),
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
behind_(p.behind_),
nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
Foam::particle::particle(const particle& p, const polyMesh& mesh)
:
mesh_(mesh),
coordinates_(p.coordinates_),
celli_(p.celli_),
tetFacei_(p.tetFacei_),
tetPti_(p.tetPti_),
facei_(p.facei_),
stepFraction_(p.stepFraction_),
behind_(p.behind_),
nBehind_(p.nBehind_),
origProc_(p.origProc_),
origId_(p.origId_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::particle::track
(
const vector& displacement,
const scalar fraction
)
{
scalar f = trackToFace(displacement, fraction);
while (onInternalFace())
{
changeCell();
f *= trackToFace(f*displacement, f*fraction);
}
return f;
}
Foam::scalar Foam::particle::trackToFace
(
const vector& displacement,
const scalar fraction
)
{
scalar f = 1;
label tetTriI = onFace() ? 0 : -1;
facei_ = -1;
// Loop the tets in the current cell
while (nBehind_ < maxNBehind_)
{
f *= trackToTri(f*displacement, f*fraction, tetTriI);
if (tetTriI == -1)
{
// The track has completed within the current tet
return 0;
}
else if (tetTriI == 0)
{
// The track has hit a face, so set the current face and return
facei_ = tetFacei_;
return f;
}
else
{
// Move to the next tet and continue the track
changeTet(tetTriI);
}
}
// Warn if stuck, and incorrectly advance the step fraction to completion
static label stuckID = -1, stuckProc = -1;
if (origId_ != stuckID && origProc_ != stuckProc)
{
WarningInFunction
<< "Particle #" << origId_ << " got stuck at " << position()
<< endl;
}
stuckID = origId_;
stuckProc = origProc_;
stepFraction_ += f*fraction;
behind_ = 0;
nBehind_ = 0;
return 0;
}
Foam::scalar Foam::particle::trackToStationaryTri
(
const vector& displacement,
const scalar fraction,
label& tetTriI
)
{
const vector x0 = position();
const vector x1 = displacement;
const barycentric y0 = coordinates_;
if (debug)
{
Pout<< "Particle " << origId() << endl << "Tracking from " << x0
<< " along " << x1 << " to " << x0 + x1 << endl;
}
// Get the tet geometry
vector centre;
scalar detA;
barycentricTensor T;
stationaryTetReverseTransform(centre, detA, T);
if (debug)
{
vector o, b, v1, v2;
stationaryTetGeometry(o, b, v1, v2);
Pout<< "Tet points o=" << o << ", b=" << b
<< ", v1=" << v1 << ", v2=" << v2 << endl
<< "Tet determinant = " << detA << endl
<< "Start local coordinates = " << y0 << endl;
}
// Calculate the local tracking displacement
barycentric Tx1(x1 & T);
if (debug)
{
Pout<< "Local displacement = " << Tx1 << "/" << detA << endl;
}
// Calculate the hit fraction
label iH = -1;
scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
for (label i = 0; i < 4; ++ i)
{
if (Tx1[i] < - detA*SMALL)
{
scalar mu = - y0[i]/Tx1[i];
if (debug)
{
Pout<< "Hit on tet face " << i << " at local coordinate "
<< y0 + mu*Tx1 << ", " << mu*detA*100 << "% of the "
<< "way along the track" << endl;
}
if (0 <= mu && mu < muH)
{
iH = i;
muH = mu;
}
}
}
// Set the new coordinates
barycentric yH = y0 + muH*Tx1;
// Clamp to zero any negative coordinates generated by round-off error
for (label i = 0; i < 4; ++ i)
{
yH.replace(i, i == iH ? 0 : max(0, yH[i]));
}
// Re-normalise if within the tet
if (iH == -1)
{
yH /= cmptSum(yH);
}
// Set the new position and hit index
coordinates_ = yH;
tetTriI = iH;
if (debug)
{
if (iH != -1)
{
Pout<< "Track hit tet face " << iH << " first" << endl;
}
else
{
Pout<< "Track hit no tet faces" << endl;
}
Pout<< "End local coordinates = " << yH << endl
<< "End global coordinates = " << position() << endl
<< "Tracking displacement = " << position() - x0 << endl
<< muH*detA*100 << "% of the step from " << stepFraction_ << " to "
<< stepFraction_ + fraction << " completed" << endl << endl;
}
// Set the proportion of the track that has been completed
stepFraction_ += fraction*muH*detA;
if (debug)
{
Pout << "Step Fraction : " << stepFraction_ << endl;
Pout << "fraction*muH*detA : " << fraction*muH*detA << endl;
Pout << "static muH : " << muH << endl;
Pout << "origId() : " << origId() << endl;
}
// Accumulate displacement behind
if (detA <= 0 || nBehind_ > 0)
{
behind_ += muH*detA*mag(displacement);
if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}
return iH != -1 ? 1 - muH*detA : 0;
}
Foam::scalar Foam::particle::trackToMovingTri
(
const vector& displacement,
const scalar fraction,
label& tetTriI
)
{
const vector x0 = position();
const vector x1 = displacement;
const barycentric y0 = coordinates_;
if (debug)
{
Pout<< "Particle " << origId() << endl << "Tracking from " << x0
<< " along " << x1 << " to " << x0 + x1 << endl;
}
// Get the tet geometry
Pair<vector> centre;
FixedList<scalar, 4> detA;
FixedList<barycentricTensor, 3> T;
movingTetReverseTransform(fraction, centre, detA, T);
if (debug)
{
Pair<vector> o, b, v1, v2;
movingTetGeometry(fraction, o, b, v1, v2);
Pout<< "Tet points o=" << o[0] << ", b=" << b[0]
<< ", v1=" << v1[0] << ", v2=" << v2[0] << endl
<< "Tet determinant = " << detA[0] << endl
<< "Start local coordinates = " << y0[0] << endl;
}
// Get the relative global position
const vector x0Rel = x0 - centre[0];
const vector x1Rel = x1 - centre[1];
// Form the determinant and hit equations
cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
const barycentric yC(1, 0, 0, 0);
const barycentric hitEqnA =
((x1Rel & T[2]) + detA[3]*yC)*sqr(detA[0]);
const barycentric hitEqnB =
((x1Rel & T[1]) + (x0Rel & T[2]) + detA[2]*yC)*detA[0];
const barycentric hitEqnC =
((x1Rel & T[0]) + (x0Rel & T[1]) + detA[1]*yC);
const barycentric hitEqnD = y0;
FixedList<cubicEqn, 4> hitEqn;
forAll(hitEqn, i)
{
hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
}
if (debug)
{
for (label i = 0; i < 4; ++ i)
{
Pout<< (i ? " " : "Hit equation ") << i << " = "
<< hitEqn[i] << endl;
}
Pout<< " DetA equation = " << detA << endl;
}
// Calculate the hit fraction
label iH = -1;
scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
for (label i = 0; i < 4; ++i)
{
const Roots<3> mu = hitEqn[i].roots();
for (label j = 0; j < 3; ++j)
{
if
(
mu.type(j) == roots::real
&& hitEqn[i].derivative(mu[j]) < - detA[0]*SMALL
)
{
if (debug)
{
const barycentric yH
(
hitEqn[0].value(mu[j]),
hitEqn[1].value(mu[j]),
hitEqn[2].value(mu[j]),
hitEqn[3].value(mu[j])
);
const scalar detAH = detAEqn.value(mu[j]);
Pout<< "Hit on tet face " << i << " at local coordinate "
<< (std::isnormal(detAH) ? name(yH/detAH) : "???")
<< ", " << mu[j]*detA[0]*100 << "% of the "
<< "way along the track" << endl;
Pout<< "derivative : " << hitEqn[i].derivative(mu[j]) << nl
<< " coord " << j << " mu[j]: " << mu[j] << nl
<< " hitEq " << i << endl;
}
if (0 <= mu[j] && mu[j] < muH)
{
iH = i;
muH = mu[j];
}
}
}
}
// Set the new coordinates
barycentric yH
(
hitEqn[0].value(muH),
hitEqn[1].value(muH),
hitEqn[2].value(muH),
hitEqn[3].value(muH)
);
// !!! <-- This fails if the tet collapses onto the particle, as detA tends
// to zero at the hit. In this instance, we can differentiate the hit and
// detA polynomials to find a limiting location, but this will not be on a
// triangle. We will then need to do a second track through the degenerate
// tet to find the final hit position. This second track is over zero
// distance and therefore can be of the static mesh type. This has not yet
// been implemented.
const scalar detAH = detAEqn.value(muH);
if (!std::isnormal(detAH))
{
FatalErrorInFunction
<< "A moving tet collapsed onto a particle. This is not supported. "
<< "The mesh is too poor, or the motion too severe, for particle "
<< "tracking to function." << exit(FatalError);
}
yH /= detAH;
// Clamp to zero any negative coordinates generated by round-off error
for (label i = 0; i < 4; ++ i)
{
yH.replace(i, i == iH ? 0 : max(0, yH[i]));
}
// Re-normalise if within the tet
if (iH == -1)
{
yH /= cmptSum(yH);
}
// Set the new position and hit index
coordinates_ = yH;
tetTriI = iH;
scalar advance = muH*detA[0];
// Set the proportion of the track that has been completed
stepFraction_ += fraction*advance;
// Accumulate displacement behind
if (detA[0] <= 0 || nBehind_ > 0)
{
behind_ += muH*detA[0]*mag(displacement);
if (behind_ > 0)
{
behind_ = 0;
nBehind_ = 0;
}
else
{
++ nBehind_;
}
}
if (debug)
{
if (iH != -1)
{
Pout<< "Track hit tet face " << iH << " first" << endl;
}
else
{
Pout<< "Track hit no tet faces" << endl;
}
// Pout<< "End local coordinates = " << yH << endl
// << "End global coordinates = " << position() << endl
// << "Tracking displacement = " << position() - x0 << endl
// << muH*detA[0]*100 << "% of the step from " << stepFraction_
// << " to " << stepFraction_ + fraction << " completed" << endl
// << endl;
}
return iH != -1 ? 1 - muH*detA[0] : 0;
}
Foam::scalar Foam::particle::trackToTri
(
const vector& displacement,
const scalar fraction,
label& tetTriI
)
{
if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
{
return trackToMovingTri(displacement, fraction, tetTriI);
}
else
{
return trackToStationaryTri(displacement, fraction, tetTriI);
}
}
Foam::vector Foam::particle::deviationFromMeshCentre() const
{
if (cmptMin(mesh_.geometricD()) == -1)
{
vector pos = position(), posC = pos;
meshTools::constrainToMeshCentre(mesh_, posC);
return pos - posC;
}
else
{
return vector::zero;
}
}
void Foam::particle::transformProperties(const tensor&)
{}
void Foam::particle::transformProperties(const vector&)
{}
void Foam::particle::prepareForParallelTransfer()
{
// Convert the face index to be local to the processor patch
facei_ = mesh_.boundaryMesh()[patch()].whichFace(facei_);
}
void Foam::particle::correctAfterParallelTransfer
(
const label patchi,
trackingData& td
)
{
const coupledPolyPatch& ppp =
refCast<const coupledPolyPatch>(mesh_.boundaryMesh()[patchi]);
if (!ppp.parallel())
{
const tensor& T =
(
ppp.forwardT().size() == 1
? ppp.forwardT()[0]
: ppp.forwardT()[facei_]
);
transformProperties(T);
}
else if (ppp.separated())
{
const vector& s =
(
(ppp.separation().size() == 1)
? ppp.separation()[0]
: ppp.separation()[facei_]
);
transformProperties(-s);
}
// Set the topology
celli_ = ppp.faceCells()[facei_];
facei_ += ppp.start();
tetFacei_ = facei_;
// Faces either side of a coupled patch are numbered in opposite directions
// as their normals both point away from their connected cells. The tet
// point therefore counts in the opposite direction from the base point.
tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
// Reflect to account for the change of triangle orientation in the new cell
reflect();
// Note that the position does not need transforming explicitly. The face-
// triangle on the receive patch is the transformation of the one on the
// send patch, so whilst the barycentric coordinates remain the same, the
// change of triangle implicitly transforms the position.
}
void Foam::particle::prepareForInteractionListReferral
(
const vectorTensorTransform& transform
)
{
// Get the transformed position
const vector pos = transform.invTransformPosition(position());
// Break the topology
celli_ = -1;
tetFacei_ = -1;
tetPti_ = -1;
facei_ = -1;
// Store the position in the barycentric data
coordinates_ = barycentric(1 - cmptSum(pos), pos.x(), pos.y(), pos.z());
// Transform the properties
transformProperties(- transform.t());
if (transform.hasR())
{
transformProperties(transform.R().T());
}
}
void Foam::particle::correctAfterInteractionListReferral(const label celli)
{
// Get the position from the barycentric data
const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
// Create some arbitrary topology for the supplied cell
celli_ = celli;
tetFacei_ = mesh_.cells()[celli_][0];
tetPti_ = 1;
facei_ = -1;
// Get the reverse transform and directly set the coordinates from the
// position. This isn't likely to be correct; the particle is probably not
// in this tet. It will, however, generate the correct vector when the
// position method is called. A referred particle should never be tracked,
// so this approximate topology is good enough. By using the nearby cell we
// minimise the error associated with the incorrect topology.
coordinates_ = barycentric(1, 0, 0, 0);
if (mesh_.moving() && stepFraction_ != 1)
{
Pair<vector> centre;
FixedList<scalar, 4> detA;
FixedList<barycentricTensor, 3> T;
movingTetReverseTransform(0, centre, detA, T);
coordinates_ += (pos - centre[0]) & T[0]/detA[0];
}
else
{
vector centre;
scalar detA;
barycentricTensor T;
stationaryTetReverseTransform(centre, detA, T);
coordinates_ += (pos - centre) & T/detA;
}
}
Foam::label Foam::particle::procTetPt
(
const polyMesh& procMesh,
const label procCell,
const label procTetFace
) const
{
// The tet point on the procMesh differs from the current tet point if the
// mesh and procMesh faces are of differing orientation. The change is the
// same as in particle::correctAfterParallelTransfer.
if
(
(mesh_.faceOwner()[tetFacei_] == celli_)
== (procMesh.faceOwner()[procTetFace] == procCell)
)
{
return tetPti_;
}
else
{
return procMesh.faces()[procTetFace].size() - 1 - tetPti_;
}
}
void Foam::particle::autoMap
(
const vector& position,
const mapPolyMesh& mapper
)
{
locate
(
position,
nullptr,
mapper.reverseCellMap()[celli_],
true,
"Particle mapped to a location outside of the mesh"
);
}
void Foam::particle::relocate(const point& position, const label celli)
{
locate
(
position,
nullptr,
celli,
true,
"Particle mapped to a location outside of the mesh"
);
}
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
bool Foam::operator==(const particle& pA, const particle& pB)
{
return (pA.origProc() == pB.origProc() && pA.origId() == pB.origId());
}
bool Foam::operator!=(const particle& pA, const particle& pB)
{
return !(pA == pB);
}
// ************************************************************************* //
Reference in New Issue View Git Blame Copy Permalink