From e2f6dae265ece7ea01cc5f9ff7c5df2006e13842 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 6 Feb 2021 20:43:11 +0100
Subject: [PATCH] MatrixFreeTools::compute_matrix() make it work for float

---
 include/deal.II/matrix_free/tools.h           | 62 ++++++++++++-
 tests/matrix_free/compute_diagonal_01.cc      | 10 ++-
 ...h_mpi=true.with_p4est=true.mpirun=1.output | 30 +++++++
 ...h_mpi=true.with_p4est=true.mpirun=3.output | 90 +++++++++++++++++++
 4 files changed, 187 insertions(+), 5 deletions(-)

diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h
index 99498540c7..55376a701b 100644
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -628,6 +628,55 @@ namespace MatrixFreeTools
       first_selected_component);
   }
 
+  namespace internal
+  {
+    /**
+     * If value type of matrix and constrains equals, return a reference
+     * to the given AffineConstraint instance.
+     */
+    template <typename MatrixType,
+              typename Number,
+              typename std::enable_if<std::is_same<
+                typename std::remove_const<typename std::remove_reference<
+                  typename MatrixType::value_type>::type>::type,
+                typename std::remove_const<typename std::remove_reference<
+                  Number>::type>::type>::value>::type * = nullptr>
+    const AffineConstraints<typename MatrixType::value_type> &
+    create_new_affine_constraints_if_needed(
+      const MatrixType &,
+      const AffineConstraints<Number> &constraints,
+      std::unique_ptr<AffineConstraints<typename MatrixType::value_type>> &)
+    {
+      return constraints;
+    }
+
+    /**
+     * If value type of matrix and constrains do not equal, a new
+     * AffineConstraint instance with the value type of the matrix is
+     * created and a reference to it is returned.
+     */
+    template <typename MatrixType,
+              typename Number,
+              typename std::enable_if<!std::is_same<
+                typename std::remove_const<typename std::remove_reference<
+                  typename MatrixType::value_type>::type>::type,
+                typename std::remove_const<typename std::remove_reference<
+                  Number>::type>::type>::value>::type * = nullptr>
+    const AffineConstraints<typename MatrixType::value_type> &
+    create_new_affine_constraints_if_needed(
+      const MatrixType &,
+      const AffineConstraints<Number> &constraints,
+      std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
+        &new_constraints)
+    {
+      new_constraints =
+        std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
+      new_constraints->copy_from(constraints);
+
+      return *new_constraints;
+    }
+  } // namespace internal
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -637,9 +686,9 @@ namespace MatrixFreeTools
             typename MatrixType>
   void
   compute_matrix(
-    const MatrixFree<dim, Number, VectorizedArrayType> &            matrix_free,
-    const AffineConstraints<Number> &                               constraints,
-    MatrixType &                                                    matrix,
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number> &                   constraints_in,
+    MatrixType &                                        matrix,
     const std::function<void(FEEvaluation<dim,
                                           fe_degree,
                                           n_q_points_1d,
@@ -650,6 +699,13 @@ namespace MatrixFreeTools
     const unsigned int                                              quad_no,
     const unsigned int first_selected_component)
   {
+    std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
+                                                              constraints_for_matrix;
+    const AffineConstraints<typename MatrixType::value_type> &constraints =
+      internal::create_new_affine_constraints_if_needed(matrix,
+                                                        constraints_in,
+                                                        constraints_for_matrix);
+
     matrix_free.template cell_loop<MatrixType, MatrixType>(
       [&](const auto &, auto &dst, const auto &, const auto range) {
         FEEvaluation<dim,
diff --git a/tests/matrix_free/compute_diagonal_01.cc b/tests/matrix_free/compute_diagonal_01.cc
index 1a9d46ba9c..ffbdedf6b1 100644
--- a/tests/matrix_free/compute_diagonal_01.cc
+++ b/tests/matrix_free/compute_diagonal_01.cc
@@ -47,8 +47,11 @@ test()
   AffineConstraints<Number> constraint;
   DoFTools::make_hanging_node_constraints(dof_handler, constraint);
 
-  VectorTools::interpolate_boundary_values(
-    dof_handler, 0, Functions::ZeroFunction<dim>(n_components), constraint);
+  VectorTools::interpolate_boundary_values(dof_handler,
+                                           0,
+                                           Functions::ZeroFunction<dim, Number>(
+                                             n_components),
+                                           constraint);
 
   constraint.close();
 
@@ -89,4 +92,7 @@ main(int argc, char **argv)
 
   test<2, 1, 2, 1>(); // scalar
   test<2, 1, 2, 2>(); // vector
+
+  test<2, 1, 2, 1, float>(); // scalar
+  test<2, 1, 2, 2, float>(); // vector
 }
diff --git a/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=1.output
index e0e5c3a0a5..60689ac866 100644
--- a/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=1.output
+++ b/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=1.output
@@ -29,3 +29,33 @@ Local range: [0, 36), global size: 36
 Vector data:
 0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
 DEAL:0::7.78888
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::7.78888
+Process #0
+Local range: [0, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::7.78888
+Process #0
+Local range: [0, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 3.000e+00 3.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:0::7.78888
diff --git a/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=3.output
index 7099f7da5c..8fdfe585ea 100644
--- a/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=3.output
+++ b/tests/matrix_free/compute_diagonal_01.with_mpi=true.with_p4est=true.mpirun=3.output
@@ -1,4 +1,34 @@
 
+Process #0
+Local range: [0, 9), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 3.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 9), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 3.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 9), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 3.000e+00 
+DEAL:0::5.50757
+Process #0
+Local range: [0, 18), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 3.000e+00 3.000e+00 
+DEAL:0::7.78888
+Process #0
+Local range: [0, 18), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 3.000e+00 3.000e+00 
+DEAL:0::7.78888
+Process #0
+Local range: [0, 18), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 3.000e+00 3.000e+00 
+DEAL:0::7.78888
 Process #0
 Local range: [0, 9), global size: 18
 Vector data:
@@ -60,6 +90,36 @@ Local range: [18, 22), global size: 36
 Vector data:
 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
 DEAL:1::7.78888
+Process #1
+Local range: [9, 11), global size: 18
+Vector data:
+0.000e+00 0.000e+00 
+DEAL:1::5.50757
+Process #1
+Local range: [9, 11), global size: 18
+Vector data:
+0.000e+00 0.000e+00 
+DEAL:1::5.50757
+Process #1
+Local range: [9, 11), global size: 18
+Vector data:
+0.000e+00 0.000e+00 
+DEAL:1::5.50757
+Process #1
+Local range: [18, 22), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:1::7.78888
+Process #1
+Local range: [18, 22), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:1::7.78888
+Process #1
+Local range: [18, 22), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:1::7.78888
 
 
 Process #2
@@ -92,4 +152,34 @@ Local range: [22, 36), global size: 36
 Vector data:
 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
 DEAL:2::7.78888
+Process #2
+Local range: [11, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::5.50757
+Process #2
+Local range: [11, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::5.50757
+Process #2
+Local range: [11, 18), global size: 18
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::5.50757
+Process #2
+Local range: [22, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::7.78888
+Process #2
+Local range: [22, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::7.78888
+Process #2
+Local range: [22, 36), global size: 36
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.667e+00 2.667e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 
+DEAL:2::7.78888
 
-- 
2.39.5