From: Peter Munch <peterrmuench@gmail.com>
Date: Wed, 7 Aug 2024 08:00:45 +0000 (+0200)
Subject: Add DoFTools::extract_elasticity_modes()
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b40548f4036f9662865ecc12baabc7e7b07e4b94;p=dealii.git

Add DoFTools::extract_elasticity_modes()
---

diff --git a/doc/news/changes/minor/20240810Munch b/doc/news/changes/minor/20240810Munch
new file mode 100644
index 0000000000..767fd3f103
--- /dev/null
+++ b/doc/news/changes/minor/20240810Munch
@@ -0,0 +1,6 @@
+New: The new functions DoFTools::extract_elasticity_modes
+and DoFTools::extract_level_elasticity_modes
+allow you to extract translational and rotational modes,
+need to set up AMG for solving elasticity problems.
+<br>
+(Marco Feder, Peter Munch, Richard Schussnig, 2024/08/10)
diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h
index 39442f70b7..a0722db5e3 100644
--- a/include/deal.II/dofs/dof_tools.h
+++ b/include/deal.II/dofs/dof_tools.h
@@ -1382,6 +1382,30 @@ namespace DoFTools
                                const DoFHandler<dim, spacedim> &dof_handler,
                                const ComponentMask             &component_mask,
                                std::vector<std::vector<bool>>  &constant_modes);
+
+  /**
+   * Same as the above but additional to the translational modes also
+   * the rotational modes are added, needed to set up AMG for solving
+   * elasticity problems.
+   */
+  template <int dim, int spacedim>
+  void
+  extract_elasticity_modes(const Mapping<dim, spacedim>     &mapping,
+                           const DoFHandler<dim, spacedim>  &dof_handler,
+                           const ComponentMask              &component_mask,
+                           std::vector<std::vector<double>> &elasticity_modes);
+
+  /**
+   * Same as above but for multigrid levels.
+   */
+  template <int dim, int spacedim>
+  void
+  extract_level_elasticity_modes(
+    const unsigned int                level,
+    const Mapping<dim, spacedim>     &mapping,
+    const DoFHandler<dim, spacedim>  &dof_handler,
+    const ComponentMask              &component_mask,
+    std::vector<std::vector<double>> &elasticity_modes);
   /** @} */
 
   /**
diff --git a/include/deal.II/lac/trilinos_precondition.h b/include/deal.II/lac/trilinos_precondition.h
index 2f99daa1d9..7d3c88ef39 100644
--- a/include/deal.II/lac/trilinos_precondition.h
+++ b/include/deal.II/lac/trilinos_precondition.h
@@ -1499,6 +1499,15 @@ namespace TrilinosWrappers
        */
       std::vector<std::vector<bool>> constant_modes;
 
+      /**
+       * Same as above, but with values instead of Booleans. This
+       * is usfule if you want to specifiy rotational modes
+       * in addition to the translational modes. See also:
+       * DoFTools::extract_elasticity_modes
+       * and DoFTools::extract_level_elasticity_modes.
+       */
+      std::vector<std::vector<double>> constant_modes_values;
+
       /**
        * Determines how many sweeps of the smoother should be performed. When
        * the flag <tt>elliptic</tt> is set to <tt>true</tt>, i.e., for
diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index 43f51b7c9f..7bf4e2ea42 100644
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -43,6 +43,8 @@
 #include <deal.II/lac/sparsity_pattern_base.h>
 #include <deal.II/lac/vector.h>
 
+#include <deal.II/numerics/vector_tools.h>
+
 #include <algorithm>
 #include <numeric>
 
@@ -1407,6 +1409,130 @@ namespace DoFTools
             }
       });
     }
+
+
+
+    /**
+     * Definition of the rigid body motions for linear elasticity.
+     */
+    template <int dim>
+    class RigidBodyMotion : public Function<dim>
+    {
+    public:
+      static constexpr unsigned int n_modes = dim * (dim + 1) / 2;
+
+      RigidBodyMotion(const unsigned int type);
+
+      virtual double
+      value(const Point<dim> &p, const unsigned int component) const override;
+
+    private:
+      const unsigned int type;
+    };
+
+
+
+    template <int dim>
+    RigidBodyMotion<dim>::RigidBodyMotion(const unsigned int type)
+      : Function<dim>(dim)
+      , type(type)
+    {
+      Assert(type < n_modes, ExcNotImplemented());
+    }
+
+
+
+    Tensor<1, 2>
+    cross_product(const Tensor<1, 2> &tensor1, const Tensor<1, 1> &tensor2)
+    {
+      // |a|   |0|   |+bc|
+      // |b| x |0| = |-ac|
+      // |0|   |c|   | 0 |
+
+      Tensor<1, 2> cproduct;
+      cproduct[0] = +tensor1[1] * tensor2[0];
+      cproduct[1] = -tensor1[0] * tensor2[0];
+      return cproduct;
+    }
+
+
+
+    Tensor<1, 3>
+    cross_product(const Tensor<1, 3> &tensor1, const Tensor<1, 3> &tensor2)
+    {
+      Tensor<1, 3> cproduct;
+      cproduct[0] = +tensor1[1] * tensor2[2] - tensor1[2] * tensor2[1];
+      cproduct[1] = +tensor1[2] * tensor2[0] - tensor1[0] * tensor2[2];
+      cproduct[2] = +tensor1[0] * tensor2[1] - tensor1[1] * tensor2[0];
+      return cproduct;
+    }
+
+
+
+    template <int dim>
+    double
+    RigidBodyMotion<dim>::value(const Point<dim>  &p,
+                                const unsigned int component) const
+    {
+      if (type < dim) // translation modes
+        return static_cast<double>(component == type);
+
+      if constexpr (dim >= 2) // rotation modes
+        {
+          Tensor<1, n_modes - dim> dir;
+          dir[type - dim] = 1.0;
+
+          return cross_product(p, dir)[component];
+        }
+      else
+        {
+          Assert(false, ExcNotImplemented());
+
+          return 0.0;
+        }
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    extract_elasticity_modes(const Mapping<dim, spacedim>     &mapping,
+                             const DoFHandler<dim, spacedim>  &dof_handler,
+                             const ComponentMask              &component_mask,
+                             const unsigned int                mg_level,
+                             std::vector<std::vector<double>> &elasticity_modes)
+    {
+      AssertDimension(dim, spacedim);
+      AssertDimension(component_mask.n_selected_components(), dim);
+
+      const unsigned int n_modes = RigidBodyMotion<dim>::n_modes;
+
+      elasticity_modes.resize(n_modes);
+
+      LinearAlgebra::distributed::Vector<double> elasticity_modes_dealii(
+        mg_level == numbers::invalid_unsigned_int ?
+          dof_handler.locally_owned_dofs() :
+          dof_handler.locally_owned_mg_dofs(mg_level),
+        mg_level == numbers::invalid_unsigned_int ?
+          DoFTools::extract_locally_active_dofs(dof_handler) :
+          DoFTools::extract_locally_active_level_dofs(dof_handler, mg_level),
+        dof_handler.get_communicator());
+
+      for (unsigned int i = 0; i < n_modes; ++i)
+        {
+          VectorTools::interpolate(mapping,
+                                   dof_handler,
+                                   RigidBodyMotion<dim>(i),
+                                   elasticity_modes_dealii,
+                                   component_mask,
+                                   mg_level);
+
+          // copy to right format
+          elasticity_modes[i].assign(elasticity_modes_dealii.begin(),
+                                     elasticity_modes_dealii.end());
+        }
+    }
+
   } // namespace internal
 
 
@@ -1440,6 +1566,37 @@ namespace DoFTools
 
 
 
+  template <int dim, int spacedim>
+  void
+  extract_elasticity_modes(const Mapping<dim, spacedim>     &mapping,
+                           const DoFHandler<dim, spacedim>  &dof_handler,
+                           const ComponentMask              &component_mask,
+                           std::vector<std::vector<double>> &constant_modes)
+  {
+    internal::extract_elasticity_modes(mapping,
+                                       dof_handler,
+                                       component_mask,
+                                       numbers::invalid_unsigned_int,
+                                       constant_modes);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  extract_level_elasticity_modes(
+    const unsigned int                level,
+    const Mapping<dim, spacedim>     &mapping,
+    const DoFHandler<dim, spacedim>  &dof_handler,
+    const ComponentMask              &component_mask,
+    std::vector<std::vector<double>> &constant_modes)
+  {
+    internal::extract_elasticity_modes(
+      mapping, dof_handler, component_mask, level, constant_modes);
+  }
+
+
+
   template <int dim, int spacedim>
   std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
            std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in
index b891f1f3a1..30732c63dd 100644
--- a/source/dofs/dof_tools.inst.in
+++ b/source/dofs/dof_tools.inst.in
@@ -199,6 +199,19 @@ for (deal_II_dimension : DIMENSIONS)
       const ComponentMask                 &selected_components,
       std::vector<std::vector<bool>>      &constant_modes);
 
+    template void DoFTools::extract_elasticity_modes<deal_II_dimension>(
+      const Mapping<deal_II_dimension>    &mapping,
+      const DoFHandler<deal_II_dimension> &dof_handler,
+      const ComponentMask                 &selected_components,
+      std::vector<std::vector<double>>    &constant_modes);
+
+    template void DoFTools::extract_level_elasticity_modes<deal_II_dimension>(
+      const unsigned int                   level,
+      const Mapping<deal_II_dimension>    &mapping,
+      const DoFHandler<deal_II_dimension> &dof_handler,
+      const ComponentMask                 &selected_components,
+      std::vector<std::vector<double>>    &constant_modes);
+
     template void DoFTools::get_active_fe_indices<deal_II_dimension>(
       const DoFHandler<deal_II_dimension> &dof_handler,
       std::vector<unsigned int>           &active_fe_indices);
diff --git a/source/lac/trilinos_precondition_ml.cc b/source/lac/trilinos_precondition_ml.cc
index 903b3392d5..af3c801540 100644
--- a/source/lac/trilinos_precondition_ml.cc
+++ b/source/lac/trilinos_precondition_ml.cc
@@ -128,57 +128,70 @@ namespace TrilinosWrappers
     std::unique_ptr<Epetra_MultiVector> &ptr_distributed_constant_modes,
     const Epetra_RowMatrix              &matrix) const
   {
-    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
-
-    const size_type constant_modes_dimension = constant_modes.size();
-    ptr_distributed_constant_modes = std::make_unique<Epetra_MultiVector>(
-      domain_map, constant_modes_dimension > 0 ? constant_modes_dimension : 1);
-    Assert(ptr_distributed_constant_modes, ExcNotInitialized());
-    Epetra_MultiVector &distributed_constant_modes =
-      *ptr_distributed_constant_modes;
-
-    if (constant_modes_dimension > 0)
+    const auto run = [&](const auto &constant_modes) {
+      const Epetra_Map &domain_map = matrix.OperatorDomainMap();
+
+      const size_type constant_modes_dimension = constant_modes.size();
+      ptr_distributed_constant_modes =
+        std::make_unique<Epetra_MultiVector>(domain_map,
+                                             constant_modes_dimension > 0 ?
+                                               constant_modes_dimension :
+                                               1);
+      Assert(ptr_distributed_constant_modes, ExcNotInitialized());
+      Epetra_MultiVector &distributed_constant_modes =
+        *ptr_distributed_constant_modes;
+
+      if (constant_modes_dimension > 0)
+        {
+          const size_type global_size = TrilinosWrappers::n_global_rows(matrix);
+          (void)global_length; // work around compiler warning about unused
+                               // function in release mode
+          Assert(global_size ==
+                   static_cast<size_type>(TrilinosWrappers::global_length(
+                     distributed_constant_modes)),
+                 ExcDimensionMismatch(global_size,
+                                      TrilinosWrappers::global_length(
+                                        distributed_constant_modes)));
+          const bool constant_modes_are_global =
+            constant_modes[0].size() == global_size;
+          const size_type my_size = domain_map.NumMyElements();
+
+          // Reshape null space as a contiguous vector of doubles so that
+          // Trilinos can read from it.
+          const size_type expected_mode_size =
+            constant_modes_are_global ? global_size : my_size;
+          for (size_type d = 0; d < constant_modes_dimension; ++d)
+            {
+              Assert(constant_modes[d].size() == expected_mode_size,
+                     ExcDimensionMismatch(constant_modes[d].size(),
+                                          expected_mode_size));
+              for (size_type row = 0; row < my_size; ++row)
+                {
+                  const TrilinosWrappers::types::int_type mode_index =
+                    constant_modes_are_global ?
+                      TrilinosWrappers::global_index(domain_map, row) :
+                      row;
+                  distributed_constant_modes[d][row] =
+                    static_cast<double>(constant_modes[d][mode_index]);
+                }
+            }
+          (void)expected_mode_size;
+
+          parameter_list.set("null space: type", "pre-computed");
+          parameter_list.set("null space: dimension",
+                             distributed_constant_modes.NumVectors());
+          parameter_list.set("null space: vectors",
+                             distributed_constant_modes.Values());
+        }
+    };
+
+    if (!constant_modes_values.empty())
       {
-        const size_type global_size = TrilinosWrappers::n_global_rows(matrix);
-        (void)global_length; // work around compiler warning about unused
-                             // function in release mode
-        Assert(global_size ==
-                 static_cast<size_type>(
-                   TrilinosWrappers::global_length(distributed_constant_modes)),
-               ExcDimensionMismatch(global_size,
-                                    TrilinosWrappers::global_length(
-                                      distributed_constant_modes)));
-        const bool constant_modes_are_global =
-          constant_modes[0].size() == global_size;
-        const size_type my_size = domain_map.NumMyElements();
-
-        // Reshape null space as a contiguous vector of doubles so that
-        // Trilinos can read from it.
-        const size_type expected_mode_size =
-          constant_modes_are_global ? global_size : my_size;
-        for (size_type d = 0; d < constant_modes_dimension; ++d)
-          {
-            Assert(constant_modes[d].size() == expected_mode_size,
-                   ExcDimensionMismatch(constant_modes[d].size(),
-                                        expected_mode_size));
-            for (size_type row = 0; row < my_size; ++row)
-              {
-                const TrilinosWrappers::types::int_type mode_index =
-                  constant_modes_are_global ?
-                    TrilinosWrappers::global_index(domain_map, row) :
-                    row;
-                distributed_constant_modes[d][row] =
-                  static_cast<double>(constant_modes[d][mode_index]);
-              }
-          }
-        (void)expected_mode_size;
-
-        parameter_list.set("null space: type", "pre-computed");
-        parameter_list.set("null space: dimension",
-                           distributed_constant_modes.NumVectors());
-        parameter_list.set("null space: vectors",
-                           distributed_constant_modes.Values());
+        AssertDimension(constant_modes.size(), 0);
+        run(constant_modes_values);
       }
+    else
+      run(constant_modes);
   }
 
 
diff --git a/tests/mpi/extract_constant_modes_05.cc b/tests/mpi/extract_constant_modes_05.cc
new file mode 100644
index 0000000000..70b738e406
--- /dev/null
+++ b/tests/mpi/extract_constant_modes_05.cc
@@ -0,0 +1,87 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// test DoFTools::extract_elasticity_modes
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+
+  MappingQ1<dim> mapping;
+
+  FESystem<dim>   fe(FE_Q<dim>(1), dim);
+  DoFHandler<dim> dofh(tr);
+  dofh.distribute_dofs(fe);
+
+  deallog << "Total dofs=" << dofh.n_dofs() << std::endl;
+
+  // extract constant modes and print
+  ComponentMask mask(fe.n_components(), true);
+
+  std::vector<std::vector<double>> constant_modes;
+  DoFTools::extract_elasticity_modes(mapping, dofh, mask, constant_modes);
+
+  for (unsigned int i = 0; i < constant_modes.size(); ++i)
+    {
+      for (unsigned int j = 0; j < constant_modes[i].size(); ++j)
+        deallog << constant_modes[i][j] << ' ';
+      deallog << std::endl;
+    }
+}
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+  deallog.push("2d");
+  test<2>();
+  deallog.pop();
+
+  deallog.push("3d");
+  test<3>();
+  deallog.pop();
+}
diff --git a/tests/mpi/extract_constant_modes_05.mpirun=4.with_p4est=true.output b/tests/mpi/extract_constant_modes_05.mpirun=4.with_p4est=true.output
new file mode 100644
index 0000000000..7d870e34b6
--- /dev/null
+++ b/tests/mpi/extract_constant_modes_05.mpirun=4.with_p4est=true.output
@@ -0,0 +1,51 @@
+
+DEAL:0:2d::Total dofs=18
+DEAL:0:2d::1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 
+DEAL:0:2d::0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 
+DEAL:0:2d::0.00000 0.00000 0.00000 -0.500000 0.500000 0.00000 0.500000 -0.500000 
+DEAL:0:3d::Total dofs=81
+DEAL:0:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:0:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:0:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.500000 -0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -0.500000 0.00000 1.00000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 
+
+DEAL:1:2d::Total dofs=18
+DEAL:1:2d::1.00000 0.00000 1.00000 0.00000 
+DEAL:1:2d::0.00000 1.00000 0.00000 1.00000 
+DEAL:1:2d::0.00000 -1.00000 0.500000 -1.00000 
+DEAL:1:3d::Total dofs=81
+DEAL:1:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:1:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:1:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:1:3d::0.00000 0.00000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 
+DEAL:1:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 
+DEAL:1:3d::1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 0.00000 
+
+
+DEAL:2:2d::Total dofs=18
+DEAL:2:2d::1.00000 0.00000 1.00000 0.00000 
+DEAL:2:2d::0.00000 1.00000 0.00000 1.00000 
+DEAL:2:2d::1.00000 0.00000 1.00000 -0.500000 
+DEAL:2:3d::Total dofs=81
+DEAL:2:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:2:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:2:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:2:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -0.500000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 
+DEAL:2:3d::-1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 1.00000 -1.00000 0.00000 1.00000 
+DEAL:2:3d::0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 
+
+
+DEAL:3:2d::Total dofs=18
+DEAL:3:2d::1.00000 0.00000 
+DEAL:3:2d::0.00000 1.00000 
+DEAL:3:2d::1.00000 -1.00000 
+DEAL:3:3d::Total dofs=81
+DEAL:3:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:3:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:3:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:3:3d::0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 
+DEAL:3:3d::-1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 1.00000 
+DEAL:3:3d::1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -1.00000 0.00000 
+
diff --git a/tests/mpi/extract_constant_modes_06.cc b/tests/mpi/extract_constant_modes_06.cc
new file mode 100644
index 0000000000..660eed0e95
--- /dev/null
+++ b/tests/mpi/extract_constant_modes_06.cc
@@ -0,0 +1,92 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// test DoFTools::extract_level_elasticity_modes
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  parallel::distributed::Triangulation<dim> tr(
+    MPI_COMM_WORLD,
+    Triangulation<dim>::none,
+    parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::subdivided_hyper_cube(tr, 2);
+
+  MappingQ1<dim> mapping;
+
+  FESystem<dim>   fe(FE_Q<dim>(1), dim);
+  DoFHandler<dim> dofh(tr);
+  dofh.distribute_dofs(fe);
+  dofh.distribute_mg_dofs();
+
+  deallog << "Total dofs=" << dofh.n_dofs() << std::endl;
+
+  // extract constant modes and print
+  ComponentMask mask(fe.n_components(), true);
+
+  std::vector<std::vector<double>> constant_modes;
+  DoFTools::extract_level_elasticity_modes(
+    0, mapping, dofh, mask, constant_modes);
+
+  for (unsigned int i = 0; i < constant_modes.size(); ++i)
+    {
+      for (unsigned int j = 0; j < constant_modes[i].size(); ++j)
+        deallog << constant_modes[i][j] << ' ';
+      deallog << std::endl;
+    }
+}
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+  deallog.push("2d");
+  test<2>();
+  deallog.pop();
+
+  deallog.push("3d");
+  test<3>();
+  deallog.pop();
+}
diff --git a/tests/mpi/extract_constant_modes_06.mpirun=4.with_p4est=true.output b/tests/mpi/extract_constant_modes_06.mpirun=4.with_p4est=true.output
new file mode 100644
index 0000000000..7d870e34b6
--- /dev/null
+++ b/tests/mpi/extract_constant_modes_06.mpirun=4.with_p4est=true.output
@@ -0,0 +1,51 @@
+
+DEAL:0:2d::Total dofs=18
+DEAL:0:2d::1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 
+DEAL:0:2d::0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 0.00000 1.00000 
+DEAL:0:2d::0.00000 0.00000 0.00000 -0.500000 0.500000 0.00000 0.500000 -0.500000 
+DEAL:0:3d::Total dofs=81
+DEAL:0:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:0:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:0:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.500000 -0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -0.500000 0.00000 1.00000 
+DEAL:0:3d::0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 
+
+DEAL:1:2d::Total dofs=18
+DEAL:1:2d::1.00000 0.00000 1.00000 0.00000 
+DEAL:1:2d::0.00000 1.00000 0.00000 1.00000 
+DEAL:1:2d::0.00000 -1.00000 0.500000 -1.00000 
+DEAL:1:3d::Total dofs=81
+DEAL:1:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:1:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:1:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:1:3d::0.00000 0.00000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 
+DEAL:1:3d::0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 
+DEAL:1:3d::1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 0.00000 
+
+
+DEAL:2:2d::Total dofs=18
+DEAL:2:2d::1.00000 0.00000 1.00000 0.00000 
+DEAL:2:2d::0.00000 1.00000 0.00000 1.00000 
+DEAL:2:2d::1.00000 0.00000 1.00000 -0.500000 
+DEAL:2:3d::Total dofs=81
+DEAL:2:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:2:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:2:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:2:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -0.500000 0.00000 1.00000 0.00000 0.00000 1.00000 -0.500000 
+DEAL:2:3d::-1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 1.00000 -1.00000 0.00000 1.00000 
+DEAL:2:3d::0.00000 0.00000 0.00000 0.00000 -0.500000 0.00000 0.500000 0.00000 0.00000 0.500000 -0.500000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 
+
+
+DEAL:3:2d::Total dofs=18
+DEAL:3:2d::1.00000 0.00000 
+DEAL:3:2d::0.00000 1.00000 
+DEAL:3:2d::1.00000 -1.00000 
+DEAL:3:3d::Total dofs=81
+DEAL:3:3d::1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
+DEAL:3:3d::0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:3:3d::0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:3:3d::0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 0.00000 1.00000 -1.00000 
+DEAL:3:3d::-1.00000 0.00000 0.00000 -1.00000 0.00000 0.500000 -1.00000 0.00000 1.00000 
+DEAL:3:3d::1.00000 0.00000 0.00000 1.00000 -0.500000 0.00000 1.00000 -1.00000 0.00000 
+