diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 44309faa5a..d4255877e4 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -531,7 +531,7 @@ private:
   void initialize() override;
 
   //! Translate a bin value to an element pointer
-  const libMesh::Elem* get_element_from_bin(int bin) const;
+  const libMesh::Elem& get_element_from_bin(int bin) const;
 
   //! Translate an element pointer to a bin value
   int get_bin_from_element(const libMesh::Elem* elem) const;
diff --git a/src/initialize.cpp b/src/initialize.cpp
index 515442f80f..c6d8250c1e 100644
--- a/src/initialize.cpp
+++ b/src/initialize.cpp
@@ -67,13 +67,18 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
   int n_threads = 1;
 #endif
 
-// initialize libmesh
+// initialize libMesh if it hasn't been initialized already
+// by a call external to OpenMC
 if (!settings::LMI && !libMesh::initialized())
 {
 #ifdef OPENMC_MPI
+  // pass command line args, empty MPI communicator, and number of threads.
+  // Because libMesh was not initialized, we assume that OpenMC is the primary
+  // application and that its main MPI comm should be used.
   settings::LMI = std::make_unique<libMesh::LibMeshInit>(argc, argv, comm, n_threads);
 #else
-  settings::LMI = std::make_unique<libMesh::LibMeshInit>(argc, argv, n_threads);
+  // pass command line args, empty MPI communicator, and number of threads
+  settings::LMI = std::make_unique<libMesh::LibMeshInit>(argc, argv, 0, n_threads);
 #endif
 }
 #endif
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 37f6a3a7ce..ccaf96efd6 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -2099,8 +2099,8 @@ LibMesh::LibMesh(const std::string& filename) {
 
 Position
 LibMesh::centroid(int bin) const {
-  auto elem = this->get_element_from_bin(bin);
-  auto centroid = elem->centroid();
+  auto& elem = this->get_element_from_bin(bin);
+  auto centroid = elem.centroid();
   return {centroid(0), centroid(1), centroid(2)};
 }
 
@@ -2145,7 +2145,7 @@ int LibMesh::n_bins() const {
 int LibMesh::n_surface_bins() const {
   int n_bins = 0;
   for (int i = 0; i < this->n_bins(); i++) {
-    const libMesh::Elem& e = m_->elem_ref(i);
+    const libMesh::Elem& e = get_element_from_bin(i);
     n_bins += e.n_faces();
     // if this is a boundary element, sure to count boundary faces
     // twice
@@ -2269,9 +2269,9 @@ std::pair<std::vector<double>, std::vector<double>>
 LibMesh::plot(Position plot_ll,
               Position plot_ur) const { return {}; }
 
-const libMesh::Elem*
+const libMesh::Elem&
 LibMesh::get_element_from_bin(int bin) const {
-  return m_->elem_ptr(bin);
+  return m_->elem_ref(bin);
 }
 
 double LibMesh::volume(int bin) const {