simplex: line_dof_identities for FE_P

diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index a869d95a6bce40787e8ecb34a05cec980ca5072d..890ef9ed136ddfc8a66dff863b29cbc730332969 100644 (file)
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -813,13 +813,44 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_line_dof_identities(
                                                [0]) < 1e-14)
             identities.emplace_back(i, j);
 
+      return identities;
+    }
+  else if (const Simplex::FE_P<dim, spacedim> *fe_p_other =
+             dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&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<unsigned int> &index_map_inverse_q =
+        this->get_poly_space_numbering_inverse();
+
+      std::vector<std::pair<unsigned int, unsigned int>> 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<const FE_Nothing<dim> *>(&fe_other) != nullptr)
     {
       // the FE_Nothing has no degrees of freedom, so there are no
       // equivalencies to be recorded
-      return std::vector<std::pair<unsigned int, unsigned int>>();
+      return {};
     }
   else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
     {
@@ -830,12 +861,12 @@ FE_Q_Base<PolynomialType, dim, spacedim>::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<std::pair<unsigned int, unsigned int>>();
+      return {};
     }
   else
     {
       Assert(false, ExcNotImplemented());
-      return std::vector<std::pair<unsigned int, unsigned int>>();
+      return {};
     }
 }
 

diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc
index f566795c62b383373cd21fdf3bfa745e038d64dc..23319c4536486a8f84cf4c899777819df8759851 100644 (file)
--- a/source/simplex/fe_lib.cc
+++ b/source/simplex/fe_lib.cc
@@ -498,19 +498,82 @@ namespace Simplex
   FE_P<dim, spacedim>::hp_line_dof_identities(
     const FiniteElement<dim, spacedim> &fe_other) const
   {
-    (void)fe_other;
-
-    Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
-           ExcNotImplemented());
     AssertDimension(dim, 2);
-    AssertDimension(this->degree, fe_other.tensor_degree());
+    Assert(this->degree <= 2, ExcNotImplemented());
 
-    std::vector<std::pair<unsigned int, unsigned int>> result;
+    if (const FE_P<dim, spacedim> *fe_p_other =
+          dynamic_cast<const FE_P<dim, spacedim> *>(&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<std::pair<unsigned int, unsigned int>> 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<dim, spacedim> *fe_q_other =
+               dynamic_cast<const FE_Q<dim, spacedim> *>(&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<unsigned int> &index_map_inverse_q_other =
+          fe_q_other->get_poly_space_numbering_inverse();
+
+        std::vector<std::pair<unsigned int, unsigned int>> 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<const FE_Nothing<dim> *>(&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 {};
+      }
   }