From bd3bf2628b0ef9c7b24f1d8a95125bb489667e44 Mon Sep 17 00:00:00 2001
From: Sebastian Grimberg <sjg@amazon.com>
Date: Thu, 14 Dec 2023 14:50:16 -0800
Subject: [PATCH 1/5] Minor cleanup of geodata helper functions and
 declarations, in preparation for Mesh class

---
 palace/drivers/basesolver.cpp                 |  15 +-
 palace/models/lumpedportoperator.cpp          |   7 +-
 palace/models/surfaceconductivityoperator.cpp |   7 +-
 palace/models/surfacecurrentoperator.cpp      |   7 +-
 palace/models/surfaceimpedanceoperator.cpp    |   7 +-
 palace/models/waveportoperator.cpp            |   7 +-
 palace/utils/configfile.cpp                   |  42 +-
 palace/utils/configfile.hpp                   |   9 +-
 palace/utils/geodata.cpp                      | 435 +++++++++---------
 palace/utils/geodata.hpp                      |  29 +-
 10 files changed, 279 insertions(+), 286 deletions(-)

diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp
index 62dd8f819..4495f41b9 100644
--- a/palace/drivers/basesolver.cpp
+++ b/palace/drivers/basesolver.cpp
@@ -157,9 +157,9 @@ void BaseSolver::SolveEstimateMarkRefine(
     constexpr bool refine = true, fix_orientation = true;
     mesh.back()->Finalize(refine, fix_orientation);
   }
+  MPI_Comm comm = mesh.back()->GetComm();
 
   // Perform initial solve and estimation.
-  MPI_Comm comm = mesh.back()->GetComm();
   auto [indicators, ntdof] = Solve(mesh);
   double err = indicators.Norml2(comm);
 
@@ -170,7 +170,7 @@ void BaseSolver::SolveEstimateMarkRefine(
     // Run out of iterations.
     ret |= (it >= refinement.max_it);
     // Run out of DOFs if a limit was set.
-    ret |= (refinement.max_size >= 1 && ntdof > refinement.max_size);
+    ret |= (refinement.max_size > 0 && ntdof > refinement.max_size);
     return ret;
   };
 
@@ -206,19 +206,20 @@ void BaseSolver::SolveEstimateMarkRefine(
         refinement.update_fraction);
 
     // Refine.
-    const auto initial_elem_count = mesh.back()->GetGlobalNE();
-    mesh.back()->GeneralRefinement(marked_elements, -1, refinement.max_nc_levels);
-    const auto final_elem_count = mesh.back()->GetGlobalNE();
+    auto &fine_mesh = *mesh.back();
+    const auto initial_elem_count = fine_mesh.GetGlobalNE();
+    fine_mesh.GeneralRefinement(marked_elements, -1, refinement.max_nc_levels);
+    const auto final_elem_count = fine_mesh.GetGlobalNE();
     Mpi::Print(" Mesh refinement added {:d} elements (initial: {}, final: {})\n",
                final_elem_count - initial_elem_count, initial_elem_count, final_elem_count);
 
     // Optionally rebalance and write the adapted mesh to file.
     const auto ratio_pre =
-        mesh::RebalanceMesh(iodata, mesh.back(), refinement.maximum_imbalance);
+        mesh::RebalanceMesh(fine_mesh, iodata, refinement.maximum_imbalance);
     if (ratio_pre > refinement.maximum_imbalance)
     {
       int min_elem, max_elem;
-      min_elem = max_elem = mesh.back()->GetNE();
+      min_elem = max_elem = fine_mesh.GetNE();
       Mpi::GlobalMin(1, &min_elem, comm);
       Mpi::GlobalMax(1, &max_elem, comm);
       const auto ratio_post = double(max_elem) / min_elem;
diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp
index eedf5a0e2..8e540b9c0 100644
--- a/palace/models/lumpedportoperator.cpp
+++ b/palace/models/lumpedportoperator.cpp
@@ -379,8 +379,7 @@ void LumpedPortOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &
           continue;
         }
         const int attr = i + 1;
-        mfem::Vector nor;
-        mesh::GetSurfaceNormal(mesh, attr, nor);
+        mfem::Vector normal = mesh::GetSurfaceNormal(mesh, attr);
         const double Rs = data.GetR() * data.GetToSquare(*elem);
         const double Ls = data.GetL() * data.GetToSquare(*elem);
         const double Cs = data.GetC() / data.GetToSquare(*elem);
@@ -418,11 +417,11 @@ void LumpedPortOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &
         }
         if (mesh.SpaceDimension() == 3)
         {
-          Mpi::Print(" n = ({:+.1f}, {:+.1f}, {:+.1f})", nor(0), nor(1), nor(2));
+          Mpi::Print(" n = ({:+.1f}, {:+.1f}, {:+.1f})", normal(0), normal(1), normal(2));
         }
         else
         {
-          Mpi::Print(" n = ({:+.1f}, {:+.1f})", nor(0), nor(1));
+          Mpi::Print(" n = ({:+.1f}, {:+.1f})", normal(0), normal(1));
         }
         Mpi::Print("\n");
       }
diff --git a/palace/models/surfaceconductivityoperator.cpp b/palace/models/surfaceconductivityoperator.cpp
index 705315f67..7c1e1591e 100644
--- a/palace/models/surfaceconductivityoperator.cpp
+++ b/palace/models/surfaceconductivityoperator.cpp
@@ -108,8 +108,7 @@ void SurfaceConductivityOperator::PrintBoundaryInfo(const IoData &iodata,
     if (conductivity_marker[i])
     {
       const int attr = i + 1;
-      mfem::Vector nor;
-      mesh::GetSurfaceNormal(mesh, attr, nor);
+      mfem::Vector normal = mesh::GetSurfaceNormal(mesh, attr);
       Mpi::Print(" {:d}: σ = {:.3e} S/m", attr,
                  iodata.DimensionalizeValue(IoData::ValueType::CONDUCTIVITY, bdr_sigma(i)));
       if (bdr_h(i) > 0.0)
@@ -119,11 +118,11 @@ void SurfaceConductivityOperator::PrintBoundaryInfo(const IoData &iodata,
       }
       if (mesh.SpaceDimension() == 3)
       {
-        Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", nor(0), nor(1), nor(2));
+        Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", normal(0), normal(1), normal(2));
       }
       else
       {
-        Mpi::Print(", n = ({:+.1f}, {:+.1f})", nor(0), nor(1));
+        Mpi::Print(", n = ({:+.1f}, {:+.1f})", normal(0), normal(1));
       }
       Mpi::Print("\n");
     }
diff --git a/palace/models/surfacecurrentoperator.cpp b/palace/models/surfacecurrentoperator.cpp
index bbcd1d17b..4c2dbea21 100644
--- a/palace/models/surfacecurrentoperator.cpp
+++ b/palace/models/surfacecurrentoperator.cpp
@@ -125,16 +125,15 @@ void SurfaceCurrentOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMe
           continue;
         }
         const int attr = i + 1;
-        mfem::Vector nor;
-        mesh::GetSurfaceNormal(mesh, attr, nor);
+        mfem::Vector normal = mesh::GetSurfaceNormal(mesh, attr);
         Mpi::Print(" {:d}: Index = {:d}", attr, idx);
         if (mesh.SpaceDimension() == 3)
         {
-          Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", nor(0), nor(1), nor(2));
+          Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", normal(0), normal(1), normal(2));
         }
         else
         {
-          Mpi::Print(", n = ({:+.1f}, {:+.1f})", nor(0), nor(1));
+          Mpi::Print(", n = ({:+.1f}, {:+.1f})", normal(0), normal(1));
         }
         Mpi::Print("\n");
       }
diff --git a/palace/models/surfaceimpedanceoperator.cpp b/palace/models/surfaceimpedanceoperator.cpp
index 4671cc6d0..8ff0eb999 100644
--- a/palace/models/surfaceimpedanceoperator.cpp
+++ b/palace/models/surfaceimpedanceoperator.cpp
@@ -109,8 +109,7 @@ void SurfaceImpedanceOperator::PrintBoundaryInfo(const IoData &iodata, mfem::Par
     if (impedance_marker[i])
     {
       const int attr = i + 1;
-      mfem::Vector nor;
-      mesh::GetSurfaceNormal(mesh, attr, nor);
+      mfem::Vector normal = mesh::GetSurfaceNormal(mesh, attr);
       bool comma = false;
       Mpi::Print(" {:d}:", attr);
       if (std::abs(Z_Rsinv(i)) > 0.0)
@@ -147,11 +146,11 @@ void SurfaceImpedanceOperator::PrintBoundaryInfo(const IoData &iodata, mfem::Par
       }
       if (mesh.SpaceDimension() == 3)
       {
-        Mpi::Print(" n = ({:+.1f}, {:+.1f}, {:+.1f})", nor(0), nor(1), nor(2));
+        Mpi::Print(" n = ({:+.1f}, {:+.1f}, {:+.1f})", normal(0), normal(1), normal(2));
       }
       else
       {
-        Mpi::Print(" n = ({:+.1f}, {:+.1f})", nor(0), nor(1));
+        Mpi::Print(" n = ({:+.1f}, {:+.1f})", normal(0), normal(1));
       }
       Mpi::Print("\n");
     }
diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp
index 54fbd1c5c..5526a2efd 100644
--- a/palace/models/waveportoperator.cpp
+++ b/palace/models/waveportoperator.cpp
@@ -1020,18 +1020,17 @@ void WavePortOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &me
         continue;
       }
       const int attr = i + 1;
-      mfem::Vector nor;
-      mesh::GetSurfaceNormal(mesh, attr, nor);
+      mfem::Vector normal = mesh::GetSurfaceNormal(mesh, attr);
       Mpi::Print(
           " {:d}: Index = {:d}, mode = {:d}, d = {:.3e} m", attr, idx, data.GetModeIndex(),
           iodata.DimensionalizeValue(IoData::ValueType::LENGTH, data.GetOffsetDistance()));
       if (mesh.SpaceDimension() == 3)
       {
-        Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", nor(0), nor(1), nor(2));
+        Mpi::Print(", n = ({:+.1f}, {:+.1f}, {:+.1f})", normal(0), normal(1), normal(2));
       }
       else
       {
-        Mpi::Print(", n = ({:+.1f}, {:+.1f})", nor(0), nor(1));
+        Mpi::Print(", n = ({:+.1f}, {:+.1f})", normal(0), normal(1));
       }
       Mpi::Print("\n");
     }
diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp
index ff887fe50..3ffc09ee4 100644
--- a/palace/utils/configfile.cpp
+++ b/palace/utils/configfile.cpp
@@ -614,8 +614,10 @@ void DomainPostData::SetUp(json &domains)
   // Store all unique postprocessing domain attributes.
   for (const auto &[idx, data] : energy)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
+  std::sort(attributes.begin(), attributes.end());
+  attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
 
   // Cleanup
   postpro->erase("Energy");
@@ -636,11 +638,14 @@ void DomainData::SetUp(json &config)
   // Store all unique domain attributes.
   for (const auto &data : materials)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
+  std::sort(attributes.begin(), attributes.end());
+  attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
   for (const auto &attr : postpro.attributes)
   {
-    MFEM_VERIFY(attributes.find(attr) != attributes.end(),
+    MFEM_VERIFY(std::lower_bound(attributes.begin(), attributes.end(), attr) !=
+                    attributes.end(),
                 "Domain postprocessing can only be enabled on domains which have a "
                 "corresponding \"Materials\" entry!");
   }
@@ -1295,19 +1300,21 @@ void BoundaryPostData::SetUp(json &boundaries)
   // Store all unique postprocessing boundary attributes.
   for (const auto &[idx, data] : capacitance)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
   for (const auto &[idx, data] : inductance)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
   for (const auto &[idx, data] : dielectric)
   {
     for (const auto &elem : data.elements)
     {
-      attributes.insert(elem.attributes.begin(), elem.attributes.end());
+      attributes.insert(attributes.end(), elem.attributes.begin(), elem.attributes.end());
     }
   }
+  std::sort(attributes.begin(), attributes.end());
+  attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
 
   // Cleanup
   postpro->erase("Capacitance");
@@ -1335,37 +1342,40 @@ void BoundaryData::SetUp(json &config)
   postpro.SetUp(*boundaries);
 
   // Store all unique boundary attributes.
-  attributes.insert(pec.attributes.begin(), pec.attributes.end());
-  attributes.insert(pmc.attributes.begin(), pmc.attributes.end());
-  attributes.insert(auxpec.attributes.begin(), auxpec.attributes.end());
-  attributes.insert(farfield.attributes.begin(), farfield.attributes.end());
+  attributes.insert(attributes.end(), pec.attributes.begin(), pec.attributes.end());
+  attributes.insert(attributes.end(), pmc.attributes.begin(), pmc.attributes.end());
+  attributes.insert(attributes.end(), auxpec.attributes.begin(), auxpec.attributes.end());
+  attributes.insert(attributes.end(), farfield.attributes.begin(),
+                    farfield.attributes.end());
   for (const auto &data : conductivity)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
   for (const auto &data : impedance)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
   for (const auto &[idx, data] : lumpedport)
   {
     for (const auto &elem : data.elements)
     {
-      attributes.insert(elem.attributes.begin(), elem.attributes.end());
+      attributes.insert(attributes.end(), elem.attributes.begin(), elem.attributes.end());
     }
   }
   for (const auto &[idx, data] : waveport)
   {
-    attributes.insert(data.attributes.begin(), data.attributes.end());
+    attributes.insert(attributes.end(), data.attributes.begin(), data.attributes.end());
   }
   for (const auto &[idx, data] : current)
   {
     for (const auto &elem : data.elements)
     {
-      attributes.insert(elem.attributes.begin(), elem.attributes.end());
+      attributes.insert(attributes.end(), elem.attributes.begin(), elem.attributes.end());
     }
   }
-  attributes.insert(postpro.attributes.begin(), postpro.attributes.end());
+  attributes.insert(attributes.end(), postpro.attributes.begin(), postpro.attributes.end());
+  std::sort(attributes.begin(), attributes.end());
+  attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
 
   // Cleanup
   boundaries->erase("PEC");
diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp
index bc7070173..5f73488bd 100644
--- a/palace/utils/configfile.hpp
+++ b/palace/utils/configfile.hpp
@@ -6,7 +6,6 @@
 
 #include <array>
 #include <map>
-#include <set>
 #include <string>
 #include <vector>
 #include <nlohmann/json_fwd.hpp>
@@ -289,7 +288,7 @@ struct DomainPostData
 {
 public:
   // Set of all postprocessing domain attributes.
-  std::set<int> attributes = {};
+  std::vector<int> attributes = {};
 
   // Domain postprocessing objects.
   DomainEnergyPostData energy;
@@ -302,7 +301,7 @@ struct DomainData
 {
 public:
   // Set of all domain attributes.
-  std::set<int> attributes = {};
+  std::vector<int> attributes = {};
 
   // Domain objects.
   DomainMaterialData materials = {};
@@ -521,7 +520,7 @@ struct BoundaryPostData
 {
 public:
   // Set of all postprocessing boundary attributes.
-  std::set<int> attributes = {};
+  std::vector<int> attributes = {};
 
   // Boundary postprocessing objects.
   CapacitancePostData capacitance = {};
@@ -535,7 +534,7 @@ struct BoundaryData
 {
 public:
   // Set of all boundary attributes.
-  std::set<int> attributes = {};
+  std::vector<int> attributes = {};
 
   // Boundary objects.
   PecBoundaryData pec = {};
diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp
index 7b5846aae..4e8b4761c 100644
--- a/palace/utils/geodata.cpp
+++ b/palace/utils/geodata.cpp
@@ -53,14 +53,9 @@ std::unique_ptr<mfem::ParMesh> DistributeMesh(MPI_Comm, const std::unique_ptr<mf
                                               const std::unique_ptr<int[]> & = nullptr,
                                               const std::string & = "");
 
-// Get list of domain and boundary attribute markers used in configuration file for mesh
-// cleaning.
-void GetUsedAttributeMarkers(const IoData &, int, int, mfem::Array<int> &,
-                             mfem::Array<int> &);
-
 // Rebalance a conformal mesh across processor ranks, using the MeshPartitioner. Gathers the
 // mesh onto the root rank before scattering the partitioned mesh.
-void RebalanceConformalMesh(std::unique_ptr<mfem::ParMesh> &, double, const std::string &);
+void RebalanceConformalMesh(mfem::ParMesh &);
 
 }  // namespace
 
@@ -257,10 +252,10 @@ void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>
       mfem::DenseMatrix pointmat;
       if (use_nodes)
       {
-        mfem::ElementTransformation *T = mesh.back()->GetElementTransformation(i);
-        mfem::Geometry::Type geo = mesh.back()->GetElementGeometry(i);
-        mfem::RefinedGeometry *RefG = mfem::GlobGeometryRefiner.Refine(geo, ref);
-        T->Transform(RefG->RefPts, pointmat);
+        mfem::ElementTransformation &T = *mesh.back()->GetElementTransformation(i);
+        mfem::Geometry::Type geom = mesh.back()->GetElementGeometry(i);
+        mfem::RefinedGeometry *RefG = mfem::GlobGeometryRefiner.Refine(geom, ref);
+        T.Transform(RefG->RefPts, pointmat);
       }
       else
       {
@@ -456,49 +451,53 @@ std::vector<mfem::Geometry::Type> ElementTypeInfo::GetGeomTypes() const
   return geom_types;
 }
 
-ElementTypeInfo CheckElements(mfem::Mesh &mesh)
+ElementTypeInfo CheckElements(const mfem::Mesh &mesh)
 {
   // MeshGenerator is reduced over the communicator. This checks for geometries on any
   // processor.
-  auto meshgen = mesh.MeshGenerator();
+  auto meshgen = const_cast<mfem::Mesh &>(mesh).MeshGenerator();
   return {bool(meshgen & 1), bool(meshgen & 2), bool(meshgen & 4), bool(meshgen & 8)};
 }
 
-void AttrToMarker(int max_attr, const mfem::Array<int> &attrs, mfem::Array<int> &marker)
+namespace
 {
-  MFEM_VERIFY(attrs.Size() == 0 || attrs.Max() <= max_attr,
-              "Invalid attribute number present (" << attrs.Max() << ")!");
-  marker.SetSize(max_attr);
-  if (attrs.Size() == 1 && attrs[0] == -1)
-  {
-    marker = 1;
-  }
-  else
-  {
-    marker = 0;
-    for (auto attr : attrs)
-    {
-      MFEM_VERIFY(attr > 0, "Attribute number less than one!");
-      MFEM_VERIFY(marker[attr - 1] == 0, "Repeate attribute in attribute list!");
-      marker[attr - 1] = 1;
-    }
-  }
+
+auto AttrListSize(const mfem::Array<int> &attr_list)
+{
+  return attr_list.Size();
 }
 
-void AttrToMarker(int max_attr, const std::vector<int> &attrs, mfem::Array<int> &marker)
+auto AttrListSize(const std::vector<int> &attr_list)
 {
-  MFEM_VERIFY(attrs.empty() || *std::max_element(attrs.begin(), attrs.end()) <= max_attr,
-              "Invalid attribute number present ("
-                  << *std::max_element(attrs.begin(), attrs.end()) << ")!");
+  return attr_list.size();
+}
+
+auto AttrListMax(const mfem::Array<int> &attr_list)
+{
+  return attr_list.Max();
+}
+
+auto AttrListMax(const std::vector<int> &attr_list)
+{
+  return *std::max_element(attr_list.begin(), attr_list.end());
+}
+
+}  // namespace
+
+template <typename T>
+void AttrToMarker(int max_attr, const T &attr_list, mfem::Array<int> &marker)
+{
+  MFEM_VERIFY(AttrListSize(attr_list) == 0 || AttrListMax(attr_list) <= max_attr,
+              "Invalid attribute number present (" << AttrListMax(attr_list) << ")!");
   marker.SetSize(max_attr);
-  if (attrs.size() == 1 && attrs[0] == -1)
+  if (AttrListSize(attr_list) == 1 && attr_list[0] == -1)
   {
     marker = 1;
   }
   else
   {
     marker = 0;
-    for (auto attr : attrs)
+    for (auto attr : attr_list)
     {
       MFEM_VERIFY(attr > 0, "Attribute number less than one!");
       MFEM_VERIFY(marker[attr - 1] == 0, "Repeate attribute in attribute list!");
@@ -507,15 +506,6 @@ void AttrToMarker(int max_attr, const std::vector<int> &attrs, mfem::Array<int>
   }
 }
 
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                               mfem::Vector &max)
-{
-  mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
-  marker = 0;
-  marker[attr - 1] = 1;
-  GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max);
-}
-
 void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                                bool bdr, mfem::Vector &min, mfem::Vector &max)
 {
@@ -529,7 +519,7 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
   }
   if (!mesh.GetNodes())
   {
-    auto BBUpdate = [&mesh, &dim, &min, &max](mfem::Array<int> &verts) -> void
+    auto BBUpdate = [&mesh, &dim, &min, &max](const mfem::Array<int> &verts) -> void
     {
       for (int j = 0; j < verts.Size(); j++)
       {
@@ -577,12 +567,12 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
   else
   {
     const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
-    auto BBUpdate = [&ref, &min, &max](mfem::ElementTransformation *T,
-                                       mfem::Geometry::Type &geo) -> void
+    auto BBUpdate = [&ref, &min, &max](mfem::ElementTransformation &T,
+                                       mfem::Geometry::Type &geom) -> void
     {
       mfem::DenseMatrix pointmat;
-      mfem::RefinedGeometry *RefG = mfem::GlobGeometryRefiner.Refine(geo, ref);
-      T->Transform(RefG->RefPts, pointmat);
+      mfem::RefinedGeometry *RefG = mfem::GlobGeometryRefiner.Refine(geom, ref);
+      T.Transform(RefG->RefPts, pointmat);
       for (int j = 0; j < pointmat.Width(); j++)
       {
         for (int d = 0; d < pointmat.Height(); d++)
@@ -606,9 +596,9 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
         {
           continue;
         }
-        mfem::ElementTransformation *T = mesh.GetBdrElementTransformation(i);
-        mfem::Geometry::Type geo = mesh.GetBdrElementGeometry(i);
-        BBUpdate(T, geo);
+        mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+        mfem::Geometry::Type geom = mesh.GetBdrElementGeometry(i);
+        BBUpdate(T, geom);
       }
     }
     else
@@ -619,16 +609,23 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
         {
           continue;
         }
-        mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
-        mfem::Geometry::Type geo = mesh.GetElementGeometry(i);
-        BBUpdate(T, geo);
+        mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+        mfem::Geometry::Type geom = mesh.GetElementGeometry(i);
+        BBUpdate(T, geom);
       }
     }
   }
-  auto *Min = min.HostReadWrite();
-  auto *Max = max.HostReadWrite();
-  Mpi::GlobalMin(dim, Min, mesh.GetComm());
-  Mpi::GlobalMax(dim, Max, mesh.GetComm());
+  Mpi::GlobalMin(dim, min.HostReadWrite(), mesh.GetComm());
+  Mpi::GlobalMax(dim, max.HostReadWrite(), mesh.GetComm());
+}
+
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
+                               mfem::Vector &max)
+{
+  mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
+  marker = 0;
+  marker[attr - 1] = 1;
+  GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max);
 }
 
 double BoundingBox::Area() const
@@ -727,8 +724,8 @@ int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
         {
           continue;
         }
-        mfem::ElementTransformation *T = mesh.GetBdrElementTransformation(i);
-        T->Transform(
+        mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+        T.Transform(
             mfem::GlobGeometryRefiner.Refine(mesh.GetBdrElementGeometry(i), ref)->RefPts,
             pointmat);
         for (int j = 0; j < pointmat.Width(); j++)
@@ -745,8 +742,8 @@ int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
         {
           continue;
         }
-        mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
-        T->Transform(
+        mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+        T.Transform(
             mfem::GlobGeometryRefiner.Refine(mesh.GetElementGeometry(i), ref)->RefPts,
             pointmat);
         for (int j = 0; j < pointmat.Width(); j++)
@@ -1091,52 +1088,72 @@ BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr)
   return GetBoundingBall(mesh, marker, bdr);
 }
 
-void GetSurfaceNormal(mfem::ParMesh &mesh, int attr, mfem::Vector &normal)
-{
-  mfem::Array<int> marker(mesh.bdr_attributes.Max());
-  marker = 0;
-  marker[attr - 1] = 1;
-  GetSurfaceNormal(mesh, marker, normal);
-}
-
-void GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
-                      mfem::Vector &normal)
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                              bool average)
 {
   int dim = mesh.SpaceDimension();
-  mfem::Vector nor(dim);
-  normal.SetSize(dim);
+  mfem::Vector loc_normal(dim), normal(dim);
   normal = 0.0;
   bool init = false;
-  for (int i = 0; i < mesh.GetNBE(); i++)
+  auto UpdateNormal = [&](mfem::ElementTransformation &T, mfem::Geometry::Type geom)
   {
-    if (!marker[mesh.GetBdrAttribute(i) - 1])
-    {
-      continue;
-    }
-    mfem::ElementTransformation *T = mesh.GetBdrElementTransformation(i);
-    const mfem::IntegrationPoint &ip =
-        mfem::Geometries.GetCenter(mesh.GetBdrElementGeometry(i));
-    T->SetIntPoint(&ip);
-    mfem::CalcOrtho(T->Jacobian(), nor);
+    const mfem::IntegrationPoint &ip = mfem::Geometries.GetCenter(geom);
+    T.SetIntPoint(&ip);
+    mfem::CalcOrtho(T.Jacobian(), loc_normal);
     if (!init)
     {
-      normal = nor;
+      normal = loc_normal;
       init = true;
     }
     else
     {
       // Check orientation and make sure consistent on this process. If a boundary has
       // conflicting normal definitions, use the first value.
-      if (nor * normal < 0.0)
+      if (loc_normal * normal < 0.0)
       {
-        normal -= nor;
+        normal -= loc_normal;
       }
       else
       {
-        normal += nor;
+        normal += loc_normal;
+      }
+    }
+  };
+  if (mesh.Dimension() == mesh.SpaceDimension())
+  {
+    // Loop over boundary elements.
+    for (int i = 0; i < mesh.GetNBE(); i++)
+    {
+      if (!marker[mesh.GetBdrAttribute(i) - 1])
+      {
+        continue;
+      }
+      mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+      UpdateNormal(T, mesh.GetBdrElementGeometry(i));
+      if (!average)
+      {
+        break;
       }
     }
   }
+  else
+  {
+    // Loop over domain elements.
+    for (int i = 0; i < mesh.GetNE(); i++)
+    {
+      if (!marker[mesh.GetAttribute(i) - 1])
+      {
+        continue;
+      }
+      mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+      UpdateNormal(T, mesh.GetElementGeometry(i));
+      if (!average)
+      {
+        break;
+      }
+    }
+  }
+
   // If different processors have different normal orientations, take that from the lowest
   // rank processor.
   MPI_Comm comm = mesh.GetComm();
@@ -1151,72 +1168,109 @@ void GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
   {
     // No boundary elements of attribute attr.
     normal = 0.0;
-    return;
+    return normal;
   }
   if (rank == Mpi::Rank(comm))
   {
     glob_normal = normal;
   }
+  Mpi::Broadcast(dim, glob_normal.HostReadWrite(), rank, comm);
+  if (average)
   {
-    auto *GlobNormal = glob_normal.HostReadWrite();
-    Mpi::Broadcast(dim, GlobNormal, rank, comm);
-  }
-  if (init && normal * glob_normal < 0.0)
-  {
-    normal.Neg();
+    if (init && normal * glob_normal < 0.0)
+    {
+      normal.Neg();
+    }
+    Mpi::GlobalSum(dim, normal.HostReadWrite(), comm);
   }
+  else
   {
-    auto *Normal = normal.HostReadWrite();
-    Mpi::GlobalSum(dim, Normal, comm);
+    normal = glob_normal;
   }
   normal /= normal.Norml2();
+
   // if (dim == 3)
   // {
   //   Mpi::Print(comm, " Surface normal {:d} = ({:+.3e}, {:+.3e}, {:+.3e})", attr,
-  //   normal(0),
-  //              normal(1), normal(2));
+  //              normal(0), normal(1), normal(2));
   // }
   // else
   // {
   //   Mpi::Print(comm, " Surface normal {:d} = ({:+.3e}, {:+.3e})", attr, normal(0),
   //              normal(1));
   // }
+
+  return normal;
 }
 
-double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh, double tol)
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, int attr, bool average)
+{
+  const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
+  mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
+  marker = 0;
+  marker[attr - 1] = 1;
+  return GetSurfaceNormal(mesh, marker, average);
+}
+
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, bool average)
+{
+  const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
+  const auto &attributes = bdr ? mesh.bdr_attributes : mesh.attributes;
+  return GetSurfaceNormal(mesh, AttrToMarker(attributes.Max(), attributes), average);
+}
+
+double RebalanceMesh(mfem::ParMesh &mesh, const IoData &iodata, double tol)
 {
   BlockTimer bt0(Timer::REBALANCE);
-  const bool save_adapt_mesh = iodata.model.refinement.save_adapt_mesh;
-  std::string serial_mesh_file;
-  if (save_adapt_mesh)
+  MPI_Comm comm = mesh.GetComm();
+  if (iodata.model.refinement.save_adapt_mesh)
   {
-    serial_mesh_file = iodata.problem.output;
-    if (serial_mesh_file.back() != '/')
+    // Create a separate serial mesh to write to disk.
+    std::string sfile = iodata.problem.output;
+    if (sfile.back() != '/')
     {
-      serial_mesh_file += '/';
+      sfile += '/';
     }
-    serial_mesh_file += "serial.mesh";
-  }
+    sfile += "serial.mesh";
 
-  MPI_Comm comm = mesh->GetComm();
-  if (Mpi::Size(comm) == 1)
-  {
-    if (save_adapt_mesh)
+    auto PrintSerial = [&](mfem::Mesh &smesh)
     {
       BlockTimer bt1(Timer::IO);
-      std::ofstream fo(serial_mesh_file);
-      fo.precision(MSH_FLT_PRECISION);
-      mesh::DimensionalizeMesh(*mesh, iodata.GetLengthScale());
-      mesh->mfem::Mesh::Print(fo);
-      mesh::NondimensionalizeMesh(*mesh, iodata.GetLengthScale());
+      if (Mpi::Root(comm))
+      {
+        std::ofstream fo(sfile);
+        // mfem::ofgzstream fo(sfile, true);  // Use zlib compression if available
+        // fo << std::fixed;
+        fo << std::scientific;
+        fo.precision(MSH_FLT_PRECISION);
+        mesh::DimensionalizeMesh(smesh, iodata.GetLengthScale());
+        smesh.Mesh::Print(fo);  // Do not need to nondimensionalize the temporary mesh
+      }
+    };
+
+    if (mesh.Nonconforming())
+    {
+      mfem::ParMesh smesh(mesh);
+      mfem::Array<int> serial_partition(mesh.GetNE());
+      serial_partition = 0;
+      smesh.Rebalance(serial_partition);
+      PrintSerial(smesh);
     }
-    return 1.0;
+    else
+    {
+      mfem::Mesh smesh(mesh.GetSerialMesh(0));
+      PrintSerial(smesh);
+    }
+    Mpi::Barrier(comm);
   }
 
   // If there is more than one processor, may perform rebalancing.
-  mesh->ExchangeFaceNbrData();
+  if (Mpi::Size(comm) == 1)
+  {
+    return 1.0;
+  }
   int min_elem, max_elem;
-  min_elem = max_elem = mesh->GetNE();
+  min_elem = max_elem = mesh.GetNE();
   Mpi::GlobalMin(1, &min_elem, comm);
   Mpi::GlobalMax(1, &max_elem, comm);
   const double ratio = double(max_elem) / min_elem;
@@ -1228,72 +1282,25 @@ double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh,
   }
   if (ratio > tol)
   {
-    if (mesh->Nonconforming() && save_adapt_mesh)
+    mesh.ExchangeFaceNbrData();
+    if (mesh.Nonconforming())
     {
-      // Do not need to duplicate the mesh, as rebalancing will undo this.
-      mfem::Array<int> serial_partition(mesh->GetNE());
-      serial_partition = 0;
-      mesh->Rebalance(serial_partition);
-      BlockTimer bt1(Timer::IO);
-      if (Mpi::Root(comm))
-      {
-        std::ofstream fo(serial_mesh_file);
-        fo.precision(MSH_FLT_PRECISION);
-        mesh::DimensionalizeMesh(*mesh, iodata.GetLengthScale());
-        mesh->Mesh::Print(fo);
-        mesh::NondimensionalizeMesh(*mesh, iodata.GetLengthScale());
-      }
-      Mpi::Barrier(comm);
-    }
-    if (mesh->Nonconforming())
-    {
-      mesh->Rebalance();
+      mesh.Rebalance();
     }
     else
     {
       // Without access to a refinement tree, partitioning must be done on the root
       // processor and then redistributed.
-      RebalanceConformalMesh(mesh, iodata.GetLengthScale(), serial_mesh_file);
+      RebalanceConformalMesh(mesh);
     }
+    mesh.ExchangeFaceNbrData();
   }
-  else if (save_adapt_mesh)
-  {
-    // Given no rebalancing will be done, need to handle the serial write more carefully.
-    // This requires creating a separate serial mesh.
-    if (mesh->Nonconforming())
-    {
-      mfem::ParMesh smesh(*mesh);
-      mfem::Array<int> serial_partition(mesh->GetNE());
-      serial_partition = 0;
-      smesh.Rebalance(serial_partition);
-      BlockTimer bt1(Timer::IO);
-      if (Mpi::Root(comm))
-      {
-        std::ofstream fo(serial_mesh_file);
-        fo.precision(MSH_FLT_PRECISION);
-        mesh::DimensionalizeMesh(smesh, iodata.GetLengthScale());
-        smesh.Mesh::Print(fo);  // Do not need to nondimensionalize the temporary mesh
-      }
-      Mpi::Barrier(comm);
-    }
-    else
-    {
-      auto smesh = std::make_unique<mfem::Mesh>(mesh->GetSerialMesh(0));
-      BlockTimer bt1(Timer::IO);
-      if (Mpi::Rank(comm) == 0)
-      {
-        std::ofstream fo(serial_mesh_file);
-        fo.precision(MSH_FLT_PRECISION);
-        mesh::DimensionalizeMesh(*smesh, iodata.GetLengthScale());
-        smesh->Print(fo);  // Do not need to nondimensionalize the temporary mesh
-      }
-      Mpi::Barrier(comm);
-    }
-  }
-  mesh->ExchangeFaceNbrData();
   return ratio;
 }
 
+template void AttrToMarker(int, const mfem::Array<int> &, mfem::Array<int> &);
+template void AttrToMarker(int, const std::vector<int> &, mfem::Array<int> &);
+
 }  // namespace mesh
 
 namespace
@@ -1477,11 +1484,12 @@ std::map<int, std::array<int, 2>> CheckMesh(mfem::Mesh &orig_mesh,
               "Nonconforming or 2D meshes have not been tested yet!");
   MFEM_VERIFY(dynamic_cast<mfem::ParMesh *>(&orig_mesh) == nullptr,
               "This function does not work for ParMesh");
-  mfem::Array<int> mat_marker, bdr_marker;
-  GetUsedAttributeMarkers(
-      iodata, orig_mesh.attributes.Size() ? orig_mesh.attributes.Max() : 0,
-      orig_mesh.bdr_attributes.Size() ? orig_mesh.bdr_attributes.Max() : 0, mat_marker,
-      bdr_marker);
+  auto mat_marker =
+      mesh::AttrToMarker(orig_mesh.attributes.Size() ? orig_mesh.attributes.Max() : 0,
+                         iodata.domains.attributes);
+  auto bdr_marker = mesh::AttrToMarker(
+      orig_mesh.bdr_attributes.Size() ? orig_mesh.bdr_attributes.Max() : 0,
+      iodata.boundaries.attributes);
   bool warn = false;
   for (int be = 0; be < orig_mesh.GetNBE(); be++)
   {
@@ -1930,41 +1938,22 @@ std::unique_ptr<mfem::ParMesh> DistributeMesh(MPI_Comm comm,
   }
 }
 
-void GetUsedAttributeMarkers(const IoData &iodata, int n_mat, int n_bdr,
-                             mfem::Array<int> &mat_marker, mfem::Array<int> &bdr_marker)
-{
-  mfem::Array<int> mat_attr, bdr_attr;
-  mat_attr.Reserve(static_cast<int>(iodata.domains.attributes.size()));
-  for (auto attr : iodata.domains.attributes)
-  {
-    mat_attr.Append(attr);
-  }
-  bdr_attr.Reserve(static_cast<int>(iodata.boundaries.attributes.size()));
-  for (auto attr : iodata.boundaries.attributes)
-  {
-    bdr_attr.Append(attr);
-  }
-  mesh::AttrToMarker(n_mat, mat_attr, mat_marker);
-  mesh::AttrToMarker(n_bdr, bdr_attr, bdr_marker);
-}
-
-void RebalanceConformalMesh(std::unique_ptr<mfem::ParMesh> &pmesh, double length_scale,
-                            const std::string &serial_mesh_file)
+void RebalanceConformalMesh(mfem::ParMesh &pmesh)
 {
   // Write the parallel mesh to a stream as a serial mesh, then read back in and partition
   // using METIS.
-  MPI_Comm comm = pmesh->GetComm();
+  MPI_Comm comm = pmesh.GetComm();
   constexpr bool generate_edges = true, refine = true, fix_orientation = true,
                  generate_bdr = false;
   std::unique_ptr<mfem::Mesh> smesh;
-  std::unique_ptr<int[]> partitioning;
   if constexpr (false)
   {
     // Write the serial mesh to a stream and read that through the Mesh constructor.
-    std::stringstream fo;
+    std::ostringstream fo(std::stringstream::out);
+    // fo << std::fixed;
+    fo << std::scientific;
     fo.precision(MSH_FLT_PRECISION);
-    pmesh->PrintAsSerial(fo);
-    pmesh.reset();
+    pmesh.PrintAsSerial(fo);
     if (Mpi::Root(comm))
     {
       smesh = std::make_unique<mfem::Mesh>(fo, generate_edges, refine, fix_orientation);
@@ -1973,34 +1962,26 @@ void RebalanceConformalMesh(std::unique_ptr<mfem::ParMesh> &pmesh, double length
   else
   {
     // Directly ingest the generated Mesh and release the no longer needed memory.
-    smesh = std::make_unique<mfem::Mesh>(pmesh->GetSerialMesh(0));
-    pmesh.reset();
-    if (!Mpi::Root(comm))
+    smesh = std::make_unique<mfem::Mesh>(pmesh.GetSerialMesh(0));
+    if (Mpi::Root(comm))
+    {
+      smesh->FinalizeTopology(generate_bdr);
+      smesh->Finalize(refine, fix_orientation);
+    }
+    else
     {
       smesh.reset();
     }
   }
+
+  // (Re)-construct the parallel mesh.
+  std::unique_ptr<int[]> partitioning;
   if (Mpi::Root(comm))
   {
-    smesh->FinalizeTopology(generate_bdr);
-    smesh->Finalize(refine, fix_orientation);
     partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm));
   }
-
-  // Construct the parallel mesh.
-  pmesh = DistributeMesh(comm, smesh, partitioning);
-  if (!serial_mesh_file.empty())
-  {
-    BlockTimer bt(Timer::IO);
-    if (Mpi::Root(comm))
-    {
-      std::ofstream fo(serial_mesh_file);
-      fo.precision(MSH_FLT_PRECISION);
-      mesh::DimensionalizeMesh(*smesh, length_scale);
-      smesh->Print(fo);  // Do not need to nondimensionalize the temporary mesh
-    }
-    Mpi::Barrier(comm);
-  }
+  auto new_pmesh = DistributeMesh(comm, smesh, partitioning);
+  pmesh = std::move(*new_pmesh);
 }
 
 }  // namespace
diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp
index cc245e5e5..54bbc454a 100644
--- a/palace/utils/geodata.hpp
+++ b/palace/utils/geodata.hpp
@@ -50,20 +50,27 @@ struct ElementTypeInfo
 };
 
 // Simplified helper for describing the element types in a (Par)Mesh.
-ElementTypeInfo CheckElements(mfem::Mesh &mesh);
+ElementTypeInfo CheckElements(const mfem::Mesh &mesh);
 
 // Helper function to convert a set of attribute numbers to a marker array. The marker array
 // will be of size max_attr and it will contain only zeroes and ones. Ones indicate which
-// attribute numbers are present in the attrs array. In the special case when attrs has a
+// attribute numbers are present in the list array. In the special case when list has a
 // single entry equal to -1 the marker array will contain all ones.
-void AttrToMarker(int max_attr, const mfem::Array<int> &attrs, mfem::Array<int> &marker);
-void AttrToMarker(int max_attr, const std::vector<int> &attrs, mfem::Array<int> &marker);
+template <typename T>
+void AttrToMarker(int max_attr, const T &attr_list, mfem::Array<int> &marker);
+template <typename T>
+mfem::Array<int> AttrToMarker(int max_attr, const T &attr_list)
+{
+  mfem::Array<int> marker;
+  AttrToMarker(max_attr, attr_list, marker);
+  return marker;
+}
 
 // Helper function to construct the bounding box for all elements with the given attribute.
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                               mfem::Vector &max);
 void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                                bool bdr, mfem::Vector &min, mfem::Vector &max);
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
+                               mfem::Vector &max);
 
 // Struct describing a bounding box in terms of the center and face normals. The normals
 // specify the direction from the center of the box.
@@ -132,14 +139,14 @@ BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr);
 
 // Helper function to compute the average surface normal for all elements with the given
 // attribute.
-void GetSurfaceNormal(mfem::ParMesh &mesh, int attr, mfem::Vector &normal);
-void GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
-                      mfem::Vector &normal);
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                              bool average = true);
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, int attr, bool average = true);
+mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, bool average = true);
 
 // Helper function responsible for rebalancing the mesh, and optionally writing meshes from
 // the intermediate stages to disk. Returns the imbalance ratio before rebalancing.
-double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh,
-                     double tol);
+double RebalanceMesh(mfem::ParMesh &mesh, const IoData &iodata, double tol);
 
 }  // namespace mesh
 

From 43e4b365afa702b226f8ca3aa199fda09626e400 Mon Sep 17 00:00:00 2001
From: Sebastian Grimberg <sjg@amazon.com>
Date: Thu, 14 Dec 2023 15:13:37 -0800
Subject: [PATCH 2/5] Const correctness for mesh and thread safety

---
 palace/models/lumpedportoperator.cpp          |  2 +-
 palace/models/lumpedportoperator.hpp          |  2 +-
 palace/models/surfaceconductivityoperator.cpp |  4 +-
 palace/models/surfaceconductivityoperator.hpp |  4 +-
 palace/models/surfacecurrentoperator.cpp      |  3 +-
 palace/models/surfacecurrentoperator.hpp      |  2 +-
 palace/models/surfaceimpedanceoperator.cpp    |  5 +-
 palace/models/surfaceimpedanceoperator.hpp    |  4 +-
 palace/models/waveportoperator.cpp            |  2 +-
 palace/models/waveportoperator.hpp            |  2 +-
 palace/utils/geodata.cpp                      | 57 ++++++++++---------
 palace/utils/geodata.hpp                      | 28 ++++-----
 12 files changed, 62 insertions(+), 53 deletions(-)

diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp
index 8e540b9c0..d7e98f6f4 100644
--- a/palace/models/lumpedportoperator.cpp
+++ b/palace/models/lumpedportoperator.cpp
@@ -360,7 +360,7 @@ void LumpedPortOperator::SetUpBoundaryProperties(const IoData &iodata,
   }
 }
 
-void LumpedPortOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh)
+void LumpedPortOperator::PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh)
 {
   // Print out BC info for all port attributes.
   if (ports.empty())
diff --git a/palace/models/lumpedportoperator.hpp b/palace/models/lumpedportoperator.hpp
index 8c5fa30d9..ec08264af 100644
--- a/palace/models/lumpedportoperator.hpp
+++ b/palace/models/lumpedportoperator.hpp
@@ -90,7 +90,7 @@ class LumpedPortOperator
   mfem::Array<int> port_marker, port_Rs_marker, port_Ls_marker, port_Cs_marker;
   void SetUpBoundaryProperties(const IoData &iodata,
                                mfem::ParFiniteElementSpace &h1_fespace);
-  void PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh);
+  void PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh);
 
 public:
   LumpedPortOperator(const IoData &iodata, mfem::ParFiniteElementSpace &h1_fespace);
diff --git a/palace/models/surfaceconductivityoperator.cpp b/palace/models/surfaceconductivityoperator.cpp
index 7c1e1591e..329d4c6c3 100644
--- a/palace/models/surfaceconductivityoperator.cpp
+++ b/palace/models/surfaceconductivityoperator.cpp
@@ -14,7 +14,7 @@ namespace palace
 using namespace std::complex_literals;
 
 SurfaceConductivityOperator::SurfaceConductivityOperator(const IoData &iodata,
-                                                         mfem::ParMesh &mesh)
+                                                         const mfem::ParMesh &mesh)
 {
   // Set up finite conductivity boundary conditions.
   SetUpBoundaryProperties(iodata, mesh);
@@ -96,7 +96,7 @@ void SurfaceConductivityOperator::SetUpBoundaryProperties(const IoData &iodata,
 }
 
 void SurfaceConductivityOperator::PrintBoundaryInfo(const IoData &iodata,
-                                                    mfem::ParMesh &mesh)
+                                                    const mfem::ParMesh &mesh)
 {
   if (conductivity_marker.Size() && conductivity_marker.Max() == 0)
   {
diff --git a/palace/models/surfaceconductivityoperator.hpp b/palace/models/surfaceconductivityoperator.hpp
index cfe40e717..7648e9647 100644
--- a/palace/models/surfaceconductivityoperator.hpp
+++ b/palace/models/surfaceconductivityoperator.hpp
@@ -23,10 +23,10 @@ class SurfaceConductivityOperator
   mfem::Vector bdr_sigma, bdr_mu, bdr_h;
   mfem::Array<int> conductivity_marker;
   void SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh);
-  void PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh);
+  void PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh);
 
 public:
-  SurfaceConductivityOperator(const IoData &iodata, mfem::ParMesh &mesh);
+  SurfaceConductivityOperator(const IoData &iodata, const mfem::ParMesh &mesh);
 
   // Returns array marking finite conductivity boundary attributes.
   const mfem::Array<int> &GetMarker() const { return conductivity_marker; }
diff --git a/palace/models/surfacecurrentoperator.cpp b/palace/models/surfacecurrentoperator.cpp
index 4c2dbea21..5c546a0a0 100644
--- a/palace/models/surfacecurrentoperator.cpp
+++ b/palace/models/surfacecurrentoperator.cpp
@@ -107,7 +107,8 @@ void SurfaceCurrentOperator::SetUpBoundaryProperties(
   }
 }
 
-void SurfaceCurrentOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh)
+void SurfaceCurrentOperator::PrintBoundaryInfo(const IoData &iodata,
+                                               const mfem::ParMesh &mesh)
 {
   if (sources.empty())
   {
diff --git a/palace/models/surfacecurrentoperator.hpp b/palace/models/surfacecurrentoperator.hpp
index 949cec765..719b91aa3 100644
--- a/palace/models/surfacecurrentoperator.hpp
+++ b/palace/models/surfacecurrentoperator.hpp
@@ -58,7 +58,7 @@ class SurfaceCurrentOperator
   mfem::Array<int> source_marker;
   void SetUpBoundaryProperties(const IoData &iodata,
                                mfem::ParFiniteElementSpace &h1_fespace);
-  void PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh);
+  void PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh);
 
 public:
   SurfaceCurrentOperator(const IoData &iodata, mfem::ParFiniteElementSpace &h1_fespace);
diff --git a/palace/models/surfaceimpedanceoperator.cpp b/palace/models/surfaceimpedanceoperator.cpp
index 8ff0eb999..07b984c0c 100644
--- a/palace/models/surfaceimpedanceoperator.cpp
+++ b/palace/models/surfaceimpedanceoperator.cpp
@@ -12,7 +12,7 @@ namespace palace
 {
 
 SurfaceImpedanceOperator::SurfaceImpedanceOperator(const IoData &iodata,
-                                                   mfem::ParMesh &mesh)
+                                                   const mfem::ParMesh &mesh)
 {
   // Set up impedance boundary conditions.
   SetUpBoundaryProperties(iodata, mesh);
@@ -97,7 +97,8 @@ void SurfaceImpedanceOperator::SetUpBoundaryProperties(const IoData &iodata,
   mesh::AttrToMarker(bdr_attr_max, impedance_Cs_bcs, impedance_Cs_marker);
 }
 
-void SurfaceImpedanceOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh)
+void SurfaceImpedanceOperator::PrintBoundaryInfo(const IoData &iodata,
+                                                 const mfem::ParMesh &mesh)
 {
   if (impedance_marker.Size() && impedance_marker.Max() == 0)
   {
diff --git a/palace/models/surfaceimpedanceoperator.hpp b/palace/models/surfaceimpedanceoperator.hpp
index 6ad90ab42..17b8620ef 100644
--- a/palace/models/surfaceimpedanceoperator.hpp
+++ b/palace/models/surfaceimpedanceoperator.hpp
@@ -24,10 +24,10 @@ class SurfaceImpedanceOperator
   mfem::Array<int> impedance_marker, impedance_Rs_marker, impedance_Ls_marker,
       impedance_Cs_marker;
   void SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh);
-  void PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh);
+  void PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh);
 
 public:
-  SurfaceImpedanceOperator(const IoData &iodata, mfem::ParMesh &mesh);
+  SurfaceImpedanceOperator(const IoData &iodata, const mfem::ParMesh &mesh);
 
   // Returns array marking surface impedance attributes.
   const mfem::Array<int> &GetMarker() const { return impedance_marker; }
diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp
index 5526a2efd..e45f0b335 100644
--- a/palace/models/waveportoperator.cpp
+++ b/palace/models/waveportoperator.cpp
@@ -1003,7 +1003,7 @@ void WavePortOperator::SetUpBoundaryProperties(
   }
 }
 
-void WavePortOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh)
+void WavePortOperator::PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh)
 {
   // Print out BC info for all port attributes.
   if (ports.empty())
diff --git a/palace/models/waveportoperator.hpp b/palace/models/waveportoperator.hpp
index c2a7b1f6a..4bb368b5a 100644
--- a/palace/models/waveportoperator.hpp
+++ b/palace/models/waveportoperator.hpp
@@ -142,7 +142,7 @@ class WavePortOperator
   void SetUpBoundaryProperties(const IoData &iodata, const MaterialOperator &mat_op,
                                const mfem::ParFiniteElementSpace &nd_fespace,
                                const mfem::ParFiniteElementSpace &h1_fespace);
-  void PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh);
+  void PrintBoundaryInfo(const IoData &iodata, const mfem::ParMesh &mesh);
 
   // Compute boundary modes for all wave port boundaries at the specified frequency.
   void Initialize(double omega);
diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp
index 4e8b4761c..f8cca97c4 100644
--- a/palace/utils/geodata.cpp
+++ b/palace/utils/geodata.cpp
@@ -506,7 +506,7 @@ void AttrToMarker(int max_attr, const T &attr_list, mfem::Array<int> &marker)
   }
 }
 
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                                bool bdr, mfem::Vector &min, mfem::Vector &max)
 {
   int dim = mesh.SpaceDimension();
@@ -566,9 +566,8 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
   }
   else
   {
-    const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
-    auto BBUpdate = [&ref, &min, &max](mfem::ElementTransformation &T,
-                                       mfem::Geometry::Type &geom) -> void
+    auto BBUpdate = [&min, &max](mfem::ElementTransformation &T, mfem::Geometry::Type &geom,
+                                 int ref) -> void
     {
       mfem::DenseMatrix pointmat;
       mfem::RefinedGeometry *RefG = mfem::GlobGeometryRefiner.Refine(geom, ref);
@@ -588,6 +587,8 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
         }
       }
     };
+    const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
+    mfem::IsoparametricTransformation T;
     if (bdr)
     {
       for (int i = 0; i < mesh.GetNBE(); i++)
@@ -596,9 +597,9 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
         {
           continue;
         }
-        mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+        mesh.GetBdrElementTransformation(i, &T);
         mfem::Geometry::Type geom = mesh.GetBdrElementGeometry(i);
-        BBUpdate(T, geom);
+        BBUpdate(T, geom, ref);
       }
     }
     else
@@ -609,9 +610,9 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
         {
           continue;
         }
-        mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+        mesh.GetElementTransformation(i, &T);
         mfem::Geometry::Type geom = mesh.GetElementGeometry(i);
-        BBUpdate(T, geom);
+        BBUpdate(T, geom, ref);
       }
     }
   }
@@ -619,8 +620,8 @@ void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &mark
   Mpi::GlobalMax(dim, max.HostReadWrite(), mesh.GetComm());
 }
 
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                               mfem::Vector &max)
+void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr,
+                               mfem::Vector &min, mfem::Vector &max)
 {
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
   marker = 0;
@@ -673,8 +674,8 @@ bool EigenLE(const Eigen::Vector3d &x, const Eigen::Vector3d &y)
 // bounding balls. Returns the dominant rank, for which the vertices argument will be
 // filled, while all other ranks will have an empty vector. Vertices are de-duplicated to a
 // certain floating point precision.
-int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
-                            std::vector<Eigen::Vector3d> &vertices)
+int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                            bool bdr, std::vector<Eigen::Vector3d> &vertices)
 {
   std::set<int> vertex_indices;
   if (!mesh.GetNodes())
@@ -716,6 +717,7 @@ int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
     // Nonlinear mesh, need to process point matrices.
     const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder();
     mfem::DenseMatrix pointmat;  // 3 x N
+    mfem::IsoparametricTransformation T;
     if (bdr)
     {
       for (int i = 0; i < mesh.GetNBE(); i++)
@@ -724,7 +726,7 @@ int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
         {
           continue;
         }
-        mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+        mesh.GetBdrElementTransformation(i, &T);
         T.Transform(
             mfem::GlobGeometryRefiner.Refine(mesh.GetBdrElementGeometry(i), ref)->RefPts,
             pointmat);
@@ -742,7 +744,7 @@ int CollectPointCloudOnRoot(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
         {
           continue;
         }
-        mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+        mesh.GetElementTransformation(i, &T);
         T.Transform(
             mfem::GlobGeometryRefiner.Refine(mesh.GetElementGeometry(i), ref)->RefPts,
             pointmat);
@@ -1041,15 +1043,15 @@ double LengthFromPointCloud(MPI_Comm comm, const std::vector<Eigen::Vector3d> &v
 
 }  // namespace
 
-double GetProjectedLength(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
-                          const std::array<double, 3> &dir)
+double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                          bool bdr, const std::array<double, 3> &dir)
 {
   std::vector<Eigen::Vector3d> vertices;
   int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
   return LengthFromPointCloud(mesh.GetComm(), vertices, dominant_rank, dir);
 }
 
-double GetProjectedLength(mfem::ParMesh &mesh, int attr, bool bdr,
+double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr,
                           const std::array<double, 3> &dir)
 {
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
@@ -1058,14 +1060,15 @@ double GetProjectedLength(mfem::ParMesh &mesh, int attr, bool bdr,
   return GetProjectedLength(mesh, marker, bdr, dir);
 }
 
-BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr)
+BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                           bool bdr)
 {
   std::vector<Eigen::Vector3d> vertices;
   int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
   return BoundingBoxFromPointCloud(mesh.GetComm(), vertices, dominant_rank);
 }
 
-BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr)
+BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr)
 {
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
   marker = 0;
@@ -1073,14 +1076,15 @@ BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr)
   return GetBoundingBox(mesh, marker, bdr);
 }
 
-BoundingBall GetBoundingBall(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr)
+BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                             bool bdr)
 {
   std::vector<Eigen::Vector3d> vertices;
   int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
   return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank);
 }
 
-BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr)
+BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr)
 {
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
   marker = 0;
@@ -1088,10 +1092,11 @@ BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr)
   return GetBoundingBall(mesh, marker, bdr);
 }
 
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                               bool average)
 {
   int dim = mesh.SpaceDimension();
+  mfem::IsoparametricTransformation T;
   mfem::Vector loc_normal(dim), normal(dim);
   normal = 0.0;
   bool init = false;
@@ -1128,7 +1133,7 @@ mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marke
       {
         continue;
       }
-      mfem::ElementTransformation &T = *mesh.GetBdrElementTransformation(i);
+      mesh.GetBdrElementTransformation(i, &T);
       UpdateNormal(T, mesh.GetBdrElementGeometry(i));
       if (!average)
       {
@@ -1145,7 +1150,7 @@ mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marke
       {
         continue;
       }
-      mfem::ElementTransformation &T = *mesh.GetElementTransformation(i);
+      mesh.GetElementTransformation(i, &T);
       UpdateNormal(T, mesh.GetElementGeometry(i));
       if (!average)
       {
@@ -1203,7 +1208,7 @@ mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marke
   return normal;
 }
 
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, int attr, bool average)
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr, bool average)
 {
   const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
@@ -1212,7 +1217,7 @@ mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, int attr, bool average)
   return GetSurfaceNormal(mesh, marker, average);
 }
 
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, bool average)
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average)
 {
   const bool bdr = (mesh.Dimension() == mesh.SpaceDimension());
   const auto &attributes = bdr ? mesh.bdr_attributes : mesh.attributes;
diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp
index 54bbc454a..4ed95f95b 100644
--- a/palace/utils/geodata.hpp
+++ b/palace/utils/geodata.hpp
@@ -67,10 +67,10 @@ mfem::Array<int> AttrToMarker(int max_attr, const T &attr_list)
 }
 
 // Helper function to construct the bounding box for all elements with the given attribute.
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                                bool bdr, mfem::Vector &min, mfem::Vector &max);
-void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                               mfem::Vector &max);
+void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr,
+                               mfem::Vector &min, mfem::Vector &max);
 
 // Struct describing a bounding box in terms of the center and face normals. The normals
 // specify the direction from the center of the box.
@@ -123,26 +123,28 @@ struct BoundingBall
 };
 
 // Helper functions for computing bounding boxes from a mesh and markers.
-BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr);
-BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr);
+BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                           bool bdr);
+BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr);
 
 // Helper function for computing the direction aligned length of a marked group.
-double GetProjectedLength(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
-                          const std::array<double, 3> &dir);
-double GetProjectedLength(mfem::ParMesh &mesh, int attr, bool bdr,
+double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                          bool bdr, const std::array<double, 3> &dir);
+double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr,
                           const std::array<double, 3> &dir);
 
 // Given a mesh and a marker, compute the diameter of a bounding circle/sphere, assuming
 // that the extrema points are in the marked group.
-BoundingBall GetBoundingBall(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr);
-BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr);
+BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                             bool bdr);
+BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr);
 
 // Helper function to compute the average surface normal for all elements with the given
 // attribute.
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array<int> &marker,
                               bool average = true);
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, int attr, bool average = true);
-mfem::Vector GetSurfaceNormal(mfem::ParMesh &mesh, bool average = true);
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr, bool average = true);
+mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average = true);
 
 // Helper function responsible for rebalancing the mesh, and optionally writing meshes from
 // the intermediate stages to disk. Returns the imbalance ratio before rebalancing.

From c586bab75836a565ce970e6d340edc7c886ed454 Mon Sep 17 00:00:00 2001
From: Sebastian Grimberg <sjg@amazon.com>
Date: Wed, 20 Dec 2023 11:22:59 -0800
Subject: [PATCH 3/5] Fix unneeded const_cast after #151

---
 palace/utils/geodata.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp
index f8cca97c4..adde903fa 100644
--- a/palace/utils/geodata.cpp
+++ b/palace/utils/geodata.cpp
@@ -455,7 +455,7 @@ ElementTypeInfo CheckElements(const mfem::Mesh &mesh)
 {
   // MeshGenerator is reduced over the communicator. This checks for geometries on any
   // processor.
-  auto meshgen = const_cast<mfem::Mesh &>(mesh).MeshGenerator();
+  auto meshgen = mesh.MeshGenerator();
   return {bool(meshgen & 1), bool(meshgen & 2), bool(meshgen & 4), bool(meshgen & 8)};
 }
 

From bfe264397f6512cd85a7e23ae301333b8914611f Mon Sep 17 00:00:00 2001
From: Sebastian Grimberg <sjg@amazon.com>
Date: Wed, 20 Dec 2023 11:25:48 -0800
Subject: [PATCH 4/5] Add shrink_to_fit for vector class members

---
 palace/utils/configfile.cpp | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp
index 3ffc09ee4..030e15525 100644
--- a/palace/utils/configfile.cpp
+++ b/palace/utils/configfile.cpp
@@ -618,6 +618,7 @@ void DomainPostData::SetUp(json &domains)
   }
   std::sort(attributes.begin(), attributes.end());
   attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
+  attributes.shrink_to_fit();
 
   // Cleanup
   postpro->erase("Energy");
@@ -642,6 +643,7 @@ void DomainData::SetUp(json &config)
   }
   std::sort(attributes.begin(), attributes.end());
   attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
+  attributes.shrink_to_fit();
   for (const auto &attr : postpro.attributes)
   {
     MFEM_VERIFY(std::lower_bound(attributes.begin(), attributes.end(), attr) !=
@@ -1315,6 +1317,7 @@ void BoundaryPostData::SetUp(json &boundaries)
   }
   std::sort(attributes.begin(), attributes.end());
   attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
+  attributes.shrink_to_fit();
 
   // Cleanup
   postpro->erase("Capacitance");
@@ -1376,6 +1379,7 @@ void BoundaryData::SetUp(json &config)
   attributes.insert(attributes.end(), postpro.attributes.begin(), postpro.attributes.end());
   std::sort(attributes.begin(), attributes.end());
   attributes.erase(unique(attributes.begin(), attributes.end()), attributes.end());
+  attributes.shrink_to_fit();
 
   // Cleanup
   boundaries->erase("PEC");

From 88746228439c860a83eeef38a643bf39a2b41482 Mon Sep 17 00:00:00 2001
From: Hugh Carson <hughcars@amazon.com>
Date: Thu, 4 Jan 2024 18:00:23 -0500
Subject: [PATCH 5/5] Patch for ParNCMesh::Update bug

---
 cmake/ExternalMFEM.cmake                        |  1 +
 extern/patch/mfem/patch_pncmesh_update_fix.diff | 13 +++++++++++++
 2 files changed, 14 insertions(+)
 create mode 100644 extern/patch/mfem/patch_pncmesh_update_fix.diff

diff --git a/cmake/ExternalMFEM.cmake b/cmake/ExternalMFEM.cmake
index 072c1fa70..8019e2167 100644
--- a/cmake/ExternalMFEM.cmake
+++ b/cmake/ExternalMFEM.cmake
@@ -363,6 +363,7 @@ set(MFEM_PATCH_FILES
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_ncmesh_interior_boundary_dev.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mesh_const_fix.diff"
+  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_pncmesh_update_fix.diff"
 )
 
 include(ExternalProject)
diff --git a/extern/patch/mfem/patch_pncmesh_update_fix.diff b/extern/patch/mfem/patch_pncmesh_update_fix.diff
new file mode 100644
index 000000000..2ca6d1d80
--- /dev/null
+++ b/extern/patch/mfem/patch_pncmesh_update_fix.diff
@@ -0,0 +1,13 @@
+diff --git a/mesh/pncmesh.cpp b/mesh/pncmesh.cpp
+index cd6625e9c..a8a88081f 100644
+--- a/mesh/pncmesh.cpp
++++ b/mesh/pncmesh.cpp
+@@ -107,6 +107,8 @@ void ParNCMesh::Update()
+       entity_owner[i].DeleteAll();
+       entity_pmat_group[i].DeleteAll();
+       entity_index_rank[i].DeleteAll();
++      entity_conf_group[i].DeleteAll();
++      entity_elem_local[i].DeleteAll();
+    }
+ 
+    shared_vertices.Clear();