From 82327c82049b96034c3e3e91b5599c1bb6c22552 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 3 Aug 2024 11:55:33 +0200
Subject: [PATCH] Fix MatrixFreeTools::compute_diagonal

---
 doc/news/changes/minor/20240803Munch         |  5 +++
 include/deal.II/matrix_free/tools.h          | 38 +++++++++++++-------
 tests/matrix_free/compute_diagonal_07.cc     | 22 +++++++-----
 tests/matrix_free/compute_diagonal_07.output | 16 +++++++++
 4 files changed, 61 insertions(+), 20 deletions(-)
 create mode 100644 doc/news/changes/minor/20240803Munch

diff --git a/doc/news/changes/minor/20240803Munch b/doc/news/changes/minor/20240803Munch
new file mode 100644
index 0000000000..106f83220e
--- /dev/null
+++ b/doc/news/changes/minor/20240803Munch
@@ -0,0 +1,5 @@
+Fixed: MatrixFreeTools::compute_diagonal() used
+to use first_selected_component also for accessing
+(block) vectors. This has been fixed.
+<br>
+(Peter Munch, Nils Much, 2024/07/30)
diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index 6e3928a1aa..8b9c1f9f64 100644
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -76,6 +76,9 @@ namespace MatrixFreeTools
    *
    * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
    * passed to the constructor of the FEEvaluation that is internally set up.
+   *
+   * The parameter @p first_vector_component is used to select the right
+   * starting block in a block vector.
    */
   template <int dim,
             int fe_degree,
@@ -97,7 +100,8 @@ namespace MatrixFreeTools
                       &cell_operation,
     const unsigned int dof_no                   = 0,
     const unsigned int quad_no                  = 0,
-    const unsigned int first_selected_component = 0);
+    const unsigned int first_selected_component = 0,
+    const unsigned int first_vector_component   = 0);
 
 
 
@@ -125,7 +129,8 @@ namespace MatrixFreeTools
     const CLASS       *owning_class,
     const unsigned int dof_no                   = 0,
     const unsigned int quad_no                  = 0,
-    const unsigned int first_selected_component = 0);
+    const unsigned int first_selected_component = 0,
+    const unsigned int first_vector_component   = 0);
 
 
 
@@ -179,7 +184,8 @@ namespace MatrixFreeTools
                       &boundary_operation,
     const unsigned int dof_no                   = 0,
     const unsigned int quad_no                  = 0,
-    const unsigned int first_selected_component = 0);
+    const unsigned int first_selected_component = 0,
+    const unsigned int first_vector_component   = 0);
 
 
 
@@ -227,7 +233,8 @@ namespace MatrixFreeTools
     const CLASS       *owning_class,
     const unsigned int dof_no                   = 0,
     const unsigned int quad_no                  = 0,
-    const unsigned int first_selected_component = 0);
+    const unsigned int first_selected_component = 0,
+    const unsigned int first_vector_component   = 0);
 
 
 
@@ -1352,7 +1359,8 @@ namespace MatrixFreeTools
                       &cell_operation,
     const unsigned int dof_no,
     const unsigned int quad_no,
-    const unsigned int first_selected_component)
+    const unsigned int first_selected_component,
+    const unsigned int first_vector_component)
   {
     compute_diagonal<dim,
                      fe_degree,
@@ -1367,7 +1375,8 @@ namespace MatrixFreeTools
                                  {},
                                  dof_no,
                                  quad_no,
-                                 first_selected_component);
+                                 first_selected_component,
+                                 first_vector_component);
   }
 
   template <typename CLASS,
@@ -1391,7 +1400,8 @@ namespace MatrixFreeTools
     const CLASS       *owning_class,
     const unsigned int dof_no,
     const unsigned int quad_no,
-    const unsigned int first_selected_component)
+    const unsigned int first_selected_component,
+    const unsigned int first_vector_component)
   {
     compute_diagonal<dim,
                      fe_degree,
@@ -1405,7 +1415,8 @@ namespace MatrixFreeTools
       [&](auto &phi) { (owning_class->*cell_operation)(phi); },
       dof_no,
       quad_no,
-      first_selected_component);
+      first_selected_component,
+      first_vector_component);
   }
 
   template <int dim,
@@ -1448,7 +1459,8 @@ namespace MatrixFreeTools
                       &boundary_operation,
     const unsigned int dof_no,
     const unsigned int quad_no,
-    const unsigned int first_selected_component)
+    const unsigned int first_selected_component,
+    const unsigned int first_vector_component)
   {
     int dummy = 0;
 
@@ -1461,7 +1473,7 @@ namespace MatrixFreeTools
     for (unsigned int d = 0; d < n_components; ++d)
       diagonal_global_components[d] = dealii::internal::
         BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::
-          get_vector_component(diagonal_global, d + first_selected_component);
+          get_vector_component(diagonal_global, d + first_vector_component);
 
     const auto &dof_info = matrix_free.get_dof_info(dof_no);
 
@@ -1745,7 +1757,8 @@ namespace MatrixFreeTools
     const CLASS       *owning_class,
     const unsigned int dof_no,
     const unsigned int quad_no,
-    const unsigned int first_selected_component)
+    const unsigned int first_selected_component,
+    const unsigned int first_vector_component)
   {
     compute_diagonal<dim,
                      fe_degree,
@@ -1763,7 +1776,8 @@ namespace MatrixFreeTools
       [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
       dof_no,
       quad_no,
-      first_selected_component);
+      first_selected_component,
+      first_vector_component);
   }
 
   namespace internal
diff --git a/tests/matrix_free/compute_diagonal_07.cc b/tests/matrix_free/compute_diagonal_07.cc
index d134e59d9d..94d72e5004 100644
--- a/tests/matrix_free/compute_diagonal_07.cc
+++ b/tests/matrix_free/compute_diagonal_07.cc
@@ -29,7 +29,7 @@ template <int dim,
           typename Number              = double,
           typename VectorizedArrayType = VectorizedArray<Number>>
 void
-test()
+test(const unsigned int first_vector_component)
 {
   Triangulation<dim> tria;
   GridGenerator::hyper_ball(tria);
@@ -104,9 +104,10 @@ test()
 
   diagonal_global.print(deallog.get_file_stream());
 
-  LinearAlgebra::distributed::BlockVector<Number> diagonal_global_block(dim);
+  LinearAlgebra::distributed::BlockVector<Number> diagonal_global_block(
+    dim + first_vector_component);
 
-  for (unsigned int comp = 0; comp < dim; ++comp)
+  for (unsigned int comp = 0; comp < dim + first_vector_component; ++comp)
     matrix_free.initialize_dof_vector(diagonal_global_block.block(comp), 1);
 
   MatrixFreeTools::compute_diagonal<dim,
@@ -114,10 +115,14 @@ test()
                                     n_points,
                                     n_components,
                                     Number,
-                                    VectorizedArrayType>(matrix_free,
-                                                         diagonal_global_block,
-                                                         kernel,
-                                                         1);
+                                    VectorizedArrayType>(
+    matrix_free,
+    diagonal_global_block,
+    kernel,
+    1,
+    0,
+    0,
+    first_vector_component);
 
   diagonal_global_block.print(deallog.get_file_stream());
 }
@@ -128,5 +133,6 @@ main(int argc, char **argv)
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
   MPILogInitAll                    all;
 
-  test<2, 1>();
+  test<2, 1>(0);
+  test<2, 1>(1 /*test first vector component*/);
 }
diff --git a/tests/matrix_free/compute_diagonal_07.output b/tests/matrix_free/compute_diagonal_07.output
index ecf8996667..fffb6018bf 100644
--- a/tests/matrix_free/compute_diagonal_07.output
+++ b/tests/matrix_free/compute_diagonal_07.output
@@ -11,3 +11,19 @@ Process #0
 Local range: [0, 8), global size: 8
 Vector data:
 0.000e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 
+Process #0
+Local range: [0, 16), global size: 16
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.621e+00 4.621e+00 4.621e+00 4.621e+00 0.000e+00 0.000e+00 4.621e+00 4.621e+00 4.621e+00 4.621e+00 0.000e+00 0.000e+00 
+Process #0
+Local range: [0, 8), global size: 8
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+Process #0
+Local range: [0, 8), global size: 8
+Vector data:
+0.000e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 
+Process #0
+Local range: [0, 8), global size: 8
+Vector data:
+0.000e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 
-- 
2.39.5