Tests for hanging nodes

diff --git a/tests/matrix_free/matrix_vector_rt_03.cc b/tests/matrix_free/matrix_vector_rt_03.cc
new file mode 100644 (file)
index 0000000..d817e02
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_rt_03.cc
@@ -0,0 +1,74 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// This test is the same as matrix_vector_rt_01.cc but with hanging nodes.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(2);
+  else
+    tria.refine_global(1);
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
+  cell = tria.begin_active(), endc = tria.end();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 0.2)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 10 - 3 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+  DoFHandler<dim>            dof(tria);
+  dof.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  constraints.close();
+
+  deallog << "Using " << dof.get_fe().get_name() << std::endl;
+  deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+          << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+          << std::endl;
+  do_test<dim, fe_degree, double>(dof, constraints, TestType::values_gradients);
+}

diff --git a/tests/matrix_free/matrix_vector_rt_03.output b/tests/matrix_free/matrix_vector_rt_03.output
new file mode 100644 (file)
index 0000000..f327e82
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_rt_03.output
@@ -0,0 +1,29 @@
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 481
+DEAL:2d::Number of degrees of freedom: 4548
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients 
+DEAL:2d::Norm of difference: 1.13503e-13
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 481
+DEAL:2d::Number of degrees of freedom: 9708
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients 
+DEAL:2d::Norm of difference: 2.05175e-13
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 85
+DEAL:3d::Number of degrees of freedom: 2392
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients 
+DEAL:3d::Norm of difference: 1.98105e-14
+DEAL:3d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Number of cells: 85
+DEAL:3d::Number of degrees of freedom: 7677
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients 
+DEAL:3d::Norm of difference: 3.46605e-14
+DEAL:3d::

diff --git a/tests/matrix_free/matrix_vector_rt_common.h b/tests/matrix_free/matrix_vector_rt_common.h
index 90ea0b41c8bc24a19408a7ae85fb5fab404e7f69..40d7022da12aae7cfe8d37562057b9377f0c2742 100644 (file)
--- a/tests/matrix_free/matrix_vector_rt_common.h
+++ b/tests/matrix_free/matrix_vector_rt_common.h
@@ -173,10 +173,11 @@ do_test(const DoFHandler<dim> &          dof,
   MatrixFree<dim, Number> mf_data;
   {
     const QGaussLobatto<1>                           quad(fe_degree + 2);
+    const MappingQ<dim>                              mapping(fe_degree + 2);
     typename MatrixFree<dim, Number>::AdditionalData data;
     data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
     data.mapping_update_flags  = update_gradients | update_hessians;
-    mf_data.reinit(dof, constraints, quad, data);
+    mf_data.reinit(mapping, dof, constraints, quad, data);
   }
 
   // create vector with random entries
@@ -191,23 +192,27 @@ do_test(const DoFHandler<dim> &          dof,
         continue;
       initial_condition[i] = random_value<Number>();
     }
+  constraints.distribute(initial_condition);
 
   // MatrixFree solution
   MatrixFreeTest<dim, -1, 0, Number> mf(mf_data, test_type);
   mf.test_functions(solution, initial_condition);
 
 
+  // Evaluation with FEValues
   SparsityPattern        sp;
   SparseMatrix<double>   system_matrix;
   DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
-  DoFTools::make_sparsity_pattern(dof, dsp);
+  DoFTools::make_sparsity_pattern(dof, dsp, constraints, false);
   sp.copy_from(dsp);
   system_matrix.reinit(sp);
 
-  FEValues<dim> fe_val(dof.get_fe(),
+  const MappingQ<dim> mapping(fe_degree + 2);
+  FEValues<dim>       fe_val(mapping,
+                       dof.get_fe(),
                        Quadrature<dim>(mf_data.get_quadrature(0)),
-                       update_values | update_gradients | update_JxW_values);
-
+                       update_values | update_gradients | update_piola |
+                         update_JxW_values);
 
   const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
@@ -248,17 +253,16 @@ do_test(const DoFHandler<dim> &          dof,
               }
         }
       cell->get_dof_indices(local_dof_indices);
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int j = 0; j < dofs_per_cell; ++j)
-          system_matrix.add(local_dof_indices[i],
-                            local_dof_indices[j],
-                            local_matrix(i, j));
+      constraints.distribute_local_to_global(local_matrix,
+                                             local_dof_indices,
+                                             system_matrix);
     }
 
   Vector<Number> ref(solution.size());
 
   // Compute reference
   system_matrix.vmult(ref, initial_condition);
+  constraints.set_zero(ref);
 
   ref -= solution;
 

diff --git a/tests/matrix_free/matrix_vector_rt_face_03.cc b/tests/matrix_free/matrix_vector_rt_face_03.cc
new file mode 100644 (file)
index 0000000..66736ef
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_rt_face_03.cc
@@ -0,0 +1,70 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// This test does the same as matrix_vector_rt_face_01.cc but also uses hanging
+// nodes.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_face_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
+  cell = tria.begin_active(), endc = tria.end();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 0.2)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 10 - 3 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+  DoFHandler<dim>            dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  constraints.close();
+
+  deallog << "Using " << dof.get_fe().get_name() << std::endl;
+  deallog << "Number of cells: " << tria.n_active_cells() << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+          << std::endl;
+  do_test<dim, fe_degree, double>(dof, constraints, TestType::values_gradients);
+  //   do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
+}

diff --git a/tests/matrix_free/matrix_vector_rt_face_03.output b/tests/matrix_free/matrix_vector_rt_face_03.output
new file mode 100644 (file)
index 0000000..31d0141
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_rt_face_03.output
@@ -0,0 +1,22 @@
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 85
+DEAL:2d::Number of degrees of freedom: 798
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients 
+DEAL:2d::Norm of difference: 1.46726e-13
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 85
+DEAL:2d::Number of degrees of freedom: 1707
+DEAL:2d::
+DEAL:2d::Testing Values and Gradients 
+DEAL:2d::Norm of difference: 1.47883e-13
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 22
+DEAL:3d::Number of degrees of freedom: 672
+DEAL:3d::
+DEAL:3d::Testing Values and Gradients 
+DEAL:3d::Norm of difference: 3.31436e-14
+DEAL:3d::

diff --git a/tests/matrix_free/matrix_vector_rt_face_common.h b/tests/matrix_free/matrix_vector_rt_face_common.h
index dd65b7f3917bfb5bd80a0c3891f82f2dd0c9483c..39a3ffa46d4b7a8d710f0605634ab46cc0c4854b 100644 (file)
--- a/tests/matrix_free/matrix_vector_rt_face_common.h
+++ b/tests/matrix_free/matrix_vector_rt_face_common.h
@@ -225,17 +225,21 @@ do_test(const DoFHandler<dim> &          dof,
 {
   deallog << "Testing " << enum_to_string(test_type) << std::endl;
 
+  constexpr unsigned int n_q_points = fe_degree + 2;
+  constexpr unsigned int m_degree   = fe_degree + 2;
+
+  const MappingQ<dim>     mapping(m_degree);
   MatrixFree<dim, Number> mf_data;
   {
-    const QGaussLobatto<1>                           quad(fe_degree + 2);
-    const MappingQ<dim>                              mapping(fe_degree);
+    const QGaussLobatto<1>                           quad(n_q_points);
     typename MatrixFree<dim, Number>::AdditionalData data;
     data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
-    data.mapping_update_flags  = update_gradients | update_JxW_values;
+    data.mapping_update_flags =
+      (update_gradients | update_JxW_values | update_hessians);
     data.mapping_update_flags_inner_faces =
-      (update_gradients | update_JxW_values);
+      (update_gradients | update_JxW_values | update_hessians);
     data.mapping_update_flags_boundary_faces =
-      (update_gradients | update_JxW_values);
+      (update_gradients | update_JxW_values | update_hessians);
     mf_data.reinit(mapping, dof, constraints, quad, data);
   }
 
@@ -251,23 +255,25 @@ do_test(const DoFHandler<dim> &          dof,
         continue;
       initial_condition[i] = random_value<Number>();
     }
+  constraints.distribute(initial_condition);
 
   MatrixFreeTest<dim, -1, 0, Number> mf(mf_data, test_type);
   mf.test_functions(solution, initial_condition);
 
 
+  // Evaluation with FEFaceValues
   SparsityPattern        sp;
   SparseMatrix<double>   system_matrix;
   DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
-  DoFTools::make_sparsity_pattern(dof, dsp);
+  DoFTools::make_sparsity_pattern(dof, dsp, constraints, false);
   sp.copy_from(dsp);
   system_matrix.reinit(sp);
 
-  FEFaceValues<dim> fe_val(
-    dof.get_fe(),
-    QGaussLobatto<dim - 1>(mf_data.get_quadrature(0).size()),
-    update_values | update_gradients | update_JxW_values);
-
+  FEFaceValues<dim> fe_val(mapping,
+                           dof.get_fe(),
+                           QGaussLobatto<dim - 1>(n_q_points),
+                           update_values | update_gradients |
+                             update_JxW_values | update_piola);
 
   const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
@@ -308,17 +314,16 @@ do_test(const DoFHandler<dim> &          dof,
                 }
           }
         cell->get_dof_indices(local_dof_indices);
-        for (unsigned int i = 0; i < dofs_per_cell; ++i)
-          for (unsigned int j = 0; j < dofs_per_cell; ++j)
-            system_matrix.add(local_dof_indices[i],
-                              local_dof_indices[j],
-                              local_matrix(i, j));
+        constraints.distribute_local_to_global(local_matrix,
+                                               local_dof_indices,
+                                               system_matrix);
       }
 
   Vector<Number> ref(solution.size());
 
   // Compute reference
   system_matrix.vmult(ref, initial_condition);
+  constraints.set_zero(ref);
 
   ref -= solution;