FEInterfaceValues: fill indices for TriaIterator<CellAccessor>

diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h
index 9093cb00d8a4b00aced320448d4420bb6bdbc27f..3b0b2371f7f386356b6b67ba933b29c34ace6bf2 100644 (file)
--- a/include/deal.II/fe/fe_interface_values.h
+++ b/include/deal.II/fe/fe_interface_values.h
@@ -2339,6 +2339,9 @@ FEInterfaceValues<dim, spacedim>::reinit(
     std::is_same_v<DoFCellAccessor<dim, spacedim, false>,
                    typename CellNeighborIteratorType::AccessorType>;
 
+  unsigned int active_fe_index          = 0;
+  unsigned int active_fe_index_neighbor = 0;
+
   if (internal_fe_face_values)
     {
       if (sub_face_no == numbers::invalid_unsigned_int)
@@ -2373,8 +2376,8 @@ FEInterfaceValues<dim, spacedim>::reinit(
     }
   else if (internal_hp_fe_face_values)
     {
-      unsigned int active_fe_index = fe_index_in;
-      unsigned int active_fe_index_neighbor =
+      active_fe_index = fe_index_in;
+      active_fe_index_neighbor =
         (fe_index_neighbor_in != numbers::invalid_unsigned_int) ?
           fe_index_neighbor_in :
           fe_index_in;
@@ -2421,37 +2424,42 @@ FEInterfaceValues<dim, spacedim>::reinit(
 
       // Third check, if the above did not already suffice. We see if we
       // can get somewhere via the dominated's finite element index.
-      const unsigned int dominated_fe_index =
-        ((used_q_index == numbers::invalid_unsigned_int) ||
-             (used_mapping_index == numbers::invalid_unsigned_int) ?
-           internal_hp_fe_face_values->get_fe_collection().find_dominated_fe(
-             {active_fe_index, active_fe_index_neighbor}) :
-           numbers::invalid_unsigned_int);
-
-      if (used_q_index == numbers::invalid_unsigned_int)
-        {
-          Assert(dominated_fe_index != numbers::invalid_fe_index,
-                 ExcMessage(
-                   "You called this function with 'q_index' left at its "
-                   "default value, but this can only work if one of "
-                   "the two finite elements adjacent to this face "
-                   "dominates the other. See the documentation "
-                   "of this function for more information of how "
-                   "to deal with this situation."));
-          used_q_index = dominated_fe_index;
-        }
-
-      if (used_mapping_index == numbers::invalid_unsigned_int)
+      if ((used_q_index == numbers::invalid_unsigned_int) ||
+          (used_mapping_index == numbers::invalid_unsigned_int))
         {
-          Assert(dominated_fe_index != numbers::invalid_fe_index,
-                 ExcMessage(
-                   "You called this function with 'mapping_index' left "
-                   "at its default value, but this can only work if one "
-                   "of the two finite elements adjacent to this face "
-                   "dominates the other. See the documentation "
-                   "of this function for more information of how "
-                   "to deal with this situation."));
-          used_mapping_index = dominated_fe_index;
+          const unsigned int dominated_fe_index =
+            ((used_q_index == numbers::invalid_unsigned_int) ||
+                 (used_mapping_index == numbers::invalid_unsigned_int) ?
+               internal_hp_fe_face_values->get_fe_collection()
+                 .find_dominated_fe(
+                   {active_fe_index, active_fe_index_neighbor}) :
+               numbers::invalid_unsigned_int);
+
+          if (used_q_index == numbers::invalid_unsigned_int)
+            {
+              Assert(dominated_fe_index != numbers::invalid_fe_index,
+                     ExcMessage(
+                       "You called this function with 'q_index' left at its "
+                       "default value, but this can only work if one of "
+                       "the two finite elements adjacent to this face "
+                       "dominates the other. See the documentation "
+                       "of this function for more information of how "
+                       "to deal with this situation."));
+              used_q_index = dominated_fe_index;
+            }
+
+          if (used_mapping_index == numbers::invalid_unsigned_int)
+            {
+              Assert(dominated_fe_index != numbers::invalid_fe_index,
+                     ExcMessage(
+                       "You called this function with 'mapping_index' left "
+                       "at its default value, but this can only work if one "
+                       "of the two finite elements adjacent to this face "
+                       "dominates the other. See the documentation "
+                       "of this function for more information of how "
+                       "to deal with this situation."));
+              used_mapping_index = dominated_fe_index;
+            }
         }
 
       // Same as if above, but when hp is enabled.
@@ -2550,6 +2558,28 @@ FEInterfaceValues<dim, spacedim>::reinit(
           ++idx;
         }
     }
+  else
+    {
+      const unsigned int n_dofs_per_cell_1 = fe_face_values->dofs_per_cell;
+      const unsigned int n_dofs_per_cell_2 =
+        fe_face_values_neighbor->dofs_per_cell;
+
+      interface_dof_indices.resize(n_dofs_per_cell_1 + n_dofs_per_cell_2);
+      dofmap.resize(n_dofs_per_cell_1 + n_dofs_per_cell_2);
+
+      for (unsigned int i = 0; i < n_dofs_per_cell_1; ++i)
+        {
+          interface_dof_indices[i] = numbers::invalid_dof_index;
+          dofmap[i]                = {i, numbers::invalid_unsigned_int};
+        }
+
+      for (unsigned int i = 0; i < n_dofs_per_cell_2; ++i)
+        {
+          interface_dof_indices[i + n_dofs_per_cell_1] =
+            numbers::invalid_dof_index;
+          dofmap[i + n_dofs_per_cell_1] = {numbers::invalid_unsigned_int, i};
+        }
+    }
 }
 
 
@@ -2589,30 +2619,24 @@ FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
       fe_face_values_neighbor = nullptr;
     }
 
+  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+
   if constexpr (std::is_same_v<typename CellIteratorType::AccessorType,
                                DoFCellAccessor<dim, spacedim, true>> ||
                 std::is_same_v<typename CellIteratorType::AccessorType,
                                DoFCellAccessor<dim, spacedim, false>>)
     {
-      if (internal_fe_face_values)
-        {
-          interface_dof_indices.resize(
-            fe_face_values->get_fe().n_dofs_per_cell());
-          cell->get_active_or_mg_dof_indices(interface_dof_indices);
-        }
-      else if (internal_hp_fe_face_values)
-        {
-          interface_dof_indices.resize(
-            fe_face_values->get_fe().n_dofs_per_cell());
-          cell->get_active_or_mg_dof_indices(interface_dof_indices);
-        }
-
-      dofmap.resize(interface_dof_indices.size());
-      for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
-        {
-          dofmap[i] = {{i, numbers::invalid_unsigned_int}};
-        }
+      cell->get_active_or_mg_dof_indices(interface_dof_indices);
     }
+  else
+    {
+      for (auto &i : interface_dof_indices)
+        i = numbers::invalid_dof_index;
+    }
+
+  dofmap.resize(interface_dof_indices.size());
+  for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
+    dofmap[i] = {{i, numbers::invalid_unsigned_int}};
 }
 
 

diff --git a/tests/feinterface/fe_interface_values_02.cc b/tests/feinterface/fe_interface_values_02.cc
index 8d62a0831376c1990e98981400f8192aa7f1dabf..4712f3e8b92ed388af9763ba4e73a3f5937ec020 100644 (file)
--- a/tests/feinterface/fe_interface_values_02.cc
+++ b/tests/feinterface/fe_interface_values_02.cc
@@ -32,31 +32,10 @@
 
 #include "../test_grids.h"
 
-
-template <int dim>
+template <int dim, typename IteratorType>
 void
-test(unsigned int fe_degree)
+test(const IteratorType &cell, FEInterfaceValues<dim> &fiv)
 {
-  Triangulation<dim> tria;
-  TestGrids::hyper_line(tria, 2);
-
-  DoFHandler<dim> dofh(tria);
-  FE_DGQ<dim>     fe(fe_degree);
-  deallog << fe.get_name() << std::endl;
-  dofh.distribute_dofs(fe);
-
-  MappingQ<dim> mapping(1);
-  UpdateFlags   update_flags = update_values | update_gradients |
-                             update_quadrature_points | update_JxW_values;
-
-  FEInterfaceValues<dim> fiv(mapping,
-                             fe,
-                             QGauss<dim - 1>(fe.degree + 1),
-                             update_flags);
-
-
-  auto cell = dofh.begin();
-
   for (const unsigned int f : GeometryInfo<dim>::face_indices())
     if (!cell->at_boundary(f))
       {
@@ -103,6 +82,32 @@ test(unsigned int fe_degree)
 }
 
 
+template <int dim>
+void
+test(unsigned int fe_degree)
+{
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(fe_degree);
+  deallog << fe.get_name() << std::endl;
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(mapping,
+                             fe,
+                             QGauss<dim - 1>(fe.degree + 1),
+                             update_flags);
+
+  test(tria.begin(), fiv);
+  test(dofh.begin(), fiv);
+}
+
+
 
 int
 main()

diff --git a/tests/feinterface/fe_interface_values_02.output b/tests/feinterface/fe_interface_values_02.output
index fe3840c85501f13500b3fbebf4ac5195c4c505ba..7360480b1234e7404737ee4158f677df52443078 100644 (file)
--- a/tests/feinterface/fe_interface_values_02.output
+++ b/tests/feinterface/fe_interface_values_02.output
@@ -4,18 +4,34 @@ DEAL::shape_value(true): 1.00000 0.00000
 DEAL::shape_value(false): 0.00000 1.00000
 DEAL::jump_in_shape_values(): 1.00000 -1.00000
 DEAL::average_of_shape_values(): 0.500000 0.500000
+DEAL::shape_value(true): 1.00000 0.00000
+DEAL::shape_value(false): 0.00000 1.00000
+DEAL::jump_in_shape_values(): 1.00000 -1.00000
+DEAL::average_of_shape_values(): 0.500000 0.500000
 DEAL::FE_DGQ<2>(1)
 DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000
 DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.500000 0.00000
 DEAL::jump_in_shape_values(): 0.00000 0.500000 0.00000 0.500000 -0.500000 0.00000 -0.500000 0.00000
 DEAL::average_of_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000
+DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.500000 0.00000
+DEAL::jump_in_shape_values(): 0.00000 0.500000 0.00000 0.500000 -0.500000 0.00000 -0.500000 0.00000
+DEAL::average_of_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000
 DEAL::FE_DGQ<3>(0)
 DEAL::shape_value(true): 1.00000 0.00000
 DEAL::shape_value(false): 0.00000 1.00000
 DEAL::jump_in_shape_values(): 1.00000 -1.00000
 DEAL::average_of_shape_values(): 0.500000 0.500000
+DEAL::shape_value(true): 1.00000 0.00000
+DEAL::shape_value(false): 0.00000 1.00000
+DEAL::jump_in_shape_values(): 1.00000 -1.00000
+DEAL::average_of_shape_values(): 0.500000 0.500000
 DEAL::FE_DGQ<3>(1)
 DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
 DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000
 DEAL::jump_in_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000
 DEAL::average_of_shape_values(): 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000
+DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000
+DEAL::jump_in_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000
+DEAL::average_of_shape_values(): 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000