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

added: conservative form support in EqualOrderOperators::Weak::Advection

Browse Source
This commit is contained in:
Arne Morten Kvarving 2020-04-03 12:06:48 +02:00
parent f5485633ab
commit 00f9fd8a03
2 changed files with 17 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Apps/Common/EqualOrderOperators.C
View File

@ -54,17 +54,24 @@ static void DivGrad(Matrix& EM, const FiniteElement& fe,
void EqualOrderOperators::Weak::Advection(Matrix& EM, const FiniteElement& fe,
const Vec3& AC,
double scale, int basis)
double scale,
WeakOperators::ConvectionForm form,
int basis)
{
Matrix C(fe.basis(basis).size(), fe.basis(basis).size());
size_t ncmp = EM.rows() / C.rows();
for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
for (size_t j = 1; j <= fe.basis(basis).size(); ++j) {
// Sum convection for each direction
for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
C(i,j) += AC[k-1]*fe.grad(basis)(j,k);
C(i,j) *= scale*fe.basis(basis)(i)*fe.detJxW;
if (form == WeakOperators::CONVECTIVE) {
// Sum convection for each direction
for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
C(i,j) += AC[k-1] * fe.grad(basis)(j,k);
C(i,j) *= scale * fe.basis(basis)(i) * fe.detJxW;
} else if (form == WeakOperators::CONSERVATIVE) {
for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
C(i,j) -= AC[k-1] * fe.grad(basis)(i,k);
C(i,j) *= scale * fe.basis(basis)(j) * fe.detJxW;
}
}
}
addComponents(EM, C, ncmp, ncmp, 0);

5
Apps/Common/EqualOrderOperators.h
View File

@ -45,9 +45,12 @@ public:
//! \param[in] fe The finite element to evaluate for
//! \param[in] AC Advecting field
//! \param[in] scale Scaling factor for contribution
//! \param[in] form Formulation to use for advection
//! \param[in] basis Basis to use
static void Advection(Matrix& EM, const FiniteElement& fe,
const Vec3& AC, double scale=1.0, int basis=1);
const Vec3& AC, double scale=1.0,
WeakOperators::ConvectionForm form = WeakOperators::CONVECTIVE,
int basis=1);
//! \brief Compute a (nonlinear) convection term.
//! \param[out] EM The element matrix to add contribution to