From 14887188a415ca016011d65ee72845a3fce4faa0 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 7 Dec 2021 18:08:58 +0100
Subject: [PATCH] MatrixFree: distinguish internally between mixed and
 non-mixed meshes

---
 .../deal.II/matrix_free/face_setup_internal.h | 35 +++++++++++++------
 .../matrix_free/mapping_info.templates.h      |  9 +++--
 2 files changed, 30 insertions(+), 14 deletions(-)

diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index 5499c18300..bcfd049961 100644
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -109,7 +109,8 @@ namespace internal
         const typename dealii::Triangulation<dim>::cell_iterator &cell,
         const unsigned int number_cell_interior,
         const typename dealii::Triangulation<dim>::cell_iterator &neighbor,
-        const unsigned int number_cell_exterior);
+        const unsigned int number_cell_exterior,
+        const bool         is_mixed_mesh);
 
       bool use_active_cells;
 
@@ -738,6 +739,8 @@ namespace internal
       const std::vector<std::pair<unsigned int, unsigned int>> &cell_levels,
       TaskInfo &                                                task_info)
     {
+      const bool is_mixed_mesh = triangulation.is_mixed_mesh();
+
       // step 1: create the inverse map between cell iterators and the
       // cell_level_index field
       std::map<std::pair<unsigned int, unsigned int>, unsigned int>
@@ -792,8 +795,10 @@ namespace internal
                         info.interior_face_no  = f;
                         info.exterior_face_no  = dcell->face(f)->boundary_id();
                         info.face_type =
-                          dcell->face(f)->reference_cell() !=
-                          dealii::ReferenceCells::get_hypercube<dim - 1>();
+                          is_mixed_mesh ?
+                            (dcell->face(f)->reference_cell() !=
+                             dealii::ReferenceCells::get_hypercube<dim - 1>()) :
+                            0;
                         info.subface_index =
                           GeometryInfo<dim>::max_children_per_cell;
                         info.face_orientation = 0;
@@ -842,7 +847,8 @@ namespace internal
                                           neighbor_c,
                                           map_to_vectorized[level_index],
                                           dcell,
-                                          cell));
+                                          cell,
+                                          is_mixed_mesh));
                                       }
                                     else if (face_is_owned[dcell->face(f)
                                                              ->child(c)
@@ -854,7 +860,8 @@ namespace internal
                                           neighbor_c,
                                           map_to_vectorized[level_index],
                                           dcell,
-                                          cell));
+                                          cell,
+                                          is_mixed_mesh));
                                       }
                                     else
                                       Assert(
@@ -901,7 +908,8 @@ namespace internal
                                       dcell,
                                       cell,
                                       neighbor,
-                                      map_to_vectorized[level_index]));
+                                      map_to_vectorized[level_index],
+                                      is_mixed_mesh));
                                   }
                                 else if (face_is_owned[dcell->face(f)
                                                          ->index()] ==
@@ -912,7 +920,8 @@ namespace internal
                                       dcell,
                                       cell,
                                       neighbor,
-                                      map_to_vectorized[level_index]));
+                                      map_to_vectorized[level_index],
+                                      is_mixed_mesh));
                                   }
                               }
                             else
@@ -933,7 +942,8 @@ namespace internal
                                               dcell,
                                               cell,
                                               neighbor,
-                                              refinement_edge_faces.size()));
+                                              refinement_edge_faces.size(),
+                                              is_mixed_mesh));
                               }
                           }
                       }
@@ -962,7 +972,8 @@ namespace internal
       const typename dealii::Triangulation<dim>::cell_iterator &cell,
       const unsigned int number_cell_interior,
       const typename dealii::Triangulation<dim>::cell_iterator &neighbor,
-      const unsigned int number_cell_exterior)
+      const unsigned int number_cell_exterior,
+      const bool         is_mixed_mesh)
     {
       FaceToCellTopology<1> info;
       info.cells_interior[0] = number_cell_interior;
@@ -973,8 +984,10 @@ namespace internal
       else
         info.exterior_face_no = cell->neighbor_face_no(face_no);
 
-      info.face_type = cell->face(face_no)->reference_cell() !=
-                       dealii::ReferenceCells::get_hypercube<dim - 1>();
+      info.face_type = is_mixed_mesh ?
+                         (cell->face(face_no)->reference_cell() !=
+                          dealii::ReferenceCells::get_hypercube<dim - 1>()) :
+                         0;
 
       info.subface_index = GeometryInfo<dim>::max_children_per_cell;
       Assert(neighbor->level() <= cell->level(), ExcInternalError());
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index d5c35a4972..fc29f2a9d4 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -124,6 +124,8 @@ namespace internal
       face_data.resize(quad.size());
       face_data_by_cells.resize(quad.size());
 
+      const bool is_mixed_mesh = tria.is_mixed_mesh();
+
       // dummy FE that is used to set up an FEValues object. Do not need the
       // actual finite element because we will only evaluate quantities for
       // the mapping that are independent of the FE
@@ -185,10 +187,11 @@ namespace internal
                   if (quad_face.second.size() > 0) // triangle
                     {
                       AssertDimension(dim, 3);
-                      face_data[my_q].descriptor[hpq * scale + 1].initialize(
-                        quad_face.second, update_default);
+                      face_data[my_q]
+                        .descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)]
+                        .initialize(quad_face.second, update_default);
                       face_data_by_cells[my_q]
-                        .descriptor[hpq * scale + 1]
+                        .descriptor[hpq * scale + (is_mixed_mesh ? 1 : 0)]
                         .initialize(quad_face.second, update_default);
                     }
 
-- 
2.39.5