From: Peter Munch Date: Sun, 4 Aug 2024 11:08:29 +0000 (+0200) Subject: Fix MatrixFreeTools::compute_matrix() X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d96ac99db76b2f5a1df9f9bd60ff392434d32006;p=dealii.git Fix MatrixFreeTools::compute_matrix() --- diff --git a/doc/news/changes/minor/20240804Munch b/doc/news/changes/minor/20240804Munch new file mode 100644 index 0000000000..09bfbe1cc0 --- /dev/null +++ b/doc/news/changes/minor/20240804Munch @@ -0,0 +1,5 @@ +Fixed: MatrixFreeTools::compute_matrix() now +correctly handles the case that individual +components are requested. +
+(Peter Munch, 2024/08/04) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index 8b9c1f9f64..e23acadee2 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -42,6 +42,7 @@ namespace MatrixFreeTools std::vector dof_numbers; std::vector quad_numbers; + std::vector n_components; std::vector first_selected_components; std::vector batch_type; static const bool is_face = is_face_; @@ -1971,6 +1972,7 @@ namespace MatrixFreeTools data_cell.dof_numbers = {dof_no}; data_cell.quad_numbers = {quad_no}; + data_cell.n_components = {n_components}; data_cell.first_selected_components = {first_selected_component}; data_cell.batch_type = {0}; @@ -2007,6 +2009,7 @@ namespace MatrixFreeTools data_face.dof_numbers = {dof_no, dof_no}; data_face.quad_numbers = {quad_no, quad_no}; + data_face.n_components = {n_components, n_components}; data_face.first_selected_components = {first_selected_component, first_selected_component}; data_face.batch_type = {1, 2}; @@ -2069,6 +2072,7 @@ namespace MatrixFreeTools data_boundary.dof_numbers = {dof_no}; data_boundary.quad_numbers = {quad_no}; + data_boundary.n_components = {n_components}; data_boundary.first_selected_components = {first_selected_component}; data_boundary.batch_type = {1}; @@ -2155,22 +2159,39 @@ namespace MatrixFreeTools for (unsigned int b = 0; b < n_blocks; ++b) { + const auto &fe = matrix_free.get_dof_handler(data.dof_numbers[b]) + .get_fe(phi[b]->get_active_fe_index()); + + const auto component_base = + matrix_free.get_dof_info(data.dof_numbers[b]) + .component_to_base_index[data.first_selected_components[b]]; + const auto component_in_base = + data.first_selected_components[b] - + matrix_free.get_dof_info(data.dof_numbers[b]) + .start_components[component_base]; + const auto &shape_info = matrix_free.get_shape_info( data.dof_numbers[b], data.quad_numbers[b], - data.first_selected_components[b], + component_base, phi[b]->get_active_fe_index(), phi[b]->get_active_quadrature_index()); dofs_per_cell[b] = - shape_info.dofs_per_component_on_cell * shape_info.n_components; + shape_info.dofs_per_component_on_cell * data.n_components[b]; - dof_indices[b].resize(dofs_per_cell[b]); + dof_indices[b].resize(fe.n_dofs_per_cell()); for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) dof_indices_mf[b][v].resize(dofs_per_cell[b]); - lexicographic_numbering[b] = shape_info.lexicographic_numbering; + lexicographic_numbering[b].insert( + lexicographic_numbering[b].begin(), + shape_info.lexicographic_numbering.begin() + + component_in_base * shape_info.dofs_per_component_on_cell, + shape_info.lexicographic_numbering.begin() + + (component_in_base + data.n_components[b]) * + shape_info.dofs_per_component_on_cell); } for (unsigned int bj = 0; bj < n_blocks; ++bj) @@ -2210,7 +2231,7 @@ namespace MatrixFreeTools else cell_iterator->get_dof_indices(dof_indices[b]); - for (unsigned int j = 0; j < dof_indices[b].size(); ++j) + for (unsigned int j = 0; j < dofs_per_cell[b]; ++j) dof_indices_mf[b][v][j] = dof_indices[b][lexicographic_numbering[b][j]]; } diff --git a/tests/matrix_free/compute_diagonal_10.cc b/tests/matrix_free/compute_diagonal_10.cc new file mode 100644 index 0000000000..2e39129ec6 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_10.cc @@ -0,0 +1,214 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Test MatrixFreeTools::compute_matrix(): +// (1) compute individual components +// (2) compute with multiple FEEvaluation instances +// + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +run(const MatrixFree &matrix_free, + const unsigned int first_selected_component) +{ + const unsigned int n_dofs = matrix_free.get_dof_handler().n_dofs(); + FullMatrix matrix(n_dofs, n_dofs); + + FullMatrix scaling(3, 3); + for (unsigned int i = 0, c = 1; i < 3; ++i) + for (unsigned int j = 0; j < 3; ++j, ++c) + scaling[i][j] = c; + + MatrixFreeTools:: + compute_matrix>( + matrix_free, + matrix_free.get_affine_constraints(), + matrix, + [&](auto &phi) { + phi.evaluate(EvaluationFlags::values); + for (const auto q : phi.quadrature_point_indices()) + { + auto value = phi.get_value(q); + + auto result = value; + result = 0.0; + + for (unsigned int i = 0; i < n_components; ++i) + for (unsigned int j = 0; j < n_components; ++j) + if constexpr (n_components > 1) + result[i] += scaling[i + first_selected_component] + [j + first_selected_component] * + value[j]; + else + result += scaling[i + first_selected_component] + [j + first_selected_component] * + value; + + phi.submit_value(result, q); + } + phi.integrate(EvaluationFlags::values); + }, + 0, + 0, + first_selected_component); + + matrix.print_formatted(deallog.get_file_stream(), 5, false, 10, "0.00000"); + deallog << std::endl; +} + + + +template +void +test_0() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(FESystem(FE_Q(2), 2, FE_Q(1), 1)); + DoFRenumbering::component_wise(dof_handler); + + QGauss quadrature(2); + + MappingQ1 mapping; + + AffineConstraints constraints; + + MatrixFree matrix_free; + + matrix_free.reinit(mapping, dof_handler, constraints, quadrature); + + run<1>(matrix_free, 0u); + run<1>(matrix_free, 1u); + run<1>(matrix_free, 2u); + + run<2>(matrix_free, 0u); +} + + + +template +void +test_1() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(FESystem(FE_Q(2), 2)); + DoFRenumbering::component_wise(dof_handler); + + QGauss quadrature(2); + + MappingQ1 mapping; + + AffineConstraints constraints; + + MatrixFree matrix_free; + + matrix_free.reinit(mapping, dof_handler, constraints, quadrature); + + MatrixFreeTools::internal:: + ComputeMatrixScratchData, false> + data_cell; + + data_cell.dof_numbers = {0, 0}; + data_cell.quad_numbers = {0, 0}; + data_cell.n_components = {1, 1}; + data_cell.first_selected_components = {0, 1}; + data_cell.batch_type = {0, 0}; + + data_cell.op_create = + [&](const std::pair &range) { + std::vector< + std::unique_ptr, false>>> + phi; + + phi.emplace_back(std::make_unique>( + matrix_free, range, 0, 0, 0)); + + phi.emplace_back(std::make_unique>( + matrix_free, range, 0, 0, 1)); + + return phi; + }; + + data_cell.op_reinit = [](auto &phi, const unsigned batch) { + static_cast &>(*phi[0]).reinit(batch); + static_cast &>(*phi[1]).reinit(batch); + }; + + data_cell.op_compute = [&](auto &phi) { + auto &phi_0 = static_cast &>(*phi[0]); + auto &phi_1 = static_cast &>(*phi[1]); + + phi_0.evaluate(EvaluationFlags::values); + phi_1.evaluate(EvaluationFlags::values); + for (const auto q : phi_0.quadrature_point_indices()) + { + auto value_0 = phi_0.get_value(q); + auto value_1 = phi_1.get_value(q); + + phi_0.submit_value(value_0 * 1.0 + value_1 * 2.0, q); + phi_1.submit_value(value_0 * 4.0 + value_1 * 5.0, q); + } + phi_0.integrate(EvaluationFlags::values); + phi_1.integrate(EvaluationFlags::values); + }; + + const unsigned int n_dofs = dof_handler.n_dofs(); + FullMatrix matrix(n_dofs, n_dofs); + + MatrixFreeTools::internal::compute_matrix( + matrix_free, constraints, data_cell, {}, {}, matrix); + + matrix.print_formatted(deallog.get_file_stream(), 5, false, 10, "0.00000"); + deallog << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + test_0<1>(); // individual components + test_1<1>(); // multiple FEEvaluation instances +} diff --git a/tests/matrix_free/compute_diagonal_10.output b/tests/matrix_free/compute_diagonal_10.output new file mode 100644 index 0000000000..5add115c38 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_10.output @@ -0,0 +1,44 @@ + +0.11111 -0.05556 0.11111 0.00000 0.00000 0.00000 0.00000 0.00000 +-0.05556 0.11111 0.11111 0.00000 0.00000 0.00000 0.00000 0.00000 +0.11111 0.11111 0.44444 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL:: +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.55556 -0.27778 0.55556 0.00000 0.00000 +0.00000 0.00000 0.00000 -0.27778 0.55556 0.55556 0.00000 0.00000 +0.00000 0.00000 0.00000 0.55556 0.55556 2.22222 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL:: +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 3.00000 1.50000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.50000 3.00000 +DEAL:: +0.11111 -0.05556 0.11111 0.22222 -0.11111 0.22222 0.00000 0.00000 +-0.05556 0.11111 0.11111 -0.11111 0.22222 0.22222 0.00000 0.00000 +0.11111 0.11111 0.44444 0.22222 0.22222 0.88889 0.00000 0.00000 +0.44444 -0.22222 0.44444 0.55556 -0.27778 0.55556 0.00000 0.00000 +-0.22222 0.44444 0.44444 -0.27778 0.55556 0.55556 0.00000 0.00000 +0.44444 0.44444 1.77778 0.55556 0.55556 2.22222 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL:: +0.11111 -0.05556 0.11111 0.22222 -0.11111 0.22222 +-0.05556 0.11111 0.11111 -0.11111 0.22222 0.22222 +0.11111 0.11111 0.44444 0.22222 0.22222 0.88889 +0.44444 -0.22222 0.44444 0.55556 -0.27778 0.55556 +-0.22222 0.44444 0.44444 -0.27778 0.55556 0.55556 +0.44444 0.44444 1.77778 0.55556 0.55556 2.22222 +DEAL::