From aabf6bb25ef08d4302e848e5472adb92c212d273 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Mon, 18 Nov 2024 14:26:27 +0100 Subject: [PATCH] Fix bug in MatrixFreeTools::compute_diagonal for FE_Nothing --- include/deal.II/matrix_free/tools.h | 26 +++-- tests/matrix_free/compute_diagonal_12.cc | 101 ++++++++++++++++++ ...iagonal_12.with_p4est=true.mpirun=1.output | 5 + 3 files changed, 124 insertions(+), 8 deletions(-) create mode 100644 tests/matrix_free/compute_diagonal_12.cc create mode 100644 tests/matrix_free/compute_diagonal_12.with_p4est=true.mpirun=1.output diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index c68490128d..597bbb2f71 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -1700,7 +1700,8 @@ namespace MatrixFreeTools }; data_cell.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); + if (phi.size() == 1) + static_cast(*phi[0]).reinit(batch); }; if (cell_operation) @@ -1761,8 +1762,11 @@ namespace MatrixFreeTools }; data_face.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); - static_cast(*phi[1]).reinit(batch); + if (phi.size() == 2) + { + static_cast(*phi[0]).reinit(batch); + static_cast(*phi[1]).reinit(batch); + } }; if (face_operation) @@ -1801,7 +1805,8 @@ namespace MatrixFreeTools }; data_boundary.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); + if (phi.size() == 1) + static_cast(*phi[0]).reinit(batch); }; if (boundary_operation) @@ -2217,7 +2222,8 @@ namespace MatrixFreeTools }; data_cell.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); + if (phi.size() == 1) + static_cast(*phi[0]).reinit(batch); }; if (cell_operation) @@ -2278,8 +2284,11 @@ namespace MatrixFreeTools }; data_face.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); - static_cast(*phi[1]).reinit(batch); + if (phi.size() == 2) + { + static_cast(*phi[0]).reinit(batch); + static_cast(*phi[1]).reinit(batch); + } }; if (face_operation) @@ -2318,7 +2327,8 @@ namespace MatrixFreeTools }; data_boundary.op_reinit = [](auto &phi, const unsigned batch) { - static_cast(*phi[0]).reinit(batch); + if (phi.size() == 1) + static_cast(*phi[0]).reinit(batch); }; if (boundary_operation) diff --git a/tests/matrix_free/compute_diagonal_12.cc b/tests/matrix_free/compute_diagonal_12.cc new file mode 100644 index 0000000000..ea88edb557 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_12.cc @@ -0,0 +1,101 @@ +// ------------------------------------------------------------------------ +// +// 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. +// +// ------------------------------------------------------------------------ + + +// Based on compute_diagonal_04 but uses FE_Nothing + +#include "compute_diagonal_util.h" + +template > +void +test() +{ + const unsigned int n_components = 1; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tria, -numbers::PI / 2, numbers::PI / 2); + + tria.refine_global(); + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_active() && cell->is_locally_owned() && + cell->center()[0] < 0.0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + + hp::FECollection fes{FE_Q(1), FE_Nothing(1)}; + + // setup dof-handlers + DoFHandler dof_handler(tria); + + unsigned int counter = 0; + + for (const auto &cell : dof_handler.active_cell_iterators()) + cell->set_active_fe_index(++counter % 2); + + dof_handler.distribute_dofs(fes); + + AffineConstraints constraint; + DoFTools::make_hanging_node_constraints(dof_handler, constraint); + + VectorTools::interpolate_boundary_values( + dof_handler, 0, Functions::ZeroFunction(n_components), constraint); + + constraint.close(); + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + + MappingQ mapping(1); + QGauss<1> quad(2); + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + LinearAlgebra::distributed::Vector diagonal; + matrix_free.initialize_dof_vector(diagonal); + + dealii::MatrixFreeTools::compute_diagonal< + dim, + -1, + 0, + n_components, + Number, + VectorizedArrayType, + LinearAlgebra::distributed::Vector>( + matrix_free, + diagonal, + [](FEEvaluation + &phi) { + phi.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + }); + + diagonal.print(deallog.get_file_stream()); +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2>(); +} diff --git a/tests/matrix_free/compute_diagonal_12.with_p4est=true.mpirun=1.output b/tests/matrix_free/compute_diagonal_12.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..4a489c4bb9 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_12.with_p4est=true.mpirun=1.output @@ -0,0 +1,5 @@ + +Process #0 +Local range: [0, 12), global size: 12 +Vector data: +2.167e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.333e+00 1.333e+00 1.333e+00 1.333e+00 0.000e+00 0.000e+00 -- 2.39.5