From 6c82f030b5d6040f6fca1a64c23bc8a63f9144ed Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 1 Nov 2022 22:22:37 +0100 Subject: [PATCH] Fast hanging-node algorithm for FE_Q_iso_Q1 --- .../matrix_free/hanging_nodes_internal.h | 11 +++++--- .../matrix_free/shape_info.templates.h | 25 +++++++++++++++---- 2 files changed, 28 insertions(+), 8 deletions(-) diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h index b5cc0cf6bb..3942abed73 100644 --- a/include/deal.II/matrix_free/hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -24,6 +24,7 @@ #include #include +#include #include #include @@ -491,9 +492,13 @@ namespace internal for (unsigned int c = 0; c < fe_collection[i].element_multiplicity(base_element_index); ++c, ++comp) - if (dim == 1 || dynamic_cast *>( - &fe_collection[i].base_element( - base_element_index)) == nullptr) + if (dim == 1 || + (dynamic_cast *>( + &fe_collection[i].base_element(base_element_index)) == + nullptr && + dynamic_cast *>( + &fe_collection[i].base_element(base_element_index)) == + nullptr)) supported_components[i][comp] = false; else supported_components[i][comp] = true; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index b9a7ba4419..3c73f22759 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -753,7 +753,8 @@ namespace internal } } - if (dim > 1 && dynamic_cast *>(&fe)) + if (dim > 1 && (dynamic_cast *>(&fe) || + dynamic_cast *>(&fe))) { auto &subface_interpolation_matrix_0 = univariate_shape_data.subface_interpolation_matrices[0]; @@ -770,19 +771,33 @@ namespace internal subface_interpolation_matrix_scalar_0.resize(nn * nn); subface_interpolation_matrix_scalar_1.resize(nn * nn); - std::vector> fe_q_points = QGaussLobatto<1>(nn).get_points(); - const std::vector> poly = + const bool is_feq = dynamic_cast *>(&fe) != nullptr; + + std::vector> fe_q_points = + is_feq ? QGaussLobatto<1>(nn).get_points() : + QIterated<1>(QTrapezoid<1>(), nn - 1).get_points(); + + const std::vector> poly_feq = Polynomials::generate_complete_Lagrange_basis(fe_q_points); + const std::vector> + poly_feq_iso_q1 = + Polynomials::generate_complete_Lagrange_basis_on_subdivisions(nn - + 1, + 1); + for (unsigned int i = 0, c = 0; i < nn; ++i) for (unsigned int j = 0; j < nn; ++j, ++c) { subface_interpolation_matrix_scalar_0[c] = - poly[j].value(0.5 * fe_q_points[i][0]); + is_feq ? poly_feq[j].value(0.5 * fe_q_points[i][0]) : + poly_feq_iso_q1[j].value(0.5 * fe_q_points[i][0]); subface_interpolation_matrix_0[c] = subface_interpolation_matrix_scalar_0[c]; subface_interpolation_matrix_scalar_1[c] = - poly[j].value(0.5 + 0.5 * fe_q_points[i][0]); + is_feq ? + poly_feq[j].value(0.5 + 0.5 * fe_q_points[i][0]) : + poly_feq_iso_q1[j].value(0.5 + 0.5 * fe_q_points[i][0]); subface_interpolation_matrix_1[c] = subface_interpolation_matrix_scalar_1[c]; } -- 2.39.5