From fa0ea629a76d012918b630221253fb7773bd1913 Mon Sep 17 00:00:00 2001 From: Ivy Weber Date: Mon, 24 Apr 2023 11:17:34 +0200 Subject: [PATCH] Added code for non-matching DOF identities between Hermite and FE_nothing --- include/deal.II/fe/fe_hermite.h | 32 +++- source/fe/fe_hermite.cc | 254 +++++++++++++++++++++++++++++++- source/fe/fe_q.cc | 16 ++ source/fe/fe_q_base.cc | 9 ++ 4 files changed, 307 insertions(+), 4 deletions(-) diff --git a/include/deal.II/fe/fe_hermite.h b/include/deal.II/fe/fe_hermite.h index 628b9e4171..eaa739c8a5 100644 --- a/include/deal.II/fe/fe_hermite.h +++ b/include/deal.II/fe/fe_hermite.h @@ -119,6 +119,36 @@ public: virtual std::unique_ptr> clone() const override; + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ +virtual std::vector> +hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ +virtual std::vector> +hp_line_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_quad_dof_identities() + */ +virtual std::vector> +hp_quad_dof_identities( + const FiniteElement &fe_other, + const unsigned int face_no = 0) const override; + + /** + * @copydoc FiniteElement::compare_for_domination() + */ + virtual FiniteElementDomination::Domination + compare_for_domination( + const FiniteElement& other_fe, + const unsigned int codim) const override; + /** * Returns the mapping between lexicographic and hierarchic numbering * schemes for Hermite. See the class documentation for diagrams of @@ -179,8 +209,6 @@ public: spacedim> &output_data) const override; - - protected: /** * Wrapper function for FE_Poly::get_data() that ensures relevant diff --git a/source/fe/fe_hermite.cc b/source/fe/fe_hermite.cc index 328eb51107..858783dd0a 100644 --- a/source/fe/fe_hermite.cc +++ b/source/fe/fe_hermite.cc @@ -21,6 +21,13 @@ #include #include +#include +#include +#include +#include +#include +#include +#include #include #include #include @@ -560,6 +567,248 @@ FE_Hermite::clone() const +/** + * A large part of the following function is copied from FE_Q_Base, the main + * difference is that some additional constraints are needed between two + * Hermite finite elements + */ +template +std::vector> +FE_Hermite::hp_vertex_dof_identities( + const FiniteElement &fe_other) const +{ + /*if (const FE_Hermite *fe_herm_other = + dynamic_cast *>(&fe_other)) + { + // In this case there will be several DoFs that need to be matched between + // the elements to ensure continuity. The number of DoFs to be matched is + // dependent on the polynomial degree of the lower order element, and dim + // + // Note: is this using hierarchical or lexicographic numbering at the vertices? + if (this->degree < fe_herm_other->degree) + { + std::vector> dof_matches; + dof_matches.reserve(this->n_dofs_per_vertex()); + + + } + else + { + + } + } + else */if (dynamic_cast *>(&fe_other) != nullptr) + { + // there should be exactly one single DoF of FE_Q_Base at a vertex, and it + // should have an identical value to the first Hermite DoF + return {{0U, 0U}}; + } + else if (dynamic_cast *>(&fe_other) != + nullptr) + { + // there should be exactly one single DoF of FE_Q_Base at a vertex, and it + // should have an identical value to the first Hermite DoF + return {{0U, 0U}}; + } + else if (dynamic_cast *>(&fe_other) != nullptr) + { + // the FE_Nothing has no degrees of freedom, so there are no + // equivalencies to be recorded + return {}; + } + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) + { + // if the other element has no elements on faces at all, + // then it would be impossible to enforce any kind of + // continuity even if we knew exactly what kind of element + // we have -- simply because the other element declares + // that it is discontinuous because it has no DoFs on + // its faces. in that case, just state that we have no + // constraints to declare + return {}; + } + else + { + Assert(false, ExcNotImplemented()); + return {}; + } +} + + + +template +std::vector> +FE_Hermite::hp_line_dof_identities( + const FiniteElement &fe_other) const +{ + if (dynamic_cast *>(&fe_other) != nullptr) + { + // the FE_Nothing has no degrees of freedom, so there are no + // equivalencies to be recorded + return {}; + } + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) + { + // if the other element has no elements on faces at all, + // then it would be impossible to enforce any kind of + // continuity even if we knew exactly what kind of element + // we have -- simply because the other element declares + // that it is discontinuous because it has no DoFs on + // its faces. in that case, just state that we have no + // constraints to declare + return {}; + } + else + { + Assert(false, ExcNotImplemented()); + return {}; + } +} + + + +template +std::vector> +FE_Hermite::hp_quad_dof_identities( + const FiniteElement &fe_other, + const unsigned int face_no) const +{ + if (dynamic_cast *>(&fe_other) != nullptr) + { + // the FE_Nothing has no degrees of freedom, so there are no + // equivalencies to be recorded + return {}; + } + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) + { + // if the other element has no elements on faces at all, + // then it would be impossible to enforce any kind of + // continuity even if we knew exactly what kind of element + // we have -- simply because the other element declares + // that it is discontinuous because it has no DoFs on + // its faces. in that case, just state that we have no + // constraints to declare + return {}; + } + else + { + Assert(false, ExcNotImplemented()); + return {}; + } +} + + + +/* +* The layout of this function is largely copied directly from FE_Q, +* however FE_Hermite can behave significantly differently in terms +* of domination due to how the function space is defined */ +template +FiniteElementDomination::Domination +FE_Hermite::compare_for_domination( + const FiniteElement& fe_other, + const unsigned int codim) const + { + Assert(codim <= dim, ExcImpossibleInDim(dim)); + + if (codim > 0) + if (dynamic_cast *>(&fe_other) != nullptr) + // there are no requirements between continuous and discontinuous elements + return FiniteElementDomination::no_requirements; + + + // vertex/line/face domination + // (if fe_other is not derived from FE_DGQ) + // & cell domination + // ---------------------------------------- + if (const FE_Hermite *fe_hermite_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_hermite_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_hermite_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + if (const FE_Q *fe_q_other = + dynamic_cast *>(&fe_other)) + { + if (fe_q_other->degree == 1) + { + if (this->degree == 1) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (this->degree <= fe_q_other->degree) + return FiniteElementDomination::this_element_dominates; + else + return FiniteElementDomination::neither_element_dominates; + } + else if (const FE_SimplexP *fe_p_other = + dynamic_cast *>(&fe_other)) + { + if (fe_p_other->degree == 1) + { + if (this->degree == 1) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (this->degree <= fe_p_other->degree) + return FiniteElementDomination::this_element_dominates; + else + return FiniteElementDomination::neither_element_dominates; + } + else if (const FE_WedgeP *fe_wp_other = + dynamic_cast *>(&fe_other)) + { + if (fe_wp_other->degree == 1) + { + if (this->degree == 1) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (this->degree <= fe_wp_other->degree) + return FiniteElementDomination::this_element_dominates; + else + return FiniteElementDomination::neither_element_dominates; + } + else if (const FE_PyramidP *fe_pp_other = + dynamic_cast *>(&fe_other)) + { + if (fe_pp_other->degree == 1) + { + if (this->degree == 1) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (this->degree <= fe_pp_other->degree) + return FiniteElementDomination::this_element_dominates; + else + return FiniteElementDomination::neither_element_dominates; + } + else if (const FE_Nothing *fe_nothing = + dynamic_cast *>(&fe_other)) + { + if (fe_nothing->is_dominating()) + return FiniteElementDomination::other_element_dominates; + else + // the FE_Nothing has no degrees of freedom and it is typically used + // in a context where we don't require any continuity along the + // interface + return FiniteElementDomination::no_requirements; + } + + Assert(false, ExcNotImplemented()); + return FiniteElementDomination::neither_element_dominates; + } + + + template std::vector FE_Hermite::get_lexicographic_to_hierarchic_numbering() const @@ -597,11 +846,12 @@ FE_Hermite::fill_fe_values( fe_internal); Assert( - (dynamic_cast::InternalData *>( + ((dynamic_cast::InternalData *>( &mapping_internal) != nullptr) || (dynamic_cast::InternalData *>( - &mapping_internal) != nullptr), + &mapping_internal) != nullptr)), ExcInternalError()); + const UpdateFlags flags(fe_data.update_each); diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index 47d5de838b..2fd438d63f 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include @@ -251,6 +252,21 @@ FE_Q::compare_for_domination( // interface return FiniteElementDomination::no_requirements; } + else if (const FE_Hermite *fe_hermite_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree == 1) + { + if (fe_hermite_other->degree > 1) + return FiniteElementDomination::this_element_dominates; + else + return FiniteElementDomination::either_element_can_dominate; + } + else if (this->degree >= fe_hermite_other->degree) + return FiniteElementDomination::other_element_dominates; + else + return FiniteElementDomination::neither_element_dominates; + } Assert(false, ExcNotImplemented()); return FiniteElementDomination::neither_element_dominates; diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 281f698e41..af66f99479 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -752,6 +753,14 @@ FE_Q_Base::hp_vertex_dof_identities( // should have identical value return {{0U, 0U}}; } + else if (dynamic_cast *>(&fe_other) != nullptr) + { + // FE_Hermite will usually have several degrees of freedom on + // each vertex, however only the first one will actually + // correspond to the shape value at the vertex, meaning it's + // the only one of interest for FE_Q_Base + return {{0U, 0U}}; + } else if (dynamic_cast *>(&fe_other) != nullptr) { // the FE_Nothing has no degrees of freedom, so there are no -- 2.39.5