From 109be21bd96fffe8d686af31e0a002762f125cca Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Thu, 24 Aug 2023 18:28:15 -0700 Subject: [PATCH] Fix runtime error when mesh has no boundary elements --- palace/models/curlcurloperator.cpp | 4 ++-- palace/models/farfieldboundaryoperator.cpp | 8 ++++---- palace/models/laplaceoperator.cpp | 4 ++-- palace/models/lumpedportoperator.cpp | 10 +++++++--- palace/models/spaceoperator.cpp | 4 ++-- palace/models/surfaceconductivityoperator.cpp | 6 +++--- palace/models/surfacecurrentoperator.cpp | 10 +++++++--- palace/models/surfaceimpedanceoperator.cpp | 10 +++++----- palace/models/surfacepostoperator.cpp | 12 ++++++------ palace/models/waveportoperator.cpp | 9 +++++++-- palace/utils/geodata.cpp | 15 ++++++++++----- 11 files changed, 55 insertions(+), 37 deletions(-) diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp index ba00d3a63..b29323293 100644 --- a/palace/models/curlcurloperator.cpp +++ b/palace/models/curlcurloperator.cpp @@ -20,7 +20,7 @@ namespace mfem::Array 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); diff --git a/palace/models/farfieldboundaryoperator.cpp b/palace/models/farfieldboundaryoperator.cpp index 12881a3d8..3a8b5e624 100644 --- a/palace/models/farfieldboundaryoperator.cpp +++ b/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 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; diff --git a/palace/models/laplaceoperator.cpp b/palace/models/laplaceoperator.cpp index 146f5cd94..b4015316c 100644 --- a/palace/models/laplaceoperator.cpp +++ b/palace/models/laplaceoperator.cpp @@ -19,7 +19,7 @@ namespace mfem::Array 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); diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp index e706ee007..eedf5a0e2 100644 --- a/palace/models/lumpedportoperator.cpp +++ b/palace/models/lumpedportoperator.cpp @@ -51,8 +51,10 @@ LumpedPortData::LumpedPortData(const config::LumpedPortData &data, for (const auto &elem : data.elements) { mfem::Array 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 bdr_attr_marker(bdr_attr_max); diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index 22641ee56..7860df797 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -22,7 +22,7 @@ namespace mfem::Array 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); diff --git a/palace/models/surfaceconductivityoperator.cpp b/palace/models/surfaceconductivityoperator.cpp index 0072bcd73..705315f67 100644 --- a/palace/models/surfaceconductivityoperator.cpp +++ b/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 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 diff --git a/palace/models/surfacecurrentoperator.cpp b/palace/models/surfacecurrentoperator.cpp index 84e680b8a..bbcd1d17b 100644 --- a/palace/models/surfacecurrentoperator.cpp +++ b/palace/models/surfacecurrentoperator.cpp @@ -20,8 +20,10 @@ SurfaceCurrentData::SurfaceCurrentData(const config::SurfaceCurrentData &data, for (const auto &elem : data.elements) { mfem::Array 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 bdr_attr_marker(bdr_attr_max); diff --git a/palace/models/surfaceimpedanceoperator.cpp b/palace/models/surfaceimpedanceoperator.cpp index 709ac2f28..4671cc6d0 100644 --- a/palace/models/surfaceimpedanceoperator.cpp +++ b/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 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; diff --git a/palace/models/surfacepostoperator.cpp b/palace/models/surfacepostoperator.cpp index 452b84e65..d1bd57896 100644 --- a/palace/models/surfacepostoperator.cpp +++ b/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 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 diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp index f1aef5e91..9ef00d249 100644 --- a/palace/models/waveportoperator.cpp +++ b/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::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 bdr_attr_marker(bdr_attr_max); diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index ab7cbc72d..fbb2b7735 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1098,6 +1098,7 @@ std::unique_ptr LoadMesh(const std::string &path, bool remove_curvat MFEM_ABORT("Unable to open translated mesh file \"" << tmp << "\"!"); } #endif + mesh = std::make_unique(fi, 1, 1, true); } else @@ -1216,8 +1217,10 @@ std::map> CheckMesh(std::unique_ptr &orig_me MFEM_VERIFY(orig_mesh->Dimension() == 3 && !orig_mesh->Nonconforming(), "Nonconforming or 2D meshes have not been tested yet!"); mfem::Array 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> CheckMesh(std::unique_ptr &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> CheckMesh(std::unique_ptr &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> CheckMesh(std::unique_ptr &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{a, b});