MatrixFree: Element-centric loops with different cell geometry compressions

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 798764922413794dc917639896e0dff226f54b98..84c26effb449f5c0f0ca292b0d420cd6dc203fe8 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -8189,8 +8189,9 @@ FEFaceEvaluation<dim,
     return;
   Assert(this->matrix_free != nullptr, ExcNotInitialized());
 
-  this->cell_type = this->matrix_free->get_mapping_info().cell_type[cell_index];
-  this->cell      = cell_index;
+  this->cell_type = this->matrix_free->get_mapping_info()
+                      .faces_by_cells_type[cell_index][face_number];
+  this->cell          = cell_index;
   this->subface_index = GeometryInfo<dim>::max_children_per_cell;
   this->dof_access_index =
     internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index 7087e5fda31a381abd4faa544d45361ce5d62ebc..bd46a1ce5fca8379968f49814fb9abaed5006178 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -67,7 +67,7 @@ namespace internal
       initialize(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-        const FaceInfo<VectorizedArrayType::size()> &             faces,
+        const FaceInfo<VectorizedArrayType::size()> &             face_info,
         const std::vector<unsigned int> &active_fe_index,
         const std::shared_ptr<dealii::hp::MappingCollection<dim>> &mapping,
         const std::vector<dealii::hp::QCollection<dim>> &          quad,
@@ -86,7 +86,7 @@ namespace internal
       update_mapping(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-        const FaceInfo<VectorizedArrayType::size()> &             faces,
+        const FaceInfo<VectorizedArrayType::size()> &             face_info,
         const std::vector<unsigned int> &active_fe_index,
         const std::shared_ptr<dealii::hp::MappingCollection<dim>> &mapping);
 
@@ -157,6 +157,15 @@ namespace internal
        */
       std::vector<GeometryType> face_type;
 
+      /**
+       * Store whether the face geometry for the face_data_by_cells data
+       * structure is represented as Cartesian (cell type 0), with constant
+       * transform data (Jacobians) (cell type 1), or a general type (cell
+       * type 3). Note that both the interior and exterior agree on the type
+       * of the data structure, using the more general of the two.
+       */
+      std::vector<std::array<GeometryType, 2 * dim>> faces_by_cells_type;
+
       /**
        * The data cache for the cells.
        */
@@ -209,15 +218,14 @@ namespace internal
        * given as a tuple of the level and index within the level as used in
        * the main initialization of the class
        *
-       * @param faces The description of the connectivity from faces to cells
-       * as filled in the MatrixFree class
+       * @param face_info The description of the connectivity from faces to
+       * cells as filled in the MatrixFree class
        */
       void
       compute_mapping_q(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
-        const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
-          &faces);
+        const FaceInfo<VectorizedArrayType::size()> &             face_info);
 
       /**
        * Computes the information in the given cells, called within
@@ -251,6 +259,7 @@ namespace internal
       initialize_faces_by_cells(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
+        const FaceInfo<VectorizedArrayType::size()> &             face_info,
         const dealii::hp::MappingCollection<dim> &                mapping);
     };
 

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 5aa4fb340922171f829a1ebe9d65ce99fe1efd20..e9178bb5e2aedce1175a2907542522bef104dee6 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -181,7 +181,7 @@ namespace internal
       // use the fast method.
       if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 &&
           dynamic_cast<const MappingQ<dim> *>(&mapping->operator[](0)))
-        compute_mapping_q(tria, cells, face_info.faces);
+        compute_mapping_q(tria, cells, face_info);
       else
         {
           // Could call these functions in parallel, but not useful because
@@ -189,7 +189,7 @@ namespace internal
           initialize_cells(tria, cells, active_fe_index, *mapping);
           initialize_faces(
             tria, cells, face_info.faces, active_fe_index, *mapping);
-          initialize_faces_by_cells(tria, cells, *mapping);
+          initialize_faces_by_cells(tria, cells, face_info, *mapping);
         }
     }
 
@@ -219,7 +219,7 @@ namespace internal
 
       if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 &&
           dynamic_cast<const MappingQ<dim> *>(&mapping->operator[](0)))
-        compute_mapping_q(tria, cells, face_info.faces);
+        compute_mapping_q(tria, cells, face_info);
       else
         {
           // Could call these functions in parallel, but not useful because
@@ -227,7 +227,7 @@ namespace internal
           initialize_cells(tria, cells, active_fe_index, *mapping);
           initialize_faces(
             tria, cells, face_info.faces, active_fe_index, *mapping);
-          initialize_faces_by_cells(tria, cells, *mapping);
+          initialize_faces_by_cells(tria, cells, face_info, *mapping);
         }
     }
 
@@ -2530,7 +2530,7 @@ namespace internal
     MappingInfo<dim, Number, VectorizedArrayType>::compute_mapping_q(
       const dealii::Triangulation<dim> &                        tria,
       const std::vector<std::pair<unsigned int, unsigned int>> &cell_array,
-      const std::vector<FaceToCellTopology<VectorizedArrayType::size()>> &faces)
+      const FaceInfo<VectorizedArrayType::size()> &             face_info)
     {
       // step 1: extract quadrature point data with the data appropriate for
       // MappingQ
@@ -2714,6 +2714,8 @@ namespace internal
                      std::size_t(2U)));
         }
 
+      const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
+        &faces = face_info.faces;
       if (faces.empty())
         return;
 
@@ -2882,7 +2884,10 @@ namespace internal
       // transitioned to extracting the information from cell quadrature
       // points but we need to figure out the correct indices of neighbors
       // within the list of arrays still
-      initialize_faces_by_cells(tria, cell_array, *this->mapping_collection);
+      initialize_faces_by_cells(tria,
+                                cell_array,
+                                face_info,
+                                *this->mapping_collection);
     }
 
 
@@ -2892,6 +2897,7 @@ namespace internal
     MappingInfo<dim, Number, VectorizedArrayType>::initialize_faces_by_cells(
       const dealii::Triangulation<dim> &                        tria,
       const std::vector<std::pair<unsigned int, unsigned int>> &cells,
+      const FaceInfo<VectorizedArrayType::size()> &             face_info,
       const dealii::hp::MappingCollection<dim> &                mapping_in)
     {
       if (update_flags_faces_by_cells == update_default)
@@ -2901,14 +2907,49 @@ namespace internal
       AssertDimension(mapping_in.size(), 1);
       const auto &mapping = mapping_in[0];
 
-      const unsigned int n_quads = face_data_by_cells.size();
-      const unsigned int n_lanes = VectorizedArrayType::size();
-      UpdateFlags        update_flags =
+      const unsigned int     n_quads = face_data_by_cells.size();
+      constexpr unsigned int n_lanes = VectorizedArrayType::size();
+      UpdateFlags            update_flags =
         (update_flags_faces_by_cells & update_quadrature_points ?
            update_quadrature_points :
            update_default) |
         update_normal_vectors | update_JxW_values | update_jacobians;
 
+      const auto compute_neighbor_index =
+        [&face_info](const unsigned int cell,
+                     const unsigned int face,
+                     const unsigned int lane) {
+          constexpr unsigned int n_lanes   = VectorizedArrayType::size();
+          const unsigned int     cell_this = cell * n_lanes + lane;
+          unsigned int           face_index =
+            face_info.cell_and_face_to_plain_faces(cell, face, lane);
+          if (face_index == numbers::invalid_unsigned_int)
+            return numbers::invalid_unsigned_int;
+
+          unsigned int cell_neighbor = face_info.faces[face_index / n_lanes]
+                                         .cells_interior[face_index % n_lanes];
+          if (cell_neighbor == cell_this)
+            cell_neighbor = face_info.faces[face_index / n_lanes]
+                              .cells_exterior[face_index % n_lanes];
+          return cell_neighbor;
+        };
+
+      faces_by_cells_type.resize(cell_type.size());
+      for (unsigned int cell = 0; cell < cell_type.size(); ++cell)
+        for (const unsigned int face : GeometryInfo<dim>::face_indices())
+          {
+            faces_by_cells_type[cell][face] = cell_type[cell];
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              {
+                const unsigned int cell_neighbor =
+                  compute_neighbor_index(cell, face, v);
+                if (cell_neighbor != numbers::invalid_unsigned_int)
+                  faces_by_cells_type[cell][face] =
+                    std::max(faces_by_cells_type[cell][face],
+                             cell_type[cell_neighbor / n_lanes]);
+              }
+          }
+
       for (unsigned int my_q = 0; my_q < n_quads; ++my_q)
         {
           // since we already know the cell type, we can pre-allocate the right
@@ -2924,7 +2965,7 @@ namespace internal
           for (unsigned int i = 0; i < cell_type.size(); ++i)
             for (const unsigned int face : GeometryInfo<dim>::face_indices())
               {
-                if (cell_type[i] <= affine)
+                if (faces_by_cells_type[i][face] <= affine)
                   {
                     face_data_by_cells[my_q].data_index_offsets
                       [i * GeometryInfo<dim>::faces_per_cell + face] =
@@ -3010,6 +3051,8 @@ namespace internal
                   .data_index_offsets[cell * GeometryInfo<dim>::faces_per_cell +
                                       face];
 
+              const GeometryType my_cell_type = faces_by_cells_type[cell][face];
+
               for (unsigned int v = 0; v < n_lanes; ++v)
                 {
                   typename dealii::Triangulation<dim>::cell_iterator cell_it(
@@ -3018,18 +3061,15 @@ namespace internal
                     cells[cell * n_lanes + v].second);
                   fe_val.reinit(cell_it, face);
 
-                  const bool is_local =
-                    (cell_it->is_active() ?
-                       cell_it->is_locally_owned() :
-                       cell_it->is_locally_owned_on_level()) &&
-                    (!cell_it->at_boundary(face) ||
-                     (cell_it->at_boundary(face) &&
-                      cell_it->has_periodic_neighbor(face)));
+                  const unsigned int cell_neighbor =
+                    compute_neighbor_index(cell, face, v);
 
-                  if (is_local)
+                  if (cell_neighbor != numbers::invalid_unsigned_int)
                     {
-                      auto cell_it_neigh =
-                        cell_it->neighbor_or_periodic_neighbor(face);
+                      typename dealii::Triangulation<dim>::cell_iterator
+                        cell_it_neigh(&tria,
+                                      cells[cell_neighbor].first,
+                                      cells[cell_neighbor].second);
                       fe_val_neigh.reinit(cell_it_neigh,
                                           cell_it->at_boundary(face) ?
                                             cell_it->periodic_neighbor_face_no(
@@ -3038,7 +3078,7 @@ namespace internal
                     }
 
                   // copy data for affine data type
-                  if (cell_type[cell] <= affine)
+                  if (my_cell_type <= affine)
                     {
                       if (update_flags & update_JxW_values)
                         face_data_by_cells[my_q].JxW_values[offset][v] =
@@ -3057,7 +3097,8 @@ namespace internal
                                   inv_jac[d][ee];
                               }
                         }
-                      if (is_local && (update_flags & update_jacobians))
+                      if (cell_neighbor != numbers::invalid_unsigned_int &&
+                          (update_flags & update_jacobians))
                         for (unsigned int q = 0; q < fe_val.n_quadrature_points;
                              ++q)
                           {
@@ -3109,7 +3150,8 @@ namespace internal
                                     inv_jac[d][ee];
                                 }
                           }
-                      if (is_local && (update_flags & update_jacobians))
+                      if (cell_neighbor != numbers::invalid_unsigned_int &&
+                          (update_flags & update_jacobians))
                         for (unsigned int q = 0; q < fe_val.n_quadrature_points;
                              ++q)
                           {
@@ -3149,9 +3191,9 @@ namespace internal
                 }
               if (update_flags & update_normal_vectors &&
                   update_flags & update_jacobians)
-                for (unsigned int q = 0; q < (cell_type[cell] <= affine ?
-                                                1 :
-                                                fe_val.n_quadrature_points);
+                for (unsigned int q = 0;
+                     q <
+                     (my_cell_type <= affine ? 1 : fe_val.n_quadrature_points);
                      ++q)
                   face_data_by_cells[my_q]
                     .normals_times_jacobians[0][offset + q] =
@@ -3159,9 +3201,9 @@ namespace internal
                     face_data_by_cells[my_q].jacobians[0][offset + q];
               if (update_flags & update_normal_vectors &&
                   update_flags & update_jacobians)
-                for (unsigned int q = 0; q < (cell_type[cell] <= affine ?
-                                                1 :
-                                                fe_val.n_quadrature_points);
+                for (unsigned int q = 0;
+                     q <
+                     (my_cell_type <= affine ? 1 : fe_val.n_quadrature_points);
                      ++q)
                   face_data_by_cells[my_q]
                     .normals_times_jacobians[1][offset + q] =

diff --git a/tests/matrix_free/compare_faces_by_cells.cc b/tests/matrix_free/compare_faces_by_cells.cc
index a30a4d3397430d0303dc0c820f44bf6ecb453523..8e23091004e103fa11aa74deadd0a9640eb17ae4 100644 (file)
--- a/tests/matrix_free/compare_faces_by_cells.cc
+++ b/tests/matrix_free/compare_faces_by_cells.cc
@@ -15,7 +15,7 @@
 
 
 
-// compares the computation of the diagonal using a face integration
+// compares the computation of the diagonal using the face integration
 // facilities with the alternative reinit(cell_index, face_number) approach
 
 #include <deal.II/base/function.h>
@@ -82,8 +82,8 @@ public:
   compute_diagonal_by_face(
     LinearAlgebra::distributed::Vector<number> &result) const
   {
-    int dummy;
-    result = 0;
+    int dummy = 0;
+    result    = 0;
     data.loop(&LaplaceOperator::local_diagonal_dummy,
               &LaplaceOperator::local_diagonal_face,
               &LaplaceOperator::local_diagonal_boundary,
@@ -96,7 +96,7 @@ public:
   compute_diagonal_by_cell(
     LinearAlgebra::distributed::Vector<number> &result) const
   {
-    int dummy;
+    int dummy = 0;
     result.zero_out_ghost_values();
     data.cell_loop(&LaplaceOperator::local_diagonal_by_cell,
                    this,