]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add MatrixFreeTools::ElementActivationAndDeactivationMatrixFree 14121/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 9 Jul 2022 08:37:56 +0000 (10:37 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 3 Sep 2022 07:05:31 +0000 (09:05 +0200)
doc/news/changes/minor/20220709Munch [new file with mode: 0644]
include/deal.II/matrix_free/tools.h
tests/matrix_free/element_birth_and_death_01.cc [new file with mode: 0644]
tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20220709Munch b/doc/news/changes/minor/20220709Munch
new file mode 100644 (file)
index 0000000..3160896
--- /dev/null
@@ -0,0 +1,7 @@
+New: The new class MatrixFreeTools::ElementActivationAndDeactivationMatrixFree
+is a wrapper around MatrixFree designed to deal with 
+DoFHandler objects involving cells without degrees of freedom, i.e.,
+cells using FENothing as element type. The class helps to implement the
+"element birth and death technique".
+<br>
+(Peter Munch, Sebastian Proell, 2022/07/09)
index c5b34249d0007aadf01ca7e637183b2179ad7798..61082de1bd006885ecd29bf0ac9ad4dcbfa5f59e 100644 (file)
@@ -157,6 +157,188 @@ namespace MatrixFreeTools
     const unsigned int first_selected_component = 0);
 
 
+
+  /**
+   * A wrapper around MatrixFree to help users to deal with DoFHandler
+   * objects involving cells without degrees of freedom, i.e.,
+   * cells using FENothing as element type.  In the following we call such cells
+   * deactivated. All other cells are activated. In contrast to MatrixFree,
+   * this class skips deactivated cells and faces between activated and
+   * deactivated cells are treated as boundary faces.
+   */
+  template <int dim,
+            typename Number,
+            typename VectorizedArrayType = VectorizedArray<Number>>
+  class ElementActivationAndDeactivationMatrixFree
+  {
+  public:
+    /**
+     * Struct that helps to configure
+     * ElementActivationAndDeactivationMatrixFree.
+     */
+    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.
+     *
+     * @note At the moment, only DoFHandler objects are accepted
+     * that are initialized with FECollection objects with at most
+     * two finite elements (i.e., `FENothing` and another finite
+     * element).
+     */
+    void
+    reinit(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+           const AdditionalData &additional_data = AdditionalData())
+    {
+      this->matrix_free = &matrix_free;
+
+      std::vector<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.push_back(i);
+
+      // TODO: relax this so that arbitrary number of valid
+      // FEs can be accepted
+      AssertDimension(valid_fe_indices.size(), 1);
+
+      fe_index_valid = *valid_fe_indices.begin();
+    }
+
+    /**
+     * Loop over all activated 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 activated cells and faces. Faces between activated and
+     * deactivated 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; // deactivated 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
diff --git a/tests/matrix_free/element_birth_and_death_01.cc b/tests/matrix_free/element_birth_and_death_01.cc
new file mode 100644 (file)
index 0000000..53171d0
--- /dev/null
@@ -0,0 +1,209 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Test MatrixFreeTools::ElementActivationAndDeactivationMatrixFree by comparing
+// the results with results obtained on fitted mesh.
+
+using namespace dealii;
+
+template <int dim>
+class Fu : public Function<dim>
+{
+public:
+  double
+  value(const Point<dim> &point, const unsigned int = 0) const
+  {
+    return std::sin(point[0] * 2 * numbers::PI) *
+           std::sin(point[1] * 2 * numbers::PI);
+  }
+};
+
+template <int dim>
+void
+test(const unsigned int n_refinements)
+{
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+  using VectorType          = LinearAlgebra::distributed::Vector<Number>;
+  using FECellIntegrator =
+    FEEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType>;
+  using FEFaceIntegrator =
+    FEFaceEvaluation<dim, -1, 0, 1, Number, VectorizedArrayType>;
+
+  for (unsigned int i = 0; i < 2; ++i)
+    {
+      parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+      if (i == 0)
+        GridGenerator::subdivided_hyper_rectangle(
+          tria, {2, 2}, {0.0, 0.0}, {1.0, 1.0}, true);
+      else
+        GridGenerator::subdivided_hyper_rectangle(
+          tria, {2, 1}, {0.0, 0.0}, {1.0, 0.5}, true);
+
+      tria.refine_global(n_refinements);
+
+      DoFHandler<dim> dof_handler(tria);
+
+      for (const auto &cell :
+           filter_iterators(dof_handler.active_cell_iterators(),
+                            IteratorFilters::LocallyOwnedCell()))
+        {
+          if (cell->center()[1] < 0.5)
+            cell->set_active_fe_index(0);
+          else
+            cell->set_active_fe_index(1);
+        }
+
+      dof_handler.distribute_dofs(
+        hp::FECollection<dim>(FE_Q<dim>(1), FE_Nothing<dim>(1)));
+
+      typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+        data;
+      data.tasks_parallel_scheme =
+        MatrixFree<dim, double>::AdditionalData::none;
+      data.mapping_update_flags                = update_gradients;
+      data.mapping_update_flags_boundary_faces = update_gradients;
+      data.mapping_update_flags_inner_faces    = update_gradients;
+
+      MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+      matrix_free.reinit(MappingQ1<dim>(),
+                         dof_handler,
+                         AffineConstraints<Number>(),
+                         QGauss<dim>(2),
+                         data);
+
+      VectorType dst, src;
+      matrix_free.initialize_dof_vector(dst);
+      matrix_free.initialize_dof_vector(src);
+
+      VectorTools::interpolate(dof_handler, Fu<dim>(), src);
+
+      const auto cell_operation = [&](const auto &matrix_free,
+                                      auto &      dst,
+                                      const auto &src,
+                                      const auto  range) {
+        FECellIntegrator phi(matrix_free);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            phi.reinit(cell);
+            phi.gather_evaluate(src, EvaluationFlags::gradients);
+
+            for (const auto q : phi.quadrature_point_indices())
+              phi.submit_gradient(phi.get_gradient(q), q);
+
+            phi.integrate_scatter(EvaluationFlags::gradients, dst);
+          }
+      };
+
+      const auto face_operation =
+        [&](const auto &, auto &, const auto &, const auto) {
+          // nothing to do here but in the DG case
+        };
+
+      const auto boundary_operation = [&](const auto &matrix_free,
+                                          auto &,
+                                          const auto &,
+                                          const auto range,
+                                          const bool is_interior_face = true) {
+        FEFaceIntegrator phi(matrix_free, is_interior_face);
+
+        for (unsigned int face = range.first; face < range.second; ++face)
+          {
+            phi.reinit(face);
+
+            const double scaling =
+              phi.at_boundary() ? (phi.boundary_id() + 1.0) : (2.0 * dim);
+
+            phi.gather_evaluate(src, EvaluationFlags::gradients);
+
+            for (const auto q : phi.quadrature_point_indices())
+              {
+                auto gradient = phi.get_gradient(q);
+                auto normal   = phi.get_normal_vector(q);
+
+                if (is_interior_face == false) // fix sign!
+                  normal *= -1.0;
+
+                phi.submit_value(gradient * normal * scaling, q);
+              }
+
+            phi.integrate_scatter(EvaluationFlags::values, dst);
+          }
+      };
+
+      if (i == 0)
+        {
+          MatrixFreeTools::ElementActivationAndDeactivationMatrixFree<
+            dim,
+            Number,
+            VectorizedArrayType>
+            element_birth_and_death;
+
+          element_birth_and_death.reinit(matrix_free);
+
+          element_birth_and_death.template cell_loop<VectorType, VectorType>(
+            cell_operation, dst, src, true);
+
+          deallog << dst.l2_norm() << std::endl;
+
+          element_birth_and_death.template loop<VectorType, VectorType>(
+            cell_operation, face_operation, boundary_operation, dst, src, true);
+
+          deallog << dst.l2_norm() << std::endl;
+        }
+      else
+        {
+          matrix_free.template cell_loop<VectorType, VectorType>(cell_operation,
+                                                                 dst,
+                                                                 src,
+                                                                 true);
+
+          deallog << dst.l2_norm() << std::endl;
+
+          matrix_free.template loop<VectorType, VectorType>(
+            cell_operation, face_operation, boundary_operation, dst, src, true);
+
+          deallog << dst.l2_norm() << std::endl;
+        }
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  mpi_initlog();
+  test<2>(2);
+}
diff --git a/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..b83d53d
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::3.59694
+DEAL::9.14748
+DEAL::3.59694
+DEAL::9.14748
diff --git a/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output b/tests/matrix_free/element_birth_and_death_01.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..b83d53d
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::3.59694
+DEAL::9.14748
+DEAL::3.59694
+DEAL::9.14748

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.