Add MatrixFreeTools::ElementBirthAndDeathMatrixFree

doc/news/changes/minor/20220709Munch

@@ -0,0 +1,7 @@
+New: The new class MatrixFreeTools::ElementBirthAndDeathMatrixFree
+is a wrapper around MatrixFree, which helps users to deal with 
+DoFHandlers, where there are cells without degrees of freedom, 
+i.e., with FENothing assigned to them. The class helps to implement 
+"element birth and death technique".
+<br>
+(Peter Munch, Sebastian Proell, 2022/07/09)

include/deal.II/matrix_free/fe_evaluation.h

@@ -2731,6 +2731,22 @@ class FEFaceEvaluation : public FEEvaluationAccess<dim,
   std_cxx20::ranges::iota_view<unsigned int, unsigned int>
   dof_indices() const;
 
+  /**
+   * Return whether the face is at the boundary.
+   */
+  bool
+  at_boundary() const;
+
+  /**
+   * Return the boundary indicator of this object.
+   *
+   * If the return value is the special value
+   * numbers::internal_face_boundary_id, then this object is in the interior of
+   * the domain.
+   */
+  types::boundary_id
+  boundary_id() const;
+
   /**
    * The number of degrees of freedom of a single component on the cell for
    * the underlying evaluation object. Usually close to
@@ -8881,6 +8897,63 @@ FEFaceEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+bool
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::at_boundary() const
+{
+  Assert(this->dof_access_index !=
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
+         ExcNotImplemented());
+
+  if (this->is_interior_face() == false)
+    return false;
+  else if (this->cell < this->matrix_free->n_inner_face_batches())
+    return false;
+  else if (this->cell < (this->matrix_free->n_inner_face_batches() +
+                         this->matrix_free->n_boundary_face_batches()))
+    return true;
+  else
+    return false;
+}
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+types::boundary_id
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::boundary_id() const
+{
+  Assert(this->dof_access_index !=
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell,
+         ExcNotImplemented());
+
+  if (at_boundary())
+    return this->matrix_free->get_boundary_id(this->cell);
+  else
+    return numbers::internal_face_boundary_id;
+}
+
+
+
 /*------------------------- end FEFaceEvaluation ------------------------- */
 
 

include/deal.II/matrix_free/tools.h

@@ -157,6 +157,183 @@ namespace MatrixFreeTools
     const unsigned int first_selected_component = 0);
 
 
+
+  /**
+   * A wrapper around MatrixFree to help users to deal with DoFHandlers,
+   * where there are cells without degrees of freedom, i.e., with
+   * FENothing assigned to them. In the following we call such cells
+   * invalid. All other cells are valid. In contrast to MatrixFree,
+   * this class skips invalid cells and faces between valid and
+   * invalid cells are treated as boundary faces.
+   *
+   * The class helps to implement "element birth and death technique",
+   * in which context valid cells are alive and invalid cells are dead.
+   */
+  template <int dim,
+            typename Number,
+            typename VectorizedArrayType = VectorizedArray<Number>>
+  class ElementBirthAndDeathMatrixFree
+  {
+  public:
+    /**
+     * Struct that helps to configure ElementBirthAndDeathMatrixFree.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor.
+       */
+      AdditionalData(const unsigned int dof_index = 0)
+        : dof_index(dof_index)
+      {}
+
+      /**
+       * Index of the DoFHandler within MatrixFree to be used.
+       */
+      unsigned int dof_index;
+    };
+
+    /**
+     * Reinitialize class based on a given MatrixFree instance.
+     * Particularly, the index of the valid FE is determined.
+     */
+    void
+    reinit(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+           const AdditionalData &additional_data = AdditionalData())
+    {
+      this->matrix_free = &matrix_free;
+
+      std::set<unsigned int> valid_fe_indices;
+
+      const auto &fe_collection =
+        matrix_free.get_dof_handler(additional_data.dof_index)
+          .get_fe_collection();
+
+      for (unsigned int i = 0; i < fe_collection.size(); ++i)
+        if (fe_collection[i].n_dofs_per_cell() > 0)
+          valid_fe_indices.insert(i);
+
+      AssertDimension(valid_fe_indices.size(), 1);
+
+      fe_index_valid = *valid_fe_indices.begin();
+    }
+
+    /**
+     * Loop over all valid cells.
+     *
+     * For the meaning of the parameters see MatrixFree::cell_loop().
+     */
+    template <typename VectorTypeOut, typename VectorTypeIn>
+    void
+    cell_loop(const std::function<void(
+                const MatrixFree<dim, Number, VectorizedArrayType> &,
+                VectorTypeOut &,
+                const VectorTypeIn &,
+                const std::pair<unsigned int, unsigned int>)> &cell_operation,
+              VectorTypeOut &                                  dst,
+              const VectorTypeIn &                             src,
+              const bool zero_dst_vector = false)
+    {
+      const auto ebd_cell_operation = [&](const auto &matrix_free,
+                                          auto &      dst,
+                                          const auto &src,
+                                          const auto  range) {
+        const auto category = matrix_free.get_cell_range_category(range);
+
+        if (category != fe_index_valid)
+          return;
+
+        cell_operation(matrix_free, dst, src, range);
+      };
+
+      matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
+        ebd_cell_operation, dst, src, zero_dst_vector);
+    }
+
+    /**
+     * Loop over all valid cells and faces. Faces between valid and
+     * invalid cells are treated as boundary faces with the boundary ID
+     * numbers::internal_face_boundary_id.
+     *
+     * For the meaning of the parameters see MatrixFree::cell_loop().
+     */
+    template <typename VectorTypeOut, typename VectorTypeIn>
+    void
+    loop(const std::function<
+           void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                VectorTypeOut &,
+                const VectorTypeIn &,
+                const std::pair<unsigned int, unsigned int>)> &cell_operation,
+         const std::function<
+           void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                VectorTypeOut &,
+                const VectorTypeIn &,
+                const std::pair<unsigned int, unsigned int>)> &face_operation,
+         const std::function<
+           void(const MatrixFree<dim, Number, VectorizedArrayType> &,
+                VectorTypeOut &,
+                const VectorTypeIn &,
+                const std::pair<unsigned int, unsigned int>,
+                const bool)> &boundary_operation,
+         VectorTypeOut &      dst,
+         const VectorTypeIn & src,
+         const bool           zero_dst_vector = false)
+    {
+      const auto ebd_cell_operation = [&](const auto &matrix_free,
+                                          auto &      dst,
+                                          const auto &src,
+                                          const auto  range) {
+        const auto category = matrix_free.get_cell_range_category(range);
+
+        if (category != fe_index_valid)
+          return;
+
+        cell_operation(matrix_free, dst, src, range);
+      };
+
+      const auto ebd_internal_or_boundary_face_operation =
+        [&](const auto &matrix_free,
+            auto &      dst,
+            const auto &src,
+            const auto  range) {
+          const auto category = matrix_free.get_face_range_category(range);
+
+          const unsigned int type =
+            static_cast<unsigned int>(category.first == fe_index_valid) +
+            static_cast<unsigned int>(category.second == fe_index_valid);
+
+          if (type == 0)
+            return; // invalid face -> nothing to do
+
+          if (type == 1) // boundary face
+            boundary_operation(
+              matrix_free, dst, src, range, category.first == fe_index_valid);
+          else if (type == 2) // internal face
+            face_operation(matrix_free, dst, src, range);
+        };
+
+      matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
+        ebd_cell_operation,
+        ebd_internal_or_boundary_face_operation,
+        ebd_internal_or_boundary_face_operation,
+        dst,
+        src,
+        zero_dst_vector);
+    }
+
+  private:
+    /**
+     * Reference to the underlying MatrixFree object.
+     */
+    SmartPointer<const MatrixFree<dim, Number, VectorizedArrayType>>
+      matrix_free;
+
+    /**
+     * Index of the valid FE. Currently, ony a single one is supported.
+     */
+    unsigned int fe_index_valid;
+  };
+
   // implementations
 
 #ifndef DOXYGEN