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

diff --git a/doc/news/changes/minor/20210729DavidWells b/doc/news/changes/minor/20210729DavidWells
new file mode 100644 (file)
index 0000000..a6189ca
--- /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.
+<br>
+(David Wells, 2021/07/29)

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index fe42e8724cca30f639ef6a126b8c455a445b47a9..6832c428a17e6e447e4eb563cbcb2ed02178f2d6 100644 (file)
--- 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 <int dim>
@@ -852,10 +862,13 @@ class QWitherdenVincentSimplex : public QSimplex<dim>
 {
 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 e617604e4fa22a524eac0107909e85fdeeabfefd..368cd980feee46a4e7afe66c741f415bea794854 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1529,7 +1529,8 @@ namespace
 
 template <int dim>
 QWitherdenVincentSimplex<dim>::QWitherdenVincentSimplex(
-  const unsigned int n_points_1D)
+  const unsigned int n_points_1D,
+  const bool         use_odd_order)
   : QSimplex<dim>(Quadrature<dim>())
 {
   Assert(1 <= dim && dim <= 3, ExcNotImplemented());
@@ -1553,6 +1554,7 @@ QWitherdenVincentSimplex<dim>::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<double, dim + 1> b_point;
@@ -1564,6 +1566,7 @@ QWitherdenVincentSimplex<dim>::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<dim>::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<double, dim + 1> b_point;
@@ -1592,128 +1596,386 @@ QWitherdenVincentSimplex<dim>::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<dim>::operator=(QWitherdenVincentSimplex<dim>(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 3a82f4ba43c85664cbc5b8fcd34d70282a552d46..2611b1a8b7e6d13fe48190a4b5c17162e42952dd 100644 (file)
--- 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 <int dim>
 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<dim>      func(monomial_powers);
-  const QWitherdenVincentSimplex<dim> quad(n_points_1D);
+  const QWitherdenVincentSimplex<dim> 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 413dac5e4c1790e596044e16cb5b77df22379b80..a05c3ed85acf792a399207e97d37647c24e70616 100644 (file)
--- 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