Add economical Gauss quadrature for square and cube (up to degree 7)

src/analysis/CMakeLists.txt

@@ -2,6 +2,7 @@
 
 add_library(
   PugsAnalysis
+  EconomicalGaussQuadrature.cpp
   SimplicialGaussQuadrature.cpp
   TensorialGaussLegendreQuadrature.cpp
   TensorialGaussLobattoQuadrature.cpp)

src/analysis/EconomicalGaussQuadrature.cpp

+#include <analysis/EconomicalGaussQuadrature.hpp>
+#include <analysis/TensorialGaussLegendreQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+template <>
+void
+EconomicalGaussQuadrature<1>::_buildPointAndWeightLists()
+{
+  TensorialGaussLegendreQuadrature<1> gauss(m_degree);
+
+  m_point_list  = gauss.pointList();
+  m_weight_list = gauss.weightList();
+}
+
+template <>
+void
+EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
+{
+  switch (m_degree) {
+  case 0:
+  case 1: {
+    constexpr size_t nb_points = 1;
+    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0]  = {0, 0};
+    weight_list[0] = 4;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 2:
+  case 3: {
+    constexpr size_t nb_points = 4;
+    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0] = {+0.577350269189626, +0.577350269189626};
+    point_list[1] = {+0.577350269189626, -0.577350269189626};
+    point_list[2] = {-0.577350269189626, +0.577350269189626};
+    point_list[3] = {-0.577350269189626, -0.577350269189626};
+
+    weight_list[0] = 1;
+    weight_list[1] = 1;
+    weight_list[2] = 1;
+    weight_list[3] = 1;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 4:
+  case 5: {
+    constexpr size_t nb_points = 8;
+    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0] = {+0.683130051063973, +0.000000000000000};
+    point_list[1] = {-0.683130051063973, +0.000000000000000};
+    point_list[2] = {+0.000000000000000, +0.683130051063973};
+    point_list[3] = {+0.000000000000000, -0.683130051063973};
+    point_list[4] = {+0.881917103688197, +0.881917103688197};
+    point_list[5] = {+0.881917103688197, -0.881917103688197};
+    point_list[6] = {-0.881917103688197, +0.881917103688197};
+    point_list[7] = {-0.881917103688197, -0.881917103688197};
+
+    weight_list[0] = 0.816326530612245;
+    weight_list[1] = 0.816326530612245;
+    weight_list[2] = 0.816326530612245;
+    weight_list[3] = 0.816326530612245;
+    weight_list[4] = 0.183673469387755;
+    weight_list[5] = 0.183673469387755;
+    weight_list[6] = 0.183673469387755;
+    weight_list[7] = 0.183673469387755;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 6:
+  case 7: {
+    constexpr size_t nb_points = 12;
+    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0]  = {+0.925820099772551, +0.000000000000000};
+    point_list[1]  = {-0.925820099772551, +0.000000000000000};
+    point_list[2]  = {+0.000000000000000, +0.925820099772551};
+    point_list[3]  = {+0.000000000000000, -0.925820099772551};
+    point_list[4]  = {+0.805979782918599, +0.805979782918599};
+    point_list[5]  = {+0.805979782918599, -0.805979782918599};
+    point_list[6]  = {-0.805979782918599, +0.805979782918599};
+    point_list[7]  = {-0.805979782918599, -0.805979782918599};
+    point_list[8]  = {+0.380554433208316, +0.380554433208316};
+    point_list[9]  = {+0.380554433208316, -0.380554433208316};
+    point_list[10] = {-0.380554433208316, +0.380554433208316};
+    point_list[11] = {-0.380554433208316, -0.380554433208316};
+
+    weight_list[0]  = 0.241975308641975;
+    weight_list[1]  = 0.241975308641975;
+    weight_list[2]  = 0.241975308641975;
+    weight_list[3]  = 0.241975308641975;
+    weight_list[4]  = 0.237431774690630;
+    weight_list[5]  = 0.237431774690630;
+    weight_list[6]  = 0.237431774690630;
+    weight_list[7]  = 0.237431774690630;
+    weight_list[8]  = 0.520592916667394;
+    weight_list[9]  = 0.520592916667394;
+    weight_list[10] = 0.520592916667394;
+    weight_list[11] = 0.520592916667394;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  default: {
+    throw NormalError("Gauss quadrature formulae handle degrees up to 7 on quadrangles");
+  }
+  }
+}
+
+template <>
+void
+EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
+{
+  switch (m_degree) {
+  case 0:
+  case 1: {
+    constexpr size_t nb_points = 1;
+    SmallArray<TinyVector<3>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0] = {0, 0, 0};
+
+    weight_list[0] = 8;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 2:
+  case 3: {
+    constexpr size_t nb_points = 6;
+    SmallArray<TinyVector<3>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0] = {+1.0000000000, +0.0000000000, +0.0000000000};
+    point_list[1] = {-1.0000000000, +0.0000000000, +0.0000000000};
+    point_list[2] = {+0.0000000000, +1.0000000000, +0.0000000000};
+    point_list[3] = {+0.0000000000, -1.0000000000, +0.0000000000};
+    point_list[4] = {+0.0000000000, +0.0000000000, +1.0000000000};
+    point_list[5] = {+0.0000000000, +0.0000000000, -1.0000000000};
+
+    weight_list[0] = 1.3333333333;
+    weight_list[1] = 1.3333333333;
+    weight_list[2] = 1.3333333333;
+    weight_list[3] = 1.3333333333;
+    weight_list[4] = 1.3333333333;
+    weight_list[5] = 1.3333333333;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 4:
+  case 5: {
+    constexpr size_t nb_points = 14;
+    SmallArray<TinyVector<3>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0]  = {+0.7958224257, +0.0000000000, +0.0000000000};
+    point_list[1]  = {-0.7958224257, +0.0000000000, +0.0000000000};
+    point_list[2]  = {+0.0000000000, +0.7958224257, +0.0000000000};
+    point_list[3]  = {+0.0000000000, -0.7958224257, +0.0000000000};
+    point_list[4]  = {+0.0000000000, +0.0000000000, +0.7958224257};
+    point_list[5]  = {+0.0000000000, +0.0000000000, -0.7958224257};
+    point_list[6]  = {+0.7587869106, +0.7587869106, +0.7587869106};
+    point_list[7]  = {+0.7587869106, -0.7587869106, +0.7587869106};
+    point_list[8]  = {+0.7587869106, +0.7587869106, -0.7587869106};
+    point_list[9]  = {+0.7587869106, -0.7587869106, -0.7587869106};
+    point_list[10] = {-0.7587869106, +0.7587869106, +0.7587869106};
+    point_list[11] = {-0.7587869106, -0.7587869106, +0.7587869106};
+    point_list[12] = {-0.7587869106, +0.7587869106, -0.7587869106};
+    point_list[13] = {-0.7587869106, -0.7587869106, -0.7587869106};
+
+    weight_list[0]  = 0.8864265927;
+    weight_list[1]  = 0.8864265927;
+    weight_list[2]  = 0.8864265927;
+    weight_list[3]  = 0.8864265927;
+    weight_list[4]  = 0.8864265927;
+    weight_list[5]  = 0.8864265927;
+    weight_list[6]  = 0.3351800554;
+    weight_list[7]  = 0.3351800554;
+    weight_list[8]  = 0.3351800554;
+    weight_list[9]  = 0.3351800554;
+    weight_list[10] = 0.3351800554;
+    weight_list[11] = 0.3351800554;
+    weight_list[12] = 0.3351800554;
+    weight_list[13] = 0.3351800554;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 6:
+  case 7: {
+    constexpr size_t nb_points = 27;
+    SmallArray<TinyVector<3>> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    point_list[0]  = {+0.0000000000, +0.0000000000, +0.0000000000};
+    point_list[1]  = {+0.8484180014, +0.0000000000, +0.0000000000};
+    point_list[2]  = {-0.8484180014, +0.0000000000, +0.0000000000};
+    point_list[3]  = {+0.0000000000, +0.8484180014, +0.0000000000};
+    point_list[4]  = {+0.0000000000, -0.8484180014, +0.0000000000};
+    point_list[5]  = {+0.0000000000, +0.0000000000, +0.8484180014};
+    point_list[6]  = {+0.0000000000, +0.0000000000, -0.8484180014};
+    point_list[7]  = {+0.6528164721, +0.6528164721, +0.6528164721};
+    point_list[8]  = {+0.6528164721, -0.6528164721, +0.6528164721};
+    point_list[9]  = {+0.6528164721, +0.6528164721, -0.6528164721};
+    point_list[10] = {+0.6528164721, -0.6528164721, -0.6528164721};
+    point_list[11] = {-0.6528164721, +0.6528164721, +0.6528164721};
+    point_list[12] = {-0.6528164721, -0.6528164721, +0.6528164721};
+    point_list[13] = {-0.6528164721, +0.6528164721, -0.6528164721};
+    point_list[14] = {-0.6528164721, -0.6528164721, -0.6528164721};
+    point_list[15] = {+0.0000000000, +1.1064128986, +1.1064128986};
+    point_list[16] = {+0.0000000000, -1.1064128986, +1.1064128986};
+    point_list[17] = {+0.0000000000, +1.1064128986, -1.1064128986};
+    point_list[18] = {+0.0000000000, -1.1064128986, -1.1064128986};
+    point_list[19] = {+1.1064128986, +0.0000000000, +1.1064128986};
+    point_list[20] = {-1.1064128986, +0.0000000000, +1.1064128986};
+    point_list[21] = {+1.1064128986, +0.0000000000, -1.1064128986};
+    point_list[22] = {-1.1064128986, +0.0000000000, -1.1064128986};
+    point_list[23] = {+1.1064128986, +1.1064128986, +0.0000000000};
+    point_list[24] = {-1.1064128986, +1.1064128986, +0.0000000000};
+    point_list[25] = {+1.1064128986, -1.1064128986, +0.0000000000};
+    point_list[26] = {-1.1064128986, -1.1064128986, +0.0000000000};
+
+    weight_list[0]  = 0.7880734827;
+    weight_list[1]  = 0.4993690023;
+    weight_list[2]  = 0.4993690023;
+    weight_list[3]  = 0.4993690023;
+    weight_list[4]  = 0.4993690023;
+    weight_list[5]  = 0.4993690023;
+    weight_list[6]  = 0.4993690023;
+    weight_list[7]  = 0.4785084494;
+    weight_list[8]  = 0.4785084494;
+    weight_list[9]  = 0.4785084494;
+    weight_list[10] = 0.4785084494;
+    weight_list[11] = 0.4785084494;
+    weight_list[12] = 0.4785084494;
+    weight_list[13] = 0.4785084494;
+    weight_list[14] = 0.4785084494;
+    weight_list[15] = 0.0323037423;
+    weight_list[16] = 0.0323037423;
+    weight_list[17] = 0.0323037423;
+    weight_list[18] = 0.0323037423;
+    weight_list[19] = 0.0323037423;
+    weight_list[20] = 0.0323037423;
+    weight_list[21] = 0.0323037423;
+    weight_list[22] = 0.0323037423;
+    weight_list[23] = 0.0323037423;
+    weight_list[24] = 0.0323037423;
+    weight_list[25] = 0.0323037423;
+    weight_list[26] = 0.0323037423;
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  default: {
+    throw NormalError("Gauss quadrature formulae handle degrees up to 7 on hexahedra");
+  }
+  }
+}

src/analysis/EconomicalGaussQuadrature.hpp

+#ifndef ECONOMICAL_GAUSS_QUADRATURE_HPP
+#define ECONOMICAL_GAUSS_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+/**
+ * Defines Economical Gauss quadrature on the reference element
+ * \f$]-1,1[^d\f$, where \f$d\f$ denotes the \a Dimension.
+ *
+ * \note formulae are extracted from High-order Finite Element Method [2004 -  Chapman & Hall]
+ */
+template <size_t Dimension>
+class EconomicalGaussQuadrature final : public QuadratureForumla<Dimension>
+{
+ private:
+  void _buildPointAndWeightLists();
+
+ public:
+  EconomicalGaussQuadrature(EconomicalGaussQuadrature&&)      = default;
+  EconomicalGaussQuadrature(const EconomicalGaussQuadrature&) = default;
+
+  explicit EconomicalGaussQuadrature(const size_t degree) : QuadratureForumla<Dimension>(degree)
+  {
+    this->_buildPointAndWeightLists();
+  }
+
+  EconomicalGaussQuadrature()  = delete;
+  ~EconomicalGaussQuadrature() = default;
+};
+
+#endif   // ECONOMICAL_GAUSS_QUADRATURE_HPP

tests/CMakeLists.txt

@@ -68,6 +68,7 @@ add_executable (unit_tests
   test_DiscreteFunctionType.cpp
   test_DiscreteFunctionUtils.cpp
   test_DoWhileProcessor.cpp
+  test_EconomicalGaussQuadrature.cpp
   test_EigenvalueSolver.cpp
   test_EmbeddedData.cpp
   test_EmbeddedIDiscreteFunctionUtils.cpp