From a14703d7f66dd2f3eb64c5959222b95bfe39d3ab Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 29 Jul 2021 12:21:54 -0400 Subject: [PATCH] Add a lot more Witherden-Vincent rules. This includes most of the even-degree rules. For posterity: here's the script I used to print quadpy rules rather than copying numbers individually: import quadpy as qp def format_16(x): return '{:.16e}'.format(float(x)) def print_scheme(scheme): sym_data = scheme.symmetry_data print("// WV-{}, {}D".format(scheme.degree, scheme.dim)) for key in sym_data.keys(): if key == 'centroid' or key == 's4': print("b_point_permutations.push_back({centroid});") print("b_weights.push_back({});".format(format_16(sym_data[key][0][0]))) elif key == 'd3_aa' or key == 's31': weights = [format_16(w) for w in sym_data[key][0]] points = [format_16(p) for p in sym_data[key][1]] assert len(weights) == len(points) for i in range(len(weights)): print("process_point_1({}, {});".format(points[i], weights[i])) elif key == 'd3_ab' or key == 's211': weights = [format_16(w) for w in sym_data[key][0]] points1 = [format_16(p) for p in sym_data[key][1]] points2 = [format_16(p) for p in sym_data[key][2]] assert len(weights) == len(points1) assert len(points2) == len(points1) for i in range(len(weights)): print(("process_point_3({}," "\n {}," "\n {});") .format(points1[i], points2[i], weights[i])) elif key == 's22': assert scheme.dim == 3 weights = [format_16(w) for w in sym_data[key][0]] points = [format_16(p) for p in sym_data[key][1]] assert len(weights) == len(points) for i in range(len(weights)): print(("process_point_2({}," "\n {});") .format(points[i], weights[i])) else: assert False --- doc/news/changes/minor/20210729DavidWells | 4 + include/deal.II/base/quadrature_lib.h | 31 +- source/base/quadrature_lib.cc | 460 +++++++++++++++----- tests/simplex/q_witherden_vincent_01.cc | 25 +- tests/simplex/q_witherden_vincent_01.output | 28 ++ 5 files changed, 437 insertions(+), 111 deletions(-) create mode 100644 doc/news/changes/minor/20210729DavidWells diff --git a/doc/news/changes/minor/20210729DavidWells b/doc/news/changes/minor/20210729DavidWells new file mode 100644 index 0000000000..a6189ca950 --- /dev/null +++ b/doc/news/changes/minor/20210729DavidWells @@ -0,0 +1,4 @@ +New: QWitherdenVincentSimplex now implements even-order rules +in addition to the standard Gauss-like odd-order rules. +
+(David Wells, 2021/07/29) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index fe42e8724c..6832c428a1 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -832,19 +832,29 @@ public: * * Like QGauss, users should specify a number `n_points_1D` as an indication * of what polynomial degree to be integrated exactly (e.g., for $n$ points, - * the rule can integrate polynomials of degree $2 n - 1$ exactly). The given - * value for n_points_1D = 1, 2, 3, 4, 5 results in the following number of - * quadrature points in 2D and 3D: - * - 2D: 1, 6, 7, 15, 19 - * - 3D: 1, 8, 14, 35, 59 + * the rule can integrate polynomials of degree $2 n - 1$ exactly). + * Additionally, since these rules were derived for simplices, there are + * also even-ordered rules (i.e., they integrate polynomials of degree $2 n$) + * available which do not have analogous 1D rules. + * + * The given value for n_points_1D = 1, 2, 3, 4, 5, 6, 7 (where the last two are + * only implemented in 2D) results in the following number of quadrature points + * in 2D and 3D: + * - 2D: odd (default): 1, 6, 7, 15, 19, 28, 37 + * - 2D: even: 3, 6, 12, 16, 25, 33, 42 + * - 3D: odd (default): 1, 8, 14, 35, 59 + * - 3D: even: 4, 14, 24, 46, 81 * * For 1D, the quadrature rule degenerates to a - * `dealii::QGauss<1>(n_points_1D)`. + * `dealii::QGauss<1>(n_points_1D)` and @p use_odd_order is ignored. * * These rules match the ones listed for Witherden-Vincent in the quadpy * @cite quadpy library and were first described in * @cite witherden2015identification. * + * @note Some rules (2D 2 odd and 3D 2 even) do not yet exist and instead a + * higher-order rule is used in their place. + * * @ingroup simplex */ template @@ -852,10 +862,13 @@ class QWitherdenVincentSimplex : public QSimplex { public: /** - * Constructor taking the number of quadrature points in 1D direction - * @p n_points_1D. + * Constructor taking the equivalent number of quadrature points in 1D + * @p n_points_1D and boolean indicating whether the rule should be order + * $2 n - 1$ or $2 n$: see the general documentation of this class for more + * information. */ - explicit QWitherdenVincentSimplex(const unsigned int n_points_1D); + explicit QWitherdenVincentSimplex(const unsigned int n_points_1D, + const bool use_odd_order = true); }; /** diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index e617604e4f..368cd980fe 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1529,7 +1529,8 @@ namespace template QWitherdenVincentSimplex::QWitherdenVincentSimplex( - const unsigned int n_points_1D) + const unsigned int n_points_1D, + const bool use_odd_order) : QSimplex(Quadrature()) { Assert(1 <= dim && dim <= 3, ExcNotImplemented()); @@ -1553,6 +1554,7 @@ QWitherdenVincentSimplex::QWitherdenVincentSimplex( // instead of writing things out explicitly. // Apply a Barycentric permutation where one point is different. + // Equivalent to d3_aa and s31 in quadpy. auto process_point_1 = [&](const double a, const double w) { const double b = 1.0 - dim * a; std::array b_point; @@ -1564,6 +1566,7 @@ QWitherdenVincentSimplex::QWitherdenVincentSimplex( }; // Apply a Barycentric permutation where two points (in 3D) are different. + // Equivalent to s22 in quadpy. auto process_point_2 = [&](const double a, const double w) { Assert(dim == 3, ExcInternalError()); const double b = (1.0 - 2.0 * a) / 2.0; @@ -1578,6 +1581,7 @@ QWitherdenVincentSimplex::QWitherdenVincentSimplex( // Apply a Barycentric permutation where three (or four) points // are different (since there are two inputs). + // Equivalent to d3_ab and s211 in quadpy. auto process_point_3 = [&](const double a, const double b, const double w) { const double c = 1.0 - (dim - 1.0) * a - b; std::array b_point; @@ -1592,128 +1596,386 @@ QWitherdenVincentSimplex::QWitherdenVincentSimplex( switch (n_points_1D) { case 1: - b_point_permutations.push_back({centroid}); - b_weights.push_back(1.0); - break; - case 2: - // This is WV-4 in 2D and WV-3 in 3D - if (dim == 2) + switch (dim) { - process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01); - process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01); + case 2: + if (use_odd_order) + { + // WV-1, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(1.0000000000000000e+00); + } + else + { + // WV-2, 2D + process_point_1(1.6666666666666669e-01, + 3.3333333333333331e-01); + } + break; + case 3: + if (use_odd_order) + { + // WV-1, 3D + b_point_permutations.push_back({centroid}); + b_weights.push_back(1.0000000000000000e+00); + } + else + { + // WV-2, 3D + process_point_1(1.3819660112501050e-01, + 2.5000000000000000e-01); + } + break; + default: + Assert(false, ExcNotImplemented()); } - else if (dim == 3) + break; + case 2: + switch (dim) { - process_point_1(3.281633025163817e-01, 1.362178425370874e-01); - process_point_1(1.080472498984286e-01, 1.137821574629126e-01); + case 2: + // WV-4 in both cases (no WV-3 in 2D) + process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01); + process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01); + break; + case 3: + if (use_odd_order) + { + // WV-3, 3D + process_point_1(3.2816330251638171e-01, + 1.3621784253708741e-01); + process_point_1(1.0804724989842859e-01, + 1.1378215746291261e-01); + } + else + { + // WV-5 (no WV-4 in 3D) + Quadrature::operator=(QWitherdenVincentSimplex(3)); + } + break; + default: + Assert(false, ExcInternalError()); } break; case 3: - // This is the WV-5 rule in both 2D and 3D - if (dim == 2) + switch (dim) { - b_weights.push_back(0.225); - b_point_permutations.push_back({centroid}); - - process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01); - process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01); - } - else if (dim == 3) - { - process_point_1(3.108859192633006e-01, 1.126879257180159e-01); - process_point_1(9.273525031089125e-02, 7.349304311636196e-02); - - process_point_2(4.550370412564964e-02, 4.254602077708147e-02); + case 2: + if (use_odd_order) + { + // WV-5, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(2.2500000000000001e-01); + process_point_1(1.0128650732345634e-01, + 1.2593918054482714e-01); + process_point_1(4.7014206410511511e-01, + 1.3239415278850619e-01); + } + else + { + // WV-6, 2D + process_point_1(6.3089014491502227e-02, + 5.0844906370206819e-02); + process_point_1(2.4928674517091043e-01, + 1.1678627572637937e-01); + process_point_3(5.3145049844816938e-02, + 3.1035245103378439e-01, + 8.2851075618373571e-02); + } + break; + case 3: + if (use_odd_order) + { + // WV-5, 3D + process_point_1(3.1088591926330061e-01, + 1.1268792571801590e-01); + process_point_1(9.2735250310891248e-02, + 7.3493043116361956e-02); + process_point_2(4.5503704125649642e-02, + 4.2546020777081472e-02); + } + else + { + // WV-6, 3D + process_point_1(4.0673958534611372e-02, + 1.0077211055320640e-02); + process_point_1(3.2233789014227548e-01, + 5.5357181543654717e-02); + process_point_1(2.1460287125915201e-01, + 3.9922750258167487e-02); + process_point_3(6.3661001875017442e-02, + 6.0300566479164919e-01, + 4.8214285714285710e-02); + } + break; + default: + Assert(false, ExcInternalError()); } break; case 4: - // This is the WV-7 rule in both 2D and 3D - if (dim == 2) + switch (dim) { - process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02); - process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02); - process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01); - process_point_3(4.7036644652595216e-02, - 1.9868331479735168e-01, - 5.5878732903199779e-02); - } - else if (dim == 3) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(9.548528946413085e-02); - - process_point_1(3.157011497782028e-01, 4.232958120996703e-02); - process_point_2(5.048982259839635e-02, 3.189692783285758e-02); - - process_point_3(1.888338310260010e-01, - 5.751716375870000e-01, - 3.720713072833462e-02); - process_point_3(2.126547254148314e-02, - 8.108302410985486e-01, - 8.110770829903342e-03); + case 2: + if (use_odd_order) + { + // WV-7, 2D + process_point_1(3.3730648554587850e-02, + 1.6545050110792131e-02); + process_point_1(4.7430969250471822e-01, + 7.7086646185986069e-02); + process_point_1(2.4157738259540357e-01, + 1.2794417123015558e-01); + process_point_3(4.7036644652595216e-02, + 1.9868331479735168e-01, + 5.5878732903199779e-02); + } + else + { + // WV-8, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(1.4431560767778717e-01); + process_point_1(5.0547228317030957e-02, + 3.2458497623198079e-02); + process_point_1(4.5929258829272313e-01, + 9.5091634267284619e-02); + process_point_1(1.7056930775176021e-01, + 1.0321737053471824e-01); + process_point_3(8.3947774099575878e-03, + 2.6311282963463811e-01, + 2.7230314174434993e-02); + } + break; + case 3: + if (use_odd_order) + { + // WV-7, 3D + b_point_permutations.push_back({centroid}); + b_weights.push_back(9.5485289464130846e-02); + process_point_1(3.1570114977820279e-01, + 4.2329581209967028e-02); + process_point_2(5.0489822598396350e-02, + 3.1896927832857580e-02); + process_point_3(1.8883383102600099e-01, + 5.7517163758699996e-01, + 3.7207130728334620e-02); + process_point_3(2.1265472541483140e-02, + 8.1083024109854862e-01, + 8.1107708299033420e-03); + } + else + { + // WV-8, 3D + process_point_1(1.0795272496221089e-01, + 2.6426650908408830e-02); + process_point_1(1.8510948778258660e-01, + 5.2031747563738531e-02); + process_point_1(4.2316543684767283e-02, + 7.5252561535401989e-03); + process_point_1(3.1418170912403898e-01, + 4.1763782856934897e-02); + process_point_2(4.3559132858383021e-01, + 3.6280930261308818e-02); + process_point_3(2.1433930127130570e-02, + 7.1746406342630831e-01, + 7.1569028908444327e-03); + process_point_3(2.0413933387602909e-01, + 5.8379737830214440e-01, + 1.5453486150960340e-02); + } + break; + default: + Assert(false, ExcInternalError()); } break; case 5: - // This is the WV-9 rule in both 2D and 3D - if (dim == 2) + switch (dim) { - b_point_permutations.push_back({centroid}); - b_weights.push_back(9.7135796282798836e-02); - - process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02); - process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02); - process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02); - process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02); - - process_point_3(3.6838412054736258e-02, - 2.2196298916076568e-01, - 4.3283539377289376e-02); - } - else if (dim == 3) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(5.801054891248025e-02); - - process_point_1(6.198169755222693e-10, 6.431928175925639e-05); - process_point_1(1.607745353952616e-01, 2.317333846242546e-02); - process_point_1(3.222765218214210e-01, 2.956291233542929e-02); - process_point_1(4.510891834541358e-02, 8.063979979616182e-03); - - process_point_2(1.122965460043761e-01, 3.813408010370246e-02); - - process_point_3(4.588714487524592e-01, - 2.554579233041310e-03, - 8.384422198298552e-03); - process_point_3(3.377587068533860e-02, - 7.183503264420745e-01, - 1.023455935274533e-02); - process_point_3(1.836413698099279e-01, - 3.441591057817528e-02, - 2.052491596798814e-02); + case 2: + if (use_odd_order) + { + // WV-9, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(9.7135796282798836e-02); + process_point_1(4.4729513394452691e-02, + 2.5577675658698031e-02); + process_point_1(4.8968251919873762e-01, + 3.1334700227139071e-02); + process_point_1(4.3708959149293664e-01, + 7.7827541004774278e-02); + process_point_1(1.8820353561903275e-01, + 7.9647738927210249e-02); + process_point_3(3.6838412054736258e-02, + 2.2196298916076568e-01, + 4.3283539377289376e-02); + } + else + { + // WV-10, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(8.1743329146285973e-02); + process_point_1(3.2055373216943517e-02, + 1.3352968813149567e-02); + process_point_1(1.4216110105656438e-01, + 4.5957963604744731e-02); + process_point_3(2.8367665339938453e-02, + 1.6370173373718250e-01, + 2.5297757707288385e-02); + process_point_3(2.9619889488729734e-02, + 3.6914678182781102e-01, + 3.4184648162959429e-02); + process_point_3(1.4813288578382056e-01, + 3.2181299528883545e-01, + 6.3904906396424044e-02); + } + break; + case 3: + if (use_odd_order) + { + // WV-9, 3D + b_point_permutations.push_back({centroid}); + b_weights.push_back(5.8010548912480253e-02); + process_point_1(6.1981697552226933e-10, + 6.4319281759256394e-05); + process_point_1(1.6077453539526160e-01, + 2.3173338462425461e-02); + process_point_1(3.2227652182142102e-01, + 2.9562912335429289e-02); + process_point_1(4.5108918345413578e-02, + 8.0639799796161822e-03); + process_point_2(1.1229654600437609e-01, + 3.8134080103702457e-02); + process_point_3(4.5887144875245922e-01, + 2.5545792330413102e-03, + 8.3844221982985519e-03); + process_point_3(3.3775870685338598e-02, + 7.1835032644207453e-01, + 1.0234559352745330e-02); + process_point_3(1.8364136980992790e-01, + 3.4415910578175279e-02, + 2.0524915967988139e-02); + } + else + { + // WV-10, 3D + b_point_permutations.push_back({centroid}); + b_weights.push_back(4.7399773556020743e-02); + process_point_1(3.1225006869518868e-01, + 2.6937059992268701e-02); + process_point_1(1.1430965385734609e-01, + 9.8691597167933822e-03); + process_point_3(4.1043073921896539e-01, + 1.6548602561961109e-01, + 1.1393881220195230e-02); + process_point_3(6.1380088247906528e-03, + 9.4298876734520487e-01, + 3.6194434433925362e-04); + process_point_3(1.2105018114558939e-01, + 4.7719037990428043e-01, + 2.5739731980456069e-02); + process_point_3(3.2779468216442620e-02, + 5.9425626948000698e-01, + 1.0135871679755789e-02); + process_point_3(3.2485281564823047e-02, + 8.0117728465834437e-01, + 6.5761472770359038e-03); + process_point_3(1.7497934218393901e-01, + 6.2807184547536599e-01, + 1.2907035798861989e-02); + } + break; + default: + Assert(false, ExcNotImplemented()); } break; case 6: // There is no WV-11 rule in 3D yet Assert(dim == 2, ExcNotImplemented()); - b_point_permutations.push_back({centroid}); - b_weights.push_back(8.5761179732224219e-02); - - process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02); - process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02); - process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02); - process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02); - process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02); - - process_point_3(7.3254276860644785e-03, - 1.4932478865208237e-01, - 1.0290289572953278e-02); - process_point_3(4.6010500165429957e-02, - 2.8958112563770588e-01, - 4.0332476640500554e-02); + if (use_odd_order) + { + // WV-11, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(8.5761179732224219e-02); + process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02); + process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02); + process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02); + process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02); + process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02); + process_point_3(7.3254276860644785e-03, + 1.4932478865208237e-01, + 1.0290289572953278e-02); + process_point_3(4.6010500165429957e-02, + 2.8958112563770588e-01, + 4.0332476640500554e-02); + } + else + { + // WV-12, 2D + process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03); + process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02); + process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02); + process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02); + process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02); + process_point_3(2.1382490256170616e-02, + 1.2727971723358933e-01, + 1.5083677576511438e-02); + process_point_3(2.3034156355267121e-02, + 2.9165567973834094e-01, + 2.1783585038607559e-02); + process_point_3(1.1629601967792658e-01, + 2.5545422863851736e-01, + 4.3227363659414209e-02); + } break; case 7: + // There is no WV-13 rule in 3D yet Assert(dim == 2, ExcNotImplemented()); + if (use_odd_order) + { + // WV-13, 2D + b_point_permutations.push_back({centroid}); + b_weights.push_back(6.7960036586831640e-02); + process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03); + process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02); + process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02); + process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02); + process_point_3(5.1263891023823893e-03, + 2.7251581777342970e-01, + 9.5906810035432631e-03); + process_point_3(2.4370186901093827e-02, + 1.1092204280346341e-01, + 1.4965401105165668e-02); + process_point_3(8.7895483032197297e-02, + 1.6359740106785048e-01, + 2.4179039811593819e-02); + process_point_3(6.8012243554206653e-02, + 3.0844176089211778e-01, + 3.4641276140848373e-02); + } + else + { + // WV-14, 2D + process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03); + process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02); + process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02); + process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02); + process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02); + process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02); + process_point_3(1.2683309328720416e-03, + 1.1897449769695684e-01, + 5.0102288385006719e-03); + process_point_3(1.4646950055654417e-02, + 2.9837288213625779e-01, + 1.4436308113533840e-02); + process_point_3(5.7124757403647919e-02, + 1.7226668782135557e-01, + 2.4665753212563674e-02); + process_point_3(9.2916249356971847e-02, + 3.3686145979634496e-01, + 3.8571510787060684e-02); + } break; + default: + Assert(false, ExcNotImplemented()); } Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError()); diff --git a/tests/simplex/q_witherden_vincent_01.cc b/tests/simplex/q_witherden_vincent_01.cc index 3a82f4ba43..2611b1a8b7 100644 --- a/tests/simplex/q_witherden_vincent_01.cc +++ b/tests/simplex/q_witherden_vincent_01.cc @@ -36,9 +36,11 @@ print(const unsigned int n_points_1D) template void -check_accuracy_1D(const unsigned int n_points_1D) +check_accuracy_1D(const unsigned int n_points_1D, + const bool use_odd_order = true) { - const unsigned int accuracy = 2 * n_points_1D - 1; + const unsigned int accuracy = + use_odd_order ? 2 * n_points_1D - 1 : 2 * n_points_1D; Tensor<1, dim> monomial_powers; unsigned int sum = 0; @@ -53,7 +55,7 @@ check_accuracy_1D(const unsigned int n_points_1D) monomial_powers[dim - 1] += accuracy - sum; const Functions::Monomial func(monomial_powers); - const QWitherdenVincentSimplex quad(n_points_1D); + const QWitherdenVincentSimplex quad(n_points_1D, use_odd_order); deallog << "Monomial powers = " << monomial_powers << std::endl; double integrand = 0.0; @@ -81,16 +83,33 @@ main() print<3>(4); deallog << std::endl << std::endl; + deallog << "check odd orders" << std::endl; check_accuracy_1D<2>(1); check_accuracy_1D<2>(2); check_accuracy_1D<2>(3); check_accuracy_1D<2>(4); check_accuracy_1D<2>(5); check_accuracy_1D<2>(6); + check_accuracy_1D<2>(7); check_accuracy_1D<3>(1); check_accuracy_1D<3>(2); check_accuracy_1D<3>(3); check_accuracy_1D<3>(4); check_accuracy_1D<3>(5); + + deallog << "check even orders" << std::endl; + check_accuracy_1D<2>(1, false); + check_accuracy_1D<2>(2, false); + check_accuracy_1D<2>(3, false); + check_accuracy_1D<2>(4, false); + check_accuracy_1D<2>(5, false); + check_accuracy_1D<2>(6, false); + check_accuracy_1D<2>(7, false); + + check_accuracy_1D<3>(1, false); + check_accuracy_1D<3>(2, false); + check_accuracy_1D<3>(3, false); + check_accuracy_1D<3>(4, false); + check_accuracy_1D<3>(5, false); } diff --git a/tests/simplex/q_witherden_vincent_01.output b/tests/simplex/q_witherden_vincent_01.output index 413dac5e4c..a05c3ed85a 100644 --- a/tests/simplex/q_witherden_vincent_01.output +++ b/tests/simplex/q_witherden_vincent_01.output @@ -107,6 +107,7 @@ DEAL::0.810830 0.0212655 0.146639 0.00135180 DEAL::0.810830 0.146639 0.0212655 0.00135180 DEAL:: DEAL:: +DEAL::check odd orders DEAL::Monomial powers = 0.00000 1.00000 DEAL::Integrand = 0.1666666666666667 DEAL::Monomial powers = 1.00000 2.00000 @@ -119,6 +120,8 @@ DEAL::Monomial powers = 4.00000 5.00000 DEAL::Integrand = 7.215007215007216e-05 DEAL::Monomial powers = 5.00000 6.00000 DEAL::Integrand = 1.387501387501388e-05 +DEAL::Monomial powers = 6.00000 7.00000 +DEAL::Integrand = 2.775002775002775e-06 DEAL::Monomial powers = 0.00000 0.00000 1.00000 DEAL::Integrand = 0.04166666666666666 DEAL::Monomial powers = 1.00000 1.00000 1.00000 @@ -129,3 +132,28 @@ DEAL::Monomial powers = 2.00000 2.00000 3.00000 DEAL::Integrand = 6.613756613756609e-06 DEAL::Monomial powers = 3.00000 3.00000 3.00000 DEAL::Integrand = 4.509379509379515e-07 +DEAL::check even orders +DEAL::Monomial powers = 1.00000 1.00000 +DEAL::Integrand = 0.04166666666666666 +DEAL::Monomial powers = 2.00000 2.00000 +DEAL::Integrand = 0.005555555555555556 +DEAL::Monomial powers = 3.00000 3.00000 +DEAL::Integrand = 0.0008928571428571427 +DEAL::Monomial powers = 4.00000 4.00000 +DEAL::Integrand = 0.0001587301587301587 +DEAL::Monomial powers = 5.00000 5.00000 +DEAL::Integrand = 3.006253006253007e-05 +DEAL::Monomial powers = 6.00000 6.00000 +DEAL::Integrand = 5.946434517863085e-06 +DEAL::Monomial powers = 7.00000 7.00000 +DEAL::Integrand = 1.214063714063714e-06 +DEAL::Monomial powers = 0.00000 0.00000 2.00000 +DEAL::Integrand = 0.01666666666666667 +DEAL::Monomial powers = 1.00000 1.00000 2.00000 +DEAL::Integrand = 0.0003968253968253969 +DEAL::Monomial powers = 2.00000 2.00000 2.00000 +DEAL::Integrand = 2.204585537918871e-05 +DEAL::Monomial powers = 2.00000 2.00000 4.00000 +DEAL::Integrand = 2.405002405002404e-06 +DEAL::Monomial powers = 3.00000 3.00000 4.00000 +DEAL::Integrand = 1.387501387501388e-07 -- 2.39.5