diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index bafef08561b7900ecda859ecaae4b54b61119c57..b25ca82accf3189fbdab58ca4161101647cb463d 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -67,16 +67,12 @@ GnuplotWriter1D::_is_node_value(const NodeValue<const DataType>&) const
   return true;
 }
 
-template <size_t Dimension>
 void
 GnuplotWriter1D::_writePreamble(const OutputNamedItemValueSet& output_named_item_value_set, std::ostream& fout) const
 {
   fout << "# list of data\n";
   fout << "# 1:x";
-  if constexpr (Dimension > 1) {
-    fout << " 2:y";
-  }
-  uint64_t i = Dimension + 1;
+  uint64_t i = 2;
   for (const auto& i_named_item_value : output_named_item_value_set) {
     const std::string name         = i_named_item_value.first;
     const auto& item_value_variant = i_named_item_value.second;
@@ -105,139 +101,135 @@ GnuplotWriter1D::_writePreamble(const OutputNamedItemValueSet& output_named_item
   fout << "\n\n";
 }
 
-template <typename DataType>
-void
-GnuplotWriter1D::_writeCellValue(const NodeValue<DataType>&, CellId, std::ostream&) const
-{}
-
-template <typename DataType>
-void
-GnuplotWriter1D::_writeCellValue(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriter1D::_itemValueNbRow(const ItemValue<DataType, item_type>&) const
 {
-  const auto& value = cell_value[cell_id];
   if constexpr (std::is_arithmetic_v<DataType>) {
-    fout << ' ' << value;
+    return 1;
   } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
-    for (size_t i = 0; i < value.dimension(); ++i) {
-      fout << ' ' << value[i];
-    }
+    return DataType::Dimension;
   } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
-    for (size_t i = 0; i < value.nbRows(); ++i) {
-      for (size_t j = 0; j < value.nbColumns(); ++j) {
-        fout << ' ' << value(i, j);
-      }
-    }
+    return DataType{}.dimension();
   } else {
     throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
   }
 }
 
-template <typename MeshType>
+template <typename MeshType, ItemType item_type>
 void
-GnuplotWriter1D::_writeCellValues(const std::shared_ptr<const MeshType>& mesh,
+GnuplotWriter1D::_writeItemValues(const std::shared_ptr<const MeshType>& mesh,
                                   const OutputNamedItemValueSet& output_named_item_value_set,
                                   std::ostream& fout) const
 {
-  if constexpr (MeshType::Dimension == 1) {
-    auto& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+  using ItemId = ItemIdT<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';
+  const size_t& number_of_columns = [&] {
+    size_t number_of_columns = 1;
+    for (auto [name, item_value] : output_named_item_value_set) {
+      std::visit([&](auto&& value) { number_of_columns += _itemValueNbRow(value); }, item_value);
     }
-  } else if constexpr (MeshType::Dimension == 2) {
-    auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+    return number_of_columns;
+  }();
 
-    for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-      std::ostringstream os;
-      os.precision(15);
-      os.setf(std::ios_base::scientific);
+  auto is_owned = mesh->connectivity().template isOwned<item_type>();
 
-      for (auto [name, item_value] : output_named_item_value_set) {
-        std::visit([&](auto&& cell_value) { _writeCellValue(cell_value, cell_id, os); }, item_value);
+  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) {
+        if (is_owned[item_id]) {
+          ++number_of_owned_items;
+        }
       }
 
-      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';
+      return number_of_owned_items;
+    } else {
+      return mesh->template numberOf<item_type>();
     }
-  } else {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-}
+  }();
 
-template <typename ValueType>
-void
-GnuplotWriter1D::_writeNodeValue(const NodeValue<ValueType>& node_value, NodeId node_id, std::ostream& fout) const
-{
-  const auto& value = node_value[node_id];
-  if constexpr (std::is_arithmetic_v<ValueType>) {
-    fout << ' ' << value;
-  } else if constexpr (is_tiny_vector_v<std::decay_t<ValueType>>) {
-    for (size_t i = 0; i < value.dimension(); ++i) {
-      fout << ' ' << value[i];
+  Array<double> values{number_of_columns * number_of_owned_lines};
+
+  if constexpr (item_type == ItemType::cell) {
+    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) {
+      if (is_owned[item_id]) {
+        values[number_of_columns * index++] = cell_center[item_id][0];
+      }
     }
-  } else if constexpr (is_tiny_matrix_v<std::decay_t<ValueType>>) {
-    for (size_t i = 0; i < value.nbRows(); ++i) {
-      for (size_t j = 0; j < value.nbColumns(); ++j) {
-        fout << ' ' << value(i, j);
+  } else if constexpr (item_type == ItemType::node) {
+    const auto& node_position = mesh->xr();
+
+    size_t index = 0;
+    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];
       }
     }
   } else {
-    throw UnexpectedError("invalid data type for cell value output: " + demangle<ValueType>());
+    throw UnexpectedError("invalid item type");
   }
-}
 
-template <typename ValueType>
-void
-GnuplotWriter1D::_writeNodeValue(const CellValue<ValueType>&, NodeId, std::ostream&) const
-{}
+  size_t column_number = 1;
+  for (auto [name, output_item_value] : output_named_item_value_set) {
+    std::visit(
+      [&](auto&& item_value) {
+        using ItemValueT = std::decay_t<decltype(item_value)>;
+        if constexpr (ItemValueT::item_t == item_type) {
+          using DataT  = std::decay_t<typename ItemValueT::data_type>;
+          size_t index = 0;
+          for (ItemId item_id = 0; item_id < item_value.size(); ++item_id) {
+            if (is_owned[item_id]) {
+              if constexpr (std::is_arithmetic_v<DataT>) {
+                values[number_of_columns * index + column_number] = item_value[item_id];
+              } else if constexpr (is_tiny_vector_v<DataT>) {
+                const size_t k = number_of_columns * index + column_number;
+                for (size_t j = 0; j < DataT::Dimension; ++j) {
+                  values[k + j] = item_value[item_id][j];
+                }
+              } else if constexpr (is_tiny_matrix_v<DataT>) {
+                size_t k = number_of_columns * index + column_number;
+                for (size_t i = 0; i < DataT{}.nbRows(); ++i) {
+                  for (size_t j = 0; j < DataT{}.nbColumns(); ++j) {
+                    values[k++] = item_value[item_id](i, j);
+                  }
+                }
+              }
+              ++index;
+            }
+          }
+        }
+        column_number += _itemValueNbRow(item_value);
+      },
+      output_item_value);
+  }
 
-template <typename MeshType>
-void
-GnuplotWriter1D::_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';
+  if (parallel::size() > 1) {
+    values = parallel::gatherVariable(values, 0);
+  }
+
+  if (parallel::rank() == 0) {
+    Assert(values.size() % number_of_columns == 0);
+
+    std::vector<size_t> line_numbers(values.size() / number_of_columns);
+    for (size_t i = 0; i < line_numbers.size(); ++i) {
+      line_numbers[i] = i;
     }
-  } 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';
+
+    std::sort(line_numbers.begin(), line_numbers.end(),
+              [&](size_t i, size_t j) { return values[i * number_of_columns] < values[j * number_of_columns]; });
+
+    for (auto i_line : line_numbers) {
+      fout << values[i_line * number_of_columns];
+      for (size_t j = 1; j < number_of_columns; ++j) {
+        fout << ' ' << values[i_line * number_of_columns + j];
       }
       fout << '\n';
     }
-  } else {
-    throw UnexpectedError("invalid mesh dimension");
   }
 }
 
@@ -263,19 +255,23 @@ GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
     throw NormalError("cannot store both node and cell values in a gnuplot file");
   }
 
-  std::ofstream fout(_getFilename());
-  fout.precision(15);
-  fout.setf(std::ios_base::scientific);
-  fout << _getDateAndVersionComment();
+  std::ofstream fout;
+
+  if (parallel::rank() == 0) {
+    fout.open(_getFilename());
+    fout.precision(15);
+    fout.setf(std::ios_base::scientific);
+    fout << _getDateAndVersionComment();
 
-  fout << "\n# time = " << time << "\n\n";
+    fout << "# time = " << time << "\n\n";
 
-  _writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
+    _writePreamble(output_named_item_value_set, fout);
+  }
 
   if (has_cell_value) {
-    this->_writeCellValues(mesh, output_named_item_value_set, fout);
+    this->_writeItemValues<MeshType, ItemType::cell>(mesh, output_named_item_value_set, fout);
   } else {   // has_node_value
-    this->_writeNodeValues(mesh, output_named_item_value_set, fout);
+    this->_writeItemValues<MeshType, ItemType::node>(mesh, output_named_item_value_set, fout);
   }
 }
 
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index a299b14eee175ea0f2928e32e0008b6421a61607..6402a6002a105e1ef8b85b60783dce7a1cba82bd 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -30,29 +30,14 @@ class GnuplotWriter1D : public WriterBase
   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, ItemType item_type>
+  size_t _itemValueNbRow(const ItemValue<DataType, item_type>&) 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,
+  template <typename MeshType, ItemType item_type>
+  void _writeItemValues(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 MeshType>