diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index 1760cf8565e9c398fd3d5275d7049ae413394c55..c818e70b3ee6e10c68b4a4c0d7000225e594ac0e 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(
PugsAnalysis
CubeGaussQuadrature.cpp
+ PrismGaussQuadrature.cpp
PyramidGaussQuadrature.cpp
SquareGaussQuadrature.cpp
TensorialGaussLegendreQuadrature.cpp
diff --git a/src/analysis/PrismGaussQuadrature.cpp b/src/analysis/PrismGaussQuadrature.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a5ee6b8902a40899b7aaa488cdd67ccae02b960
--- /dev/null
+++ b/src/analysis/PrismGaussQuadrature.cpp
@@ -0,0 +1,801 @@
+#include <analysis/PrismGaussQuadrature.hpp>
+#include <utils/Exceptions.hpp>
+
+void
+PrismGaussQuadrature::_buildPointAndWeightLists()
+{
+ using R2 = TinyVector<2>;
+ using R3 = TinyVector<3>;
+
+ struct Descriptor
+ {
+ int id;
+ double weight;
+ std::vector<double> lambda_list;
+ };
+
+ auto fill_quadrature_points = [](auto descriptor_list, auto& point_list, auto& weight_list) {
+ auto to_R3 = [](R2 X, double z) { return R3{X[0], X[1], z}; };
+
+ Assert(point_list.size() == weight_list.size());
+
+ const R2 A = {0, 0};
+ const R2 B = {1, 0};
+ const R2 C = {0, 1};
+
+ size_t k = 0;
+ for (size_t i = 0; i < descriptor_list.size(); ++i) {
+ const auto [id, w, value_list] = descriptor_list[i];
+
+ switch (id) {
+ case 1: {
+ Assert(value_list.size() == 0);
+
+ point_list[k] = {1. / 3, 1. / 3, 0.};
+ weight_list[k] = w;
+
+ k += 1;
+ break;
+ }
+ case 2: {
+ Assert(value_list.size() == 1);
+ const double c = value_list[0];
+ point_list[k + 0] = {1. / 3, 1. / 3, +c};
+ point_list[k + 1] = {1. / 3, 1. / 3, -c};
+
+ for (size_t l = 0; l < 2; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 2;
+ break;
+ }
+ case 3: {
+ Assert(value_list.size() == 1);
+ const double l0 = value_list[0];
+ const double l1 = 1 - 2 * l0;
+
+ point_list[k + 0] = to_R3(l0 * A + l0 * B + l1 * C, 0);
+ point_list[k + 1] = to_R3(l0 * A + l1 * B + l0 * C, 0);
+ point_list[k + 2] = to_R3(l1 * A + l0 * B + l0 * C, 0);
+
+ for (size_t l = 0; l < 3; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 3;
+ break;
+ }
+ case 4: {
+ Assert(value_list.size() == 2);
+ const double l0 = value_list[0];
+ const double l1 = 1 - 2 * l0;
+ const double z = value_list[1];
+
+ point_list[k + 0] = to_R3(l0 * A + l0 * B + l1 * C, +z);
+ point_list[k + 1] = to_R3(l0 * A + l1 * B + l0 * C, +z);
+ point_list[k + 2] = to_R3(l1 * A + l0 * B + l0 * C, +z);
+ point_list[k + 3] = to_R3(l0 * A + l0 * B + l1 * C, -z);
+ point_list[k + 4] = to_R3(l0 * A + l1 * B + l0 * C, -z);
+ point_list[k + 5] = to_R3(l1 * A + l0 * B + l0 * C, -z);
+
+ for (size_t l = 0; l < 6; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 6;
+ break;
+ }
+ case 5: {
+ Assert(value_list.size() == 2);
+ const double l0 = value_list[0];
+ const double l1 = value_list[1];
+ const double l2 = 1 - l0 - l1;
+
+ point_list[k + 0] = to_R3(l0 * A + l1 * B + l2 * C, 0);
+ point_list[k + 1] = to_R3(l0 * A + l2 * B + l1 * C, 0);
+ point_list[k + 2] = to_R3(l1 * A + l0 * B + l2 * C, 0);
+ point_list[k + 3] = to_R3(l1 * A + l2 * B + l0 * C, 0);
+ point_list[k + 4] = to_R3(l2 * A + l0 * B + l1 * C, 0);
+ point_list[k + 5] = to_R3(l2 * A + l1 * B + l0 * C, 0);
+
+ for (size_t l = 0; l < 6; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 6;
+ break;
+ }
+ case 6: {
+ Assert(value_list.size() == 3);
+ const double l0 = value_list[0];
+ const double l1 = value_list[1];
+ const double l2 = 1 - l0 - l1;
+ const double z = value_list[2];
+
+ point_list[k + 0] = to_R3(l0 * A + l1 * B + l2 * C, +z);
+ point_list[k + 1] = to_R3(l0 * A + l2 * B + l1 * C, +z);
+ point_list[k + 2] = to_R3(l1 * A + l0 * B + l2 * C, +z);
+ point_list[k + 3] = to_R3(l1 * A + l2 * B + l0 * C, +z);
+ point_list[k + 4] = to_R3(l2 * A + l0 * B + l1 * C, +z);
+ point_list[k + 5] = to_R3(l2 * A + l1 * B + l0 * C, +z);
+ point_list[k + 6] = to_R3(l0 * A + l1 * B + l2 * C, -z);
+ point_list[k + 7] = to_R3(l0 * A + l2 * B + l1 * C, -z);
+ point_list[k + 8] = to_R3(l1 * A + l0 * B + l2 * C, -z);
+ point_list[k + 9] = to_R3(l1 * A + l2 * B + l0 * C, -z);
+ point_list[k + 10] = to_R3(l2 * A + l0 * B + l1 * C, -z);
+ point_list[k + 11] = to_R3(l2 * A + l1 * B + l0 * C, -z);
+
+ for (size_t l = 0; l < 12; ++l) {
+ weight_list[k + l] = w;
+ }
+
+ k += 12;
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid quadrature id");
+ }
+ }
+ }
+ };
+
+ switch (m_degree) {
+ case 0:
+ case 1: {
+ constexpr size_t nb_points = 1;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{1, 1.000000000000000e+00, {}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 2: {
+ constexpr size_t nb_points = 5;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 2.344355869392759e-01, {-8.431650688665664e-01}},
+ Descriptor{3, 1.770429420404827e-01, {+1.046424703769979e-01}}};
+
+ 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 = 8;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.803034341765672e-01, {-9.614404179888022e-01}},
+ Descriptor{3, 1.134313729984015e-01, {+4.890588576053607e-01}},
+ Descriptor{3, 9.969967088388698e-02, {+7.783177802730620e-02}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 4: {
+ constexpr size_t nb_points = 11;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.079119748155355e-01, {-8.668619740090348e-01}},
+ Descriptor{3, 1.364146126054776e-01, {+4.686558098619952e-01}},
+ Descriptor{4, 6.248870209208267e-02, {+1.007404057989106e-01, +6.756398236822597e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 5: {
+ constexpr size_t nb_points = 16;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{1, 2.071428343483058e-01, {}}, Descriptor{3, 3.807558903099764e-02, {+5.176461782716475e-02}},
+ Descriptor{4, 7.637426139226190e-02, {+1.663967696311171e-01, +8.071634863884439e-01}},
+ Descriptor{4, 3.673080503418831e-02, {+4.976649895838920e-01, -3.972616744496609e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 6: {
+ constexpr size_t nb_points = 28;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 3.614462938209505e-02, {+9.850385924415599e-01}},
+ Descriptor{2, 5.544690202422076e-02, {-5.039978952117941e-01}},
+ Descriptor{3, 1.164296357658436e-02, {+1.750424465865124e-02}},
+ Descriptor{3, 7.688064192655711e-02, {+1.617417813899514e-01}},
+ Descriptor{4, 4.962685514962321e-02, {+4.656535513495914e-01, +4.811426008466984e-01}},
+ Descriptor{6, 2.112374914835039e-02, {+3.459486985245699e-02, +2.025039451729335e-01, -8.094634904091534e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 7: {
+ constexpr size_t nb_points = 35;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 2.844568111268929e-02, {+9.802280695908916e-01}},
+ Descriptor{3, 6.118242173095039e-03, {+8.337595466021473e-03}},
+ Descriptor{4, 2.155086408622646e-02, {+4.815219753291366e-01, -8.413873542065260e-01}},
+ Descriptor{4, 2.917852490209846e-02, {+9.548324837148936e-02, -7.958490905869831e-01}},
+ Descriptor{5, 2.551485633514925e-02, {+7.429966820728956e-01, +1.214913159837829e-02}},
+ Descriptor{6, 3.894070327620761e-02, {+1.529845984247976e-01, +3.051562164322261e-01, +4.039345605321099e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 8: {
+ constexpr size_t nb_points = 46;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 2.073493664285546e-02, {+2.244338661059906e-01}},
+ Descriptor{2, 5.108942711218155e-02, {-6.818254415708658e-01}},
+ Descriptor{3, 5.269974490650717e-02, {+4.600889628137106e-01}},
+ Descriptor{3, 1.815224878414972e-02, {+5.341235093690714e-02}},
+ Descriptor{4, 7.145357753538616e-03, {+4.723878583976938e-02, +8.288154430586802e-01}},
+ Descriptor{4, 4.062926265898933e-02, {+1.740616079243704e-01, +4.060842293545903e-01}},
+ Descriptor{4, 1.114302734846531e-02, {+1.597492639425890e-01, +9.658651665406544e-01}},
+ Descriptor{4, 2.108650929358649e-02, {+4.585690687909513e-01, +8.761959100002542e-01}},
+ Descriptor{6, 1.364752909087307e-02, {+8.588127507759013e-03, +7.285980718010000e-01, +5.773502691896257e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 9: {
+ constexpr size_t nb_points = 59;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 3.428347129092372e-02, {-2.750184190077768e-01}},
+ Descriptor{3, 1.639279770445363e-02, {+4.938630638969568e-01}},
+ Descriptor{4, 2.347178069870322e-02, {+2.362540169543293e-01, +8.712387576289576e-01}},
+ Descriptor{4, 5.761848158278895e-03, {+7.223833941638236e-02, +9.551937402361604e-01}},
+ Descriptor{4, 3.492479305778020e-02, {+4.383137607101617e-01, +4.364117817130079e-01}},
+ Descriptor{4, 7.373899623212321e-03, {+3.643401649407788e-02, +5.030056360566418e-01}},
+ Descriptor{4, 5.722460846448911e-03, {+4.828779929693860e-01, -9.817457230015865e-01}},
+ Descriptor{4, 2.433874301444529e-02, {+1.628698857202373e-01, -1.748668736196710e-01}},
+ Descriptor{5, 9.412963663135166e-03, {+8.213377527237301e-01, +1.626087609745086e-01}},
+ Descriptor{6, 1.801797749439730e-02, {+4.249444950639278e-02, +2.561926710584905e-01, +6.836158559766136e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 10: {
+ constexpr size_t nb_points = 84;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{1, 3.219608823836797e-02, {}},
+ Descriptor{2, 1.323959532338741e-02, {+9.650866816169476e-01}},
+ Descriptor{3, 1.663071726413940e-02, {+4.300104731739727e-01}},
+ Descriptor{3, 2.237117375679376e-02, {+1.095696683547513e-01}},
+ Descriptor{3, 2.965144268028371e-03, {+1.581483663138682e-02}},
+ Descriptor{4, 2.041220712674361e-02, {+4.356018236697312e-01, +7.829035904287732e-01}},
+ Descriptor{4, 1.414434969765846e-02, {+1.485976165307029e-01, +7.673257532061416e-01}},
+ Descriptor{4, 3.545708824271026e-03, {+5.196648358223575e-02, +4.892913949246038e-01}},
+ Descriptor{4, 6.058080482992846e-03, {+4.990219694442680e-01, +6.916716292769725e-01}},
+ Descriptor{4, 3.555888039958573e-03, {+4.155250160277021e-02, +8.910010843120678e-01}},
+ Descriptor{4, 3.038102466521747e-02, {+2.343945196820784e-01, +4.651342767282841e-01}},
+ Descriptor{6, 6.211856757092111e-03, {+5.281151684656214e-02, +2.716521744885937e-01, +9.566215690575565e-01}},
+ Descriptor{6, 6.618419166146100e-03, {+1.675738337212976e-01, +1.019125209869291e-02, +5.391797458851759e-01}},
+ Descriptor{6, 1.607306259567184e-02, {+3.291869417398026e-01, +5.090816276695182e-02, -2.642847151350913e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 11: {
+ constexpr size_t nb_points = 99;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 2.168014406730469e-02, {-5.549549808499665e-01}},
+ Descriptor{2, 1.045760798587228e-02, {+9.704322100480244e-01}},
+ Descriptor{2, 3.493462507713402e-03, {-8.225054517645587e-02}},
+ Descriptor{3, 1.610610418171415e-02, {+4.269221881685347e-01}},
+ Descriptor{4, 2.203310285466430e-02, {+4.396327659106838e-01, +6.759979260851990e-01}},
+ Descriptor{4, 1.274418696180960e-02, {+2.009128761138035e-01, +8.360535436001461e-01}},
+ Descriptor{4, 2.842526101831001e-03, {+2.158509804317693e-02, +5.861905524880673e-01}},
+ Descriptor{4, 1.509857722810566e-03, {+4.941546402080231e-01, +9.747001655491775e-01}},
+ Descriptor{4, 3.650427064786428e-03, {+5.996027726544360e-02, +9.450272677498353e-01}},
+ Descriptor{4, 2.100290671698337e-02, {+2.328102236807494e-01, +2.635224843921006e-01}},
+ Descriptor{4, 1.876195853667291e-02, {+1.213201391549184e-01, +4.439791020520624e-01}},
+ Descriptor{4, 5.194819129164601e-03, {+4.987543911906001e-01, +5.750153420806698e-01}},
+ Descriptor{5, 7.050587938770526e-03, {+1.046140524481813e-01, +1.651849633425114e-02}},
+ Descriptor{6, 5.901269948294771e-03, {+7.196325697558484e-02, +2.979854667459965e-01, +9.371332380619902e-01}},
+ Descriptor{6, 6.212665465530579e-03, {+1.998026706474004e-01, +1.720948255102629e-02, +7.302645437140292e-01}},
+ Descriptor{6, 1.385914960018441e-02, {+3.197359768880742e-01, +4.628523865251271e-02, -2.345966255449601e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 12: {
+ constexpr size_t nb_points = 136;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.122540095264984e-02, {+9.249605281775015e-01}},
+ Descriptor{2, 2.250978701634355e-02, {+3.174565048121913e-01}},
+ Descriptor{3, 5.489424083937044e-03, {+4.853372007888398e-02}},
+ Descriptor{3, 4.964406537953573e-04, {+3.340541016130764e-04}},
+ Descriptor{4, 1.020904001099484e-02, {+4.849074297077183e-01, +5.587032406470354e-01}},
+ Descriptor{4, 2.100649550809441e-03, {+1.497566246841529e-01, -9.990255354639939e-01}},
+ Descriptor{4, 1.520942066705075e-02, {+2.344536517846724e-01, +3.308248202440709e-01}},
+ Descriptor{4, 1.365828920018060e-02, {+4.409030162469282e-01, -7.263193646044468e-02}},
+ Descriptor{4, 2.788739322048965e-03, {+4.870412467653055e-01, +9.373420953931805e-01}},
+ Descriptor{4, 2.252095085593022e-03, {+2.486381930021292e-02, +6.814367870790086e-01}},
+ Descriptor{4, 9.917448258967848e-03, {+1.118542147928236e-01, +7.546508799087877e-01}},
+ Descriptor{5, 9.524594103832613e-03, {+1.207867185816364e-01, +7.097236881695401e-01}},
+ Descriptor{5, 6.688135429090100e-03, {+6.481099336610571e-01, +3.419306029008594e-01}},
+ Descriptor{6, 5.667322471331417e-03, {+2.009036502771764e-02, +1.335404714654308e-01, +3.889378443355619e-01}},
+ Descriptor{6, 1.258664545850983e-02, {+3.214917379706315e-01, +1.648190492804087e-01, +7.118012620024268e-01}},
+ Descriptor{6, 1.143701180003211e-02, {+6.899565054914573e-02, +6.636404691861656e-01, -4.223089212637487e-01}},
+ Descriptor{6, 4.495597556144068e-03, {+1.996263664334138e-02, +2.615529990296567e-01, +8.442893497066150e-01}},
+ Descriptor{6, 1.189774020910150e-03, {+9.002977238916547e-02, +1.552185633457250e-02, +9.525314828501794e-01}},
+ Descriptor{6, 4.663778699523011e-03, {+1.172214513692467e-01, +5.641016399337317e-01, -9.443120180277191e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 13: {
+ constexpr size_t nb_points = 162;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{4, 2.629059706793623e-03, {+4.898439392870481e-01, +9.257498273505042e-01}},
+ Descriptor{4, 5.973423749312596e-03, {+4.652402431593082e-01, -6.807809122309588e-01}},
+ Descriptor{4, 8.093549424155690e-03, {+1.189161740974999e-01, +8.113054084779335e-01}},
+ Descriptor{4, 6.819116259035672e-03, {+4.065656818638698e-01, +1.925054913726048e-02}},
+ Descriptor{4, 4.933451368773049e-03, {+3.731243598486834e-01, +9.244211981221383e-01}},
+ Descriptor{4, 2.014708013288427e-03, {+2.222577118367550e-02, +2.903696330212439e-01}},
+ Descriptor{4, 1.233577073707530e-02, {+4.252710267490136e-01, +6.716647847389562e-01}},
+ Descriptor{4, 1.136479433767108e-02, {+2.896929731593648e-01, +3.940165471189255e-01}},
+ Descriptor{4, 6.666864960682452e-03, {+1.072868932263340e-01, +1.612419567102727e-01}},
+ Descriptor{4, 1.194522013639762e-02, {+2.080125480772502e-01, +2.621034383425459e-01}},
+ Descriptor{4, 1.068788094061122e-02, {+2.230030950549982e-01, +7.517293034088793e-01}},
+ Descriptor{4, 1.666344327463221e-03, {+2.617439160849268e-02, +8.530545934139608e-01}},
+ Descriptor{5, 4.512051617065152e-03, {+6.914170159250816e-03, +2.742659465495736e-01}},
+ Descriptor{5, 3.488706839443290e-03, {+1.210652398684976e-01, +2.342377584496339e-02}},
+ Descriptor{5, 1.118035069441576e-02, {+9.353601775463583e-02, +5.415308571572887e-01}},
+ Descriptor{6, 4.205260232450528e-03, {+1.908095378196189e-02, +4.425451813256520e-01, +3.662111224258147e-01}},
+ Descriptor{6, 4.376791749742736e-03, {+2.923396969545124e-01, +1.817765226028588e-02, +7.893982266689565e-01}},
+ Descriptor{6, 4.406755041936957e-03, {+2.123093816715971e-02, +8.615902483468726e-01, +5.566323460176222e-01}},
+ Descriptor{6, 1.203503267726805e-02, {+2.541452627735283e-01, +7.525507194012332e-02, -4.323213872909485e-01}},
+ Descriptor{6, 1.398692188464922e-03, {+1.414034499801619e-01, +2.485309321657718e-02, +9.754978010051167e-01}},
+ Descriptor{6, 4.755154887378063e-03, {+5.920519581168684e-01, +1.079442302776815e-01, +9.602045342477669e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 14: {
+ constexpr size_t nb_points = 194;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.177233712035789e-02, {-5.431508197665302e-01}},
+ Descriptor{3, 1.353619130586776e-02, {+1.531409494078462e-01}},
+ Descriptor{3, 2.863624459018961e-03, {+1.094642266754876e-01}},
+ Descriptor{4, 2.831710195314422e-03, {+4.852139749601931e-01, +9.499413855584002e-01}},
+ Descriptor{4, 4.354649559352045e-03, {+4.811857851081537e-01, -3.416804188824072e-01}},
+ Descriptor{4, 9.145278212110571e-04, {+8.542525491697950e-02, +9.993425247803482e-01}},
+ Descriptor{4, 1.640665865351787e-03, {+4.949498254464672e-01, +4.043764192878801e-01}},
+ Descriptor{4, 1.217329578934041e-02, {+2.006722689888837e-01, +6.704323691216523e-01}},
+ Descriptor{4, 1.062128023171532e-03, {+1.346314023891985e-02, +4.402219142928553e-01}},
+ Descriptor{4, 8.146219031353574e-03, {+3.888073656098609e-01, +8.664942125245655e-01}},
+ Descriptor{4, 1.188198250865211e-02, {+2.589597694248802e-01, -3.805332165701859e-01}},
+ Descriptor{4, 4.510606073543659e-03, {+6.566746391585167e-02, +3.180500426640964e-01}},
+ Descriptor{4, 1.201058503487853e-02, {+4.462504387663027e-01, +6.737378781532631e-01}},
+ Descriptor{4, 4.766438302156675e-03, {+1.155506883468244e-01, +8.527706525272908e-01}},
+ Descriptor{4, 9.777024948606015e-04, {+2.579397983791752e-02, +9.172594869558284e-01}},
+ Descriptor{5, 5.111880702605887e-03, {+2.313133650799548e-02, +3.737411598567404e-01}},
+ Descriptor{5, 2.171762107022176e-03, {+7.808376567245785e-02, +8.716783743732383e-03}},
+ Descriptor{5, 1.040227457387526e-02, {+1.932090397525293e-01, +4.408189862889988e-01}},
+ Descriptor{6, 3.590362689299532e-03, {+1.106956102435793e-02, +3.495104935335414e-01, +6.888113400872247e-01}},
+ Descriptor{6, 5.210848278113965e-03, {+2.813786819004687e-01, +7.252407174904094e-02, +8.818281117320524e-01}},
+ Descriptor{6, 2.902421092721445e-03, {+1.360674597305564e-02, +7.694543610610796e-01, +2.082229797704095e-01}},
+ Descriptor{6, 1.093709739870311e-02, {+3.160642637515643e-01, +9.806129213500263e-02, -2.960531739682007e-01}},
+ Descriptor{6, 1.070623154200375e-03, {+2.010219519914552e-01, +1.015338057469439e-02, +9.567691001326991e-01}},
+ Descriptor{6, 2.366054231962321e-03, {+5.404862370418329e-01, +1.675040220811526e-01, +9.843823434581401e-01}},
+ Descriptor{6, 2.619672090010859e-03, {+1.016331118955575e-01, +1.682509189588544e-02, +7.176031930976530e-01}},
+ Descriptor{6, 7.096030229028862e-03, {+7.553391402298588e-01, +5.528506804337090e-02, -4.978572721822902e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 15: {
+ constexpr size_t nb_points = 238;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.097348853586182e-02, {+6.802218229609156e-01}},
+ Descriptor{2, 3.261395830734888e-03, {+9.049666192679098e-01}},
+ Descriptor{3, 9.394140546438022e-03, {+2.180179116942069e-01}},
+ Descriptor{3, 7.472192538440717e-03, {+1.819996803717089e-01}},
+ Descriptor{3, 1.125148661004542e-03, {+1.740140116312610e-02}},
+ Descriptor{3, 4.027894259445648e-03, {+5.427115018599665e-02}},
+ Descriptor{4, 4.698740527412616e-03, {+3.702896314191779e-01, +1.017726699824298e-01}},
+ Descriptor{4, 4.423106873185891e-03, {+6.837923242190364e-02, +7.328812049735380e-01}},
+ Descriptor{4, 1.190871191679164e-02, {+4.041029865296590e-01, -3.417820508269636e-01}},
+ Descriptor{4, 3.419058355265174e-03, {+4.924009437532922e-01, -1.425907114174295e-01}},
+ Descriptor{4, 7.607869583069166e-03, {+4.679592294726652e-01, +6.612792268372575e-01}},
+ Descriptor{4, 3.487721814344312e-03, {+2.538613877625208e-01, +9.653983029752667e-01}},
+ Descriptor{4, 4.126571123768199e-03, {+1.490233678882002e-01, +9.207186791179688e-01}},
+ Descriptor{4, 2.055229043939243e-03, {+4.943016520035221e-01, +8.528797591642616e-01}},
+ Descriptor{5, 5.333436981117544e-03, {+2.380227290434666e-01, +1.919120975700847e-02}},
+ Descriptor{6, 4.565836294062908e-03, {+2.658819381343187e-01, +9.470242762478344e-02, +8.051851008291565e-01}},
+ Descriptor{6, 5.459132755230557e-04, {+2.824880313716481e-02, +6.797646253429619e-03, +5.550126100536413e-01}},
+ Descriptor{6, 4.319194829106329e-03, {+6.410939576871652e-01, +1.511219453234863e-02, +4.965601017444452e-01}},
+ Descriptor{6, 5.377414350253174e-03, {+1.570345607673937e-01, +7.920488512684705e-02, +2.807402393928705e-01}},
+ Descriptor{6, 7.161694599840893e-03, {+2.618681572601245e-01, +1.978262294209248e-01, +5.539260799071125e-01}},
+ Descriptor{6, 5.425741614460213e-04, {+6.532458874680514e-03, +4.868056225266830e-02, +9.130997730502041e-01}},
+ Descriptor{6, 1.437669184909563e-03, {+3.702624031957258e-01, +9.129771330289427e-02, +9.857586619210927e-01}},
+ Descriptor{6, 1.116039038725027e-03, {+1.377079768509805e-01, +3.730083695338125e-02, +9.761262490766817e-01}},
+ Descriptor{6, 4.855981285749960e-03, {+4.784950101267365e-01, +1.635654688089934e-01, +8.308622970412985e-01}},
+ Descriptor{6, 1.797153116855535e-03, {+3.278240979071815e-01, +2.496158166292168e-02, +9.410500001533124e-01}},
+ Descriptor{6, 9.605058598669360e-03, {+3.475735063324100e-01, +8.713022883913982e-02, -2.293965609348112e-01}},
+ Descriptor{6, 2.702127731284175e-03, {+1.096327891311607e-01, +1.283956773909193e-02, +3.801848607953903e-01}},
+ Descriptor{6, 5.504446034785920e-03, {+2.078784801536998e-01, +8.490456755292054e-02, +5.633005891537184e-01}},
+ Descriptor{6, 2.394682993576167e-03, {+1.953330193360589e-01, +1.430457871757475e-02, +7.767703733759026e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 16: {
+ constexpr size_t nb_points = 287;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 1.260828009387497e-02, {+6.806235853917343e-01}},
+ Descriptor{3, 5.804416586274102e-04, {+8.954084031214650e-03}},
+ Descriptor{3, 1.451527527703537e-03, {+5.539081004402164e-02}},
+ Descriptor{3, 1.103925194575912e-02, {+3.916317627335951e-01}},
+ Descriptor{4, 2.389660788445647e-03, {+5.610007788583466e-02, +7.685450672389825e-01}},
+ Descriptor{4, 1.118529592275642e-02, {+2.669967593271997e-01, +3.187166719009041e-01}},
+ Descriptor{4, 4.452710626148950e-03, {+1.619375672046858e-01, +7.425213340738127e-01}},
+ Descriptor{4, 9.931042920523256e-03, {+4.259891438745869e-01, -5.391540839796620e-01}},
+ Descriptor{4, 2.798331509865767e-03, {+4.943682456985981e-01, +2.910077027212657e-01}},
+ Descriptor{4, 4.595079144000768e-03, {+4.745817664118026e-01, +7.702843820193086e-01}},
+ Descriptor{4, 2.649741982735800e-03, {+2.785739189962955e-01, +9.637659513314421e-01}},
+ Descriptor{4, 2.379174945004374e-03, {+1.564426776295505e-01, +9.700236378009609e-01}},
+ Descriptor{4, 1.290211537994051e-03, {+4.933597406415689e-01, +9.402630634429575e-01}},
+ Descriptor{5, 7.174970445767767e-03, {+5.254267651000835e-01, +9.931494549031464e-02}},
+ Descriptor{5, 3.904196912906489e-03, {+3.482743570815052e-01, +2.127155614722863e-02}},
+ Descriptor{5, 5.896847921349438e-03, {+1.654756461617898e-01, +8.608509799836524e-02}},
+ Descriptor{6, 2.434940232310643e-03, {+1.639763763638563e-01, +5.517222908705197e-02, +8.775190479419079e-01}},
+ Descriptor{6, 7.867926281754163e-04, {+3.703825228193449e-02, +1.050366659141256e-02, +5.241361844419883e-01}},
+ Descriptor{6, 1.973345895093291e-03, {+6.350080118967869e-01, +5.563959422105249e-03, +6.546900713490693e-01}},
+ Descriptor{6, 3.902804476486658e-03, {+1.297894012222233e-01, +7.248681972075718e-02, +4.803772456633624e-01}},
+ Descriptor{6, 2.865194848993297e-03, {+3.059321030072774e-01, +6.212018826128094e-02, +8.420714610259229e-01}},
+ Descriptor{6, 4.265862973729131e-03, {+1.677554820715698e-01, +2.612303617350463e-01, +6.466653853149219e-01}},
+ Descriptor{6, 2.763846307136280e-04, {+1.529621929633006e-03, +3.376869229909243e-02, +8.957708274465785e-01}},
+ Descriptor{6, 1.594720494514084e-03, {+3.532796606132861e-01, +9.492422215886979e-02, +9.832783024215900e-01}},
+ Descriptor{6, 5.967018601525356e-04, {+9.592623259204197e-02, +2.488171702893576e-02, +9.797358781454657e-01}},
+ Descriptor{6, 5.089010857368133e-03, {+4.963305263335406e-01, +1.757029526928022e-01, +8.611424584617942e-01}},
+ Descriptor{6, 1.023252728565292e-03, {+2.738501035919119e-01, +1.407861220922605e-02, +9.516023826285763e-01}},
+ Descriptor{6, 4.833611837936114e-03, {+3.828844998895772e-01, +5.651716232003427e-02, -3.926210681927463e-01}},
+ Descriptor{6, 2.256570449498587e-03, {+8.993653239530038e-02, +1.924125702016291e-02, +2.277368351329287e-01}},
+ Descriptor{6, 4.394217441895532e-03, {+2.516074150213084e-01, +6.113701845424014e-02, +5.771276473820092e-01}},
+ Descriptor{6, 1.441280420852857e-03, {+1.524244609819593e-01, +5.018269956194812e-03, +7.029237083896506e-01}},
+ Descriptor{6, 7.638903866482450e-03, {+2.475011786624472e-01, +1.365910574355268e-01, -2.555004528315441e-01}},
+ Descriptor{6, 3.266920063147978e-03, {+7.591852986956938e-01, +1.735437181095786e-02, +2.659716737332686e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 17: {
+ constexpr size_t nb_points = 338;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 7.521652592973654e-03, {+7.332980977058475e-01}},
+ Descriptor{4, 3.952340142081542e-03, {+4.665556153120432e-01, +7.407048594232563e-01}},
+ Descriptor{4, 2.658797102389008e-03, {+4.923809674098066e-01, +1.925784688493710e-01}},
+ Descriptor{4, 2.303261680132425e-03, {+4.466963397682847e-01, +9.724914617208265e-01}},
+ Descriptor{4, 4.300609620564292e-03, {+2.058899494758731e-01, +2.386810048778012e-01}},
+ Descriptor{4, 3.022106015301819e-03, {+2.305779360323416e-01, +3.012229853956593e-01}},
+ Descriptor{4, 3.076949950084475e-04, {+7.072650234259644e-03, +6.517254037629554e-01}},
+ Descriptor{4, 3.073440319150984e-03, {+5.248385455049231e-02, -5.044383522897413e-01}},
+ Descriptor{4, 4.905529550079285e-03, {+2.780743897936196e-01, -4.157283981367805e-01}},
+ Descriptor{4, 2.706669648409243e-03, {+3.893995813144771e-01, -8.036253483433273e-01}},
+ Descriptor{4, 3.457774927613496e-03, {+1.132999955833483e-01, +4.250397278118490e-01}},
+ Descriptor{4, 5.170517454197127e-03, {+4.687774749283292e-01, -2.967208107336694e-01}},
+ Descriptor{4, 1.402258311530629e-03, {+4.942337889964747e-01, +8.919484290177192e-01}},
+ Descriptor{4, 4.240102604381064e-03, {+2.906999098818323e-01, +1.317224598994015e-01}},
+ Descriptor{4, 3.284148741170156e-03, {+9.376930399184973e-02, +8.547449509711846e-01}},
+ Descriptor{4, 1.134225985312425e-03, {+1.632812980104028e-01, +9.979247130558654e-01}},
+ Descriptor{4, 7.470854816360809e-03, {+1.774626418569745e-01, -6.791385858496140e-01}},
+ Descriptor{4, 3.434038461352071e-03, {+3.930786561208929e-01, +2.690582708305222e-01}},
+ Descriptor{4, 3.146260969412623e-03, {+2.711153947647855e-01, +9.623986887607461e-01}},
+ Descriptor{4, 5.752613942102343e-04, {+2.430622100156484e-02, +9.425587161003264e-01}},
+ Descriptor{5, 6.555279553032400e-03, {+1.273859447314594e-01, +2.765687474993688e-01}},
+ Descriptor{6, 2.945409561987205e-03, {+2.450510833431427e-01, +8.524685018280553e-02, +9.106045966534763e-01}},
+ Descriptor{6, 2.453063041563170e-03, {+3.192911229101512e-01, +2.764498521125525e-02, +5.755372475087557e-02}},
+ Descriptor{6, 2.870396640498966e-03, {+4.577741357714497e-01, +1.385428807092584e-01, -9.235677741171985e-02}},
+ Descriptor{6, 1.826827790065983e-03, {+2.421105196372634e-01, +1.245683910267135e-02, -2.750816491221573e-01}},
+ Descriptor{6, 4.684287720797846e-03, {+7.269945541994681e-01, +8.880930240992388e-02, +2.353034707366664e-01}},
+ Descriptor{6, 6.198783236777368e-03, {+4.870517086755036e-01, +1.725779658472295e-01, -5.739332375158235e-01}},
+ Descriptor{6, 4.197863049050683e-03, {+1.582123188291182e-01, +3.108975273704328e-01, +8.562321460031186e-01}},
+ Descriptor{6, 2.552814201463011e-03, {+5.973833093262283e-01, +4.329678563553625e-02, +7.520938435967152e-01}},
+ Descriptor{6, 2.090115170917154e-03, {+1.111929938785087e-01, +3.682608320557917e-02, +2.115412627288322e-02}},
+ Descriptor{6, 7.568671403396213e-04, {+3.878650062905668e-02, +7.045410253068108e-03, +2.249484087437713e-01}},
+ Descriptor{6, 1.300065728839461e-03, {+1.417441449927015e-01, +6.656135005864013e-03, +3.824194976983090e-01}},
+ Descriptor{6, 8.479873802967710e-04, {+1.343472140817544e-01, +2.698141838221961e-02, +9.741848074979857e-01}},
+ Descriptor{6, 6.404613040101877e-03, {+3.195031824862086e-01, +8.134405613722046e-02, +4.381759877407309e-01}},
+ Descriptor{6, 1.132352094196614e-03, {+8.263084385883397e-02, +1.026354650397643e-02, +7.684631026560178e-01}},
+ Descriptor{6, 2.064410253144382e-03, {+3.827807429336583e-01, +6.749013415056980e-03, +5.465728592769658e-01}},
+ Descriptor{6, 1.048193458332381e-03, {+2.446543644502634e-01, +5.650199208693200e-03, +8.553667006980621e-01}},
+ Descriptor{6, 9.390285812952761e-04, {+2.637071367985937e-02, +6.264598409479015e-01, +9.755523282064171e-01}},
+ Descriptor{6, 4.216060332324916e-03, {+4.340745420159236e-02, +1.947680793017931e-01, +6.465743617002622e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 18: {
+ constexpr size_t nb_points = 396;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{4, 2.257692848096014e-03, {+4.728631477307805e-01, +9.068774432765224e-01}},
+ Descriptor{4, 3.522000993611865e-03, {+1.090318158369727e-01, +5.191310054808369e-01}},
+ Descriptor{4, 4.461723497786998e-03, {+2.141954002982068e-01, +6.975267033041567e-01}},
+ Descriptor{4, 3.707999359392404e-04, {+8.415088753342596e-03, +4.704963597796460e-01}},
+ Descriptor{4, 2.423661142721756e-03, {+4.941564402940962e-02, -3.771046294090744e-01}},
+ Descriptor{4, 2.814996700016239e-03, {+4.574258170355346e-01, -2.709290831593626e-01}},
+ Descriptor{4, 5.337642166957648e-03, {+2.945570490820424e-01, -1.757797989950063e-01}},
+ Descriptor{4, 8.574266664464696e-03, {+4.148583469418042e-01, -2.903119874076979e-01}},
+ Descriptor{4, 2.265725871486650e-03, {+1.801567120271234e-01, -2.996313713953939e-01}},
+ Descriptor{4, 2.955157959707394e-03, {+4.092592655375933e-01, -8.687914069087144e-01}},
+ Descriptor{4, 3.919466496845992e-03, {+4.106876572727113e-01, -7.102916194693503e-01}},
+ Descriptor{4, 6.030674708196229e-04, {+4.942616172602019e-01, +9.613370670800333e-01}},
+ Descriptor{4, 4.577980994542588e-03, {+1.455821883151533e-01, +2.580797174139677e-01}},
+ Descriptor{4, 2.430663642304682e-03, {+7.759286964629900e-02, +8.574527537341975e-01}},
+ Descriptor{4, 2.159969925839523e-03, {+1.688466480848478e-01, +9.442832226314590e-01}},
+ Descriptor{4, 7.901051089233524e-03, {+2.737821968608845e-01, -5.261705496722796e-01}},
+ Descriptor{4, 2.173845391531516e-03, {+3.193703110432156e-01, +7.838756187326278e-01}},
+ Descriptor{4, 2.248646315174777e-03, {+2.789278090499041e-01, +9.684120338438452e-01}},
+ Descriptor{4, 3.442548195872391e-04, {+1.473433086781356e-02, +9.055493494371784e-01}},
+ Descriptor{5, 8.069935581612784e-03, {+2.803931989784251e-01, +1.730708483947652e-01}},
+ Descriptor{6, 2.387723163590078e-03, {+3.067402740370590e-01, +4.725833161993727e-02, +8.686462758628247e-01}},
+ Descriptor{6, 3.430528221040949e-03, {+5.587472044782793e-01, +6.374864935475613e-02, +1.027445357789138e-01}},
+ Descriptor{6, 6.029545367216524e-03, {+1.295467219204834e-01, +2.834002542935580e-01, -4.512744270430775e-01}},
+ Descriptor{6, 3.276823879158102e-03, {+2.943165499688927e-01, +3.236061609245209e-02, -4.364849670272578e-01}},
+ Descriptor{6, 3.663549210540714e-03, {+7.045618883179275e-01, +6.827271159006545e-02, +1.782321853488997e-01}},
+ Descriptor{6, 4.424488262807732e-03, {+5.336092397180819e-01, +8.232278479338200e-02, -6.571528038169470e-01}},
+ Descriptor{6, 2.760063523214631e-03, {+1.659343236267790e-01, +2.958566157742775e-01, +8.788603720074640e-01}},
+ Descriptor{6, 1.392109889443924e-03, {+5.873912217089957e-01, +9.153855225433397e-03, +7.488835368067833e-01}},
+ Descriptor{6, 2.138435872642536e-03, {+1.216404492275597e-01, +4.924556835776799e-02, -8.479900295560096e-02}},
+ Descriptor{6, 5.166067073488927e-04, {+4.208095332260363e-02, +7.753734519664009e-03, +9.396055358943239e-02}},
+ Descriptor{6, 1.953953141133318e-03, {+4.516351625773616e-01, +2.282191697171249e-02, +4.454084028264259e-01}},
+ Descriptor{6, 1.255008598363820e-03, {+1.252580641579170e-01, +7.396327404727600e-03, +3.172925739593950e-01}},
+ Descriptor{6, 1.042334114181810e-03, {+1.765722195894258e-01, +6.530021157253776e-02, +9.730112939053892e-01}},
+ Descriptor{6, 1.279342041886178e-03, {+2.404267015692479e-01, +8.587901550681881e-03, +1.654571949278815e-01}},
+ Descriptor{6, 1.087772234137440e-03, {+6.440585522631566e-02, +1.226780343665597e-02, +7.058192837170146e-01}},
+ Descriptor{6, 8.506443220554228e-04, {+2.681069831840066e-01, +4.479326561367245e-03, +6.566582555690533e-01}},
+ Descriptor{6, 1.071738436276040e-03, {+1.727205252244992e-01, +1.173345352446697e-02, +8.740347679939118e-01}},
+ Descriptor{6, 1.222389838372311e-03, {+1.118329953359960e-01, +5.446775610503505e-01, +9.879516615930261e-01}},
+ Descriptor{6, 1.426482164900413e-03, {+5.899536365475773e-01, +8.017587417702520e-03, +1.572003893989445e-01}},
+ Descriptor{6, 3.834177899773436e-03, {+9.615405978617643e-02, +2.067498834662211e-01, +7.510430318557424e-01}},
+ Descriptor{6, 4.758021594541443e-04, {+3.068968355307726e-01, +1.177465532391100e-02, +9.786522981024399e-01}},
+ Descriptor{6, 3.711298716255058e-04, {+7.984864478527595e-02, +1.740633561957181e-02, +9.787643170060012e-01}},
+ Descriptor{6, 2.736409660029034e-03, {+3.607777906277682e-02, +1.618562004775664e-01, +5.932565552690645e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 19: {
+ constexpr size_t nb_points = 420;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 6.766510972429800e-03, {+5.244830821459505e-01}},
+ Descriptor{2, 2.666997686915707e-03, {+7.453333881549420e-02}},
+ Descriptor{2, 2.198296683612197e-03, {+9.566097292149883e-01}},
+ Descriptor{3, 1.803557661024045e-03, {+4.935398626809867e-01}},
+ Descriptor{3, 6.354430508374619e-03, {+2.767659367008514e-01}},
+ Descriptor{3, 1.636834721720784e-03, {+5.457779485028427e-02}},
+ Descriptor{3, 4.893654007837226e-04, {+1.399697352188861e-02}},
+ Descriptor{4, 1.846418513040312e-03, {+1.418740786751734e-01, +9.552594993861951e-01}},
+ Descriptor{4, 2.773402798745861e-04, {+8.395599010104145e-03, +6.731987542402844e-01}},
+ Descriptor{4, 2.823634437215422e-03, {+1.035973359771996e-01, -2.005717856238193e-01}},
+ Descriptor{4, 7.726075998933063e-03, {+3.956080996828434e-01, +3.080298478277679e-01}},
+ Descriptor{4, 2.193983773241012e-04, {+1.400193119774357e-02, +2.354401887617481e-01}},
+ Descriptor{4, 4.366650366888420e-03, {+4.801897764603656e-01, +4.828040748109914e-01}},
+ Descriptor{4, 4.029293922819735e-03, {+1.623484102426076e-01, -2.002124012363609e-01}},
+ Descriptor{4, 1.130768692395794e-03, {+4.972236736522475e-01, +7.211608708255114e-01}},
+ Descriptor{4, 4.287863606686208e-03, {+1.990216655966096e-01, +7.769851170153003e-01}},
+ Descriptor{4, 1.592256072404633e-03, {+1.202889101503504e-01, -6.997989522584932e-01}},
+ Descriptor{4, 2.052668995583038e-03, {+5.321915260988313e-02, +6.322324604037224e-01}},
+ Descriptor{4, 3.894795325529755e-03, {+4.012796612843944e-01, +8.985968033126313e-01}},
+ Descriptor{4, 1.803402089659100e-03, {+1.089883669604236e-01, +6.414118273921905e-01}},
+ Descriptor{4, 4.906673814096024e-03, {+4.472748897677901e-01, +1.418065900083249e-01}},
+ Descriptor{4, 8.385245625735418e-04, {+6.923354485715784e-02, +9.592516295656434e-01}},
+ Descriptor{4, 1.344422701677694e-03, {+2.475145234973103e-01, +9.911400740550351e-01}},
+ Descriptor{4, 1.212351202429953e-03, {+4.533377451558121e-01, +9.881823361248698e-01}},
+ Descriptor{4, 2.704967330298844e-04, {+1.808374112681072e-02, +9.537067614804079e-01}},
+ Descriptor{5, 1.194495738162939e-03, {+1.297897131077726e-01, +6.917496569505294e-03}},
+ Descriptor{5, 3.103142225008695e-03, {+3.779418721990778e-01, +3.100617423599103e-02}},
+ Descriptor{5, 5.131660654552337e-03, {+6.446313119552312e-02, +2.203251202566757e-01}},
+ Descriptor{6, 4.803074338815450e-03, {+1.625899734623090e-01, +2.824754920851798e-01, +1.568351483023571e-01}},
+ Descriptor{6, 8.150376242849885e-04, {+8.118414809439785e-03, +5.478835112834492e-02, -3.967511256013889e-01}},
+ Descriptor{6, 4.498489872233412e-03, {+6.582273496659641e-01, +1.178005364736766e-01, -4.885151647510187e-01}},
+ Descriptor{6, 2.731777351387361e-03, {+2.876874815978630e-01, +3.552685414544682e-02, +5.867894747260645e-01}},
+ Descriptor{6, 3.917995492165165e-03, {+2.839484400308330e-01, +2.073917255967352e-01, +4.831782153742011e-01}},
+ Descriptor{6, 8.946420147168489e-04, {+2.824570038349533e-02, +8.622300269568527e-02, -2.019344195912693e-01}},
+ Descriptor{6, 1.272619581981183e-03, {+3.915043117238656e-01, +5.619651603686126e-03, +3.317594047785773e-01}},
+ Descriptor{6, 4.539117516816274e-03, {+7.376676137243446e-02, +3.389265718024346e-01, -3.045944670757652e-01}},
+ Descriptor{6, 1.228077139941574e-03, {+1.582419524500696e-01, +8.020518897879757e-03, +6.242276968562098e-01}},
+ Descriptor{6, 2.481352686408294e-03, {+2.270465393157224e-01, +7.611830639869993e-02, +7.779863298475371e-01}},
+ Descriptor{6, 6.742435916065792e-04, {+1.074023202728770e-02, +6.748601969235382e-02, +8.305354371679037e-01}},
+ Descriptor{6, 3.085688768414200e-03, {+4.506297544312898e-02, +1.495289831435912e-01, +3.936943128499052e-01}},
+ Descriptor{6, 2.827744882828859e-03, {+3.899263834124326e-01, +5.365999012062311e-02, +8.362140373651296e-01}},
+ Descriptor{6, 7.880156766981593e-04, {+1.426168210785968e-02, +4.039638109964547e-01, +9.545849701585185e-01}},
+ Descriptor{6, 9.734356783437631e-04, {+3.314062155351160e-01, +8.833083116611448e-03, +7.096795149078843e-01}},
+ Descriptor{6, 2.344473397349397e-04, {+1.223644156202090e-01, +1.039237935936368e-02, +9.900280695189960e-01}},
+ Descriptor{6, 2.754928134962740e-03, {+2.929652829983945e-01, +1.324448044285512e-01, +9.178743530773790e-01}},
+ Descriptor{6, 8.530955496579550e-04, {+2.534748953911949e-01, +5.351470915991541e-02, +9.845698380380822e-01}},
+ Descriptor{6, 1.720791066604888e-03, {+4.407858720850005e-02, +1.461188761758480e-01, +8.569909006217998e-01}},
+ Descriptor{6, 5.238060160822010e-03, {+5.031457633505478e-01, +1.264604222462178e-01, +6.487185269932558e-01}},
+ Descriptor{6, 7.780035531900878e-04, {+8.749997984117155e-03, +2.433623618774898e-01, +9.028600070589583e-01}},
+ Descriptor{6, 2.799719934133974e-03, {+4.433279475253553e-01, +2.569581578007861e-01, +7.536500284803950e-01}},
+ Descriptor{6, 1.887126925839933e-03, {+2.485412372956145e-01, +1.127715403373798e-02, +2.418109616381699e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ case 20: {
+ constexpr size_t nb_points = 518;
+ SmallArray<R3> point_list(nb_points);
+ SmallArray<double> weight_list(nb_points);
+
+ std::array descriptor_list = //
+ {Descriptor{2, 4.747044232976959e-03, {+8.420166009462728e-01}},
+ Descriptor{3, 6.640322548036050e-04, {+1.441108736989831e-02}},
+ Descriptor{3, 9.620914054915759e-04, {+5.949741342252451e-02}},
+ Descriptor{4, 2.742173179580797e-03, {+4.716249531173324e-01, +8.152882860873994e-01}},
+ Descriptor{4, 2.814324102526247e-03, {+4.084460043805324e-01, +5.308959921730525e-01}},
+ Descriptor{4, 3.864759174126348e-03, {+2.544280410411358e-01, +7.145164808360277e-01}},
+ Descriptor{4, 3.774538054936480e-03, {+4.774877095488475e-01, +3.302023087843083e-01}},
+ Descriptor{4, 2.695893315348362e-04, {+8.029263668218678e-03, +5.253966144526213e-01}},
+ Descriptor{4, 1.117541442536487e-03, {+4.968613511179721e-01, +1.678823330502635e-01}},
+ Descriptor{4, 1.168183967320753e-03, {+4.952615769053836e-01, +7.059801252317832e-01}},
+ Descriptor{4, 9.039109611475474e-04, {+3.535780971713121e-02, +3.380124871900346e-01}},
+ Descriptor{4, 3.369833201547859e-03, {+2.045970282680664e-01, -9.824414827677468e-02}},
+ Descriptor{4, 2.639764171119988e-03, {+1.406131665045265e-01, +6.840978690759630e-01}},
+ Descriptor{4, 2.599301978924752e-03, {+9.825610533777711e-02, -2.793583988322257e-01}},
+ Descriptor{4, 1.282233014437331e-03, {+5.059516645881403e-02, +5.440523359665514e-01}},
+ Descriptor{4, 1.290196339749456e-03, {+4.494558038329708e-01, +9.661789228499832e-01}},
+ Descriptor{4, 3.634034328176076e-03, {+2.188296334498393e-01, +5.161030390440184e-01}},
+ Descriptor{4, 3.865555241109362e-03, {+3.839280551265914e-01, +6.414153617240698e-01}},
+ Descriptor{4, 1.054906754865126e-03, {+5.997065917075951e-02, +8.943918441529238e-01}},
+ Descriptor{4, 1.562311721738531e-03, {+2.735494152095392e-01, +9.793184961065893e-01}},
+ Descriptor{4, 4.729101647325770e-04, {+4.907393281824849e-01, +9.886251338305795e-01}},
+ Descriptor{4, 1.740638131578173e-04, {+1.149301083025444e-02, +9.267723967676791e-01}},
+ Descriptor{4, 2.895959965319792e-03, {+4.048678509173179e-01, +9.164077711364235e-01}},
+ Descriptor{4, 5.234744278809832e-03, {+2.917111538437177e-01, -3.721772106793935e-01}},
+ Descriptor{4, 3.176639444983119e-03, {+4.459340787080981e-01, +3.782961145575000e-02}},
+ Descriptor{5, 2.716223123318332e-03, {+2.107382124681030e-01, +3.137124037930802e-02}},
+ Descriptor{5, 3.281825797738248e-03, {+2.031438036613264e-01, +1.114400380412596e-01}},
+ Descriptor{5, 3.187894698489262e-03, {+3.734168270486775e-01, +3.629012071225871e-02}},
+ Descriptor{6, 3.590341132235904e-03, {+5.331515253263366e-01, +1.335223252058013e-01, +7.802503826507077e-01}},
+ Descriptor{6, 5.184823629819776e-04, {+3.105220296426577e-03, +6.506617704521390e-02, -2.884770252636346e-01}},
+ Descriptor{6, 2.044053800301851e-03, {+7.387825606686489e-01, +9.320636776575365e-02, -2.442855072953742e-01}},
+ Descriptor{6, 2.454384702584994e-03, {+2.402209181076414e-01, +3.428169043211691e-02, +4.260690982676672e-01}},
+ Descriptor{6, 4.119406732850509e-03, {+8.161298386181808e-02, +3.087274419444517e-01, +2.420309868035705e-01}},
+ Descriptor{6, 5.762946980041836e-03, {+1.728159430622389e-01, +3.228866908185781e-01, -2.377845513467044e-01}},
+ Descriptor{6, 3.254649507017313e-03, {+2.247127800382396e-01, +1.327427633967232e-01, +4.683664775890787e-01}},
+ Descriptor{6, 1.376952861188114e-03, {+2.850103512086270e-02, +1.058138309132070e-01, -1.189934946219634e-01}},
+ Descriptor{6, 1.390144266230514e-03, {+3.896100790636586e-01, +1.097122964947396e-02, +4.812562527542947e-01}},
+ Descriptor{6, 1.973482620122281e-03, {+2.493335646225524e-01, +3.397868885786173e-01, -3.494229691424887e-02}},
+ Descriptor{6, 1.112029053786251e-03, {+7.042023600039281e-03, +3.177207979690921e-01, +2.078247495083118e-01}},
+ Descriptor{6, 1.196110133680204e-03, {+1.229469656661441e-01, +1.566757860788512e-02, +6.094751351099792e-01}},
+ Descriptor{6, 1.888982141486211e-03, {+2.259160262165568e-01, +5.215188464418734e-02, +8.360830877198899e-01}},
+ Descriptor{6, 5.659037184544154e-04, {+9.471493326637022e-03, +4.964414152642235e-02, +7.613841735998289e-01}},
+ Descriptor{6, 1.832188797223636e-03, {+5.243114241272343e-02, +1.452721772236918e-01, +4.996792895702993e-01}},
+ Descriptor{6, 8.460253945737610e-04, {+8.872202555510618e-03, +3.816302433543118e-01, +8.877001109429070e-01}},
+ Descriptor{6, 1.597513197343773e-03, {+2.592518129005248e-01, +5.551397100853696e-01, +8.790486279875863e-01}},
+ Descriptor{6, 2.678157870344808e-03, {+3.527097683072426e-01, +4.546097143409199e-02, +6.369991526576181e-01}},
+ Descriptor{6, 7.373212126545897e-04, {+1.450709676542589e-01, +1.181702543851609e-02, +8.952124765779996e-01}},
+ Descriptor{6, 1.390059413378896e-03, {+1.874143309040665e-01, +1.371981166173774e-01, +9.134247047008914e-01}},
+ Descriptor{6, 7.964285504550415e-04, {+6.814016086596048e-02, +1.498311852045099e-01, +9.751784156687211e-01}},
+ Descriptor{6, 1.460116442574025e-03, {+3.340585216441213e-01, +5.496962861346901e-02, +9.356631230196827e-01}},
+ Descriptor{6, 1.348373684024378e-03, {+6.785157976700519e-02, +1.205147270278029e-01, +7.655747942493478e-01}},
+ Descriptor{6, 3.286050878128491e-03, {+5.023136055537912e-01, +1.076311573233394e-01, +5.182173603742516e-01}},
+ Descriptor{6, 1.124720118777690e-03, {+7.819399743424233e-03, +2.556542053755758e-01, +7.166820851155357e-01}},
+ Descriptor{6, 2.127107274890790e-03, {+6.609230621965876e-01, +9.839569585073644e-02, +6.803594853699388e-01}},
+ Descriptor{6, 7.135219223491392e-04, {+1.784736485903538e-01, +1.796539994312394e-03, +2.866535225594761e-01}},
+ Descriptor{6, 3.719588400835986e-04, {+2.532227359651285e-01, +1.312175908589084e-02, +9.829888431964284e-01}},
+ Descriptor{6, 2.151829924432428e-04, {+1.443505367133109e-02, +6.867387590589011e-02, +9.834928544568246e-01}},
+ Descriptor{6, 8.163226512583331e-04, {+2.794904874495157e-01, +1.274003485132245e-01, +9.907920620663775e-01}}};
+
+ fill_quadrature_points(descriptor_list, point_list, weight_list);
+
+ m_point_list = point_list;
+ m_weight_list = weight_list;
+ break;
+ }
+ default: {
+ throw NormalError("Gauss quadrature formulae handle degrees up to 20 on prisms");
+ }
+ }
+}
diff --git a/src/analysis/PrismGaussQuadrature.hpp b/src/analysis/PrismGaussQuadrature.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..96ae6fbee84fc22ac404153e44a8e1d49a846924
--- /dev/null
+++ b/src/analysis/PrismGaussQuadrature.hpp
@@ -0,0 +1,32 @@
+#ifndef PRISM_GAUSS_QUADRATURE_HPP
+#define PRISM_GAUSS_QUADRATURE_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+
+/**
+ * Defines Gauss quadrature on the reference prism element
+ *
+ * \note formulae are provided by
+ *
+ * 'High-order symmetric cubature rules for tetrahedra and pyramids'
+ * Jan JasĖkowiec & N. Sukumar (2020).
+ */
+class PrismGaussQuadrature final : public QuadratureForumla<3>
+{
+ private:
+ void _buildPointAndWeightLists();
+
+ public:
+ PrismGaussQuadrature(PrismGaussQuadrature&&) = default;
+ PrismGaussQuadrature(const PrismGaussQuadrature&) = default;
+
+ explicit PrismGaussQuadrature(const size_t degree) : QuadratureForumla<3>(degree)
+ {
+ this->_buildPointAndWeightLists();
+ }
+
+ PrismGaussQuadrature() = delete;
+ ~PrismGaussQuadrature() = default;
+};
+
+#endif // PRISM_GAUSS_QUADRATURE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3ae6ff750ecb45cf1f493e400fa452ed88e960f3..bbe9ee92269704f85964748817663619ffb21a1e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -95,6 +95,7 @@ add_executable (unit_tests
test_OStream.cpp
test_ParseError.cpp
test_PETScUtils.cpp
+ test_PrismGaussQuadrature.cpp
test_PugsAssert.cpp
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
diff --git a/tests/test_PrismGaussQuadrature.cpp b/tests/test_PrismGaussQuadrature.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c94f63c2d851b6299171fd651b111d8cdc397647
--- /dev/null
+++ b/tests/test_PrismGaussQuadrature.cpp
@@ -0,0 +1,684 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <analysis/PrismGaussQuadrature.hpp>
+#include <geometry/AffineTransformation.hpp>
+#include <utils/Exceptions.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("PrismGaussQuadrature", "[analysis]")
+{
+ auto integrate = [](auto f, auto quadrature_formula) {
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+
+ auto value = weight_list[0] * f(point_list[0]);
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(point_list[i]);
+ }
+
+ return value;
+ };
+
+ auto integrate_on_prism = [](auto f, auto quadrature_formula, const std::array<TinyVector<3>, 4>& tetrahedron) {
+ const auto& A = tetrahedron[0];
+ const auto& B = tetrahedron[1];
+ const auto& C = tetrahedron[2];
+ const auto& D = tetrahedron[3];
+
+ TinyMatrix<3> J;
+ for (size_t i = 0; i < 3; ++i) {
+ J(i, 0) = B[i] - A[i];
+ J(i, 1) = C[i] - A[i];
+ J(i, 2) = 0.5 * (D[i] - A[i]);
+ }
+ TinyVector s = 0.5 * (A + D);
+
+ auto point_list = quadrature_formula.pointList();
+ auto weight_list = quadrature_formula.weightList();
+ auto value = weight_list[0] * f(J * (point_list[0]) + s);
+ for (size_t i = 1; i < weight_list.size(); ++i) {
+ value += weight_list[i] * f(J * (point_list[i]) + s);
+ }
+
+ return det(J) * value;
+ };
+
+ auto get_order = [&integrate, &integrate_on_prism](auto f, auto quadrature_formula, const double exact_value) {
+ using R3 = TinyVector<3>;
+
+ const R3 A{+0, +0, -1};
+ const R3 B{+1, +0, -1};
+ const R3 C{+0, +1, -1};
+ const R3 D{+0, +0, +1};
+ const R3 E{+1, +0, +1};
+ const R3 F{+0, +1, +1};
+
+ const R3 M_AB = 0.5 * (A + B);
+ const R3 M_AC = 0.5 * (A + C);
+ const R3 M_BC = 0.5 * (B + C);
+ const R3 M_DE = 0.5 * (D + E);
+ const R3 M_AD = 0.5 * (A + D);
+ const R3 M_BE = 0.5 * (B + E);
+ const R3 M_CF = 0.5 * (C + F);
+ const R3 M_ABED = 0.25 * (A + B + E + D);
+ const R3 M_BCFE = 0.25 * (B + C + E + F);
+ const R3 M_ADFC = 0.25 * (A + D + F + C);
+
+ const double int_T_hat = integrate(f, quadrature_formula);
+ const double int_refined //
+ = integrate_on_prism(f, quadrature_formula, {A, M_AB, M_AC, M_AD}) +
+ integrate_on_prism(f, quadrature_formula, {B, M_BC, M_AB, M_BE}) +
+ integrate_on_prism(f, quadrature_formula, {C, M_AC, M_BC, M_CF}) +
+ integrate_on_prism(f, quadrature_formula, {M_AB, M_BC, M_AC, M_ABED}) +
+ integrate_on_prism(f, quadrature_formula, {M_AD, M_ABED, M_ADFC, D}) +
+ integrate_on_prism(f, quadrature_formula, {M_BE, M_BCFE, M_ABED, E}) +
+ integrate_on_prism(f, quadrature_formula, {M_CF, M_ADFC, M_BCFE, F}) +
+ integrate_on_prism(f, quadrature_formula, {M_ABED, M_BCFE, M_ADFC, M_DE});
+
+ return -std::log(std::abs(int_refined - exact_value) / std::abs(int_T_hat - exact_value)) / std::log(2);
+ };
+
+ auto p0 = [](const TinyVector<3>&) { return 4; };
+ auto p1 = [](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return 2 * x + 3 * y + z - 1;
+ };
+ auto p2 = [&p1](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p1(X) * (2.5 * x - 3 * y + z + 3);
+ };
+ auto p3 = [&p2](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p2(X) * (3 * x + 2 * y - 3 * z - 1);
+ };
+ auto p4 = [&p3](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p3(X) * (2 * x - 0.5 * y - 1.3 * z + 1);
+ };
+ auto p5 = [&p4](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p4(X) * (-0.1 * x + 1.3 * y - 3 * z + 1);
+ };
+ auto p6 = [&p5](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p5(X) * 7875. / 143443 * (2 * x - y + 4 * z + 1);
+ };
+ auto p7 = [&p6](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p6(X) * (0.7 * x - 2.7 * y + 1.3 * z - 2);
+ };
+ auto p8 = [&p7](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p7(X) * (0.3 * x + 1.2 * y - 0.7 * z + 0.2);
+ };
+ auto p9 = [&p8](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p8(X) * (-0.2 * x + 1.1 * y - 0.5 * z + 0.6);
+ };
+ auto p10 = [&p9](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p9(X) * (0.7 * x - 0.6 * y - 0.7 * z - 0.2);
+ };
+ auto p11 = [&p10](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p10(X) * (-1.3 * x + 0.6 * y - 1.3 * z + 0.7);
+ };
+ auto p12 = [&p11](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p11(X) * (0.3 * x - 0.7 * y + 0.3 * z + 0.7);
+ };
+ auto p13 = [&p12](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p12(X) * (0.9 * x + 0.2 * y - 0.4 * z + 0.5);
+ };
+ auto p14 = [&p13](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p13(X) * (0.6 * x - 1.2 * y + 0.7 * z - 0.4);
+ };
+ auto p15 = [&p14](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p14(X) * (-0.2 * x - 0.7 * y + 0.9 * z + 0.8);
+ };
+ auto p16 = [&p15](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p15(X) * (0.7 * x + 0.2 * y - 0.6 * z + 0.4);
+ };
+ auto p17 = [&p16](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p16(X) * (-0.1 * x + 0.8 * y + 0.3 * z - 0.2);
+ };
+ auto p18 = [&p17](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p17(X) * (0.7 * x - 0.2 * y - 0.3 * z + 0.8);
+ };
+ auto p19 = [&p18](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p18(X) * (-0.7 * x + 1.2 * y + 1.3 * z + 0.8);
+ };
+ auto p20 = [&p19](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p19(X) * (0.7 * x - 1.2 * y + 0.3 * z - 0.6);
+ };
+ auto p21 = [&p20](const TinyVector<3>& X) {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return p20(X) * (0.7 * x - 1.2 * y + 0.3 * z - 0.6);
+ };
+
+ SECTION("degree 1")
+ {
+ PrismGaussQuadrature l1(1);
+
+ REQUIRE(l1.numberOfPoints() == 1);
+
+ REQUIRE(integrate(p0, l1) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l1) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l1) != Catch::Approx(47. / 24));
+
+ REQUIRE(get_order(p2, l1, 47. / 24) == Catch::Approx(2));
+ }
+
+ SECTION("degree 2")
+ {
+ PrismGaussQuadrature l2(2);
+
+ REQUIRE(l2.numberOfPoints() == 5);
+
+ REQUIRE(integrate(p0, l2) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l2) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l2) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l2) != Catch::Approx(-427. / 360));
+
+ // In a weird way, one gets 4th order on this test
+ REQUIRE(get_order(p3, l2, -427. / 360) == Catch::Approx(4));
+ }
+
+ SECTION("degree 3")
+ {
+ PrismGaussQuadrature l3(3);
+
+ REQUIRE(l3.numberOfPoints() == 8);
+
+ REQUIRE(integrate(p0, l3) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l3) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l3) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l3) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l3) != Catch::Approx(1003. / 3600));
+
+ REQUIRE(get_order(p4, l3, 1003. / 3600) == Catch::Approx(4));
+ }
+
+ SECTION("degree 4")
+ {
+ PrismGaussQuadrature l4(4);
+
+ REQUIRE(l4.numberOfPoints() == 11);
+
+ REQUIRE(integrate(p0, l4) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l4) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l4) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l4) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l4) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l4) != Catch::Approx(34031. / 18000));
+
+ REQUIRE(get_order(p5, l4, 34031. / 18000) >= Catch::Approx(5));
+ }
+
+ SECTION("degree 5")
+ {
+ PrismGaussQuadrature l5(5);
+
+ REQUIRE(l5.numberOfPoints() == 16);
+
+ REQUIRE(integrate(p0, l5) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l5) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l5) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l5) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l5) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l5) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l5) != Catch::Approx(36346369. / 55082112));
+
+ REQUIRE(get_order(p6, l5, 36346369. / 55082112) == Catch::Approx(6));
+ }
+
+ SECTION("degree 6")
+ {
+ PrismGaussQuadrature l6(6);
+
+ REQUIRE(l6.numberOfPoints() == 28);
+
+ REQUIRE(integrate(p0, l6) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l6) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l6) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l6) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l6) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l6) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l6) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l6) != Catch::Approx(-1128861017. / 918035200));
+
+ // In a weird way, one gets 8th order on this test
+ REQUIRE(get_order(p7, l6, -1128861017. / 918035200) == Catch::Approx(8));
+ }
+
+ SECTION("degree 7")
+ {
+ PrismGaussQuadrature l7(7);
+
+ REQUIRE(l7.numberOfPoints() == 35);
+
+ REQUIRE(integrate(p0, l7) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l7) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l7) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l7) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l7) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l7) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l7) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l7) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l7) != Catch::Approx(-21178319419. / 27541056000));
+
+ REQUIRE(get_order(p8, l7, -21178319419. / 27541056000) == Catch::Approx(8));
+ }
+
+ SECTION("degree 8")
+ {
+ PrismGaussQuadrature l8(8);
+
+ REQUIRE(l8.numberOfPoints() == 46);
+
+ REQUIRE(integrate(p0, l8) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l8) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l8) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l8) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l8) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l8) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l8) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l8) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l8) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l8) != Catch::Approx(-5483758803191. / 9088548480000));
+
+ // In a weird way, one gets 10th order on this test
+ REQUIRE(get_order(p9, l8, -5483758803191. / 9088548480000) == Catch::Approx(10));
+ }
+
+ SECTION("degree 9")
+ {
+ PrismGaussQuadrature l9(9);
+
+ REQUIRE(l9.numberOfPoints() == 59);
+
+ REQUIRE(integrate(p0, l9) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l9) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l9) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l9) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l9) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l9) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l9) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l9) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l9) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l9) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l9) != Catch::Approx(-9456848221657. / 22721371200000));
+
+ REQUIRE(get_order(p10, l9, -9456848221657. / 22721371200000) == Catch::Approx(10));
+ }
+
+ SECTION("degree 10")
+ {
+ PrismGaussQuadrature l10(10);
+
+ REQUIRE(l10.numberOfPoints() == 84);
+
+ REQUIRE(integrate(p0, l10) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l10) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l10) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l10) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l10) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l10) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l10) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l10) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l10) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l10) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l10) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l10) != Catch::Approx(-4571362439539697. / 7518708288000000));
+
+ // In a weird way, one gets 12th order on this test
+ REQUIRE(get_order(p11, l10, -4571362439539697. / 7518708288000000) == Catch::Approx(12));
+ }
+
+ SECTION("degree 11")
+ {
+ PrismGaussQuadrature l11(11);
+
+ REQUIRE(l11.numberOfPoints() == 99);
+
+ REQUIRE(integrate(p0, l11) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l11) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l11) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l11) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l11) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l11) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l11) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l11) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l11) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l11) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l11) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l11) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l11) != Catch::Approx(-491755535075074133. / 1378429852800000000));
+
+ REQUIRE(get_order(p12, l11, -491755535075074133. / 1378429852800000000) == Catch::Approx(12));
+ }
+
+ SECTION("degree 12")
+ {
+ PrismGaussQuadrature l12(12);
+
+ REQUIRE(l12.numberOfPoints() == 136);
+
+ REQUIRE(integrate(p0, l12) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l12) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l12) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l12) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l12) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l12) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l12) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l12) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l12) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l12) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l12) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l12) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l12) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l12) != Catch::Approx(-1620413117251976393. / 4135289558400000000));
+
+ // In a weird way, one gets almost 14th order on this test
+ REQUIRE(get_order(p13, l12, -1620413117251976393. / 4135289558400000000) == Catch::Approx(14).margin(0.01));
+ }
+
+ SECTION("degree 13")
+ {
+ PrismGaussQuadrature l13(13);
+
+ REQUIRE(l13.numberOfPoints() == 162);
+
+ REQUIRE(integrate(p0, l13) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l13) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l13) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l13) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l13) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l13) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l13) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l13) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l13) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l13) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l13) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l13) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l13) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l13) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l13) != Catch::Approx(296918520496968826367. / 827057911680000000000.));
+
+ REQUIRE(get_order(p14, l13, 296918520496968826367. / 827057911680000000000.) == Catch::Approx(14));
+ }
+
+ SECTION("degree 14")
+ {
+ PrismGaussQuadrature l14(14);
+
+ REQUIRE(l14.numberOfPoints() == 194);
+
+ REQUIRE(integrate(p0, l14) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l14) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l14) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l14) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l14) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l14) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l14) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l14) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l14) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l14) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l14) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l14) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l14) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l14) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l14) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l14) != Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+
+ // In a weird way, one gets almost 16th order on this test
+ REQUIRE(get_order(p15, l14, -7727953154629829488841. / 210899767478400000000000.) ==
+ Catch::Approx(16).margin(0.01));
+ }
+
+ SECTION("degree 15")
+ {
+ PrismGaussQuadrature l15(15);
+
+ REQUIRE(l15.numberOfPoints() == 238);
+
+ REQUIRE(integrate(p0, l15) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l15) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l15) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l15) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l15) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l15) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l15) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l15) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l15) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l15) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l15) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l15) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l15) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l15) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l15) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l15) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l15) != Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+
+ REQUIRE(get_order(p16, l15, -18153283669186101815689. / 527249418696000000000000.) == Catch::Approx(16));
+ }
+
+ SECTION("degree 16")
+ {
+ PrismGaussQuadrature l16(16);
+
+ REQUIRE(l16.numberOfPoints() == 287);
+
+ REQUIRE(integrate(p0, l16) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l16) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l16) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l16) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l16) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l16) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l16) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l16) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l16) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l16) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l16) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l16) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l16) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l16) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l16) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l16) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l16) == Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+ REQUIRE(integrate(p17, l16) != Catch::Approx(5157361121064510230030071. / 200354779104480000000000000.));
+
+ // In a weird way, one gets almost 18th order on this test
+ REQUIRE(get_order(p17, l16, 5157361121064510230030071. / 200354779104480000000000000.) ==
+ Catch::Approx(18).margin(0.1));
+ }
+
+ SECTION("degree 17")
+ {
+ PrismGaussQuadrature l17(17);
+
+ REQUIRE(l17.numberOfPoints() == 338);
+
+ REQUIRE(integrate(p0, l17) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l17) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l17) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l17) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l17) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l17) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l17) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l17) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l17) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l17) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l17) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l17) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l17) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l17) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l17) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l17) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l17) == Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+ REQUIRE(integrate(p17, l17) == Catch::Approx(5157361121064510230030071. / 200354779104480000000000000.));
+ REQUIRE(integrate(p18, l17) !=
+ Catch::Approx(195337148397715128549413507. / 5609933814925440000000000000.).epsilon(1E-10));
+
+ REQUIRE(get_order(p18, l17, 195337148397715128549413507. / 5609933814925440000000000000.) == Catch::Approx(18));
+ }
+
+ SECTION("degree 18")
+ {
+ PrismGaussQuadrature l18(18);
+
+ REQUIRE(l18.numberOfPoints() == 396);
+
+ REQUIRE(integrate(p0, l18) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l18) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l18) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l18) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l18) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l18) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l18) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l18) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l18) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l18) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l18) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l18) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l18) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l18) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l18) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l18) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l18) == Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+ REQUIRE(integrate(p17, l18) == Catch::Approx(5157361121064510230030071. / 200354779104480000000000000.));
+ REQUIRE(integrate(p18, l18) == Catch::Approx(195337148397715128549413507. / 5609933814925440000000000000.));
+ REQUIRE(integrate(p19, l18) !=
+ Catch::Approx(-417563570921497136922189149. / 14024834537313600000000000000.).epsilon(1E-10));
+
+ // In a weird way, one gets almost 20th order on this test
+ REQUIRE(get_order(p19, l18, -417563570921497136922189149. / 14024834537313600000000000000.) ==
+ Catch::Approx(20).margin(0.01));
+ }
+
+ SECTION("degree 19")
+ {
+ PrismGaussQuadrature l19(19);
+
+ REQUIRE(l19.numberOfPoints() == 420);
+
+ REQUIRE(integrate(p0, l19) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l19) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l19) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l19) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l19) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l19) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l19) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l19) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l19) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l19) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l19) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l19) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l19) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l19) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l19) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l19) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l19) == Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+ REQUIRE(integrate(p17, l19) == Catch::Approx(5157361121064510230030071. / 200354779104480000000000000.));
+ REQUIRE(integrate(p18, l19) == Catch::Approx(195337148397715128549413507. / 5609933814925440000000000000.));
+ REQUIRE(integrate(p19, l19) == Catch::Approx(-417563570921497136922189149. / 14024834537313600000000000000.));
+ REQUIRE(integrate(p20, l19) !=
+ Catch::Approx(4740816174053415637444760963. / 205697573213932800000000000000.).epsilon(1E-10));
+
+ REQUIRE(get_order(p20, l19, 4740816174053415637444760963. / 205697573213932800000000000000.) ==
+ Catch::Approx(20).margin(0.01));
+ }
+
+ SECTION("degree 20")
+ {
+ PrismGaussQuadrature l20(20);
+
+ REQUIRE(l20.numberOfPoints() == 518);
+
+ REQUIRE(integrate(p0, l20) == Catch::Approx(4));
+ REQUIRE(integrate(p1, l20) == Catch::Approx(2. / 3));
+ REQUIRE(integrate(p2, l20) == Catch::Approx(47. / 24));
+ REQUIRE(integrate(p3, l20) == Catch::Approx(-427. / 360));
+ REQUIRE(integrate(p4, l20) == Catch::Approx(1003. / 3600));
+ REQUIRE(integrate(p5, l20) == Catch::Approx(34031. / 18000));
+ REQUIRE(integrate(p6, l20) == Catch::Approx(36346369. / 55082112));
+ REQUIRE(integrate(p7, l20) == Catch::Approx(-1128861017. / 918035200));
+ REQUIRE(integrate(p8, l20) == Catch::Approx(-21178319419. / 27541056000));
+ REQUIRE(integrate(p9, l20) == Catch::Approx(-5483758803191. / 9088548480000));
+ REQUIRE(integrate(p10, l20) == Catch::Approx(-9456848221657. / 22721371200000));
+ REQUIRE(integrate(p11, l20) == Catch::Approx(-4571362439539697. / 7518708288000000));
+ REQUIRE(integrate(p12, l20) == Catch::Approx(-491755535075074133. / 1378429852800000000));
+ REQUIRE(integrate(p13, l20) == Catch::Approx(-1620413117251976393. / 4135289558400000000));
+ REQUIRE(integrate(p14, l20) == Catch::Approx(296918520496968826367. / 827057911680000000000.));
+ REQUIRE(integrate(p15, l20) == Catch::Approx(-7727953154629829488841. / 210899767478400000000000.));
+ REQUIRE(integrate(p16, l20) == Catch::Approx(-18153283669186101815689. / 527249418696000000000000.));
+ REQUIRE(integrate(p17, l20) == Catch::Approx(5157361121064510230030071. / 200354779104480000000000000.));
+ REQUIRE(integrate(p18, l20) == Catch::Approx(195337148397715128549413507. / 5609933814925440000000000000.));
+ REQUIRE(integrate(p19, l20) == Catch::Approx(-417563570921497136922189149. / 14024834537313600000000000000.));
+ REQUIRE(integrate(p20, l20) == Catch::Approx(4740816174053415637444760963. / 205697573213932800000000000000.));
+ REQUIRE(integrate(p21, l20) !=
+ Catch::Approx(-164372186128198750911065614811351. / 7096566275880681600000000000000000.).epsilon(1E-10));
+
+ // In a weird way, one gets almost 22th order on this test
+ REQUIRE(get_order(p21, l20, -164372186128198750911065614811351. / 7096566275880681600000000000000000.) ==
+ Catch::Approx(22).margin(0.01));
+ }
+
+ SECTION("not implemented formulae")
+ {
+ REQUIRE_THROWS_WITH(PrismGaussQuadrature(21), "error: Gauss quadrature formulae handle degrees up to 20 on prisms");
+ }
+}