From 679e5554d83972ce10fe8358d67ba7c49d041351 Mon Sep 17 00:00:00 2001
From: Leila Ghaffari <Leila.Ghaffari@colorado.edu>
Date: Wed, 17 Jul 2024 14:47:35 -0700
Subject: [PATCH] Get face and element point matrices

---
 palace/fem/libceed/restriction.cpp | 56 +++++++++++++++++++++++++++---
 1 file changed, 52 insertions(+), 4 deletions(-)

diff --git a/palace/fem/libceed/restriction.cpp b/palace/fem/libceed/restriction.cpp
index 12d5d655f..bb9324df7 100644
--- a/palace/fem/libceed/restriction.cpp
+++ b/palace/fem/libceed/restriction.cpp
@@ -102,8 +102,8 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
     if (!fe)
     {
       int elem_id, face_info;
-      int FaceNo = fespace.GetMesh()->GetBdrElementFaceIndex(indices[0]);
-      fespace.GetMesh()->GetBdrElementAdjacentElement(FaceNo, elem_id, face_info);
+      int FaceNo = fespace.GetMesh()->GetBdrElementFaceIndex(indices[0]); // TODO: I think we don't this line
+      fespace.GetMesh()->GetBdrElementAdjacentElement(FaceNo, elem_id, face_info); // I think the first argument should be the boundary element index
       mfem::Geometry::Type face_geom = fespace.GetMesh()->GetBdrElementGeometry(FaceNo);
       fe = fespace.GetTraceElement(elem_id, face_geom);
       face_flg = true;
@@ -163,8 +163,56 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
         }
         else
         {
-          int f = fespace.GetMesh()->GetBdrElementFaceIndex(e);
-          fespace.GetFaceVDofs(f, dofs);
+          // Get coordinates of face dofs
+          int elem_id, face_info;
+          fespace.GetMesh()->GetBdrElementAdjacentElement(e, elem_id, face_info);
+          mfem::Geometry::Type face_geom = fespace.GetMesh()->GetBdrElementGeometry(e);
+          face_info = fespace.GetMesh()->EncodeFaceInfo(
+                         fespace.GetMesh()->DecodeFaceInfoLocalIndex(face_info),
+                         mfem::Geometry::GetInverseOrientation(
+                            face_geom, fespace.GetMesh()->DecodeFaceInfoOrientation(face_info))
+                      );
+          mfem::IntegrationPointTransformation Loc1;
+          fespace.GetMesh()->GetLocalFaceTransformation(fespace.GetMesh()->GetBdrElementType(e),
+                                     fespace.GetMesh()->GetElementType(elem_id),
+                                     Loc1.Transf, face_info);
+          const mfem::FiniteElement *face_el = fespace.GetTraceElement(elem_id, face_geom);
+          MFEM_VERIFY(dynamic_cast<const mfem::NodalFiniteElement*>(face_el),
+                      "Mesh requires nodal Finite Element.");
+          mfem::IntegrationRule face_ir(face_el->GetDof());
+          Loc1.Transf.ElementNo = elem_id;
+          Loc1.Transf.mesh = fespace.GetMesh();
+          Loc1.Transf.ElementType = mfem::ElementTransformation::ELEMENT;
+          Loc1.Transform(face_el->GetNodes(), face_ir);
+          mfem::IsoparametricTransformation face_tr;
+          face_tr.ElementNo = e;
+          face_tr.ElementType = mfem::ElementTransformation::BDR_ELEMENT;
+          face_tr.mesh = fespace.GetMesh();
+          face_tr.Attribute = fespace.GetMesh()->GetBdrAttribute(e);
+          mfem::DenseMatrix &face_pm = face_tr.GetPointMat(); // dim x dof
+          face_tr.Reset();
+          fespace.GetMesh()->GetNodes()->GetVectorValues(Loc1.Transf, face_ir, face_pm);
+
+          std::cout << std::fixed << std::setprecision(2);
+          for (int m=0; m<face_pm.Height(); m++) {
+            for (int n=0; n<face_pm.Width(); n++) std::cout << "face[" <<m<<","<<n<<"] = " << face_pm(m, n) << "\t";
+            std::cout << std::endl;
+          }
+
+          // Get coordinates of element dofs
+          mfem::DenseMatrix elem_pm;
+          const mfem::FiniteElement *fe_elem = fespace.GetFE(elem_id);
+          const mfem::IntegrationRule &fe_elem_nodes = fe_elem->GetNodes();
+          mfem::ElementTransformation *T = fespace.GetMesh()->GetElementTransformation(elem_id);
+          T->Transform(fe_elem_nodes, elem_pm);
+
+          for (int m=0; m<elem_pm.Height(); m++) {
+            for (int n=0; n<elem_pm.Width(); n++) std::cout << "elem[" <<m<<","<<n<<"] = " << elem_pm(m, n) << "\t";
+            std::cout << std::endl;
+          }
+
+          // Find the dofs
+
         }
       }
       else