const unsigned int first_selected_component = 0);
+ /**
+ * Compute the matrix representation of a linear operator (@p matrix), given
+ * @p matrix_free and the local cell integral operation @p local_vmult. The
+ * vector is initialized to the right size in the function.
+ *
+ * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+ * passed to the constructor of the FEEvaluation that is internally set up.
+ */
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> & matrix_free,
+ const AffineConstraints<Number> & constraints,
+ MatrixType & matrix,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)> &local_vmult,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Same as above but with a class and a function pointer.
+ */
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> & constraints,
+ MatrixType & matrix,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ const CLASS * owning_class,
+ const unsigned int dof_no = 0,
+ const unsigned int quad_no = 0,
+ const unsigned int first_selected_component = 0);
+
+
// implementations
false);
}
+ template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> & matrix_free,
+ const AffineConstraints<Number> & constraints,
+ MatrixType & matrix,
+ const std::function<void(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &)> &local_vmult,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ matrix_free.template cell_loop<MatrixType, MatrixType>(
+ [&](const auto &, auto &dst, const auto &, const auto range) {
+ FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType>
+ integrator(
+ matrix_free, range, dof_no, quad_no, first_selected_component);
+
+ unsigned int const dofs_per_cell = integrator.dofs_per_cell;
+
+ std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+ std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
+
+ std::array<FullMatrix<typename MatrixType::value_type>,
+ VectorizedArrayType::size()>
+ matrices;
+
+ std::fill_n(matrices.begin(),
+ VectorizedArrayType::size(),
+ FullMatrix<typename MatrixType::value_type>(dofs_per_cell,
+ dofs_per_cell));
+
+ const auto lexicographic_numbering =
+ matrix_free
+ .get_shape_info(dof_no,
+ quad_no,
+ first_selected_component,
+ integrator.get_active_fe_index(),
+ integrator.get_active_quadrature_index())
+ .lexicographic_numbering;
+
+ for (auto cell = range.first; cell < range.second; ++cell)
+ {
+ integrator.reinit(cell);
+
+ unsigned int const n_filled_lanes =
+ matrix_free.n_active_entries_per_cell_batch(cell);
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ matrices[v] = 0.0;
+
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ integrator.begin_dof_values()[i] =
+ static_cast<Number>(i == j);
+
+ local_vmult(integrator);
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ matrices[v](i, j) = integrator.begin_dof_values()[i][v];
+ }
+
+ for (unsigned int v = 0; v < n_filled_lanes; ++v)
+ {
+ const auto cell_v = matrix_free.get_cell_iterator(cell, v);
+
+ if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
+ cell_v->get_mg_dof_indices(dof_indices);
+ else
+ cell_v->get_dof_indices(dof_indices);
+
+ for (unsigned int j = 0; j < dof_indices.size(); ++j)
+ dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
+
+ constraints.distribute_local_to_global(matrices[v],
+ dof_indices_mf,
+ dof_indices_mf,
+ dst);
+ }
+ }
+ },
+ matrix,
+ matrix);
+
+ matrix.compress(VectorOperation::add);
+
+ for (unsigned int i = 0; i < matrix.m(); ++i)
+ if (matrix.get_row_length(i) > 0 && matrix(i, i) == 0.0 &&
+ constraints.is_constrained(i))
+ matrix.add(i, i, 1);
+ }
+
+ template <typename CLASS,
+ int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ int n_components,
+ typename Number,
+ typename VectorizedArrayType,
+ typename MatrixType>
+ void
+ compute_matrix(
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> & constraints,
+ MatrixType & matrix,
+ void (CLASS::*cell_operation)(FEEvaluation<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType> &) const,
+ const CLASS * owning_class,
+ const unsigned int dof_no,
+ const unsigned int quad_no,
+ const unsigned int first_selected_component)
+ {
+ compute_matrix<dim,
+ fe_degree,
+ n_q_points_1d,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ MatrixType>(matrix_free,
+ constraints,
+ matrix,
+ [&](auto &feeval) {
+ (owning_class->*cell_operation)(feeval);
+ },
+ dof_no,
+ quad_no,
+ first_selected_component);
+ }
+
} // namespace MatrixFreeTools
DEAL_II_NAMESPACE_CLOSE
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_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>
using VectorType = LinearAlgebra::distributed::Vector<Number>;
public:
- Test(const MatrixFree<dim, Number, VectorizedArrayType> matrix_free,
+ Test(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+ const AffineConstraints<Number> & constraints,
const std::function<void(FEEvaluation<dim,
fe_degree,
n_points,
VectorizedArrayType> &)>
&cell_operation)
: matrix_free(matrix_free)
+ , constraints(constraints)
, cell_operation(cell_operation)
{}
do_test()
{
// compute diagonal with the new function
+ VectorType diagonal_global, diagonal_global_reference;
+
+ SparseMatrix<Number> A1, A2, A_ref;
+ SparsityPattern sparsity_pattern;
+
+ const bool test_matrix = Utilities::MPI::n_mpi_processes(
+ matrix_free.get_task_info().communicator) == 1;
+
+ if (test_matrix)
+ {
+ DynamicSparsityPattern dsp(matrix_free.get_dof_handler().n_dofs());
+ DoFTools::make_sparsity_pattern(matrix_free.get_dof_handler(),
+ dsp,
+ constraints);
+ sparsity_pattern.copy_from(dsp);
+ A1.reinit(sparsity_pattern);
+ A2.reinit(sparsity_pattern);
+ A_ref.reinit(sparsity_pattern);
+ }
+
{
- VectorType diagonal_global;
MatrixFreeTools::compute_diagonal<dim,
fe_degree,
n_points,
deallog << diagonal_global.l2_norm() << std::endl;
}
+ if (test_matrix)
+ {
+ MatrixFreeTools::compute_matrix<dim,
+ fe_degree,
+ n_points,
+ n_components,
+ Number,
+ VectorizedArrayType,
+ SparseMatrix<Number>>(
+ matrix_free, constraints, A1, [&](auto &phi) {
+ this->cell_operation(phi);
+ });
+ }
+
+ if (test_matrix)
+ {
+ MatrixFreeTools::compute_matrix(
+ matrix_free, constraints, A2, &Test::cell_function, this);
+ }
+
// compute diagonal globally
{
- VectorType src, dst, temp;
+ VectorType src, temp;
matrix_free.initialize_dof_vector(src);
- matrix_free.initialize_dof_vector(dst);
+ matrix_free.initialize_dof_vector(diagonal_global_reference);
matrix_free.initialize_dof_vector(temp);
for (unsigned int i = 0; i < src.size(); i++)
if (src.get_partitioner()->in_local_range(i))
{
- dst[i] = temp[i];
- src[i] = 0.0;
+ diagonal_global_reference[i] = temp[i];
+ src[i] = 0.0;
}
+
+ if (test_matrix)
+ {
+ for (unsigned int j = 0; j < src.size(); j++)
+ if (temp[j] != 0.0)
+ A_ref(j, i) = temp[j];
+ else if (i == j)
+ A_ref(j, i) = 1.0;
+ }
+
temp = 0.0;
}
- dst.print(deallog.get_file_stream());
- deallog << dst.l2_norm() << std::endl;
+ diagonal_global_reference.print(deallog.get_file_stream());
+ deallog << diagonal_global_reference.l2_norm() << std::endl;
}
+
+ if (test_matrix)
+ {
+ auto a1 = A1.begin();
+ auto a2 = A2.begin();
+ auto a_ref = A_ref.begin();
+
+ for (; a1 != A1.end(); ++a1, ++a2, ++a_ref)
+ {
+ Assert(std::abs(a1->value() - a_ref->value()) < 1e-6,
+ ExcNotImplemented());
+ Assert(std::abs(a2->value() - a_ref->value()) < 1e-6,
+ ExcNotImplemented());
+ }
+ }
}
void
}
}
- const MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+ void
+ cell_function(FEEvaluation<dim,
+ fe_degree,
+ n_points,
+ n_components,
+ Number,
+ VectorizedArrayType> &phi) const
+ {
+ this->cell_operation(phi);
+ }
+
+ const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
+ const AffineConstraints<Number> & constraints;
const std::function<void(FEEvaluation<dim,
fe_degree,
n_points,
projects / dealii.git / commitdiff