From: Maximilian Bergbauer Date: Fri, 17 Jan 2025 16:12:01 +0000 (+0100) Subject: MatrixFree: enable zeroing in loops if assignment operator is available for VectorType X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cc36aae4c506a12128a01f11a50ede5efaa36873;p=dealii.git MatrixFree: enable zeroing in loops if assignment operator is available for VectorType --- diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 4eecd5539a..c4802d352d 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -3940,8 +3940,9 @@ namespace internal */ template , VectorType> - * = nullptr, - typename VectorType::value_type * = nullptr> + * = nullptr, + std::enable_if_t, 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)")); } diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h index 9f28fcbd34..b70ecf6d0e 100644 --- a/include/deal.II/matrix_free/type_traits.h +++ b/include/deal.II/matrix_free/type_traits.h @@ -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 + using assignment_operator_t = + decltype(std::declval().operator=(typename T::value_type())); + + template + constexpr bool has_assignment_operator = + is_supported_operation; + + + // a helper type-trait that leverage SFINAE to figure out if type T has // T::communication_block_size template diff --git a/source/lac/affine_constraints.inst.in b/source/lac/affine_constraints.inst.in index bffe8f17bc..d5de48a6ff 100644 --- a/source/lac/affine_constraints.inst.in +++ b/source/lac/affine_constraints.inst.in @@ -138,6 +138,12 @@ for (S : REAL_AND_COMPLEX_SCALARS; T : DEAL_II_VEC_TEMPLATES) const std::vector &, DiagonalMatrix> &) const; + template void AffineConstraints::distribute_local_to_global< + DiagonalMatrix>>(const FullMatrix &, + const std::vector &, + const std::vector &, + DiagonalMatrix> &) const; + template void AffineConstraints::distribute_local_to_global< DiagonalMatrix>>( const FullMatrix &, diff --git a/tests/matrix_free/compute_diagonal_08.cc b/tests/matrix_free/compute_diagonal_08.cc index 709a18b791..1510a8e0f2 100644 --- a/tests/matrix_free/compute_diagonal_08.cc +++ b/tests/matrix_free/compute_diagonal_08.cc @@ -20,6 +20,8 @@ #include +#include + #include #include #include @@ -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 diagonal_matrix(diagonal_2); + + MatrixFreeTools::compute_matrix(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; }