MatrixFree: store refinement configuration once per cell

diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h
index 7b0583ea9e9af29397fc1b16ca9f06dae4cac215..1dc6024ff10ebec3891b1d0a3b104779aecef6e6 100644 (file)
--- a/include/deal.II/matrix_free/dof_info.h
+++ b/include/deal.II/matrix_free/dof_info.h
@@ -526,6 +526,12 @@ namespace internal
        */
       std::vector<unsigned int> dof_indices;
 
+      /**
+       * Supported components of all entries of the hp::FECollection object of
+       * the given DoFHandler.
+       */
+      std::vector<std::vector<bool>> hanging_node_constraint_masks_comp;
+
       /**
        * Masks indicating for each cell and component if the optimized
        * hanging-node constraint is applicable and if yes which type.

diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h
index 63ff5f22a39dee0b3cba881744355a740012cc3e..9e3571a1544cf400e1f5fe11327dea035258c72f 100644 (file)
--- a/include/deal.II/matrix_free/dof_info.templates.h
+++ b/include/deal.II/matrix_free/dof_info.templates.h
@@ -283,13 +283,38 @@ namespace internal
       const TriaIterator<DoFCellAccessor<dim, dim, false>> &cell,
       std::vector<types::global_dof_index> &                dof_indices)
     {
-      const ArrayView<ConstraintKinds> mask_view(
-        hanging_node_constraint_masks.data() +
-          cell_number * cell->get_fe().n_components(),
-        cell->get_fe().n_components());
-
-      return hanging_nodes.setup_constraints(
-        cell, {}, lexicographic_mapping, dof_indices, mask_view);
+      // 1) check if finite elements support fast hanging-node algorithm
+      this->hanging_node_constraint_masks_comp =
+        hanging_nodes.compute_supported_components(
+          cell->get_dof_handler().get_fe_collection());
+
+      if ([](const auto &supported_components) {
+            return std::none_of(supported_components.begin(),
+                                supported_components.end(),
+                                [](const auto &a) {
+                                  return *std::max_element(a.begin(), a.end());
+                                });
+          }(hanging_node_constraint_masks_comp))
+        return false;
+
+      // 2) determine the refinement configuration of the cell
+      const auto refinement_configuration =
+        hanging_nodes.compute_refinement_configuration(cell);
+
+      if (refinement_configuration == ConstraintKinds::unconstrained)
+        return false;
+
+      // 3) update DoF indices of cell for specified components
+      hanging_nodes.update_dof_indices(cell,
+                                       {},
+                                       lexicographic_mapping,
+                                       hanging_node_constraint_masks_comp,
+                                       refinement_configuration,
+                                       dof_indices);
+
+      hanging_node_constraint_masks[cell_number] = refinement_configuration;
+
+      return true;
     }
 
 

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 647edf920880123aab6362d3c79cb48045172ce8..e5f58db72788b1e8b2901a08677bcc9d39938be4 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3325,10 +3325,10 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       for (unsigned int v = 0; v < n_lanes; ++v)
         if (cells[v] != numbers::invalid_unsigned_int &&
             this->dof_info->hanging_node_constraint_masks.size() > 0 &&
-            this->dof_info->hanging_node_constraint_masks
-                [cells[v] * this->n_fe_components +
-                 this->first_selected_component] !=
-              internal::MatrixFreeFunctions::ConstraintKinds::unconstrained)
+            this->dof_info->hanging_node_constraint_masks[cells[v]] !=
+              internal::MatrixFreeFunctions::ConstraintKinds::unconstrained &&
+            this->dof_info->hanging_node_constraint_masks_comp
+              [this->active_fe_index][this->first_selected_component])
           has_hn_constraints = true;
     }
 
@@ -3427,10 +3427,10 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             has_constraints = true;
 
           if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
-              this->dof_info->hanging_node_constraint_masks
-                  [cells[v] * this->n_fe_components +
-                   this->first_selected_component] !=
-                internal::MatrixFreeFunctions::ConstraintKinds::unconstrained)
+              this->dof_info->hanging_node_constraint_masks[cells[v]] !=
+                internal::MatrixFreeFunctions::ConstraintKinds::unconstrained &&
+              this->dof_info->hanging_node_constraint_masks_comp
+                [this->active_fe_index][this->first_selected_component])
             has_hn_constraints = true;
 
           Assert(my_index_start[n_components_read].first ==
@@ -3514,9 +3514,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           (this->dof_info->row_starts[cell_dof_index].second !=
              this->dof_info->row_starts[cell_dof_index + n_components_read]
                .second ||
-           (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
-            this->dof_info->hanging_node_constraint_masks[cell_dof_index] !=
-              internal::MatrixFreeFunctions::ConstraintKinds::unconstrained)))
+           ((this->dof_info->hanging_node_constraint_masks.size() > 0 &&
+             this->dof_info->hanging_node_constraint_masks[cell_index] !=
+               internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) &&
+            this->dof_info->hanging_node_constraint_masks_comp
+              [this->active_fe_index][this->first_selected_component])))
         {
           Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
                    numbers::invalid_unsigned_int,
@@ -4152,7 +4154,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   apply_hanging_node_constraints(const bool transpose) const
 {
   if (this->dof_info == nullptr ||
-      this->dof_info->hanging_node_constraint_masks.size() == 0)
+      this->dof_info->hanging_node_constraint_masks.size() == 0 ||
+      this->dof_info->hanging_node_constraint_masks_comp
+          [this->active_fe_index][this->first_selected_component] == false)
     return; // nothing to do with faces
 
   constexpr unsigned int n_lanes = VectorizedArrayType::size();
@@ -4174,11 +4178,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         }
 
       const unsigned int cell_index = cells[v];
-      const unsigned int cell_dof_index =
-        cell_index * this->n_fe_components + this->first_selected_component;
-
-      const auto mask =
-        this->dof_info->hanging_node_constraint_masks[cell_dof_index];
+      const auto         mask =
+        this->dof_info->hanging_node_constraint_masks[cell_index];
       constraint_mask[v] = mask;
 
       hn_available |=

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 7a341ee05e227edef5ff3d1ec9cc75c1b7382abd..f254e9ed87055fadf1b1df11b7d4eba7b66f95c6 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -1236,8 +1236,7 @@ namespace internal
         hanging_nodes = std::make_unique<
           dealii::internal::MatrixFreeFunctions::HangingNodes<dim>>(tria);
         for (unsigned int no = 0; no < n_dof_handlers; ++no)
-          dof_info[no].hanging_node_constraint_masks.resize(
-            n_active_cells * dof_handler[no]->get_fe().n_components());
+          dof_info[no].hanging_node_constraint_masks.resize(n_active_cells);
       }
 
     for (unsigned int counter = 0; counter < n_active_cells; ++counter)

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index a2429989e6f91e9aaec1eef25d9aea50b9bc0abe..39d52486a36f173cc54a2bc6e9423374cb40e865 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -447,13 +447,13 @@ namespace MatrixFreeTools
               locally_relevant_constrains.end());
 
             // STEP 2c: apply hanging-node constraints
-            if (dof_info.hanging_node_constraint_masks.size() > 0)
+            if (dof_info.hanging_node_constraint_masks.size() > 0 &&
+                dof_info
+                  .hanging_node_constraint_masks_comp[phi.get_active_fe_index()]
+                                                     [first_selected_component])
               {
                 const auto mask =
-                  dof_info
-                    .hanging_node_constraint_masks[(cell * n_lanes + v) *
-                                                     n_fe_components +
-                                                   first_selected_component];
+                  dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
 
                 // cell has hanging nodes
                 if (mask != dealii::internal::MatrixFreeFunctions::

diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc
index 6bd7669f1149892bda8ffa202e7a65ad765d2946..c9e2717eab504fe0b33d30ed85bd9560a9a1fbc5 100644 (file)
--- a/source/matrix_free/dof_info.cc
+++ b/source/matrix_free/dof_info.cc
@@ -104,8 +104,9 @@ namespace internal
           // one
           const bool has_constraints =
             (hanging_node_constraint_masks.size() != 0 &&
-             hanging_node_constraint_masks[ib] !=
-               ConstraintKinds::unconstrained) ||
+             hanging_node_constraint_masks[cell * n_vectorization + v] !=
+               ConstraintKinds::unconstrained &&
+             hanging_node_constraint_masks_comp[fe_index][0 /*TODO*/]) ||
             (row_starts[ib].second != row_starts[ib + n_fe_components].second);
 
           auto do_copy = [&](const unsigned int *begin,
@@ -226,13 +227,19 @@ namespace internal
                 {
                   bool has_hanging_nodes = false;
 
-                  if (hanging_node_constraint_masks.size() > 0)
+                  const unsigned int fe_index =
+                    (cell_active_fe_index.size() == 0 ||
+                     dofs_per_cell.size() == 1) ?
+                      0 :
+                      cell_active_fe_index[boundary_cells[i]];
+                  AssertIndexRange(fe_index, dofs_per_cell.size());
+
+                  if (hanging_node_constraint_masks.size() > 0 &&
+                      hanging_node_constraint_masks[boundary_cells[i]] !=
+                        ConstraintKinds::unconstrained)
                     for (unsigned int comp = 0; comp < n_components; ++comp)
                       has_hanging_nodes |=
-                        hanging_node_constraint_masks[boundary_cells[i] *
-                                                        n_components +
-                                                      comp] !=
-                        ConstraintKinds::unconstrained;
+                        hanging_node_constraint_masks_comp[fe_index][comp];
 
                   if (has_hanging_nodes ||
                       row_starts[boundary_cells[i] * n_components].second !=
@@ -242,12 +249,6 @@ namespace internal
                       unsigned int *data_ptr =
                         plain_dof_indices.data() +
                         row_starts_plain_indices[boundary_cells[i]];
-                      const unsigned int fe_index =
-                        (cell_active_fe_index.size() == 0 ||
-                         dofs_per_cell.size() == 1) ?
-                          0 :
-                          cell_active_fe_index[boundary_cells[i]];
-                      AssertIndexRange(fe_index, dofs_per_cell.size());
                       const unsigned int *row_end =
                         data_ptr + dofs_per_cell[fe_index];
                       for (; data_ptr != row_end; ++data_ptr)
@@ -340,7 +341,7 @@ namespace internal
 
       std::vector<ConstraintKinds> new_hanging_node_constraint_masks;
       new_hanging_node_constraint_masks.reserve(
-        new_hanging_node_constraint_masks.size());
+        hanging_node_constraint_masks.size());
 
       if (store_plain_indices == true)
         {
@@ -373,6 +374,19 @@ namespace internal
 
               bool has_hanging_nodes = false;
 
+              if (hanging_node_constraint_masks.size() > 0)
+                {
+                  const auto mask =
+                    hanging_node_constraint_masks[renumbering[position_cell +
+                                                              j]];
+                  new_hanging_node_constraint_masks.push_back(mask);
+
+                  if (mask != ConstraintKinds::unconstrained)
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      has_hanging_nodes |= hanging_node_constraint_masks_comp
+                        [have_hp ? cell_active_fe_index[i] : 0][comp];
+                }
+
               for (unsigned int comp = 0; comp < n_components; ++comp)
                 {
                   new_row_starts[(i * vectorization_length + j) * n_components +
@@ -382,15 +396,6 @@ namespace internal
                                  comp]
                     .second = new_constraint_indicator.size();
 
-                  if (hanging_node_constraint_masks.size() > 0)
-                    {
-                      const auto mask =
-                        hanging_node_constraint_masks[cell_no + comp];
-                      new_hanging_node_constraint_masks.push_back(mask);
-                      has_hanging_nodes |=
-                        mask != ConstraintKinds::unconstrained;
-                    }
-
                   new_dof_indices.insert(
                     new_dof_indices.end(),
                     dof_indices.data() + row_starts[cell_no + comp].first,
@@ -426,11 +431,13 @@ namespace internal
                 new_row_starts[(i * vectorization_length + j) * n_components +
                                comp]
                   .second = new_constraint_indicator.size();
-
-                if (hanging_node_constraint_masks.size() > 0)
-                  new_hanging_node_constraint_masks.push_back(
-                    ConstraintKinds::unconstrained);
               }
+
+          for (unsigned int j = n_vect; j < vectorization_length; ++j)
+            if (hanging_node_constraint_masks.size() > 0)
+              new_hanging_node_constraint_masks.push_back(
+                ConstraintKinds::unconstrained);
+
           position_cell += n_vect;
         }
       AssertDimension(position_cell * n_components + 1, row_starts.size());