From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 4 Aug 2024 11:08:29 +0000 (+0200)
Subject: Fix MatrixFreeTools::compute_matrix()
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d96ac99db76b2f5a1df9f9bd60ff392434d32006;p=dealii.git

Fix MatrixFreeTools::compute_matrix()
---

diff --git a/doc/news/changes/minor/20240804Munch b/doc/news/changes/minor/20240804Munch
new file mode 100644
index 0000000000..09bfbe1cc0
--- /dev/null
+++ b/doc/news/changes/minor/20240804Munch
@@ -0,0 +1,5 @@
+Fixed: MatrixFreeTools::compute_matrix() now
+correctly handles the case that individual
+components are requested.
+<br>
+(Peter Munch, 2024/08/04)
diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index 8b9c1f9f64..e23acadee2 100644
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -42,6 +42,7 @@ namespace MatrixFreeTools
 
       std::vector<unsigned int> dof_numbers;
       std::vector<unsigned int> quad_numbers;
+      std::vector<unsigned int> n_components;
       std::vector<unsigned int> first_selected_components;
       std::vector<unsigned int> batch_type;
       static const bool         is_face = is_face_;
@@ -1971,6 +1972,7 @@ namespace MatrixFreeTools
 
     data_cell.dof_numbers               = {dof_no};
     data_cell.quad_numbers              = {quad_no};
+    data_cell.n_components              = {n_components};
     data_cell.first_selected_components = {first_selected_component};
     data_cell.batch_type                = {0};
 
@@ -2007,6 +2009,7 @@ namespace MatrixFreeTools
 
     data_face.dof_numbers               = {dof_no, dof_no};
     data_face.quad_numbers              = {quad_no, quad_no};
+    data_face.n_components              = {n_components, n_components};
     data_face.first_selected_components = {first_selected_component,
                                            first_selected_component};
     data_face.batch_type                = {1, 2};
@@ -2069,6 +2072,7 @@ namespace MatrixFreeTools
 
     data_boundary.dof_numbers               = {dof_no};
     data_boundary.quad_numbers              = {quad_no};
+    data_boundary.n_components              = {n_components};
     data_boundary.first_selected_components = {first_selected_component};
     data_boundary.batch_type                = {1};
 
@@ -2155,22 +2159,39 @@ namespace MatrixFreeTools
 
           for (unsigned int b = 0; b < n_blocks; ++b)
             {
+              const auto &fe = matrix_free.get_dof_handler(data.dof_numbers[b])
+                                 .get_fe(phi[b]->get_active_fe_index());
+
+              const auto component_base =
+                matrix_free.get_dof_info(data.dof_numbers[b])
+                  .component_to_base_index[data.first_selected_components[b]];
+              const auto component_in_base =
+                data.first_selected_components[b] -
+                matrix_free.get_dof_info(data.dof_numbers[b])
+                  .start_components[component_base];
+
               const auto &shape_info = matrix_free.get_shape_info(
                 data.dof_numbers[b],
                 data.quad_numbers[b],
-                data.first_selected_components[b],
+                component_base,
                 phi[b]->get_active_fe_index(),
                 phi[b]->get_active_quadrature_index());
 
               dofs_per_cell[b] =
-                shape_info.dofs_per_component_on_cell * shape_info.n_components;
+                shape_info.dofs_per_component_on_cell * data.n_components[b];
 
-              dof_indices[b].resize(dofs_per_cell[b]);
+              dof_indices[b].resize(fe.n_dofs_per_cell());
 
               for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
                 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
 
-              lexicographic_numbering[b] = shape_info.lexicographic_numbering;
+              lexicographic_numbering[b].insert(
+                lexicographic_numbering[b].begin(),
+                shape_info.lexicographic_numbering.begin() +
+                  component_in_base * shape_info.dofs_per_component_on_cell,
+                shape_info.lexicographic_numbering.begin() +
+                  (component_in_base + data.n_components[b]) *
+                    shape_info.dofs_per_component_on_cell);
             }
 
           for (unsigned int bj = 0; bj < n_blocks; ++bj)
@@ -2210,7 +2231,7 @@ namespace MatrixFreeTools
                     else
                       cell_iterator->get_dof_indices(dof_indices[b]);
 
-                    for (unsigned int j = 0; j < dof_indices[b].size(); ++j)
+                    for (unsigned int j = 0; j < dofs_per_cell[b]; ++j)
                       dof_indices_mf[b][v][j] =
                         dof_indices[b][lexicographic_numbering[b][j]];
                   }
diff --git a/tests/matrix_free/compute_diagonal_10.cc b/tests/matrix_free/compute_diagonal_10.cc
new file mode 100644
index 0000000000..2e39129ec6
--- /dev/null
+++ b/tests/matrix_free/compute_diagonal_10.cc
@@ -0,0 +1,214 @@
+// ------------------------------------------------------------------------
+//
+// 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 MatrixFreeTools::compute_matrix():
+//   (1) compute individual components
+//   (2) compute with multiple FEEvaluation instances
+//
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.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/tria.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <unsigned int n_components,
+          int          dim,
+          typename Number,
+          typename VectorizedArrayType>
+void
+run(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const unsigned int first_selected_component)
+{
+  const unsigned int n_dofs = matrix_free.get_dof_handler().n_dofs();
+  FullMatrix<double> matrix(n_dofs, n_dofs);
+
+  FullMatrix<double> scaling(3, 3);
+  for (unsigned int i = 0, c = 1; i < 3; ++i)
+    for (unsigned int j = 0; j < 3; ++j, ++c)
+      scaling[i][j] = c;
+
+  MatrixFreeTools::
+    compute_matrix<dim, -1, 0, n_components, double, VectorizedArray<double>>(
+      matrix_free,
+      matrix_free.get_affine_constraints(),
+      matrix,
+      [&](auto &phi) {
+        phi.evaluate(EvaluationFlags::values);
+        for (const auto q : phi.quadrature_point_indices())
+          {
+            auto value = phi.get_value(q);
+
+            auto result = value;
+            result      = 0.0;
+
+            for (unsigned int i = 0; i < n_components; ++i)
+              for (unsigned int j = 0; j < n_components; ++j)
+                if constexpr (n_components > 1)
+                  result[i] += scaling[i + first_selected_component]
+                                      [j + first_selected_component] *
+                               value[j];
+                else
+                  result += scaling[i + first_selected_component]
+                                   [j + first_selected_component] *
+                            value;
+
+            phi.submit_value(result, q);
+          }
+        phi.integrate(EvaluationFlags::values);
+      },
+      0,
+      0,
+      first_selected_component);
+
+  matrix.print_formatted(deallog.get_file_stream(), 5, false, 10, "0.00000");
+  deallog << std::endl;
+}
+
+
+
+template <int dim>
+void
+test_0()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FESystem<dim>(FE_Q<dim>(2), 2, FE_Q<dim>(1), 1));
+  DoFRenumbering::component_wise(dof_handler);
+
+  QGauss<dim> quadrature(2);
+
+  MappingQ1<dim> mapping;
+
+  AffineConstraints<double> constraints;
+
+  MatrixFree<dim, double> matrix_free;
+
+  matrix_free.reinit(mapping, dof_handler, constraints, quadrature);
+
+  run<1>(matrix_free, 0u);
+  run<1>(matrix_free, 1u);
+  run<1>(matrix_free, 2u);
+
+  run<2>(matrix_free, 0u);
+}
+
+
+
+template <int dim>
+void
+test_1()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FESystem<dim>(FE_Q<dim>(2), 2));
+  DoFRenumbering::component_wise(dof_handler);
+
+  QGauss<dim> quadrature(2);
+
+  MappingQ1<dim> mapping;
+
+  AffineConstraints<double> constraints;
+
+  MatrixFree<dim, double> matrix_free;
+
+  matrix_free.reinit(mapping, dof_handler, constraints, quadrature);
+
+  MatrixFreeTools::internal::
+    ComputeMatrixScratchData<dim, VectorizedArray<double>, false>
+      data_cell;
+
+  data_cell.dof_numbers               = {0, 0};
+  data_cell.quad_numbers              = {0, 0};
+  data_cell.n_components              = {1, 1};
+  data_cell.first_selected_components = {0, 1};
+  data_cell.batch_type                = {0, 0};
+
+  data_cell.op_create =
+    [&](const std::pair<unsigned int, unsigned int> &range) {
+      std::vector<
+        std::unique_ptr<FEEvaluationData<dim, VectorizedArray<double>, false>>>
+        phi;
+
+      phi.emplace_back(std::make_unique<FEEvaluation<dim, -1, 0, 1, double>>(
+        matrix_free, range, 0, 0, 0));
+
+      phi.emplace_back(std::make_unique<FEEvaluation<dim, -1, 0, 1, double>>(
+        matrix_free, range, 0, 0, 1));
+
+      return phi;
+    };
+
+  data_cell.op_reinit = [](auto &phi, const unsigned batch) {
+    static_cast<FEEvaluation<dim, -1, 0, 1, double> &>(*phi[0]).reinit(batch);
+    static_cast<FEEvaluation<dim, -1, 0, 1, double> &>(*phi[1]).reinit(batch);
+  };
+
+  data_cell.op_compute = [&](auto &phi) {
+    auto &phi_0 = static_cast<FEEvaluation<dim, -1, 0, 1, double> &>(*phi[0]);
+    auto &phi_1 = static_cast<FEEvaluation<dim, -1, 0, 1, double> &>(*phi[1]);
+
+    phi_0.evaluate(EvaluationFlags::values);
+    phi_1.evaluate(EvaluationFlags::values);
+    for (const auto q : phi_0.quadrature_point_indices())
+      {
+        auto value_0 = phi_0.get_value(q);
+        auto value_1 = phi_1.get_value(q);
+
+        phi_0.submit_value(value_0 * 1.0 + value_1 * 2.0, q);
+        phi_1.submit_value(value_0 * 4.0 + value_1 * 5.0, q);
+      }
+    phi_0.integrate(EvaluationFlags::values);
+    phi_1.integrate(EvaluationFlags::values);
+  };
+
+  const unsigned int n_dofs = dof_handler.n_dofs();
+  FullMatrix<double> matrix(n_dofs, n_dofs);
+
+  MatrixFreeTools::internal::compute_matrix(
+    matrix_free, constraints, data_cell, {}, {}, matrix);
+
+  matrix.print_formatted(deallog.get_file_stream(), 5, false, 10, "0.00000");
+  deallog << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  test_0<1>(); // individual components
+  test_1<1>(); // multiple FEEvaluation instances
+}
diff --git a/tests/matrix_free/compute_diagonal_10.output b/tests/matrix_free/compute_diagonal_10.output
new file mode 100644
index 0000000000..5add115c38
--- /dev/null
+++ b/tests/matrix_free/compute_diagonal_10.output
@@ -0,0 +1,44 @@
+
+0.11111    -0.05556   0.11111    0.00000    0.00000    0.00000    0.00000    0.00000    
+-0.05556   0.11111    0.11111    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.11111    0.11111    0.44444    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+DEAL::
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.55556    -0.27778   0.55556    0.00000    0.00000    
+0.00000    0.00000    0.00000    -0.27778   0.55556    0.55556    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.55556    0.55556    2.22222    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+DEAL::
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    3.00000    1.50000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    1.50000    3.00000    
+DEAL::
+0.11111    -0.05556   0.11111    0.22222    -0.11111   0.22222    0.00000    0.00000    
+-0.05556   0.11111    0.11111    -0.11111   0.22222    0.22222    0.00000    0.00000    
+0.11111    0.11111    0.44444    0.22222    0.22222    0.88889    0.00000    0.00000    
+0.44444    -0.22222   0.44444    0.55556    -0.27778   0.55556    0.00000    0.00000    
+-0.22222   0.44444    0.44444    -0.27778   0.55556    0.55556    0.00000    0.00000    
+0.44444    0.44444    1.77778    0.55556    0.55556    2.22222    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    0.00000    
+DEAL::
+0.11111    -0.05556   0.11111    0.22222    -0.11111   0.22222    
+-0.05556   0.11111    0.11111    -0.11111   0.22222    0.22222    
+0.11111    0.11111    0.44444    0.22222    0.22222    0.88889    
+0.44444    -0.22222   0.44444    0.55556    -0.27778   0.55556    
+-0.22222   0.44444    0.44444    -0.27778   0.55556    0.55556    
+0.44444    0.44444    1.77778    0.55556    0.55556    2.22222    
+DEAL::