diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 6a2995687cbff74907766d694b5443ed95551b9a..e555ca224c1d02f11e11696542774f2e470b6432 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -70,7 +70,6 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -80,7 +79,16 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
@@ -95,8 +103,17 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -157,8 +174,17 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
@@ -229,7 +255,15 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
       face_list[j] = ref_face_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_face_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -239,7 +273,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list, is_boundary});
+    descriptor.addRefItemList(RefFaceList{ref_id, face_list, is_boundary});
   }
 
   Array<bool> is_boundary_node(descriptor.node_number_vector.size());
@@ -265,7 +299,15 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -274,7 +316,7 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -372,8 +414,17 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
       cell_list[j] = ref_cell_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
-    descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list, false));
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_cell_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_cell_list.first};
+      }
+    }();
+
+    descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
   }
 
   ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
@@ -481,7 +532,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
         face_list[j] = ref_face_list.second[j];
       }
-      const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+
+      auto ref_id = [&] {
+        if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_face_list.first);
+            i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+          return i_physical_ref->second.refId();
+        } else {
+          return RefId{ref_face_list.first};
+        }
+      }();
 
       bool is_boundary = true;
       for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -490,8 +549,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
           break;
         }
       }
-
-      descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list, is_boundary});
+      descriptor.addRefItemList(RefFaceList(ref_id, face_list, is_boundary));
     }
   }
 
@@ -566,7 +624,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       for (size_t j = 0; j < ref_edge_list.second.size(); ++j) {
         edge_list[j] = ref_edge_list.second[j];
       }
-      const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_edge_list.first);
+
+      auto ref_id = [&] {
+        if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_edge_list.first);
+            i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+          return i_physical_ref->second.refId();
+        } else {
+          return RefId{ref_edge_list.first};
+        }
+      }();
 
       bool is_boundary = true;
       for (size_t i_node = 0; i_node < edge_list.size(); ++i_node) {
@@ -575,7 +641,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
         }
       }
 
-      descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list, is_boundary});
+      descriptor.addRefItemList(RefEdgeList(ref_id, edge_list, is_boundary));
     }
   }
 
@@ -602,7 +668,15 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
     for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
       point_list[j] = ref_point_list.second[j];
     }
-    const GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_point_list.first);
+
+    auto ref_id = [&] {
+      if (auto i_physical_ref = gmsh_data.m_physical_ref_map.find(ref_point_list.first);
+          i_physical_ref != gmsh_data.m_physical_ref_map.end()) {
+        return i_physical_ref->second.refId();
+      } else {
+        return RefId{ref_point_list.first};
+      }
+    }();
 
     bool is_boundary = true;
     for (size_t i_node = 0; i_node < point_list.size(); ++i_node) {
@@ -611,7 +685,7 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
       }
     }
 
-    descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list, is_boundary));
+    descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
   }
 
   descriptor.cell_owner_vector.resize(nb_cells);
@@ -740,6 +814,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
     this->__proceedData();
   }
   std::cout << std::flush;
+
   if (parallel::size() > 1) {
     std::cout << "Sequential mesh read! Need to be dispatched\n" << std::flush;
 
@@ -763,6 +838,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
 
       return *begin(dimension_set);
     }();
+
     switch (mesh_dimension) {
     case 1: {
       this->_dispatch<1>();
@@ -1440,6 +1516,7 @@ GmshReader::__proceedData()
       xr[i][2] = m_mesh_data.__vertices[i][2];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<3>();
 
   } else if (dot(dimension2_mask, elementNumber) > 0) {
     const size_t nb_cells = dot(dimension2_mask, elementNumber);
@@ -1459,6 +1536,7 @@ GmshReader::__proceedData()
       xr[i][1] = m_mesh_data.__vertices[i][1];
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<2>();
 
   } else if (dot(dimension1_mask, elementNumber) > 0) {
     const size_t nb_cells = dot(dimension1_mask, elementNumber);
@@ -1478,6 +1556,7 @@ GmshReader::__proceedData()
     }
 
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
+    this->_checkMesh<1>();
 
   } else {
     throw NormalError("using a dimension 0 mesh is forbidden");
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 9e62baf56e7915379d65a5c88c21fe7411fca253..a3a8440bd25499d8aea34a08e7c836b66a4c3a32 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -5,12 +5,14 @@
 #include <mesh/ConnectivityDispatcher.hpp>
 #include <mesh/ItemId.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 
 #include <vector>
 
-template <int Dimension>
+template <size_t Dimension>
 void
 MeshBuilderBase::_dispatch()
 {
@@ -28,8 +30,8 @@ MeshBuilderBase::_dispatch()
     NodeValue<Rd> xr;
     m_mesh = std::make_shared<MeshType>(connectivity, xr);
   }
-  const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);
 
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*m_mesh);
   ConnectivityDispatcher<Dimension> dispatcher(mesh.connectivity());
 
   std::shared_ptr dispatched_connectivity = dispatcher.dispatchedConnectivity();
@@ -41,3 +43,165 @@ MeshBuilderBase::_dispatch()
 template void MeshBuilderBase::_dispatch<1>();
 template void MeshBuilderBase::_dispatch<2>();
 template void MeshBuilderBase::_dispatch<3>();
+
+template <size_t Dimension>
+void
+MeshBuilderBase::_checkMesh() const
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+
+  if (not m_mesh) {
+    throw UnexpectedError("mesh is not built yet");
+  }
+
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*m_mesh);
+
+  const ConnectivityType& connectivity = mesh.connectivity();
+
+  if constexpr (Dimension > 2) {   // check for duplicated edges
+    auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+
+    std::vector<std::vector<EdgeId>> node_edges(mesh.numberOfNodes());
+    for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
+      auto node_list = edge_to_node_matrix[edge_id];
+      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+        node_edges[node_list[i_node]].push_back(edge_id);
+      }
+    }
+
+    for (auto&& edge_list : node_edges) {
+      std::sort(edge_list.begin(), edge_list.end());
+    }
+
+    for (EdgeId edge_id = 0; edge_id < mesh.numberOfEdges(); ++edge_id) {
+      auto node_list                   = edge_to_node_matrix[edge_id];
+      std::vector<EdgeId> intersection = node_edges[node_list[0]];
+      for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+        std::vector<EdgeId> local_intersection;
+        std::set_intersection(intersection.begin(), intersection.end(),   //
+                              node_edges[node_list[i_node]].begin(), node_edges[node_list[i_node]].end(),
+                              std::back_inserter(local_intersection));
+        std::swap(local_intersection, intersection);
+        if (intersection.size() < 2) {
+          break;
+        }
+      }
+
+      if (intersection.size() > 1) {
+        std::ostringstream error_msg;
+        error_msg << "invalid mesh.\n\tFollowing edges\n";
+        for (EdgeId edge_id : intersection) {
+          error_msg << "\t - id=" << edge_id << " number=" << connectivity.edgeNumber()[edge_id] << '\n';
+        }
+        error_msg << "\tare duplicated";
+        throw NormalError(error_msg.str());
+      }
+    }
+  }
+
+  if constexpr (Dimension > 1) {   // check for duplicated faces
+    auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+
+    std::vector<std::vector<FaceId>> node_faces(mesh.numberOfNodes());
+    for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
+      auto node_list = face_to_node_matrix[face_id];
+      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+        node_faces[node_list[i_node]].push_back(face_id);
+      }
+    }
+
+    for (auto&& face_list : node_faces) {
+      std::sort(face_list.begin(), face_list.end());
+    }
+
+    for (FaceId face_id = 0; face_id < mesh.numberOfFaces(); ++face_id) {
+      auto node_list                   = face_to_node_matrix[face_id];
+      std::vector<FaceId> intersection = node_faces[node_list[0]];
+      for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+        std::vector<FaceId> local_intersection;
+        std::set_intersection(intersection.begin(), intersection.end(),   //
+                              node_faces[node_list[i_node]].begin(), node_faces[node_list[i_node]].end(),
+                              std::back_inserter(local_intersection));
+        std::swap(local_intersection, intersection);
+        if (intersection.size() < 2) {
+          break;
+        }
+      }
+
+      if (intersection.size() > 1) {
+        std::ostringstream error_msg;
+        error_msg << "invalid mesh.\n\tFollowing faces\n";
+        for (FaceId face_id : intersection) {
+          error_msg << "\t - id=" << face_id << " number=" << connectivity.faceNumber()[face_id] << '\n';
+        }
+        error_msg << "\tare duplicated";
+        throw NormalError(error_msg.str());
+      }
+    }
+  }
+
+  auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+  {   // check for duplicated cells
+    std::vector<std::vector<CellId>> node_cells(mesh.numberOfNodes());
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto node_list = cell_to_node_matrix[cell_id];
+      for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+        node_cells[node_list[i_node]].push_back(cell_id);
+      }
+    }
+
+    for (auto&& cell_list : node_cells) {
+      std::sort(cell_list.begin(), cell_list.end());
+    }
+
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto node_list                   = cell_to_node_matrix[cell_id];
+      std::vector<CellId> intersection = node_cells[node_list[0]];
+      for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+        std::vector<CellId> local_intersection;
+        std::set_intersection(intersection.begin(), intersection.end(),   //
+                              node_cells[node_list[i_node]].begin(), node_cells[node_list[i_node]].end(),
+                              std::back_inserter(local_intersection));
+        std::swap(local_intersection, intersection);
+        if (intersection.size() < 2) {
+          break;
+        }
+      }
+
+      if (intersection.size() > 1) {
+        std::ostringstream error_msg;
+        error_msg << "invalid mesh.\n\tFollowing cells\n";
+        for (CellId cell_id : intersection) {
+          error_msg << "\t - id=" << cell_id << " number=" << connectivity.cellNumber()[cell_id] << '\n';
+        }
+        error_msg << "\tare duplicated";
+        throw NormalError(error_msg.str());
+      }
+    }
+  }
+
+  const auto& Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+  const auto& xr  = mesh.xr();
+
+  for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+    double cell_volume = 0;
+    auto cell_nodes    = cell_to_node_matrix[cell_id];
+    for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+      cell_volume += dot(Cjr(cell_id, i_node), xr[cell_nodes[i_node]]);
+    }
+
+    if (cell_volume <= 0) {
+      std::ostringstream error_msg;
+      error_msg << "invalid mesh.\n\tThe following cell\n";
+      error_msg << "\t - id=" << cell_id << " number=" << connectivity.cellNumber()[cell_id] << '\n';
+      error_msg << "\thas non-positive volume: " << cell_volume / Dimension;
+      throw NormalError(error_msg.str());
+    }
+  }
+}
+
+template void MeshBuilderBase::_checkMesh<1>() const;
+template void MeshBuilderBase::_checkMesh<2>() const;
+template void MeshBuilderBase::_checkMesh<3>() const;
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 6c0328bb8ee624e19d42fb4660ba3665ff7164de..63c55ec565b1f48023906f5aa602db13ef76d62c 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -10,9 +10,12 @@ class MeshBuilderBase
  protected:
   std::shared_ptr<const IMesh> m_mesh;
 
-  template <int Dimension>
+  template <size_t Dimension>
   void _dispatch();
 
+  template <size_t Dimension>
+  void _checkMesh() const;
+
  public:
   std::shared_ptr<const IMesh>
   mesh() const
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 9b0c89c8a6df4c4aab44b4b2f1aadc0615265fe8..a448862a01391a7c51bd1e7e9d75aff9e3df996e 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -9,6 +9,8 @@
 #include <utils/Messenger.hpp>
 #include <utils/PugsUtils.hpp>
 
+#include <output/VTKWriter.hpp>
+
 #include <map>
 
 template <size_t Dimension>
@@ -535,7 +537,12 @@ class MeshData : public IMeshData
     }();
 
     if (not parallel::allReduceAnd(is_valid)) {
-      throw NormalError("mesh contains cells of non-positive volume");
+      VTKWriter writer("bad_mesh");
+      writer.writeMesh(m_mesh);
+      std::ostringstream error_msg;
+      error_msg << "mesh contains cells of non-positive volume (see " << rang::fgB::yellow << "bad_mesh.pvd"
+                << rang::fg::reset << " file).";
+      throw NormalError(error_msg.str());
     }
   }
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index ef976be76e21b8d28ecc307a58a77161ffa2804c..4a346cf33e4c7f475a0abb8f0c3136ed88a72cc0 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -147,20 +147,20 @@ GnuplotWriter::_writeData(const ItemDataT& item_data,
 
 template <typename MeshType>
 void
-GnuplotWriter::_writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
+GnuplotWriter::_writeDataAtNodes(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
                                  std::ostream& fout) const
 {
   if constexpr (MeshType::Dimension == 1) {
-    auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-    auto cell_is_owned       = mesh->connectivity().cellIsOwned();
+    auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    auto cell_is_owned       = mesh.connectivity().cellIsOwned();
 
-    for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
       if (cell_is_owned[cell_id]) {
         const auto& cell_nodes = cell_to_node_matrix[cell_id];
         for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
           const NodeId& node_id                     = cell_nodes[i_node];
-          const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
+          const TinyVector<MeshType::Dimension>& xr = mesh.xr()[node_id];
           fout << xr[0];
           for (auto [name, item_data_variant] : output_named_item_data_set) {
             std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
@@ -172,15 +172,15 @@ GnuplotWriter::_writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
       }
     }
   } else if constexpr (MeshType::Dimension == 2) {
-    auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-    auto cell_is_owned       = mesh->connectivity().cellIsOwned();
+    auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    auto cell_is_owned       = mesh.connectivity().cellIsOwned();
 
-    for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
       if (cell_is_owned[cell_id]) {
         const auto& cell_nodes = cell_to_node_matrix[cell_id];
         for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
           const NodeId& node_id                     = cell_nodes[i_node];
-          const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
+          const TinyVector<MeshType::Dimension>& xr = mesh.xr()[node_id];
           fout << xr[0] << ' ' << xr[1];
           for (auto [name, item_data_variant] : output_named_item_data_set) {
             std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
@@ -188,7 +188,7 @@ GnuplotWriter::_writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
           fout << '\n';
         }
         const NodeId& node_id                     = cell_nodes[0];
-        const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
+        const TinyVector<MeshType::Dimension>& xr = mesh.xr()[node_id];
         fout << xr[0] << ' ' << xr[1];
         for (auto [name, item_data_variant] : output_named_item_data_set) {
           std::visit([&](auto&& item_data) { _writeData(item_data, cell_id, node_id, fout); }, item_data_variant);
@@ -239,7 +239,7 @@ GnuplotWriter::_writeNodeData(const NodeArray<DataType>& node_array, NodeId node
 
 template <typename MeshType>
 void
-GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
+GnuplotWriter::_write(const MeshType& mesh,
                       const OutputNamedItemDataSet& output_named_item_data_set,
                       std::optional<double> time) const
 {
@@ -296,64 +296,64 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 }
 
 void
-GnuplotWriter::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+GnuplotWriter::_writeAtTime(const IMesh& mesh,
                             const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                             double time) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh.dimension()));
   }
   }
 }
 
 void
-GnuplotWriter::_writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
+GnuplotWriter::_writeMesh(const IMesh& mesh) const
 {
   OutputNamedItemDataSet output_named_item_data_set{};
 
-  switch (p_mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + stringify(p_mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh.dimension()));
   }
   }
 }
 
 void
-GnuplotWriter::_write(const std::shared_ptr<const IMesh>& mesh,
+GnuplotWriter::_write(const IMesh& mesh,
                       const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh.dimension()));
   }
   }
 }
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 7d4c0c8a430df6dec0eb21a8246b848dcc41cc39..6a4c6bd5b4765c6402f56fdc5b700cd2b584cac7 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -41,23 +41,23 @@ class GnuplotWriter final : public WriterBase
                   std::ostream& fout) const;
 
   template <typename MeshType>
-  void _writeDataAtNodes(const std::shared_ptr<const MeshType>& mesh,
+  void _writeDataAtNodes(const MeshType& mesh,
                          const OutputNamedItemDataSet& output_named_item_data_set,
                          std::ostream& fout) const;
 
   template <typename MeshType>
-  void _write(const std::shared_ptr<const MeshType>& mesh,
+  void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
 
-  void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+  void _writeAtTime(const IMesh& mesh,
                     const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                     double time) const final;
 
-  void _write(const std::shared_ptr<const IMesh>& mesh,
+  void _write(const IMesh& mesh,
               const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
 
-  void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+  void _writeMesh(const IMesh& mesh) const final;
 
  public:
   GnuplotWriter(const std::string& base_filename) : WriterBase(base_filename) {}
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index e98ebdd17780b385783eff189538c0ae90f6dcc6..ad398126c4f075bf8c092ebb0aaac6ed1e4c0fd9 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -141,7 +141,7 @@ GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array
 
 template <typename MeshType, ItemType item_type>
 void
-GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
+GnuplotWriter1D::_writeItemDatas(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
                                  std::ostream& fout) const
 {
@@ -155,12 +155,12 @@ GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
     return number_of_columns;
   }();
 
-  auto is_owned = mesh->connectivity().template isOwned<item_type>();
+  auto is_owned = mesh.connectivity().template isOwned<item_type>();
 
   const size_t& number_of_owned_lines = [&]() {
     if (parallel::size() > 1) {
       size_t number_of_owned_items = 0;
-      for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
+      for (ItemId item_id = 0; item_id < mesh.template numberOf<item_type>(); ++item_id) {
         if (is_owned[item_id]) {
           ++number_of_owned_items;
         }
@@ -168,27 +168,27 @@ GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
 
       return number_of_owned_items;
     } else {
-      return mesh->template numberOf<item_type>();
+      return mesh.template numberOf<item_type>();
     }
   }();
 
   Array<double> values{number_of_columns * number_of_owned_lines};
 
   if constexpr (item_type == ItemType::cell) {
-    auto& mesh_data         = MeshDataManager::instance().getMeshData(*mesh);
+    auto& mesh_data         = MeshDataManager::instance().getMeshData(mesh);
     const auto& cell_center = mesh_data.xj();
 
     size_t index = 0;
-    for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
+    for (ItemId item_id = 0; item_id < mesh.template numberOf<item_type>(); ++item_id) {
       if (is_owned[item_id]) {
         values[number_of_columns * index++] = cell_center[item_id][0];
       }
     }
   } else if constexpr (item_type == ItemType::node) {
-    const auto& node_position = mesh->xr();
+    const auto& node_position = mesh.xr();
 
     size_t index = 0;
-    for (ItemId item_id = 0; item_id < mesh->template numberOf<item_type>(); ++item_id) {
+    for (ItemId item_id = 0; item_id < mesh.template numberOf<item_type>(); ++item_id) {
       if (is_owned[item_id]) {
         values[number_of_columns * index++] = node_position[item_id][0];
       }
@@ -274,7 +274,7 @@ GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
 
 template <typename MeshType>
 void
-GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
+GnuplotWriter1D::_write(const MeshType& mesh,
                         const OutputNamedItemDataSet& output_named_item_data_set,
                         std::optional<double> time) const
 {
@@ -341,7 +341,7 @@ GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
 }
 
 void
-GnuplotWriter1D::_writeMesh(const std::shared_ptr<const IMesh>&) const
+GnuplotWriter1D::_writeMesh(const IMesh&) const
 {
   std::ostringstream errorMsg;
   errorMsg << "gnuplot_1d_writer does not write meshes\n"
@@ -351,47 +351,47 @@ GnuplotWriter1D::_writeMesh(const std::shared_ptr<const IMesh>&) const
 }
 
 void
-GnuplotWriter1D::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+GnuplotWriter1D::_writeAtTime(const IMesh& mesh,
                               const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                               double time) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
     std::ostringstream errorMsg;
-    errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(mesh->dimension()) << '\n'
+    errorMsg << "gnuplot_1d_writer is not available in dimension " << stringify(mesh.dimension()) << '\n'
              << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
              << "gnuplot_writer" << rang::fg::reset << " in dimension 2";
     throw NormalError(errorMsg.str());
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh.dimension()));
   }
   }
 }
 
 void
-GnuplotWriter1D::_write(const std::shared_ptr<const IMesh>& mesh,
+GnuplotWriter1D::_write(const IMesh& mesh,
                         const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   default: {
-    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh->dimension()));
+    throw NormalError("gnuplot format is not available in dimension " + stringify(mesh.dimension()));
   }
   }
 }
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 93f81404659e8ef9aa869ad59b9a78885381a70a..fd0cac99ba56322100536483d14e78335015265d 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -38,25 +38,25 @@ class GnuplotWriter1D final : public WriterBase
   size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
 
   template <typename MeshType, ItemType item_type>
-  void _writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
+  void _writeItemDatas(const MeshType& mesh,
                        const OutputNamedItemDataSet& output_named_item_data_set,
                        std::ostream& fout) const;
 
   void _writePreamble(const OutputNamedItemDataSet& output_named_item_value_set, std::ostream& fout) const;
 
   template <typename MeshType>
-  void _write(const std::shared_ptr<const MeshType>& mesh,
+  void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_value_set,
               std::optional<double> time) const;
 
-  void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+  void _writeAtTime(const IMesh& mesh,
                     const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                     double time) const final;
 
-  void _write(const std::shared_ptr<const IMesh>& mesh,
+  void _write(const IMesh& mesh,
               const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
 
-  void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+  void _writeMesh(const IMesh& mesh) const final;
 
  public:
   GnuplotWriter1D(const std::string& base_filename) : WriterBase(base_filename) {}
diff --git a/src/output/IWriter.hpp b/src/output/IWriter.hpp
index a277bdea3c7d33051ddf228b3ce95b052655b454..14c77330065fbf35b7edf707f110ae621e4b22bf 100644
--- a/src/output/IWriter.hpp
+++ b/src/output/IWriter.hpp
@@ -12,6 +12,7 @@ class IWriter
 {
  public:
   virtual void writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
+  virtual void writeMesh(const IMesh& mesh) const                        = 0;
 
   virtual void write(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const = 0;
 
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 342a384e7f0ea4989b53fd9fd39ede1d0973171c..f68607d7da6daf5afe8980a36e669e87b013dd3e 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -355,14 +355,14 @@ VTKWriter::_write_node_data(std::ofstream& os,
 
 template <typename MeshType>
 void
-VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
+VTKWriter::_write(const MeshType& mesh,
                   const OutputNamedItemDataSet& given_output_named_item_data_set,
                   std::optional<double> time) const
 {
   OutputNamedItemDataSet output_named_item_data_set{given_output_named_item_data_set};
   // Adding basic mesh information
-  output_named_item_data_set.add(NamedItemData{"cell_number", mesh->connectivity().cellNumber()});
-  output_named_item_data_set.add(NamedItemData{"node_number", mesh->connectivity().nodeNumber()});
+  output_named_item_data_set.add(NamedItemData{"cell_number", mesh.connectivity().cellNumber()});
+  output_named_item_data_set.add(NamedItemData{"node_number", mesh.connectivity().nodeNumber()});
 
   if (parallel::rank() == 0) {   // write PVTK file
     std::ofstream fout(_getFilenamePVTU());
@@ -445,7 +445,7 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 #endif
     fout << ">\n";
     fout << "<UnstructuredGrid>\n";
-    fout << "<Piece NumberOfPoints=\"" << mesh->numberOfNodes() << "\" NumberOfCells=\"" << mesh->numberOfCells()
+    fout << "<Piece NumberOfPoints=\"" << mesh.numberOfNodes() << "\" NumberOfCells=\"" << mesh.numberOfCells()
          << "\">\n";
     fout << "<CellData>\n";
     for (const auto& [name, item_value_variant] : output_named_item_data_set) {
@@ -464,10 +464,10 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << "<Points>\n";
     {
       using Rd                      = TinyVector<MeshType::Dimension>;
-      const NodeValue<const Rd>& xr = mesh->xr();
-      Array<TinyVector<3>> positions(mesh->numberOfNodes());
+      const NodeValue<const Rd>& xr = mesh.xr();
+      Array<TinyVector<3>> positions(mesh.numberOfNodes());
       parallel_for(
-        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
           for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
             positions[r][i] = xr[r][i];
           }
@@ -481,16 +481,16 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 
     fout << "<Cells>\n";
     {
-      const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+      const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
       _write_array(fout, "connectivity", cell_to_node_matrix.values(), serialize_data_list);
     }
 
     {
-      const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-      Array<unsigned int> offsets(mesh->numberOfCells());
+      const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+      Array<unsigned int> offsets(mesh.numberOfCells());
       unsigned int offset = 0;
-      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+      for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
         offset += cell_nodes.size();
         offsets[j] = offset;
@@ -499,12 +499,12 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     }
 
     {
-      Array<int8_t> types(mesh->numberOfCells());
-      const auto& cell_type           = mesh->connectivity().cellType();
-      const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+      Array<int8_t> types(mesh.numberOfCells());
+      const auto& cell_type           = mesh.connectivity().cellType();
+      const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
       parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+        mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
           switch (cell_type[j]) {
           case CellType::Line: {
             types[j] = 3;
@@ -563,14 +563,14 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
           return false;
         }();
         if (has_general_polyhedron) {
-          const auto& cell_to_face_matrix   = mesh->connectivity().cellToFaceMatrix();
-          const auto& face_to_node_matrix   = mesh->connectivity().faceToNodeMatrix();
-          const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+          const auto& cell_to_face_matrix   = mesh.connectivity().cellToFaceMatrix();
+          const auto& face_to_node_matrix   = mesh.connectivity().faceToNodeMatrix();
+          const auto& cell_face_is_reversed = mesh.connectivity().cellFaceIsReversed();
 
-          Array<size_t> faces_offsets(mesh->numberOfCells());
+          Array<size_t> faces_offsets(mesh.numberOfCells());
           size_t next_offset = 0;
           fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
-          for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
             const auto& cell_nodes = cell_to_node_matrix[cell_id];
             std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
             for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
@@ -658,21 +658,21 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
 }
 
 void
-VTKWriter::_writeMesh(const std::shared_ptr<const IMesh>& mesh) const
+VTKWriter::_writeMesh(const IMesh& mesh) const
 {
   OutputNamedItemDataSet output_named_item_data_set;
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<3>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   default: {
@@ -682,23 +682,23 @@ VTKWriter::_writeMesh(const std::shared_ptr<const IMesh>& mesh) const
 }
 
 void
-VTKWriter::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+VTKWriter::_writeAtTime(const IMesh& mesh,
                         const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                         double time) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_data_set, time);
+    this->_write(dynamic_cast<const Mesh<Connectivity<3>>&>(mesh), output_named_item_data_set, time);
     break;
   }
   default: {
@@ -708,22 +708,22 @@ VTKWriter::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
 }
 
 void
-VTKWriter::_write(const std::shared_ptr<const IMesh>& mesh,
+VTKWriter::_write(const IMesh& mesh,
                   const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
 {
   OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
 
-  switch (mesh->dimension()) {
+  switch (mesh.dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<1>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<2>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_data_set, {});
+    this->_write(dynamic_cast<const Mesh<Connectivity<3>>&>(mesh), output_named_item_data_set, {});
     break;
   }
   default: {
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 1b6fc35f6ee02dc7a441f9dc956165af32c33276..1fa439e94c8f1b48930df2835d68caa4030d7061 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -96,18 +96,18 @@ class VTKWriter final : public WriterBase
                         SerializedDataList& serialized_data_list) const;
 
   template <typename MeshType>
-  void _write(const std::shared_ptr<const MeshType>& mesh,
+  void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
 
-  void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+  void _writeAtTime(const IMesh& mesh,
                     const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                     double time) const final;
 
-  void _write(const std::shared_ptr<const IMesh>& mesh,
+  void _write(const IMesh& mesh,
               const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
 
-  void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+  void _writeMesh(const IMesh& mesh) const final;
 
  public:
   VTKWriter(const std::string& base_filename) : WriterBase(base_filename) {}
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 0ce785cd9e6bd3f3985466d6f2f90e38fc49b802..fcb2cf78f21f3bebc878b0343826c3abab4636c1 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -370,7 +370,7 @@ WriterBase::write(const std::vector<std::shared_ptr<const INamedDiscreteData>>&
     throw NormalError("this writer requires time value");
   } else {
     std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
-    this->_write(mesh, named_discrete_data_list);
+    this->_write(*mesh, named_discrete_data_list);
   }
 }
 
@@ -385,7 +385,7 @@ WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const INamedDiscrete
 
     if (time >= m_period_manager->nextTime()) {
       std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
-      this->_writeAtTime(mesh, named_discrete_data_list, time);
+      this->_writeAtTime(*mesh, named_discrete_data_list, time);
       m_period_manager->setSaveTime(time);
     }
   } else {
@@ -401,7 +401,7 @@ WriterBase::writeForced(const std::vector<std::shared_ptr<const INamedDiscreteDa
     if (time == m_period_manager->getLastTime())
       return;   // output already performed
     std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
-    this->_writeAtTime(mesh, named_discrete_data_list, time);
+    this->_writeAtTime(*mesh, named_discrete_data_list, time);
     m_period_manager->setSaveTime(time);
   } else {
     throw NormalError("this writer does not allow time value");
@@ -417,7 +417,7 @@ WriterBase::writeOnMesh(const std::shared_ptr<const IMesh>& mesh,
   } else {
     _checkMesh(mesh, named_discrete_data_list);
     _checkConnectivity(mesh, named_discrete_data_list);
-    this->_write(mesh, named_discrete_data_list);
+    this->_write(*mesh, named_discrete_data_list);
   }
 }
 
@@ -431,7 +431,7 @@ WriterBase::writeOnMeshIfNeeded(const std::shared_ptr<const IMesh>& mesh,
       return;   // output already performed
     _checkMesh(mesh, named_discrete_data_list);
     _checkConnectivity(mesh, named_discrete_data_list);
-    this->_writeAtTime(mesh, named_discrete_data_list, time);
+    this->_writeAtTime(*mesh, named_discrete_data_list, time);
     m_period_manager->setSaveTime(time);
   } else {
     throw NormalError("this writer does not allow time value");
@@ -448,7 +448,7 @@ WriterBase::writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
       return;   // output already performed
     _checkMesh(mesh, named_discrete_data_list);
     _checkConnectivity(mesh, named_discrete_data_list);
-    this->_writeAtTime(mesh, named_discrete_data_list, time);
+    this->_writeAtTime(*mesh, named_discrete_data_list, time);
     m_period_manager->setSaveTime(time);
   } else {
     throw NormalError("this writer does not allow time value");
@@ -457,6 +457,12 @@ WriterBase::writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
 
 void
 WriterBase::writeMesh(const std::shared_ptr<const IMesh>& mesh) const
+{
+  writeMesh(*mesh);
+}
+
+void
+WriterBase::writeMesh(const IMesh& mesh) const
 {
   if (m_period_manager.has_value()) {
     throw NormalError("write_mesh requires a writer without time period");
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
index 6f6a8e143b5b181bed10d76223f809851cda4229..339cd7645c8cbdc903357cee6000d1c6ab45dda3 100644
--- a/src/output/WriterBase.hpp
+++ b/src/output/WriterBase.hpp
@@ -109,14 +109,14 @@ class WriterBase : public IWriter
   OutputNamedItemDataSet _getOutputNamedItemDataSet(
     const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const;
 
-  virtual void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+  virtual void _writeAtTime(const IMesh& mesh,
                             const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
                             double time) const = 0;
 
-  virtual void _write(const std::shared_ptr<const IMesh>& mesh,
+  virtual void _write(const IMesh& mesh,
                       const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const = 0;
 
-  virtual void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
+  virtual void _writeMesh(const IMesh& mesh) const = 0;
 
  public:
   void write(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
@@ -138,6 +138,7 @@ class WriterBase : public IWriter
                          double time) const final;
 
   void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+  void writeMesh(const IMesh& mesh) const final;
 
   WriterBase() = delete;