From c9bc665861e4ccf584562e05750effc11b4e2d95 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Oct 2016 15:32:44 +0200 Subject: [PATCH] changed: make WeakOperators a class with static functions allows for templating over operator classes --- Apps/Common/EqualOrderOperators.C | 238 ++++++++++++++++++ Apps/Common/EqualOrderOperators.h | 191 ++++++++++++++ Apps/Common/ResidualOperators.h | 114 --------- ...kOperators.C => TestEqualOrderOperators.C} | 52 ++-- Apps/Common/WeakOperators.C | 170 ------------- Apps/Common/WeakOperators.h | 145 ----------- 6 files changed, 455 insertions(+), 455 deletions(-) create mode 100644 Apps/Common/EqualOrderOperators.C create mode 100644 Apps/Common/EqualOrderOperators.h delete mode 100644 Apps/Common/ResidualOperators.h rename Apps/Common/Test/{TestWeakOperators.C => TestEqualOrderOperators.C} (78%) delete mode 100644 Apps/Common/WeakOperators.C delete mode 100644 Apps/Common/WeakOperators.h diff --git a/Apps/Common/EqualOrderOperators.C b/Apps/Common/EqualOrderOperators.C new file mode 100644 index 00000000..f0933f3b --- /dev/null +++ b/Apps/Common/EqualOrderOperators.C @@ -0,0 +1,238 @@ +//============================================================================== +//! +//! \file EqualOrderOperators.C +//! +//! \date Jul 22 2015 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Various discrete equal-ordered operators. +//! +//============================================================================== + +#include "EqualOrderOperators.h" +#include "FiniteElement.h" +#include "Vec3.h" + +//! \brief Helper for adding an element matrix to several components. +//! \param[out] EM The element matrix to add to. +//! \param[in] A The scalar element matrix to add. +//! \param[in] cmp Number of components to add matrix to +//! \param[in] nf Number of components in total matrix. +//! \param[in] scmp Index of first component to add matrix to. +static void addComponents(Matrix& EM, const Matrix& A, + size_t cmp, size_t nf, size_t scmp) +{ + if (cmp == 1 && nf == 1) + EM += A; + else + for (size_t i = 1; i <= A.rows(); ++i) + for (size_t j = 1; j <= A.cols(); ++j) + for (size_t k = 1; k <= cmp; ++k) + EM(nf*(i-1)+k+scmp,nf*(j-1)+k+scmp) += A(i, j); +} + + +//! \brief Helper applying a divergence (1) or a gradient (2) operation +template +static void DivGrad(Matrix& EM, const FiniteElement& fe, + double scale, int basis, int tbasis) +{ + 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 EqualOrderOperators::Weak::Advection(Matrix& EM, const FiniteElement& fe, + const Vec3& AC, + double scale, 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; + } + } + addComponents(EM, C, ncmp, ncmp, 0); +} + + +void EqualOrderOperators::Weak::Convection(Matrix& EM, const FiniteElement& fe, + const Vec3& U, const Tensor& dUdX, + double scale, bool conservative, int basis) +{ + size_t cmp = EM.rows() / fe.basis(basis).size(); + if (conservative) { + 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)*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, 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)*cmp+l,(i-1)*cmp+k) += scale*dUdX(l,k)*fe.basis(basis)(j)*fe.basis(basis)(i)*fe.detJxW; + } + } + } +} + + +void EqualOrderOperators::Weak::Divergence(Matrix& EM, const FiniteElement& fe, + double scale, int basis, int tbasis) +{ + DivGrad<1>(EM,fe,scale,basis,tbasis); +} + + +void EqualOrderOperators::Weak::Gradient(Matrix& EM, const FiniteElement& fe, + double scale, int basis, int tbasis) +{ + DivGrad<2>(EM,fe,scale,basis,tbasis); +} + + +void EqualOrderOperators::Weak::Divergence(Vector& EV, const FiniteElement& fe, + const Vec3& D, double scale, int basis) +{ + 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 += D[k-1]*fe.grad(basis)(i,k); + EV(i) += scale*div*fe.detJxW; + } +} + + +void EqualOrderOperators::Weak::Gradient(Vector& EV, const FiniteElement& fe, + double scale, int basis) +{ + 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 EqualOrderOperators::Weak::Laplacian(Matrix& EM, const FiniteElement& fe, + 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, cmp, 0); + if (stress) + 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(cmp*(j-1)+k,cmp*(i-1)+l) += scale*fe.grad(basis)(i,k)*fe.grad(basis)(j,l)*fe.detJxW; +} + + +void EqualOrderOperators::Weak::LaplacianCoeff(Matrix& EM, const Matrix& K, + const FiniteElement& fe, + double scale, int basis) +{ + Matrix KB; + KB.multiply(K,fe.grad(basis),false,true).multiply(scale*fe.detJxW); + EM.multiply(fe.grad(basis),KB,false,false,true); +} + + +void EqualOrderOperators::Weak::Mass(Matrix& EM, const FiniteElement& fe, + 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, ncmp, ncmp, 0); +} + + +void EqualOrderOperators::Weak::Source(Vector& EV, const FiniteElement& fe, + double scale, int cmp, int basis) +{ + 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 = (cmp == 0 ? 1: cmp); + k <= (cmp == 0 ? ncmp : cmp); ++k) + EV(ncmp*(i-1)+k) += scale*fe.basis(basis)(i)*fe.detJxW; + } +} + + +void EqualOrderOperators::Weak::Source(Vector& EV, const FiniteElement& fe, + const Vec3& f, double scale, int basis) +{ + 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) + EV(cmp*(i-1)+k) += scale*f[k-1]*fe.basis(basis)(i)*fe.detJxW; +} + + +void EqualOrderOperators::Residual::Convection(Vector& EV, const FiniteElement& fe, + const Vec3& U, const Tensor& dUdX, + const Vec3& UC, double scale, + bool conservative, int basis) +{ + size_t cmp = EV.size() / fe.basis(basis).size(); + if (conservative) { + 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)*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++) { + // Convection + double conv = 0.0; + for (size_t l = 1;l <= cmp;l++) + conv += UC[l-1]*dUdX(k,l); + conv *= scale*fe.detJxW; + + for (size_t i = 1;i <= fe.basis(basis).size();i++) + EV((i-1)*cmp + k) -= conv*fe.basis(basis)(i); + } + } +} + + +void EqualOrderOperators::Residual::Divergence(Vector& EV, const FiniteElement& fe, + const Tensor& dUdX, + double scale, size_t basis) +{ + 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(i) += scale*div*fe.basis(basis)(i)*fe.detJxW; + } +} diff --git a/Apps/Common/EqualOrderOperators.h b/Apps/Common/EqualOrderOperators.h new file mode 100644 index 00000000..111e224d --- /dev/null +++ b/Apps/Common/EqualOrderOperators.h @@ -0,0 +1,191 @@ +//============================================================================== +//! +//! \file EqualOrderOperators.h +//! +//! \date Jul 22 2015 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Various equal-ordered discrete operators. +//! +//============================================================================== + +#ifndef EQUAL_ORDER_OPERATORS_H +#define EQUAL_ORDER_OPERATORS_H + +class Vec3; + +#include "FiniteElement.h" +#include "MatVec.h" + +/*! \brief Common discrete operators using equal-ordered discretizations. + */ + +class EqualOrderOperators +{ +public: + //! \brief Common weak operators using equal-ordered discretizations. + class Weak { + public: + //! \brief Compute an advection term. + //! \param[out] EM The element matrix to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] AC Advecting field + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use + static void Advection(Matrix& EM, const FiniteElement& fe, + const Vec3& AC, double scale=1.0, int basis=1); + + //! \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(Matrix& EM, const FiniteElement& fe, + const Vec3& U, const Tensor& dUdX, double scale, + bool conservative=false, int basis=1); + + //! \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] scale Scaling factor for contribution + //! \param[in] basis Basis for field + //! \param[in] tbasis Test function basis + static void Divergence(Matrix& EM, const FiniteElement& fe, + 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] basis Test function basis + static void Divergence(Vector& EV, const FiniteElement& fe, + const Vec3& D, double scale=1.0, int basis=1); + + //! \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 + static void Gradient(Matrix& EM, const FiniteElement& fe, + 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] tbasis Test function basis + static void Gradient(Vector& EV, const FiniteElement& fe, + 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] stress Whether to add extra stress formulation terms + //! \param[in] basis Basis to use + static void Laplacian(Matrix& EM, const FiniteElement& fe, + 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 + static void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe, + 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] basis Basis to use + static void Mass(Matrix& EM, const FiniteElement& fe, + 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] basis Basis to use + //! \param[in] cmp Component to add (0 for all) + static void Source(Vector& EV, const FiniteElement& fe, + 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] basis Basis to use + static void Source(Vector& EV, const FiniteElement& fe, + const Vec3& f, double scale=1.0, int basis=1); + }; + + //! \brief Common weak residual operators using equal-ordered discretizations. + 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(Vector& EV, const FiniteElement& fe, + const Vec3& U, const Tensor& dUdX, const Vec3& UC, + double scale, bool conservative=false, int basis=1); + + //! \brief Compute a divergence 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] dUdX Gradient of field + //! \param[in] scale Scaling factor for contribution + //! \param[in] basis Basis to use + static void Divergence(Vector& EV, const FiniteElement& fe, + const Tensor& dUdX, double scale=1.0, + size_t basis=1); + + //! \brief Compute a laplacian 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] 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 + static void Laplacian(Vector& EV, const FiniteElement& fe, + const T& dUdX, double scale=1.0, + bool stress=false, int basis=1) + { + 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.grad(basis)(i,l); + diff *= scale*fe.detJxW; + + // Add negative residual to rhs of momentum equation + EV((i-1)*cmp + k) -= diff; + } + } + + // Use stress formulation + if (stress) { + 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)*cmp + k) -= scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW; + } + } + }; +}; + +#endif diff --git a/Apps/Common/ResidualOperators.h b/Apps/Common/ResidualOperators.h deleted file mode 100644 index 245e0299..00000000 --- a/Apps/Common/ResidualOperators.h +++ /dev/null @@ -1,114 +0,0 @@ -//============================================================================== -//! -//! \file ResidualOperators.h -//! -//! \date Aug 19 2015 -//! -//! \author Arne Morten Kvarving / SINTEF -//! -//! \brief Various weak, discrete residual operators -//! -//============================================================================== - -#ifndef RESIDUALOPERATORS_H_ -#define RESIDUALOPERATORS_H_ - -class Vec3; - -#include "FiniteElement.h" -#include "MatVec.h" - -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] 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, bool conservative=false, int basis=1) - { - size_t cmp = EV.size() / fe.basis(basis).size(); - if (conservative) { - 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)*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++) { - // Convection - double conv = 0.0; - for (size_t l = 1;l <= cmp;l++) - conv += UC[l-1]*dUdX(k,l); - conv *= scale*fe.detJxW; - - 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] basis Basis to use - template - void Divergence(Vector& EV, const FiniteElement& fe, - const T& dUdX, double scale=1.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(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] 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, - bool stress=false, int basis=1) - { - 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.grad(basis)(i,l); - diff *= scale*fe.detJxW; - - // Add negative residual to rhs of momentum equation - EV((i-1)*cmp + k) -= diff; - } - } - - // Use stress formulation - if (stress) { - 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)*cmp + k) -= scale*dUdX(l,k)*fe.grad(basis)(i,l)*fe.detJxW; - } - } -} - -#endif diff --git a/Apps/Common/Test/TestWeakOperators.C b/Apps/Common/Test/TestEqualOrderOperators.C similarity index 78% rename from Apps/Common/Test/TestWeakOperators.C rename to Apps/Common/Test/TestEqualOrderOperators.C index fac7e017..c92ef41c 100644 --- a/Apps/Common/Test/TestWeakOperators.C +++ b/Apps/Common/Test/TestEqualOrderOperators.C @@ -1,16 +1,16 @@ //============================================================================== //! -//! \file TestWeakoperators.C +//! \file TestEqualOrderOperators.C //! //! \date Feb 16 2015 //! //! \author Arne Morten Kvarving / SINTEF //! -//! \brief Tests for various mesh quality indicators +//! \brief Tests for various discrete equal-order operators. //! //============================================================================== -#include "WeakOperators.h" +#include "EqualOrderOperators.h" #include "gtest/gtest.h" #include "FiniteElement.h" @@ -36,7 +36,7 @@ static FiniteElement getFE() return fe; } -TEST(TestWeakOperators, Advection) +TEST(TestEqualOrderOperators, Advection) { FiniteElement fe = getFE(); @@ -44,14 +44,14 @@ TEST(TestWeakOperators, Advection) Vec3 U; U[0] = 1.0; U[1] = 2.0; // First component only - WeakOperators::Advection(EM_scalar, fe, U); + EqualOrderOperators::Weak::Advection(EM_scalar, fe, U); const DoubleVec EM_scalar_ref = {{ 7.0, 10}, {14.0, 20.0}}; check_matrix_equal(EM_scalar, EM_scalar_ref); Matrix EM_vec(2*2, 2*2); U[0] = 3.0; U[1] = 4.0; - WeakOperators::Advection(EM_vec, fe, U, 2.0); + EqualOrderOperators::Weak::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}, {60.0, 0.0, 88.0, 0.0}, @@ -60,12 +60,12 @@ TEST(TestWeakOperators, Advection) } -TEST(TestWeakOperators, Divergence) +TEST(TestEqualOrderOperators, Divergence) { FiniteElement fe = getFE(); Matrix EM_scalar(2, 2*2); - WeakOperators::Divergence(EM_scalar, fe); + EqualOrderOperators::Weak::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); @@ -73,7 +73,7 @@ TEST(TestWeakOperators, Divergence) Vector EV_scalar(4); Vec3 D; D[0] = 1.0; D[1] = 2.0; - WeakOperators::Divergence(EV_scalar, fe, D, 1.0); + EqualOrderOperators::Weak::Divergence(EV_scalar, fe, D, 1.0); ASSERT_NEAR(EV_scalar(1), 7.0, 1e-13); ASSERT_NEAR(EV_scalar(2), 10.0, 1e-13); ASSERT_NEAR(EV_scalar(3), 0.0, 1e-13); @@ -81,13 +81,13 @@ TEST(TestWeakOperators, Divergence) } -TEST(TestWeakOperators, Gradient) +TEST(TestEqualOrderOperators, Gradient) { FiniteElement fe = getFE(); Matrix EM_scalar(2*2,2); // single component + pressure - WeakOperators::Gradient(EM_scalar, fe); + EqualOrderOperators::Weak::Gradient(EM_scalar, fe); const DoubleVec EM_scalar_ref = {{-1.0,-2.0}, {-3.0,-6.0}, {-2.0,-4.0}, @@ -95,7 +95,7 @@ TEST(TestWeakOperators, Gradient) check_matrix_equal(EM_scalar, EM_scalar_ref); Vector EV_scalar(4); - WeakOperators::Gradient(EV_scalar, fe); + EqualOrderOperators::Weak::Gradient(EV_scalar, fe); 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); @@ -103,13 +103,13 @@ TEST(TestWeakOperators, Gradient) } -TEST(TestWeakOperators, Laplacian) +TEST(TestEqualOrderOperators, Laplacian) { FiniteElement fe = getFE(); // single scalar block Matrix EM_scalar(2,2); - WeakOperators::Laplacian(EM_scalar, fe); + EqualOrderOperators::Weak::Laplacian(EM_scalar, fe); const DoubleVec EM_scalar_ref = {{10.0, 14.0}, {14.0, 20.0}}; @@ -118,7 +118,7 @@ TEST(TestWeakOperators, Laplacian) // multiple (2) blocks in 3 component element matrix Matrix EM_multi(2*2,2*2); - WeakOperators::Laplacian(EM_multi, fe, 1.0); + EqualOrderOperators::Weak::Laplacian(EM_multi, fe, 1.0); const DoubleVec EM_multi_ref = {{10.0, 0.0, 14.0, 0.0}, { 0.0, 10.0, 0.0, 14.0}, @@ -129,7 +129,7 @@ TEST(TestWeakOperators, Laplacian) // stress formulation Matrix EM_stress(2*2,2*2); - WeakOperators::Laplacian(EM_stress, fe, 1.0, true); + EqualOrderOperators::Weak::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}, @@ -138,25 +138,25 @@ TEST(TestWeakOperators, Laplacian) check_matrix_equal(EM_stress, EM_stress_ref); Matrix EM_coeff(2, 2); - WeakOperators::LaplacianCoeff(EM_coeff, fe.dNdX, fe, 2.0); + EqualOrderOperators::Weak::LaplacianCoeff(EM_coeff, fe.dNdX, fe, 2.0); const DoubleVec EM_coeff_ref = {{104.0, 148.0}, {152.0, 216.0}}; check_matrix_equal(EM_coeff, EM_coeff_ref); } -TEST(TestWeakOperators, Mass) +TEST(TestEqualOrderOperators, Mass) { FiniteElement fe = getFE(); Matrix EM_scalar(2,2); - WeakOperators::Mass(EM_scalar, fe, 2.0); + EqualOrderOperators::Weak::Mass(EM_scalar, fe, 2.0); const DoubleVec EM_scalar_ref = {{2.0, 4.0}, {4.0, 8.0}}; check_matrix_equal(EM_scalar, EM_scalar_ref); Matrix EM_vec(2*2,2*2); - WeakOperators::Mass(EM_vec, fe); + EqualOrderOperators::Weak::Mass(EM_vec, fe); const DoubleVec EM_vec_ref = {{1.0, 0.0, 2.0, 0.0}, {0.0, 1.0, 0.0, 2.0}, @@ -166,37 +166,37 @@ TEST(TestWeakOperators, Mass) } -TEST(TestWeakOperators, Source) +TEST(TestEqualOrderOperators, Source) { FiniteElement fe = getFE(); Vector EV_scalar(2); - WeakOperators::Source(EV_scalar, fe, 2.0); + EqualOrderOperators::Weak::Source(EV_scalar, fe, 2.0); ASSERT_NEAR(EV_scalar(1), 2.0, 1e-13); ASSERT_NEAR(EV_scalar(2), 4.0, 1e-13); Vector EV_vec(4); Vec3 f; f[0] = 1.0; f[1] = 2.0; - WeakOperators::Source(EV_vec, fe, f, 1.0); + EqualOrderOperators::Weak::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), 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); + EqualOrderOperators::Weak::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); + EqualOrderOperators::Weak::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); + EqualOrderOperators::Weak::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); diff --git a/Apps/Common/WeakOperators.C b/Apps/Common/WeakOperators.C deleted file mode 100644 index 2cf505e9..00000000 --- a/Apps/Common/WeakOperators.C +++ /dev/null @@ -1,170 +0,0 @@ -//============================================================================== -//! -//! \file WeakOperators.C -//! -//! \date Jul 22 2015 -//! -//! \author Arne Morten Kvarving / SINTEF -//! -//! \brief Various weak, discrete operators. -//! -//============================================================================== - -#include "WeakOperators.h" -#include "FiniteElement.h" -#include "Vec3.h" - -namespace WeakOperators { - //! \brief Helper for adding an element matrix to several components. - //! \param[out] EM The element matrix to add to. - //! \param[in] A The scalar element matrix to add. - //! \param[in] cmp Number of components to add matrix to - //! \param[in] nf Number of components in total matrix. - //! \param[in] scmp Index of first component to add matrix to. - static void addComponents(Matrix& EM, const Matrix& A, - size_t cmp, size_t nf, size_t scmp) - { - if (cmp == 1 && nf == 1) - EM += A; - else - for (size_t i = 1; i <= A.rows(); ++i) - for (size_t j = 1; j <= A.cols(); ++j) - for (size_t k = 1; k <= cmp; ++k) - EM(nf*(i-1)+k+scmp,nf*(j-1)+k+scmp) += A(i, j); - } - - - //! \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, int basis, int tbasis) - { - 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, 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.dNdX(j,k); - - C(i,j) *= scale*fe.N(i)*fe.detJxW; - } - } - addComponents(EM, C, ncmp, ncmp, 0); - } - - - void Divergence(Matrix& EM, const FiniteElement& fe, double scale, - int basis, int tbasis) - { - DivGrad<1>(EM,fe,scale,basis,tbasis); - } - - - void Gradient(Matrix& EM, const FiniteElement& fe, double scale, - int basis, int tbasis) - { - DivGrad<2>(EM,fe,scale,basis,tbasis); - } - - - void Divergence(Vector& EV, const FiniteElement& fe, - const Vec3& D, double scale, int basis) - { - 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 += D[k-1]*fe.grad(basis)(i,k); - EV(i) += scale*div*fe.detJxW; - } - } - - - void Gradient(Vector& EV, const FiniteElement& fe, - double scale, int basis) - { - 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, 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, cmp, 0); - if (stress) - 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(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, int basis) - { - Matrix KB; - 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, 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, ncmp, ncmp, 0); - } - - - void Source(Vector& EV, const FiniteElement& fe, - double scale, int cmp, int basis) - { - 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 = (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, int basis) - { - 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) - EV(cmp*(i-1)+k) += scale*f[k-1]*fe.basis(basis)(i)*fe.detJxW; - } -} diff --git a/Apps/Common/WeakOperators.h b/Apps/Common/WeakOperators.h deleted file mode 100644 index 96fb7f0e..00000000 --- a/Apps/Common/WeakOperators.h +++ /dev/null @@ -1,145 +0,0 @@ -//============================================================================== -//! -//! \file WeakOperators.h -//! -//! \date Jul 22 2015 -//! -//! \author Arne Morten Kvarving / SINTEF -//! -//! \brief Various weak, discrete operators. -//! -//============================================================================== - -#ifndef WEAKOPERATORS_H_ -#define WEAKOPERATORS_H_ - -class Vec3; - -#include "FiniteElement.h" -#include "MatVec.h" - -namespace WeakOperators -{ - //! \brief Compute an advection term. - //! \param[out] EM The element matrix to add contribution to - //! \param[in] fe The finite element to evaluate for - //! \param[in] AC Advecting field - //! \param[in] scale Scaling factor for contribution - //! \param[in] basis Basis to use - void Advection(Matrix& EM, const FiniteElement& fe, - const Vec3& AC, double scale=1.0, int basis=1); - - //! \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 - template - void Convection(Matrix& EM, const FiniteElement& fe, - const Vec3& U, const T& dUdX, double scale, - bool conservative=false, int basis=1) - { - size_t cmp = EM.rows() / fe.basis(basis).size(); - if (conservative) { - 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)*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, 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)*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] scale Scaling factor for contribution - //! \param[in] basis Basis for field - //! \param[in] tbasis Test function basis - void Divergence(Matrix& EM, const FiniteElement& fe, - 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] basis Test function basis - void Divergence(Vector& EV, const FiniteElement& fe, - const Vec3& D, double scale=1.0, int basis=1); - - //! \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, - 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] tbasis Test function basis - void Gradient(Vector& EV, const FiniteElement& fe, - 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] 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, 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 - void LaplacianCoeff(Matrix& EM, const Matrix& K, const FiniteElement& fe, - 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] basis Basis to use - void Mass(Matrix& EM, const FiniteElement& fe, - 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] basis Basis to use - //! \param[in] cmp Component to add (0 for all) - void Source(Vector& EV, const FiniteElement& fe, - 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] basis Basis to use - void Source(Vector& EV, const FiniteElement& fe, - const Vec3& f, double scale=1.0, int basis=1); -} - -#endif