month={Jan},
pages={A2830–A2857}
}
+
+@article{witherden2015identification,
+ title={On the identification of symmetric quadrature rules for finite element methods},
+ author={Witherden, Freddie D and Vincent, Peter E},
+ journal={Computers \& Mathematics with Applications},
+ volume={69},
+ number={10},
+ pages={1232--1241},
+ year={2015},
+ publisher={Elsevier}
+}
+
+@online{quadpy,
+ author = {Schl{\"o}mer, Nico},
+ title = {quadpy: Your one-stop shop for numerical integration in Python},
+ year = 2021,
+ url = {https://github.com/nschloe/quadpy/},
+ urldate = {2021-02-08}
+}
explicit QGauss(const unsigned int n_points_1D);
};
+ /**
+ * Witherden-Vincent rules for simplex entities.
+ *
+ * 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
+ *
+ * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`.
+ *
+ * These rules match the ones listed for Witherden-Vincent in the quadpy
+ * @cite quadpy library and were first described in
+ * @cite witherden2015identification.
+ *
+ * @ingroup simplex
+ */
+ template <int dim>
+ class QWitherdenVincent : public QSimplex<dim>
+ {
+ public:
+ /**
+ * Constructor taking the number of quadrature points in 1D direction
+ * @p n_points_1D.
+ */
+ explicit QWitherdenVincent(const unsigned int n_points_1D);
+ };
+
/**
* Integration rule for wedge entities.
*/
Utilities::to_string(n_points_1D)));
}
+ namespace
+ {
+ template <std::size_t b_dim>
+ std::vector<std::array<double, b_dim>>
+ all_permutations(const std::array<double, b_dim> &b_point)
+ {
+ std::vector<std::array<double, b_dim>> output;
+
+ // We want all possible permutations of the barycentric coordinates.
+ // The easiest way to get all of them is to sort the input first and
+ // then use next_permutation to cycle through them all.
+ std::array<double, b_dim> temp = b_point;
+ std::sort(temp.begin(), temp.end());
+ do
+ {
+ output.push_back(temp);
+ }
+ while (std::next_permutation(temp.begin(), temp.end()));
+
+ return output;
+ }
+ } // namespace
+
+
+
+ template <int dim>
+ QWitherdenVincent<dim>::QWitherdenVincent(const unsigned int n_points_1D)
+ : QSimplex<dim>(Quadrature<dim>())
+ {
+ Assert(1 <= dim && dim <= 3, ExcNotImplemented());
+ // Just use Gauss in 1D: this is a high-order open rule so this is a
+ // reasonable equivalent for generic programming.
+ if (dim == 1)
+ {
+ Quadrature<dim>::operator=(dealii::QGauss<dim>(n_points_1D));
+ return;
+ }
+
+ std::array<double, dim + 1> centroid;
+ std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
+ std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
+ std::vector<double> b_weights;
+
+ // We can simplify the implementation of these quadrature rules
+ // by quite a bit by exploiting symmetry - we do essentially the
+ // same thing for each barycentric coordinate, so we can express
+ // our quadrature rule as permutations of barycentric points
+ // instead of writing things out explicitly.
+
+ // Apply a Barycentric permutation where one point is different.
+ auto process_point_1 = [&](const double a, const double w) {
+ const double b = 1.0 - dim * a;
+ std::array<double, dim + 1> b_point;
+ std::fill(b_point.begin(), b_point.begin() + dim, a);
+ b_point[dim] = b;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ // Apply a Barycentric permutation where two points (in 3D) are different.
+ auto process_point_2 = [&](const double a, const double w) {
+ Assert(dim == 3, ExcInternalError());
+ const double b = (1.0 - 2.0 * a) / 2.0;
+ std::array<double, dim + 1> b_point;
+ std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
+ b_point[dim - 1] = b;
+ b_point[dim] = b;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ // Apply a Barycentric permutation where three (or four) points
+ // are different (since there are two inputs).
+ 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;
+ std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
+ b_point[dim - 1] = b;
+ b_point[dim] = c;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ if (n_points_1D == 1)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(1.0);
+ }
+ else if (n_points_1D == 2)
+ {
+ // This is WV-4 in 2D and WV-3 in 3D
+ if (dim == 2)
+ {
+ process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
+ process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
+ }
+ else if (dim == 3)
+ {
+ process_point_1(3.281633025163817e-01, 1.362178425370874e-01);
+ process_point_1(1.080472498984286e-01, 1.137821574629126e-01);
+ }
+ }
+ else if (n_points_1D == 3)
+ {
+ // This is the WV-5 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ 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);
+ }
+ }
+ else if (n_points_1D == 4)
+ {
+ // This is the WV-7 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ 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);
+ }
+ }
+ else if (n_points_1D == 5)
+ {
+ // This is the WV-9 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ 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);
+ }
+ }
+ else if (n_points_1D == 6)
+ {
+ // There is no WV-11 rule in 3D yet
+ if (dim == 2)
+ {
+ 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 if (dim == 3)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
+ for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
+ ++permutation_n)
+ {
+ for (const std::array<double, dim + 1> &b_point :
+ b_point_permutations[permutation_n])
+ {
+ const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
+ this->weights.emplace_back(volume * b_weights[permutation_n]);
+ Point<dim> c_point;
+ std::copy(b_point.begin(),
+ b_point.begin() + dim,
+ c_point.begin_raw());
+ this->quadrature_points.emplace_back(c_point);
+ }
+ }
+ }
+
template <int dim>
template class Simplex::QGaussPyramid<2>;
template class Simplex::QGaussPyramid<3>;
+template class Simplex::QWitherdenVincent<1>;
+template class Simplex::QWitherdenVincent<2>;
+template class Simplex::QWitherdenVincent<3>;
+
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+print(const unsigned int n_points_1D)
+{
+ deallog << "n_points_1D = " << n_points_1D << std::endl;
+ const Simplex::QWitherdenVincent<dim> quad(n_points_1D);
+
+ deallog << "quad size = " << quad.size() << std::endl;
+ for (unsigned int q = 0; q < quad.size(); ++q)
+ {
+ deallog << quad.point(q) << " ";
+ deallog << quad.weight(q) << " ";
+ deallog << std::endl;
+ }
+}
+
+template <int dim>
+void
+check_accuracy_1D(const unsigned int n_points_1D)
+{
+ const unsigned int accuracy = 2 * n_points_1D - 1;
+
+ Tensor<1, dim> monomial_powers;
+ unsigned int sum = 0;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ monomial_powers[d] += accuracy / dim;
+ sum += accuracy / dim;
+ }
+
+ // if we aren't at the correct degree then add the rest to the final
+ // component
+ monomial_powers[dim - 1] += accuracy - sum;
+
+ const Functions::Monomial<dim> func(monomial_powers);
+ const Simplex::QWitherdenVincent<dim> quad(n_points_1D);
+
+ deallog << "Monomial powers = " << monomial_powers << std::endl;
+ double integrand = 0.0;
+ for (unsigned int q = 0; q < quad.size(); ++q)
+ integrand += quad.weight(q) * func.value(quad.point(q));
+ auto old_precision = deallog.precision(16);
+ deallog << "Integrand = " << integrand << std::endl;
+ deallog.precision(old_precision);
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog << "2D:" << std::endl;
+ print<2>(1);
+ print<2>(2);
+ print<2>(3);
+ print<2>(4);
+ deallog << std::endl << "3D:" << std::endl;
+ print<3>(1);
+ print<3>(2);
+ print<3>(3);
+ print<3>(4);
+
+ deallog << std::endl << 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<3>(1);
+ check_accuracy_1D<3>(2);
+ check_accuracy_1D<3>(3);
+ check_accuracy_1D<3>(4);
+ check_accuracy_1D<3>(5);
+}
--- /dev/null
+
+DEAL::2D:
+DEAL::n_points_1D = 1
+DEAL::quad size = 1
+DEAL::0.333333 0.333333 0.500000
+DEAL::n_points_1D = 2
+DEAL::quad size = 6
+DEAL::0.0915762 0.0915762 0.0549759
+DEAL::0.0915762 0.816848 0.0549759
+DEAL::0.816848 0.0915762 0.0549759
+DEAL::0.108103 0.445948 0.111691
+DEAL::0.445948 0.108103 0.111691
+DEAL::0.445948 0.445948 0.111691
+DEAL::n_points_1D = 3
+DEAL::quad size = 7
+DEAL::0.333333 0.333333 0.112500
+DEAL::0.101287 0.101287 0.0629696
+DEAL::0.101287 0.797427 0.0629696
+DEAL::0.797427 0.101287 0.0629696
+DEAL::0.0597159 0.470142 0.0661971
+DEAL::0.470142 0.0597159 0.0661971
+DEAL::0.470142 0.470142 0.0661971
+DEAL::n_points_1D = 4
+DEAL::quad size = 15
+DEAL::0.0337306 0.0337306 0.00827253
+DEAL::0.0337306 0.932539 0.00827253
+DEAL::0.932539 0.0337306 0.00827253
+DEAL::0.0513806 0.474310 0.0385433
+DEAL::0.474310 0.0513806 0.0385433
+DEAL::0.474310 0.474310 0.0385433
+DEAL::0.241577 0.241577 0.0639721
+DEAL::0.241577 0.516845 0.0639721
+DEAL::0.516845 0.241577 0.0639721
+DEAL::0.0470366 0.198683 0.0279394
+DEAL::0.0470366 0.754280 0.0279394
+DEAL::0.198683 0.0470366 0.0279394
+DEAL::0.198683 0.754280 0.0279394
+DEAL::0.754280 0.0470366 0.0279394
+DEAL::0.754280 0.198683 0.0279394
+DEAL::
+DEAL::3D:
+DEAL::n_points_1D = 1
+DEAL::quad size = 1
+DEAL::0.250000 0.250000 0.250000 0.166667
+DEAL::n_points_1D = 2
+DEAL::quad size = 8
+DEAL::0.0155101 0.328163 0.328163 0.0227030
+DEAL::0.328163 0.0155101 0.328163 0.0227030
+DEAL::0.328163 0.328163 0.0155101 0.0227030
+DEAL::0.328163 0.328163 0.328163 0.0227030
+DEAL::0.108047 0.108047 0.108047 0.0189637
+DEAL::0.108047 0.108047 0.675858 0.0189637
+DEAL::0.108047 0.675858 0.108047 0.0189637
+DEAL::0.675858 0.108047 0.108047 0.0189637
+DEAL::n_points_1D = 3
+DEAL::quad size = 14
+DEAL::0.0673422 0.310886 0.310886 0.0187813
+DEAL::0.310886 0.0673422 0.310886 0.0187813
+DEAL::0.310886 0.310886 0.0673422 0.0187813
+DEAL::0.310886 0.310886 0.310886 0.0187813
+DEAL::0.0927353 0.0927353 0.0927353 0.0122488
+DEAL::0.0927353 0.0927353 0.721794 0.0122488
+DEAL::0.0927353 0.721794 0.0927353 0.0122488
+DEAL::0.721794 0.0927353 0.0927353 0.0122488
+DEAL::0.0455037 0.0455037 0.454496 0.00709100
+DEAL::0.0455037 0.454496 0.0455037 0.00709100
+DEAL::0.0455037 0.454496 0.454496 0.00709100
+DEAL::0.454496 0.0455037 0.0455037 0.00709100
+DEAL::0.454496 0.0455037 0.454496 0.00709100
+DEAL::0.454496 0.454496 0.0455037 0.00709100
+DEAL::n_points_1D = 4
+DEAL::quad size = 35
+DEAL::0.250000 0.250000 0.250000 0.0159142
+DEAL::0.0528966 0.315701 0.315701 0.00705493
+DEAL::0.315701 0.0528966 0.315701 0.00705493
+DEAL::0.315701 0.315701 0.0528966 0.00705493
+DEAL::0.315701 0.315701 0.315701 0.00705493
+DEAL::0.0504898 0.0504898 0.449510 0.00531615
+DEAL::0.0504898 0.449510 0.0504898 0.00531615
+DEAL::0.0504898 0.449510 0.449510 0.00531615
+DEAL::0.449510 0.0504898 0.0504898 0.00531615
+DEAL::0.449510 0.0504898 0.449510 0.00531615
+DEAL::0.449510 0.449510 0.0504898 0.00531615
+DEAL::0.0471607 0.188834 0.188834 0.00620119
+DEAL::0.0471607 0.188834 0.575172 0.00620119
+DEAL::0.0471607 0.575172 0.188834 0.00620119
+DEAL::0.188834 0.0471607 0.188834 0.00620119
+DEAL::0.188834 0.0471607 0.575172 0.00620119
+DEAL::0.188834 0.188834 0.0471607 0.00620119
+DEAL::0.188834 0.188834 0.575172 0.00620119
+DEAL::0.188834 0.575172 0.0471607 0.00620119
+DEAL::0.188834 0.575172 0.188834 0.00620119
+DEAL::0.575172 0.0471607 0.188834 0.00620119
+DEAL::0.575172 0.188834 0.0471607 0.00620119
+DEAL::0.575172 0.188834 0.188834 0.00620119
+DEAL::0.0212655 0.0212655 0.146639 0.00135180
+DEAL::0.0212655 0.0212655 0.810830 0.00135180
+DEAL::0.0212655 0.146639 0.0212655 0.00135180
+DEAL::0.0212655 0.146639 0.810830 0.00135180
+DEAL::0.0212655 0.810830 0.0212655 0.00135180
+DEAL::0.0212655 0.810830 0.146639 0.00135180
+DEAL::0.146639 0.0212655 0.0212655 0.00135180
+DEAL::0.146639 0.0212655 0.810830 0.00135180
+DEAL::0.146639 0.810830 0.0212655 0.00135180
+DEAL::0.810830 0.0212655 0.0212655 0.00135180
+DEAL::0.810830 0.0212655 0.146639 0.00135180
+DEAL::0.810830 0.146639 0.0212655 0.00135180
+DEAL::
+DEAL::
+DEAL::Monomial powers = 0.00000 1.00000
+DEAL::Integrand = 0.1666666666666667
+DEAL::Monomial powers = 1.00000 2.00000
+DEAL::Integrand = 0.01666666666666667
+DEAL::Monomial powers = 2.00000 3.00000
+DEAL::Integrand = 0.002380952380952381
+DEAL::Monomial powers = 3.00000 4.00000
+DEAL::Integrand = 0.0003968253968253969
+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 = 0.00000 0.00000 1.00000
+DEAL::Integrand = 0.04166666666666666
+DEAL::Monomial powers = 1.00000 1.00000 1.00000
+DEAL::Integrand = 0.001388888888888889
+DEAL::Monomial powers = 1.00000 1.00000 3.00000
+DEAL::Integrand = 0.0001488095238095239
+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
projects / dealii.git / commitdiff