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;