Implement QIteratedSimplex.

A limited version of QIterated for simplex meshes.

diff --git a/doc/news/changes/minor/20210715DavidWells b/doc/news/changes/minor/20210715DavidWells
new file mode 100644 (file)
index 0000000..ae20692
--- /dev/null
+++ b/doc/news/changes/minor/20210715DavidWells
@@ -0,0 +1,3 @@
+New: Added a class QIteratedSimplex for building composite simplex quadrature rules.
+<br>
+(David Wells, 2021/07/15)

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index 0e17569e45d4acbcbb4237c672f7085e751400b4..fe42e8724cca30f639ef6a126b8c455a445b47a9 100644 (file)
--- 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 <int dim>
+class QIteratedSimplex : public Quadrature<dim>
+{
+public:
+  QIteratedSimplex(const Quadrature<dim> &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 083018184c3f42103254d96877c160c09b0b35ff..dd213ecb89734b411deb54a7c3673d3e53003116 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -17,6 +17,13 @@
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/reference_cell.h>
+#include <deal.II/grid/tria.h>
+
 #include <algorithm>
 #include <cmath>
 #include <functional>
@@ -1739,6 +1746,81 @@ QWitherdenVincentSimplex<dim>::QWitherdenVincentSimplex(
 
 
 
+namespace
+{
+  template <int dim>
+  Quadrature<dim>
+  setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
+  {
+    Assert(false, ExcInternalError());
+    return Quadrature<dim>();
+  }
+
+
+
+  Quadrature<1>
+  setup_qiterated_1D(const Quadrature<1> &base_quad,
+                     const unsigned int   n_copies)
+  {
+    return QIterated<1>(base_quad, n_copies);
+  }
+} // namespace
+
+
+
+template <int dim>
+QIteratedSimplex<dim>::QIteratedSimplex(const Quadrature<dim> &base_quad,
+                                        const unsigned int     n_copies)
+{
+  switch (dim)
+    {
+      case 1:
+        static_cast<Quadrature<dim> &>(*this) =
+          setup_qiterated_1D(base_quad, n_copies);
+        break;
+      case 2:
+      case 3:
+        {
+          const auto n_refinements =
+            static_cast<unsigned int>(std::round(std::log2(n_copies)));
+          Assert((1u << n_refinements) == n_copies,
+                 ExcMessage("The number of copies must be a power of 2."));
+          Triangulation<dim> tria;
+          const auto reference_cell = ReferenceCells::get_simplex<dim>();
+          GridGenerator::reference_cell(tria, reference_cell);
+          tria.refine_global(n_refinements);
+          const Mapping<dim> &mapping =
+            reference_cell.template get_default_linear_mapping<dim>();
+          FE_Nothing<dim> fe(reference_cell);
+
+          FEValues<dim>           fe_values(mapping,
+                                  fe,
+                                  base_quad,
+                                  update_quadrature_points | update_JxW_values);
+          std::vector<Point<dim>> points;
+          std::vector<double>     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<Quadrature<dim> &>(*this) =
+            Quadrature<dim>(points, weights);
+
+          break;
+        }
+      default:
+        Assert(false, ExcNotImplemented());
+    }
+}
+
+
+
 template <int dim>
 QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
   : Quadrature<dim>()
@@ -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 (file)
index 0000000..e990469
--- /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 <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include "../tests.h"
+
+// Test QIteratedSimplex accuracy
+
+template <int dim>
+void
+print(const Quadrature<dim> &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 <int dim>
+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<dim> func(monomial_powers);
+  const QIteratedSimplex<dim> quad(QWitherdenVincentSimplex<dim>(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 (file)
index 0000000..98d806d
--- /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