From: David Wells <drwells@email.unc.edu>
Date: Sun, 29 Jan 2023 13:47:17 +0000 (-0500)
Subject: FE_Simplex_FE_{DG}P: Add a basic convergence rate test.
X-Git-Tag: relicensing~816^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bac21fd728d791f6b01e42617765676733dc0306;p=dealii.git

FE_Simplex_FE_{DG}P: Add a basic convergence rate test.
---

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 <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_creator.h>
+#include <deal.II/numerics/vector_tools_integrate_difference.h>
+#include <deal.II/numerics/vector_tools_project.h>
+#include <deal.II/numerics/vector_tools_rhs.h>
+
+#include "../tests.h"
+
+#define SIMPLEX
+
+template <int dim>
+class Solution : public Function<dim>
+{
+public:
+  double
+  value(const Point<dim> &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 <int dim>
+class ForcingH1 : public Function<dim>
+{
+public:
+  double
+  value(const Point<dim> &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 <int dim>
+void
+test(const unsigned int degree)
+{
+#ifdef SIMPLEX
+  FE_SimplexP<dim>   fe(degree);
+  QGaussSimplex<dim> quadrature(degree + 2);
+#else
+  FE_Q<dim>   fe(degree);
+  QGauss<dim> 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<dim> 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<dim> dof_handler(tria);
+      dof_handler.distribute_dofs(fe);
+
+      Vector<double>            cell_errors(tria.n_active_cells());
+      Vector<double>            solution(dof_handler.n_dofs());
+      Solution<dim>             function;
+      AffineConstraints<double> dummy;
+      const auto &              mapping =
+        reference_cell.template get_default_linear_mapping<dim>();
+      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<double> h1_matrix(sparsity_pattern);
+          SparseMatrix<double> 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<double> rhs(solution.size());
+          VectorTools::create_right_hand_side(
+            mapping, dof_handler, quadrature, ForcingH1<dim>(), rhs);
+
+          SolverControl            solver_control(1000,
+                                       1e-12 * rhs.l2_norm(),
+                                       false,
+                                       false);
+          SolverCG<Vector<double>> cg(solver_control);
+
+          cg.solve(h1_matrix, solution, rhs, PreconditionIdentity());
+        }
+
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        function,
+                                        cell_errors,
+                                        Quadrature<dim>(
+                                          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<dim> 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 <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools_integrate_difference.h>
+#include <deal.II/numerics/vector_tools_project.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  FE_SimplexDGP<dim> fe(degree);
+  deallog << "FE = " << fe.get_name() << std::endl;
+  QGaussSimplex<dim> quadrature(degree + 1);
+
+  double previous_error = 1.0;
+
+  for (unsigned int r = 0; r < 4; ++r)
+    {
+      Triangulation<dim> 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<dim> dof_handler(tria);
+      dof_handler.distribute_dofs(fe);
+
+      Vector<double>                 cell_errors(tria.n_active_cells());
+      Vector<double>                 solution(dof_handler.n_dofs());
+      Functions::CosineFunction<dim> function;
+      AffineConstraints<double>      dummy;
+      const auto &                   mapping =
+        reference_cell.template get_default_linear_mapping<dim>();
+      dummy.close();
+      VectorTools::project(
+        mapping, dof_handler, dummy, quadrature, function, solution);
+
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        function,
+                                        cell_errors,
+                                        Quadrature<dim>(
+                                          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<dim> 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