From 277affa69e8dc1a74dd8fc8ad26809c36a4f9a66 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 15 Jul 2021 19:42:00 -0400 Subject: [PATCH] Implement QIteratedSimplex. A limited version of QIterated for simplex meshes. --- doc/news/changes/minor/20210715DavidWells | 3 + include/deal.II/base/quadrature_lib.h | 15 +++ source/base/quadrature_lib.cc | 86 ++++++++++++ tests/simplex/q_iterated_simplex_01.cc | 99 ++++++++++++++ tests/simplex/q_iterated_simplex_01.output | 148 +++++++++++++++++++++ 5 files changed, 351 insertions(+) create mode 100644 doc/news/changes/minor/20210715DavidWells create mode 100644 tests/simplex/q_iterated_simplex_01.cc create mode 100644 tests/simplex/q_iterated_simplex_01.output diff --git a/doc/news/changes/minor/20210715DavidWells b/doc/news/changes/minor/20210715DavidWells new file mode 100644 index 0000000000..ae20692bcd --- /dev/null +++ b/doc/news/changes/minor/20210715DavidWells @@ -0,0 +1,3 @@ +New: Added a class QIteratedSimplex for building composite simplex quadrature rules. +
+(David Wells, 2021/07/15) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 0e17569e45..fe42e8724c 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -858,6 +858,21 @@ public: explicit QWitherdenVincentSimplex(const unsigned int n_points_1D); }; +/** + * Iterated quadrature for simplices. Since simplex cannot be described as + * tensor products the base quadrature has equal dimension. + * + * At the moment @p n_copies must be a power of 2 due to the complexity of + * subdividing a simplex. + */ +template +class QIteratedSimplex : public Quadrature +{ +public: + QIteratedSimplex(const Quadrature &base_quadrature, + const unsigned int n_copies); +}; + /** * Integration rule for wedge entities. */ diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 083018184c..dd213ecb89 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -17,6 +17,13 @@ #include #include +#include +#include + +#include +#include +#include + #include #include #include @@ -1739,6 +1746,81 @@ QWitherdenVincentSimplex::QWitherdenVincentSimplex( +namespace +{ + template + Quadrature + setup_qiterated_1D(const Quadrature &, const unsigned int) + { + Assert(false, ExcInternalError()); + return Quadrature(); + } + + + + Quadrature<1> + setup_qiterated_1D(const Quadrature<1> &base_quad, + const unsigned int n_copies) + { + return QIterated<1>(base_quad, n_copies); + } +} // namespace + + + +template +QIteratedSimplex::QIteratedSimplex(const Quadrature &base_quad, + const unsigned int n_copies) +{ + switch (dim) + { + case 1: + static_cast &>(*this) = + setup_qiterated_1D(base_quad, n_copies); + break; + case 2: + case 3: + { + const auto n_refinements = + static_cast(std::round(std::log2(n_copies))); + Assert((1u << n_refinements) == n_copies, + ExcMessage("The number of copies must be a power of 2.")); + Triangulation tria; + const auto reference_cell = ReferenceCells::get_simplex(); + GridGenerator::reference_cell(tria, reference_cell); + tria.refine_global(n_refinements); + const Mapping &mapping = + reference_cell.template get_default_linear_mapping(); + FE_Nothing fe(reference_cell); + + FEValues fe_values(mapping, + fe, + base_quad, + update_quadrature_points | update_JxW_values); + std::vector> points; + std::vector weights; + for (const auto &cell : tria.active_cell_iterators()) + { + fe_values.reinit(cell); + for (unsigned int qp = 0; qp < base_quad.size(); ++qp) + { + points.push_back(fe_values.quadrature_point(qp)); + weights.push_back(fe_values.JxW(qp)); + } + } + + static_cast &>(*this) = + Quadrature(points, weights); + + break; + } + default: + Assert(false, ExcNotImplemented()); + } +} + + + template QGaussWedge::QGaussWedge(const unsigned int n_points) : Quadrature() @@ -1850,6 +1932,10 @@ template class QSimplex<1>; template class QSimplex<2>; template class QSimplex<3>; +template class QIteratedSimplex<1>; +template class QIteratedSimplex<2>; +template class QIteratedSimplex<3>; + template class QSplit<1>; template class QSplit<2>; template class QSplit<3>; diff --git a/tests/simplex/q_iterated_simplex_01.cc b/tests/simplex/q_iterated_simplex_01.cc new file mode 100644 index 0000000000..e9904690be --- /dev/null +++ b/tests/simplex/q_iterated_simplex_01.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// 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" + +// Test QIteratedSimplex accuracy + +template +void +print(const Quadrature &quad) +{ + 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 n_copies) +{ + 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 QIteratedSimplex quad(QWitherdenVincentSimplex(n_points_1D), + n_copies); + + 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(QIteratedSimplex<2>(QWitherdenVincentSimplex<2>(2), 1)); + print(QIteratedSimplex<2>(QWitherdenVincentSimplex<2>(2), 2)); + deallog << std::endl << "3D:" << std::endl; + print(QIteratedSimplex<3>(QWitherdenVincentSimplex<3>(2), 1)); + print(QIteratedSimplex<3>(QWitherdenVincentSimplex<3>(2), 2)); + + deallog << std::endl << std::endl; + check_accuracy_1D<2>(1, 1); + check_accuracy_1D<2>(1, 2); + check_accuracy_1D<2>(1, 4); + check_accuracy_1D<2>(2, 1); + check_accuracy_1D<2>(2, 2); + check_accuracy_1D<2>(2, 4); + check_accuracy_1D<2>(6, 1); + check_accuracy_1D<2>(6, 2); + check_accuracy_1D<2>(6, 4); + + check_accuracy_1D<3>(1, 1); + check_accuracy_1D<3>(1, 2); + check_accuracy_1D<3>(1, 4); + check_accuracy_1D<3>(3, 1); + check_accuracy_1D<3>(3, 2); + check_accuracy_1D<3>(3, 4); + check_accuracy_1D<3>(5, 1); + check_accuracy_1D<3>(5, 2); + check_accuracy_1D<3>(5, 4); +} diff --git a/tests/simplex/q_iterated_simplex_01.output b/tests/simplex/q_iterated_simplex_01.output new file mode 100644 index 0000000000..98d806df7d --- /dev/null +++ b/tests/simplex/q_iterated_simplex_01.output @@ -0,0 +1,148 @@ + +DEAL::2D: +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::quad size = 24 +DEAL::0.0457881 0.0457881 0.0137440 +DEAL::0.0457881 0.408424 0.0137440 +DEAL::0.408424 0.0457881 0.0137440 +DEAL::0.0540515 0.222974 0.0279227 +DEAL::0.222974 0.0540515 0.0279227 +DEAL::0.222974 0.222974 0.0279227 +DEAL::0.545788 0.0457881 0.0137440 +DEAL::0.545788 0.408424 0.0137440 +DEAL::0.908424 0.0457881 0.0137440 +DEAL::0.554052 0.222974 0.0279227 +DEAL::0.722974 0.0540515 0.0279227 +DEAL::0.722974 0.222974 0.0279227 +DEAL::0.0457881 0.545788 0.0137440 +DEAL::0.0457881 0.908424 0.0137440 +DEAL::0.408424 0.545788 0.0137440 +DEAL::0.0540515 0.722974 0.0279227 +DEAL::0.222974 0.554052 0.0279227 +DEAL::0.222974 0.722974 0.0279227 +DEAL::0.454212 0.0915762 0.0137440 +DEAL::0.0915762 0.454212 0.0137440 +DEAL::0.454212 0.454212 0.0137440 +DEAL::0.277026 0.277026 0.0279227 +DEAL::0.445948 0.277026 0.0279227 +DEAL::0.277026 0.445948 0.0279227 +DEAL:: +DEAL::3D: +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::quad size = 64 +DEAL::0.00775505 0.164082 0.164082 0.00283787 +DEAL::0.164082 0.00775505 0.164082 0.00283787 +DEAL::0.164082 0.164082 0.00775505 0.00283787 +DEAL::0.164082 0.164082 0.164082 0.00283787 +DEAL::0.0540236 0.0540236 0.0540236 0.00237046 +DEAL::0.0540236 0.0540236 0.337929 0.00237046 +DEAL::0.0540236 0.337929 0.0540236 0.00237046 +DEAL::0.337929 0.0540236 0.0540236 0.00237046 +DEAL::0.507755 0.164082 0.164082 0.00283787 +DEAL::0.664082 0.00775505 0.164082 0.00283787 +DEAL::0.664082 0.164082 0.00775505 0.00283787 +DEAL::0.664082 0.164082 0.164082 0.00283787 +DEAL::0.554024 0.0540236 0.0540236 0.00237046 +DEAL::0.554024 0.0540236 0.337929 0.00237046 +DEAL::0.554024 0.337929 0.0540236 0.00237046 +DEAL::0.837929 0.0540236 0.0540236 0.00237046 +DEAL::0.00775505 0.664082 0.164082 0.00283787 +DEAL::0.164082 0.507755 0.164082 0.00283787 +DEAL::0.164082 0.664082 0.00775505 0.00283787 +DEAL::0.164082 0.664082 0.164082 0.00283787 +DEAL::0.0540236 0.554024 0.0540236 0.00237046 +DEAL::0.0540236 0.554024 0.337929 0.00237046 +DEAL::0.0540236 0.837929 0.0540236 0.00237046 +DEAL::0.337929 0.554024 0.0540236 0.00237046 +DEAL::0.00775505 0.164082 0.664082 0.00283787 +DEAL::0.164082 0.00775505 0.664082 0.00283787 +DEAL::0.164082 0.164082 0.507755 0.00283787 +DEAL::0.164082 0.164082 0.664082 0.00283787 +DEAL::0.0540236 0.0540236 0.554024 0.00237046 +DEAL::0.0540236 0.0540236 0.837929 0.00237046 +DEAL::0.0540236 0.337929 0.554024 0.00237046 +DEAL::0.337929 0.0540236 0.554024 0.00237046 +DEAL::0.335918 0.171837 0.164082 0.00283787 +DEAL::0.492245 0.171837 0.164082 0.00283787 +DEAL::0.335918 0.328163 0.00775505 0.00283787 +DEAL::0.335918 0.328163 0.164082 0.00283787 +DEAL::0.445976 0.108047 0.0540236 0.00237046 +DEAL::0.445976 0.108047 0.337929 0.00237046 +DEAL::0.162071 0.391953 0.0540236 0.00237046 +DEAL::0.445976 0.391953 0.0540236 0.00237046 +DEAL::0.328163 0.164082 0.171837 0.00283787 +DEAL::0.171837 0.164082 0.171837 0.00283787 +DEAL::0.328163 0.00775505 0.328163 0.00283787 +DEAL::0.171837 0.164082 0.328163 0.00283787 +DEAL::0.391953 0.0540236 0.108047 0.00237046 +DEAL::0.108047 0.337929 0.108047 0.00237046 +DEAL::0.391953 0.0540236 0.391953 0.00237046 +DEAL::0.108047 0.0540236 0.391953 0.00237046 +DEAL::0.164082 0.171837 0.335918 0.00283787 +DEAL::0.164082 0.328163 0.335918 0.00283787 +DEAL::0.00775505 0.328163 0.335918 0.00283787 +DEAL::0.164082 0.171837 0.492245 0.00283787 +DEAL::0.0540236 0.391953 0.162071 0.00237046 +DEAL::0.337929 0.108047 0.445976 0.00237046 +DEAL::0.0540236 0.108047 0.445976 0.00237046 +DEAL::0.0540236 0.391953 0.445976 0.00237046 +DEAL::0.171837 0.492245 0.171837 0.00283787 +DEAL::0.328163 0.335918 0.171837 0.00283787 +DEAL::0.328163 0.335918 0.328163 0.00283787 +DEAL::0.171837 0.335918 0.328163 0.00283787 +DEAL::0.391953 0.445976 0.108047 0.00237046 +DEAL::0.108047 0.445976 0.108047 0.00237046 +DEAL::0.108047 0.445976 0.391953 0.00237046 +DEAL::0.391953 0.162071 0.391953 0.00237046 +DEAL:: +DEAL:: +DEAL::Monomial powers = 0.00000 1.00000 +DEAL::Integrand = 0.1666666666666667 +DEAL::Monomial powers = 0.00000 1.00000 +DEAL::Integrand = 0.1666666666666667 +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 = 1.00000 2.00000 +DEAL::Integrand = 0.01666666666666667 +DEAL::Monomial powers = 1.00000 2.00000 +DEAL::Integrand = 0.01666666666666667 +DEAL::Monomial powers = 5.00000 6.00000 +DEAL::Integrand = 1.387501387501388e-05 +DEAL::Monomial powers = 5.00000 6.00000 +DEAL::Integrand = 1.387501387501388e-05 +DEAL::Monomial powers = 5.00000 6.00000 +DEAL::Integrand = 1.387501387501387e-05 +DEAL::Monomial powers = 0.00000 0.00000 1.00000 +DEAL::Integrand = 0.04166666666666666 +DEAL::Monomial powers = 0.00000 0.00000 1.00000 +DEAL::Integrand = 0.04166666666666666 +DEAL::Monomial powers = 0.00000 0.00000 1.00000 +DEAL::Integrand = 0.04166666666666667 +DEAL::Monomial powers = 1.00000 1.00000 3.00000 +DEAL::Integrand = 0.0001488095238095239 +DEAL::Monomial powers = 1.00000 1.00000 3.00000 +DEAL::Integrand = 0.0001488095238095238 +DEAL::Monomial powers = 1.00000 1.00000 3.00000 +DEAL::Integrand = 0.0001488095238095236 +DEAL::Monomial powers = 3.00000 3.00000 3.00000 +DEAL::Integrand = 4.509379509379515e-07 +DEAL::Monomial powers = 3.00000 3.00000 3.00000 +DEAL::Integrand = 4.509379509379510e-07 +DEAL::Monomial powers = 3.00000 3.00000 3.00000 +DEAL::Integrand = 4.509379509379522e-07 -- 2.39.5