Fix race condition in GetFaceDofsFromAdjacentElement from mfem internal IsoparametricTransformation usage

Author: hughcars

CHANGELOG.md

@@ -108,6 +108,8 @@ The format of this changelog is based on
   - Fix bug in using `"MakeSimplex"` which would cause undefined behaviour for higher order
     meshes.
   - Fix bug when combining OpenMP and GPU builds in reduction operations over `Vector`.
+  - Fix race condition that would affect OpenMP parallelism with periodic boundaries (exact
+    and Floquet).
 
 ## [0.13.0] - 2024-05-20
 

palace/fem/libceed/restriction.cpp

@@ -51,20 +51,25 @@ mfem::Array<int> GetFaceDofsFromAdjacentElement(const mfem::FiniteElementSpace &
   // Get coordinates of element dofs.
   mfem::DenseMatrix elem_pm;
   const mfem::FiniteElement *fe_elem = fespace.GetFE(elem_id);
-  mfem::ElementTransformation *T = fespace.GetMesh()->GetElementTransformation(elem_id);
-  T->Transform(fe_elem->GetNodes(), elem_pm);
+  mfem::IsoparametricTransformation T;
+  fespace.GetMesh()->GetElementTransformation(elem_id, &T);
+  T.Transform(fe_elem->GetNodes(), elem_pm);
 
   // Find the dofs.
   double tol = 1E-5;
   mfem::Array<int> elem_dofs, dofs(P);
+  dofs = -123456789;
   fespace.GetElementDofs(elem_id, elem_dofs, dof_trans);
+
   for (int l = 0; l < P; l++)
   {
     double norm2_f = 0.0;
     for (int m = 0; m < face_pm.Height(); m++)
     {
       norm2_f += face_pm(m, l) * face_pm(m, l);
     }
+
+    bool found_match = false;
     for (int m = 0; m < elem_pm.Width(); m++)
     {
       double norm2_e = 0.0;
@@ -81,9 +86,28 @@ mfem::Array<int> GetFaceDofsFromAdjacentElement(const mfem::FiniteElementSpace &
       if (diff <= relative_tol)
       {
         dofs[l] = elem_dofs[m];
+        found_match = true;
+        break;
       }
     }
+
+    MFEM_ASSERT(found_match,
+                [&]()
+                {
+                  std::stringstream msg;
+                  msg << "l " << l << '\n';
+                  msg << "elem_dofs\n";
+                  for (auto x : elem_dofs)
+                    msg << x << ' ';
+
+                  msg << "\ndofs\n";
+                  for (auto x : dofs)
+                    msg << x << ' ';
+                  msg << '\n';
+                  return msg.str();
+                }());
   }
+
   return dofs;
 };