diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index a1fcfd4dc0165b4eec1a15ae366903bee82f4879..436ae7ffa22d7f3aaa20cbd44e2331393958823c 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -39,34 +39,6 @@ GnuplotWriter::_getFilename() const
   return sout.str();
 }
 
-template <typename DataType>
-bool
-GnuplotWriter::_is_cell_value(const CellValue<const DataType>&) const
-{
-  return true;
-}
-
-template <typename DataType>
-bool
-GnuplotWriter::_is_cell_value(const NodeValue<const DataType>&) const
-{
-  return false;
-}
-
-template <typename DataType>
-bool
-GnuplotWriter::_is_node_value(const CellValue<const DataType>&) const
-{
-  return false;
-}
-
-template <typename DataType>
-bool
-GnuplotWriter::_is_node_value(const NodeValue<const DataType>&) const
-{
-  return true;
-}
-
 template <size_t Dimension>
 void
 GnuplotWriter::_writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const
@@ -105,11 +77,6 @@ GnuplotWriter::_writePreamble(const OutputNamedItemValueSet& output_named_item_v
   fout << "\n\n";
 }
 
-template <typename DataType>
-void
-GnuplotWriter::_writeCellValue(const NodeValue<DataType>&, CellId, std::ostream&) const
-{}
-
 template <typename DataType>
 void
 GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const
@@ -132,42 +99,50 @@ GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cel
   }
 }
 
-template <typename MeshType>
+template <typename DataType, ItemType item_type>
 void
-GnuplotWriter::_writeCellValues(const std::shared_ptr<const MeshType>& mesh,
-                                const OutputNamedItemValueSet& output_named_item_value_set,
-                                std::ostream& fout) const
+GnuplotWriter::_writeValue(const ItemValue<DataType, item_type>& item_value,
+                           CellId cell_id,
+                           NodeId node_id,
+                           std::ostream& fout) const
 {
-  if constexpr (MeshType::Dimension == 1) {
-    auto& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+  if constexpr (item_type == ItemType::cell) {
+    this->_writeCellValue(item_value, cell_id, fout);
+  } else if constexpr (item_type == ItemType::node) {
+    this->_writeNodeValue(item_value, node_id, fout);
+  } else {
+    throw UnexpectedError{"invalid item type"};
+  }
+}
 
-    const auto& cell_center = mesh_data.xj();
-    for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-      fout << cell_center[cell_id][0];
-      for (auto [name, item_value] : output_named_item_value_set) {
-        std::visit([&](auto&& cell_value) { _writeCellValue(cell_value, cell_id, fout); }, item_value);
-      }
-      fout << '\n';
-    }
-  } else if constexpr (MeshType::Dimension == 2) {
+template <typename MeshType>
+void
+GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
+                            const OutputNamedItemValueSet& output_named_item_value_set,
+                            std::ostream& fout) const
+{
+  if constexpr (MeshType::Dimension <= 2) {
     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) {
-      std::ostringstream os;
-      os.precision(15);
-      os.setf(std::ios_base::scientific);
-
-      for (auto [name, item_value] : output_named_item_value_set) {
-        std::visit([&](auto&& cell_value) { _writeCellValue(cell_value, cell_id, os); }, item_value);
-      }
+      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];
+          fout << xr[0];
+          if (MeshType::Dimension == 2) {
+            fout << ' ' << xr[1];
+          }
+          for (auto [name, item_value] : output_named_item_value_set) {
+            std::visit([&](auto&& cell_value) { _writeValue(cell_value, cell_id, node_id, fout); }, item_value);
+          }
 
-      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<2>& xr = mesh->xr()[node_id];
-        fout << xr[0] << ' ' << xr[1] << ' ' << os.str() << '\n';
+          fout << '\n';
+        }
+        fout << '\n';
       }
-      fout << '\n';
     }
   } else {
     throw UnexpectedError("invalid mesh dimension");
@@ -196,86 +171,32 @@ GnuplotWriter::_writeNodeValue(const NodeValue<ValueType>& node_value, NodeId no
   }
 }
 
-template <typename ValueType>
-void
-GnuplotWriter::_writeNodeValue(const CellValue<ValueType>&, NodeId, std::ostream&) const
-{}
-
-template <typename MeshType>
-void
-GnuplotWriter::_writeNodeValues(const std::shared_ptr<const MeshType>& mesh,
-                                const OutputNamedItemValueSet& output_named_item_value_set,
-                                std::ostream& fout) const
-{
-  if constexpr (MeshType::Dimension == 1) {
-    const auto& xr = mesh->xr();
-    for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-      fout << xr[node_id][0];
-      for (auto [name, item_value] : output_named_item_value_set) {
-        std::visit([&](auto&& node_value) { _writeNodeValue(node_value, node_id, fout); }, item_value);
-      }
-      fout << '\n';
-    }
-  } else if constexpr (MeshType::Dimension == 2) {
-    auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
-
-    for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-      std::ostringstream os;
-      os.precision(15);
-      os.setf(std::ios_base::scientific);
-
-      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<2>& xr = mesh->xr()[node_id];
-        fout << xr[0] << ' ' << xr[1] << ' ' << os.str();
-        for (auto [name, item_value] : output_named_item_value_set) {
-          std::visit([&](auto&& node_value) { _writeNodeValue(node_value, node_id, os); }, item_value);
-        }
-        fout << '\n';
-      }
-      fout << '\n';
-    }
-  } else {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-}
-
 template <typename MeshType>
 void
 GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
                       const OutputNamedItemValueSet& output_named_item_value_set,
                       double time) const
 {
-  bool has_cell_value = false;
-  for (const auto& [name, item_value_variant] : output_named_item_value_set) {
-    has_cell_value |=
-      std::visit([&](auto&& item_value) { return this->_is_cell_value(item_value); }, item_value_variant);
-  }
+  if (parallel::rank() == 0) {
+    std::ofstream fout{_getFilename()};
+    fout.precision(15);
+    fout.setf(std::ios_base::scientific);
 
-  bool has_node_value = false;
-  for (const auto& [name, item_value_variant] : output_named_item_value_set) {
-    has_node_value |=
-      std::visit([&](auto&& item_value) { return this->_is_node_value(item_value); }, item_value_variant);
-  }
+    fout << this->_getDateAndVersionComment();
+    fout << "\n# time = " << time << "\n\n";
 
-  if (has_cell_value and has_node_value) {
-    throw NormalError("cannot store both node and cell values in a gnuplot file");
+    this->_writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
   }
 
-  std::ofstream fout(_getFilename());
-  fout.precision(15);
-  fout.setf(std::ios_base::scientific);
-  fout << _getDateAndVersionComment();
-
-  fout << "\n# time = " << time << "\n\n";
+  for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+    if (i_rank == parallel::rank()) {
+      std::ofstream fout{_getFilename(), std::ios_base::app};
+      fout.precision(15);
+      fout.setf(std::ios_base::scientific);
 
-  _writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
-
-  if (has_cell_value) {
-    this->_writeCellValues(mesh, output_named_item_value_set, fout);
-  } else {   // has_node_value
-    this->_writeNodeValues(mesh, output_named_item_value_set, fout);
+      this->_writeValues(mesh, output_named_item_value_set, fout);
+    }
+    parallel::barrier();
   }
 }
 
@@ -291,15 +212,15 @@ GnuplotWriter::writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
 
   switch (p_mesh->dimension()) {
   case 1: {
-    _writePreamble<1>(output_named_item_value_set, fout);
-    this->_writeNodeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_value_set,
-                           fout);
+    this->_writePreamble<1>(output_named_item_value_set, fout);
+    this->_writeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_value_set,
+                       fout);
     break;
   }
   case 2: {
-    _writePreamble<2>(output_named_item_value_set, fout);
-    this->_writeNodeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_value_set,
-                           fout);
+    this->_writePreamble<2>(output_named_item_value_set, fout);
+    this->_writeValues(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_value_set,
+                       fout);
     break;
   }
   default: {
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 0006d5ccfde75c50f0ca8d73c0479fa2d6d6f00e..bf6487c998c29d8a5c12177a574eea587a1fb6ea 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -18,43 +18,26 @@ class GnuplotWriter : public WriterBase
 
   std::string _getFilename() const;
 
-  template <typename DataType>
-  bool _is_cell_value(const CellValue<const DataType>&) const;
-
-  template <typename DataType>
-  bool _is_cell_value(const NodeValue<const DataType>&) const;
-
-  template <typename DataType>
-  bool _is_node_value(const CellValue<const DataType>&) const;
-
-  template <typename DataType>
-  bool _is_node_value(const NodeValue<const DataType>&) const;
-
   template <typename DataType>
   void _writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
 
-  template <typename DataType>
-  void _writeCellValue(const NodeValue<DataType>& node_value, CellId cell_id, std::ostream& fout) const;
-
-  template <typename MeshType>
-  void _writeCellValues(const std::shared_ptr<const MeshType>& mesh,
-                        const OutputNamedItemValueSet& output_named_item_value_set,
-                        std::ostream& fout) const;
-
-  template <typename DataType>
-  void _writeNodeValue(const CellValue<DataType>& cell_value, NodeId cell_id, std::ostream& fout) const;
-
   template <typename DataType>
   void _writeNodeValue(const NodeValue<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
 
-  template <typename MeshType>
-  void _writeNodeValues(const std::shared_ptr<const MeshType>& mesh,
-                        const OutputNamedItemValueSet& output_named_item_value_set,
-                        std::ostream& fout) const;
-
   template <size_t Dimension>
   void _writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const;
 
+  template <typename DataType, ItemType item_type>
+  void _writeValue(const ItemValue<DataType, item_type>& item_value,
+                   CellId cell_id,
+                   NodeId node_id,
+                   std::ostream& fout) const;
+
+  template <typename MeshType>
+  void _writeValues(const std::shared_ptr<const MeshType>& mesh,
+                    const OutputNamedItemValueSet& output_named_item_value_set,
+                    std::ostream& fout) const;
+
   template <typename MeshType>
   void _write(const std::shared_ptr<const MeshType>& mesh,
               const OutputNamedItemValueSet& output_named_item_value_set,