From 2a3c73318298c47bdb7b04841de6d86af6498214 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 11 Jul 2021 21:53:55 +0200 Subject: [PATCH] Introduce ShapeInfo::subface_interpolation_matrix --- include/deal.II/matrix_free/shape_info.h | 6 ++ .../matrix_free/shape_info.templates.h | 62 +++++++++++++++++++ 2 files changed, 68 insertions(+) diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 51d61df5cf..ac641152f8 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -261,6 +261,12 @@ namespace internal */ std::array, 2> hessians_within_subface; + /** + * A 1D subface interpolation matrix to the first quadrant. This data + * structure is only set up for FE_Q for dim > 1. + */ + AlignedVector subface_interpolation_matrix; + /** * We store a copy of the one-dimensional quadrature formula * used for initialization. diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 2a9b360878..3c716a4fe4 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -30,9 +30,11 @@ #include #include #include +#include #include #include #include +#include #include #include @@ -171,6 +173,29 @@ namespace internal + template + std::unique_ptr> + create_fe(const FiniteElement &fe) + { + std::string fe_name = fe.get_name(); + + Assert( + fe_name.find("FESystem") == std::string::npos, + ExcMessage( + "This function can not accept FESystem but only base elements.")); + + { + const std::size_t template_starts = fe_name.find_first_of('<'); + Assert(fe_name[template_starts + 1] == + (dim == 1 ? '1' : (dim == 2 ? '2' : '3')), + ExcInternalError()); + fe_name[template_starts + 1] = std::to_string(dim_to).c_str()[0]; + } + return FETools::get_fe_by_name(fe_name); + } + + + template ShapeInfo::ShapeInfo() : element_type(tensor_general) @@ -600,6 +625,43 @@ namespace internal } } + if (dim > 1 && dynamic_cast *>(&fe)) + { + auto &subface_interpolation_matrix = + univariate_shape_data.subface_interpolation_matrix; + + const auto fe_1d = create_fe<1>(fe); + const auto fe_2d = create_fe<2>(fe); + + FullMatrix interpolation_matrix(fe_2d->n_dofs_per_face(0), + fe_2d->n_dofs_per_face(0)); + + fe_2d->get_subface_interpolation_matrix(*fe_2d, + 0, + interpolation_matrix, + 0); + + ElementType element_type; + std::vector scalar_lexicographic; + std::vector lexicographic_numbering; + + get_element_type_specific_information(*fe_1d, + *fe_1d, + 0, + element_type, + scalar_lexicographic, + lexicographic_numbering); + + subface_interpolation_matrix.resize(fe_1d->n_dofs_per_cell() * + fe_1d->n_dofs_per_cell()); + + for (unsigned int i = 0, c = 0; i < fe_1d->n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < fe_1d->n_dofs_per_cell(); ++j, ++c) + subface_interpolation_matrix[c] = + interpolation_matrix(scalar_lexicographic[i], + scalar_lexicographic[j]); + } + // get gradient and Hessian transformation matrix for the polynomial // space associated with the quadrature rule (collocation space). We // need to avoid the case with more than a few hundreds of quadrature -- 2.39.5