diff --git a/src/Omega_h_gmsh.cpp b/src/Omega_h_gmsh.cpp
index 74191c17d..8d2de612d 100644
--- a/src/Omega_h_gmsh.cpp
+++ b/src/Omega_h_gmsh.cpp
@@ -144,7 +144,7 @@ static void read(
 }
 
 static void read_internal_entities_section(Mesh& mesh, Real format,
-    std::vector<std::string>& physical_names, std::istream& stream,
+    std::unordered_map<int, std::string>& physical_names, std::istream& stream,
     bool is_binary, bool needs_swapping) {
   Int num_points, num_curves, num_surfaces, num_volumes;
   read(stream, num_points, is_binary, needs_swapping);
@@ -171,8 +171,9 @@ static void read_internal_entities_section(Mesh& mesh, Real format,
       read(stream, physical, is_binary, needs_swapping);
       OMEGA_H_CHECK(physical != 0);
       if (physical > 0) {
-        const auto& physicalname = physical_names[physical - 1];
-        mesh.class_sets[physicalname].emplace_back(0, tag);
+        const auto& physicalname = physical_names.find(physical);
+        OMEGA_H_CHECK(physicalname != physical_names.end());
+        mesh.class_sets[physicalname->second].emplace_back(0, tag);
       }
     }
   }
@@ -198,8 +199,9 @@ static void read_internal_entities_section(Mesh& mesh, Real format,
         read(stream, physical, is_binary, needs_swapping);
         OMEGA_H_CHECK(physical != 0);
         if (physical > 0) {
-          const auto& physical_name = physical_names[physical - 1];
-          mesh.class_sets[physical_name].emplace_back(dim, tag);
+          auto physical_name = physical_names.find(physical);
+          OMEGA_H_CHECK(physical_name != physical_names.end());
+          mesh.class_sets[physical_name->second].emplace_back(dim, tag);
         }
       }
       Int num_bounding_points;
@@ -232,20 +234,18 @@ void read_internal(std::istream& stream, Mesh* mesh) {
     }
   }
   OMEGA_H_CHECK(data_size == sizeof(Real));
-  std::vector<std::string> physical_names;
+  std::unordered_map<int, std::string> physical_names;
   if (seek_optional_section(stream, "$PhysicalNames")) {
     Int num_physicals;
     read(stream, num_physicals, is_binary, needs_swapping);
-    physical_names.reserve(static_cast<std::size_t>(num_physicals));
     eat_newlines(stream);
     for (auto i = 0; i < num_physicals; ++i) {
       Int dim, number;
       read(stream, dim, is_binary, needs_swapping);
       read(stream, number, is_binary, needs_swapping);
-      OMEGA_H_CHECK(number == i + 1);
       std::string name;
       stream >> name;
-      physical_names.push_back(name.substr(1, name.size() - 2));
+      physical_names.emplace(number, name.substr(1, name.size() - 2));
     }
   }
   if (seek_optional_section(stream, "$Entities")) {
@@ -465,8 +465,9 @@ void read_internal(std::istream& stream, Mesh* mesh) {
         const auto entity = pair.first;
         const auto physical = pair.second;
         std::string physical_name(std::to_string(physical));
-        if (physical <= static_cast<Int>(physical_names.size())) {
-          physical_name = physical_names[physical - 1];
+        auto physical_name_it = physical_names.find(physical);
+        if (physical_name_it != physical_names.end()) {
+          physical_name = physical_name_it->second;
         }
         mesh->class_sets[physical_name].emplace_back(dim, entity);
       }
diff --git a/src/unit_io.cpp b/src/unit_io.cpp
index aba6700fe..5fe4ceed3 100644
--- a/src/unit_io.cpp
+++ b/src/unit_io.cpp
@@ -34,6 +34,36 @@ Mesh 2;
 )GMSH";
 #endif
 
+static const char* GMSH_PHYSICAL_ENTITIES_MSH2 = R"GMSH(
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+4
+2 13 "patch1"
+2 14 "patchInBetween"
+3 15 "comp1"
+3 16 "comp2"
+$EndPhysicalNames
+$Nodes
+6
+1 0 0 0
+2 1 0 0
+3 0 1 0
+4 0 0 1
+5 -1 0 0
+6 0 -1 0
+$EndNodes
+$Elements
+5
+1 2 2 14 2 1 4 2
+2 2 2 13 4 1 5 3
+3 4 2 15 1 2 4 3 1
+4 4 2 15 1 5 3 4 1
+5 4 2 16 2 1 2 4 6
+$EndElements
+)GMSH";
+
 static const char* GMSH_SQUARE_MSH2 = R"GMSH(
 $MeshFormat
 2.2 0 8
@@ -1393,6 +1423,26 @@ static void test_gmsh(Library* lib) {
       OMEGA_H_CHECK(mesh.class_sets.empty());
     }
   }
+  {
+    if (lib->world()->rank() == 0) {
+      {
+        std::ofstream oss("physical-entities.msh");
+        oss << GMSH_PHYSICAL_ENTITIES_MSH2;
+      }
+      lib->world()->barrier();
+      auto mesh = Omega_h::gmsh::read("physical-entities.msh", lib->world());
+      const auto num_class_sets = mesh.class_sets.size();
+      OMEGA_H_CHECK(num_class_sets == 4);
+      OMEGA_H_CHECK(mesh.class_sets["patch1"].size() == 1);
+      OMEGA_H_CHECK(mesh.class_sets["patch1"][0].id == 4);
+      OMEGA_H_CHECK(mesh.class_sets["patchInBetween"].size() == 1);
+      OMEGA_H_CHECK(mesh.class_sets["patchInBetween"][0].id == 2);
+      OMEGA_H_CHECK(mesh.class_sets["comp1"].size() == 1);
+      OMEGA_H_CHECK(mesh.class_sets["comp1"][0].id == 1);
+      OMEGA_H_CHECK(mesh.class_sets["comp2"].size() == 1);
+      OMEGA_H_CHECK(mesh.class_sets["comp2"][0].id == 2);
+    }
+  }
   {
     const std::vector<const char*> meshes{
         GMSH_PHYSICAL_MSH2, GMSH_PHYSICAL_MSH40, GMSH_PHYSICAL_MSH41};