diff --git a/src/analysis/EconomicalGaussQuadrature.cpp b/src/analysis/EconomicalGaussQuadrature.cpp
index 04560888cdd38a3a5ce8299fa65185d84c2ceebc..efc83e10070ea381662c8a69fc46eacebb946cba 100644
--- a/src/analysis/EconomicalGaussQuadrature.cpp
+++ b/src/analysis/EconomicalGaussQuadrature.cpp
@@ -24,7 +24,7 @@ EconomicalGaussQuadrature<2>::_buildPointAndWeightLists()
 
     size_t k = 0;
     for (size_t i = 0; i < quadrature_group_list.size(); ++i) {
-      auto [w, xi, eta] = quadrature_group_list[i];
+      const auto& [w, xi, eta] = quadrature_group_list[i];
       if (xi != eta) {
         Assert(xi != 0, "quadrature group was not defined correctly");
         if (eta == 0) {
@@ -321,7 +321,7 @@ EconomicalGaussQuadrature<3>::_buildPointAndWeightLists()
 
     size_t k = 0;
     for (size_t i = 0; i < quadrature_group_list.size(); ++i) {
-      auto [w, xi, eta, zeta] = quadrature_group_list[i];
+      const auto& [w, xi, eta, zeta] = quadrature_group_list[i];
       if (xi == eta) {
         if (xi == zeta) {
           if (xi == 0) {
diff --git a/src/analysis/SimplicialGaussQuadrature.cpp b/src/analysis/SimplicialGaussQuadrature.cpp
index 03ccc138e1071f977a28b03ed29550c034bf7cc6..d759050c3ace24e700bffe30e16462e8f272851a 100644
--- a/src/analysis/SimplicialGaussQuadrature.cpp
+++ b/src/analysis/SimplicialGaussQuadrature.cpp
@@ -16,15 +16,71 @@ template <>
 void
 SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
 {
+  using R2              = TinyVector<2>;
+  using QuadratureGroup = std::array<double, 4>;
+
+  auto fill_quadrature_points = [](auto quadrature_group_list, auto& point_list, auto& weight_list) {
+    Assert(point_list.size() == weight_list.size());
+
+    const R2 A = {-1, -1};
+    const R2 B = {+1, -1};
+    const R2 C = {-1, +1};
+
+    size_t k = 0;
+    for (size_t i = 0; i < quadrature_group_list.size(); ++i) {
+      const auto [half_weight, lambda_1, lambda_2, lambda_3] = quadrature_group_list[i];
+
+      const double w = 2 * half_weight;
+
+      if (lambda_2 == lambda_3) {
+        if (lambda_1 == lambda_2) {
+          point_list[k]  = lambda_1 * A + lambda_2 * B + lambda_3 * C;
+          weight_list[k] = w;
+
+          k += 1;
+        } else {
+          Assert(lambda_1 != lambda_2);
+
+          point_list[k + 0] = lambda_1 * A + lambda_2 * B + lambda_3 * C;
+          point_list[k + 1] = lambda_2 * A + lambda_3 * B + lambda_1 * C;
+          point_list[k + 2] = lambda_3 * A + lambda_1 * B + lambda_2 * C;
+
+          for (size_t i = 0; i < 3; ++i) {
+            weight_list[k + i] = w;
+          }
+
+          k += 3;
+        }
+      } else {
+        point_list[k + 0] = lambda_1 * A + lambda_2 * B + lambda_3 * C;
+        point_list[k + 1] = lambda_1 * A + lambda_3 * B + lambda_2 * C;
+        point_list[k + 2] = lambda_2 * A + lambda_1 * B + lambda_3 * C;
+        point_list[k + 3] = lambda_2 * A + lambda_3 * B + lambda_1 * C;
+        point_list[k + 4] = lambda_3 * A + lambda_1 * B + lambda_2 * C;
+        point_list[k + 5] = lambda_3 * A + lambda_2 * B + lambda_1 * C;
+
+        for (size_t i = 0; i < 6; ++i) {
+          weight_list[k + i] = w;
+        }
+
+        k += 6;
+      }
+    }
+
+    Assert(k == point_list.size(), "invalid number of quadrature points");
+  };
+
   switch (m_degree) {
   case 0:
   case 1: {
     constexpr size_t nb_points = 1;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0]  = {-0.333333333333333, -0.333333333333333};
-    weight_list[0] = 2;
+    std::array<QuadratureGroup, 1> quadrature_group_list =   //
+      {QuadratureGroup{1.000000000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -32,16 +88,13 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 2: {
     constexpr size_t nb_points = 3;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0] = {-0.666666666666667, -0.666666666666667};
-    point_list[1] = {-0.666666666666667, +0.333333333333333};
-    point_list[2] = {+0.333333333333333, -0.666666666666667};
+    std::array<QuadratureGroup, 1> quadrature_group_list =   //
+      {QuadratureGroup{0.333333333333333, 0.666666666666667, 0.166666666666667, 0.166666666666667}};
 
-    weight_list[0] = 0.666666666666667;
-    weight_list[1] = 0.666666666666667;
-    weight_list[2] = 0.666666666666667;
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -49,18 +102,14 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 3: {
     constexpr size_t nb_points = 4;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0] = {-0.333333333333333, -0.333333333333333};
-    point_list[1] = {-0.600000000000000, -0.600000000000000};
-    point_list[2] = {-0.600000000000000, +0.200000000000000};
-    point_list[3] = {+0.200000000000000, -0.600000000000000};
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{-0.562500000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{+0.520833333333333, 0.600000000000000, 0.200000000000000, 0.200000000000000}};
 
-    weight_list[0] = -1.125;
-    weight_list[1] = 1.041666666666667;
-    weight_list[2] = 1.041666666666667;
-    weight_list[3] = 1.041666666666667;
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -68,22 +117,14 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 4: {
     constexpr size_t nb_points = 6;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0] = {-0.108103018168070, -0.108103018168070};
-    point_list[1] = {-0.108103018168070, -0.783793963663860};
-    point_list[2] = {-0.783793963663860, -0.108103018168070};
-    point_list[3] = {-0.816847572980458, -0.816847572980458};
-    point_list[4] = {-0.816847572980458, +0.633695145960918};
-    point_list[5] = {+0.633695145960918, -0.816847572980458};
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.223381589678011, 0.108103018168070, 0.445948490915965, 0.445948490915965},
+       QuadratureGroup{0.109951743655322, 0.816847572980459, 0.091576213509771, 0.091576213509771}};
 
-    weight_list[0] = 0.446763179356022;
-    weight_list[1] = 0.446763179356022;
-    weight_list[2] = 0.446763179356022;
-    weight_list[3] = 0.219903487310644;
-    weight_list[4] = 0.219903487310644;
-    weight_list[5] = 0.219903487310644;
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -91,24 +132,15 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 5: {
     constexpr size_t nb_points = 7;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0] = {-0.333333333333333, -0.333333333333333};
-    point_list[1] = {-0.059715871789770, -0.059715871789770};
-    point_list[2] = {-0.059715871789770, -0.880568256420460};
-    point_list[3] = {-0.880568256420460, -0.059715871789770};
-    point_list[4] = {-0.797426985353088, -0.797426985353088};
-    point_list[5] = {-0.797426985353088, +0.594853970706174};
-    point_list[6] = {+0.594853970706174, -0.797426985353088};
-
-    weight_list[0] = 0.45;
-    weight_list[1] = 0.264788305577012;
-    weight_list[2] = 0.264788305577012;
-    weight_list[3] = 0.264788305577012;
-    weight_list[4] = 0.251878361089654;
-    weight_list[5] = 0.251878361089654;
-    weight_list[6] = 0.251878361089654;
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.225000000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.132394152788506, 0.059715871789770, 0.470142064105115, 0.470142064105115},
+       QuadratureGroup{0.125939180544827, 0.797426985353087, 0.101286507323456, 0.101286507323456}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -116,34 +148,15 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 6: {
     constexpr size_t nb_points = 12;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0]  = {-0.501426509658180, -0.501426509658180};
-    point_list[1]  = {-0.501426509658180, +0.002853019316358};
-    point_list[2]  = {+0.002853019316358, -0.501426509658180};
-    point_list[3]  = {-0.873821971016996, -0.873821971016996};
-    point_list[4]  = {-0.873821971016996, +0.747643942033992};
-    point_list[5]  = {+0.747643942033992, -0.873821971016996};
-    point_list[6]  = {-0.379295097932432, +0.273004998242798};
-    point_list[7]  = {+0.273004998242798, -0.893709900310366};
-    point_list[8]  = {-0.893709900310366, -0.379295097932432};
-    point_list[9]  = {-0.379295097932432, -0.893709900310366};
-    point_list[10] = {+0.273004998242798, -0.379295097932432};
-    point_list[11] = {-0.893709900310366, +0.273004998242798};
-
-    weight_list[0]  = 0.233572551452758;
-    weight_list[1]  = 0.233572551452758;
-    weight_list[2]  = 0.233572551452758;
-    weight_list[3]  = 0.101689812740414;
-    weight_list[4]  = 0.101689812740414;
-    weight_list[5]  = 0.101689812740414;
-    weight_list[6]  = 0.165702151236748;
-    weight_list[7]  = 0.165702151236748;
-    weight_list[8]  = 0.165702151236748;
-    weight_list[9]  = 0.165702151236748;
-    weight_list[10] = 0.165702151236748;
-    weight_list[11] = 0.165702151236748;
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.116786275726379, 0.501426509658179, 0.249286745170910, 0.249286745170910},
+       QuadratureGroup{0.050844906370207, 0.873821971016996, 0.063089014491502, 0.063089014491502},
+       QuadratureGroup{0.082851075618374, 0.053145049844817, 0.310352451033784, 0.636502499121399}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -151,36 +164,329 @@ SimplicialGaussQuadrature<2>::_buildPointAndWeightLists()
   }
   case 7: {
     constexpr size_t nb_points = 13;
-    SmallArray<TinyVector<2>> point_list(nb_points);
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{-0.149570044467682, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{+0.175615257433208, 0.479308067841920, 0.260345966079040, 0.260345966079040},
+       QuadratureGroup{+0.053347235608838, 0.869739794195568, 0.065130102902216, 0.065130102902216},
+       QuadratureGroup{+0.077113760890257, 0.048690315425316, 0.312865496004874, 0.638444188569810}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 8: {
+    constexpr size_t nb_points = 16;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.144315607677787, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.095091634267285, 0.081414823414554, 0.459292588292723, 0.459292588292723},
+       QuadratureGroup{0.103217370534718, 0.658861384496480, 0.170569307751760, 0.170569307751760},
+       QuadratureGroup{0.032458497623198, 0.898905543365938, 0.050547228317031, 0.050547228317031},
+       QuadratureGroup{0.027230314174435, 0.008394777409958, 0.263112829634638, 0.728492392955404}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 9: {
+    constexpr size_t nb_points = 19;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.097135796282799, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.031334700227139, 0.020634961602525, 0.489682519198738, 0.489682519198738},
+       QuadratureGroup{0.077827541004774, 0.125820817014127, 0.437089591492937, 0.437089591492937},
+       QuadratureGroup{0.079647738927210, 0.623592928761935, 0.188203535619033, 0.188203535619033},
+       QuadratureGroup{0.025577675658698, 0.910540973211095, 0.044729513394453, 0.044729513394453},
+       QuadratureGroup{0.043283539377289, 0.036838412054736, 0.221962989160766, 0.741198598784498}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 10: {
+    constexpr size_t nb_points = 25;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.090817990382754, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.036725957756467, 0.028844733232685, 0.485577633383657, 0.485577633383657},
+       QuadratureGroup{0.045321059435528, 0.781036849029926, 0.109481575485037, 0.109481575485037},
+       QuadratureGroup{0.072757916845420, 0.141707219414880, 0.307939838764121, 0.550352941820999},
+       QuadratureGroup{0.028327242531057, 0.025003534762686, 0.246672560639903, 0.728323904597411},
+       QuadratureGroup{0.009421666963733, 0.009540815400299, 0.066803251012200, 0.923655933587500}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 11: {
+    constexpr size_t nb_points = 27;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.000927006328961, -0.069222096541517, +0.534611048270758, +0.534611048270758},
+       QuadratureGroup{0.077149534914813, +0.202061394068290, +0.398969302965855, +0.398969302965855},
+       QuadratureGroup{0.059322977380774, +0.593380199137435, +0.203309900431282, +0.203309900431282},
+       QuadratureGroup{0.036184540503418, +0.761298175434837, +0.119350912282581, +0.119350912282581},
+       QuadratureGroup{0.013659731002678, +0.935270103777448, +0.032364948111276, +0.032364948111276},
+       QuadratureGroup{0.052337111962204, +0.050178138310495, +0.356620648261293, +0.593201213428213},
+       QuadratureGroup{0.020707659639141, +0.021022016536166, +0.171488980304042, +0.807489003159792}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 12: {
+    constexpr size_t nb_points = 33;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.025731066440455, 0.023565220452390, 0.488217389773805, 0.488217389773805},
+       QuadratureGroup{0.043692544538038, 0.120551215411079, 0.439724392294460, 0.439724392294460},
+       QuadratureGroup{0.062858224217885, 0.457579229975768, 0.271210385012116, 0.271210385012116},
+       QuadratureGroup{0.034796112930709, 0.744847708916828, 0.127576145541586, 0.127576145541586},
+       QuadratureGroup{0.006166261051559, 0.957365299093579, 0.021317350453210, 0.021317350453210},
+       QuadratureGroup{0.040371557766381, 0.115343494534698, 0.275713269685514, 0.608943235779788},
+       QuadratureGroup{0.022356773202303, 0.022838332222257, 0.281325580989940, 0.695836086787803},
+       QuadratureGroup{0.017316231108659, 0.025734050548330, 0.116251915907597, 0.858014033544073}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 13: {
+    constexpr size_t nb_points = 37;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.052520923400802, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.011280145209330, 0.009903630120591, 0.495048184939705, 0.495048184939705},
+       QuadratureGroup{0.031423518362454, 0.062566729780852, 0.468716635109574, 0.468716635109574},
+       QuadratureGroup{0.047072502504194, 0.170957326397447, 0.414521336801277, 0.414521336801277},
+       QuadratureGroup{0.047363586536355, 0.541200855914337, 0.229399572042831, 0.229399572042831},
+       QuadratureGroup{0.031167529045794, 0.771151009607340, 0.114424495196330, 0.114424495196330},
+       QuadratureGroup{0.007975771465074, 0.950377217273082, 0.024811391363459, 0.024811391363459},
+       QuadratureGroup{0.036848402728732, 0.094853828379579, 0.268794997058761, 0.636351174561660},
+       QuadratureGroup{0.017401463303822, 0.018100773278807, 0.291730066734288, 0.690169159986905},
+       QuadratureGroup{0.015521786839045, 0.022233076674090, 0.126357385491669, 0.851409537834241}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 14: {
+    constexpr size_t nb_points = 42;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.021883581369429, 0.022072179275643, 0.488963910362179, 0.488963910362179},
+       QuadratureGroup{0.032788353544125, 0.164710561319092, 0.417644719340454, 0.417644719340454},
+       QuadratureGroup{0.051774104507292, 0.453044943382323, 0.273477528308839, 0.273477528308839},
+       QuadratureGroup{0.042162588736993, 0.645588935174913, 0.177205532412543, 0.177205532412543},
+       QuadratureGroup{0.014433699669777, 0.876400233818255, 0.061799883090873, 0.061799883090873},
+       QuadratureGroup{0.004923403602400, 0.961218077502598, 0.019390961248701, 0.019390961248701},
+       QuadratureGroup{0.024665753212564, 0.057124757403648, 0.172266687821356, 0.770608554774996},
+       QuadratureGroup{0.038571510787061, 0.092916249356972, 0.336861459796345, 0.570222290846683},
+       QuadratureGroup{0.014436308113534, 0.014646950055654, 0.298372882136258, 0.686980167808088},
+       QuadratureGroup{0.005010228838501, 0.001268330932872, 0.118974497696957, 0.879757171370171}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 15: {
+    constexpr size_t nb_points = 48;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.001916875642849, -0.013945833716486, +0.506972916858243, +0.506972916858243},
+       QuadratureGroup{0.044249027271145, +0.137187291433955, +0.431406354283023, +0.431406354283023},
+       QuadratureGroup{0.051186548718852, +0.444612710305711, +0.277693644847144, +0.277693644847144},
+       QuadratureGroup{0.023687735870688, +0.747070217917492, +0.126464891041254, +0.126464891041254},
+       QuadratureGroup{0.013289775690021, +0.858383228050628, +0.070808385974686, +0.070808385974686},
+       QuadratureGroup{0.004748916608192, +0.962069659517853, +0.018965170241073, +0.018965170241073},
+       QuadratureGroup{0.038550072599593, +0.133734161966621, +0.261311371140087, +0.604954466893291},
+       QuadratureGroup{0.027215814320624, +0.036366677396917, +0.388046767090269, +0.575586555512814},
+       QuadratureGroup{0.002182077366797, -0.010174883126571, +0.285712220049916, +0.724462663076655},
+       QuadratureGroup{0.021505319847731, +0.036843869875878, +0.215599664072284, +0.747556466051838},
+       QuadratureGroup{0.007673942631049, +0.012459809331199, +0.103575616576386, +0.883964574092416}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 16: {
+    constexpr size_t nb_points = 52;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.046875697427642, +0.333333333333333, +0.333333333333333, +0.333333333333333},
+       QuadratureGroup{0.006405878578585, +0.005238916103123, +0.497380541948438, +0.497380541948438},
+       QuadratureGroup{0.041710296739387, +0.173061122901295, +0.413469438549352, +0.413469438549352},
+       QuadratureGroup{0.026891484250064, +0.059082801866017, +0.470458599066991, +0.470458599066991},
+       QuadratureGroup{0.042132522761650, +0.518892500060958, +0.240553749969521, +0.240553749969521},
+       QuadratureGroup{0.030000266842773, +0.704068411554854, +0.147965794222573, +0.147965794222573},
+       QuadratureGroup{0.014200098925024, +0.849069624685052, +0.075465187657474, +0.075465187657474},
+       QuadratureGroup{0.003582462351273, +0.966807194753950, +0.016596402623025, +0.016596402623025},
+       QuadratureGroup{0.032773147460627, +0.103575692245252, +0.296555596579887, +0.599868711174861},
+       QuadratureGroup{0.015298306248441, +0.020083411655416, +0.337723063403079, +0.642193524941505},
+       QuadratureGroup{0.002386244192839, -0.004341002614139, +0.204748281642812, +0.799592720971327},
+       QuadratureGroup{0.019084792755899, +0.041941786468010, +0.189358492130623, +0.768699721401368},
+       QuadratureGroup{0.006850054546542, +0.014317320230681, +0.085283615682657, +0.900399064086661}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 17: {
+    constexpr size_t nb_points = 61;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.033437199290803, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.005093415440507, 0.005658918886452, 0.497170540556774, 0.497170540556774},
+       QuadratureGroup{0.014670864527638, 0.035647354750751, 0.482176322624625, 0.482176322624625},
+       QuadratureGroup{0.024350878353672, 0.099520061958437, 0.450239969020782, 0.450239969020782},
+       QuadratureGroup{0.031107550868969, 0.199467521245206, 0.400266239377397, 0.400266239377397},
+       QuadratureGroup{0.031257111218620, 0.495717464058095, 0.252141267970953, 0.252141267970953},
+       QuadratureGroup{0.024815654339665, 0.675905990683077, 0.162047004658461, 0.162047004658461},
+       QuadratureGroup{0.014056073070557, 0.848248235478508, 0.075875882260746, 0.075875882260746},
+       QuadratureGroup{0.003194676173779, 0.968690546064356, 0.015654726967822, 0.015654726967822},
+       QuadratureGroup{0.008119655318993, 0.010186928826919, 0.334319867363658, 0.655493203809423},
+       QuadratureGroup{0.026805742283163, 0.135440871671036, 0.292221537796944, 0.572337590532020},
+       QuadratureGroup{0.018459993210822, 0.054423924290583, 0.319574885423190, 0.626001190286228},
+       QuadratureGroup{0.008476868534328, 0.012868560833637, 0.190704224192292, 0.796427214974071},
+       QuadratureGroup{0.018292796770025, 0.067165782413524, 0.180483211648746, 0.752351005937729},
+       QuadratureGroup{0.006665632004165, 0.014663182224828, 0.080711313679564, 0.904625504095608}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 18: {
+    constexpr size_t nb_points = 70;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{+0.030809939937647, +0.333333333333333, +0.333333333333333, +0.333333333333333},
+       QuadratureGroup{+0.009072436679404, +0.013310382738157, +0.493344808630921, +0.493344808630921},
+       QuadratureGroup{+0.018761316939594, +0.061578811516086, +0.469210594241957, +0.469210594241957},
+       QuadratureGroup{+0.019441097985477, +0.127437208225989, +0.436281395887006, +0.436281395887006},
+       QuadratureGroup{+0.027753948610810, +0.210307658653168, +0.394846170673416, +0.394846170673416},
+       QuadratureGroup{+0.032256225351457, +0.500410862393686, +0.249794568803157, +0.249794568803157},
+       QuadratureGroup{+0.025074032616922, +0.677135612512315, +0.161432193743843, +0.161432193743843},
+       QuadratureGroup{+0.015271927971832, +0.846803545029257, +0.076598227485371, +0.076598227485371},
+       QuadratureGroup{+0.006793922022963, +0.951495121293100, +0.024252439353450, +0.024252439353450},
+       QuadratureGroup{-0.002223098729920, +0.913707265566071, +0.043146367216965, +0.043146367216965},
+       QuadratureGroup{+0.006331914076406, +0.008430536202420, +0.358911494940944, +0.632657968856636},
+       QuadratureGroup{+0.027257538049138, +0.131186551737188, +0.294402476751957, +0.574410971510855},
+       QuadratureGroup{+0.017676785649465, +0.050203151565675, +0.325017801641814, +0.624779046792512},
+       QuadratureGroup{+0.018379484638070, +0.066329263810916, +0.184737559666046, +0.748933176523037},
+       QuadratureGroup{+0.008104732808192, +0.011996194566236, +0.218796800013321, +0.769207005420443},
+       QuadratureGroup{+0.007634129070725, +0.014858100590125, +0.101179597136408, +0.883962302273467},
+       QuadratureGroup{+0.000046187660794, -0.035222015287949, +0.020874755282586, +1.014347260005363}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 19: {
+    constexpr size_t nb_points = 73;
+    SmallArray<R2> point_list(nb_points);
+    SmallArray<double> weight_list(nb_points);
+
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{0.032906331388919, 0.333333333333333, 0.333333333333333, 0.333333333333333},
+       QuadratureGroup{0.010330731891272, 0.020780025853987, 0.489609987073006, 0.489609987073006},
+       QuadratureGroup{0.022387247263016, 0.090926214604215, 0.454536892697893, 0.454536892697893},
+       QuadratureGroup{0.030266125869468, 0.197166638701138, 0.401416680649431, 0.401416680649431},
+       QuadratureGroup{0.030490967802198, 0.488896691193805, 0.255551654403098, 0.255551654403098},
+       QuadratureGroup{0.024159212741641, 0.645844115695741, 0.177077942152130, 0.177077942152130},
+       QuadratureGroup{0.016050803586801, 0.779877893544096, 0.110061053227952, 0.110061053227952},
+       QuadratureGroup{0.008084580261784, 0.888942751496321, 0.055528624251840, 0.055528624251840},
+       QuadratureGroup{0.002079362027485, 0.974756272445543, 0.012621863777229, 0.012621863777229},
+       QuadratureGroup{0.003884876904981, 0.003611417848412, 0.395754787356943, 0.600633794794645},
+       QuadratureGroup{0.025574160612022, 0.134466754530780, 0.307929983880436, 0.557603261588784},
+       QuadratureGroup{0.008880903573338, 0.014446025776115, 0.264566948406520, 0.720987025817365},
+       QuadratureGroup{0.016124546761731, 0.046933578838178, 0.358539352205951, 0.594527068955871},
+       QuadratureGroup{0.002491941817491, 0.002861120350567, 0.157807405968595, 0.839331473680839},
+       QuadratureGroup{0.018242840118951, 0.223861424097916, 0.075050596975911, 0.701087978926173},
+       QuadratureGroup{0.010258563736199, 0.034647074816760, 0.142421601113383, 0.822931324069857},
+       QuadratureGroup{0.003799928855302, 0.010161119296278, 0.065494628082938, 0.924344252620784}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+
+    m_point_list  = point_list;
+    m_weight_list = weight_list;
+    break;
+  }
+  case 20: {
+    constexpr size_t nb_points = 79;
+    SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    point_list[0]  = {-0.333333333333333, -0.333333333333333};
-    point_list[1]  = {-0.479308067841920, -0.479308067841920};
-    point_list[2]  = {-0.479308067841920, -0.041383864316160};
-    point_list[3]  = {-0.041383864316160, -0.479308067841920};
-    point_list[4]  = {-0.869739794195568, -0.869739794195568};
-    point_list[5]  = {-0.869739794195568, +0.739479588391136};
-    point_list[6]  = {+0.739479588391136, -0.869739794195568};
-    point_list[7]  = {-0.374269007990252, +0.276888377139620};
-    point_list[8]  = {+0.276888377139620, -0.902619369149368};
-    point_list[9]  = {-0.902619369149368, -0.374269007990252};
-    point_list[10] = {-0.374269007990252, -0.902619369149368};
-    point_list[11] = {+0.276888377139620, -0.374269007990252};
-    point_list[12] = {-0.902619369149368, +0.276888377139620};
-
-    weight_list[0]  = -0.299140088935364;
-    weight_list[1]  = +0.351230514866416;
-    weight_list[2]  = +0.351230514866416;
-    weight_list[3]  = +0.351230514866416;
-    weight_list[4]  = +0.106694471217676;
-    weight_list[5]  = +0.106694471217676;
-    weight_list[6]  = +0.106694471217676;
-    weight_list[7]  = +0.154227521780514;
-    weight_list[8]  = +0.154227521780514;
-    weight_list[9]  = +0.154227521780514;
-    weight_list[10] = +0.154227521780514;
-    weight_list[11] = +0.154227521780514;
-    weight_list[12] = +0.154227521780514;
+    std::array quadrature_group_list =   //
+      {QuadratureGroup{+0.033057055541624, +0.333333333333333, +0.333333333333333, +0.333333333333333},
+       QuadratureGroup{+0.000867019185663, -0.001900928704400, +0.500950464352200, +0.500950464352200},
+       QuadratureGroup{+0.011660052716448, +0.023574084130543, +0.488212957934729, +0.488212957934729},
+       QuadratureGroup{+0.022876936356421, +0.089726636099435, +0.455136681950283, +0.455136681950283},
+       QuadratureGroup{+0.030448982673938, +0.196007481363421, +0.401996259318289, +0.401996259318289},
+       QuadratureGroup{+0.030624891725355, +0.488214180481157, +0.255892909759421, +0.255892909759421},
+       QuadratureGroup{+0.024368057676800, +0.647023488009788, +0.176488255995106, +0.176488255995106},
+       QuadratureGroup{+0.015997432032024, +0.791658289326483, +0.104170855336758, +0.104170855336758},
+       QuadratureGroup{+0.007698301815602, +0.893862072318140, +0.053068963840930, +0.053068963840930},
+       QuadratureGroup{-0.000632060497488, +0.916762569607942, +0.041618715196029, +0.041618715196029},
+       QuadratureGroup{+0.001751134301193, +0.976836157186356, +0.011581921406822, +0.011581921406822},
+       QuadratureGroup{+0.016465839189576, +0.048741583664839, +0.344855770229001, +0.606402646106160},
+       QuadratureGroup{+0.004839033540485, +0.006314115948605, +0.377843269594854, +0.615842614456541},
+       QuadratureGroup{+0.025804906534650, +0.134316520547348, +0.306635479062357, +0.559048000390295},
+       QuadratureGroup{+0.008471091054441, +0.013973893962392, +0.249419362774742, +0.736606743262866},
+       QuadratureGroup{+0.018354914106280, +0.075549132909764, +0.212775724802802, +0.711675142287434},
+       QuadratureGroup{+0.000704404677908, -0.008368153208227, +0.146965436053239, +0.861402717154987},
+       QuadratureGroup{+0.010112684927462, +0.026686063258714, +0.137726978828923, +0.835586957912363},
+       QuadratureGroup{+0.003573909385950, +0.010547719294141, +0.059696109149007, +0.929756171556853}};
+
+    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;