Fix runtime error when mesh has no boundary elements

palace/models/curlcurloperator.cpp

@@ -20,7 +20,7 @@ namespace
 
 mfem::Array<int> SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh)
 {
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.pec.empty())
   {
     // Check that boundary attributes have been specified correctly.
@@ -92,7 +92,7 @@ CurlCurlOperator::CurlCurlOperator(const IoData &iodata,
   CheckBoundaryProperties();
 
   // Print essential BC information.
-  if (dbc_marker.Max() > 0)
+  if (dbc_marker.Size() && dbc_marker.Max() > 0)
   {
     Mpi::Print("\nConfiguring Dirichlet BC at attributes:\n");
     utils::PrettyPrintMarker(dbc_marker);

palace/models/farfieldboundaryoperator.cpp

@@ -22,7 +22,7 @@ FarfieldBoundaryOperator::FarfieldBoundaryOperator(const IoData &iodata,
   SetUpBoundaryProperties(iodata, mesh);
 
   // Print out BC info for all farfield attributes.
-  if (farfield_marker.Max() > 0)
+  if (farfield_marker.Size() && farfield_marker.Max() > 0)
   {
     Mpi::Print("\nConfiguring Robin absorbing BC (order {:d}) at attributes:\n", order);
     utils::PrettyPrintMarker(farfield_marker);
@@ -33,7 +33,7 @@ void FarfieldBoundaryOperator::SetUpBoundaryProperties(const IoData &iodata,
                                                        const mfem::ParMesh &mesh)
 {
   // Check that impedance boundary attributes have been specified correctly.
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.farfield.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);
@@ -67,7 +67,7 @@ void FarfieldBoundaryOperator::AddDampingBdrCoefficients(double coef,
                                                          SumMatrixCoefficient &fb)
 {
   // First-order absorbing boundary condition.
-  if (farfield_marker.Max() > 0)
+  if (farfield_marker.Size() && farfield_marker.Max() > 0)
   {
     constexpr auto MatType = MaterialPropertyType::INV_Z0;
     constexpr auto ElemType = MeshElementType::BDR_ELEMENT;
@@ -87,7 +87,7 @@ void FarfieldBoundaryOperator::AddExtraSystemBdrCoefficients(double omega,
   // purely imaginary. Multiplying through by μ⁻¹ we get the material coefficient to ω as
   // 1 / (μ √με). Also, this implementation ignores the divergence term ∇⋅Eₜ, as COMSOL
   // does as well.
-  if (farfield_marker.Max() > 0 && order > 1)
+  if (farfield_marker.Size() && farfield_marker.Max() > 0 && order > 1)
   {
     constexpr auto MatType = MaterialPropertyType::INV_PERMEABILITY_C0;
     constexpr auto ElemType = MeshElementType::BDR_ELEMENT;

palace/models/laplaceoperator.cpp

@@ -19,7 +19,7 @@ namespace
 
 mfem::Array<int> SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh)
 {
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.pec.empty() || !iodata.boundaries.lumpedport.empty())
   {
     // Check that boundary attributes have been specified correctly.
@@ -126,7 +126,7 @@ LaplaceOperator::LaplaceOperator(const IoData &iodata,
     source_attr_lists(ConstructSources(iodata))
 {
   // Print essential BC information.
-  if (dbc_marker.Max() > 0)
+  if (dbc_marker.Size() && dbc_marker.Max() > 0)
   {
     Mpi::Print("\nConfiguring Dirichlet BC at attributes:\n");
     utils::PrettyPrintMarker(dbc_marker);

palace/models/lumpedportoperator.cpp

@@ -51,8 +51,10 @@ LumpedPortData::LumpedPortData(const config::LumpedPortData &data,
   for (const auto &elem : data.elements)
   {
     mfem::Array<int> attr_marker;
-    mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Max(), elem.attributes,
-                       attr_marker);
+    mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Size()
+                           ? h1_fespace.GetParMesh()->bdr_attributes.Max()
+                           : 0,
+                       elem.attributes, attr_marker);
     switch (elem.coordinate_system)
     {
       case config::internal::ElementData::CoordinateSystem::CYLINDRICAL:
@@ -290,7 +292,9 @@ void LumpedPortOperator::SetUpBoundaryProperties(const IoData &iodata,
                                                  mfem::ParFiniteElementSpace &h1_fespace)
 {
   // Check that lumped port boundary attributes have been specified correctly.
-  int bdr_attr_max = h1_fespace.GetParMesh()->bdr_attributes.Max();
+  int bdr_attr_max = h1_fespace.GetParMesh()->bdr_attributes.Size()
+                         ? h1_fespace.GetParMesh()->bdr_attributes.Max()
+                         : 0;
   if (!iodata.boundaries.lumpedport.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);

palace/models/spaceoperator.cpp

@@ -22,7 +22,7 @@ namespace
 
 mfem::Array<int> SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMesh &mesh)
 {
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.pec.empty())
   {
     // Check that boundary attributes have been specified correctly.
@@ -100,7 +100,7 @@ SpaceOperator::SpaceOperator(const IoData &iodata,
   CheckBoundaryProperties();
 
   // Print essential BC information.
-  if (dbc_marker.Max() > 0)
+  if (dbc_marker.Size() && dbc_marker.Max() > 0)
   {
     Mpi::Print("\nConfiguring Dirichlet PEC BC at attributes:\n");
     utils::PrettyPrintMarker(dbc_marker);

palace/models/surfaceconductivityoperator.cpp

@@ -25,7 +25,7 @@ void SurfaceConductivityOperator::SetUpBoundaryProperties(const IoData &iodata,
                                                           const mfem::ParMesh &mesh)
 {
   // Check that conductivity boundary attributes have been specified correctly.
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.conductivity.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);
@@ -98,7 +98,7 @@ void SurfaceConductivityOperator::SetUpBoundaryProperties(const IoData &iodata,
 void SurfaceConductivityOperator::PrintBoundaryInfo(const IoData &iodata,
                                                     mfem::ParMesh &mesh)
 {
-  if (conductivity_marker.Max() == 0)
+  if (conductivity_marker.Size() && conductivity_marker.Max() == 0)
   {
     return;
   }
@@ -134,7 +134,7 @@ void SurfaceConductivityOperator::AddExtraSystemBdrCoefficients(double omega,
                                                                 SumMatrixCoefficient &fbr,
                                                                 SumMatrixCoefficient &fbi)
 {
-  if (conductivity_marker.Max() > 0)
+  if (conductivity_marker.Size() && conductivity_marker.Max() > 0)
   {
     // If the provided conductor thickness is empty (zero), prescribe a surface impedance
     // (1+i)/σδ, where δ is the skin depth. If it is nonzero, use a finite thickness

palace/models/surfacecurrentoperator.cpp

@@ -20,8 +20,10 @@ SurfaceCurrentData::SurfaceCurrentData(const config::SurfaceCurrentData &data,
   for (const auto &elem : data.elements)
   {
     mfem::Array<int> attr_marker;
-    mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Max(), elem.attributes,
-                       attr_marker);
+    mesh::AttrToMarker(h1_fespace.GetParMesh()->bdr_attributes.Size()
+                           ? h1_fespace.GetParMesh()->bdr_attributes.Max()
+                           : 0,
+                       elem.attributes, attr_marker);
     switch (elem.coordinate_system)
     {
       case config::internal::ElementData::CoordinateSystem::CYLINDRICAL:
@@ -54,7 +56,9 @@ void SurfaceCurrentOperator::SetUpBoundaryProperties(
     const IoData &iodata, mfem::ParFiniteElementSpace &h1_fespace)
 {
   // Check that surface current boundary attributes have been specified correctly.
-  int bdr_attr_max = h1_fespace.GetParMesh()->bdr_attributes.Max();
+  int bdr_attr_max = h1_fespace.GetParMesh()->bdr_attributes.Size()
+                         ? h1_fespace.GetParMesh()->bdr_attributes.Max()
+                         : 0;
   if (!iodata.boundaries.current.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);

palace/models/surfaceimpedanceoperator.cpp

@@ -23,7 +23,7 @@ void SurfaceImpedanceOperator::SetUpBoundaryProperties(const IoData &iodata,
                                                        const mfem::ParMesh &mesh)
 {
   // Check that impedance boundary attributes have been specified correctly.
-  int bdr_attr_max = mesh.bdr_attributes.Max();
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   if (!iodata.boundaries.impedance.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);
@@ -99,7 +99,7 @@ void SurfaceImpedanceOperator::SetUpBoundaryProperties(const IoData &iodata,
 
 void SurfaceImpedanceOperator::PrintBoundaryInfo(const IoData &iodata, mfem::ParMesh &mesh)
 {
-  if (impedance_marker.Max() == 0)
+  if (impedance_marker.Size() && impedance_marker.Max() == 0)
   {
     return;
   }
@@ -162,7 +162,7 @@ void SurfaceImpedanceOperator::AddStiffnessBdrCoefficients(double coef,
                                                            SumMatrixCoefficient &fb)
 {
   // Lumped inductor boundaries.
-  if (impedance_Ls_marker.Max() > 0)
+  if (impedance_Ls_marker.Size() && impedance_Ls_marker.Max() > 0)
   {
     mfem::Vector v(Z_Lsinv);
     v *= coef;
@@ -174,7 +174,7 @@ void SurfaceImpedanceOperator::AddStiffnessBdrCoefficients(double coef,
 void SurfaceImpedanceOperator::AddMassBdrCoefficients(double coef, SumMatrixCoefficient &fb)
 {
   // Lumped capacitor boundaries.
-  if (impedance_Cs_marker.Max() > 0)
+  if (impedance_Cs_marker.Size() && impedance_Cs_marker.Max() > 0)
   {
     mfem::Vector v(Z_Cs);
     v *= coef;
@@ -186,7 +186,7 @@ void SurfaceImpedanceOperator::AddDampingBdrCoefficients(double coef,
                                                          SumMatrixCoefficient &fb)
 {
   // Lumped resistor boundaries.
-  if (impedance_Rs_marker.Max() > 0)
+  if (impedance_Rs_marker.Size() && impedance_Rs_marker.Max() > 0)
   {
     mfem::Vector v(Z_Rsinv);
     v *= coef;

palace/models/surfacepostoperator.cpp

@@ -72,8 +72,8 @@ SurfacePostOperator::InterfaceDielectricData::InterfaceDielectricData(
     }
 
     // Store markers for this element of the postprocessing boundary.
-    mesh::AttrToMarker(mesh.bdr_attributes.Max(), elem.attributes,
-                       attr_markers.emplace_back());
+    mesh::AttrToMarker(mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0,
+                       elem.attributes, attr_markers.emplace_back());
   }
 }
 
@@ -103,8 +103,8 @@ SurfacePostOperator::InterfaceDielectricData::GetCoefficient(
 SurfacePostOperator::SurfaceChargeData::SurfaceChargeData(
     const config::CapacitanceData &data, mfem::ParMesh &mesh)
 {
-  mesh::AttrToMarker(mesh.bdr_attributes.Max(), data.attributes,
-                     attr_markers.emplace_back());
+  mesh::AttrToMarker(mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0,
+                     data.attributes, attr_markers.emplace_back());
 }
 
 std::unique_ptr<mfem::Coefficient> SurfacePostOperator::SurfaceChargeData::GetCoefficient(
@@ -124,8 +124,8 @@ SurfacePostOperator::SurfaceFluxData::SurfaceFluxData(const config::InductanceDa
 
   // Construct the coefficient for this postprocessing boundary (copies the direction
   // vector).
-  mesh::AttrToMarker(mesh.bdr_attributes.Max(), data.attributes,
-                     attr_markers.emplace_back());
+  mesh::AttrToMarker(mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0,
+                     data.attributes, attr_markers.emplace_back());
 }
 
 std::unique_ptr<mfem::Coefficient>

palace/models/waveportoperator.cpp

@@ -507,7 +507,10 @@ WavePortData::WavePortData(const config::WavePortData &data, const MaterialOpera
   {
     attr_list.Append(attr);
   }
-  mesh::AttrToMarker(nd_fespace.GetParMesh()->bdr_attributes.Max(), attr_list, attr_marker);
+  mesh::AttrToMarker(nd_fespace.GetParMesh()->bdr_attributes.Size()
+                         ? nd_fespace.GetParMesh()->bdr_attributes.Max()
+                         : 0,
+                     attr_list, attr_marker);
   port_mesh = std::make_unique<mfem::ParSubMesh>(
       mfem::ParSubMesh::CreateFromBoundary(mesh, attr_list));
 
@@ -955,7 +958,9 @@ void WavePortOperator::SetUpBoundaryProperties(const IoData &iodata,
                                                mfem::ParFiniteElementSpace &h1_fespace)
 {
   // Check that wave port boundary attributes have been specified correctly.
-  int bdr_attr_max = nd_fespace.GetParMesh()->bdr_attributes.Max();
+  int bdr_attr_max = nd_fespace.GetParMesh()->bdr_attributes.Size()
+                         ? nd_fespace.GetParMesh()->bdr_attributes.Max()
+                         : 0;
   if (!iodata.boundaries.waveport.empty())
   {
     mfem::Array<int> bdr_attr_marker(bdr_attr_max);

palace/utils/geodata.cpp

@@ -1098,6 +1098,7 @@ std::unique_ptr<mfem::Mesh> LoadMesh(const std::string &path, bool remove_curvat
       MFEM_ABORT("Unable to open translated mesh file \"" << tmp << "\"!");
     }
 #endif
+
     mesh = std::make_unique<mfem::Mesh>(fi, 1, 1, true);
   }
   else
@@ -1216,8 +1217,10 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
   MFEM_VERIFY(orig_mesh->Dimension() == 3 && !orig_mesh->Nonconforming(),
               "Nonconforming or 2D meshes have not been tested yet!");
   mfem::Array<int> mat_marker, bdr_marker;
-  GetUsedAttributeMarkers(iodata, orig_mesh->attributes.Max(),
-                          orig_mesh->bdr_attributes.Max(), 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);
   bool warn = false;
   for (int be = 0; be < orig_mesh->GetNBE(); be++)
   {
@@ -1394,7 +1397,7 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
     }
     if (new_nbdr > new_nbdr_step2)
     {
-      Mpi::Print("Added boundary elements for subdomain interfaces to the mesh\n",
+      Mpi::Print("Added {:d} boundary elements for subdomain interfaces to the mesh\n",
                  new_nbdr - new_nbdr_step2);
     }
   }
@@ -1454,7 +1457,8 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
     // 1-based, some boundary attributes may be empty since they were removed from the
     // original mesh, but to keep indices the same as config file we don't compact the
     // list.
-    int max_bdr_attr = orig_mesh->bdr_attributes.Max();
+    int max_bdr_attr =
+        orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0;
     for (int f = 0; f < orig_mesh->GetNumFaces(); f++)
     {
       if (add_bdr_faces[f] > 0)
@@ -1485,7 +1489,8 @@ std::map<int, std::array<int, 2>> CheckMesh(std::unique_ptr<mfem::Mesh> &orig_me
           b = 0;
         }
         MFEM_VERIFY(a + b > 0, "Invalid new boundary element attribute!");
-        int new_attr = max_bdr_attr + (a * (a - 1)) / 2 + b;  // At least max_bdr_attr + 1
+        int new_attr = max_bdr_attr +
+                       (b > 0 ? (a * (a - 1)) / 2 + b : a);  // At least max_bdr_attr + 1
         if (new_attr_map.find(new_attr) == new_attr_map.end())
         {
           new_attr_map.emplace(new_attr, std::array<int, 2>{a, b});