Commit 6b316ba3 authored by Will Bainbridge's avatar Will Bainbridge Committed by Andrew Heather
Browse files

interpolation: Optimise by using particle local coordinates

This change changes the point-tetIndices-face interpolation function
method to take barycentric-tetIndices-face arguments instead. This
function is, at present, only used for interpolating Eulerian data to
Lagrangian particles.

This change prevents an inefficiency in cellPointInterpolation whereby
the position of the particle is calculated from it's barycentric
coordinates, before immediately being converted back to barycentric
coordinates to perform the interpolation.
parent 7a02a507
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -135,18 +135,24 @@ public:
const label facei = -1
) const = 0;
//- Interpolate field to the given point in the tetrahedron
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices. Calls interpolate function
// above here execpt where overridden by derived
// interpolation types.
virtual Type interpolate
(
const vector& position,
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const
{
return interpolate(position, tetIs.cell(), facei);
return
interpolate
(
tetIs.tet(pMesh_).barycentricToPoint(coordinates),
tetIs.cell(),
facei
);
}
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -75,6 +75,19 @@ public:
const label celli,
const label facei = -1
) const;
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices. This is an optimisation which skips
// calculating the position, as cell interpolation doesn't need it.
inline Type interpolate
(
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const
{
return interpolate(vector::zero, tetIs.cell(), facei);
}
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -77,6 +77,19 @@ public:
const label celli,
const label facei = -1
) const;
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices. This is an optimisation which skips
// calculating the position, as cell interpolation doesn't need it.
inline Type interpolate
(
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const
{
return interpolate(vector::zero, tetIs.cell(), facei);
}
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -86,11 +86,11 @@ public:
const label facei = -1
) const;
//- Interpolate field to the given point in the tetrahedron
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices.
inline Type interpolate
(
const vector& position,
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const;
......
......@@ -58,7 +58,7 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
template<class Type>
inline Type Foam::interpolationCellPoint<Type>::interpolate
(
const vector& position,
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei
) const
......@@ -79,21 +79,13 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
}
}
FixedList<scalar, 4> weights;
tetIs.tet(this->pMesh_).barycentric(position, weights);
// Order of weights is the same as that of the vertices of the tet, i.e.
// cellCentre, faceBasePt, facePtA, facePtB.
Type t = this->psi_[tetIs.cell()]*weights[0];
const triFace triIs = tetIs.faceTriIs(this->pMesh_);
t += psip_[triIs[0]]*weights[1];
t += psip_[triIs[1]]*weights[2];
t += psip_[triIs[2]]*weights[3];
return t;
return
this->psi_[tetIs.cell()]*coordinates[0]
+ psip_[triIs[0]]*coordinates[1]
+ psip_[triIs[1]]*coordinates[2]
+ psip_[triIs[2]]*coordinates[3];
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -78,11 +78,11 @@ public:
const label facei = -1
) const;
//- Interpolate field to the given point in the tetrahedron
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices.
inline Type interpolate
(
const vector& position,
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const;
......
......@@ -67,7 +67,7 @@ inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
template<class Type>
inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
(
const vector& position,
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei
) const
......@@ -101,7 +101,7 @@ inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
return interpolationCellPoint<Type>::interpolate
(
position,
coordinates,
tetIs,
facei
);
......
......@@ -48,7 +48,7 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
{
tetIndices tetIs = this->currentTetIndices();
rhoc_ = td.rhoInterp().interpolate(this->position(), tetIs);
rhoc_ = td.rhoInterp().interpolate(this->coordinates(), tetIs);
if (rhoc_ < td.cloud().constProps().rhoMin())
{
......@@ -62,9 +62,9 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
rhoc_ = td.cloud().constProps().rhoMin();
}
Uc_ = td.UInterp().interpolate(this->position(), tetIs);
Uc_ = td.UInterp().interpolate(this->coordinates(), tetIs);
muc_ = td.muInterp().interpolate(this->position(), tetIs);
muc_ = td.muInterp().interpolate(this->coordinates(), tetIs);
// Apply dispersion components to carrier phase velocity
Uc_ = td.cloud().dispersion().update
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -211,7 +211,7 @@ void Foam::ReactingParcel<ParcelType>::setCellValues
pc_ = td.pInterp().interpolate
(
this->position(),
this->coordinates(),
this->currentTetIndices()
);
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -43,9 +43,9 @@ void Foam::ThermoParcel<ParcelType>::setCellValues
tetIndices tetIs = this->currentTetIndices();
Cpc_ = td.CpInterp().interpolate(this->position(), tetIs);
Cpc_ = td.CpInterp().interpolate(this->coordinates(), tetIs);
Tc_ = td.TInterp().interpolate(this->position(), tetIs);
Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
if (Tc_ < td.cloud().constProps().TMin())
{
......@@ -124,8 +124,8 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
rhos = this->rhoc_*TRatio;
tetIndices tetIs = this->currentTetIndices();
mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
kappas = td.kappaInterp().interpolate(this->position(), tetIs)/TRatio;
mus = td.muInterp().interpolate(this->coordinates(), tetIs)/TRatio;
kappas = td.kappaInterp().interpolate(this->coordinates(), tetIs)/TRatio;
Pr = Cpc_*mus/kappas;
Pr = max(ROOTVSMALL, Pr);
......@@ -286,7 +286,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
if (td.cloud().radiation())
{
tetIndices tetIs = this->currentTetIndices();
const scalar Gc = td.GInterp().interpolate(this->position(), tetIs);
const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs);
const scalar sigma = physicoChemical::sigma.value();
const scalar epsilon = td.cloud().constProps().epsilon0();
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -137,7 +137,7 @@ Foam::forceSuSp Foam::LiftForce<CloudType>::calcCoupled
forceSuSp value(Zero, 0.0);
vector curlUc =
curlUcInterp().interpolate(p.position(), p.currentTetIndices());
curlUcInterp().interpolate(p.coordinates(), p.currentTetIndices());
scalar Cl = this->Cl(p, curlUc, Re, muc);
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -110,7 +110,7 @@ Foam::forceSuSp Foam::ParamagneticForce<CloudType>::calcNonCoupled
value.Su()=
mass*3.0*constant::electromagnetic::mu0.value()/p.rho()
*magneticSusceptibility_/(magneticSusceptibility_ + 3)
*HdotGradHInterp.interpolate(p.position(), p.currentTetIndices());
*HdotGradHInterp.interpolate(p.coordinates(), p.currentTetIndices());
return value;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -128,7 +128,7 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
forceSuSp value(Zero, 0.0);
vector DUcDt =
DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
DUcDtInterp().interpolate(p.coordinates(), p.currentTetIndices());
value.Su() = mass*p.rhoc()/p.rho()*DUcDt;
......
......@@ -55,7 +55,6 @@ bool Foam::solidParticle::move
}
const label celli = cell();
const scalar sfrac = stepFraction();
const scalar f = 1 - stepFraction();
......@@ -63,10 +62,10 @@ bool Foam::solidParticle::move
const scalar dt = (stepFraction() - sfrac)*trackTime;
cellPointWeight cpw(mesh(), position(), celli, face());
scalar rhoc = td.rhoInterp().interpolate(cpw);
vector Uc = td.UInterp().interpolate(cpw);
scalar nuc = td.nuInterp().interpolate(cpw);
const tetIndices tetIs = this->currentTetIndices();
scalar rhoc = td.rhoInterp().interpolate(this->coordinates(), tetIs);
vector Uc = td.UInterp().interpolate(this->coordinates(), tetIs);
scalar nuc = td.nuInterp().interpolate(this->coordinates(), tetIs);
scalar rhop = td.cloud().rhop();
scalar magUr = mag(Uc - U_);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment