diff --git a/Apps/Common/CompatibleOperators.C b/Apps/Common/CompatibleOperators.C new file mode 100644 index 00000000..df94934c --- /dev/null +++ b/Apps/Common/CompatibleOperators.C @@ -0,0 +1,108 @@ +//============================================================================== +//! +//! \file CompatibleOperators.C +//! +//! \date Oct 9 2016 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Various weak, discrete div-compatible operators. +//! +//============================================================================== + +#include "CompatibleOperators.h" +#include "EqualOrderOperators.h" +#include "FiniteElement.h" +#include "Vec3.h" + + +void CompatibleOperators::Weak::Advection(std::vector& EM, + const FiniteElement& fe, + const Vec3& AC, double scale) +{ + size_t nsd = fe.grad(1).cols(); + for (size_t n = 1; n <= nsd; ++n) + for (size_t i = 1; i <= fe.basis(n).size(); ++i) + for (size_t j = 1; j <= fe.basis(n).size(); ++j) + for (size_t k = 1; k <= nsd; ++k) + EM[n](i,j) += scale*AC[k]*fe.grad(n)(j,k)*fe.basis(n)(i)*fe.detJxW; +} + +void CompatibleOperators::Weak::Gradient(std::vector& EM, + const FiniteElement& fe, + double scale) +{ + size_t nsd = fe.grad(1).cols(); + for (size_t n = 1; n <= nsd; ++n) + for (size_t i=1; i <= fe.basis(n).size(); ++i) + for (size_t j=1; j <= fe.basis(nsd+1).size(); ++j) + EM[8+4*(n-1)](i,j) += -scale*fe.grad(n)(i,n)*fe.basis(nsd+1)(j)*fe.detJxW; +} + + +void CompatibleOperators::Weak::Gradient(Vectors& EV, const FiniteElement& fe, + double scale) +{ + for (size_t n = 1; n <= fe.grad(n).cols(); ++n) + for (size_t i = 1; i <= fe.basis(n).size(); ++i) + EV[n](i) += scale*fe.grad(n)(i,n)*fe.detJxW; +} + + +void CompatibleOperators::Weak::Laplacian(std::vector& EM, + const FiniteElement& fe, + double scale, bool stress) +{ + size_t nsd = fe.grad(1).cols(); + for (size_t n = 1; n <= nsd; ++n) + EqualOrderOperators::Weak::Laplacian(EM[n], fe, scale, false, n); + + for (size_t m = 1; m <= nsd && stress; m++) + for (size_t i = 1; i <= fe.basis(m).size(); ++i) + for (size_t n = m; n <= nsd; n++) + for (size_t j = 1; j <= fe.basis(n).size(); ++j) { + int idx = m == n ? m : (m == 1 ? 5+n-m : 10+n-m); + EM[idx](i,j) += scale* fe.grad(n)(j,m) * fe.grad(m)(i,n) * fe.detJxW; + } +} + + +void CompatibleOperators::Weak::Mass(std::vector& EM, + const FiniteElement& fe, double scale) +{ + for (size_t k = 1; k <= fe.grad(1).cols(); ++k) + EqualOrderOperators::Weak::Mass(EM[k], fe, scale, k); +} + +void CompatibleOperators::Weak::Source(Vectors& EV, + const FiniteElement& fe, + const Vec3& f, double scale) +{ + for (size_t k = 1; k <= fe.grad(1).cols(); ++k) + EqualOrderOperators::Weak::Source(EV[k], fe, scale*f[k-1], 1, k); +} + + +void CompatibleOperators::Weak::Source(Vectors& EV, + const FiniteElement& fe, + double scale) +{ + for (size_t k = 1; k <= fe.grad(1).cols(); ++k) + EqualOrderOperators::Weak::Source(EV[k], fe, scale, 1, k); +} + + +void CompatibleOperators::Residual::Laplacian(Vectors& EV, + const FiniteElement& fe, + const Tensor& dUdX, + double scale, bool stress) +{ + size_t nsd = fe.grad(1).cols(); + for (size_t k = 1; k <= nsd; ++k) + for (size_t i = 1; i <= fe.basis(k).size(); ++i) { + double diff = 0.0; + for (size_t m = 1; m <= nsd; ++m) + diff += fe.grad(k)(i,m)*dUdX(k,m); + EV[k](i) += scale*diff*fe.detJxW; + } +} diff --git a/Apps/Common/CompatibleOperators.h b/Apps/Common/CompatibleOperators.h new file mode 100644 index 00000000..6ca0926a --- /dev/null +++ b/Apps/Common/CompatibleOperators.h @@ -0,0 +1,102 @@ +//============================================================================== +//! +//! \file CompatibleOperators.h +//! +//! \date Oct 9 2016 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Various weak, div-compatible discrete operators. +//! +//============================================================================== + +#ifndef COMPATIBLE_OPERATORS_H_ +#define COMPATIBLE_OPERATORS_H_ + +class FiniteElement; +class Tensor; +class Vec3; + +#include "MatVec.h" + +/*! \brief Common operators using div-compatible discretizations. + * \details The operators use the block ordering used in the BlockElmMats class. + */ + +class CompatibleOperators +{ +public: + //! \brief Common weak operators using div-compatible discretizations. + class Weak { + public: + //! \brief Compute an advection term. + //! \param[out] EM The element matrices 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 + static void Advection(std::vector& EM, const FiniteElement& fe, + const Vec3& AC, 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] scale Scaling factor for contribution + static void Gradient(std::vector& EM, + const FiniteElement& fe, double scale=1.0); + + //! \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 + static void Gradient(Vectors& EV, const FiniteElement& fe, + double scale=1.0); + + //! \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 + //! \details Only the upper blocks are added with the stress formulation + static void Laplacian(std::vector& EM, const FiniteElement& fe, + double scale=1.0, bool stress=false); + + //! \brief Compute a mass term. + //! \param[out] EM The element matrices to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] scale Scaling factor for contribution + static void Mass(std::vector& EM, + const FiniteElement& fe, double scale=1.0); + + //! \brief Compute a vector-source term. + //! \param[out] EV The element vectors to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] f Vector with contributions + //! \param[in] scale Scaling factor for contribution + static void Source(Vectors& EV, const FiniteElement& fe, double scale=1.0); + + //! \brief Compute a vector-source term. + //! \param[out] EV The element vectors to add contribution to + //! \param[in] fe The finite element to evaluate for + //! \param[in] f Vector with contributions + //! \param[in] scale Scaling factor for contribution + static void Source(Vectors& EV, const FiniteElement& fe, + const Vec3& f, double scale=1.0); + }; + + //! \brief Common weak residual operators using div-compatible discretizations. + class Residual + { + public: + //! \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 + static void Laplacian(Vectors& EV, const FiniteElement& fe, + const Tensor& dUdX, double scale=1.0, + bool stress=false); + }; +}; + +#endif diff --git a/Apps/Common/Test/TestCompatibleOperators.C b/Apps/Common/Test/TestCompatibleOperators.C new file mode 100644 index 00000000..ec0aa118 --- /dev/null +++ b/Apps/Common/Test/TestCompatibleOperators.C @@ -0,0 +1,239 @@ +//============================================================================== +//! +//! \file TestCompatibleOperators.C +//! +//! \date Oct 17 2016 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Tests for various discrete div-compatible operators. +//! +//============================================================================== + +#include "CompatibleOperators.h" +#include "FiniteElement.h" +#include "gtest/gtest.h" + +typedef std::vector> DoubleVec; +const auto&& check_matrix_equal = [](const Matrix& A, const DoubleVec& B) + { + for (size_t i=1;i<=A.rows();++i) + for (size_t j=1;j<=A.cols();++j) + ASSERT_NEAR(A(i,j), B[i-1][j-1], 1e-13); + }; + +static MxFiniteElement getFE() +{ + MxFiniteElement fe({6,6,4}); + for (size_t i = 1; i <= 2; ++i) { + fe.grad(i).resize(6,2); + fe.basis(i).resize(6); + for (size_t j = 1; j <= 6; ++j) { + fe.grad(i)(j,1) = 1.0+12*(i-1)+2*(j-1); + fe.grad(i)(j,2) = 2.0+12*(i-1)+2*(j-1); + fe.basis(i)(j) = j + 6*(i-1); + } + } + + fe.grad(3).resize(4,2); + fe.basis(3).resize(4); + for (size_t j = 1; j <= 4; ++j) { + fe.grad(3)(j,1) = 1.0 + 14.0 + 2*(j-1); + fe.grad(3)(j,2) = 2.0 + 14.0 + 2*(j-1); + fe.basis(3)(j) = 12 + j; + } + + return fe; +} + + +/* +TESTI(TestEqualOrderOperators, Advection) +{ + FiniteElement fe = getFE(); + + Matrix EM_scalar(2, 2); + Vec3 U; + U[0] = 1.0; U[1] = 2.0; + // First component only + 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; + 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}, + { 0.0, 60.0, 0.0, 88.0}}; + check_matrix_equal(EM_vec, EM_vec_ref); +} + + +TESTI(TestEqualOrderOperators, Divergence) +{ + FiniteElement fe = getFE(); + + Matrix EM_scalar(2, 2*2); + 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); + + Vector EV_scalar(4); + Vec3 D; + D[0] = 1.0; D[1] = 2.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); + ASSERT_NEAR(EV_scalar(4), 0.0, 1e-13); +} + + +TESTI(TestEqualOrderOperators, Gradient) +{ + FiniteElement fe = getFE(); + + Matrix EM_scalar(2*2,2); + // single component + pressure + EqualOrderOperators::Weak::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); + 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); + ASSERT_NEAR(EV_scalar(4), fe.dNdX(2,2), 1e-13); +} +*/ + + +TEST(TestCompatibleOperators, Laplacian) +{ + MxFiniteElement fe = getFE(); + + std::vector EM(3); + EM[1].resize(6,6); + EM[2].resize(6,6); + + CompatibleOperators::Weak::Laplacian(EM, fe); + + const DoubleVec EM_1_ref = {{ 5.0, 11.0, 17.0, 23.0, 29.0, 35.0}, + {11.0, 25.0, 39.0, 53.0, 67.0, 81.0}, + {17.0, 39.0, 61.0, 83.0, 105.0, 127.0}, + {23.0, 53.0, 83.0, 113.0, 143.0, 173.0}, + {29.0, 67.0, 105.0, 143.0, 181.0, 219.0}, + {35.0, 81.0, 127.0, 173.0, 219.0, 265.0}}; + + check_matrix_equal(EM[1], EM_1_ref); + + const DoubleVec EM_2_ref = {{365.0, 419.0, 473.0, 527.0, 581.0, 635.0}, + {419.0, 481.0, 543.0, 605.0, 667.0, 729.0}, + {473.0, 543.0, 613.0, 683.0, 753.0, 823.0}, + {527.0, 605.0, 683.0, 761.0, 839.0, 917.0}, + {581.0, 667.0, 753.0, 839.0, 925.0, 1011.0}, + {635.0, 729.0, 823.0, 917.0, 1011.0, 1105.0}}; + + check_matrix_equal(EM[2], EM_2_ref); + + // stress formulation + std::vector EM_stress(21); + + EM_stress[1].resize(6,6); + EM_stress[2].resize(6,6); + EM_stress[6].resize(6,6); + EM_stress[10].resize(6,6); + + CompatibleOperators::Weak::Laplacian(EM_stress, fe, 1.0, true); + std::cout << EM_stress[1] << std::endl; + +/* + 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}, + { 6.0, 26.0, 8.0, 36.0}}; + + check_matrix_equal(EM_stress, EM_stress_ref); + + Matrix EM_coeff(2, 2); + 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(TestCompatibleOperators, Mass) +{ + MxFiniteElement fe = getFE(); + + std::vector EM_vec(3); + EM_vec[1].resize(6,6); + EM_vec[2].resize(6,6); + + CompatibleOperators::Weak::Mass(EM_vec, fe); + + const DoubleVec EM_1_ref = {{1.0, 2.0, 3.0, 4.0, 5.0, 6.0}, + {2.0, 4.0, 6.0, 8.0, 10.0, 12.0}, + {3.0, 6.0, 9.0, 12.0, 15.0, 18.0}, + {4.0, 8.0, 12.0, 16.0, 20.0, 24.0}, + {5.0, 10.0, 15.0, 20.0, 25.0, 30.0}, + {6.0, 12.0, 18.0, 24.0, 30.0, 36.0}}; + + check_matrix_equal(EM_vec[1], EM_1_ref); + + const DoubleVec EM_2_ref = {{49.0, 56.0, 63.0, 70.0, 77.0, 84.0}, + {56.0, 64.0, 72.0, 80.0, 88.0, 96.0}, + {63.0, 72.0, 81.0, 90.0, 99.0, 108.0}, + {70.0, 80.0, 90.0, 100.0, 110.0, 120.0}, + {77.0, 88.0, 99.0, 110.0, 121.0, 132.0}, + {84.0, 96.0, 108.0, 120.0, 132.0, 144.0}}; + + check_matrix_equal(EM_vec[2], EM_2_ref); +} + + +TEST(TestCompatibleOperators, Source) +{ + MxFiniteElement fe = getFE(); + + Vectors EV_scalar(3); + EV_scalar[1].resize(6); + EV_scalar[2].resize(6); + + CompatibleOperators::Weak::Source(EV_scalar, fe, 2.0); + + for (size_t i = 1; i <= 6; ++i) { + ASSERT_FLOAT_EQ(EV_scalar[1](i), 2.0*fe.basis(1)(i)); + ASSERT_FLOAT_EQ(EV_scalar[2](i), 2.0*fe.basis(2)(i)); + } + + EV_scalar[1].fill(0.0); + EV_scalar[2].fill(0.0); + Vec3 f; + f[0] = 1.0; + f[1] = 2.0; + + CompatibleOperators::Weak::Source(EV_scalar, fe, f, 1.0); + for (size_t i = 1; i <= 6; ++i) { + ASSERT_FLOAT_EQ(EV_scalar[1](i), fe.basis(1)(i)); + ASSERT_FLOAT_EQ(EV_scalar[2](i), 2.0*fe.basis(2)(i)); + } + + EV_scalar[1].fill(0.0); + EV_scalar[2].fill(0.0); + CompatibleOperators::Weak::Source(EV_scalar, fe, f, 2.0); + for (size_t i = 1; i <= 6; ++i) { + ASSERT_FLOAT_EQ(EV_scalar[1](i), 2.0*fe.basis(1)(i)); + ASSERT_FLOAT_EQ(EV_scalar[2](i), 4.0*fe.basis(2)(i)); + } +}