Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Implement QIteratedSimplex.
A limited version of QIterated for simplex meshes.
  • Loading branch information
drwells committed Jul 16, 2021
1 parent a8d2d28 commit 277affa
Show file tree
Hide file tree
Showing 5 changed files with 351 additions and 0 deletions.
3 changes: 3 additions & 0 deletions 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)
15 changes: 15 additions & 0 deletions include/deal.II/base/quadrature_lib.h
Expand Up @@ -858,6 +858,21 @@ class QWitherdenVincentSimplex : public QSimplex<dim>
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.
*/
Expand Down
86 changes: 86 additions & 0 deletions source/base/quadrature_lib.cc
Expand Up @@ -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>
Expand Down Expand Up @@ -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>()
Expand Down Expand Up @@ -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>;
Expand Down
99 changes: 99 additions & 0 deletions 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);
}
148 changes: 148 additions & 0 deletions 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

0 comments on commit 277affa

Please sign in to comment.