From bac21fd728d791f6b01e42617765676733dc0306 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 29 Jan 2023 08:47:17 -0500 Subject: [PATCH] FE_Simplex_FE_{DG}P: Add a basic convergence rate test. --- tests/simplex/simplex_project_cg.cc | 205 ++++++++++++++++++++++++ tests/simplex/simplex_project_cg.output | 37 +++++ tests/simplex/simplex_project_dg.cc | 111 +++++++++++++ tests/simplex/simplex_project_dg.output | 37 +++++ 4 files changed, 390 insertions(+) create mode 100644 tests/simplex/simplex_project_cg.cc create mode 100644 tests/simplex/simplex_project_cg.output create mode 100644 tests/simplex/simplex_project_dg.cc create mode 100644 tests/simplex/simplex_project_dg.output diff --git a/tests/simplex/simplex_project_cg.cc b/tests/simplex/simplex_project_cg.cc new file mode 100644 index 0000000000..e7d3ac6e93 --- /dev/null +++ b/tests/simplex/simplex_project_cg.cc @@ -0,0 +1,205 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2022 - 2022 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. + * + * --------------------------------------------------------------------- + */ + +// Verify convergence rates for continuous simplex elements. + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "../tests.h" + +#define SIMPLEX + +template +class Solution : public Function +{ +public: + double + value(const Point &p, const unsigned int = 0) const override + { + double u = 1.0; + for (int d = 0; d < dim; ++d) + u *= std::sin(numbers::PI * p(d)); + + return u; + } +}; + +template +class ForcingH1 : public Function +{ +public: + double + value(const Point &p, const unsigned int = 0) const override + { + double u = 1.0; + for (int d = 0; d < dim; ++d) + u *= std::sin(numbers::PI * p(d)); + + return (1 + 0.0 * dim * numbers::PI * numbers::PI) * u; + } +}; + +template +void +test(const unsigned int degree) +{ +#ifdef SIMPLEX + FE_SimplexP fe(degree); + QGaussSimplex quadrature(degree + 2); +#else + FE_Q fe(degree); + QGauss quadrature(degree + 2); +#endif + deallog << "FE = " << fe.get_name() << std::endl; + + double previous_error = 1.0; + + for (unsigned int r = 0; r < 4; ++r) + { + Triangulation tria_hex, tria; + GridGenerator::hyper_cube(tria_hex); + tria_hex.refine_global(r); +#ifdef SIMPLEX + GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria); +#else + tria.copy_triangulation(tria_hex); +#endif + + ReferenceCell reference_cell = tria.begin_active()->reference_cell(); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector cell_errors(tria.n_active_cells()); + Vector solution(dof_handler.n_dofs()); + Solution function; + AffineConstraints dummy; + const auto & mapping = + reference_cell.template get_default_linear_mapping(); + dummy.close(); + + const bool l2_projection = false; + + if (l2_projection) + { + VectorTools::project( + mapping, dof_handler, dummy, quadrature, function, solution); + } + else + { + SparsityPattern sparsity_pattern; + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + } + + SparseMatrix h1_matrix(sparsity_pattern); + SparseMatrix laplace_matrix(sparsity_pattern); + + MatrixCreator::create_mass_matrix(mapping, + dof_handler, + quadrature, + h1_matrix); + MatrixCreator::create_laplace_matrix(mapping, + dof_handler, + quadrature, + laplace_matrix); + + h1_matrix.add(0.0, laplace_matrix); + + Vector rhs(solution.size()); + VectorTools::create_right_hand_side( + mapping, dof_handler, quadrature, ForcingH1(), rhs); + + SolverControl solver_control(1000, + 1e-12 * rhs.l2_norm(), + false, + false); + SolverCG> cg(solver_control); + + cg.solve(h1_matrix, solution, rhs, PreconditionIdentity()); + } + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + function, + cell_errors, + Quadrature( + fe.get_unit_support_points()), + VectorTools::Linfty_norm); + + const double max_error = + *std::max_element(cell_errors.begin(), cell_errors.end()); + deallog << "max error = " << max_error << std::endl; + if (max_error != 0.0) + deallog << "ratio = " << previous_error / max_error << std::endl; + previous_error = max_error; + +#if 0 + if (dim == 2) + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u"); + data_out.build_patches(2); + + std::ofstream output("out-" + std::to_string(degree) + "-" + + std::to_string(r) + ".vtu"); + data_out.write_vtu(output); + } +#endif + } +} + +int +main() +{ + initlog(); + + test<2>(1); + test<2>(2); + + test<3>(1); + test<3>(2); +} diff --git a/tests/simplex/simplex_project_cg.output b/tests/simplex/simplex_project_cg.output new file mode 100644 index 0000000000..287c3090a0 --- /dev/null +++ b/tests/simplex/simplex_project_cg.output @@ -0,0 +1,37 @@ + +DEAL::FE = FE_SimplexP<2>(1) +DEAL::max error = 0.332362 +DEAL::ratio = 3.00877 +DEAL::max error = 0.100426 +DEAL::ratio = 3.30953 +DEAL::max error = 0.0271946 +DEAL::ratio = 3.69285 +DEAL::max error = 0.00747939 +DEAL::ratio = 3.63594 +DEAL::FE = FE_SimplexP<2>(2) +DEAL::max error = 0.0430400 +DEAL::ratio = 23.2342 +DEAL::max error = 0.0116055 +DEAL::ratio = 3.70859 +DEAL::max error = 0.00232162 +DEAL::ratio = 4.99887 +DEAL::max error = 0.000336599 +DEAL::ratio = 6.89729 +DEAL::FE = FE_SimplexP<3>(1) +DEAL::max error = 0.620494 +DEAL::ratio = 1.61162 +DEAL::max error = 0.259138 +DEAL::ratio = 2.39445 +DEAL::max error = 0.0619310 +DEAL::ratio = 4.18431 +DEAL::max error = 0.0151022 +DEAL::ratio = 4.10080 +DEAL::FE = FE_SimplexP<3>(2) +DEAL::max error = 0.0550999 +DEAL::ratio = 18.1489 +DEAL::max error = 0.0466852 +DEAL::ratio = 1.18024 +DEAL::max error = 0.00502563 +DEAL::ratio = 9.28942 +DEAL::max error = 0.000810078 +DEAL::ratio = 6.20388 diff --git a/tests/simplex/simplex_project_dg.cc b/tests/simplex/simplex_project_dg.cc new file mode 100644 index 0000000000..1ec5d28a6f --- /dev/null +++ b/tests/simplex/simplex_project_dg.cc @@ -0,0 +1,111 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2022 - 2022 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. + * + * --------------------------------------------------------------------- + */ + +// Verify convergence rates for discontinuous simplex elements. + +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +template +void +test(const unsigned int degree) +{ + FE_SimplexDGP fe(degree); + deallog << "FE = " << fe.get_name() << std::endl; + QGaussSimplex quadrature(degree + 1); + + double previous_error = 1.0; + + for (unsigned int r = 0; r < 4; ++r) + { + Triangulation tria_hex, tria; + GridGenerator::hyper_cube(tria_hex); + tria_hex.refine_global(r); + GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria); + + ReferenceCell reference_cell = tria.begin_active()->reference_cell(); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector cell_errors(tria.n_active_cells()); + Vector solution(dof_handler.n_dofs()); + Functions::CosineFunction function; + AffineConstraints dummy; + const auto & mapping = + reference_cell.template get_default_linear_mapping(); + dummy.close(); + VectorTools::project( + mapping, dof_handler, dummy, quadrature, function, solution); + + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + function, + cell_errors, + Quadrature( + fe.get_unit_support_points()), + VectorTools::Linfty_norm); + + const double max_error = + *std::max_element(cell_errors.begin(), cell_errors.end()); + deallog << "max error = " << max_error << std::endl; + if (max_error != 0.0) + deallog << "ratio = " << previous_error / max_error << std::endl; + previous_error = max_error; + +#if 0 + if (dim == 2) + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u"); + data_out.build_patches(2); + + std::ofstream output("out-" + std::to_string(degree) + "-" + + std::to_string(r) + ".vtu"); + data_out.write_vtu(output); + } +#endif + } +} + +int +main() +{ + initlog(); + + test<2>(1); + test<2>(2); + + test<3>(1); + test<3>(2); +} diff --git a/tests/simplex/simplex_project_dg.output b/tests/simplex/simplex_project_dg.output new file mode 100644 index 0000000000..4bb911968d --- /dev/null +++ b/tests/simplex/simplex_project_dg.output @@ -0,0 +1,37 @@ + +DEAL::FE = FE_SimplexDGP<2>(1) +DEAL::max error = 0.113889 +DEAL::ratio = 8.78048 +DEAL::max error = 0.0302392 +DEAL::ratio = 3.76627 +DEAL::max error = 0.00767275 +DEAL::ratio = 3.94112 +DEAL::max error = 0.00192529 +DEAL::ratio = 3.98525 +DEAL::FE = FE_SimplexDGP<2>(2) +DEAL::max error = 0.0174164 +DEAL::ratio = 57.4171 +DEAL::max error = 0.00227404 +DEAL::ratio = 7.65880 +DEAL::max error = 0.000287343 +DEAL::ratio = 7.91404 +DEAL::max error = 3.60148e-05 +DEAL::ratio = 7.97847 +DEAL::FE = FE_SimplexDGP<3>(1) +DEAL::max error = 0.208217 +DEAL::ratio = 4.80269 +DEAL::max error = 0.0629409 +DEAL::ratio = 3.30813 +DEAL::max error = 0.0164605 +DEAL::ratio = 3.82376 +DEAL::max error = 0.00416117 +DEAL::ratio = 3.95573 +DEAL::FE = FE_SimplexDGP<3>(2) +DEAL::max error = 0.0415553 +DEAL::ratio = 24.0643 +DEAL::max error = 0.00641441 +DEAL::ratio = 6.47844 +DEAL::max error = 0.000842371 +DEAL::ratio = 7.61471 +DEAL::max error = 0.000106584 +DEAL::ratio = 7.90336 -- 2.39.5