zhyang-dev/OpenFOAM-6
1
0
Fork 0
mirror of https://github.com/OpenFOAM/OpenFOAM-6.git synced 2025-12-08 06:57:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

particle: Optimisation by inlining key methods

Browse Source 
Resolves bug report https://bugs.openfoam.org/view.php?id=2871
This commit is contained in:
Will Bainbridge
2018-03-13 08:22:33 +00:00
parent be604abcb9
commit 81947c8085
4 changed files with 135 additions and 136 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -42,34 +42,6 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::particle::stationaryTetGeometry
(
vector& centre,
vector& base,
vector& vertex1,
vector& vertex2
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
const vectorField& ccs = mesh_.cellCentres();
const pointField& pts = mesh_.points();
centre = ccs[celli_];
base = pts[triIs[0]];
vertex1 = pts[triIs[1]];
vertex2 = pts[triIs[2]];
}
Foam::barycentricTensor Foam::particle::stationaryTetTransform() const
{
vector centre, base, vertex1, vertex2;
stationaryTetGeometry(centre, base, vertex1, vertex2);
return barycentricTensor(centre, base, vertex1, vertex2);
}
void Foam::particle::stationaryTetReverseTransform
(
vector& centre,
@ -99,61 +71,6 @@ void Foam::particle::stationaryTetReverseTransform
}
void Foam::particle::movingTetGeometry
(
const scalar fraction,
Pair<vector>& centre,
Pair<vector>& base,
Pair<vector>& vertex1,
Pair<vector>& vertex2
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
const pointField& ptsOld = mesh_.oldPoints();
const pointField& ptsNew = mesh_.points();
// !!! <-- We would be better off using mesh_.cellCentres() here. However,
// we need to put a mesh_.oldCellCentres() method in for this to work. The
// values obtained from the mesh and those obtained from the cell do not
// necessarily match. See mantis #1993.
const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
// Old and new points and cell centres are not sub-cycled. If we are sub-
// cycling, then we have to account for the timestep change here by
// modifying the fractions that we take of the old and new geometry.
const Pair<scalar> s = stepFractionSpan();
const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
centre[0] = ccOld + f0*(ccNew - ccOld);
base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
centre[1] = f1*(ccNew - ccOld);
base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
}
Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
(
const scalar fraction
) const
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(fraction, centre, base, vertex1, vertex2);
return
Pair<barycentricTensor>
(
barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
);
}
void Foam::particle::movingTetReverseTransform
(
const scalar fraction,
@ -961,47 +878,6 @@ Foam::vector Foam::particle::deviationFromMeshCentre() const
}
void Foam::particle::patchData(vector& n, vector& U) const
{
if (!onBoundaryFace())
{
FatalErrorInFunction
<< "Patch data was requested for a particle that isn't on a patch"
<< exit(FatalError);
}
if (mesh_.moving())
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);
n = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
n /= mag(n);
// Interpolate the motion of the three face vertices to the current
// coordinates
U =
coordinates_.b()*base[1]
+ coordinates_.c()*vertex1[1]
+ coordinates_.d()*vertex2[1];
// The movingTetGeometry method gives the motion as a displacement
// across the time-step, so we divide by the time-step to get velocity
U /= mesh_.time().deltaTValue();
}
else
{
vector centre, base, vertex1, vertex2;
stationaryTetGeometry(centre, base, vertex1, vertex2);
n = triPointRef(base, vertex1, vertex2).normal();
n /= mag(n);
U = Zero;
}
}
void Foam::particle::transformProperties(const tensor&)
{}

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

@ -168,7 +168,7 @@ private:
// Tetrahedra functions
//- Get the vertices of the current tet
void stationaryTetGeometry
inline void stationaryTetGeometry
(
vector& centre,
vector& base,
@ -181,7 +181,7 @@ private:
// cartesian position in the global coordinate system. The
// conversion is x = A & y, where x is the cartesian position, y is
// the barycentric position and A is the transformation tensor.
barycentricTensor stationaryTetTransform() const;
inline barycentricTensor stationaryTetTransform() const;
//- Get the reverse transform associated with the current tet. The
// conversion is detA*y = (x - centre) & T. The variables x, y and
@ -200,7 +200,7 @@ private:
//- Get the vertices of the current moving tet. Two values are
// returned for each vertex. The first is a constant, and the
// second is a linear coefficient of the track fraction.
void movingTetGeometry
inline void movingTetGeometry
(
const scalar endStepFraction,
Pair<vector>& centre,
@ -213,7 +213,7 @@ private:
// This is of the same form as for the static case. As with the
// moving geometry, a linear function of the tracking fraction is
// returned for each component.
Pair<barycentricTensor> movingTetTransform
inline Pair<barycentricTensor> movingTetTransform
(
const scalar endStepFraction
) const;
@ -570,7 +570,7 @@ public:
// Patch data
//- Get the normal and velocity of the current patch location
void patchData(vector& n, vector& U) const;
inline void patchData(vector& n, vector& U) const;
// Transformations

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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,6 +26,91 @@ License
#include "polyMesh.H"
#include "Time.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::particle::stationaryTetGeometry
(
vector& centre,
vector& base,
vector& vertex1,
vector& vertex2
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
const vectorField& ccs = mesh_.cellCentres();
const pointField& pts = mesh_.points();
centre = ccs[celli_];
base = pts[triIs[0]];
vertex1 = pts[triIs[1]];
vertex2 = pts[triIs[2]];
}
inline Foam::barycentricTensor Foam::particle::stationaryTetTransform() const
{
vector centre, base, vertex1, vertex2;
stationaryTetGeometry(centre, base, vertex1, vertex2);
return barycentricTensor(centre, base, vertex1, vertex2);
}
inline void Foam::particle::movingTetGeometry
(
const scalar fraction,
Pair<vector>& centre,
Pair<vector>& base,
Pair<vector>& vertex1,
Pair<vector>& vertex2
) const
{
const triFace triIs(currentTetIndices().faceTriIs(mesh_));
const pointField& ptsOld = mesh_.oldPoints();
const pointField& ptsNew = mesh_.points();
// !!! <-- We would be better off using mesh_.cellCentres() here. However,
// we need to put a mesh_.oldCellCentres() method in for this to work. The
// values obtained from the mesh and those obtained from the cell do not
// necessarily match. See mantis #1993.
const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
// Old and new points and cell centres are not sub-cycled. If we are sub-
// cycling, then we have to account for the timestep change here by
// modifying the fractions that we take of the old and new geometry.
const Pair<scalar> s = stepFractionSpan();
const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
centre[0] = ccOld + f0*(ccNew - ccOld);
base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
centre[1] = f1*(ccNew - ccOld);
base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
}
inline Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
(
const scalar fraction
) const
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(fraction, centre, base, vertex1, vertex2);
return
Pair<barycentricTensor>
(
barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::particle::getNewParticleID() const
@ -207,4 +292,45 @@ inline Foam::vector Foam::particle::position() const
}
void Foam::particle::patchData(vector& n, vector& U) const
{
if (!onBoundaryFace())
{
FatalErrorInFunction
<< "Patch data was requested for a particle that isn't on a patch"
<< exit(FatalError);
}
if (mesh_.moving())
{
Pair<vector> centre, base, vertex1, vertex2;
movingTetGeometry(1, centre, base, vertex1, vertex2);
n = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
n /= mag(n);
// Interpolate the motion of the three face vertices to the current
// coordinates
U =
coordinates_.b()*base[1]
+ coordinates_.c()*vertex1[1]
+ coordinates_.d()*vertex2[1];
// The movingTetGeometry method gives the motion as a displacement
// across the time-step, so we divide by the time-step to get velocity
U /= mesh_.time().deltaTValue();
}
else
{
vector centre, base, vertex1, vertex2;
stationaryTetGeometry(centre, base, vertex1, vertex2);
n = triPointRef(base, vertex1, vertex2).normal();
n /= mag(n);
U = Zero;
}
}
// ************************************************************************* //

9
src/lagrangian/basic/particle/particleTemplates.C
View File

@ -122,7 +122,9 @@ void Foam::particle::hitFace
}
else if (onBoundaryFace())
{
if(!p.hitPatch(cloud, ttd))
changeToMasterPatch();
if (!p.hitPatch(cloud, ttd))
{
const polyPatch& patch = mesh_.boundaryMesh()[p.patch()];
@ -178,11 +180,6 @@ void Foam::particle::trackToAndHitFace
{
trackToFace(direction, fraction);
if (onBoundaryFace())
{
changeToMasterPatch();
}
hitFace(direction, cloud, td);
}