From: David Wells <drwells@email.unc.edu>
Date: Fri, 30 Jun 2023 22:21:26 +0000 (-0400)
Subject: Add FE_SimplexP(3) and FE_SimplexDGP(3).
X-Git-Tag: relicensing~395^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d792d234b3d71a2f957a8d374e00fcfd0327f016;p=dealii.git

Add FE_SimplexP(3) and FE_SimplexDGP(3).
---

diff --git a/doc/news/changes/major/20231011DavidWellsPeterMunch b/doc/news/changes/major/20231011DavidWellsPeterMunch
new file mode 100644
index 0000000000..58709a9890
--- /dev/null
+++ b/doc/news/changes/major/20231011DavidWellsPeterMunch
@@ -0,0 +1,3 @@
+New: FE_SimplexP and FE_SimplexDGP now support cubic elements.
+<br>
+(David Wells and Peter Munch, 2023/10/11)
diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc
index 32f42a2167..767d5b1ab5 100644
--- a/source/fe/fe_simplex_p.cc
+++ b/source/fe/fe_simplex_p.cc
@@ -28,6 +28,31 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace
 {
+  // TODO: replace this with QProjector once QProjector learns how to project
+  // quadrature points onto separate faces of simplices
+  template <int dim>
+  Point<dim>
+  face_midpoint(const unsigned int face_no)
+  {
+    const auto reference_cell = ReferenceCells::get_simplex<dim>();
+    const auto face_reference_cell =
+      reference_cell.face_reference_cell(face_no);
+
+    Point<dim> midpoint;
+    for (const auto face_vertex_no : face_reference_cell.vertex_indices())
+      {
+        const auto vertex_no = reference_cell.face_to_cell_vertices(
+          face_no,
+          face_vertex_no,
+          ReferenceCell::default_combined_face_orientation());
+
+        midpoint += reference_cell.template vertex<dim>(vertex_no);
+      }
+
+    midpoint /= reference_cell.face_reference_cell(0).n_vertices();
+    return midpoint;
+  }
+
   /**
    * Helper function to set up the dpo vector of FE_SimplexP for a given @p dim and
    * @p degree.
@@ -44,6 +69,8 @@ namespace
                 return {1, 0};
               case 2:
                 return {1, 1};
+              case 3:
+                return {1, 2};
               default:
                 Assert(false, ExcNotImplemented());
             }
@@ -54,6 +81,8 @@ namespace
                 return {1, 0, 0};
               case 2:
                 return {1, 1, 0};
+              case 3:
+                return {1, 2, 1};
               default:
                 Assert(false, ExcNotImplemented());
             }
@@ -64,6 +93,8 @@ namespace
                 return {1, 0, 0, 0};
               case 2:
                 return {1, 1, 0, 0};
+              case 3:
+                return {1, 2, 1, 0};
               default:
                 Assert(false, ExcNotImplemented());
             }
@@ -83,15 +114,11 @@ namespace
   std::vector<Point<dim>>
   unit_support_points_fe_p(const unsigned int degree)
   {
+    Assert(dim != 0, ExcInternalError());
+    Assert(degree <= 3, ExcNotImplemented());
     std::vector<Point<dim>> unit_points;
-    // If we do dim - 1 we can get here in dim = 0:
-    if (dim == 0)
-      {
-        unit_points.emplace_back();
-        return unit_points;
-      }
+    const auto              reference_cell = ReferenceCells::get_simplex<dim>();
 
-    const auto reference_cell = ReferenceCells::get_simplex<dim>();
     // Piecewise constants are a special case: use a support point at the
     // centroid and only the centroid
     if (degree == 0)
@@ -100,23 +127,45 @@ namespace
         return unit_points;
       }
 
-    Assert(degree <= 2, ExcNotImplemented());
-    for (const auto &vertex_no : reference_cell.vertex_indices())
-      unit_points.emplace_back(reference_cell.template vertex<dim>(vertex_no));
-
-    if (degree == 2)
-      for (const auto &line_no : reference_cell.line_indices())
-        {
-          const auto v0 = reference_cell.template vertex<dim>(
-            reference_cell.line_to_cell_vertices(line_no, 0));
-          const auto v1 = reference_cell.template vertex<dim>(
-            reference_cell.line_to_cell_vertices(line_no, 1));
-          unit_points.emplace_back((v0 + v1) / 2.0);
-        }
+    // otherwise write everything as linear combinations of vertices
+    const auto dpo = get_dpo_vector_fe_p(dim, degree);
+    Assert(dpo.size() == dim + 1, ExcInternalError());
+    Assert(dpo[0] == 1, ExcNotImplemented());
+    // vertices:
+    for (const unsigned int d : reference_cell.vertex_indices())
+      unit_points.push_back(reference_cell.template vertex<dim>(d));
+    // lines:
+    for (const unsigned int l : reference_cell.line_indices())
+      {
+        const Point<dim> p0 =
+          unit_points[reference_cell.line_to_cell_vertices(l, 0)];
+        const Point<dim> p1 =
+          unit_points[reference_cell.line_to_cell_vertices(l, 1)];
+        for (unsigned int p = 0; p < dpo[1]; ++p)
+          unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
+                                (double(p + 1) / (dpo[1] + 1)) * p1);
+      }
+    // quads:
+    if (dim > 1 && dpo[2] > 0)
+      {
+        Assert(dpo[2] == 1, ExcNotImplemented());
+        if (dim == 2)
+          unit_points.push_back(reference_cell.template barycenter<dim>());
+        if (dim == 3)
+          for (const unsigned int f : reference_cell.face_indices())
+            unit_points.push_back(face_midpoint<dim>(f));
+      }
 
     return unit_points;
   }
 
+  template <>
+  std::vector<Point<0>>
+  unit_support_points_fe_p(const unsigned int /*degree*/)
+  {
+    return {Point<0>()};
+  }
+
   /**
    * Set up a vector that contains the unit (reference) cell's faces support
    * points for FE_SimplexP and sufficiently similar elements.
@@ -166,9 +215,8 @@ namespace
   FullMatrix<double>
   constraints_fe_p<2>(const unsigned int degree)
   {
-    const unsigned int dim = 2;
-
-    Assert(degree <= 2, ExcNotImplemented());
+    constexpr int dim = 2;
+    Assert(degree <= 3, ExcNotImplemented());
 
     // the following implements the 2d case
     // (the 3d case is not implemented yet)
@@ -575,7 +623,14 @@ FE_SimplexP<dim, spacedim>::FE_SimplexP(const unsigned int degree)
       unit_support_points_fe_p<dim>(degree),
       unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
       constraints_fe_p<dim>(degree))
-{}
+{
+  if (degree > 2)
+    for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
+      this->adjust_line_dof_index_for_line_orientation_table[i] =
+        this->n_dofs_per_line() - 1 - i - i;
+  // We do not support multiple DoFs per quad yet
+  Assert(degree <= 3, ExcNotImplemented());
+}
 
 
 
diff --git a/tests/simplex/simplex_project_cg.cc b/tests/simplex/simplex_project_cg.cc
index 31132a4752..b0377afbf6 100644
--- a/tests/simplex/simplex_project_cg.cc
+++ b/tests/simplex/simplex_project_cg.cc
@@ -84,12 +84,13 @@ test(const unsigned int degree)
 {
 #ifdef SIMPLEX
   FE_SimplexP<dim>   fe(degree);
-  QGaussSimplex<dim> quadrature(degree + 2);
+  QGaussSimplex<dim> quadrature(degree + 1);
 #else
   FE_Q<dim>   fe(degree);
-  QGauss<dim> quadrature(degree + 2);
+  QGauss<dim> quadrature(degree + 1);
 #endif
   deallog << "FE = " << fe.get_name() << std::endl;
+  deallog << std::setprecision(4);
 
   double previous_error = 1.0;
 
@@ -116,65 +117,52 @@ test(const unsigned int degree)
         reference_cell.template get_default_linear_mapping<dim>();
       dummy.close();
 
-      const bool l2_projection = false;
+      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);
+      }
 
-      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());
-        }
+      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;
+                                        quadrature,
+                                        VectorTools::L2_norm);
+      const double L2_error =
+        VectorTools::compute_global_error(tria,
+                                          cell_errors,
+                                          VectorTools::L2_norm);
+
+      deallog << "L2 error = " << L2_error << std::endl;
+      if (L2_error != 0.0)
+        deallog << "ratio = " << previous_error / L2_error << std::endl;
+      previous_error = L2_error;
 
 #if 0
       if (dim == 2)
@@ -199,7 +187,9 @@ main()
 
   test<2>(1);
   test<2>(2);
+  test<2>(3);
 
   test<3>(1);
   test<3>(2);
+  test<3>(3);
 }
diff --git a/tests/simplex/simplex_project_cg.output b/tests/simplex/simplex_project_cg.output
index 287c3090a0..dfd98e33ac 100644
--- a/tests/simplex/simplex_project_cg.output
+++ b/tests/simplex/simplex_project_cg.output
@@ -1,37 +1,55 @@
 
 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::L2 error = 0.05224
+DEAL::ratio = 19.14
+DEAL::L2 error = 0.02949
+DEAL::ratio = 1.772
+DEAL::L2 error = 0.007768
+DEAL::ratio = 3.796
+DEAL::L2 error = 0.001951
+DEAL::ratio = 3.981
 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::L2 error = 0.004078
+DEAL::ratio = 245.2
+DEAL::L2 error = 0.002064
+DEAL::ratio = 1.975
+DEAL::L2 error = 0.0003512
+DEAL::ratio = 5.879
+DEAL::L2 error = 5.040e-05
+DEAL::ratio = 6.967
+DEAL::FE = FE_SimplexP<2>(3)
+DEAL::L2 error = 0.004061
+DEAL::ratio = 246.3
+DEAL::L2 error = 0.0002464
+DEAL::ratio = 16.48
+DEAL::L2 error = 1.603e-05
+DEAL::ratio = 15.38
+DEAL::L2 error = 1.009e-06
+DEAL::ratio = 15.88
 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::L2 error = 0.09826
+DEAL::ratio = 10.18
+DEAL::L2 error = 0.02938
+DEAL::ratio = 3.344
+DEAL::L2 error = 0.007387
+DEAL::ratio = 3.977
+DEAL::L2 error = 0.001763
+DEAL::ratio = 4.191
 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
+DEAL::L2 error = 0.01087
+DEAL::ratio = 91.95
+DEAL::L2 error = 0.004092
+DEAL::ratio = 2.657
+DEAL::L2 error = 0.0007005
+DEAL::ratio = 5.842
+DEAL::L2 error = 0.0001018
+DEAL::ratio = 6.880
+DEAL::FE = FE_SimplexP<3>(3)
+DEAL::L2 error = 0.006416
+DEAL::ratio = 155.9
+DEAL::L2 error = 0.0003557
+DEAL::ratio = 18.03
+DEAL::L2 error = 2.324e-05
+DEAL::ratio = 15.31
+DEAL::L2 error = 1.447e-06
+DEAL::ratio = 16.06
diff --git a/tests/simplex/simplex_project_dg.cc b/tests/simplex/simplex_project_dg.cc
index 5edb6fae80..4d277490c6 100644
--- a/tests/simplex/simplex_project_dg.cc
+++ b/tests/simplex/simplex_project_dg.cc
@@ -20,18 +20,29 @@
 #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"
 
@@ -41,6 +52,7 @@ test(const unsigned int degree)
 {
   FE_SimplexDGP<dim> fe(degree);
   deallog << "FE = " << fe.get_name() << std::endl;
+  deallog << std::setprecision(4);
   QGaussSimplex<dim> quadrature(degree + 1);
 
   double previous_error = 1.0;
@@ -63,24 +75,44 @@ test(const unsigned int degree)
       const auto                    &mapping =
         reference_cell.template get_default_linear_mapping<dim>();
       dummy.close();
-      VectorTools::project(
-        mapping, dof_handler, dummy, quadrature, function, solution);
+
+      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> mass_matrix(sparsity_pattern);
+      MatrixCreator::create_mass_matrix(mapping,
+                                        dof_handler,
+                                        quadrature,
+                                        mass_matrix);
+      Vector<double> rhs(solution.size());
+      VectorTools::create_right_hand_side(
+        mapping, dof_handler, quadrature, function, rhs);
+
+      SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false);
+      SolverCG<Vector<double>> cg(solver_control);
+
+      cg.solve(mass_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;
+                                        quadrature,
+                                        VectorTools::L2_norm);
+      const double L2_error =
+        VectorTools::compute_global_error(tria,
+                                          cell_errors,
+                                          VectorTools::L2_norm);
+
+      deallog << "L2 error = " << L2_error << std::endl;
+      if (L2_error != 0.0)
+        deallog << "ratio = " << previous_error / L2_error << std::endl;
+      previous_error = L2_error;
 
 #if 0
       if (dim == 2)
@@ -105,7 +137,9 @@ main()
 
   test<2>(1);
   test<2>(2);
+  test<2>(3);
 
   test<3>(1);
   test<3>(2);
+  test<3>(3);
 }
diff --git a/tests/simplex/simplex_project_dg.output b/tests/simplex/simplex_project_dg.output
index 4bb911968d..1b38184f1a 100644
--- a/tests/simplex/simplex_project_dg.output
+++ b/tests/simplex/simplex_project_dg.output
@@ -1,37 +1,55 @@
 
 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::L2 error = 0.009794
+DEAL::ratio = 102.1
+DEAL::L2 error = 0.002479
+DEAL::ratio = 3.952
+DEAL::L2 error = 0.0006215
+DEAL::ratio = 3.988
+DEAL::L2 error = 0.0001555
+DEAL::ratio = 3.997
 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::L2 error = 0.0002017
+DEAL::ratio = 4957.
+DEAL::L2 error = 2.334e-05
+DEAL::ratio = 8.644
+DEAL::L2 error = 2.854e-06
+DEAL::ratio = 8.178
+DEAL::L2 error = 3.547e-07
+DEAL::ratio = 8.046
+DEAL::FE = FE_SimplexDGP<2>(3)
+DEAL::L2 error = 0.0001932
+DEAL::ratio = 5175.
+DEAL::L2 error = 1.224e-05
+DEAL::ratio = 15.79
+DEAL::L2 error = 7.677e-07
+DEAL::ratio = 15.95
+DEAL::L2 error = 4.802e-08
+DEAL::ratio = 15.99
 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::L2 error = 0.01241
+DEAL::ratio = 80.59
+DEAL::L2 error = 0.003150
+DEAL::ratio = 3.939
+DEAL::L2 error = 0.0007879
+DEAL::ratio = 3.998
+DEAL::L2 error = 0.0001968
+DEAL::ratio = 4.004
 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
+DEAL::L2 error = 0.001633
+DEAL::ratio = 612.3
+DEAL::L2 error = 0.0002072
+DEAL::ratio = 7.883
+DEAL::L2 error = 2.599e-05
+DEAL::ratio = 7.971
+DEAL::L2 error = 3.252e-06
+DEAL::ratio = 7.993
+DEAL::FE = FE_SimplexDGP<3>(3)
+DEAL::L2 error = 0.0002064
+DEAL::ratio = 4845.
+DEAL::L2 error = 1.309e-05
+DEAL::ratio = 15.77
+DEAL::L2 error = 8.214e-07
+DEAL::ratio = 15.94
+DEAL::L2 error = 5.138e-08
+DEAL::ratio = 15.99