Merge pull request #154 from awslabs/sjg/geodata-updates

Cleanup in `geodata.hpp/cpp`
awslabs · Jan 5, 2024 · 9b9cb88 · 9b9cb88
2 parents bc2b658 + 9ff441d
commit 9b9cb88
Show file tree

Hide file tree

Showing 17 changed files with 338 additions and 318 deletions.
diff --git 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
@@ -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();
diff --git 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
@@ -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())
@@ -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/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
@@ -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)
   {
@@ -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/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
@@ -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())
   {
@@ -125,16 +126,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/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
@@ -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)
   {
@@ -109,8 +110,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 +147,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/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
@@ -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())
@@ -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/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);