From f00eb8c9c042fff45f146ea348e6c1f5297687eb Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 22 Jun 2021 14:00:07 -0400 Subject: [PATCH] Add a test with poor output. This test passes on master, but shouldn't, since several of the printed integrands are wrong. This will be corrected in the next commit. --- tests/simplex/qgauss_01.cc | 82 ++++++++++++++++++++++++++++++++++ tests/simplex/qgauss_01.output | 19 ++++++++ 2 files changed, 101 insertions(+) create mode 100644 tests/simplex/qgauss_01.cc create mode 100644 tests/simplex/qgauss_01.output diff --git a/tests/simplex/qgauss_01.cc b/tests/simplex/qgauss_01.cc new file mode 100644 index 0000000000..a517fbb910 --- /dev/null +++ b/tests/simplex/qgauss_01.cc @@ -0,0 +1,82 @@ +// --------------------------------------------------------------------- +// +// 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 QGaussSimplex 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 QGaussSimplex 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 << 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<3>(1); + check_accuracy_1D<3>(2); + check_accuracy_1D<3>(3); + check_accuracy_1D<3>(4); +} diff --git a/tests/simplex/qgauss_01.output b/tests/simplex/qgauss_01.output new file mode 100644 index 0000000000..d69117871d --- /dev/null +++ b/tests/simplex/qgauss_01.output @@ -0,0 +1,19 @@ + +DEAL:: +DEAL:: +DEAL::Monomial powers = 0.00000 1.00000 +DEAL::Integrand = 0.1666666666666667 +DEAL::Monomial powers = 1.00000 2.00000 +DEAL::Integrand = 0.01620370370370370 +DEAL::Monomial powers = 2.00000 3.00000 +DEAL::Integrand = 0.002380952380952380 +DEAL::Monomial powers = 3.00000 4.00000 +DEAL::Integrand = 0.0003968253968253969 +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.001507514161979123 +DEAL::Monomial powers = 1.00000 1.00000 3.00000 +DEAL::Integrand = 0.0001578591272523996 +DEAL::Monomial powers = 2.00000 2.00000 3.00000 +DEAL::Integrand = 6.613756613756609e-06 -- 2.39.5