From 9e78b76efd5c2fe67ac4936af4fa4a6b68159335 Mon Sep 17 00:00:00 2001
From: Leila Ghaffari <Leila.Ghaffari@colorado.edu>
Date: Tue, 25 Jun 2024 13:28:18 -0700
Subject: [PATCH] WIP: get fe from trace elements on the boundaries when there
 is some periodicity in the mesh

There are no dofs on the faces to constrain the desired boundaries
---
 palace/fem/libceed/restriction.cpp | 59 +++++++++++++++++++++++++-----
 1 file changed, 50 insertions(+), 9 deletions(-)

diff --git a/palace/fem/libceed/restriction.cpp b/palace/fem/libceed/restriction.cpp
index 130308397..12d5d655f 100644
--- a/palace/fem/libceed/restriction.cpp
+++ b/palace/fem/libceed/restriction.cpp
@@ -90,9 +90,26 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
                      Ceed ceed, CeedElemRestriction *restr)
 {
   const std::size_t num_elem = indices.size();
-  const mfem::FiniteElement &fe =
-      use_bdr ? *fespace.GetBE(indices[0]) : *fespace.GetFE(indices[0]);
-  const int P = fe.GetDof();
+  const mfem::FiniteElement *fe;
+  bool face_flg = false;
+  if (!use_bdr)
+  {
+    fe = fespace.GetFE(indices[0]);
+  }
+  else
+  {
+    fe = fespace.GetBE(indices[0]);
+    if (!fe)
+    {
+      int elem_id, face_info;
+      int FaceNo = fespace.GetMesh()->GetBdrElementFaceIndex(indices[0]);
+      fespace.GetMesh()->GetBdrElementAdjacentElement(FaceNo, elem_id, face_info);
+      mfem::Geometry::Type face_geom = fespace.GetMesh()->GetBdrElementGeometry(FaceNo);
+      fe = fespace.GetTraceElement(elem_id, face_geom);
+      face_flg = true;
+    }
+  }
+  const int P = fe->GetDof();
   const CeedInt comp_stride =
       (fespace.GetVDim() == 1 || fespace.GetOrdering() == mfem::Ordering::byVDIM)
           ? 1
@@ -105,7 +122,7 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
     {
       return false;
     }
-    const auto geom = fe.GetGeomType();
+    const auto geom = fe->GetGeomType();
     const auto *dof_trans = fespace.FEColl()->DofTransformationForGeometry(geom);
     return (dof_trans && !dof_trans->IsIdentity());
   }();
@@ -140,7 +157,15 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
       const auto e = indices[i];
       if (use_bdr)
       {
-        fespace.GetBdrElementDofs(e, dofs, dof_trans);
+        if (!face_flg)
+        {
+          fespace.GetBdrElementDofs(e, dofs, dof_trans);
+        }
+        else
+        {
+          int f = fespace.GetMesh()->GetBdrElementFaceIndex(e);
+          fespace.GetFaceVDofs(f, dofs);
+        }
       }
       else
       {
@@ -216,6 +241,7 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
           }
         }
       }
+
     }
     use_el_orients += use_el_orients_loc;
   }
@@ -258,10 +284,25 @@ void InitRestriction(const mfem::FiniteElementSpace &fespace,
               << indices[0] << ", " << use_bdr << ", " << is_interp << ", "
               << is_interp_range << ")\n";
   }
-  const mfem::FiniteElement &fe =
-      use_bdr ? *fespace.GetBE(indices[0]) : *fespace.GetFE(indices[0]);
-  const mfem::TensorBasisElement *tfe = dynamic_cast<const mfem::TensorBasisElement *>(&fe);
-  const bool vector = fe.GetRangeType() == mfem::FiniteElement::VECTOR;
+  const mfem::FiniteElement *fe;
+  if (!use_bdr)
+  {
+    fe = fespace.GetFE(indices[0]);
+  }
+  else
+  {
+    fe = fespace.GetBE(indices[0]);
+    if (!fe)
+    {
+      int elem_id, face_info;
+      int FaceNo = fespace.GetMesh()->GetBdrElementFaceIndex(indices[0]);
+      fespace.GetMesh()->GetBdrElementAdjacentElement(FaceNo, elem_id, face_info);
+      mfem::Geometry::Type face_geom = fespace.GetMesh()->GetBdrElementGeometry(FaceNo);
+      fe = fespace.GetTraceElement(elem_id, face_geom);
+    }
+  }
+  const mfem::TensorBasisElement *tfe = dynamic_cast<const mfem::TensorBasisElement *>(fe);
+  const bool vector = fe->GetRangeType() == mfem::FiniteElement::VECTOR;
   const bool lexico = (tfe && tfe->GetDofMap().Size() > 0 && !vector && !is_interp);
   if (lexico)
   {