Add face_to_cell_index_nodal for simplex elements

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 9858856a05a1755bcfbaea6111f8eb9036b9b4eb..8d12c435a7870d9d6d2d8e03d0dde9a31dc74ff9 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -496,6 +496,84 @@ namespace internal
 
           univariate_shape_data.nodal_at_cell_boundaries = true;
 
+          const ReferenceCell reference_cell = fe.reference_cell();
+          if (reference_cell.is_simplex())
+            {
+              if (dim == 2)
+                dofs_per_component_on_face = fe.degree + 1;
+              else
+                dofs_per_component_on_face =
+                  (fe.degree + 1) * (fe.degree + 2) / 2;
+
+              face_to_cell_index_nodal.reinit(reference_cell.n_faces(),
+                                              dofs_per_component_on_face);
+
+
+              for (unsigned int face = 0; face < reference_cell.n_faces();
+                   ++face)
+                {
+                  // first get info from reference cell, i.e. the linear case
+                  unsigned int d = 0;
+                  for (; d < dim; ++d)
+                    face_to_cell_index_nodal[face][d] =
+                      reference_cell.face_to_cell_vertices(face, d, 1);
+
+                  // now fill the rest of the indices, start with the lines
+                  if (fe.degree == 2)
+                    for (; d < dofs_per_component_on_face; ++d)
+                      face_to_cell_index_nodal[face][d] =
+                        reference_cell.n_vertices() +
+                        reference_cell.face_to_cell_lines(face, d - dim, 1);
+
+                  // in the cubic case it is more complicated as more DoFs are
+                  // on the lines
+                  else if (fe.degree == 3)
+                    {
+                      for (unsigned int line = 0;
+                           d < dofs_per_component_on_face - 1;
+                           ++line, d += 2)
+                        {
+                          const unsigned int face_to_cell_lines =
+                            reference_cell.face_to_cell_lines(face, line, 1);
+                          // check the direction of the line
+                          // is it 0 -> 1 or 1 -> 0
+                          // as DoFs on the line are ordered differently
+                          if (reference_cell.line_to_cell_vertices(
+                                face_to_cell_lines, 0) ==
+                              reference_cell.face_to_cell_vertices(face,
+                                                                   line,
+                                                                   1))
+                            {
+                              face_to_cell_index_nodal[face][d] =
+                                reference_cell.n_vertices() +
+                                2 * face_to_cell_lines;
+                              face_to_cell_index_nodal[face][d + 1] =
+                                reference_cell.n_vertices() +
+                                2 * face_to_cell_lines + 1;
+                            }
+                          else
+                            {
+                              face_to_cell_index_nodal[face][d + 1] =
+                                reference_cell.n_vertices() +
+                                2 * face_to_cell_lines;
+                              face_to_cell_index_nodal[face][d] =
+                                reference_cell.n_vertices() +
+                                2 * face_to_cell_lines + 1;
+                            }
+                        }
+                      //  in 3D we also need the DoFs on the quads
+                      if (dim == 3)
+                        {
+                          face_to_cell_index_nodal
+                            [face][dofs_per_component_on_face - 1] =
+                              reference_cell.n_vertices() +
+                              2 * reference_cell.n_lines() + face;
+                        }
+                    }
+                  else if (fe.degree > 3)
+                    DEAL_II_NOT_IMPLEMENTED();
+                }
+            }
           // TODO: set up face_to_cell_index_nodal, face_to_cell_index_hermite,
           //  face_orientations
 

diff --git a/tests/matrix_free/face_to_cell_index_nodal_simplex.cc b/tests/matrix_free/face_to_cell_index_nodal_simplex.cc
new file mode 100644 (file)
index 0000000..1a8805d
--- /dev/null
+++ b/tests/matrix_free/face_to_cell_index_nodal_simplex.cc
@@ -0,0 +1,102 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2020 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// test the correctness of face_to_cell_index_nodal
+// in internal::MatrixFreeFunctions::ShapeInfo
+// for simplex elements
+
+#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/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/shape_info.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(const FiniteElement<dim> &fe, const Quadrature<dim> &quad)
+{
+  for (unsigned int i = 0; i < fe.n_base_elements(); ++i)
+    {
+      internal::MatrixFreeFunctions::ShapeInfo<double> shape_info;
+      shape_info.reinit(quad, fe, i);
+      deallog << "Testing: " << fe.get_name() << std::endl;
+
+      Triangulation<dim> tria;
+      GridGenerator::reference_cell(tria, fe.reference_cell());
+
+      FE_SimplexP<dim> fe_continous(fe.degree);
+
+      DoFHandler<dim> dof_handler(tria);
+      dof_handler.distribute_dofs(fe_continous);
+
+      std::vector<types::global_dof_index> dof_indices(
+        fe_continous.n_dofs_per_face());
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          unsigned int face_counter = 0;
+          for (const auto &face : cell->face_iterators())
+            {
+              face->get_dof_indices(dof_indices);
+
+              unsigned int dof_counter = 0;
+              for (const auto i : dof_indices)
+                {
+                  if (i == shape_info.face_to_cell_index_nodal[face_counter]
+                                                              [dof_counter])
+                    deallog << "Ok ";
+                  else
+                    deallog << "Not ok";
+                  ++dof_counter;
+                }
+              deallog << std::endl;
+              ++face_counter;
+            }
+        }
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+  for (unsigned int degree = 1; degree < 4; ++degree)
+    {
+      FE_SimplexP<2>   fe(degree);
+      FE_SimplexDGP<2> fe_dg(degree);
+      QGaussSimplex<2> quad(1);
+      test(fe, quad);
+      test(fe_dg, quad);
+
+
+      FE_SimplexP<3>   fe_3(degree);
+      FE_SimplexDGP<3> fe_dg_3(degree);
+      QGaussSimplex<3> quad_3(1);
+      test(fe_3, quad_3);
+      test(fe_dg_3, quad_3);
+    }
+}

diff --git a/tests/matrix_free/face_to_cell_index_nodal_simplex.output b/tests/matrix_free/face_to_cell_index_nodal_simplex.output
new file mode 100644 (file)
index 0000000..422b0e5
--- /dev/null
+++ b/tests/matrix_free/face_to_cell_index_nodal_simplex.output
@@ -0,0 +1,55 @@
+
+DEAL::Testing: FE_SimplexP<2>(1)
+DEAL::Ok Ok 
+DEAL::Ok Ok 
+DEAL::Ok Ok 
+DEAL::Testing: FE_SimplexDGP<2>(1)
+DEAL::Ok Ok 
+DEAL::Ok Ok 
+DEAL::Ok Ok 
+DEAL::Testing: FE_SimplexP<3>(1)
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Testing: FE_SimplexDGP<3>(1)
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Testing: FE_SimplexP<2>(2)
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Testing: FE_SimplexDGP<2>(2)
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Ok Ok Ok 
+DEAL::Testing: FE_SimplexP<3>(2)
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Testing: FE_SimplexDGP<3>(2)
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok 
+DEAL::Testing: FE_SimplexP<2>(3)
+DEAL::Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok 
+DEAL::Testing: FE_SimplexDGP<2>(3)
+DEAL::Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok 
+DEAL::Testing: FE_SimplexP<3>(3)
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Testing: FE_SimplexDGP<3>(3)
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok 
+DEAL::Ok Ok Ok Ok Ok Ok Ok Ok Ok Ok