From ca9c1e5e44b8db51de00e330626196397200aa9a Mon Sep 17 00:00:00 2001
From: Leila Ghaffari <Leila.Ghaffari@colorado.edu>
Date: Mon, 12 Aug 2024 17:04:05 -0700
Subject: [PATCH] WIP - Removed all boundary elements that have two elments on
 both sides (not considering interior BEs).

---
 palace/fem/libceed/restriction.cpp |  6 +--
 palace/utils/geodata.cpp           | 72 ++++++++++++++++++++----------
 2 files changed, 50 insertions(+), 28 deletions(-)

diff --git a/palace/fem/libceed/restriction.cpp b/palace/fem/libceed/restriction.cpp
index c014e4464..6eaee312d 100644
--- a/palace/fem/libceed/restriction.cpp
+++ b/palace/fem/libceed/restriction.cpp
@@ -202,7 +202,7 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
           // Find the dofs
           mfem::real_t tol = 1E-5;
           mfem::Array<int> elem_dofs;
-          fespace.GetElementDofs(elem_id, elem_dofs, dof_trans); // TODO: Check if passing dof_trans is OK
+          fespace.GetElementDofs(elem_id, elem_dofs, dof_trans);
           dofs.SetSize(P);
           for (int l=0; l< P; l++)
           {
@@ -226,10 +226,6 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace,
                 {
                   dofs[l] = elem_dofs[m];
                 }
-                else
-                {
-                  continue;
-                }
               }
               else
               {
diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp
index 7af80e8c0..980dd1b21 100644
--- a/palace/utils/geodata.cpp
+++ b/palace/utils/geodata.cpp
@@ -2175,14 +2175,31 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
       bdr_elem_delete(orig_mesh->GetNBE(), false);
   std::unordered_map<int, int> face_to_be, new_face_be;
   face_to_be.reserve(orig_mesh->GetNBE());
-  //for (int be = 0; be < orig_mesh->GetNBE(); be++) // TODO: need to fix this!
-  //{
-  //  int f, o;
-  //  orig_mesh->GetBdrElementFace(be, &f, &o);
-  //  MFEM_VERIFY(face_to_be.find(f) == face_to_be.end(),
-  //              "Mesh should not define boundary elements multiple times!");
-  //  face_to_be[f] = be;
-  //}
+  for (int be = 0; be < orig_mesh->GetNBE(); be++)
+  {
+    int attr = orig_mesh->GetBdrAttribute(be);
+    bool periodic_be = false;
+    for (auto &data : iodata.boundaries.periodic)
+    {
+      for (auto d_attr : data.donor_attributes) if (d_attr == attr) periodic_be = true;
+      for (auto r_attr : data.receiver_attributes) if (r_attr == attr) periodic_be = true;
+    }
+    int face, o, e1, e2;
+    orig_mesh->GetBdrElementFace(be, &face, &o);
+    orig_mesh->GetFaceElements(face, &e1, &e2);
+    // If there are two elements associated with the face on the boundary, then the mesh is periodic.
+    if ((e1 >= 0) && (e2 >= 0))
+    {
+      periodic_be = true;
+    }
+    if (!periodic_be) {
+      int f, o;
+      orig_mesh->GetBdrElementFace(be, &f, &o);
+      MFEM_VERIFY(face_to_be.find(f) == face_to_be.end(),
+                  "Mesh should not define boundary elements multiple times!");
+      face_to_be[f] = be;
+    }
+  }
 
   if (clean_elem)
   {
@@ -2215,23 +2232,32 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
       }
     }
     // Remove any boundary element that is on the periodic surfaces.
-    if (!iodata.boundaries.periodic.empty()) {
-      for (int be = 0; be < orig_mesh->GetNBE(); be++)
+    for (int be = 0; be < orig_mesh->GetNBE(); be++)
+    {
+      bool periodic_be = false;
+      int attr = orig_mesh->GetBdrAttribute(be);
+      for (auto &data : iodata.boundaries.periodic)
       {
-        int attr = orig_mesh->GetBdrAttribute(be);
-        for (auto &data : iodata.boundaries.periodic)
-        {
-          bool periodic_be = false;
-          for (auto d_attr : data.donor_attributes) if (d_attr == attr) periodic_be = true;
-          for (auto r_attr : data.receiver_attributes) if (r_attr == attr) periodic_be = true;
-          if (periodic_be)
-          {
-            bdr_elem_delete[be] = true;
-            new_nbe--;
-          }
-        }
+        for (auto d_attr : data.donor_attributes) if (d_attr == attr) periodic_be = true;
+        for (auto r_attr : data.receiver_attributes) if (r_attr == attr) periodic_be = true;
+      }
+      int face, o, e1, e2;
+      orig_mesh->GetBdrElementFace(be, &face, &o);
+      orig_mesh->GetFaceElements(face, &e1, &e2);
+      bool elem1 = !(e1 < 0 || elem_delete[e1]);
+      bool elem2 = !(e2 < 0 || elem_delete[e2]);
+      // If there are two elements associated with the face on the boundary, then the mesh is periodic.
+      if (elem1 && elem2)
+      {
+        periodic_be = true;
+      }
+      if (periodic_be)
+      {
+        bdr_elem_delete[be] = true;
+        new_nbe--;
       }
     }
+
     if (new_ne < orig_mesh->GetNE())
     {
       Mpi::Print("Removed {:d} unmarked domain elements from the mesh\n",
@@ -2255,7 +2281,7 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
     int new_nbe_ext = 0, new_nbe_int = 0;
     for (int f = 0; f < orig_mesh->GetNumFaces(); f++)
     {
-      if (face_to_be.find(f) != face_to_be.end() || !iodata.boundaries.periodic.empty()) // TODO: I might need to fix this
+      if (face_to_be.find(f) != face_to_be.end())
       {
         continue;
       }