*/
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
{
/**
* 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)"));
}
+ // 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>
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> &,
#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>
deallog << std::endl;
}
- // Compute diagonal via MatrixFreeTools
+ // Compute diagonal via MatrixFreeTools::compute_diagonal
VectorType diagonal_mft;
matrix_free.initialize_dof_vector(diagonal_mft);
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;
}