diff --git a/src/analysis/TriangleGaussQuadrature.cpp b/src/analysis/TriangleGaussQuadrature.cpp
index 2086712713df036601dcd694b7c481fbfa490dc3..13ae4693d69ac2cb0497950cbcece93cd4505170 100644
--- a/src/analysis/TriangleGaussQuadrature.cpp
+++ b/src/analysis/TriangleGaussQuadrature.cpp
@@ -4,8 +4,14 @@
 void
 TriangleGaussQuadrature::_buildPointAndWeightLists()
 {
-  using R2         = TinyVector<2>;
-  using Descriptor = std::array<double, 4>;
+  using R2 = TinyVector<2>;
+
+  struct Descriptor
+  {
+    int id;
+    double weight;
+    std::vector<double> lambda_list;
+  };
 
   auto fill_quadrature_points = [](auto descriptor_list, auto& point_list, auto& weight_list) {
     Assert(point_list.size() == weight_list.size());
@@ -16,42 +22,57 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
 
     size_t k = 0;
     for (size_t i = 0; i < descriptor_list.size(); ++i) {
-      const auto [unit_weight, lambda_1, lambda_2, lambda_3] = descriptor_list[i];
-
-      const double w = 0.5 * unit_weight;
+      const auto [id, w, position_list] = descriptor_list[i];
 
-      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;
+      switch (id) {
+      case 1: {
+        Assert(position_list.size() == 0);
 
-          k += 1;
-        } else {
-          Assert(lambda_1 != lambda_2);
+        point_list[k]  = (1. / 3) * (A + B + C);
+        weight_list[k] = w;
 
-          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;
+        k += 1;
+        break;
+      }
+      case 2: {
+        Assert(position_list.size() == 1);
+        const double& l0 = position_list[0];
+        const double& l1 = 1 - 2 * l0;
 
-          for (size_t i = 0; i < 3; ++i) {
-            weight_list[k + i] = w;
-          }
+        point_list[k + 0] = l0 * A + l0 * B + l1 * C;
+        point_list[k + 1] = l0 * A + l1 * B + l0 * C;
+        point_list[k + 2] = l1 * A + l0 * B + l0 * C;
 
-          k += 3;
+        for (size_t l = 0; l < 3; ++l) {
+          weight_list[k + l] = w;
         }
-      } 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 += 3;
+        break;
+      }
+      case 3: {
+        Assert(position_list.size() == 2);
+        const double& l0 = position_list[0];
+        const double& l1 = position_list[1];
+        const double& l2 = 1 - l0 - l1;
+
+        point_list[k + 0] = l0 * A + l1 * B + l2 * C;
+        point_list[k + 1] = l0 * A + l2 * B + l1 * C;
+        point_list[k + 2] = l1 * A + l0 * B + l2 * C;
+        point_list[k + 3] = l1 * A + l2 * B + l0 * C;
+        point_list[k + 4] = l2 * A + l0 * B + l1 * C;
+        point_list[k + 5] = l2 * A + l1 * B + l0 * C;
+
+        for (size_t l = 0; l < 6; ++l) {
+          weight_list[k + l] = w;
         }
 
         k += 6;
+        break;
+      }
+      default: {
+        throw UnexpectedError("invalid quadrature id");
+      }
       }
     }
 
@@ -65,8 +86,8 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array<Descriptor, 1> descriptor_list =   //
-      {Descriptor{1.000000000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333}};
+    std::array descriptor_list =   //
+      {Descriptor{1, 5.000000000000000e-01, {}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -79,23 +100,8 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array<Descriptor, 1> descriptor_list =   //
-      {Descriptor{0.333333333333333, 0.666666666666667, 0.166666666666667, 0.166666666666667}};
-
-    fill_quadrature_points(descriptor_list, point_list, weight_list);
-
-    m_point_list  = point_list;
-    m_weight_list = weight_list;
-    break;
-  }
-  case 3: {
-    constexpr size_t nb_points = 4;
-    SmallArray<R2> point_list(nb_points);
-    SmallArray<double> weight_list(nb_points);
-
     std::array descriptor_list =   //
-      {Descriptor{-0.562500000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{+0.520833333333333, 0.600000000000000, 0.200000000000000, 0.200000000000000}};
+      {Descriptor{2, 1.666666666666667e-01, {1.666666666666667e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -103,17 +109,17 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     m_weight_list = weight_list;
     break;
   }
+  case 3:
   case 4: {
     constexpr size_t nb_points = 6;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.223381589678011, 0.108103018168070, 0.445948490915965, 0.445948490915965},
-       Descriptor{0.109951743655322, 0.816847572980459, 0.091576213509771, 0.091576213509771}};
+      {Descriptor{2, 1.116907948390057e-01, {4.459484909159649e-01}},
+       Descriptor{2, 5.497587182766094e-02, {9.157621350977074e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
-
     m_point_list  = point_list;
     m_weight_list = weight_list;
     break;
@@ -124,9 +130,8 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.225000000000000, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.132394152788506, 0.059715871789770, 0.470142064105115, 0.470142064105115},
-       Descriptor{0.125939180544827, 0.797426985353087, 0.101286507323456, 0.101286507323456}};
+      {Descriptor{1, 1.125000000000000e-01, {}}, Descriptor{2, 6.296959027241357e-02, {1.012865073234563e-01}},
+       Descriptor{2, 6.619707639425310e-02, {4.701420641051151e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -140,9 +145,9 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.116786275726379, 0.501426509658179, 0.249286745170910, 0.249286745170910},
-       Descriptor{0.050844906370207, 0.873821971016996, 0.063089014491502, 0.063089014491502},
-       Descriptor{0.082851075618374, 0.053145049844817, 0.310352451033784, 0.636502499121399}};
+      {Descriptor{2, 2.542245318510341e-02, {6.308901449150223e-02}},
+       Descriptor{2, 5.839313786318968e-02, {2.492867451709104e-01}},
+       Descriptor{3, 4.142553780918679e-02, {3.103524510337844e-01, 5.314504984481695e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -151,15 +156,15 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 7: {
-    constexpr size_t nb_points = 13;
+    constexpr size_t nb_points = 15;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{-0.149570044467682, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{+0.175615257433208, 0.479308067841920, 0.260345966079040, 0.260345966079040},
-       Descriptor{+0.053347235608838, 0.869739794195568, 0.065130102902216, 0.065130102902216},
-       Descriptor{+0.077113760890257, 0.048690315425316, 0.312865496004874, 0.638444188569810}};
+      {Descriptor{2, 8.272525055396066e-03, {3.373064855458785e-02}},
+       Descriptor{2, 6.397208561507779e-02, {2.415773825954036e-01}},
+       Descriptor{2, 3.854332309299303e-02, {4.743096925047182e-01}},
+       Descriptor{3, 2.793936645159989e-02, {1.986833147973516e-01, 4.703664465259523e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -173,11 +178,10 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.144315607677787, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.095091634267285, 0.081414823414554, 0.459292588292723, 0.459292588292723},
-       Descriptor{0.103217370534718, 0.658861384496480, 0.170569307751760, 0.170569307751760},
-       Descriptor{0.032458497623198, 0.898905543365938, 0.050547228317031, 0.050547228317031},
-       Descriptor{0.027230314174435, 0.008394777409958, 0.263112829634638, 0.728492392955404}};
+      {Descriptor{1, 7.215780383889359e-02, {}}, Descriptor{2, 4.754581713364231e-02, {4.592925882927232e-01}},
+       Descriptor{2, 5.160868526735912e-02, {1.705693077517602e-01}},
+       Descriptor{2, 1.622924881159904e-02, {5.054722831703098e-02}},
+       Descriptor{3, 1.361515708721750e-02, {2.631128296346381e-01, 8.394777409957605e-03}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -191,12 +195,12 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.097135796282799, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.031334700227139, 0.020634961602525, 0.489682519198738, 0.489682519198738},
-       Descriptor{0.077827541004774, 0.125820817014127, 0.437089591492937, 0.437089591492937},
-       Descriptor{0.079647738927210, 0.623592928761935, 0.188203535619033, 0.188203535619033},
-       Descriptor{0.025577675658698, 0.910540973211095, 0.044729513394453, 0.044729513394453},
-       Descriptor{0.043283539377289, 0.036838412054736, 0.221962989160766, 0.741198598784498}};
+      {Descriptor{1, 4.856789814139942e-02, {}},
+       Descriptor{2, 3.891377050238714e-02, {4.370895914929366e-01}},
+       Descriptor{2, 3.982386946360512e-02, {1.882035356190327e-01}},
+       Descriptor{2, 1.566735011356954e-02, {4.896825191987376e-01}},
+       Descriptor{2, 1.278883782934902e-02, {4.472951339445271e-02}},
+       Descriptor{3, 2.164176968864469e-02, {2.219629891607657e-01, 3.683841205473629e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -210,12 +214,12 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.090817990382754, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.036725957756467, 0.028844733232685, 0.485577633383657, 0.485577633383657},
-       Descriptor{0.045321059435528, 0.781036849029926, 0.109481575485037, 0.109481575485037},
-       Descriptor{0.072757916845420, 0.141707219414880, 0.307939838764121, 0.550352941820999},
-       Descriptor{0.028327242531057, 0.025003534762686, 0.246672560639903, 0.728323904597411},
-       Descriptor{0.009421666963733, 0.009540815400299, 0.066803251012200, 0.923655933587500}};
+      {Descriptor{1, 4.087166457314299e-02, {}},
+       Descriptor{2, 6.676484406574783e-03, {3.205537321694351e-02}},
+       Descriptor{2, 2.297898180237237e-02, {1.421611010565644e-01}},
+       Descriptor{3, 3.195245319821202e-02, {1.481328857838206e-01, 3.218129952888354e-01}},
+       Descriptor{3, 1.709232408147971e-02, {3.691467818278110e-01, 2.961988948872977e-02}},
+       Descriptor{3, 1.264887885364419e-02, {1.637017337371825e-01, 2.836766533993844e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -224,18 +228,19 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 11: {
-    constexpr size_t nb_points = 27;
+    constexpr size_t nb_points = 28;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.000927006328961, -0.069222096541517, +0.534611048270758, +0.534611048270758},
-       Descriptor{0.077149534914813, +0.202061394068290, +0.398969302965855, +0.398969302965855},
-       Descriptor{0.059322977380774, +0.593380199137435, +0.203309900431282, +0.203309900431282},
-       Descriptor{0.036184540503418, +0.761298175434837, +0.119350912282581, +0.119350912282581},
-       Descriptor{0.013659731002678, +0.935270103777448, +0.032364948111276, +0.032364948111276},
-       Descriptor{0.052337111962204, +0.050178138310495, +0.356620648261293, +0.593201213428213},
-       Descriptor{0.020707659639141, +0.021022016536166, +0.171488980304042, +0.807489003159792}};
+      {Descriptor{1, 4.288058986611211e-02, {}},
+       Descriptor{2, 5.215935256447348e-03, {2.848541761437191e-02}},
+       Descriptor{2, 3.525784205585829e-02, {2.102199567031783e-01}},
+       Descriptor{2, 1.931537961850966e-02, {1.026354827122464e-01}},
+       Descriptor{2, 8.303136527292684e-03, {4.958919009658909e-01}},
+       Descriptor{2, 3.365807703973415e-02, {4.384659267643522e-01}},
+       Descriptor{3, 5.145144786476639e-03, {7.325427686064452e-03, 1.493247886520824e-01}},
+       Descriptor{3, 2.016623832025028e-02, {2.895811256377058e-01, 4.601050016542996e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -249,14 +254,14 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.025731066440455, 0.023565220452390, 0.488217389773805, 0.488217389773805},
-       Descriptor{0.043692544538038, 0.120551215411079, 0.439724392294460, 0.439724392294460},
-       Descriptor{0.062858224217885, 0.457579229975768, 0.271210385012116, 0.271210385012116},
-       Descriptor{0.034796112930709, 0.744847708916828, 0.127576145541586, 0.127576145541586},
-       Descriptor{0.006166261051559, 0.957365299093579, 0.021317350453210, 0.021317350453210},
-       Descriptor{0.040371557766381, 0.115343494534698, 0.275713269685514, 0.608943235779788},
-       Descriptor{0.022356773202303, 0.022838332222257, 0.281325580989940, 0.695836086787803},
-       Descriptor{0.017316231108659, 0.025734050548330, 0.116251915907597, 0.858014033544073}};
+      {Descriptor{2, 1.213341904072602e-02, {4.882037509455415e-01}},
+       Descriptor{2, 1.424302603443877e-02, {1.092578276593543e-01}},
+       Descriptor{2, 3.127060659795138e-02, {2.714625070149261e-01}},
+       Descriptor{2, 3.965821254986819e-03, {2.464636343633559e-02}},
+       Descriptor{2, 2.495916746403047e-02, {4.401116486585931e-01}},
+       Descriptor{3, 1.089179251930378e-02, {2.303415635526714e-02, 2.916556797383409e-01}},
+       Descriptor{3, 2.161368182970710e-02, {2.554542286385174e-01, 1.162960196779266e-01}},
+       Descriptor{3, 7.541838788255719e-03, {2.138249025617059e-02, 8.513377925102400e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -270,16 +275,15 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.052520923400802, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.011280145209330, 0.009903630120591, 0.495048184939705, 0.495048184939705},
-       Descriptor{0.031423518362454, 0.062566729780852, 0.468716635109574, 0.468716635109574},
-       Descriptor{0.047072502504194, 0.170957326397447, 0.414521336801277, 0.414521336801277},
-       Descriptor{0.047363586536355, 0.541200855914337, 0.229399572042831, 0.229399572042831},
-       Descriptor{0.031167529045794, 0.771151009607340, 0.114424495196330, 0.114424495196330},
-       Descriptor{0.007975771465074, 0.950377217273082, 0.024811391363459, 0.024811391363459},
-       Descriptor{0.036848402728732, 0.094853828379579, 0.268794997058761, 0.636351174561660},
-       Descriptor{0.017401463303822, 0.018100773278807, 0.291730066734288, 0.690169159986905},
-       Descriptor{0.015521786839045, 0.022233076674090, 0.126357385491669, 0.851409537834241}};
+      {Descriptor{1, 3.398001829341582e-02, {}},
+       Descriptor{2, 1.199720096444737e-02, {4.890769464525394e-01}},
+       Descriptor{2, 2.913924255959999e-02, {2.213722862918329e-01}},
+       Descriptor{2, 2.780098376522666e-02, {4.269414142598004e-01}},
+       Descriptor{2, 3.026168551769586e-03, {2.150968110884318e-02}},
+       Descriptor{3, 1.208951990579691e-02, {1.635974010678505e-01, 8.789548303219732e-02}},
+       Descriptor{3, 7.482700552582834e-03, {2.437018690109383e-02, 1.109220428034634e-01}},
+       Descriptor{3, 1.732063807042419e-02, {6.801224355420665e-02, 3.084417608921178e-01}},
+       Descriptor{3, 4.795340501771632e-03, {2.725158177734296e-01, 5.126389102382369e-03}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -293,16 +297,16 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.021883581369429, 0.022072179275643, 0.488963910362179, 0.488963910362179},
-       Descriptor{0.032788353544125, 0.164710561319092, 0.417644719340454, 0.417644719340454},
-       Descriptor{0.051774104507292, 0.453044943382323, 0.273477528308839, 0.273477528308839},
-       Descriptor{0.042162588736993, 0.645588935174913, 0.177205532412543, 0.177205532412543},
-       Descriptor{0.014433699669777, 0.876400233818255, 0.061799883090873, 0.061799883090873},
-       Descriptor{0.004923403602400, 0.961218077502598, 0.019390961248701, 0.019390961248701},
-       Descriptor{0.024665753212564, 0.057124757403648, 0.172266687821356, 0.770608554774996},
-       Descriptor{0.038571510787061, 0.092916249356972, 0.336861459796345, 0.570222290846683},
-       Descriptor{0.014436308113534, 0.014646950055654, 0.298372882136258, 0.686980167808088},
-       Descriptor{0.005010228838501, 0.001268330932872, 0.118974497696957, 0.879757171370171}};
+      {Descriptor{2, 2.108129436849651e-02, {1.772055324125434e-01}},
+       Descriptor{2, 1.639417677206267e-02, {4.176447193404539e-01}},
+       Descriptor{2, 7.216849834888334e-03, {6.179988309087260e-02}},
+       Descriptor{2, 1.094179068471444e-02, {4.889639103621786e-01}},
+       Descriptor{2, 2.588705225364579e-02, {2.734775283088386e-01}},
+       Descriptor{2, 2.461701801200041e-03, {1.939096124870105e-02}},
+       Descriptor{3, 7.218154056766920e-03, {1.464695005565441e-02, 2.983728821362577e-01}},
+       Descriptor{3, 1.233287660628184e-02, {1.722666878213556e-01, 5.712475740364794e-02}},
+       Descriptor{3, 1.928575539353034e-02, {9.291624935697182e-02, 3.368614597963450e-01}},
+       Descriptor{3, 2.505114419250336e-03, {1.189744976969569e-01, 1.268330932872025e-03}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -311,22 +315,23 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 15: {
-    constexpr size_t nb_points = 48;
+    constexpr size_t nb_points = 49;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.001916875642849, -0.013945833716486, +0.506972916858243, +0.506972916858243},
-       Descriptor{0.044249027271145, +0.137187291433955, +0.431406354283023, +0.431406354283023},
-       Descriptor{0.051186548718852, +0.444612710305711, +0.277693644847144, +0.277693644847144},
-       Descriptor{0.023687735870688, +0.747070217917492, +0.126464891041254, +0.126464891041254},
-       Descriptor{0.013289775690021, +0.858383228050628, +0.070808385974686, +0.070808385974686},
-       Descriptor{0.004748916608192, +0.962069659517853, +0.018965170241073, +0.018965170241073},
-       Descriptor{0.038550072599593, +0.133734161966621, +0.261311371140087, +0.604954466893291},
-       Descriptor{0.027215814320624, +0.036366677396917, +0.388046767090269, +0.575586555512814},
-       Descriptor{0.002182077366797, -0.010174883126571, +0.285712220049916, +0.724462663076655},
-       Descriptor{0.021505319847731, +0.036843869875878, +0.215599664072284, +0.747556466051838},
-       Descriptor{0.007673942631049, +0.012459809331199, +0.103575616576386, +0.883964574092416}};
+      {Descriptor{1, 2.216769369109204e-02, {}},
+       Descriptor{2, 2.135689078573028e-02, {4.053622141339755e-01}},
+       Descriptor{2, 8.222368781312581e-03, {7.017355289998602e-02}},
+       Descriptor{2, 8.698074000381707e-03, {4.741706814380198e-01}},
+       Descriptor{2, 2.339168086435481e-02, {2.263787134203496e-01}},
+       Descriptor{2, 4.786923091230043e-03, {4.949969567691262e-01}},
+       Descriptor{2, 1.480387318952688e-03, {1.581172625098864e-02}},
+       Descriptor{3, 7.801286415287982e-03, {3.146482428124508e-01, 1.837611238568109e-02}},
+       Descriptor{3, 2.014926686009050e-03, {7.094860523645553e-02, 9.139237037308396e-03}},
+       Descriptor{3, 1.436029346260067e-02, {9.424205359215536e-02, 1.905355894763939e-01}},
+       Descriptor{3, 5.836310590787923e-03, {1.863871372816638e-02, 1.680686452224144e-01}},
+       Descriptor{3, 1.565773814248464e-02, {9.579672364760859e-02, 3.389506114752772e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -335,24 +340,24 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 16: {
-    constexpr size_t nb_points = 52;
+    constexpr size_t nb_points = 55;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.046875697427642, +0.333333333333333, +0.333333333333333, +0.333333333333333},
-       Descriptor{0.006405878578585, +0.005238916103123, +0.497380541948438, +0.497380541948438},
-       Descriptor{0.041710296739387, +0.173061122901295, +0.413469438549352, +0.413469438549352},
-       Descriptor{0.026891484250064, +0.059082801866017, +0.470458599066991, +0.470458599066991},
-       Descriptor{0.042132522761650, +0.518892500060958, +0.240553749969521, +0.240553749969521},
-       Descriptor{0.030000266842773, +0.704068411554854, +0.147965794222573, +0.147965794222573},
-       Descriptor{0.014200098925024, +0.849069624685052, +0.075465187657474, +0.075465187657474},
-       Descriptor{0.003582462351273, +0.966807194753950, +0.016596402623025, +0.016596402623025},
-       Descriptor{0.032773147460627, +0.103575692245252, +0.296555596579887, +0.599868711174861},
-       Descriptor{0.015298306248441, +0.020083411655416, +0.337723063403079, +0.642193524941505},
-       Descriptor{0.002386244192839, -0.004341002614139, +0.204748281642812, +0.799592720971327},
-       Descriptor{0.019084792755899, +0.041941786468010, +0.189358492130623, +0.768699721401368},
-       Descriptor{0.006850054546542, +0.014317320230681, +0.085283615682657, +0.900399064086661}};
+      {Descriptor{1, 2.263228303690940e-02, {}},
+       Descriptor{2, 2.054646157184948e-02, {2.459900704671417e-01}},
+       Descriptor{2, 2.035591665621268e-02, {4.155848968854205e-01}},
+       Descriptor{2, 7.390817345112202e-03, {8.535556658670032e-02}},
+       Descriptor{2, 1.470920484949405e-02, {1.619186441912712e-01}},
+       Descriptor{2, 2.209273156075285e-03, {5.000000000000000e-01}},
+       Descriptor{2, 1.298716664913858e-02, {4.752807275459421e-01}},
+       Descriptor{3, 9.469136232207850e-03, {1.910747636405291e-01, 5.475517491470312e-02}},
+       Descriptor{3, 8.272333574175241e-04, {8.552204200227611e-03, 2.320342776881371e-02}},
+       Descriptor{3, 7.504300892142903e-03, {3.317645234741477e-01, 1.893177828040591e-02}},
+       Descriptor{3, 3.973796966696249e-03, {8.069616698587292e-02, 1.903012974369745e-02}},
+       Descriptor{3, 1.599180503968503e-02, {3.082449691963540e-01, 1.026061902393981e-01}},
+       Descriptor{3, 2.695593558424406e-03, {1.874417824837821e-01, 5.936350016822270e-03}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -361,26 +366,24 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 17: {
-    constexpr size_t nb_points = 61;
+    constexpr size_t nb_points = 60;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.033437199290803, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.005093415440507, 0.005658918886452, 0.497170540556774, 0.497170540556774},
-       Descriptor{0.014670864527638, 0.035647354750751, 0.482176322624625, 0.482176322624625},
-       Descriptor{0.024350878353672, 0.099520061958437, 0.450239969020782, 0.450239969020782},
-       Descriptor{0.031107550868969, 0.199467521245206, 0.400266239377397, 0.400266239377397},
-       Descriptor{0.031257111218620, 0.495717464058095, 0.252141267970953, 0.252141267970953},
-       Descriptor{0.024815654339665, 0.675905990683077, 0.162047004658461, 0.162047004658461},
-       Descriptor{0.014056073070557, 0.848248235478508, 0.075875882260746, 0.075875882260746},
-       Descriptor{0.003194676173779, 0.968690546064356, 0.015654726967822, 0.015654726967822},
-       Descriptor{0.008119655318993, 0.010186928826919, 0.334319867363658, 0.655493203809423},
-       Descriptor{0.026805742283163, 0.135440871671036, 0.292221537796944, 0.572337590532020},
-       Descriptor{0.018459993210822, 0.054423924290583, 0.319574885423190, 0.626001190286228},
-       Descriptor{0.008476868534328, 0.012868560833637, 0.190704224192292, 0.796427214974071},
-       Descriptor{0.018292796770025, 0.067165782413524, 0.180483211648746, 0.752351005937729},
-       Descriptor{0.006665632004165, 0.014663182224828, 0.080711313679564, 0.904625504095608}};
+      {Descriptor{2, 1.365546326405105e-02, {4.171034443615992e-01}},
+       Descriptor{2, 1.386943788818821e-03, {1.475549166075395e-02}},
+       Descriptor{2, 1.250972547524868e-02, {4.655978716188903e-01}},
+       Descriptor{2, 1.315631529400899e-02, {1.803581162663706e-01}},
+       Descriptor{2, 6.229500401152721e-03, {6.665406347959693e-02}},
+       Descriptor{2, 1.885811857639764e-02, {2.857065024365866e-01}},
+       Descriptor{3, 3.989150102964797e-03, {1.591922874727927e-01, 1.601764236211930e-02}},
+       Descriptor{3, 1.124388627334553e-02, {6.734937786736120e-02, 3.062815917461865e-01}},
+       Descriptor{3, 5.199219977919768e-03, {4.154754592952291e-01, 1.322967276008689e-02}},
+       Descriptor{3, 1.027894916022726e-02, {1.687225134952595e-01, 7.804234056828242e-02}},
+       Descriptor{3, 4.346107250500596e-03, {2.717918700553548e-01, 1.313587083400269e-02}},
+       Descriptor{3, 2.292174200867934e-03, {7.250547079900242e-02, 1.157517590318062e-02}},
+       Descriptor{3, 1.308581296766849e-02, {2.992189424769703e-01, 1.575054779268699e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -389,28 +392,26 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     break;
   }
   case 18: {
-    constexpr size_t nb_points = 70;
+    constexpr size_t nb_points = 67;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{+0.030809939937647, +0.333333333333333, +0.333333333333333, +0.333333333333333},
-       Descriptor{+0.009072436679404, +0.013310382738157, +0.493344808630921, +0.493344808630921},
-       Descriptor{+0.018761316939594, +0.061578811516086, +0.469210594241957, +0.469210594241957},
-       Descriptor{+0.019441097985477, +0.127437208225989, +0.436281395887006, +0.436281395887006},
-       Descriptor{+0.027753948610810, +0.210307658653168, +0.394846170673416, +0.394846170673416},
-       Descriptor{+0.032256225351457, +0.500410862393686, +0.249794568803157, +0.249794568803157},
-       Descriptor{+0.025074032616922, +0.677135612512315, +0.161432193743843, +0.161432193743843},
-       Descriptor{+0.015271927971832, +0.846803545029257, +0.076598227485371, +0.076598227485371},
-       Descriptor{+0.006793922022963, +0.951495121293100, +0.024252439353450, +0.024252439353450},
-       Descriptor{-0.002223098729920, +0.913707265566071, +0.043146367216965, +0.043146367216965},
-       Descriptor{+0.006331914076406, +0.008430536202420, +0.358911494940944, +0.632657968856636},
-       Descriptor{+0.027257538049138, +0.131186551737188, +0.294402476751957, +0.574410971510855},
-       Descriptor{+0.017676785649465, +0.050203151565675, +0.325017801641814, +0.624779046792512},
-       Descriptor{+0.018379484638070, +0.066329263810916, +0.184737559666046, +0.748933176523037},
-       Descriptor{+0.008104732808192, +0.011996194566236, +0.218796800013321, +0.769207005420443},
-       Descriptor{+0.007634129070725, +0.014858100590125, +0.101179597136408, +0.883962302273467},
-       Descriptor{+0.000046187660794, -0.035222015287949, +0.020874755282586, +1.014347260005363}};
+      {Descriptor{1, 1.817786765071333e-02, {}},
+       Descriptor{2, 1.665223501669507e-02, {3.999556280675762e-01}},
+       Descriptor{2, 6.023323816999855e-03, {4.875803015748695e-01}},
+       Descriptor{2, 9.474585753389433e-03, {4.618095064064492e-01}},
+       Descriptor{2, 1.823754470447182e-02, {2.422647025142720e-01}},
+       Descriptor{2, 3.564663009859485e-03, {3.883025608868559e-02}},
+       Descriptor{2, 8.279579976001624e-03, {9.194774212164319e-02}},
+       Descriptor{3, 6.879808117471103e-03, {1.838227079254640e-01, 4.580491585986078e-02}},
+       Descriptor{3, 1.189095545007642e-02, {1.226967573719275e-01, 2.063492574338379e-01}},
+       Descriptor{3, 2.265267251128533e-03, {3.956834343322697e-01, 3.897611033473383e-03}},
+       Descriptor{3, 3.420055059803591e-03, {1.081957937910333e-01, 1.346201674144499e-02}},
+       Descriptor{3, 8.873744551010202e-03, {3.197516245253774e-01, 4.026028346990806e-02}},
+       Descriptor{3, 2.505330437289861e-03, {2.357721849581917e-01, 5.298335186609765e-03}},
+       Descriptor{3, 6.114740634805449e-04, {2.709091099516201e-02, 5.483600420423190e-04}},
+       Descriptor{3, 1.274108765591222e-02, {3.334935294498808e-01, 1.205876951639246e-01}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -424,23 +425,23 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{0.032906331388919, 0.333333333333333, 0.333333333333333, 0.333333333333333},
-       Descriptor{0.010330731891272, 0.020780025853987, 0.489609987073006, 0.489609987073006},
-       Descriptor{0.022387247263016, 0.090926214604215, 0.454536892697893, 0.454536892697893},
-       Descriptor{0.030266125869468, 0.197166638701138, 0.401416680649431, 0.401416680649431},
-       Descriptor{0.030490967802198, 0.488896691193805, 0.255551654403098, 0.255551654403098},
-       Descriptor{0.024159212741641, 0.645844115695741, 0.177077942152130, 0.177077942152130},
-       Descriptor{0.016050803586801, 0.779877893544096, 0.110061053227952, 0.110061053227952},
-       Descriptor{0.008084580261784, 0.888942751496321, 0.055528624251840, 0.055528624251840},
-       Descriptor{0.002079362027485, 0.974756272445543, 0.012621863777229, 0.012621863777229},
-       Descriptor{0.003884876904981, 0.003611417848412, 0.395754787356943, 0.600633794794645},
-       Descriptor{0.025574160612022, 0.134466754530780, 0.307929983880436, 0.557603261588784},
-       Descriptor{0.008880903573338, 0.014446025776115, 0.264566948406520, 0.720987025817365},
-       Descriptor{0.016124546761731, 0.046933578838178, 0.358539352205951, 0.594527068955871},
-       Descriptor{0.002491941817491, 0.002861120350567, 0.157807405968595, 0.839331473680839},
-       Descriptor{0.018242840118951, 0.223861424097916, 0.075050596975911, 0.701087978926173},
-       Descriptor{0.010258563736199, 0.034647074816760, 0.142421601113383, 0.822931324069857},
-       Descriptor{0.003799928855302, 0.010161119296278, 0.065494628082938, 0.924344252620784}};
+      {Descriptor{1, 1.723469885200617e-02, {}},
+       Descriptor{2, 3.554628298899065e-03, {5.252389035120897e-02}},
+       Descriptor{2, 5.160877571472141e-03, {4.925126750413369e-01}},
+       Descriptor{2, 7.617175546509150e-03, {1.114488733230214e-01}},
+       Descriptor{2, 1.149179501337080e-02, {4.591942010395437e-01}},
+       Descriptor{2, 1.576876744657749e-02, {4.039697225519012e-01}},
+       Descriptor{2, 1.232595742409543e-02, {1.781701047817643e-01}},
+       Descriptor{2, 8.826613882214238e-04, {1.163946118378945e-02}},
+       Descriptor{2, 1.587650968300154e-02, {2.551616329136077e-01}},
+       Descriptor{3, 4.847742243427523e-03, {3.914585933169222e-02, 1.306976762680324e-01}},
+       Descriptor{3, 1.317316098869537e-02, {1.293125644701578e-01, 3.113176298095413e-01}},
+       Descriptor{3, 1.641038275917910e-03, {3.646177809746111e-01, 2.068925896604807e-03}},
+       Descriptor{3, 9.053972465606226e-03, {2.214348854323312e-01, 7.456029460162668e-02}},
+       Descriptor{3, 1.463157551735100e-03, {1.424257573657563e-01, 5.007288257354491e-03}},
+       Descriptor{3, 8.051081382012054e-03, {3.540280097352752e-01, 4.088801119601688e-02}},
+       Descriptor{3, 4.227943749768248e-03, {1.492405208198407e-02, 2.418945789605796e-01}},
+       Descriptor{3, 1.663600681429694e-03, {9.776025800888155e-03, 6.008627532230670e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
@@ -454,25 +455,24 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<double> weight_list(nb_points);
 
     std::array descriptor_list =   //
-      {Descriptor{+0.033057055541624, +0.333333333333333, +0.333333333333333, +0.333333333333333},
-       Descriptor{+0.000867019185663, -0.001900928704400, +0.500950464352200, +0.500950464352200},
-       Descriptor{+0.011660052716448, +0.023574084130543, +0.488212957934729, +0.488212957934729},
-       Descriptor{+0.022876936356421, +0.089726636099435, +0.455136681950283, +0.455136681950283},
-       Descriptor{+0.030448982673938, +0.196007481363421, +0.401996259318289, +0.401996259318289},
-       Descriptor{+0.030624891725355, +0.488214180481157, +0.255892909759421, +0.255892909759421},
-       Descriptor{+0.024368057676800, +0.647023488009788, +0.176488255995106, +0.176488255995106},
-       Descriptor{+0.015997432032024, +0.791658289326483, +0.104170855336758, +0.104170855336758},
-       Descriptor{+0.007698301815602, +0.893862072318140, +0.053068963840930, +0.053068963840930},
-       Descriptor{-0.000632060497488, +0.916762569607942, +0.041618715196029, +0.041618715196029},
-       Descriptor{+0.001751134301193, +0.976836157186356, +0.011581921406822, +0.011581921406822},
-       Descriptor{+0.016465839189576, +0.048741583664839, +0.344855770229001, +0.606402646106160},
-       Descriptor{+0.004839033540485, +0.006314115948605, +0.377843269594854, +0.615842614456541},
-       Descriptor{+0.025804906534650, +0.134316520547348, +0.306635479062357, +0.559048000390295},
-       Descriptor{+0.008471091054441, +0.013973893962392, +0.249419362774742, +0.736606743262866},
-       Descriptor{+0.018354914106280, +0.075549132909764, +0.212775724802802, +0.711675142287434},
-       Descriptor{+0.000704404677908, -0.008368153208227, +0.146965436053239, +0.861402717154987},
-       Descriptor{+0.010112684927462, +0.026686063258714, +0.137726978828923, +0.835586957912363},
-       Descriptor{+0.003573909385950, +0.010547719294141, +0.059696109149007, +0.929756171556853}};
+      {Descriptor{1, 1.391011070145312e-02, {}},
+       Descriptor{2, 1.408320130752025e-02, {2.545792676733391e-01}},
+       Descriptor{2, 7.988407910666199e-04, {1.097614102839776e-02}},
+       Descriptor{2, 7.830230776074533e-03, {1.093835967117146e-01}},
+       Descriptor{2, 9.173462974252915e-03, {1.862949977445409e-01}},
+       Descriptor{2, 9.452399933232448e-03, {4.455510569559248e-01}},
+       Descriptor{2, 2.161275410665577e-03, {3.731088059888470e-02}},
+       Descriptor{2, 1.378805062907046e-02, {3.934253478170999e-01}},
+       Descriptor{2, 7.101825303408441e-03, {4.762456115404990e-01}},
+       Descriptor{3, 2.202897418558497e-03, {1.591337076570672e-01, 7.570780504696529e-03}},
+       Descriptor{3, 5.986398578954690e-03, {1.985181322287882e-01, 4.656036490766432e-02}},
+       Descriptor{3, 1.129869602125866e-03, {4.854937607623754e-03, 6.409058560843406e-02}},
+       Descriptor{3, 8.667225567219333e-03, {3.331348173095875e-01, 5.498747914298681e-02}},
+       Descriptor{3, 4.145711527613858e-03, {3.836368477537459e-02, 9.995229628813866e-02}},
+       Descriptor{3, 7.722607822099230e-03, {2.156070573900944e-01, 1.062272047202700e-01}},
+       Descriptor{3, 3.695681500255298e-03, {9.831548292802561e-03, 4.200237588162241e-01}},
+       Descriptor{3, 1.169174573182774e-02, {1.398080719917999e-01, 3.178601238357720e-01}},
+       Descriptor{3, 3.578200238457685e-03, {2.805814114236652e-01, 1.073721285601109e-02}}};
 
     fill_quadrature_points(descriptor_list, point_list, weight_list);
 
diff --git a/src/analysis/TriangleGaussQuadrature.hpp b/src/analysis/TriangleGaussQuadrature.hpp
index bf11f9ca3e7c5ce0f6805e982563365a9229d9d1..58b0be03e2efb7b23e6247b2564e238ac1f01734 100644
--- a/src/analysis/TriangleGaussQuadrature.hpp
+++ b/src/analysis/TriangleGaussQuadrature.hpp
@@ -6,7 +6,10 @@
 /**
  * Defines Gauss quadrature on the P1 reference triangle
  *
- * formulae are extracted from Dunavant (https://doi.org/10.1002/nme.1620210612)
+ * \note formulae are provided by
+ *
+ * ‘On the identification of symmetric quadrature rules for finite
+ * element methods‘ by F.D. Witherden & P.E. Vincent (2015).
  */
 class TriangleGaussQuadrature final : public QuadratureForumla<2>
 {
diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp
index 642ce4ae428c8da74ff786b7a291ad91741d69a8..a02a8f7285f162c416c38e9dcb80fb54a730884d 100644
--- a/tests/test_TriangleGaussQuadrature.cpp
+++ b/tests/test_TriangleGaussQuadrature.cpp
@@ -191,24 +191,19 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
     REQUIRE(get_order(p3, l2, 146. / 5) == Catch::Approx(4));
   }
 
-  SECTION("degree 3")
+  SECTION("degree 3 and 4")
   {
     TriangleGaussQuadrature l3(3);
+    TriangleGaussQuadrature l4(4);
 
-    REQUIRE(l3.numberOfPoints() == 4);
-
-    REQUIRE(integrate(p0, l3) == Catch::Approx(4));
-    REQUIRE(integrate(p1, l3) == Catch::Approx(-4. / 3));
-    REQUIRE(integrate(p2, l3) == Catch::Approx(-19. / 3));
-    REQUIRE(integrate(p3, l3) == Catch::Approx(146. / 5));
-    REQUIRE(integrate(p4, l3) != Catch::Approx(-25. / 6));
-
-    REQUIRE(get_order(p4, l3, -25. / 6) == Catch::Approx(4));
-  }
+    bool is_same = true;
+    is_same &= (l3.numberOfPoints() == l4.numberOfPoints());
+    for (size_t i = 0; i < l3.pointList().size(); ++i) {
+      is_same &= (l3.point(i) == l3.point(i));
+      is_same &= (l3.weight(i) == l3.weight(i));
+    }
 
-  SECTION("degree 4")
-  {
-    TriangleGaussQuadrature l4(4);
+    REQUIRE(is_same);
 
     REQUIRE(l4.numberOfPoints() == 6);
 
@@ -262,7 +257,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l7(7);
 
-    REQUIRE(l7.numberOfPoints() == 13);
+    REQUIRE(l7.numberOfPoints() == 15);
 
     REQUIRE(integrate(p0, l7) == Catch::Approx(4));
     REQUIRE(integrate(p1, l7) == Catch::Approx(-4. / 3));
@@ -346,7 +341,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l11(11);
 
-    REQUIRE(l11.numberOfPoints() == 27);
+    REQUIRE(l11.numberOfPoints() == 28);
 
     REQUIRE(integrate(p0, l11) == Catch::Approx(4));
     REQUIRE(integrate(p1, l11) == Catch::Approx(-4. / 3));
@@ -446,7 +441,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l15(15);
 
-    REQUIRE(l15.numberOfPoints() == 48);
+    REQUIRE(l15.numberOfPoints() == 49);
 
     REQUIRE(integrate(p0, l15) == Catch::Approx(4));
     REQUIRE(integrate(p1, l15) == Catch::Approx(-4. / 3));
@@ -473,7 +468,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l16(16);
 
-    REQUIRE(l16.numberOfPoints() == 52);
+    REQUIRE(l16.numberOfPoints() == 55);
 
     REQUIRE(integrate(p0, l16) == Catch::Approx(4));
     REQUIRE(integrate(p1, l16) == Catch::Approx(-4. / 3));
@@ -502,7 +497,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l17(17);
 
-    REQUIRE(l17.numberOfPoints() == 61);
+    REQUIRE(l17.numberOfPoints() == 60);
 
     REQUIRE(integrate(p0, l17) == Catch::Approx(4));
     REQUIRE(integrate(p1, l17) == Catch::Approx(-4. / 3));
@@ -531,7 +526,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
   {
     TriangleGaussQuadrature l18(18);
 
-    REQUIRE(l18.numberOfPoints() == 70);
+    REQUIRE(l18.numberOfPoints() == 67);
 
     REQUIRE(integrate(p0, l18) == Catch::Approx(4));
     REQUIRE(integrate(p1, l18) == Catch::Approx(-4. / 3));