diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index c6a870f8d..f6908be43 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -123,6 +123,9 @@ public:
   //! \param[in] bin Mesh bin to generate a label for
   virtual std::string bin_label(int bin) const = 0;
 
+  //! Return the mesh type
+  virtual std::string get_mesh_type() const = 0;
+
   // Data members
   int id_ {-1};     //!< User-specified ID
   int n_dimension_; //!< Number of dimensions
@@ -135,9 +138,12 @@ public:
   virtual ~StructuredMesh() = default;
 
   using MeshIndex = std::array<int, 3>;
+
   struct MeshDistance {
-    MeshDistance() {};
-    MeshDistance(int _index, bool _maxSurface, double _distance): nextIndex{_index}, maxSurface(_maxSurface), distance{_distance} {};
+    MeshDistance() = default;
+    MeshDistance(int _index, bool _maxSurface, double _distance): 
+      nextIndex{_index}, maxSurface{_maxSurface}, distance{_distance}
+      { }
     int nextIndex { -1 };
     double distance { INFTY };
     bool maxSurface { true };
@@ -185,13 +191,13 @@ public:
   //
   //! \param[in] r Position to get indices for
   //! \param[out] in_mesh Whether position is in mesh
-  //! \return ijk Array of mesh indices
+  //! \return Array of mesh indices
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
 
   //! Get mesh indices corresponding to a mesh bin
   //
   //! \param[in] bin Mesh bin
-  //! \param[out] ijk Mesh indices
+  //! \return ijk Mesh indices
   virtual MeshIndex get_indices_from_bin(int bin) const;
 
   //! Get mesh index in a particular direction
@@ -200,19 +206,30 @@ public:
   //! \param[in] i Direction index
   virtual int get_index_in_direction(double r, int i) const = 0;
 
-  //! Get the closest distance from to coordinate r to the grid for the mesh grid boundary in the negative direction
+  //! Get the closest distance from the coordinate r to the grid surface  
+  //! in i direction  that bounds mesh cell ijk and that is larger than l
+  //! The coordinate r does not have to be inside the mesh cell ijk. In 
+  //! curved coordinates, multiple crossings of the same surface can happen,
+  //! these are selected by the parameter l 
   //!
   //! \param[in] ijk Array of mesh indices
-  //! \param[in] i Direction index  
+  //! \param[in] i direction index of grid surface
+  //! \param[in] r0 position, from where to calculate the distance
+  //! \param[in] u direction of flight. actual position is r0 + l * u
+  //! \param[in] l actual chord length
+  //! \return MeshDistance struct with closest distance, next cell index in i-direction and min/max surface indicator
   virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0;
 
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
+  //! Get shape as xt::xtensor
+  xt::xtensor<int, 1> get_x_shape() const;
+
   // Data members
   xt::xtensor<double, 1> lower_left_;  //!< Lower-left coordinates of mesh
   xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
-  xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
+  std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
 
 protected:
 
@@ -231,7 +248,9 @@ public:
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
 
+  virtual std::string get_mesh_type() const override;
 
+  static const std::string mesh_type;
 
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
 
@@ -274,6 +293,10 @@ public:
 
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
+
+  virtual std::string get_mesh_type() const override;
+
+  static const std::string mesh_type;
     
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
 
@@ -296,7 +319,7 @@ public:
   double negative_grid_boundary(const MeshIndex& ijk, int i) const;
 
 
-  vector<vector<double>> grid_;
+  array<vector<double>, 3> grid_;
 
   int set_grid();
 };
@@ -312,6 +335,10 @@ public:
 
   int get_index_in_direction(double r, int i) const override;
 
+  virtual std::string get_mesh_type() const override;
+
+  static const std::string mesh_type;
+
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
 
   std::pair<vector<double>, vector<double>> plot(
@@ -319,7 +346,7 @@ public:
 
   void to_hdf5(hid_t group) const override;
 
-  vector<vector<double>> grid_;
+  array<vector<double>, 3> grid_;
 
   int set_grid();
 
@@ -358,6 +385,10 @@ public:
 
   int get_index_in_direction(double r, int i) const override;
 
+  virtual std::string get_mesh_type() const override;
+
+  static const std::string mesh_type;
+
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
 
   std::pair<vector<double>, vector<double>> plot(
@@ -365,7 +396,7 @@ public:
 
   void to_hdf5(hid_t group) const override;
 
-  vector<vector<double>> grid_;
+  array<vector<double>, 3> grid_;
 
   int set_grid();
 
@@ -405,6 +436,8 @@ public:
   UnstructuredMesh(pugi::xml_node node);
   UnstructuredMesh(const std::string& filename);
 
+  static const std::string mesh_type;
+
   // Overridden Methods
   void surface_bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins) const override;
@@ -476,6 +509,8 @@ public:
   MOABMesh(const std::string& filename, double length_multiplier = 1.0);
   MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
 
+  static const std::string mesh_lib_type;
+
   // Overridden Methods
 
   void bins_crossed(Position r0, Position r1, const Direction& u,
@@ -623,6 +658,8 @@ public:
   LibMesh(pugi::xml_node node);
   LibMesh(const std::string& filename, double length_multiplier = 1.0);
 
+  static const std::string mesh_lib_type;
+
   // Overridden Methods
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 82545e2bf..43f66d7a1 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -471,8 +471,8 @@ class SphericalMesh(Mesh):
 _MESH_TYPE_MAP = {
     'regular': RegularMesh,
     'rectilinear': RectilinearMesh,
-    'cylindrical': RectilinearMesh,
-    'spherical': RectilinearMesh
+    'cylindrical': CylindricalMesh,
+    'spherical': SphericalMesh
 }
 
 
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index b35281f4a..5584210a1 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -577,7 +577,7 @@ void ufs_count_sites()
 
 #ifdef OPENMC_MPI
     // Send source fraction to all processors
-    int n_bins = xt::prod(simulation::ufs_mesh->shape_)();
+    int n_bins = simulation::ufs_mesh->n_bins();
     MPI_Bcast(
       simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
 #endif
diff --git a/src/mesh.cpp b/src/mesh.cpp
index d4948fb86..bd139e94b 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -32,8 +32,6 @@
 #include "openmc/tallies/tally.h"
 #include "openmc/xml_interface.h"
 
-#include <iostream>
-
 #ifdef LIBMESH
 #include "libmesh/mesh_tools.h"
 #include "libmesh/numeric_vector.h"
@@ -158,6 +156,10 @@ std::string StructuredMesh::bin_label(int bin) const
   }
 }
 
+xt::xtensor<int, 1> StructuredMesh::get_x_shape() const {
+  return xt::adapt(shape_, { n_dimension_});
+}
+
 //==============================================================================
 // Unstructured Mesh implementation
 //==============================================================================
@@ -168,7 +170,7 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   // check the mesh type
   if (check_for_node(node, "type")) {
     auto temp = get_node_value(node, "type", true, true);
-    if (temp != "unstructured") {
+    if (temp != mesh_type) {
       fatal_error(fmt::format("Invalid mesh type: {}", temp));
     }
   }
@@ -197,6 +199,8 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   }
 }
 
+const std::string UnstructuredMesh::mesh_type = "unstructured";
+
 void UnstructuredMesh::surface_bins_crossed(
   Position r0, Position r1, const Direction& u, vector<int>& bins) const
 {
@@ -212,7 +216,7 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
 {
   hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_));
 
-  write_dataset(mesh_group, "type", "unstructured");
+  write_dataset(mesh_group, "type", mesh_type);
   write_dataset(mesh_group, "filename", filename_);
   write_dataset(mesh_group, "library", this->library());
   // write volume of each element
@@ -298,7 +302,7 @@ int StructuredMesh::get_bin(Position r) const
 
 int StructuredMesh::n_bins() const
 {
-  return xt::prod(shape_)();
+  return std::accumulate(shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
 }
 
 int StructuredMesh::n_surface_bins() const
@@ -428,7 +432,7 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u,
       ijk[k] = distances[k].nextIndex;
       distances[k] = distance_to_grid_boundary(ijk, k, r0, u, l);
 
-      // Check if we have left the interiour of the mesh 
+      // Check if we have left the interior of the mesh 
       in_mesh = ((ijk[k]>=1) and (ijk[k]<=shape_[k]));
 
       // If we are still inside the mesh, tally inward current for the next cell
@@ -524,14 +528,15 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
     fatal_error("Must specify <dimension> on a regular mesh.");
   }
 
-  shape_ = get_node_xarray<int>(node, "dimension");
-  int n = n_dimension_ = shape_.size();
+  xt::xtensor<int, 1> shape = get_node_xarray<int>(node, "dimension");
+  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());
 
   // Check that dimensions are all greater than zero
-  if (xt::any(shape_ <= 0)) {
+  if (xt::any(shape <= 0)) {
     fatal_error("All entries on the <dimension> element for a tally "
                 "mesh must be positive.");
   }
@@ -545,7 +550,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 
   // Make sure lower_left and dimension match
-  if (shape_.size() != lower_left_.size()) {
+  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>.");
   }
@@ -571,7 +576,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
     }
 
     // Set width and upper right coordinate
-    upper_right_ = xt::eval(lower_left_ + shape_ * width_);
+    upper_right_ = xt::eval(lower_left_ + shape * width_);
 
   } else if (check_for_node(node, "upper_right")) {
     upper_right_ = get_node_xarray<double>(node, "upper_right");
@@ -590,13 +595,13 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
     }
 
     // Set width
-    width_ = xt::eval((upper_right_ - lower_left_) / shape_);
+    width_ = xt::eval((upper_right_ - lower_left_) / shape);
   } else {
     fatal_error("Must specify either <upper_right> or <width> on a mesh.");
   }
 
   // Set volume fraction
-  volume_frac_ = 1.0 / xt::prod(shape_)();
+  volume_frac_ = 1.0 / xt::prod(shape)();
 }
 
 int RegularMesh::get_index_in_direction(double r, int i) const
@@ -604,6 +609,12 @@ int RegularMesh::get_index_in_direction(double r, int i) const
   return std::ceil((r - lower_left_[i]) / width_[i]);
 }
 
+const std::string RegularMesh::mesh_type = "regular";
+
+std::string RegularMesh::get_mesh_type() const {
+  return mesh_type;
+}
+
 double RegularMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const
 {
   return lower_left_[i] + ijk[i] * width_[i];
@@ -679,7 +690,7 @@ void RegularMesh::to_hdf5(hid_t group) const
   hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
 
   write_dataset(mesh_group, "type", "regular");
-  write_dataset(mesh_group, "dimension", shape_);
+  write_dataset(mesh_group, "dimension", get_x_shape());
   write_dataset(mesh_group, "lower_left", lower_left_);
   write_dataset(mesh_group, "upper_right", upper_right_);
   write_dataset(mesh_group, "width", width_);
@@ -752,7 +763,6 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
 {
   n_dimension_ = 3;
 
-  grid_.resize(n_dimension_);
   grid_[0] = get_node_array<double>(node, "x_grid");
   grid_[1] = get_node_array<double>(node, "y_grid");
   grid_[2] = get_node_array<double>(node, "z_grid");
@@ -762,6 +772,12 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+const std::string RectilinearMesh::mesh_type = "rectilinear";
+
+std::string RectilinearMesh::get_mesh_type() const {
+  return mesh_type;
+}
+
 double RectilinearMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const
 {
   return grid_[i][ijk[i]];
@@ -802,12 +818,10 @@ int RectilinearMesh::set_grid()
                  "must each have at least 2 points");
       return OPENMC_E_INVALID_ARGUMENT;
     }
-    for (int i = 1; i < g.size(); ++i) {
-      if (g[i] <= g[i - 1]) {
-        set_errmsg("Values in for x-, y-, and z- grids for "
-                   "rectilinear meshes must be sorted and unique.");
-        return OPENMC_E_INVALID_ARGUMENT;
-      }
+    if (not std::is_sorted(g.begin(), g.end(), std::less_equal<double>())) {
+      set_errmsg("Values in for x-, y-, and z- grids for "
+                 "rectilinear meshes must be sorted and unique.");
+      return OPENMC_E_INVALID_ARGUMENT;
     }
   }
 
@@ -817,8 +831,6 @@ int RectilinearMesh::set_grid()
   return 0;
 }
 
-
-
 int RectilinearMesh::get_index_in_direction(double r, int i) const
 {
   return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
@@ -874,7 +886,6 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node}
 {
   n_dimension_ = 3;
 
-  grid_.resize(n_dimension_);
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "p_grid");
   grid_[2] = get_node_array<double>(node, "z_grid");
@@ -884,6 +895,12 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+const std::string CylindricalMesh::mesh_type = "cylindrical";
+
+std::string CylindricalMesh::get_mesh_type() const {
+  return mesh_type;
+}
+
 StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh) const
 {
   Position mapped_r;
@@ -1020,12 +1037,10 @@ int CylindricalMesh::set_grid()
                  "must each have at least 2 points");
       return OPENMC_E_INVALID_ARGUMENT;
     }
-    for (int i = 1; i < g.size(); ++i) {
-      if (g[i] <= g[i - 1]) {
-        set_errmsg("Values in for r-, phi-, and z- grids for "
-                   "cylindrical meshes must be sorted and unique.");
-        return OPENMC_E_INVALID_ARGUMENT;
-      }
+    if (not std::is_sorted(g.begin(), g.end(), std::less_equal<double>())) {
+      set_errmsg("Values in for r-, phi-, and z- grids for "
+                 "cylindrical meshes must be sorted and unique.");
+      return OPENMC_E_INVALID_ARGUMENT;
     }
   }
   if (grid_[0].front() < 0.0) {
@@ -1093,7 +1108,6 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node}
 {
   n_dimension_ = 3;
 
-  grid_.resize(n_dimension_);
   grid_[0] = get_node_array<double>(node, "r_grid");
   grid_[1] = get_node_array<double>(node, "t_grid");
   grid_[2] = get_node_array<double>(node, "p_grid");
@@ -1103,6 +1117,12 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+const std::string SphericalMesh::mesh_type = "spherical";
+
+std::string SphericalMesh::get_mesh_type() const {
+  return mesh_type;
+}
+
 StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) const
 {
   Position mapped_r;
@@ -1128,8 +1148,6 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh)
 double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const 
 {
   
-  ////std::cout << "    r: " << shell << " " << grid_[0].size() << " " << shape_[0] << "\n";
-
   if ((shell < 0) || (shell >= shape_[0])) return INFTY;
   
   // solve |r+s*u| = r0
@@ -1275,12 +1293,10 @@ int SphericalMesh::set_grid()
                  "must each have at least 2 points");
       return OPENMC_E_INVALID_ARGUMENT;
     }
-    for (int i = 1; i < g.size(); ++i) {
-      if (g[i] <= g[i - 1]) {
-        set_errmsg("Values in for r-, theta-, and phi- grids for "
-                   "spherical meshes must be sorted and unique.");
-        return OPENMC_E_INVALID_ARGUMENT;
-      }
+    if (not std::is_sorted(g.begin(), g.end(), std::less_equal<double>())) {
+      set_errmsg("Values in for r-, theta-, and phi- grids for "
+                 "spherical meshes must be sorted and unique.");
+      return OPENMC_E_INVALID_ARGUMENT;
     }
     if (g.front() < 0.0) {
         set_errmsg("r-, theta-, and phi- grids for "
@@ -1333,7 +1349,7 @@ void SphericalMesh::to_hdf5(hid_t group) const
 {
   hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
 
-  write_dataset(mesh_group, "type", "spherical");
+  write_dataset(mesh_group, "type", SphericalMesh::mesh_type);
   write_dataset(mesh_group, "r_grid", grid_[0]);
   write_dataset(mesh_group, "t_grid", grid_[1]);
   write_dataset(mesh_group, "p_grid", grid_[2]);
@@ -1388,15 +1404,8 @@ extern "C" int openmc_mesh_get_type(int32_t index, char* type)
   if (int err = check_mesh(index))
     return err;
 
-  if (is_mesh_type<RegularMesh>(index)) {
-    std::strcpy(type, "regular");
-  } else if (is_mesh_type<RectilinearMesh>(index)) {
-    std::strcpy(type, "rectilinear");
-  } else if (is_mesh_type<CylindricalMesh>(index)) {
-    std::strcpy(type, "cylindrical");
-  } else if (is_mesh_type<SphericalMesh>(index)) {
-    std::strcpy(type, "spherical");
-  }
+  std::strcpy(type, model::meshes[index].get()->get_mesh_type().c_str());
+
   return 0;
 }
 
@@ -1409,13 +1418,14 @@ extern "C" int openmc_extend_meshes(
   std::string mesh_type;
 
   for (int i = 0; i < n; ++i) {
-    if (std::strcmp(type, "regular") == 0) {
+    //if (std::strcmp(type, RegularMesh::mesh_type.c_str()) == 0) {
+    if (RegularMesh::mesh_type == type)  {
       model::meshes.push_back(make_unique<RegularMesh>());
-    } else if (std::strcmp(type, "rectilinear") == 0) {
+    } else if (RectilinearMesh::mesh_type == type) {
       model::meshes.push_back(make_unique<RectilinearMesh>());
-    } else if (std::strcmp(type, "cylindrical") == 0) {
+    } else if (CylindricalMesh::mesh_type == type) {
       model::meshes.push_back(make_unique<CylindricalMesh>());
-    } else if (std::strcmp(type, "spherical") == 0) {
+    } else if (SphericalMesh::mesh_type == type) {
       model::meshes.push_back(make_unique<SphericalMesh>());
     } else {
       throw std::runtime_error {"Unknown mesh type: " + std::string(type)};
@@ -1436,14 +1446,14 @@ extern "C" int openmc_add_unstructured_mesh(
   bool valid_lib = false;
 
 #ifdef DAGMC
-  if (lib_name == "moab") {
+  if (lib_name == MOABMesh::mesh_lib_type) {
     model::meshes.push_back(std::move(make_unique<MOABMesh>(mesh_file)));
     valid_lib = true;
   }
 #endif
 
 #ifdef LIBMESH
-  if (lib_name == "libmesh") {
+  if (lib_name == LibMesh::mesh_lib_type) {
     model::meshes.push_back(std::move(make_unique<LibMesh>(mesh_file)));
     valid_lib = true;
   }
@@ -1515,9 +1525,8 @@ extern "C" int openmc_regular_mesh_set_dimension(
   RegularMesh* mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
 
   // Copy dimension
-  vector<std::size_t> shape = {static_cast<std::size_t>(n)};
-  mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
-  mesh->n_dimension_ = mesh->shape_.size();
+  mesh->n_dimension_ = n;
+  std::copy(dims, dims + n, mesh->shape_.begin());
   return 0;
 }
 
@@ -1553,15 +1562,15 @@ extern "C" int openmc_regular_mesh_set_params(
   if (ll && ur) {
     m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
     m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
-    m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_;
+    m->width_ = (m->upper_right_ - m->lower_left_) / m->get_x_shape();
   } else if (ll && width) {
     m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
     m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
-    m->upper_right_ = m->lower_left_ + m->shape_ * m->width_;
+    m->upper_right_ = m->lower_left_ + m->get_x_shape() * m->width_;
   } else if (ur && width) {
     m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
     m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
-    m->lower_left_ = m->upper_right_ - m->shape_ * m->width_;
+    m->lower_left_ = m->upper_right_ - m->get_x_shape() * m->width_;
   } else {
     set_errmsg("At least two parameters must be specified.");
     return OPENMC_E_INVALID_ARGUMENT;
@@ -1570,42 +1579,23 @@ extern "C" int openmc_regular_mesh_set_params(
   return 0;
 }
 
-//! Get the rectilinear mesh grid
-extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
-  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
-{
-  if (int err = check_mesh_type<RectilinearMesh>(index))
-    return err;
-  RectilinearMesh* m =
-    dynamic_cast<RectilinearMesh*>(model::meshes[index].get());
 
-  if (m->lower_left_.dimension() == 0) {
-    set_errmsg("Mesh parameters have not been set.");
-    return OPENMC_E_ALLOCATE;
-  }
-
-  *grid_x = m->grid_[0].data();
-  *nx = m->grid_[0].size();
-  *grid_y = m->grid_[1].data();
-  *ny = m->grid_[1].size();
-  *grid_z = m->grid_[2].data();
-  *nz = m->grid_[2].size();
-
-  return 0;
-}
-
-//! Set the rectilienar mesh parameters
-extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index,
+//! Set the mesh parameters for rectilinear, cylindrical and spharical meshes
+template <class C>
+int openmc_structured_mesh_set_grid_impl(int32_t index,
   const double* grid_x, const int nx, const double* grid_y, const int ny,
   const double* grid_z, const int nz)
 {
-  if (int err = check_mesh_type<RectilinearMesh>(index))
+  if (int err = check_mesh_type<C>(index))
     return err;
-  RectilinearMesh* m =
-    dynamic_cast<RectilinearMesh*>(model::meshes[index].get());
+
+  C* m = dynamic_cast<C*>(model::meshes[index].get());
 
   m->n_dimension_ = 3;
-  m->grid_.resize(m->n_dimension_);
+
+  m->grid_[0].reserve(nx);
+  m->grid_[1].reserve(ny);
+  m->grid_[2].reserve(nz);
 
   for (int i = 0; i < nx; i++) {
     m->grid_[0].push_back(grid_x[i]);
@@ -1618,17 +1608,17 @@ extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index,
   }
 
   int err = m->set_grid();
-  return err;
+  return err;  
 }
 
-//! Get the cylindrical mesh grid
-extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x,
+//! Get the mesh parameters for rectilinear, cylindrical and spharical meshes
+template <class C>
+int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x,
   int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
 {
-  if (int err = check_mesh_type<CylindricalMesh>(index))
+  if (int err = check_mesh_type<C>(index))
     return err;
-  CylindricalMesh* m =
-    dynamic_cast<CylindricalMesh*>(model::meshes[index].get());
+  C* m = dynamic_cast<C*>(model::meshes[index].get());
 
   if (m->lower_left_.dimension() == 0) {
     set_errmsg("Mesh parameters have not been set.");
@@ -1645,55 +1635,44 @@ extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x,
   return 0;
 }
 
+
+//! Get the rectilinear mesh grid
+extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
+  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
+{
+  return openmc_structured_mesh_get_grid_impl<RectilinearMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
+
+
+//! Set the rectilienar mesh parameters
+extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index,
+  const double* grid_x, const int nx, const double* grid_y, const int ny,
+  const double* grid_z, const int nz)
+{
+  return openmc_structured_mesh_set_grid_impl<RectilinearMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
+
+//! Get the cylindrical mesh grid
+extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x,
+  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
+{
+  return openmc_structured_mesh_get_grid_impl<CylindricalMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
+
 //! Set the cylindrical mesh parameters
 extern "C" int openmc_cylindrical_mesh_set_grid(int32_t index,
   const double* grid_x, const int nx, const double* grid_y, const int ny,
   const double* grid_z, const int nz)
 {
-  if (int err = check_mesh_type<CylindricalMesh>(index))
-    return err;
-  CylindricalMesh* m =
-    dynamic_cast<CylindricalMesh*>(model::meshes[index].get());
-
-  m->n_dimension_ = 3;
-  m->grid_.resize(m->n_dimension_);
-
-  for (int i = 0; i < nx; i++) {
-    m->grid_[0].push_back(grid_x[i]);
-  }
-  for (int i = 0; i < ny; i++) {
-    m->grid_[1].push_back(grid_y[i]);
-  }
-  for (int i = 0; i < nz; i++) {
-    m->grid_[2].push_back(grid_z[i]);
-  }
-
-  int err = m->set_grid();
-  return err;
+  return openmc_structured_mesh_set_grid_impl<CylindricalMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);
 }
 
 //! Get the spherical mesh grid
 extern "C" int openmc_spherical_mesh_get_grid(int32_t index, double** grid_x,
   int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
 {
-  if (int err = check_mesh_type<SphericalMesh>(index))
-    return err;
-  SphericalMesh* m =
-    dynamic_cast<SphericalMesh*>(model::meshes[index].get());
 
-  if (m->lower_left_.dimension() == 0) {
-    set_errmsg("Mesh parameters have not been set.");
-    return OPENMC_E_ALLOCATE;
-  }
-
-  *grid_x = m->grid_[0].data();
-  *nx = m->grid_[0].size();
-  *grid_y = m->grid_[1].data();
-  *ny = m->grid_[1].size();
-  *grid_z = m->grid_[2].data();
-  *nz = m->grid_[2].size();
-
-  return 0;
+  return openmc_structured_mesh_get_grid_impl<SphericalMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);;
 }
 
 //! Set the spherical mesh parameters
@@ -1701,30 +1680,13 @@ extern "C" int openmc_spherical_mesh_set_grid(int32_t index,
   const double* grid_x, const int nx, const double* grid_y, const int ny,
   const double* grid_z, const int nz)
 {
-  if (int err = check_mesh_type<SphericalMesh>(index))
-    return err;
-  SphericalMesh* m =
-    dynamic_cast<SphericalMesh*>(model::meshes[index].get());
-
-  m->n_dimension_ = 3;
-  m->grid_.resize(m->n_dimension_);
-
-  for (int i = 0; i < nx; i++) {
-    m->grid_[0].push_back(grid_x[i]);
-  }
-  for (int i = 0; i < ny; i++) {
-    m->grid_[1].push_back(grid_y[i]);
-  }
-  for (int i = 0; i < nz; i++) {
-    m->grid_[2].push_back(grid_z[i]);
-  }
-
-  int err = m->set_grid();
-  return err;
+  return openmc_structured_mesh_set_grid_impl<SphericalMesh>(index, grid_x, nx, grid_y, ny, grid_z, nz);
 }
 
 #ifdef DAGMC
 
+const std::string MOABMesh::mesh_lib_type = "moab";
+
 MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
 {
   initialize();
@@ -1994,7 +1956,7 @@ double MOABMesh::volume(int bin) const
 
 std::string MOABMesh::library() const
 {
-  return "moab";
+  return MOABMesh::mesh_lib_type;
 }
 
 double MOABMesh::tet_volume(moab::EntityHandle tet) const
@@ -2313,6 +2275,8 @@ void MOABMesh::write(const std::string& base_filename) const
 
 #ifdef LIBMESH
 
+const std::string LibMesh::mesh_lib_type = "libmesh";
+
 LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
 {
   initialize();
@@ -2387,7 +2351,7 @@ Position LibMesh::centroid(int bin) const
 
 std::string LibMesh::library() const
 {
-  return "libmesh";
+  return LibMesh::mesh_lib_type;
 }
 
 int LibMesh::n_bins() const
@@ -2564,23 +2528,23 @@ void read_meshes(pugi::xml_node root)
     }
 
     // Read mesh and add to vector
-    if (mesh_type == "regular") {
+    if (mesh_type == RegularMesh::mesh_type) {
       model::meshes.push_back(make_unique<RegularMesh>(node));
-    } else if (mesh_type == "rectilinear") {
+    } else if (mesh_type == RectilinearMesh::mesh_type) {
       model::meshes.push_back(make_unique<RectilinearMesh>(node));
-    } else if (mesh_type == "cylindrical") {
+    } else if (mesh_type == CylindricalMesh::mesh_type) {
       model::meshes.push_back(make_unique<CylindricalMesh>(node));
-    } else if (mesh_type == "spherical") {
+    } else if (mesh_type == SphericalMesh::mesh_type) {
       model::meshes.push_back(make_unique<SphericalMesh>(node));
 #ifdef DAGMC
-    } else if (mesh_type == "unstructured" && mesh_lib == "moab") {
+    } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_lib == MOABMesh::mesh_lib_type) {
       model::meshes.push_back(make_unique<MOABMesh>(node));
 #endif
 #ifdef LIBMESH
-    } else if (mesh_type == "unstructured" && mesh_lib == "libmesh") {
+    } 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 == "unstructured") {
+    } else if (mesh_type == UnstructuredMesh::mesh_type) {
       fatal_error("Unstructured mesh support is not enabled or the mesh "
                   "library is invalid.");
     } else {
diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat
index 60b3f1928..4247f6ed1 100644
--- a/tests/regression_tests/filter_mesh/results_true.dat
+++ b/tests/regression_tests/filter_mesh/results_true.dat
@@ -1 +1 @@
-79f82f287ff427b896ab8aedc81f4671949571316b90e89806f0753ed0628dd5ac32f8d0e4a3f5f5f5c65a9e8a06a6869eb5fc896b4bb423870877afca2a7d4f
\ No newline at end of file
+6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c