FE_SimplexP: support hanging nodes with cubic elements.

This more-or-less just copies the same code used to place interpolation points
on child lines we use in FE_Q_Base.

diff --git a/doc/news/changes/minor/20241015DavidWells b/doc/news/changes/minor/20241015DavidWells
new file mode 100644 (file)
index 0000000..2eca181
--- /dev/null
+++ b/doc/news/changes/minor/20241015DavidWells
@@ -0,0 +1,3 @@
+Fixed: FESystem now works correctly with cubic FE_SimplexP elements.
+<br>
+(David Wells, 2024/10/15)

diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc
index d3cede7d321003f41f9e07dabbdf7d9c20a52af2..1522a514208c9834b62ff2897fc92038259ea408 100644 (file)
--- a/source/fe/fe_simplex_p.cc
+++ b/source/fe/fe_simplex_p.cc
@@ -218,7 +218,6 @@ namespace
   constraints_fe_p<2>(const unsigned int degree)
   {
     constexpr int dim = 2;
-    Assert(degree <= 3, ExcNotImplemented());
 
     // the following implements the 2d case
     // (the 3d case is not implemented yet)
@@ -229,13 +228,16 @@ namespace
     std::vector<Point<dim - 1>> constraint_points;
     // midpoint
     constraint_points.emplace_back(0.5);
-    if (degree == 2)
-      {
-        // midpoint on subface 0
-        constraint_points.emplace_back(0.25);
-        // midpoint on subface 1
-        constraint_points.emplace_back(0.75);
-      }
+    // subface 0
+    for (unsigned int i = 1; i < degree; ++i)
+      constraint_points.push_back(
+        GeometryInfo<dim - 1>::child_to_cell_coordinates(
+          Point<dim - 1>(i / double(degree)), 0));
+    // subface 1
+    for (unsigned int i = 1; i < degree; ++i)
+      constraint_points.push_back(
+        GeometryInfo<dim - 1>::child_to_cell_coordinates(
+          Point<dim - 1>(i / double(degree)), 1));
 
     // Now construct relation between destination (child) and source (mother)
     // dofs.

diff --git a/tests/simplex/simplex_project_cg.cc b/tests/simplex/simplex_project_cg.cc
index 710e2f610a0db27e705972aa29fbf2090b21fc05..f79fd1d25acdb53ff1ff4700b54cf13490b29a9a 100644 (file)
--- a/tests/simplex/simplex_project_cg.cc
+++ b/tests/simplex/simplex_project_cg.cc
@@ -79,7 +79,7 @@ public:
 
 template <int dim>
 void
-test(const unsigned int degree)
+test(const unsigned int degree, const bool use_amr = false)
 {
 #ifdef SIMPLEX
   FE_SimplexP<dim>   fe(degree);
@@ -90,10 +90,12 @@ test(const unsigned int degree)
 #endif
   deallog << "FE = " << fe.get_name() << std::endl;
   deallog << std::setprecision(4);
+  if (use_amr)
+    deallog << "using AMR" << std::endl;
 
   double previous_error = 1.0;
 
-  for (unsigned int r = 0; r < 4; ++r)
+  for (unsigned int r = 0; r < (use_amr ? 3 : 4); ++r)
     {
       Triangulation<dim> tria_hex, tria;
       GridGenerator::hyper_cube(tria_hex);
@@ -103,6 +105,13 @@ test(const unsigned int degree)
 #else
       tria.copy_triangulation(tria_hex);
 #endif
+      if (use_amr)
+        {
+          for (auto &cell : tria.active_cell_iterators())
+            if (cell->index() % 3 == 0)
+              cell->set_refine_flag();
+          tria.execute_coarsening_and_refinement();
+        }
 
       ReferenceCell   reference_cell = tria.begin_active()->reference_cell();
       DoFHandler<dim> dof_handler(tria);
@@ -111,40 +120,38 @@ test(const unsigned int degree)
       Vector<double>            cell_errors(tria.n_active_cells());
       Vector<double>            solution(dof_handler.n_dofs());
       Solution<dim>             function;
-      AffineConstraints<double> dummy;
+      AffineConstraints<double> constraints;
       const auto               &mapping =
         reference_cell.template get_default_linear_mapping<dim>();
-      dummy.close();
+      DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+      constraints.close();
 
       SparsityPattern sparsity_pattern;
       {
         DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
-        DoFTools::make_sparsity_pattern(dof_handler, dsp);
+        DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, true);
         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);
+      MatrixCreator::create_mass_matrix(
+        mapping, dof_handler, quadrature, h1_matrix, {}, constraints);
+      MatrixCreator::create_laplace_matrix(
+        mapping, dof_handler, quadrature, laplace_matrix, {}, constraints);
 
       h1_matrix.add(0.0, laplace_matrix);
 
       Vector<double> rhs(solution.size());
       VectorTools::create_right_hand_side(
-        mapping, dof_handler, quadrature, ForcingH1<dim>(), rhs);
+        mapping, dof_handler, quadrature, ForcingH1<dim>(), rhs, constraints);
 
       SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false);
       SolverCG<Vector<double>> cg(solver_control);
 
       cg.solve(h1_matrix, solution, rhs, PreconditionIdentity());
+      constraints.distribute(solution);
 
       VectorTools::integrate_difference(mapping,
                                         dof_handler,
@@ -191,4 +198,8 @@ main()
   test<3>(1);
   test<3>(2);
   test<3>(3);
+
+  test<2>(1, true);
+  test<2>(2, true);
+  test<2>(3, true);
 }

diff --git a/tests/simplex/simplex_project_cg.output b/tests/simplex/simplex_project_cg.output
index dfd98e33ac75787637c828dff5bb709f38d9fb7b..97c1f6475babf244265344705810cc9a94dd9057 100644 (file)
--- a/tests/simplex/simplex_project_cg.output
+++ b/tests/simplex/simplex_project_cg.output
@@ -53,3 +53,27 @@ DEAL::L2 error = 2.324e-05
 DEAL::ratio = 15.31
 DEAL::L2 error = 1.447e-06
 DEAL::ratio = 16.06
+DEAL::FE = FE_SimplexP<2>(1)
+DEAL::using AMR
+DEAL::L2 error = 0.07022
+DEAL::ratio = 14.24
+DEAL::L2 error = 0.03106
+DEAL::ratio = 2.261
+DEAL::L2 error = 0.008068
+DEAL::ratio = 3.850
+DEAL::FE = FE_SimplexP<2>(2)
+DEAL::using AMR
+DEAL::L2 error = 0.004448
+DEAL::ratio = 224.8
+DEAL::L2 error = 0.002054
+DEAL::ratio = 2.165
+DEAL::L2 error = 0.0003406
+DEAL::ratio = 6.030
+DEAL::FE = FE_SimplexP<2>(3)
+DEAL::using AMR
+DEAL::L2 error = 0.003640
+DEAL::ratio = 274.7
+DEAL::L2 error = 0.0002153
+DEAL::ratio = 16.91
+DEAL::L2 error = 1.432e-05
+DEAL::ratio = 15.03