From 30153aa43afd5f312539e4bc79949fed46e5e456 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 13 Oct 2020 08:59:26 +0200 Subject: [PATCH] Fix MatrixFree setup for hp::FECollection::size() neq hp::QCollection::size() --- .../matrix_free/mapping_info.templates.h | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 58255da351..85d747d18c 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -799,8 +799,15 @@ namespace internal // FEValues std::vector>>> fe_values(mapping_info.cell_data.size()); + + const unsigned int n_active_fe_indices = + active_fe_index.size() > 0 ? + (*std::max_element(active_fe_index.begin(), active_fe_index.end()) + + 1) : + 1; + for (unsigned int i = 0; i < fe_values.size(); ++i) - fe_values[i].resize(mapping_info.cell_data[i].descriptor.size()); + fe_values[i].resize(n_active_fe_indices); const UpdateFlags update_flags = mapping_info.update_flags_cells; const UpdateFlags update_flags_feval = (update_flags & update_jacobians ? update_jacobians : @@ -823,15 +830,19 @@ namespace internal // fill this data in const unsigned int fe_index = active_fe_index.size() > 0 ? active_fe_index[cell] : 0; - const unsigned int n_q_points = - mapping_info.cell_data[my_q].descriptor[fe_index].n_q_points; + const unsigned int hp_quad_index = + mapping_info.cell_data[my_q].descriptor.size() == 1 ? 0 : + fe_index; + const unsigned int n_q_points = mapping_info.cell_data[my_q] + .descriptor[hp_quad_index] + .n_q_points; if (fe_values[my_q][fe_index].get() == nullptr) fe_values[my_q][fe_index] = std::make_shared>( mapping, dummy_fe, mapping_info.cell_data[my_q] - .descriptor[fe_index] + .descriptor[hp_quad_index] .quadrature, update_flags_feval); dealii::FEValues &fe_val = *fe_values[my_q][fe_index]; -- 2.39.5