From 382733684f2a68ff236f13ec6005987dd930a31c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 1 Dec 2020 16:04:18 +0100 Subject: [PATCH] Fix lexicographic_numbering for vector-valued Simplex::FE_Q and FE_DGQ --- .../matrix_free/shape_info.templates.h | 192 ++++++++++-------- 1 file changed, 110 insertions(+), 82 deletions(-) diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 90f7cea9c3..7b617551c1 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -77,6 +77,96 @@ namespace internal } + template + void + get_element_type_specific_information( + const FiniteElement &fe_in, + const FiniteElement *fe, + const unsigned int base_element_number, + ElementType & element_type, + std::vector & scalar_lexicographic, + std::vector & lexicographic_numbering) + { + element_type = tensor_general; + + const auto fe_poly = dynamic_cast *>(fe); + + if (dynamic_cast *>(fe) != nullptr || + dynamic_cast *>(fe) != nullptr) + { + scalar_lexicographic.resize(fe->n_dofs_per_cell()); + for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i) + scalar_lexicographic[i] = i; + element_type = tensor_none; + } + else if (fe_poly != nullptr && + (dynamic_cast *>( + &fe_poly->get_poly_space()) != nullptr || + dynamic_cast> *>( + &fe_poly->get_poly_space()) != nullptr)) + scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); + else if (const auto fe_dgp = dynamic_cast *>(fe)) + { + scalar_lexicographic.resize(fe_dgp->n_dofs_per_cell()); + for (unsigned int i = 0; i < fe_dgp->n_dofs_per_cell(); ++i) + scalar_lexicographic[i] = i; + element_type = truncated_tensor; + } + else if (const auto fe_q_dg0 = dynamic_cast *>(fe)) + { + scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); + element_type = tensor_symmetric_plus_dg0; + } + else if (fe->n_dofs_per_cell() == 0) + { + // FE_Nothing case -> nothing to do here + } + else + Assert(false, ExcNotImplemented()); + + // Finally store the renumbering into the member variable of this + // class + if (fe_in.n_components() == 1) + lexicographic_numbering = scalar_lexicographic; + else + { + // have more than one component, get the inverse + // permutation, invert it, sort the components one after one, + // and invert back + std::vector scalar_inv = + Utilities::invert_permutation(scalar_lexicographic); + std::vector lexicographic( + fe_in.n_dofs_per_cell(), numbers::invalid_unsigned_int); + unsigned int components_before = 0; + for (unsigned int e = 0; e < base_element_number; ++e) + components_before += fe_in.element_multiplicity(e); + for (unsigned int comp = 0; + comp < fe_in.element_multiplicity(base_element_number); + ++comp) + for (unsigned int i = 0; i < scalar_inv.size(); ++i) + lexicographic[fe_in.component_to_system_index( + comp + components_before, i)] = + scalar_inv.size() * comp + scalar_inv[i]; + + // invert numbering again. Need to do it manually because we might + // have undefined blocks + lexicographic_numbering.resize(fe_in.element_multiplicity( + base_element_number) * + fe->n_dofs_per_cell(), + numbers::invalid_unsigned_int); + for (unsigned int i = 0; i < lexicographic.size(); ++i) + if (lexicographic[i] != numbers::invalid_unsigned_int) + { + AssertIndexRange(lexicographic[i], + lexicographic_numbering.size()); + lexicographic_numbering[lexicographic[i]] = i; + } + } + } + + template ShapeInfo::ShapeInfo() @@ -150,8 +240,10 @@ namespace internal { #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT if (quad_in.is_tensor_product() == false || - dynamic_cast *>(&fe_in) || - dynamic_cast *>(&fe_in)) + dynamic_cast *>( + &fe_in.base_element(base_element_number)) || + dynamic_cast *>( + &fe_in.base_element(base_element_number))) { // specialization for arbitrary finite elements and quadrature rules // as needed in the context, e.g., of simplices @@ -282,17 +374,19 @@ namespace internal // shape_gradients_collocation_eo, shape_hessians_collocation_eo, // inverse_shape_values_eo cannot be filled - // indicate that no tensor product properties could be exploited - element_type = tensor_none; + std::vector scalar_lexicographic; + internal::MatrixFreeFunctions::get_element_type_specific_information( + fe_in, + fe, + base_element_number, + element_type, + scalar_lexicographic, + lexicographic_numbering); + univariate_shape_data.element_type = this->element_type; univariate_shape_data.nodal_at_cell_boundaries = true; - lexicographic_numbering.resize(n_dofs); - - for (unsigned int i = 0; i < n_dofs; ++i) - lexicographic_numbering[i] = i; - // TODO: setup face_to_cell_index_nodal, face_to_cell_index_hermite, // face_orientations @@ -358,79 +452,13 @@ namespace internal Assert(fe->n_components() == 1, ExcMessage("Expected a scalar element")); - const FE_Poly *fe_poly = - dynamic_cast *>(fe); - - const FE_DGP *fe_dgp = dynamic_cast *>(fe); - - const FE_Q_DG0 *fe_q_dg0 = dynamic_cast *>(fe); - - element_type = tensor_general; - if (fe_poly != nullptr && - (dynamic_cast *>( - &fe_poly->get_poly_space()) != nullptr || - dynamic_cast> *>( - &fe_poly->get_poly_space()) != nullptr)) - scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); - else if (fe_dgp != nullptr) - { - scalar_lexicographic.resize(fe_dgp->n_dofs_per_cell()); - for (unsigned int i = 0; i < fe_dgp->n_dofs_per_cell(); ++i) - scalar_lexicographic[i] = i; - element_type = truncated_tensor; - } - else if (fe_q_dg0 != nullptr) - { - scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); - element_type = tensor_symmetric_plus_dg0; - } - else if (fe->n_dofs_per_cell() == 0) - { - // FE_Nothing case -> nothing to do here - } - else - Assert(false, ExcNotImplemented()); - - // Finally store the renumbering into the member variable of this - // class - if (fe_in.n_components() == 1) - lexicographic_numbering = scalar_lexicographic; - else - { - // have more than one component, get the inverse - // permutation, invert it, sort the components one after one, - // and invert back - std::vector scalar_inv = - Utilities::invert_permutation(scalar_lexicographic); - std::vector lexicographic( - fe_in.n_dofs_per_cell(), numbers::invalid_unsigned_int); - unsigned int components_before = 0; - for (unsigned int e = 0; e < base_element_number; ++e) - components_before += fe_in.element_multiplicity(e); - for (unsigned int comp = 0; - comp < fe_in.element_multiplicity(base_element_number); - ++comp) - for (unsigned int i = 0; i < scalar_inv.size(); ++i) - lexicographic[fe_in.component_to_system_index( - comp + components_before, i)] = - scalar_inv.size() * comp + scalar_inv[i]; - - // invert numbering again. Need to do it manually because we might - // have undefined blocks - lexicographic_numbering.resize(fe_in.element_multiplicity( - base_element_number) * - fe->n_dofs_per_cell(), - numbers::invalid_unsigned_int); - for (unsigned int i = 0; i < lexicographic.size(); ++i) - if (lexicographic[i] != numbers::invalid_unsigned_int) - { - AssertIndexRange(lexicographic[i], - lexicographic_numbering.size()); - lexicographic_numbering[lexicographic[i]] = i; - } - } + internal::MatrixFreeFunctions::get_element_type_specific_information( + fe_in, + fe, + base_element_number, + element_type, + scalar_lexicographic, + lexicographic_numbering); // to evaluate 1D polynomials, evaluate along the line with the first // unit support point, assuming that fe.shape_value(0,unit_point) == -- 2.39.5