diff --git a/src/analysis/SquareGaussQuadrature.cpp b/src/analysis/SquareGaussQuadrature.cpp
index b8762716f9c012e4961aeff62df04a0a261aaffe..99cbc79608140a5564b083f33c71f92faa64de13 100644
--- a/src/analysis/SquareGaussQuadrature.cpp
+++ b/src/analysis/SquareGaussQuadrature.cpp
@@ -4,63 +4,88 @@
 void
 SquareGaussQuadrature::_buildPointAndWeightLists()
 {
-  using R2              = TinyVector<2>;
-  using QuadratureGroup = std::array<double, 3>;
+  using R2 = TinyVector<2>;
 
-  auto fill_quadrature_points = [](auto quadrature_group_list, auto& point_list, auto& weight_list) {
+  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());
 
     size_t k = 0;
-    for (size_t i = 0; i < quadrature_group_list.size(); ++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) {
-          point_list[k + 0] = {xi, eta};
-          point_list[k + 1] = {-xi, eta};
-          point_list[k + 2] = {eta, xi};
-          point_list[k + 3] = {eta, -xi};
-
-          for (size_t i = 0; i < 4; ++i) {
-            weight_list[k + i] = w;
-          }
-
-          k += 4;
-        } else {
-          point_list[k + 0] = {xi, eta};
-          point_list[k + 1] = {-xi, eta};
-          point_list[k + 2] = {xi, -eta};
-          point_list[k + 3] = {-xi, -eta};
-          point_list[k + 4] = {eta, xi};
-          point_list[k + 5] = {-eta, xi};
-          point_list[k + 6] = {eta, -xi};
-          point_list[k + 7] = {-eta, -xi};
-
-          for (size_t i = 0; i < 8; ++i) {
-            weight_list[k + i] = w;
-          }
-
-          k += 8;
+    for (size_t i = 0; i < descriptor_list.size(); ++i) {
+      const auto [id, w, position_list] = descriptor_list[i];
+
+      switch (id) {
+      case 1: {
+        Assert(position_list.size() == 0);
+
+        point_list[k]  = {0, 0};
+        weight_list[k] = w;
+
+        k += 1;
+        break;
+      }
+      case 2: {
+        Assert(position_list.size() == 1);
+        const double& a = position_list[0];
+
+        point_list[k + 0] = {+a, 0};
+        point_list[k + 1] = {-a, 0};
+        point_list[k + 2] = {0, +a};
+        point_list[k + 3] = {0, -a};
+
+        for (size_t l = 0; l < 4; ++l) {
+          weight_list[k + l] = w;
         }
-      } else {
-        if (xi == 0) {
-          point_list[k] = {xi, eta};
 
-          weight_list[k] = w;
+        k += 4;
+        break;
+      }
+      case 3: {
+        Assert(position_list.size() == 1);
+        const double& a = position_list[0];
 
-          k += 1;
-        } else {
-          point_list[k + 0] = {xi, eta};
-          point_list[k + 1] = {-xi, eta};
-          point_list[k + 2] = {xi, -eta};
-          point_list[k + 3] = {-xi, -eta};
+        point_list[k + 0] = {+a, +a};
+        point_list[k + 1] = {+a, -a};
+        point_list[k + 2] = {-a, +a};
+        point_list[k + 3] = {-a, -a};
 
-          for (size_t i = 0; i < 4; ++i) {
-            weight_list[k + i] = w;
-          }
+        for (size_t l = 0; l < 4; ++l) {
+          weight_list[k + l] = w;
+        }
 
-          k += 4;
+        k += 4;
+        break;
+      }
+      case 4: {
+        Assert(position_list.size() == 2);
+        const double& a = position_list[0];
+        const double& b = position_list[1];
+
+        point_list[k + 0] = {+a, +b};
+        point_list[k + 1] = {+a, -b};
+        point_list[k + 2] = {-a, +b};
+        point_list[k + 3] = {-a, -b};
+        point_list[k + 4] = {+b, +a};
+        point_list[k + 5] = {+b, -a};
+        point_list[k + 6] = {-b, +a};
+        point_list[k + 7] = {-b, -a};
+
+        for (size_t l = 0; l < 8; ++l) {
+          weight_list[k + l] = w;
         }
+
+        k += 8;
+        break;
+      }
+      default: {
+        throw UnexpectedError("invalid quadrature id");
+      }
       }
     }
 
@@ -74,10 +99,10 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array<QuadratureGroup, 1> quadrature_group_list =   //
-      {QuadratureGroup{4.000000000000000, 0.000000000000000, 0.000000000000000}};
+    std::array descriptor_list =   //
+      {Descriptor{1, 4.000000000000000e+00, {}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -89,10 +114,10 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array<QuadratureGroup, 1> quadrature_group_list =   //
-      {QuadratureGroup{1.000000000000000, 0.577350269189626, 0.577350269189626}};
+    std::array descriptor_list =   //
+      {Descriptor{3, 1.000000000000000e+00, {5.773502691896257e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -104,11 +129,11 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.816326530612245, 0.683130051063973, 0.000000000000000},
-       QuadratureGroup{0.183673469387755, 0.881917103688197, 0.881917103688197}};
+    std::array descriptor_list =   //
+      {Descriptor{2, 8.163265306122449e-01, {6.831300510639732e-01}},
+       Descriptor{3, 1.836734693877551e-01, {8.819171036881969e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -120,12 +145,12 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.241975308641975, 0.925820099772551, 0.000000000000000},
-       QuadratureGroup{0.237431774690630, 0.805979782918599, 0.805979782918599},
-       QuadratureGroup{0.520592916667394, 0.380554433208316, 0.380554433208316}};
+    std::array descriptor_list =   //
+      {Descriptor{2, 2.419753086419753e-01, {9.258200997725514e-01}},
+       Descriptor{3, 2.374317746906302e-01, {8.059797829185987e-01}},
+       Descriptor{3, 5.205929166673945e-01, {3.805544332083157e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -137,13 +162,13 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.018475842507491, 1.121225763866564, 0.000000000000000},
-       QuadratureGroup{0.390052939160735, 0.451773049920657, 0.000000000000000},
-       QuadratureGroup{0.083095178026482, 0.891849420851512, 0.891849420851512},
-       QuadratureGroup{0.254188020152646, 0.824396370149276, 0.411623426336542}};
+    std::array descriptor_list =   //
+      {Descriptor{2, 4.541639606867490e-01, {4.889268569743691e-01}},
+       Descriptor{3, 4.273123186577577e-02, {9.396552580968377e-01}},
+       Descriptor{3, 2.142003609268616e-01, {6.908805504863439e-01}},
+       Descriptor{4, 1.444522232603068e-01, {9.186204410567222e-01, 3.448720253644036e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -151,19 +176,18 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
   }
   case 10:
   case 11: {
-    constexpr size_t nb_points = 25;
+    constexpr size_t nb_points = 28;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.365379525585903, 0.000000000000000, 0.000000000000000},
-       QuadratureGroup{0.027756165564204, 1.044402915409813, 0.000000000000000},
-       QuadratureGroup{0.244272057751754, 0.769799068396649, 0.000000000000000},
-       QuadratureGroup{0.034265103851229, 0.935787012440540, 0.935787012440540},
-       QuadratureGroup{0.308993036133713, 0.413491953449114, 0.413491953449114},
-       QuadratureGroup{0.146684377651312, 0.883025508525690, 0.575653595840465}};
+    std::array descriptor_list =   //
+      {Descriptor{2, 2.174004398687120e-01, {7.146178296646060e-01}},
+       Descriptor{3, 2.772741029838511e-01, {2.736572101714596e-01}},
+       Descriptor{3, 2.139336378782481e-01, {6.366039322123010e-01}},
+       Descriptor{4, 4.407456911498309e-02, {9.516303887840335e-01, 8.155654336896384e-01}},
+       Descriptor{4, 1.016213405196113e-01, {3.462072000476454e-01, 9.355678714875911e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -171,20 +195,21 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
   }
   case 12:
   case 13: {
-    constexpr size_t nb_points = 36;
+    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.005656169693764, 1.086056158573971, 0.000000000000000},
-       QuadratureGroup{0.192443867470396, 0.658208197042585, 0.000000000000000},
-       QuadratureGroup{0.005166832979773, 1.001300602991729, 1.001300602991729},
-       QuadratureGroup{0.200302559622138, 0.584636168775946, 0.584636168775946},
-       QuadratureGroup{0.228125175912536, 0.246795612720261, 0.246795612720261},
-       QuadratureGroup{0.117496926974491, 0.900258815287201, 0.304720678579870},
-       QuadratureGroup{0.066655770186205, 0.929866705560780, 0.745052720131169}};
+    std::array descriptor_list =   //
+      {Descriptor{1, 2.999933443589096e-01, {}},
+       Descriptor{2, 3.812862534985274e-02, {9.834613326132455e-01}},
+       Descriptor{2, 1.845354689698072e-01, {6.398614183671097e-01}},
+       Descriptor{3, 3.950714064745244e-02, {9.187784880797114e-01}},
+       Descriptor{3, 2.313985194006854e-01, {3.796285051867438e-01}},
+       Descriptor{3, 1.372420967130035e-01, {6.995542165133511e-01}},
+       Descriptor{4, 3.351978005038143e-02, {6.435973749966181e-01, 9.750688839059836e-01}},
+       Descriptor{4, 1.135751263643542e-01, {3.335398811647831e-01, 8.644276092670619e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -192,22 +217,21 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
   }
   case 14:
   case 15: {
-    constexpr size_t nb_points = 45;
+    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.001768979827207, 0.000000000000000, 0.000000000000000},
-       QuadratureGroup{+0.012816726617512, 1.027314357719367, 0.000000000000000},
-       QuadratureGroup{+0.119897873101347, 0.856766776147643, 0.000000000000000},
-       QuadratureGroup{+0.210885452208801, 0.327332998189723, 0.000000000000000},
-       QuadratureGroup{+0.006392720128215, 0.967223740028505, 0.967223740028505},
-       QuadratureGroup{+0.104415680788580, 0.732168901749711, 0.732168901749711},
-       QuadratureGroup{+0.168053047203816, 0.621974427996805, 0.321696694921009},
-       QuadratureGroup{+0.076169694452294, 0.928618480068352, 0.455124178121179},
-       QuadratureGroup{+0.028794154400064, 0.960457474887516, 0.809863684081217}};
+    std::array descriptor_list =   //
+      {Descriptor{2, 1.149329272635254e-01, {7.980989878970071e-01}},
+       Descriptor{2, 1.816857589672718e-01, {3.038530263984597e-01}},
+       Descriptor{3, 4.123848378876581e-02, {8.824222701068853e-01}},
+       Descriptor{3, 5.933746485839221e-03, {9.777897995399027e-01}},
+       Descriptor{4, 1.011491434774324e-01, {8.087121358195035e-01, 5.672135733965907e-01}},
+       Descriptor{4, 1.478367216328812e-01, {3.046547990370526e-01, 5.787361940358066e-01}},
+       Descriptor{4, 2.327319414467321e-02, {9.805048437245319e-01, 6.974636319909671e-01}},
+       Descriptor{4, 5.584548249231202e-02, {2.648441558723162e-01, 9.557425198095117e-01}}};
 
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -219,19 +243,20 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.020614915919991, 0.989353074512600, 0.000000000000000},
-       QuadratureGroup{0.128025716179910, 0.376285207157973, 0.000000000000000},
-       QuadratureGroup{0.005511739534032, 0.978848279262233, 0.978848279262233},
-       QuadratureGroup{0.039207712457142, 0.885794729164116, 0.885794729164116},
-       QuadratureGroup{0.076396945079863, 0.171756123838348, 0.171756123838348},
-       QuadratureGroup{0.141513729949972, 0.590499273806002, 0.319505036634574},
-       QuadratureGroup{0.083903279363798, 0.799079131916863, 0.597972451929457},
-       QuadratureGroup{0.060394163649685, 0.803743962958745, 0.058344481776551},
-       QuadratureGroup{0.057387752969213, 0.936506276127495, 0.347386316166203},
-       QuadratureGroup{0.021922559481864, 0.981321179805452, 0.706000287798646}};
-
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    std::array descriptor_list =   //
+      {Descriptor{2, 1.388043776872199e-01, {5.450429872091483e-01}},
+       Descriptor{2, 6.534782380746491e-02, {9.174971882191532e-01}},
+       Descriptor{3, 1.428564196221996e-01, {1.925542426140469e-01}},
+       Descriptor{3, 3.671590848525151e-02, {8.713286846954432e-01}},
+       Descriptor{3, 8.625569963814381e-02, {6.943880900895507e-01}},
+       Descriptor{3, 1.345562732119169e-01, {4.576829928831555e-01}},
+       Descriptor{3, 7.655633414909936e-03, {9.692042560915087e-01}},
+       Descriptor{4, 1.031767247849354e-01, {7.620882760801982e-01, 2.794777667151921e-01}},
+       Descriptor{4, 5.482640023729144e-02, {9.073943521982236e-01, 5.468986704143101e-01}},
+       Descriptor{4, 1.511981038825686e-02, {7.612486664390827e-01, 9.837879715015495e-01}},
+       Descriptor{4, 2.078099665596304e-02, {9.870947070447679e-01, 3.028074817888614e-01}}};
+
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -243,21 +268,21 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.038205406871462, 0.943962831808239, 0.000000000000000},
-       QuadratureGroup{0.135368502976521, 0.536918434376013, 0.000000000000000},
-       QuadratureGroup{0.005773708558664, 0.973981076394170, 0.973981076394170},
-       QuadratureGroup{0.067460759759473, 0.742995535327609, 0.742995535327609},
-       QuadratureGroup{0.140899115227892, 0.285010052188916, 0.285010052188916},
-       QuadratureGroup{0.047466627685662, 0.068354569272491, 0.068354569272491},
-       QuadratureGroup{0.078619467342982, 0.802952004398543, 0.203345534163332},
-       QuadratureGroup{0.094979169511394, 0.634244672807882, 0.426572172992877},
-       QuadratureGroup{0.022331162356015, 0.978350706908227, 0.295830776620995},
-       QuadratureGroup{0.055594877793785, 0.901672714410389, 0.541983037327871},
-       QuadratureGroup{0.006049054506376, 1.007018449383116, 0.669414798783936},
-       QuadratureGroup{0.024839207949609, 0.945161453573471, 0.829501421477824}};
-
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    std::array descriptor_list =   //
+      {Descriptor{2, 9.773736469282875e-02, {7.143442595727942e-01}},
+       Descriptor{2, 1.393076129224823e-01, {2.656720521209637e-01}},
+       Descriptor{2, 3.486958349188791e-02, {9.644342692579673e-01}},
+       Descriptor{3, 4.780454716279579e-02, {8.033797294676850e-01}},
+       Descriptor{3, 1.617715177911761e-03, {9.921654058971348e-01}},
+       Descriptor{3, 1.744104803435576e-02, {9.294496027996094e-01}},
+       Descriptor{4, 1.177258328400561e-01, {5.102782573693434e-01, 2.666403145945622e-01}},
+       Descriptor{4, 1.617957761165785e-02, {3.907342057752498e-01, 9.878929331417353e-01}},
+       Descriptor{4, 8.284661898340490e-02, {7.171679213097452e-01, 5.124918772160977e-01}},
+       Descriptor{4, 6.498908149259605e-02, {2.654400078112960e-01, 8.689024341545042e-01}},
+       Descriptor{4, 3.898055239641054e-02, {6.200353986932564e-01, 9.263029558071293e-01}},
+       Descriptor{4, 9.889400934743484e-03, {8.016715847185969e-01, 9.884465306839737e-01}}};
+
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
@@ -265,27 +290,28 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
   }
   case 20:
   case 21: {
-    constexpr size_t nb_points = 88;
+    constexpr size_t nb_points = 85;
     SmallArray<R2> point_list(nb_points);
     SmallArray<double> weight_list(nb_points);
 
-    std::array quadrature_group_list =   //
-      {QuadratureGroup{0.019503841092684, 0.980883148832881, 0.000000000000000},
-       QuadratureGroup{0.089012127744268, 0.678152700336576, 0.000000000000000},
-       QuadratureGroup{0.114568584702749, 0.240599282275864, 0.000000000000000},
-       QuadratureGroup{0.007463627359106, 0.965176994929162, 0.965176994929162},
-       QuadratureGroup{0.050585943594705, 0.749698539312765, 0.749698539312765},
-       QuadratureGroup{0.074613865184212, 0.568983925500818, 0.568983925500818},
-       QuadratureGroup{0.023501091310143, 0.971086142843168, 0.355832132274584},
-       QuadratureGroup{0.011588562644144, 0.983453947854968, 0.645588139196562},
-       QuadratureGroup{0.023073245798171, 0.933927707027213, 0.821920249234369},
-       QuadratureGroup{0.001570221774472, 1.014086498915039, 0.862185099566557},
-       QuadratureGroup{0.049102258016277, 0.877914842155496, 0.168914072450263},
-       QuadratureGroup{0.042512352239126, 0.882246882640128, 0.568113580166780},
-       QuadratureGroup{0.067270936863160, 0.741324453314596, 0.371360260002223},
-       QuadratureGroup{0.103507336515645, 0.469570217710647, 0.237333359193547}};
-
-    fill_quadrature_points(quadrature_group_list, point_list, weight_list);
+    std::array descriptor_list =   //
+      {Descriptor{1, 1.351928403561428e-01, {}},
+       Descriptor{2, 1.067941587859404e-01, {4.733510743582242e-01}},
+       Descriptor{2, 6.625720800494439e-02, {8.352784297558683e-01}},
+       Descriptor{3, 1.198844836463405e-01, {2.573719006072290e-01}},
+       Descriptor{3, 8.229461289220009e-03, {9.633378321156234e-01}},
+       Descriptor{3, 3.060364092877565e-02, {8.624519253796515e-01}},
+       Descriptor{3, 9.679179189359521e-02, {4.968979625193457e-01}},
+       Descriptor{3, 6.151634148131656e-02, {7.043321751954006e-01}},
+       Descriptor{4, 8.672849951320823e-02, {2.418781854767020e-01, 6.741462199553178e-01}},
+       Descriptor{4, 5.587911740412735e-02, {4.801569663127951e-01, 8.246473752709207e-01}},
+       Descriptor{4, 3.712882744091382e-02, {9.322419359217540e-01, 2.706991841016649e-01}},
+       Descriptor{4, 2.877560085102053e-02, {9.349023258240106e-01, 6.750395674370753e-01}},
+       Descriptor{4, 1.150952044249317e-02, {9.904554675295242e-01, 4.889162511771669e-01}},
+       Descriptor{4, 7.528482130401646e-03, {8.311352459759014e-01, 9.889606113865724e-01}},
+       Descriptor{4, 1.051230415825108e-02, {9.146960092229130e-02, 9.840171484895990e-01}}};
+
+    fill_quadrature_points(descriptor_list, point_list, weight_list);
 
     m_point_list  = point_list;
     m_weight_list = weight_list;
diff --git a/src/analysis/SquareGaussQuadrature.hpp b/src/analysis/SquareGaussQuadrature.hpp
index c699bcd214f1a05cd0a516676258ee4c0a0f221f..8dd1a00225658c0088c677702458e1353d903ad8 100644
--- a/src/analysis/SquareGaussQuadrature.hpp
+++ b/src/analysis/SquareGaussQuadrature.hpp
@@ -9,8 +9,8 @@
  *
  * \note formulae are provided by
  *
- * ‘Square symmetrical quadrature rules for complete polynomials over
- * a square domain‘ D. A. Dunavant (1985).
+ * ‘On the identification of symmetric quadrature rules for finite
+ * element methods‘ by F.D. Witherden & P.E. Vincent (2015).
  */
 class SquareGaussQuadrature final : public QuadratureForumla<2>
 {
diff --git a/tests/test_SquareGaussQuadrature.cpp b/tests/test_SquareGaussQuadrature.cpp
index 75146af1f3936c1e7df43b05eacba5dd5978fb4c..ee6d8d6d6346f647f7361ae6b7d4533677f812b0 100644
--- a/tests/test_SquareGaussQuadrature.cpp
+++ b/tests/test_SquareGaussQuadrature.cpp
@@ -297,7 +297,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
 
     REQUIRE(same_formulas(l10, l11));
 
-    REQUIRE(l11.numberOfPoints() == 25);
+    REQUIRE(l11.numberOfPoints() == 28);
 
     REQUIRE(integrate(p0, l11) == Catch::Approx(8));
     REQUIRE(integrate(p1, l11) == Catch::Approx(4));
@@ -323,7 +323,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
 
     REQUIRE(same_formulas(l12, l13));
 
-    REQUIRE(l13.numberOfPoints() == 36);
+    REQUIRE(l13.numberOfPoints() == 37);
 
     REQUIRE(integrate(p0, l13) == Catch::Approx(8));
     REQUIRE(integrate(p1, l13) == Catch::Approx(4));
@@ -351,7 +351,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
 
     REQUIRE(same_formulas(l14, l15));
 
-    REQUIRE(l15.numberOfPoints() == 45);
+    REQUIRE(l15.numberOfPoints() == 48);
 
     REQUIRE(integrate(p0, l15) == Catch::Approx(8));
     REQUIRE(integrate(p1, l15) == Catch::Approx(4));
@@ -435,7 +435,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
     REQUIRE(integrate(p17, l19) == Catch::Approx(16706355156097391. / 12085937500000));
     REQUIRE(integrate(p18, l19) == Catch::Approx(495866230514635109957. / 397838847656250000.));
     REQUIRE(integrate(p19, l19) == Catch::Approx(703712939580204375319. / 132612949218750000.));
-    REQUIRE(integrate(p20, l19) != Catch::Approx(-891851528496270127477. / 884086328125000000.));
+    REQUIRE(integrate(p20, l19) != Catch::Approx(-891851528496270127477. / 884086328125000000.).epsilon(1E-8));
 
     REQUIRE(get_order(p20, l19, -891851528496270127477. / 884086328125000000.) == Catch::Approx(20).margin(0.001));
   }
@@ -447,7 +447,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
 
     REQUIRE(same_formulas(l20, l21));
 
-    REQUIRE(l21.numberOfPoints() == 88);
+    REQUIRE(l21.numberOfPoints() == 85);
 
     REQUIRE(integrate(p0, l21) == Catch::Approx(8));
     REQUIRE(integrate(p1, l21) == Catch::Approx(4));