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

simplify WeakOperators and ResidualOperators interfaces

Browse Source 
assumes that everything is always split by field. thus
we do not need the cmp/nf/scmp mess.
This commit is contained in:
Arne Morten Kvarving 2016-10-06 12:36:35 +02:00
parent 7275a4fecc
commit 5dac04e420
4 changed files with 218 additions and 247 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
Apps/Common/ResidualOperators.h
View File

@ -21,27 +21,26 @@ class Vec3;
namespace ResidualOperators namespace ResidualOperators
{ {
//! \brief Compute a convection term in a residual vector. //! \brief Compute a convection term in a residual vector.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] U Advected field. //! \param[in] U Advected field
//! \param[in] dUdX Advected field gradient. //! \param[in] dUdX Advected field gradient
//! \param[in] UC Advecting field. //! \param[in] UC Advecting field
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] conservative True to use the conservative formulation
//! \param[in] nf Number of fields in basis. //! \param[in] basis Basis to use
//! \param[in] conservative True to use the conservative formulation.
template<class T> template<class T>
void Convection(Vector& EV, const FiniteElement& fe, void Convection(Vector& EV, const FiniteElement& fe,
const Vec3& U, const T& dUdX, const Vec3& UC, const Vec3& U, const T& dUdX, const Vec3& UC,
double scale, size_t cmp, size_t nf, double scale, bool conservative=false, int basis=1)
bool conservative)
{ {
size_t cmp = EV.size() / fe.basis(basis).size();
if (conservative) { if (conservative) {
for (size_t i = 1;i <= fe.N.size();i++) for (size_t i = 1;i <= fe.basis(basis).size();i++)
for (size_t k = 1;k <= cmp;k++) for (size_t k = 1;k <= cmp;k++)
for (size_t l = 1;l <= cmp;l++) for (size_t l = 1;l <= cmp;l++)
// Convection // Convection
EV((i-1)*nf + k) += scale*U[k-1]*UC[l-1]*fe.dNdX(i,l)*fe.detJxW; EV((i-1)*cmp + k) += scale*U[k-1]*UC[l-1]*fe.grad(basis)(i,l)*fe.detJxW;
} }
else { else {
for (size_t k = 1;k <= cmp;k++) { for (size_t k = 1;k <= cmp;k++) {
@ -51,65 +50,63 @@ namespace ResidualOperators
conv += UC[l-1]*dUdX(k,l); conv += UC[l-1]*dUdX(k,l);
conv *= scale*fe.detJxW; conv *= scale*fe.detJxW;
for (size_t i = 1;i <= fe.N.size();i++) for (size_t i = 1;i <= fe.basis(basis).size();i++)
EV((i-1)*nf + k) -= conv*fe.N(i); EV((i-1)*cmp + k) -= conv*fe.basis(basis)(i);
} }
} }
} }
//! \brief Compute a divergence term in a residual vector. //! \brief Compute a divergence term in a residual vector.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] dUdX Gradient of field. //! \param[in] dUdX Gradient of field
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] nf Number of fields in basis. //! \param[in] basis Basis to use
template<class T> template<class T>
void Divergence(Vector& EV, const FiniteElement& fe, void Divergence(Vector& EV, const FiniteElement& fe,
const T& dUdX, double scale=1.0, const T& dUdX, double scale=1.0,
size_t cmp=1, size_t nf=1, size_t basis=1, size_t ofs=0) size_t basis=1)
{ {
for (size_t i = 1; i <= fe.basis(basis).size(); ++i) { for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
double div=0.0; double div=0.0;
for (size_t k = 1; k <= fe.grad(basis).cols(); ++k) for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
div += dUdX(k,k); div += dUdX(k,k);
EV(ofs+(i-1)*nf+cmp) += scale*div*fe.basis(basis)(i)*fe.detJxW; EV(i) += scale*div*fe.basis(basis)(i)*fe.detJxW;
} }
} }
//! \brief Compute a laplacian term in a residual vector. //! \brief Compute a laplacian term in a residual vector.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] dUdX Current solution gradient. //! \param[in] dUdX Current solution gradient
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] stress Whether to add extra stress formulation terms
//! \param[in] nf Number of fields in basis. //! \param[in] basis Basis to use
//! \param[in] stress Whether to add extra stress formulation terms.
//! \param[in] scmp Starting component.
template<class T> template<class T>
void Laplacian(Vector& EV, const FiniteElement& fe, void Laplacian(Vector& EV, const FiniteElement& fe,
const T& dUdX, double scale=1.0, const T& dUdX, double scale=1.0,
size_t cmp=1, size_t nf=1, bool stress=false, int basis=1)
bool stress=false, size_t scmp=0)
{ {
for (size_t i = 1;i <= fe.N.size();i++) { size_t cmp = EV.size() / fe.basis(basis).size();
for (size_t i = 1;i <= fe.basis(basis).size();i++) {
for (size_t k = 1;k <= cmp;k++) { for (size_t k = 1;k <= cmp;k++) {
double diff = 0.0; double diff = 0.0;
for (size_t l = 1;l <= cmp;l++) for (size_t l = 1;l <= cmp;l++)
diff += dUdX(k,l)*fe.dNdX(i,l); diff += dUdX(k,l)*fe.grad(basis)(i,l);
diff *= scale*fe.detJxW; diff *= scale*fe.detJxW;
// Add negative residual to rhs of momentum equation // Add negative residual to rhs of momentum equation
EV((i-1)*nf + k+scmp) -= diff; EV((i-1)*cmp + k) -= diff;
} }
} }
// Use stress formulation // Use stress formulation
if (stress) { if (stress) {
for (size_t i = 1;i <= fe.N.size();i++) for (size_t i = 1;i <= fe.basis(basis).size();i++)
for (size_t k = 1;k <= cmp;k++) for (size_t k = 1;k <= cmp;k++)
for (size_t l = 1;l <= cmp;l++) for (size_t l = 1;l <= cmp;l++)
// Diffusion // Diffusion
EV((i-1)*nf + k) -= scale*dUdX(l,k)*fe.dNdX(i,l)*fe.detJxW; EV((i-1)*cmp + k) -= scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW;
} }
} }
} }

109
Apps/Common/Test/TestWeakOperators.C
View File

@ -1,6 +1,6 @@
//============================================================================== //==============================================================================
//! //!
//! \file TestMeshUtils.C //! \file TestWeakoperators.C
//! //!
//! \date Feb 16 2015 //! \date Feb 16 2015
//! //!
@ -51,13 +51,10 @@ TEST(TestWeakOperators, Advection)
Matrix EM_vec(2*2, 2*2); Matrix EM_vec(2*2, 2*2);
U[0] = 3.0; U[1] = 4.0; U[0] = 3.0; U[1] = 4.0;
// First component only WeakOperators::Advection(EM_vec, fe, U, 2.0);
WeakOperators::Advection(EM_vec, fe, U, 1.0, 1, 2); const DoubleVec EM_vec_ref = {{30.0, 0.0, 44.0, 0.0},
// Second component only
WeakOperators::Advection(EM_vec, fe, U, 2.0, 1, 2, 1);
const DoubleVec EM_vec_ref = {{15.0, 0.0, 22.0, 0.0},
{ 0.0, 30.0, 0.0, 44.0}, { 0.0, 30.0, 0.0, 44.0},
{30.0, 0.0, 44.0, 0.0}, {60.0, 0.0, 88.0, 0.0},
{ 0.0, 60.0, 0.0, 88.0}}; { 0.0, 60.0, 0.0, 88.0}};
check_matrix_equal(EM_vec, EM_vec_ref); check_matrix_equal(EM_vec, EM_vec_ref);
} }
@ -67,12 +64,9 @@ TEST(TestWeakOperators, Divergence)
{ {
FiniteElement fe = getFE(); FiniteElement fe = getFE();
Matrix EM_scalar(2*2, 2*2); Matrix EM_scalar(2, 2*2);
// single component + pressure WeakOperators::Divergence(EM_scalar, fe);
WeakOperators::Divergence(EM_scalar, fe, 2); const DoubleVec EM_scalar_ref = {{1.0, 3.0, 2.0, 4.0},
const DoubleVec EM_scalar_ref = {{0.0, 0.0, 0.0, 0.0},
{1.0, 3.0, 2.0, 4.0},
{0.0, 0.0, 0.0, 0.0},
{2.0, 6.0, 4.0, 8.0}}; {2.0, 6.0, 4.0, 8.0}};
check_matrix_equal(EM_scalar, EM_scalar_ref); check_matrix_equal(EM_scalar, EM_scalar_ref);
@ -91,21 +85,21 @@ TEST(TestWeakOperators, Gradient)
{ {
FiniteElement fe = getFE(); FiniteElement fe = getFE();
Matrix EM_scalar(2*2, 2*2); Matrix EM_scalar(2*2,2);
// single component + pressure // single component + pressure
WeakOperators::Gradient(EM_scalar, fe, 2); WeakOperators::Gradient(EM_scalar, fe);
const DoubleVec EM_scalar_ref = {{0.0,-1.0, 0.0,-2.0}, const DoubleVec EM_scalar_ref = {{-1.0,-2.0},
{0.0,-3.0, 0.0,-6.0}, {-3.0,-6.0},
{0.0,-2.0, 0.0,-4.0}, {-2.0,-4.0},
{0.0,-4.0, 0.0,-8.0}}; {-4.0,-8.0}};
check_matrix_equal(EM_scalar, EM_scalar_ref); check_matrix_equal(EM_scalar, EM_scalar_ref);
Vector EV_scalar(4); Vector EV_scalar(4);
WeakOperators::Gradient(EV_scalar, fe); WeakOperators::Gradient(EV_scalar, fe);
ASSERT_NEAR(EV_scalar(1), 1.0, 1e-13); ASSERT_NEAR(EV_scalar(1), fe.dNdX(1,1), 1e-13);
ASSERT_NEAR(EV_scalar(2), 5.0, 1e-13); ASSERT_NEAR(EV_scalar(2), fe.dNdX(1,2), 1e-13);
ASSERT_NEAR(EV_scalar(3), 4.0, 1e-13); ASSERT_NEAR(EV_scalar(3), fe.dNdX(2,1), 1e-13);
ASSERT_NEAR(EV_scalar(4), 0.0, 1e-13); ASSERT_NEAR(EV_scalar(4), fe.dNdX(2,2), 1e-13);
} }
@ -115,7 +109,7 @@ TEST(TestWeakOperators, Laplacian)
// single scalar block // single scalar block
Matrix EM_scalar(2,2); Matrix EM_scalar(2,2);
WeakOperators::Laplacian(EM_scalar, fe, 1.0, 1, 1); WeakOperators::Laplacian(EM_scalar, fe);
const DoubleVec EM_scalar_ref = {{10.0, 14.0}, const DoubleVec EM_scalar_ref = {{10.0, 14.0},
{14.0, 20.0}}; {14.0, 20.0}};
@ -123,22 +117,19 @@ TEST(TestWeakOperators, Laplacian)
check_matrix_equal(EM_scalar, EM_scalar_ref); check_matrix_equal(EM_scalar, EM_scalar_ref);
// multiple (2) blocks in 3 component element matrix // multiple (2) blocks in 3 component element matrix
Matrix EM_vec(2*3,2*3); Matrix EM_multi(2*2,2*2);
WeakOperators::Laplacian(EM_vec, fe, 1.0, 2, 3); WeakOperators::Laplacian(EM_multi, fe, 1.0);
// last block separately
WeakOperators::Laplacian(EM_vec, fe, 2.0, 1, 3, false, 2);
const DoubleVec EM_vec_ref = {{10.0, 0.0, 0.0, 14.0, 0.0, 0.0},
{ 0.0, 10.0, 0.0, 0.0, 14.0, 0.0},
{ 0.0, 0.0, 20.0, 0.0, 0.0, 28.0},
{14.0, 0.0, 0.0, 20.0, 0.0, 0.0},
{ 0.0, 14.0, 0.0, 0.0, 20.0, 0.0},
{ 0.0, 0.0, 28.0, 0.0, 0.0, 40.0}};
check_matrix_equal(EM_vec, EM_vec_ref); const DoubleVec EM_multi_ref = {{10.0, 0.0, 14.0, 0.0},
{ 0.0, 10.0, 0.0, 14.0},
{14.0, 0.0, 20.0, 0.0},
{0.0, 14.0, 0.0, 20.0}};
check_matrix_equal(EM_multi, EM_multi_ref);
// stress formulation // stress formulation
Matrix EM_stress(2*2,2*2); Matrix EM_stress(2*2,2*2);
WeakOperators::Laplacian(EM_stress, fe, 1.0, 2, 2, true); WeakOperators::Laplacian(EM_stress, fe, 1.0, true);
const DoubleVec EM_stress_ref = {{11.0, 3.0, 16.0, 6.0}, const DoubleVec EM_stress_ref = {{11.0, 3.0, 16.0, 6.0},
{ 3.0, 19.0, 4.0, 26.0}, { 3.0, 19.0, 4.0, 26.0},
{16.0, 4.0, 24.0, 8.0}, {16.0, 4.0, 24.0, 8.0},
@ -164,18 +155,13 @@ TEST(TestWeakOperators, Mass)
{4.0, 8.0}}; {4.0, 8.0}};
check_matrix_equal(EM_scalar, EM_scalar_ref); check_matrix_equal(EM_scalar, EM_scalar_ref);
Matrix EM_vec(2*3,2*3); Matrix EM_vec(2*2,2*2);
// Two first components WeakOperators::Mass(EM_vec, fe);
WeakOperators::Mass(EM_vec, fe, 1.0, 2, 3);
// Last component const DoubleVec EM_vec_ref = {{1.0, 0.0, 2.0, 0.0},
WeakOperators::Mass(EM_vec, fe, 1.0, 1, 3, 2); {0.0, 1.0, 0.0, 2.0},
const DoubleVec EM_vec_ref = {{1.0, 0.0, 0.0, 2.0, 0.0, 0.0}, {2.0, 0.0, 4.0, 0.0},
{0.0, 1.0, 0.0, 0.0, 2.0, 0.0}, {0.0, 2.0, 0.0, 4.0}};
{0.0, 0.0, 1.0, 0.0, 0.0, 2.0},
{2.0, 0.0, 0.0, 4.0, 0.0, 0.0},
{0.0, 2.0, 0.0, 0.0, 4.0, 0.0},
{0.0, 0.0, 2.0, 0.0, 0.0, 4.0}};
check_matrix_equal(EM_vec, EM_vec_ref); check_matrix_equal(EM_vec, EM_vec_ref);
} }
@ -189,15 +175,30 @@ TEST(TestWeakOperators, Source)
ASSERT_NEAR(EV_scalar(1), 2.0, 1e-13); ASSERT_NEAR(EV_scalar(1), 2.0, 1e-13);
ASSERT_NEAR(EV_scalar(2), 4.0, 1e-13); ASSERT_NEAR(EV_scalar(2), 4.0, 1e-13);
Vector EV_vec(2*3); Vector EV_vec(4);
Vec3 f; Vec3 f;
f[0] = 1.0; f[1] = 2.0; f[0] = 1.0; f[1] = 2.0;
WeakOperators::Source(EV_vec, fe, f, 1.0, 2, 3); WeakOperators::Source(EV_vec, fe, f, 1.0);
WeakOperators::Source(EV_vec, fe, 5.0, 1, 3, 2);
ASSERT_NEAR(EV_vec(1), 1.0, 1e-13); ASSERT_NEAR(EV_vec(1), 1.0, 1e-13);
ASSERT_NEAR(EV_vec(2), 2.0, 1e-13); ASSERT_NEAR(EV_vec(2), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(3), 5.0, 1e-13); ASSERT_NEAR(EV_vec(3), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(4), 2.0, 1e-13); ASSERT_NEAR(EV_vec(4), 4.0, 1e-13);
ASSERT_NEAR(EV_vec(5), 4.0, 1e-13); EV_vec.fill(0.0);
ASSERT_NEAR(EV_vec(6), 10.0, 1e-13); WeakOperators::Source(EV_vec, fe, 2.0, 0);
ASSERT_NEAR(EV_vec(1), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(2), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(3), 4.0, 1e-13);
ASSERT_NEAR(EV_vec(4), 4.0, 1e-13);
EV_vec.fill(0.0);
WeakOperators::Source(EV_vec, fe, 2.0, 1);
ASSERT_NEAR(EV_vec(1), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(2), 0.0, 1e-13);
ASSERT_NEAR(EV_vec(3), 4.0, 1e-13);
ASSERT_NEAR(EV_vec(4), 0.0, 1e-13);
EV_vec.fill(0.0);
WeakOperators::Source(EV_vec, fe, 2.0, 2);
ASSERT_NEAR(EV_vec(1), 0.0, 1e-13);
ASSERT_NEAR(EV_vec(2), 2.0, 1e-13);
ASSERT_NEAR(EV_vec(3), 0.0, 1e-13);
ASSERT_NEAR(EV_vec(4), 4.0, 1e-13);
} }

114
Apps/Common/WeakOperators.C
View File

@ -37,136 +37,134 @@ namespace WeakOperators {
//! \brief Helper applying a divergence (1) or a gradient (2) or both (0) operation //! \brief Helper applying a divergence (1) or a gradient (2) or both (0) operation
template<int Operation> template<int Operation>
static void DivGrad(Matrix& EM, const FiniteElement& fe, static void DivGrad(Matrix& EM, const FiniteElement& fe,
double scale=1.0, size_t nf=1) double scale, int basis, int tbasis)
{ {
for (size_t i = 1; i <= fe.N.size();++i) size_t nsd = fe.grad(basis).cols();
for (size_t j = 1; j <= fe.N.size();++j) for (size_t i = 1; i <= fe.basis(tbasis).size();++i)
for (size_t k = 1; k <= fe.dNdX.cols(); ++k) { for (size_t j = 1; j <= fe.basis(basis).size();++j)
double div = fe.N(j)*fe.dNdX(i,k)*fe.detJxW; for (size_t k = 1; k <= nsd; ++k) {
if (Operation == 2 || Operation == 0) double div = fe.basis(basis)(j)*fe.grad(tbasis)(i,k)*fe.detJxW;
EM(nf*(i-1)+k,nf*j) += -scale*div; if (Operation == 2)
if (Operation == 1 || Operation == 0) EM((i-1)*nsd+k,j) += -scale*div;
EM(nf*j, nf*(i-1)+k) += scale*div; if (Operation == 1)
EM(j, (i-1)*nsd+k) += scale*div;
} }
} }
void Advection(Matrix& EM, const FiniteElement& fe, void Advection(Matrix& EM, const FiniteElement& fe,
const Vec3& AC, double scale, const Vec3& AC, double scale, int basis)
size_t cmp, size_t nf, size_t scmp)
{ {
Matrix C(fe.N.size(), fe.N.size()); Matrix C(fe.basis(basis).size(), fe.basis(basis).size());
for (size_t i = 1; i <= fe.N.size(); ++i) { size_t ncmp = EM.rows() / C.rows();
for (size_t j = 1; j <= fe.N.size(); ++j) { 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 // Sum convection for each direction
for (size_t k = 1; k <= fe.dNdX.cols(); ++k) for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
C(i,j) += AC[k-1]*fe.dNdX(j,k); C(i,j) += AC[k-1]*fe.dNdX(j,k);
C(i,j) *= scale*fe.N(i)*fe.detJxW; C(i,j) *= scale*fe.N(i)*fe.detJxW;
} }
} }
addComponents(EM, C, cmp, nf, scmp); addComponents(EM, C, ncmp, ncmp, 0);
} }
void Divergence(Matrix& EM, const FiniteElement& fe, void Divergence(Matrix& EM, const FiniteElement& fe, double scale,
size_t nf, double scale) int basis, int tbasis)
{ {
DivGrad<1>(EM,fe,scale,nf); DivGrad<1>(EM,fe,scale,basis,tbasis);
} }
void PressureDiv(Matrix& EM, const FiniteElement& fe, void Gradient(Matrix& EM, const FiniteElement& fe, double scale,
size_t nf, double scale) int basis, int tbasis)
{ {
DivGrad<0>(EM,fe,scale,nf); DivGrad<2>(EM,fe,scale,basis,tbasis);
}
void Gradient(Matrix& EM, const FiniteElement& fe,
size_t nf, double scale)
{
DivGrad<2>(EM,fe,scale,nf);
} }
void Divergence(Vector& EV, const FiniteElement& fe, void Divergence(Vector& EV, const FiniteElement& fe,
const Vec3& D, double scale, size_t cmp, size_t nf) const Vec3& D, double scale, int basis)
{ {
for (size_t i = 1; i <= fe.N.size(); ++i) { for (size_t i = 1; i <= fe.basis(basis).size(); ++i) {
double div=0.0; double div=0.0;
for (size_t k = 1; k <= fe.dNdX.cols(); ++k) for (size_t k = 1; k <= fe.grad(basis).cols(); ++k)
div += D[k-1]*fe.dNdX(i,k); div += D[k-1]*fe.grad(basis)(i,k);
EV((i-1)*nf+cmp) += scale*div*fe.detJxW; EV(i) += scale*div*fe.detJxW;
} }
} }
void Gradient(Vector& EV, const FiniteElement& fe, void Gradient(Vector& EV, const FiniteElement& fe,
double scale, size_t nf) double scale, int basis)
{ {
for (size_t i = 1; i <= fe.N.size(); ++i) size_t nsd = fe.grad(basis).cols();
for (size_t k = 1; k <= fe.dNdX.cols(); ++k) for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
EV((i-1)*nf+k) += scale*fe.dNdX(i,k)*fe.detJxW; for (size_t k = 1; k <= nsd; ++k)
EV((i-1)*nsd+k) += scale*fe.grad(basis)(i,k)*fe.detJxW;
} }
void Laplacian(Matrix& EM, const FiniteElement& fe, void Laplacian(Matrix& EM, const FiniteElement& fe,
double scale, size_t cmp, size_t nf, double scale, bool stress, int basis)
bool stress, size_t scmp, unsigned char basis)
{ {
size_t cmp = EM.rows() / fe.basis(basis).size();
Matrix A; Matrix A;
A.multiply(fe.grad(basis),fe.grad(basis),false,true); A.multiply(fe.grad(basis),fe.grad(basis),false,true);
A *= scale*fe.detJxW; A *= scale*fe.detJxW;
addComponents(EM, A, cmp, nf, scmp); addComponents(EM, A, cmp, cmp, 0);
if (stress) if (stress)
for (size_t i = 1; i <= fe.N.size(); i++) for (size_t i = 1; i <= fe.basis(basis).size(); i++)
for (size_t j = 1; j <= fe.N.size(); j++) for (size_t j = 1; j <= fe.basis(basis).size(); j++)
for (size_t k = 1; k <= cmp; k++) for (size_t k = 1; k <= cmp; k++)
for (size_t l = 1; l <= cmp; l++) for (size_t l = 1; l <= cmp; l++)
EM(nf*(j-1)+k+scmp,nf*(i-1)+l+scmp) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW; EM(cmp*(j-1)+k,cmp*(i-1)+l) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW;
} }
void LaplacianCoeff(Matrix& EM, const Matrix& K, void LaplacianCoeff(Matrix& EM, const Matrix& K,
const FiniteElement& fe, const FiniteElement& fe,
double scale) double scale, int basis)
{ {
Matrix KB; Matrix KB;
KB.multiply(K,fe.dNdX,false,true).multiply(scale*fe.detJxW); KB.multiply(K,fe.grad(basis),false,true).multiply(scale*fe.detJxW);
EM.multiply(fe.dNdX,KB,false,false,true); EM.multiply(fe.grad(basis),KB,false,false,true);
} }
void Mass(Matrix& EM, const FiniteElement& fe, void Mass(Matrix& EM, const FiniteElement& fe,
double scale, size_t cmp, size_t nf, size_t scmp, double scale, int basis)
unsigned char basis)
{ {
size_t ncmp = EM.rows()/fe.basis(basis).size();
Matrix A; Matrix A;
A.outer_product(fe.basis(basis),fe.basis(basis),false); A.outer_product(fe.basis(basis),fe.basis(basis),false);
A *= scale*fe.detJxW; A *= scale*fe.detJxW;
addComponents(EM, A, cmp, nf, scmp); addComponents(EM, A, ncmp, ncmp, 0);
} }
void Source(Vector& EV, const FiniteElement& fe, void Source(Vector& EV, const FiniteElement& fe,
double scale, size_t cmp, size_t nf, size_t scmp, double scale, int cmp, int basis)
unsigned char basis)
{ {
if (cmp == 1 && nf == 1) size_t ncmp = EV.size() / fe.basis(basis).size();
if (cmp == 1 && ncmp == 1)
EV.add(fe.basis(basis), scale*fe.detJxW); EV.add(fe.basis(basis), scale*fe.detJxW);
else { else {
for (size_t i = 1; i <= fe.basis(basis).size(); ++i) for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
for (size_t k = 1; k <= cmp; ++k) for (size_t k = (cmp == 0 ? 1: cmp);
EV(nf*(i-1)+k+scmp) += scale*fe.basis(basis)(i)*fe.detJxW; k <= (cmp == 0 ? ncmp : cmp); ++k)
EV(ncmp*(i-1)+k) += scale*fe.basis(basis)(i)*fe.detJxW;
} }
} }
void Source(Vector& EV, const FiniteElement& fe, void Source(Vector& EV, const FiniteElement& fe,
const Vec3& f, double scale, size_t cmp, size_t nf, size_t scmp) const Vec3& f, double scale, int basis)
{ {
for (size_t i=0;i<cmp;++i) size_t cmp = EV.size() / fe.basis(basis).size();
Source(EV, fe, f[i+scmp]*scale, 1, nf, scmp+i); for (size_t i = 1; i <= fe.basis(basis).size(); ++i)
for (size_t k = 1; k <= cmp; ++k)
EV(cmp*(i-1)+k) += scale*f[k-1]*fe.basis(basis)(i)*fe.detJxW;
} }
} }

171
Apps/Common/WeakOperators.h
View File

@ -21,150 +21,125 @@ class Vec3;
namespace WeakOperators namespace WeakOperators
{ {
//! \brief Compute an advection term. //! \brief Compute an advection term.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] AC Advecting field. //! \param[in] AC Advecting field
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] basis Basis to use
//! \param[in] nf Number of fields in basis.
//! \param[in] scmp Starting component.
void Advection(Matrix& EM, const FiniteElement& fe, void Advection(Matrix& EM, const FiniteElement& fe,
const Vec3& AC, double scale=1.0, const Vec3& AC, double scale=1.0, int basis=1);
size_t cmp=1, size_t nf=1, size_t scmp=0);
//! \brief Compute a (nonlinear) convection term. //! \brief Compute a (nonlinear) convection term.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] U Advecting field. //! \param[in] U Advecting field
//! \param[in] scale Scaling factor for contribution. //! \param[in] conservative True to use the conservative formulation
//! \param[in] cmp Number of components to add. //! \param[in] basis Basis to use
//! \param[in] nf Number of fields in basis.
//! \param[in] conservative True to use the conservative formulation.
template<class T> template<class T>
void Convection(Matrix& EM, const FiniteElement& fe, void Convection(Matrix& EM, const FiniteElement& fe,
const Vec3& U, const T& dUdX, double scale, const Vec3& U, const T& dUdX, double scale,
size_t cmp, size_t nf, bool conservative) bool conservative=false, int basis=1)
{ {
size_t cmp = EM.rows() / fe.basis(basis).size();
if (conservative) { if (conservative) {
Advection(EM, fe, U, -scale, cmp, nf); Advection(EM, fe, U, -scale, basis);
for (size_t i = 1;i <= fe.N.size();i++) for (size_t i = 1;i <= fe.basis(basis).size();i++)
for (size_t j = 1;j <= fe.N.size();j++) { for (size_t j = 1;j <= fe.basis(basis).size();j++) {
for (size_t k = 1;k <= cmp;k++) { for (size_t k = 1;k <= cmp;k++) {
for (size_t l = 1;l <= cmp;l++) for (size_t l = 1;l <= cmp;l++)
EM((j-1)*nf+l,(i-1)*nf+k) -= scale*U[l-1]*fe.N(i)*fe.dNdX(j,k)*fe.detJxW; EM((j-1)*cmp+l,(i-1)*cmp+k) -= scale*U[l-1]*fe.basis(basis)(i)*fe.grad(basis)(j,k)*fe.detJxW;
} }
} }
} }
else { else {
Advection(EM, fe, U, scale, cmp, nf); Advection(EM, fe, U, scale, basis);
for (size_t i = 1;i <= fe.N.size();i++) for (size_t i = 1;i <= fe.basis(basis).size();i++)
for (size_t j = 1;j <= fe.N.size();j++) { for (size_t j = 1;j <= fe.basis(basis).size();j++) {
for (size_t k = 1;k <= cmp;k++) { for (size_t k = 1;k <= cmp;k++) {
for (size_t l = 1;l <= cmp;l++) for (size_t l = 1;l <= cmp;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += scale*dUdX(l,k)*fe.N(j)*fe.N(i)*fe.detJxW; EM((j-1)*cmp+l,(i-1)*cmp+k) += scale*dUdX(l,k)*fe.basis(basis)(j)*fe.basis(basis)(i)*fe.detJxW;
} }
} }
} }
} }
//! \brief Compute a divergence term. //! \brief Compute a divergence term.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] nf Number of fields in basis. //! \param[in] scale Scaling factor for contribution
//! \param[in] scale Scaling factor for contribution. //! \param[in] basis Basis for field
//! \param[in] tbasis Test function basis
void Divergence(Matrix& EM, const FiniteElement& fe, void Divergence(Matrix& EM, const FiniteElement& fe,
size_t nf, double scale=1.0); double scale=1.0, int basis=1, int tbasis=1);
//! \brief Compute a divergence term. //! \brief Compute a divergence term.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] D Divergence of field. //! \param[in] D Divergence of field
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] nf Number of fields in basis. //! \param[in] basis Test function basis
//
void Divergence(Vector& EV, const FiniteElement& fe, void Divergence(Vector& EV, const FiniteElement& fe,
const Vec3& D, double scale=1.0, const Vec3& D, double scale=1.0, int basis=1);
size_t cmp=1, size_t nf=1);
//! \brief Compute a divergence/gradient term. //! \brief Compute a gradient term for a (potentially mixed) vector/scalar field.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] nf Number of fields in basis. //! \param[in] scale Scaling factor for contribution
//! \param[in] scale Scaling factor for contribution. //! \param[in] basis Basis for field
void PressureDiv(Matrix& EM, const FiniteElement& fe, //! \param[in] tbasis Test function basis
size_t nf, double scale=1.0);
//! \brief Compute a gradient term.
//! \param[out] EM The element matrix to add contribution to.
//! \param[in] fe The finite element to evaluate for.
//! \param[in] nf Number of fields in basis.
//! \param[in] scale Scaling factor for contribution.
void Gradient(Matrix& EM, const FiniteElement& fe, void Gradient(Matrix& EM, const FiniteElement& fe,
size_t nf, double scale=1.0); double scale=1.0, int basis=1, int tbasis=1);
//! \brief Compute a gradient term. //! \brief Compute a gradient term.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] nf Number of fields in basis. //! \param[in] tbasis Test function basis
void Gradient(Vector& EV, const FiniteElement& fe, void Gradient(Vector& EV, const FiniteElement& fe,
double scale=1.0, size_t nf=1); double scale=1.0, int basis=1);
//! \brief Compute a laplacian. //! \brief Compute a laplacian.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] stress Whether to add extra stress formulation terms
//! \param[in] nf Number of fields in basis. //! \param[in] basis Basis to use
//! \param[in] stress Whether to add extra stress formulation terms.
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Laplacian(Matrix& EM, const FiniteElement& fe, void Laplacian(Matrix& EM, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1, double scale=1.0, bool stress=false, int basis=1);
bool stress=false, size_t scmp=0, unsigned char basis=1);
//! \brief Compute a heteregenous coefficient laplacian. //! \brief Compute a heteregenous coefficient laplacian.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[out] K The coefficient matrix. //! \param[out] K The coefficient matrix
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe, void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe,
double scale=1.0); double scale=1.0, int basis=1);
//! \brief Compute a mass term. //! \brief Compute a mass term.
//! \param[out] EM The element matrix to add contribution to. //! \param[out] EM The element matrix to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] basis Basis to use
//! \param[in] nf Number of fields in basis.
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Mass(Matrix& EM, const FiniteElement& fe, void Mass(Matrix& EM, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, double scale=1.0, int basis=1);
unsigned char basis=1);
//! \brief Compute a source term. //! \brief Compute a source term.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] basis Basis to use
//! \param[in] nf Number of fields in basis. //! \param[in] cmp Component to add (0 for all)
//! \param[in] scmp Starting component.
//! \param[in] basis Basis to use.
void Source(Vector& EV, const FiniteElement& fe, void Source(Vector& EV, const FiniteElement& fe,
double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, double scale=1.0, int cmp=1, int basis=1);
unsigned char basis=1);
//! \brief Compute a vector-source term. //! \brief Compute a vector-source term.
//! \param[out] EV The element vector to add contribution to. //! \param[out] EV The element vector to add contribution to
//! \param[in] fe The finite element to evaluate for. //! \param[in] fe The finite element to evaluate for
//! \param[in] scale Vector with contributions. //! \param[in] scale Vector with contributions
//! \param[in] scale Scaling factor for contribution. //! \param[in] scale Scaling factor for contribution
//! \param[in] cmp Number of components to add. //! \param[in] basis Basis to use
//! \param[in] nf Number of fields in basis.
//! \param[in] scmp Starting component.
void Source(Vector& EV, const FiniteElement& fe, void Source(Vector& EV, const FiniteElement& fe,
const Vec3& f, double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0); const Vec3& f, double scale=1.0, int basis=1);
} }
#endif #endif