From 5c665959975b0c0a11e4b60738d0db70f5a073b1 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 9 Feb 2024 19:53:23 +0100 Subject: [PATCH] FEEvaluation + DG + hp: fix quadrature indexing --- .../20240209BergbauerMunchRitthalerSchreter | 3 + include/deal.II/matrix_free/fe_evaluation.h | 9 +- tests/matrix_free/hp_dg_01.cc | 98 +++++++++++++++++++ tests/matrix_free/hp_dg_01.output | 2 + 4 files changed, 109 insertions(+), 3 deletions(-) create mode 100644 doc/news/changes/minor/20240209BergbauerMunchRitthalerSchreter create mode 100644 tests/matrix_free/hp_dg_01.cc create mode 100644 tests/matrix_free/hp_dg_01.output diff --git a/doc/news/changes/minor/20240209BergbauerMunchRitthalerSchreter b/doc/news/changes/minor/20240209BergbauerMunchRitthalerSchreter new file mode 100644 index 0000000000..65a61f92b5 --- /dev/null +++ b/doc/news/changes/minor/20240209BergbauerMunchRitthalerSchreter @@ -0,0 +1,3 @@ +Fixed: Fix initialization of FEFaceEvaluation for hp in 3D. +
+(Maximilian Bergbauer, Peter Munch, Andreas Ritthaler, Magdalena Schreter-Fleischhacker, 2024/02/09) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 4cd421eadd..ad32c9e623 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2858,13 +2858,16 @@ namespace internal (active_fe_index_given != numbers::invalid_unsigned_int ? active_fe_index_given : 0); + init_data.active_quad_index = fe_degree == numbers::invalid_unsigned_int ? (active_quad_index_given != numbers::invalid_unsigned_int ? active_quad_index_given : - std::min(init_data.active_fe_index, - init_data.mapping_data->descriptor.size() - - 1)) : + std::min( + init_data.active_fe_index, + init_data.mapping_data->descriptor.size() / + (is_face ? std::max(1, dim - 1) : 1) - + 1)) : init_data.mapping_data->quad_index_from_n_q_points(n_q_points); init_data.shape_info = &matrix_free.get_shape_info( diff --git a/tests/matrix_free/hp_dg_01.cc b/tests/matrix_free/hp_dg_01.cc new file mode 100644 index 0000000000..00fbc13d01 --- /dev/null +++ b/tests/matrix_free/hp_dg_01.cc @@ -0,0 +1,98 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check if intialization of FEFaceEvaluation works for hp + DG in 3D. + +#include +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test() +{ + using VectorType = Vector; + + Triangulation tria; + + std::vector repetitions(dim, 1); + repetitions[0] = 2; + + Point p1; + Point p2; + + p2[0] = 2.0; + for (unsigned int i = 1; i < dim; ++i) + p2[i] = 1.0; + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); + + + hp::FECollection fe; + fe.push_back(FE_DGQ(1)); + fe.push_back(FE_DGQ(1)); + + hp::QCollection quad; + quad.push_back(QGauss(2)); + + MappingQ1 mapping; + + DoFHandler dof_handler(tria); + + for (const auto &cell : dof_handler.active_cell_iterators()) + cell->set_active_fe_index(cell->active_cell_index()); + + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + + typename MatrixFree::AdditionalData ad; + + ad.mapping_update_flags_inner_faces = update_gradients; + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraints, quad, ad); + + VectorType dst(dof_handler.n_dofs()); + VectorType src(dof_handler.n_dofs()); + + matrix_free.template loop( + [](const auto &, auto &, const auto &, const auto) {}, + [](const auto &matrix_free, auto &, const auto &, const auto range) { + FEFaceEvaluation eval(matrix_free, range); + }, + [](const auto &, auto &, const auto &, const auto) {}, + dst, + src); +} + +int +main() +{ + initlog(); + + test<3>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/matrix_free/hp_dg_01.output b/tests/matrix_free/hp_dg_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/matrix_free/hp_dg_01.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5