diff --git a/Apps/Common/ResidualOperators.h b/Apps/Common/ResidualOperators.h index 8f21985d..245e0299 100644 --- a/Apps/Common/ResidualOperators.h +++ b/Apps/Common/ResidualOperators.h @@ -21,27 +21,26 @@ class Vec3; namespace ResidualOperators { //! \brief Compute a convection term in a residual vector. - //! \param[out] 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] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] conservative True to use the conservative formulation. + //! \param[out] 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 template void Convection(Vector& EV, const FiniteElement& fe, const Vec3& U, const T& dUdX, const Vec3& UC, - double scale, size_t cmp, size_t nf, - bool conservative) + double scale, bool conservative=false, int basis=1) { + size_t cmp = EV.size() / fe.basis(basis).size(); 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 l = 1;l <= cmp;l++) // 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 { for (size_t k = 1;k <= cmp;k++) { @@ -51,65 +50,63 @@ namespace ResidualOperators conv += UC[l-1]*dUdX(k,l); conv *= scale*fe.detJxW; - for (size_t i = 1;i <= fe.N.size();i++) - EV((i-1)*nf + k) -= conv*fe.N(i); + for (size_t i = 1;i <= fe.basis(basis).size();i++) + EV((i-1)*cmp + k) -= conv*fe.basis(basis)(i); } } } //! \brief Compute a divergence term in a residual vector. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] dUdX Gradient of field. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] nf Number of fields in basis. + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] dUdX Gradient of field + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use template void Divergence(Vector& EV, const FiniteElement& fe, 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) { double div=0.0; for (size_t k = 1; k <= fe.grad(basis).cols(); ++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. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] dUdX Current solution gradient. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] stress Whether to add extra stress formulation terms. - //! \param[in] scmp Starting component. + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] dUdX Current solution gradient + //! \param[in] scale Scaling factor for contribution + //! \param[in] stress Whether to add extra stress formulation terms + //! \param[in] basis Basis to use template void Laplacian(Vector& EV, const FiniteElement& fe, const T& dUdX, double scale=1.0, - size_t cmp=1, size_t nf=1, - bool stress=false, size_t scmp=0) + bool stress=false, int basis=1) { - 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++) { double diff = 0.0; 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; // Add negative residual to rhs of momentum equation - EV((i-1)*nf + k+scmp) -= diff; + EV((i-1)*cmp + k) -= diff; } } // Use stress formulation 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 l = 1;l <= cmp;l++) // 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; } } } diff --git a/Apps/Common/Test/TestWeakOperators.C b/Apps/Common/Test/TestWeakOperators.C index 34c748d2..fac7e017 100644 --- a/Apps/Common/Test/TestWeakOperators.C +++ b/Apps/Common/Test/TestWeakOperators.C @@ -1,6 +1,6 @@ //============================================================================== //! -//! \file TestMeshUtils.C +//! \file TestWeakoperators.C //! //! \date Feb 16 2015 //! @@ -51,13 +51,10 @@ TEST(TestWeakOperators, Advection) Matrix EM_vec(2*2, 2*2); U[0] = 3.0; U[1] = 4.0; - // First component only - WeakOperators::Advection(EM_vec, fe, U, 1.0, 1, 2); - // 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}, + WeakOperators::Advection(EM_vec, fe, U, 2.0); + const DoubleVec EM_vec_ref = {{30.0, 0.0, 44.0, 0.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}}; check_matrix_equal(EM_vec, EM_vec_ref); } @@ -67,12 +64,9 @@ TEST(TestWeakOperators, Divergence) { FiniteElement fe = getFE(); - Matrix EM_scalar(2*2, 2*2); - // single component + pressure - WeakOperators::Divergence(EM_scalar, fe, 2); - 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}, + Matrix EM_scalar(2, 2*2); + WeakOperators::Divergence(EM_scalar, fe); + const DoubleVec EM_scalar_ref = {{1.0, 3.0, 2.0, 4.0}, {2.0, 6.0, 4.0, 8.0}}; check_matrix_equal(EM_scalar, EM_scalar_ref); @@ -91,21 +85,21 @@ TEST(TestWeakOperators, Gradient) { FiniteElement fe = getFE(); - Matrix EM_scalar(2*2, 2*2); + Matrix EM_scalar(2*2,2); // single component + pressure - WeakOperators::Gradient(EM_scalar, fe, 2); - const DoubleVec EM_scalar_ref = {{0.0,-1.0, 0.0,-2.0}, - {0.0,-3.0, 0.0,-6.0}, - {0.0,-2.0, 0.0,-4.0}, - {0.0,-4.0, 0.0,-8.0}}; + WeakOperators::Gradient(EM_scalar, fe); + const DoubleVec EM_scalar_ref = {{-1.0,-2.0}, + {-3.0,-6.0}, + {-2.0,-4.0}, + {-4.0,-8.0}}; check_matrix_equal(EM_scalar, EM_scalar_ref); Vector EV_scalar(4); WeakOperators::Gradient(EV_scalar, fe); - ASSERT_NEAR(EV_scalar(1), 1.0, 1e-13); - ASSERT_NEAR(EV_scalar(2), 5.0, 1e-13); - ASSERT_NEAR(EV_scalar(3), 4.0, 1e-13); - ASSERT_NEAR(EV_scalar(4), 0.0, 1e-13); + ASSERT_NEAR(EV_scalar(1), fe.dNdX(1,1), 1e-13); + ASSERT_NEAR(EV_scalar(2), fe.dNdX(1,2), 1e-13); + ASSERT_NEAR(EV_scalar(3), fe.dNdX(2,1), 1e-13); + ASSERT_NEAR(EV_scalar(4), fe.dNdX(2,2), 1e-13); } @@ -115,7 +109,7 @@ TEST(TestWeakOperators, Laplacian) // single scalar block 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}, {14.0, 20.0}}; @@ -123,22 +117,19 @@ TEST(TestWeakOperators, Laplacian) check_matrix_equal(EM_scalar, EM_scalar_ref); // multiple (2) blocks in 3 component element matrix - Matrix EM_vec(2*3,2*3); - WeakOperators::Laplacian(EM_vec, fe, 1.0, 2, 3); - // 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}}; + Matrix EM_multi(2*2,2*2); + WeakOperators::Laplacian(EM_multi, fe, 1.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 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}, { 3.0, 19.0, 4.0, 26.0}, {16.0, 4.0, 24.0, 8.0}, @@ -164,18 +155,13 @@ TEST(TestWeakOperators, Mass) {4.0, 8.0}}; check_matrix_equal(EM_scalar, EM_scalar_ref); - Matrix EM_vec(2*3,2*3); - // Two first components - WeakOperators::Mass(EM_vec, fe, 1.0, 2, 3); + Matrix EM_vec(2*2,2*2); + WeakOperators::Mass(EM_vec, fe); - // Last component - WeakOperators::Mass(EM_vec, fe, 1.0, 1, 3, 2); - const DoubleVec EM_vec_ref = {{1.0, 0.0, 0.0, 2.0, 0.0, 0.0}, - {0.0, 1.0, 0.0, 0.0, 2.0, 0.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}}; + const DoubleVec EM_vec_ref = {{1.0, 0.0, 2.0, 0.0}, + {0.0, 1.0, 0.0, 2.0}, + {2.0, 0.0, 4.0, 0.0}, + {0.0, 2.0, 0.0, 4.0}}; 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(2), 4.0, 1e-13); - Vector EV_vec(2*3); + Vector EV_vec(4); Vec3 f; f[0] = 1.0; f[1] = 2.0; - WeakOperators::Source(EV_vec, fe, f, 1.0, 2, 3); - WeakOperators::Source(EV_vec, fe, 5.0, 1, 3, 2); + WeakOperators::Source(EV_vec, fe, f, 1.0); ASSERT_NEAR(EV_vec(1), 1.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(4), 2.0, 1e-13); - ASSERT_NEAR(EV_vec(5), 4.0, 1e-13); - ASSERT_NEAR(EV_vec(6), 10.0, 1e-13); + ASSERT_NEAR(EV_vec(3), 2.0, 1e-13); + ASSERT_NEAR(EV_vec(4), 4.0, 1e-13); + EV_vec.fill(0.0); + 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); } diff --git a/Apps/Common/WeakOperators.C b/Apps/Common/WeakOperators.C index d60b54d2..2cf505e9 100644 --- a/Apps/Common/WeakOperators.C +++ b/Apps/Common/WeakOperators.C @@ -37,136 +37,134 @@ namespace WeakOperators { //! \brief Helper applying a divergence (1) or a gradient (2) or both (0) operation template 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) - for (size_t j = 1; j <= fe.N.size();++j) - for (size_t k = 1; k <= fe.dNdX.cols(); ++k) { - double div = fe.N(j)*fe.dNdX(i,k)*fe.detJxW; - if (Operation == 2 || Operation == 0) - EM(nf*(i-1)+k,nf*j) += -scale*div; - if (Operation == 1 || Operation == 0) - EM(nf*j, nf*(i-1)+k) += scale*div; + size_t nsd = fe.grad(basis).cols(); + for (size_t i = 1; i <= fe.basis(tbasis).size();++i) + for (size_t j = 1; j <= fe.basis(basis).size();++j) + for (size_t k = 1; k <= nsd; ++k) { + double div = fe.basis(basis)(j)*fe.grad(tbasis)(i,k)*fe.detJxW; + if (Operation == 2) + EM((i-1)*nsd+k,j) += -scale*div; + if (Operation == 1) + EM(j, (i-1)*nsd+k) += scale*div; } } void Advection(Matrix& EM, const FiniteElement& fe, - const Vec3& AC, double scale, - size_t cmp, size_t nf, size_t scmp) + const Vec3& AC, double scale, int basis) { - Matrix C(fe.N.size(), fe.N.size()); - for (size_t i = 1; i <= fe.N.size(); ++i) { - for (size_t j = 1; j <= fe.N.size(); ++j) { + 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.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) *= 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, - size_t nf, double scale) + void Divergence(Matrix& EM, const FiniteElement& fe, 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, - size_t nf, double scale) + void Gradient(Matrix& EM, const FiniteElement& fe, double scale, + int basis, int tbasis) { - DivGrad<0>(EM,fe,scale,nf); - } - - - void Gradient(Matrix& EM, const FiniteElement& fe, - size_t nf, double scale) - { - DivGrad<2>(EM,fe,scale,nf); + DivGrad<2>(EM,fe,scale,basis,tbasis); } 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; - for (size_t k = 1; k <= fe.dNdX.cols(); ++k) - div += D[k-1]*fe.dNdX(i,k); - EV((i-1)*nf+cmp) += scale*div*fe.detJxW; + for (size_t k = 1; k <= fe.grad(basis).cols(); ++k) + div += D[k-1]*fe.grad(basis)(i,k); + EV(i) += scale*div*fe.detJxW; } } 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) - for (size_t k = 1; k <= fe.dNdX.cols(); ++k) - EV((i-1)*nf+k) += scale*fe.dNdX(i,k)*fe.detJxW; + size_t nsd = fe.grad(basis).cols(); + for (size_t i = 1; i <= fe.basis(basis).size(); ++i) + 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, - double scale, size_t cmp, size_t nf, - bool stress, size_t scmp, unsigned char basis) + double scale, bool stress, int basis) { + size_t cmp = EM.rows() / fe.basis(basis).size(); Matrix A; A.multiply(fe.grad(basis),fe.grad(basis),false,true); A *= scale*fe.detJxW; - addComponents(EM, A, cmp, nf, scmp); + addComponents(EM, A, cmp, cmp, 0); if (stress) - for (size_t i = 1; i <= fe.N.size(); i++) - 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++) for (size_t k = 1; k <= cmp; k++) 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, const FiniteElement& fe, - double scale) + double scale, int basis) { Matrix KB; - KB.multiply(K,fe.dNdX,false,true).multiply(scale*fe.detJxW); - EM.multiply(fe.dNdX,KB,false,false,true); + KB.multiply(K,fe.grad(basis),false,true).multiply(scale*fe.detJxW); + EM.multiply(fe.grad(basis),KB,false,false,true); } void Mass(Matrix& EM, const FiniteElement& fe, - double scale, size_t cmp, size_t nf, size_t scmp, - unsigned char basis) + double scale, int basis) { + size_t ncmp = EM.rows()/fe.basis(basis).size(); Matrix A; A.outer_product(fe.basis(basis),fe.basis(basis),false); A *= scale*fe.detJxW; - addComponents(EM, A, cmp, nf, scmp); + addComponents(EM, A, ncmp, ncmp, 0); } void Source(Vector& EV, const FiniteElement& fe, - double scale, size_t cmp, size_t nf, size_t scmp, - unsigned char basis) + double scale, int cmp, int 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); else { for (size_t i = 1; i <= fe.basis(basis).size(); ++i) - for (size_t k = 1; k <= cmp; ++k) - EV(nf*(i-1)+k+scmp) += scale*fe.basis(basis)(i)*fe.detJxW; + for (size_t k = (cmp == 0 ? 1: cmp); + 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, - 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 void Convection(Matrix& EM, const FiniteElement& fe, 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) { - Advection(EM, fe, U, -scale, cmp, nf); - for (size_t i = 1;i <= fe.N.size();i++) - for (size_t j = 1;j <= fe.N.size();j++) { + Advection(EM, fe, U, -scale, basis); + for (size_t i = 1;i <= fe.basis(basis).size();i++) + for (size_t j = 1;j <= fe.basis(basis).size();j++) { for (size_t k = 1;k <= cmp;k++) { 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 { - Advection(EM, fe, U, scale, cmp, nf); - for (size_t i = 1;i <= fe.N.size();i++) - for (size_t j = 1;j <= fe.N.size();j++) { + Advection(EM, fe, U, scale, basis); + for (size_t i = 1;i <= fe.basis(basis).size();i++) + for (size_t j = 1;j <= fe.basis(basis).size();j++) { for (size_t k = 1;k <= cmp;k++) { 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. - //! \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. + //! \param[out] EM The element matrix to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \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, - size_t nf, double scale=1.0); + double scale=1.0, int basis=1, int tbasis=1); //! \brief Compute a divergence term. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] D Divergence of field. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] nf Number of fields in basis. - // + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] D Divergence of field + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Test function basis void Divergence(Vector& EV, const FiniteElement& fe, - const Vec3& D, double scale=1.0, - size_t cmp=1, size_t nf=1); + const Vec3& D, double scale=1.0, int basis=1); - //! \brief Compute a divergence/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 PressureDiv(Matrix& EM, const FiniteElement& fe, - 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. + //! \brief Compute a gradient term for a (potentially mixed) vector/scalar field. + //! \param[out] EM The element matrix to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis for field + //! \param[in] tbasis Test function basis 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. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] nf Number of fields in basis. + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + //! \param[in] tbasis Test function basis 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. - //! \param[out] EM The element matrix to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] stress Whether to add extra stress formulation terms. - //! \param[in] scmp Starting component. - //! \param[in] basis Basis to use. + //! \param[out] EM The element matrix to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + //! \param[in] stress Whether to add extra stress formulation terms + //! \param[in] basis Basis to use void Laplacian(Matrix& EM, const FiniteElement& fe, - double scale=1.0, size_t cmp=1, size_t nf=1, - bool stress=false, size_t scmp=0, unsigned char basis=1); + double scale=1.0, bool stress=false, int basis=1); //! \brief Compute a heteregenous coefficient laplacian. - //! \param[out] EM The element matrix to add contribution to. - //! \param[out] K The coefficient matrix. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Scaling factor for contribution. + //! \param[out] EM The element matrix to add contribution to + //! \param[out] K The coefficient matrix + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution 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. - //! \param[out] EM The element matrix to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] scmp Starting component. - //! \param[in] basis Basis to use. + //! \param[out] EM The element matrix to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use void Mass(Matrix& EM, const FiniteElement& fe, - double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, - unsigned char basis=1); + double scale=1.0, int basis=1); //! \brief Compute a source term. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] scmp Starting component. - //! \param[in] basis Basis to use. + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use + //! \param[in] cmp Component to add (0 for all) void Source(Vector& EV, const FiniteElement& fe, - double scale=1.0, size_t cmp=1, size_t nf=1, size_t scmp=0, - unsigned char basis=1); + double scale=1.0, int cmp=1, int basis=1); //! \brief Compute a vector-source term. - //! \param[out] EV The element vector to add contribution to. - //! \param[in] fe The finite element to evaluate for. - //! \param[in] scale Vector with contributions. - //! \param[in] scale Scaling factor for contribution. - //! \param[in] cmp Number of components to add. - //! \param[in] nf Number of fields in basis. - //! \param[in] scmp Starting component. + //! \param[out] EV The element vector to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Vector with contributions + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use 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