From 33cdee0fe32c33ac3aef0ac8676a3af57020ffbd Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Wed, 13 Jan 2021 22:40:41 -0700 Subject: [PATCH] simplex: line_dof_identities for FE_P --- source/fe/fe_q_base.cc | 37 ++++++++++++++++-- source/simplex/fe_lib.cc | 81 +++++++++++++++++++++++++++++++++++----- 2 files changed, 106 insertions(+), 12 deletions(-) diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index a869d95a6b..890ef9ed13 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -813,13 +813,44 @@ FE_Q_Base::hp_line_dof_identities( [0]) < 1e-14) identities.emplace_back(i, j); + return identities; + } + else if (const Simplex::FE_P *fe_p_other = + dynamic_cast *>(&fe_other)) + { + // DoFs are located along lines, so two dofs are identical if they are + // located at identical positions. If we had only equidistant points, we + // could simply check for similarity like (i+1)*q == (j+1)*p, but we + // might have other support points (e.g. Gauss-Lobatto + // points). Therefore, read the points in unit_support_points for the + // first coordinate direction. For FE_Q, we take the lexicographic + // ordering of the line support points in the first direction (i.e., + // x-direction), which we access between index 1 and p-1 (index 0 and p + // are vertex dofs). For FE_P, they are currently hard-coded and we + // iterate over points on the first line which begin after the 3 vertex + // points in the complete list of unit support points + + Assert(fe_p_other->degree <= 2, ExcNotImplemented()); + + const std::vector &index_map_inverse_q = + this->get_poly_space_numbering_inverse(); + + std::vector> identities; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j) + if (std::fabs( + this->unit_support_points[index_map_inverse_q[i + 1]][0] - + fe_p_other->get_unit_support_points()[j + 3][0]) < 1e-14) + identities.emplace_back(i, j); + return identities; } else if (dynamic_cast *>(&fe_other) != nullptr) { // the FE_Nothing has no degrees of freedom, so there are no // equivalencies to be recorded - return std::vector>(); + return {}; } else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) { @@ -830,12 +861,12 @@ FE_Q_Base::hp_line_dof_identities( // 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 std::vector>(); + return {}; } else { Assert(false, ExcNotImplemented()); - return std::vector>(); + return {}; } } diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index f566795c62..23319c4536 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -498,19 +498,82 @@ namespace Simplex FE_P::hp_line_dof_identities( const FiniteElement &fe_other) const { - (void)fe_other; - - Assert((dynamic_cast *>(&fe_other)), - ExcNotImplemented()); AssertDimension(dim, 2); - AssertDimension(this->degree, fe_other.tensor_degree()); + Assert(this->degree <= 2, ExcNotImplemented()); - std::vector> result; + if (const FE_P *fe_p_other = + dynamic_cast *>(&fe_other)) + { + // dofs are located along lines, so two dofs are identical if they are + // located at identical positions. + // Therefore, read the points in unit_support_points for the + // first coordinate direction. For FE_P, they are currently hard-coded + // and we iterate over points on the first line which begin after the 3 + // vertex points in the complete list of unit support points - for (unsigned int i = 0; i < this->degree - 1; ++i) - result.emplace_back(i, i); + Assert(fe_p_other->degree <= 2, ExcNotImplemented()); - return result; + std::vector> identities; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j) + if (std::fabs(this->unit_support_points[i + 3][0] - + fe_p_other->unit_support_points[i + 3][0]) < 1e-14) + identities.emplace_back(i, j); + + return identities; + } + else if (const FE_Q *fe_q_other = + dynamic_cast *>(&fe_other)) + { + // dofs are located along lines, so two dofs are identical if they are + // located at identical positions. if we had only equidistant points, we + // could simply check for similarity like (i+1)*q == (j+1)*p, but we + // might have other support points (e.g. Gauss-Lobatto + // points). Therefore, read the points in unit_support_points for the + // first coordinate direction. For FE_Q, we take the lexicographic + // ordering of the line support points in the first direction (i.e., + // x-direction), which we access between index 1 and p-1 (index 0 and p + // are vertex dofs). For FE_P, they are currently hard-coded and we + // iterate over points on the first line which begin after the 3 vertex + // points in the complete list of unit support points + + const std::vector &index_map_inverse_q_other = + fe_q_other->get_poly_space_numbering_inverse(); + + std::vector> identities; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j) + if (std::fabs(this->unit_support_points[i + 3][0] - + fe_q_other->get_unit_support_points() + [index_map_inverse_q_other[j + 1]][0]) < 1e-14) + identities.emplace_back(i, j); + + return identities; + } + 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 {}; + } } -- 2.39.5