Add DoFTools::extract_elasticity_modes()

diff --git a/doc/news/changes/minor/20240810Munch b/doc/news/changes/minor/20240810Munch
new file mode 100644 (file)
index 0000000..767fd3f
--- /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 39442f70b76c2ef1ac96a936b24953ad7b6a26e6..a0722db5e37effe05913d8252a6ac2558717af4e 100644 (file)
--- 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 2f99daa1d9f4af6f50679b5424c9a994af1c95c5..7d3c88ef392270303dc294c2a4a311e4f8fda474 100644 (file)
--- 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 43f51b7c9f7f7e93b5b3bd51b7776b3aa3e73099..7bf4e2ea424876106e7aa631d01f01348ac80045 100644 (file)
--- 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 b891f1f3a1a0de56bc4d4d4e7fc6a5f19338c04e..30732c63dd3121879cc27aab823498acb8c351ce 100644 (file)
--- 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 903b3392d5793f9fa1d9110ed67dcb7cce938911..af3c801540c3cbb5ff3940402d792575cd234faa 100644 (file)
--- 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 (file)
index 0000000..70b738e
--- /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 (file)
index 0000000..7d870e3
--- /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 (file)
index 0000000..660eed0
--- /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 (file)
index 0000000..7d870e3
--- /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 
+