From 9a2329651ec2a1c629a8b8a011e4dc1dea552c24 Mon Sep 17 00:00:00 2001
From: Marc Fehling <marc.fehling@gmx.net>
Date: Wed, 13 Jan 2021 22:19:32 -0700
Subject: [PATCH] simplex: vertex_dof_identities for FE_P

---
 source/fe/fe_q_base.cc   | 24 ++++++++++++++---------
 source/simplex/fe_lib.cc | 41 ++++++++++++++++++++++++++++++++++------
 2 files changed, 50 insertions(+), 15 deletions(-)

diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index 9f444f7aa0..a869d95a6b 100644
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -29,6 +29,8 @@
 #include <deal.II/fe/fe_q_base.h>
 #include <deal.II/fe/fe_tools.h>
 
+#include <deal.II/simplex/fe_lib.h>
+
 #include <memory>
 #include <sstream>
 #include <vector>
@@ -734,21 +736,25 @@ std::vector<std::pair<unsigned int, unsigned int>>
 FE_Q_Base<PolynomialType, dim, spacedim>::hp_vertex_dof_identities(
   const FiniteElement<dim, spacedim> &fe_other) const
 {
-  // we can presently only compute these identities if both FEs are FE_Qs or
-  // if the other one is an FE_Nothing. in the first case, there should be
-  // exactly one single DoF of each FE at a vertex, and they should have
-  // identical value
   if (dynamic_cast<const FE_Q_Base<PolynomialType, dim, spacedim> *>(
         &fe_other) != nullptr)
     {
-      return std::vector<std::pair<unsigned int, unsigned int>>(
-        1, std::make_pair(0U, 0U));
+      // there should be exactly one single DoF of each FE at a vertex, and they
+      // should have identical value
+      return {{0U, 0U}};
+    }
+  else if (dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other) !=
+           nullptr)
+    {
+      // there should be exactly one single DoF of each FE at a vertex, and they
+      // should have identical value
+      return {{0U, 0U}};
     }
   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)
     {
@@ -759,12 +765,12 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_vertex_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 ca3141b1cd..038b76a292 100644
--- a/source/simplex/fe_lib.cc
+++ b/source/simplex/fe_lib.cc
@@ -16,6 +16,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_tools.h>
 
@@ -423,14 +424,42 @@ namespace Simplex
   FE_P<dim, spacedim>::hp_vertex_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());
 
-    return {{0, 0}};
+    if (dynamic_cast<const FE_P<dim, spacedim> *>(&fe_other) != nullptr)
+      {
+        // there should be exactly one single DoF of each FE at a vertex, and
+        // they should have identical value
+        return {{0U, 0U}};
+      }
+    else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
+      {
+        // there should be exactly one single DoF of each FE at a vertex, and
+        // they should have identical value
+        return {{0U, 0U}};
+      }
+    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 {};
+      }
   }
 
 
-- 
2.39.5