From: David Wells Date: Thu, 4 Feb 2021 21:12:37 +0000 (-0500) Subject: Add Simplex::QWitherdenVincent. X-Git-Tag: v9.3.0-rc1~481^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6d8c12bfff5e8185ee8bb66a1b461601704da46b;p=dealii.git Add Simplex::QWitherdenVincent. This simplex quadrature formula has a clean implementation in terms of symmetry and permits construction of higher-order rules. --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index bf632316a6..2694ad09b9 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1089,3 +1089,22 @@ year = {2008}, 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} +} diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h index 7962091ca5..93f351f54a 100644 --- a/include/deal.II/simplex/quadrature_lib.h +++ b/include/deal.II/simplex/quadrature_lib.h @@ -51,6 +51,36 @@ namespace Simplex 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 + class QWitherdenVincent : public QSimplex + { + 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. */ diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc index 55d8c9eb69..624c3e8125 100644 --- a/source/simplex/quadrature_lib.cc +++ b/source/simplex/quadrature_lib.cc @@ -263,6 +263,247 @@ namespace Simplex Utilities::to_string(n_points_1D))); } + namespace + { + template + std::vector> + all_permutations(const std::array &b_point) + { + std::vector> 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 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 + QWitherdenVincent::QWitherdenVincent(const unsigned int n_points_1D) + : QSimplex(Quadrature()) + { + 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::operator=(dealii::QGauss(n_points_1D)); + return; + } + + std::array centroid; + std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0)); + std::vector>> b_point_permutations; + std::vector 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 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 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 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 &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 c_point; + std::copy(b_point.begin(), + b_point.begin() + dim, + c_point.begin_raw()); + this->quadrature_points.emplace_back(c_point); + } + } + } + template @@ -345,4 +586,8 @@ template class Simplex::QGaussPyramid<1>; 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 diff --git a/tests/simplex/q_witherden_vincent_01.cc b/tests/simplex/q_witherden_vincent_01.cc new file mode 100644 index 0000000000..67010b0a31 --- /dev/null +++ b/tests/simplex/q_witherden_vincent_01.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#include "../tests.h" + +template +void +print(const unsigned int n_points_1D) +{ + deallog << "n_points_1D = " << n_points_1D << std::endl; + const Simplex::QWitherdenVincent 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 +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 func(monomial_powers); + const Simplex::QWitherdenVincent 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); +} diff --git a/tests/simplex/q_witherden_vincent_01.with_simplex_support=on.output b/tests/simplex/q_witherden_vincent_01.with_simplex_support=on.output new file mode 100644 index 0000000000..413dac5e4c --- /dev/null +++ b/tests/simplex/q_witherden_vincent_01.with_simplex_support=on.output @@ -0,0 +1,131 @@ + +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