OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

add CompatibleOperators::(Residual|Weak)::Convection

Browse Source
This commit is contained in:
Arne Morten Kvarving 2016-12-01 12:26:25 +01:00 committed by timovanopstal
parent 85a7024501
commit faa0e9adf0
2 changed files with 67 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
Apps/Common/CompatibleOperators.C
View File

@ -28,6 +28,32 @@ void CompatibleOperators::Weak::Advection(std::vector<Matrix>& EM,
EM[n](i,j) += scale*AC[k]*fe.grad(n)(j,k)*fe.basis(n)(i)*fe.detJxW;
}
void CompatibleOperators::Weak::Convection(std::vector<Matrix>& EM,
const FiniteElement& fe,
const Vec3& U,
const Tensor& dUdX,
double scale,
bool conservative)
{
// Convection
static const double vidx[3][3] = {{1, 6, 7},
{10, 2, 11},
{14, 15, 3}};
size_t nsd = fe.grad(1).cols();
for (size_t m=1; m <= nsd; ++m)
for (size_t i=1; i <= fe.basis(m).size(); ++i)
for (size_t n = 1; n <= nsd; ++n)
for (size_t j=1; j <= fe.basis(n).size(); ++j) {
double conv = fe.basis(n)(j) * dUdX(m,n);
if (m == n)
for (size_t k = 1; k <= nsd; ++k)
conv += U[k-1] * fe.grad(n)(j,k);
EM[vidx[m-1][n-1]](i,j) += scale * conv * fe.basis(m)(i) * fe.detJxW;
}
}
void CompatibleOperators::Weak::Gradient(std::vector<Matrix>& EM,
const FiniteElement& fe,
double scale)
@ -92,6 +118,22 @@ void CompatibleOperators::Weak::Source(Vectors& EV,
}
void CompatibleOperators::Residual::Convection(Vectors& EV, const FiniteElement& fe,
const Vec3& U, const Tensor& dUdX,
const Vec3& UC, double scale,
bool conservative)
{
size_t nsd = fe.grad(1).cols();
for (size_t n=1; n <= nsd; ++n) {
double conv = 0.0;
for (size_t k = 1; k<= nsd; ++k)
conv += U[k-1]*dUdX(n, k);
for (size_t i=1; i <= fe.basis(n).size(); ++i)
EV[n](i) -= scale * conv * fe.basis(n)(i) * fe.detJxW;
}
}
void CompatibleOperators::Residual::Laplacian(Vectors& EV,
const FiniteElement& fe,
const Tensor& dUdX,

23
Apps/Common/CompatibleOperators.h
View File

@ -37,6 +37,16 @@ public:
static void Advection(std::vector<Matrix>& EM, const FiniteElement& fe,
const Vec3& AC, double scale=1.0);
//! \brief Compute a (nonlinear) convection term.
//! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for
//! \param[in] U Advecting field
//! \param[in] conservative True to use the conservative formulation
//! \param[in] basis Basis to use
static void Convection(std::vector<Matrix>& EM, const FiniteElement& fe,
const Vec3& U, const Tensor& dUdX, double scale,
bool conservative=false);
//! \brief Compute a gradient term.
//! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for
@ -87,6 +97,19 @@ public:
class Residual
{
public:
//! \brief Compute a convection term in a residual vector.
//! \param EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for
//! \param[in] U Advected field
//! \param[in] dUdX Advected field gradient
//! \param[in] UC Advecting field
//! \param[in] scale Scaling factor for contribution
//! \param[in] conservative True to use the conservative formulation
//! \param[in] basis Basis to use
static void Convection(Vectors& EV, const FiniteElement& fe,
const Vec3& U, const Tensor& dUdX, const Vec3& UC,
double scale, bool conservative=false);
//! \brief Compute a laplacian term in a residual vector.
//! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for