diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index a56705c9ec5..5c9272e93b9 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -132,8 +132,14 @@ class Mesh {
   // Constructors and destructor
   Mesh() = default;
   Mesh(pugi::xml_node node);
+  Mesh(hid_t group);
   virtual ~Mesh() = default;
 
+  // Factory method for creating meshes from either an XML node or HDF5 group
+  template<typename T>
+  static const std::unique_ptr<Mesh>& create(
+    T dataset, const std::string& mesh_type, const std::string& mesh_library);
+
   // Methods
   //! Perform any preparation needed to support point location within the mesh
   virtual void prepare_for_point_location() {};
@@ -258,6 +264,7 @@ class StructuredMesh : public Mesh {
 public:
   StructuredMesh() = default;
   StructuredMesh(pugi::xml_node node) : Mesh {node} {};
+  StructuredMesh(hid_t group) : Mesh {group} {};
   virtual ~StructuredMesh() = default;
 
   using MeshIndex = std::array<int, 3>;
@@ -423,6 +430,7 @@ class PeriodicStructuredMesh : public StructuredMesh {
 public:
   PeriodicStructuredMesh() = default;
   PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
+  PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {};
 
   Position local_coords(const Position& r) const override
   {
@@ -442,6 +450,7 @@ class RegularMesh : public StructuredMesh {
   // Constructors
   RegularMesh() = default;
   RegularMesh(pugi::xml_node node);
+  RegularMesh(hid_t group);
 
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
@@ -481,6 +490,8 @@ class RegularMesh : public StructuredMesh {
   //! Return the volume for a given mesh index
   double volume(const MeshIndex& ijk) const override;
 
+  int set_grid();
+
   // Data members
   double volume_frac_;           //!< Volume fraction of each mesh element
   double element_volume_;        //!< Volume of each mesh element
@@ -492,6 +503,7 @@ class RectilinearMesh : public StructuredMesh {
   // Constructors
   RectilinearMesh() = default;
   RectilinearMesh(pugi::xml_node node);
+  RectilinearMesh(hid_t group);
 
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
@@ -534,6 +546,7 @@ class CylindricalMesh : public PeriodicStructuredMesh {
   // Constructors
   CylindricalMesh() = default;
   CylindricalMesh(pugi::xml_node node);
+  CylindricalMesh(hid_t group);
 
   // Overridden methods
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@@ -598,6 +611,7 @@ class SphericalMesh : public PeriodicStructuredMesh {
   // Constructors
   SphericalMesh() = default;
   SphericalMesh(pugi::xml_node node);
+  SphericalMesh(hid_t group);
 
   // Overridden methods
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@@ -668,6 +682,7 @@ class UnstructuredMesh : public Mesh {
   // Constructors
   UnstructuredMesh() { n_dimension_ = 3; };
   UnstructuredMesh(pugi::xml_node node);
+  UnstructuredMesh(hid_t group);
 
   static const std::string mesh_type;
   virtual std::string get_mesh_type() const override;
@@ -774,6 +789,7 @@ class MOABMesh : public UnstructuredMesh {
   // Constructors
   MOABMesh() = default;
   MOABMesh(pugi::xml_node);
+  MOABMesh(hid_t group);
   MOABMesh(const std::string& filename, double length_multiplier = 1.0);
   MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
 
@@ -943,6 +959,7 @@ class LibMesh : public UnstructuredMesh {
 public:
   // Constructors
   LibMesh(pugi::xml_node node);
+  LibMesh(hid_t group);
   LibMesh(const std::string& filename, double length_multiplier = 1.0);
   LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
 
@@ -1069,6 +1086,11 @@ class AdaptiveLibMesh : public LibMesh {
 //! \param[in] root XML node
 void read_meshes(pugi::xml_node root);
 
+//! Read meshes from an HDF5 file
+//
+//! \param[in] group HDF5 group ("meshes" group)
+void read_meshes(hid_t group);
+
 //! Write mesh data to an HDF5 group
 //
 //! \param[in] group HDF5 group
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 58d218b9cad..610d057cf08 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -230,6 +230,42 @@ void MaterialVolumes::add_volume_unsafe(
 // Mesh implementation
 //==============================================================================
 
+template<typename T>
+const std::unique_ptr<Mesh>& Mesh::create(
+  T dataset, const std::string& mesh_type, const std::string& mesh_library)
+{
+  // Determine mesh type. Add to model vector and map
+  if (mesh_type == RegularMesh::mesh_type) {
+    model::meshes.push_back(make_unique<RegularMesh>(dataset));
+  } else if (mesh_type == RectilinearMesh::mesh_type) {
+    model::meshes.push_back(make_unique<RectilinearMesh>(dataset));
+  } else if (mesh_type == CylindricalMesh::mesh_type) {
+    model::meshes.push_back(make_unique<CylindricalMesh>(dataset));
+  } else if (mesh_type == SphericalMesh::mesh_type) {
+    model::meshes.push_back(make_unique<SphericalMesh>(dataset));
+#ifdef OPENMC_DAGMC_ENABLED
+  } else if (mesh_type == UnstructuredMesh::mesh_type &&
+             mesh_library == MOABMesh::mesh_lib_type) {
+    model::meshes.push_back(make_unique<MOABMesh>(dataset));
+#endif
+#ifdef OPENMC_LIBMESH_ENABLED
+  } else if (mesh_type == UnstructuredMesh::mesh_type &&
+             mesh_library == LibMesh::mesh_lib_type) {
+    model::meshes.push_back(make_unique<LibMesh>(dataset));
+#endif
+  } else if (mesh_type == UnstructuredMesh::mesh_type) {
+    fatal_error("Unstructured mesh support is not enabled or the mesh "
+                "library is invalid.");
+  } else {
+    fatal_error(fmt::format("Invalid mesh type: {}", mesh_type));
+  }
+
+  // Map ID to position in vector
+  model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
+
+  return model::meshes.back();
+}
+
 Mesh::Mesh(pugi::xml_node node)
 {
   // Read mesh id
@@ -238,6 +274,17 @@ Mesh::Mesh(pugi::xml_node node)
     name_ = get_node_value(node, "name");
 }
 
+Mesh::Mesh(hid_t group)
+{
+  // Read mesh ID
+  read_attribute(group, "id", id_);
+
+  // Read mesh name
+  if (object_exists(group, "name")) {
+    read_dataset(group, "name", name_);
+  }
+}
+
 void Mesh::set_id(int32_t id)
 {
   assert(id >= 0 || id == C_NONE);
@@ -265,7 +312,13 @@ void Mesh::set_id(int32_t id)
 
   // Update ID and entry in the mesh map
   id_ = id;
-  model::mesh_map[id] = model::meshes.size() - 1;
+
+  // find the index of this mesh in the model::meshes vector
+  // (search in reverse because this mesh was likely just added to the vector)
+  auto it = std::find_if(model::meshes.rbegin(), model::meshes.rend(),
+    [this](const std::unique_ptr<Mesh>& mesh) { return mesh.get() == this; });
+
+  model::mesh_map[id] = std::distance(model::meshes.begin(), it.base()) - 1;
 }
 
 vector<double> Mesh::volumes() const
@@ -627,6 +680,46 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   }
 }
 
+UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group)
+{
+  n_dimension_ = 3;
+
+  // check the mesh type
+  if (object_exists(group, "type")) {
+    std::string temp;
+    read_dataset(group, "type", temp);
+    if (temp != mesh_type) {
+      fatal_error(fmt::format("Invalid mesh type: {}", temp));
+    }
+  }
+
+  // check if a length unit multiplier was specified
+  if (object_exists(group, "length_multiplier")) {
+    read_dataset(group, "length_multiplier", length_multiplier_);
+  }
+
+  // get the filename of the unstructured mesh to load
+  if (object_exists(group, "filename")) {
+    read_dataset(group, "filename", filename_);
+    if (!file_exists(filename_)) {
+      fatal_error("Mesh file '" + filename_ + "' does not exist!");
+    }
+  } else {
+    fatal_error(fmt::format(
+      "No filename supplied for unstructured mesh with ID: {}", id_));
+  }
+
+  if (attribute_exists(group, "options")) {
+    read_attribute(group, "options", options_);
+  }
+
+  // check if mesh tally data should be written with
+  // statepoint files
+  if (attribute_exists(group, "output")) {
+    read_attribute(group, "output", output_);
+  }
+}
+
 void UnstructuredMesh::determine_bounds()
 {
   double xmin = INFTY;
@@ -1086,6 +1179,72 @@ void StructuredMesh::surface_bins_crossed(
 // RegularMesh implementation
 //==============================================================================
 
+int RegularMesh::set_grid()
+{
+  auto shape = xt::adapt(shape_, {n_dimension_});
+
+  // Check that dimensions are all greater than zero
+  if (xt::any(shape <= 0)) {
+    set_errmsg("All entries for a regular mesh dimensions "
+               "must be positive.");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  // Make sure lower_left and dimension match
+  if (lower_left_.size() != n_dimension_) {
+    set_errmsg("Number of entries in lower_left must be the same "
+               "as the regular mesh dimensions.");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+  if (width_.size() > 0) {
+
+    // Check to ensure width has same dimensions
+    if (width_.size() != n_dimension_) {
+      set_errmsg("Number of entries on width must be the same as "
+                 "the regular mesh dimensions.");
+      return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // Check for negative widths
+    if (xt::any(width_ < 0.0)) {
+      set_errmsg("Cannot have a negative width on a regular mesh.");
+      return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // Set width and upper right coordinate
+    upper_right_ = xt::eval(lower_left_ + shape * width_);
+
+  } else if (upper_right_.size() > 0) {
+
+    // Check to ensure upper_right_ has same dimensions
+    if (upper_right_.size() != n_dimension_) {
+      set_errmsg("Number of entries on upper_right must be the "
+                 "same as the regular mesh dimensions.");
+      return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // Check that upper-right is above lower-left
+    if (xt::any(upper_right_ < lower_left_)) {
+      set_errmsg(
+        "The upper_right coordinates of a regular mesh must be greater than "
+        "the lower_left coordinates.");
+      return OPENMC_E_INVALID_ARGUMENT;
+    }
+
+    // Set width
+    width_ = xt::eval((upper_right_ - lower_left_) / shape);
+  }
+
+  // Set material volumes
+  volume_frac_ = 1.0 / xt::prod(shape)();
+
+  element_volume_ = 1.0;
+  for (int i = 0; i < n_dimension_; i++) {
+    element_volume_ *= width_[i];
+  }
+  return 0;
+}
+
 RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
 {
   // Determine number of dimensions for mesh
@@ -1100,12 +1259,6 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
   }
   std::copy(shape.begin(), shape.end(), shape_.begin());
 
-  // Check that dimensions are all greater than zero
-  if (xt::any(shape <= 0)) {
-    fatal_error("All entries on the <dimension> element for a tally "
-                "mesh must be positive.");
-  }
-
   // Check for lower-left coordinates
   if (check_for_node(node, "lower_left")) {
     // Read mesh lower-left corner location
@@ -1114,12 +1267,6 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
     fatal_error("Must specify <lower_left> on a mesh.");
   }
 
-  // Make sure lower_left and dimension match
-  if (shape.size() != lower_left_.size()) {
-    fatal_error("Number of entries on <lower_left> must be the same "
-                "as the number of entries on <dimension>.");
-  }
-
   if (check_for_node(node, "width")) {
     // Make sure one of upper-right or width were specified
     if (check_for_node(node, "upper_right")) {
@@ -1128,49 +1275,52 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
 
     width_ = get_node_xarray<double>(node, "width");
 
-    // Check to ensure width has same dimensions
-    auto n = width_.size();
-    if (n != lower_left_.size()) {
-      fatal_error("Number of entries on <width> must be the same as "
-                  "the number of entries on <lower_left>.");
-    }
+  } else if (check_for_node(node, "upper_right")) {
 
-    // Check for negative widths
-    if (xt::any(width_ < 0.0)) {
-      fatal_error("Cannot have a negative <width> on a tally mesh.");
-    }
+    upper_right_ = get_node_xarray<double>(node, "upper_right");
 
-    // Set width and upper right coordinate
-    upper_right_ = xt::eval(lower_left_ + shape * width_);
+  } else {
+    fatal_error("Must specify either <upper_right> or <width> on a mesh.");
+  }
 
-  } else if (check_for_node(node, "upper_right")) {
-    upper_right_ = get_node_xarray<double>(node, "upper_right");
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
 
-    // Check to ensure width has same dimensions
-    auto n = upper_right_.size();
-    if (n != lower_left_.size()) {
-      fatal_error("Number of entries on <upper_right> must be the "
-                  "same as the number of entries on <lower_left>.");
-    }
+RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
+{
+  // Determine number of dimensions for mesh
+  if (!object_exists(group, "dimension")) {
+    fatal_error("Must specify <dimension> on a regular mesh.");
+  }
 
-    // Check that upper-right is above lower-left
-    if (xt::any(upper_right_ < lower_left_)) {
-      fatal_error("The <upper_right> coordinates must be greater than "
-                  "the <lower_left> coordinates on a tally mesh.");
-    }
+  xt::xtensor<int, 1> shape;
+  read_dataset(group, "dimension", shape);
+  int n = n_dimension_ = shape.size();
+  if (n != 1 && n != 2 && n != 3) {
+    fatal_error("Mesh must be one, two, or three dimensions.");
+  }
+  std::copy(shape.begin(), shape.end(), shape_.begin());
 
-    // Set width
-    width_ = xt::eval((upper_right_ - lower_left_) / shape);
+  // Check for lower-left coordinates
+  if (object_exists(group, "lower_left")) {
+    // Read mesh lower-left corner location
+    read_dataset(group, "lower_left", lower_left_);
   } else {
-    fatal_error("Must specify either <upper_right> or <width> on a mesh.");
+    fatal_error("Must specify lower_left dataset on a mesh.");
   }
 
-  // Set material volumes
-  volume_frac_ = 1.0 / xt::prod(shape)();
+  if (object_exists(group, "upper_right")) {
 
-  element_volume_ = 1.0;
-  for (int i = 0; i < n_dimension_; i++) {
-    element_volume_ *= width_[i];
+    read_dataset(group, "upper_right", upper_right_);
+
+  } else {
+    fatal_error("Must specify either upper_right dataset on a mesh.");
+  }
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
   }
 }
 
@@ -1343,6 +1493,19 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+RectilinearMesh::RectilinearMesh(hid_t group) : StructuredMesh {group}
+{
+  n_dimension_ = 3;
+
+  read_dataset(group, "x_grid", grid_[0]);
+  read_dataset(group, "y_grid", grid_[1]);
+  read_dataset(group, "z_grid", grid_[2]);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string RectilinearMesh::mesh_type = "rectilinear";
 
 std::string RectilinearMesh::get_mesh_type() const
@@ -1478,6 +1641,19 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node)
   }
 }
 
+CylindricalMesh::CylindricalMesh(hid_t group) : PeriodicStructuredMesh {group}
+{
+  n_dimension_ = 3;
+  read_dataset(group, "r_grid", grid_[0]);
+  read_dataset(group, "phi_grid", grid_[1]);
+  read_dataset(group, "z_grid", grid_[2]);
+  read_dataset(group, "origin", origin_);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string CylindricalMesh::mesh_type = "cylindrical";
 
 std::string CylindricalMesh::get_mesh_type() const
@@ -1756,6 +1932,20 @@ SphericalMesh::SphericalMesh(pugi::xml_node node)
   }
 }
 
+SphericalMesh::SphericalMesh(hid_t group) : PeriodicStructuredMesh {group}
+{
+  n_dimension_ = 3;
+
+  read_dataset(group, "r_grid", grid_[0]);
+  read_dataset(group, "theta_grid", grid_[1]);
+  read_dataset(group, "phi_grid", grid_[2]);
+  read_dataset(group, "origin", origin_);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string SphericalMesh::mesh_type = "spherical";
 
 std::string SphericalMesh::get_mesh_type() const
@@ -2520,6 +2710,11 @@ MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
+MOABMesh::MOABMesh(hid_t group) : UnstructuredMesh(group)
+{
+  initialize();
+}
+
 MOABMesh::MOABMesh(const std::string& filename, double length_multiplier)
   : UnstructuredMesh()
 {
@@ -3228,6 +3423,15 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
+LibMesh::LibMesh(hid_t group) : UnstructuredMesh(group)
+{
+  // filename_ and length_multiplier_ will already be set by the
+  // UnstructuredMesh constructor
+  set_mesh_pointer_from_filename(filename_);
+  set_length_multiplier(length_multiplier_);
+  initialize();
+}
+
 // create the mesh from a pointer to a libMesh Mesh
 LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
 {
@@ -3618,34 +3822,51 @@ void read_meshes(pugi::xml_node root)
       mesh_lib = get_node_value(node, "library", true, true);
     }
 
-    // Read mesh and add to vector
-    if (mesh_type == RegularMesh::mesh_type) {
-      model::meshes.push_back(make_unique<RegularMesh>(node));
-    } else if (mesh_type == RectilinearMesh::mesh_type) {
-      model::meshes.push_back(make_unique<RectilinearMesh>(node));
-    } else if (mesh_type == CylindricalMesh::mesh_type) {
-      model::meshes.push_back(make_unique<CylindricalMesh>(node));
-    } else if (mesh_type == SphericalMesh::mesh_type) {
-      model::meshes.push_back(make_unique<SphericalMesh>(node));
-#ifdef OPENMC_DAGMC_ENABLED
-    } else if (mesh_type == UnstructuredMesh::mesh_type &&
-               mesh_lib == MOABMesh::mesh_lib_type) {
-      model::meshes.push_back(make_unique<MOABMesh>(node));
-#endif
-#ifdef OPENMC_LIBMESH_ENABLED
-    } else if (mesh_type == UnstructuredMesh::mesh_type &&
-               mesh_lib == LibMesh::mesh_lib_type) {
-      model::meshes.push_back(make_unique<LibMesh>(node));
-#endif
-    } else if (mesh_type == UnstructuredMesh::mesh_type) {
-      fatal_error("Unstructured mesh support is not enabled or the mesh "
-                  "library is invalid.");
+    Mesh::create(node, mesh_type, mesh_lib);
+  }
+}
+
+void read_meshes(hid_t group)
+{
+  std::unordered_set<int> mesh_ids;
+
+  std::vector<int> ids;
+  read_attribute(group, "ids", ids);
+
+  for (auto id : ids) {
+
+    // Check to make sure multiple meshes in the same file don't share IDs
+    if (contains(mesh_ids, id)) {
+      fatal_error(fmt::format("Two or more meshes use the same unique ID "
+                              "'{}' in the same HDF5 input file",
+        id));
+    }
+    mesh_ids.insert(id);
+
+    // If we've already read a mesh with the same ID in a *different* file,
+    // assume it is the same here
+    if (model::mesh_map.find(id) != model::mesh_map.end()) {
+      warning(fmt::format("Mesh with ID={} appears in multiple files.", id));
+      continue;
+    }
+
+    std::string name = fmt::format("mesh {}", id);
+    hid_t mesh_group = open_group(group, name.c_str());
+
+    std::string mesh_type;
+    if (object_exists(mesh_group, "type")) {
+      read_dataset(mesh_group, "type", mesh_type);
     } else {
-      fatal_error("Invalid mesh type: " + mesh_type);
+      mesh_type = "regular";
+    }
+
+    // determine the mesh library to use
+    std::string mesh_lib;
+    if (object_exists(mesh_group, "library")) {
+      read_dataset(mesh_group, "library", mesh_lib);
     }
 
-    // Map ID to position in vector
-    model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
+    Mesh::create(mesh_group, mesh_type, mesh_lib);
   }
 }
 
diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp
index 26762ad18a6..0d648d3335e 100644
--- a/src/weight_windows.cpp
+++ b/src/weight_windows.cpp
@@ -1328,6 +1328,10 @@ extern "C" int openmc_weight_windows_import(const char* filename)
 
   hid_t weight_windows_group = open_group(ww_file, "weight_windows");
 
+  hid_t mesh_group = open_group(ww_file, "meshes");
+
+  read_meshes(mesh_group);
+
   std::vector<std::string> names = group_names(weight_windows_group);
 
   for (const auto& name : names) {
diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt
index 8fedc2daa57..5f87db9eac2 100644
--- a/tests/cpp_unit_tests/CMakeLists.txt
+++ b/tests/cpp_unit_tests/CMakeLists.txt
@@ -5,6 +5,7 @@ set(TEST_NAMES
   test_interpolate
   test_math
   test_mcpl_stat_sum
+  test_mesh
   # Add additional unit test files here
 )
 
diff --git a/tests/cpp_unit_tests/test_mesh.cpp b/tests/cpp_unit_tests/test_mesh.cpp
new file mode 100644
index 00000000000..24c4f77373f
--- /dev/null
+++ b/tests/cpp_unit_tests/test_mesh.cpp
@@ -0,0 +1,257 @@
+#include <cstdio>
+#include <iostream>
+#include <string>
+
+#include <catch2/catch_test_macros.hpp>
+#include <pugixml.hpp>
+
+#include "openmc/hdf5_interface.h"
+#include "openmc/mesh.h"
+
+using namespace openmc;
+
+TEST_CASE("Test mesh hdf5 roundtrip - regular")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1">
+            <dimension>3 4 5</dimension>
+            <lower_left>-2 -3 -5</lower_left>
+            <upper_right>2 3 5</upper_right>
+       </mesh>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = RegularMesh(root);
+
+  hid_t file_id = file_open("mesh.h5", 'w');
+
+  mesh.to_hdf5(file_id);
+
+  file_close(file_id);
+
+  hid_t file_id2 = file_open("mesh.h5", 'r');
+
+  hid_t group = open_group(file_id2, "mesh 1");
+
+  auto mesh2 = RegularMesh(group);
+
+  file_close(file_id2);
+
+  remove("mesh.h5");
+
+  REQUIRE(mesh2.shape_ == mesh.shape_);
+
+  REQUIRE(mesh2.lower_left() == mesh.lower_left());
+
+  REQUIRE(mesh2.upper_right() == mesh.upper_right());
+}
+
+TEST_CASE("Test mesh hdf5 roundtrip - rectilinear")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1" type="rectilinear">
+            <x_grid>0.0 1.0 5.0 10.0</x_grid>
+            <y_grid>-10.0 -5.0 0.0</y_grid>
+            <z_grid>-100.0 0.0 100.0</z_grid>
+        </mesh>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = RectilinearMesh(root);
+
+  hid_t file_id = file_open("mesh.h5", 'w');
+
+  mesh.to_hdf5(file_id);
+
+  file_close(file_id);
+
+  hid_t file_id2 = file_open("mesh.h5", 'r');
+
+  hid_t group = open_group(file_id2, "mesh 1");
+
+  auto mesh2 = RectilinearMesh(group);
+
+  file_close(file_id2);
+
+  remove("mesh.h5");
+
+  REQUIRE(mesh2.shape_ == mesh.shape_);
+
+  REQUIRE(mesh2.grid_ == mesh.grid_);
+}
+
+TEST_CASE("Test mesh hdf5 roundtrip - cylindrical")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1" type="cylindrical">
+            <r_grid>0.1 0.2 0.5 1.0</r_grid>
+            <phi_grid>0.0 6.283185307179586</phi_grid>
+            <z_grid>0.1 0.2 0.4 0.6 1.0</z_grid>
+            <origin>0 0 0</origin>
+        </mesh>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = CylindricalMesh(root);
+
+  hid_t file_id = file_open("mesh.h5", 'w');
+
+  mesh.to_hdf5(file_id);
+
+  file_close(file_id);
+
+  hid_t file_id2 = file_open("mesh.h5", 'r');
+
+  hid_t group = open_group(file_id2, "mesh 1");
+
+  auto mesh2 = CylindricalMesh(group);
+
+  file_close(file_id2);
+
+  remove("mesh.h5");
+
+  REQUIRE(mesh2.shape_ == mesh.shape_);
+
+  REQUIRE(mesh2.grid_ == mesh.grid_);
+}
+
+TEST_CASE("Test mesh hdf5 roundtrip - spherical")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1" type="spherical">
+            <r_grid>0.1 0.2 0.5 1.0</r_grid>
+            <theta_grid>0.0 3.141592653589793</theta_grid>
+            <phi_grid>0.0 6.283185307179586</phi_grid>
+            <origin>0.0 0.0 0.0</origin>
+        </mesh>'
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = SphericalMesh(root);
+
+  hid_t file_id = file_open("mesh.h5", 'w');
+
+  mesh.to_hdf5(file_id);
+
+  file_close(file_id);
+
+  hid_t file_id2 = file_open("mesh.h5", 'r');
+
+  hid_t group = open_group(file_id2, "mesh 1");
+
+  auto mesh2 = SphericalMesh(group);
+
+  file_close(file_id2);
+
+  remove("mesh.h5");
+
+  REQUIRE(mesh2.shape_ == mesh.shape_);
+
+  REQUIRE(mesh2.grid_ == mesh.grid_);
+}
+
+TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <meshes>
+        <mesh id="1" type="spherical">
+            <r_grid>0.1 0.2 0.5 1.0</r_grid>
+            <theta_grid>0.0 3.141592653589793</theta_grid>
+            <phi_grid>0.0 6.283185307179586</phi_grid>
+            <origin>0.0 0.0 0.0</origin>
+        </mesh>
+         <mesh id="2">
+            <dimension>3 4 5</dimension>
+            <lower_left>-2 -3 -5</lower_left>
+            <upper_right>2 3 5</upper_right>
+       </mesh>
+       </meshes>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("meshes");
+
+  read_meshes(root);
+
+  const auto spherical_mesh_xml =
+    dynamic_cast<SphericalMesh*>(model::meshes[0].get());
+  const auto regular_mesh_xml =
+    dynamic_cast<RegularMesh*>(model::meshes[1].get());
+
+  hid_t file_id = file_open("meshes.h5", 'w');
+
+  hid_t root_group = create_group(file_id, "root");
+
+  open_group(file_id, "root");
+
+  meshes_to_hdf5(root_group);
+
+  close_group(root_group);
+
+  file_close(file_id);
+
+  hid_t file_id2 = file_open("meshes.h5", 'r');
+
+  hid_t root_group_read = open_group(file_id2, "root");
+
+  hid_t mesh_group_read = open_group(root_group_read, "meshes");
+
+  read_meshes(mesh_group_read);
+
+  // increment mesh IDs to avoid collision during read
+  for (auto& mesh : model::meshes) {
+    mesh->set_id(mesh->id() + 10);
+  }
+
+  const auto spherical_mesh_hdf5 = dynamic_cast<SphericalMesh*>(
+    model::meshes[model::mesh_map[spherical_mesh_xml->id_]].get());
+  const auto regular_mesh_hdf5 = dynamic_cast<RegularMesh*>(
+    model::meshes[model::mesh_map[regular_mesh_xml->id_]].get());
+
+  remove("meshes.h5");
+
+  REQUIRE(spherical_mesh_hdf5->shape_ == spherical_mesh_xml->shape_);
+  REQUIRE(spherical_mesh_hdf5->grid_ == spherical_mesh_xml->grid_);
+
+  REQUIRE(regular_mesh_hdf5->shape_ == regular_mesh_xml->shape_);
+  REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left());
+  REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right());
+}