]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MatrixFree: enable zeroing in loops if assignment operator is available for VectorType 18008/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Fri, 17 Jan 2025 16:12:01 +0000 (17:12 +0100)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 28 Jan 2025 13:04:29 +0000 (14:04 +0100)
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/type_traits.h
source/lac/affine_constraints.inst.in
tests/matrix_free/compute_diagonal_08.cc

index 4eecd5539a14f71efe924ebdaa66e4288a87f2cb..c4802d352d4f3732f262c839cf2abc7db2689c0d 100644 (file)
@@ -3940,8 +3940,9 @@ namespace internal
      */
     template <typename VectorType,
               std::enable_if_t<!has_exchange_on_subset<VectorType>, VectorType>
-                                              * = nullptr,
-              typename VectorType::value_type * = nullptr>
+                * = nullptr,
+              std::enable_if_t<has_assignment_operator<VectorType>, VectorType>
+                * = nullptr>
     void
     zero_vector_region(const unsigned int range_index, VectorType &vec) const
     {
@@ -3963,14 +3964,15 @@ namespace internal
 
     /**
      * Zero out vector region for non-vector types, i.e., classes that do not
-     * have VectorType::value_type
+     * have operator=(const VectorType::value_type)
      */
     void
     zero_vector_region(const unsigned int, ...) const
     {
       Assert(false,
-             ExcNotImplemented("Zeroing is only implemented for vector types "
-                               "which provide VectorType::value_type"));
+             ExcNotImplemented(
+               "Zeroing is only implemented for vector types "
+               "which provide operator=(const VectorType::value_type)"));
     }
 
 
index 9f28fcbd345d1d190d3384b6c53da465a81ba7a1..b70ecf6d0e34e5b13564ef9c7d085b44f2d70405 100644 (file)
@@ -146,6 +146,18 @@ namespace internal
 
 
 
+  // a helper type-trait that leverage SFINAE to figure out if type T has
+  // T & T::operator=(const T::value_type) const
+  template <typename T>
+  using assignment_operator_t =
+    decltype(std::declval<T>().operator=(typename T::value_type()));
+
+  template <typename T>
+  constexpr bool has_assignment_operator =
+    is_supported_operation<assignment_operator_t, T>;
+
+
+
   // a helper type-trait that leverage SFINAE to figure out if type T has
   // T::communication_block_size
   template <typename T>
index bffe8f17bc9c11f3d173b9b7541f6707e784ebf1..d5de48a6ffb5caa8dbefe7a69429ce8f73690703 100644 (file)
@@ -138,6 +138,12 @@ for (S : REAL_AND_COMPLEX_SCALARS; T : DEAL_II_VEC_TEMPLATES)
       const std::vector<size_type> &,
       DiagonalMatrix<LinearAlgebra::distributed::T<S>> &) const;
 
+    template void AffineConstraints<S>::distribute_local_to_global<
+      DiagonalMatrix<T<S>>>(const FullMatrix<S> &,
+                            const std::vector<size_type> &,
+                            const std::vector<size_type> &,
+                            DiagonalMatrix<T<S>> &) const;
+
     template void AffineConstraints<S>::distribute_local_to_global<
       DiagonalMatrix<LinearAlgebra::distributed::T<S>>>(
       const FullMatrix<S> &,
index 709a18b791136d688cf0de62a56235281332b9b0..1510a8e0f2e24b00359d9d9630e5142fa7078ffd 100644 (file)
@@ -20,6 +20,8 @@
 
 #include <deal.II/grid/grid_generator.h>
 
+#include <deal.II/lac/diagonal_matrix.h>
+
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 #include <deal.II/matrix_free/tools.h>
@@ -238,7 +240,7 @@ test()
           deallog << std::endl;
         }
 
-  // Compute diagonal via MatrixFreeTools
+  // Compute diagonal via MatrixFreeTools::compute_diagonal
   VectorType diagonal_mft;
   matrix_free.initialize_dof_vector(diagonal_mft);
 
@@ -263,6 +265,35 @@ test()
         deallog << std::endl;
       }
 
+  // Compute diagonal via MatrixFreeTools::compute_matrix
+  VectorType diagonal_2;
+  matrix_free.initialize_dof_vector(diagonal_2);
+  DiagonalMatrix<VectorType> diagonal_matrix(diagonal_2);
+
+  MatrixFreeTools::compute_matrix<dim,
+                                  fe_degree,
+                                  n_points,
+                                  n_components,
+                                  Number,
+                                  VectorizedArrayType>(matrix_free,
+                                                       constraints,
+                                                       diagonal_matrix,
+                                                       cell_operation,
+                                                       face_operation,
+                                                       boundary_operation);
+
+  diagonal_2 = diagonal_matrix.get_vector();
+
+  for (unsigned int i = 0; i < diagonal_2.size(); ++i)
+    if (std::abs(diagonal[i] - diagonal_2[i]) > 1e-6)
+      {
+        diagonal.print(deallog.get_file_stream());
+        deallog << std::endl;
+
+        diagonal_2.print(deallog.get_file_stream());
+        deallog << std::endl;
+      }
+
   deallog << "OK!" << std::endl;
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.