From 0e509c3cc44af92e88a9aac0999a3abba4bf5616 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 21 May 2020 12:30:54 -0400 Subject: [PATCH] Fix numbering --- include/deal.II/fe/fe_poly.templates.h | 17 +++++- source/fe/fe_dgq.cc | 12 ++-- source/fe/fe_q_base.cc | 78 ++++++++++++++------------ 3 files changed, 63 insertions(+), 44 deletions(-) diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index fea8899e25..2f2e1776b7 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -579,13 +580,23 @@ FE_Poly::get_poly_space_numbering() const { auto *const space_tensor_prod = dynamic_cast *>(this->poly_space.get()); - if (space_tensor_prod != nullptr) return space_tensor_prod->get_numbering(); - auto *const space_tensor_prod_const = - dynamic_cast *>(this->poly_space.get()); + auto *const space_tensor_prod_piecewise = dynamic_cast< + TensorProductPolynomials> *>( + this->poly_space.get()); + if (space_tensor_prod_piecewise != nullptr) + return space_tensor_prod_piecewise->get_numbering(); + auto *const space_tensor_prod_bubbles = + dynamic_cast *>( + this->poly_space.get()); + if (space_tensor_prod_bubbles != nullptr) + return space_tensor_prod_bubbles->get_numbering(); + + auto *const space_tensor_prod_const = + dynamic_cast *>(this->poly_space.get()); if (space_tensor_prod_const != nullptr) return space_tensor_prod_const->get_numbering(); diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index f54154595d..ad26b3049d 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -835,11 +835,11 @@ FE_DGQArbitraryNodes::get_name() const bool equidistant = true; std::vector points(this->degree + 1); - auto *const polynomial_space_p = + auto *const polynomial_space = dynamic_cast *>(this->poly_space.get()); - Assert(polynomial_space_p != nullptr, ExcInternalError()); + Assert(polynomial_space != nullptr, ExcInternalError()); std::vector lexicographic = - polynomial_space_p->get_numbering_inverse(); + polynomial_space->get_numbering_inverse(); for (unsigned int j = 0; j <= this->degree; j++) points[j] = this->unit_support_points[lexicographic[j]][0]; @@ -952,11 +952,11 @@ FE_DGQArbitraryNodes::clone() const { // Construct a dummy quadrature formula containing the FE's nodes: std::vector> qpoints(this->degree + 1); - auto *const polynomial_space_p = + auto *const polynomial_space = dynamic_cast *>(this->poly_space.get()); - Assert(polynomial_space_p != nullptr, ExcInternalError()); + Assert(polynomial_space != nullptr, ExcInternalError()); std::vector lexicographic = - polynomial_space_p->get_numbering_inverse(); + polynomial_space->get_numbering_inverse(); for (unsigned int i = 0; i <= this->degree; ++i) qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]); Quadrature<1> pquadrature(qpoints); diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index dc183f5151..1dff72f2b4 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -180,10 +180,8 @@ struct FE_Q_Base::Implementation // use that the element evaluates to 1 at index 0 and along the line at // zero - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(fe.poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); + fe.get_poly_space_numbering_inverse(); const std::vector face_index_map = FETools::lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs( @@ -317,10 +315,8 @@ struct FE_Q_Base::Implementation // use that the element evaluates to 1 at index 0 and along the line at // zero - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(fe.poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); + fe.get_poly_space_numbering_inverse(); const std::vector face_index_map = FETools::lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs( @@ -444,15 +440,44 @@ FE_Q_Base::initialize( q_dofs_per_cell + dim == this->dofs_per_cell, ExcInternalError()); - { + [this, q_dofs_per_cell]() { std::vector renumber = FETools::hierarchic_to_lexicographic_numbering(q_degree); for (unsigned int i = q_dofs_per_cell; i < this->dofs_per_cell; ++i) renumber.push_back(i); - TensorProductPolynomials *poly_space_derived_ptr = + auto *tensor_poly_space_ptr = dynamic_cast *>(this->poly_space.get()); - poly_space_derived_ptr->set_numbering(renumber); - } + if (tensor_poly_space_ptr != nullptr) + { + tensor_poly_space_ptr->set_numbering(renumber); + return; + } + auto *tensor_piecewise_poly_space_ptr = dynamic_cast< + TensorProductPolynomials> + *>(this->poly_space.get()); + if (tensor_piecewise_poly_space_ptr != nullptr) + { + tensor_piecewise_poly_space_ptr->set_numbering(renumber); + return; + } + auto *tensor_bubbles_poly_space_ptr = + dynamic_cast *>( + this->poly_space.get()); + if (tensor_bubbles_poly_space_ptr != nullptr) + { + tensor_bubbles_poly_space_ptr->set_numbering(renumber); + return; + } + auto *tensor_const_poly_space_ptr = + dynamic_cast *>( + this->poly_space.get()); + if (tensor_const_poly_space_ptr != nullptr) + { + tensor_const_poly_space_ptr->set_numbering(renumber); + return; + } + Assert(false, ExcNotImplemented()); + }(); // Finally fill in support points on cell and face and initialize // constraints. All of this can happen in parallel @@ -754,15 +779,10 @@ FE_Q_Base::hp_line_dof_identities( std::vector> identities; - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); - TensorProductPolynomials *poly_space_derived_ptr_other = - dynamic_cast *>( - fe_q_other->poly_space.get()); + this->get_poly_space_numbering_inverse(); const std::vector &index_map_inverse_other = - poly_space_derived_ptr_other->get_numbering_inverse(); + fe_q_other->get_poly_space_numbering_inverse(); for (unsigned int i = 0; i < p - 1; ++i) for (unsigned int j = 0; j < q - 1; ++j) @@ -821,15 +841,10 @@ FE_Q_Base::hp_quad_dof_identities( std::vector> identities; - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); - TensorProductPolynomials *poly_space_derived_ptr_other = - dynamic_cast *>( - fe_q_other->poly_space.get()); + this->get_poly_space_numbering_inverse(); const std::vector &index_map_inverse_other = - poly_space_derived_ptr_other->get_numbering_inverse(); + fe_q_other->get_poly_space_numbering_inverse(); for (unsigned int i1 = 0; i1 < p - 1; ++i1) for (unsigned int i2 = 0; i2 < p - 1; ++i2) @@ -886,10 +901,8 @@ void FE_Q_Base::initialize_unit_support_points( const std::vector> &points) { - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); + this->get_poly_space_numbering_inverse(); // We can compute the support points by computing the tensor // product of the 1d set of points. We could do this by hand, but it's @@ -1223,10 +1236,8 @@ FE_Q_Base::get_prolongation_matrix( std::vector> subcell_evaluations( dim, Table<2, double>(dofs1d, dofs1d)); - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); + this->get_poly_space_numbering_inverse(); // helper value: step size how to walk through diagonal and how many // points we have left apart from the first dimension @@ -1400,12 +1411,9 @@ FE_Q_Base::get_restriction_matrix( // assumption that whenever a row makes a non-zero contribution to the // mother's residual, the correct value is interpolated. - const double eps = 1e-15 * q_degree * dim; - - TensorProductPolynomials *poly_space_derived_ptr = - dynamic_cast *>(this->poly_space.get()); + const double eps = 1e-15 * q_degree * dim; const std::vector &index_map_inverse = - poly_space_derived_ptr->get_numbering_inverse(); + this->get_poly_space_numbering_inverse(); const unsigned int dofs1d = q_degree + 1; std::vector> evaluations1d(dofs1d); -- 2.39.5