diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 4b1bf84b80752bb8e3e2110fd3526e523658c60d..78b1e25ce8277888819395cafe92118a46d992fb 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -121,7 +121,7 @@ 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) {
+  if constexpr (MeshType::Dimension == 1) {
     auto cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
     auto cell_is_owned       = mesh->connectivity().cellIsOwned();
 
@@ -132,16 +132,38 @@ GnuplotWriter::_writeValues(const std::shared_ptr<const MeshType>& mesh,
           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);
           }
+          fout << '\n';
+        }
+
+        fout << "\n\n";
+      }
+    }
+  } else 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) {
+      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] << ' ' << 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);
+          }
           fout << '\n';
         }
-        fout << '\n';
+        const NodeId& node_id                     = cell_nodes[0];
+        const TinyVector<MeshType::Dimension>& xr = mesh->xr()[node_id];
+        fout << xr[0] << ' ' << 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);
+        }
+        fout << "\n\n\n";
       }
     }
   } else {
@@ -177,22 +199,26 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
                       const OutputNamedItemValueSet& output_named_item_value_set,
                       double time) const
 {
-  if (parallel::size() == 1) {
+  if (parallel::rank() == 0) {
     std::ofstream fout{_getFilename()};
     fout.precision(15);
     fout.setf(std::ios_base::scientific);
 
     fout << this->_getDateAndVersionComment();
-    fout << "\n# time = " << time << "\n\n";
+    fout << "# time = " << time << "\n\n";
 
     this->_writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
+  }
 
-    fout.precision(15);
-    fout.setf(std::ios_base::scientific);
+  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);
 
-    this->_writeValues(mesh, output_named_item_value_set, fout);
-  } else {
-    throw NotImplementedError("not implemented in parallel");
+      this->_writeValues(mesh, output_named_item_value_set, fout);
+    }
+    parallel::barrier();
   }
 }
 
@@ -201,22 +227,13 @@ GnuplotWriter::writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
 {
   OutputNamedItemValueSet output_named_item_value_set{};
 
-  std::ofstream fout(_getFilename());
-  fout.precision(15);
-  fout.setf(std::ios_base::scientific);
-  fout << _getDateAndVersionComment();
-
   switch (p_mesh->dimension()) {
   case 1: {
-    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);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh), output_named_item_value_set, 0);
     break;
   }
   case 2: {
-    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);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh), output_named_item_value_set, 0);
     break;
   }
   default: {