diff --git a/doc/userdoc.org b/doc/userdoc.org
index eaac7da46edb913aff2df4d527d7c9ac20456cfe..8270c6076e9ef1f65e217c06da6722f0d8a0c02c 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -4740,11 +4740,14 @@ Running this code produces the gnuplot file displayed in Figure
 
 ***** ~gnuplot~ writers <<gnuplot-writers>>
 
-There are two ~gnuplot~ writers. One is dedicated to output of dimension
-1 results (~gnuplot_1d_writer~) and the other allows post processing in
-dimension 1 and 2 (~gnuplot_writer~).
-
-Both of these writers can be used for single output or time series
+There are three ~gnuplot~ writers. One is dedicated to output of
+dimension 1 results (~gnuplot_1d_writer~), the second allows post
+processing in dimension 1 and 2 (~gnuplot_writer~) and the third can be
+used to write arbitrary data in the ~gnuplot~ format, whichever the
+space dimension is (~gnuplot_raw_writer~). The last one is useful to
+create scatter plots for instance.
+
+All of these writers can be used for single output or time series
 outputs. In the case of single output, the filename is completed by
 adding the extension ~.gnu~, in the case of time series, the filename is
 extending by adding ~.abcd.gnu~, where ~abcd~ is the number in the output
@@ -5006,6 +5009,58 @@ The gnuplot result is displayed on Figure [[fig:writer-gp-2d-cos-sin]].
 This is the time series function in the case of the ~gnuplot_writer~. It
 behaves similarly to [[gnuplot-1d-series]].
 
+****** ~gnuplot_raw_writer~ functions <<gnuplot-raw-writer-function>>
+
+These writers are used to store raw data to the ~gnuplot~ format. Doing
+so, it is the user responsibility to store coordinates in the file if
+needed. Therefore it can be used in any space dimension.
+
+******* ~gnuplot_raw_writer: string -> writer~
+
+Here is an illustrating example in dimension 3, see the result on
+Figure [[fig:writer-raw-sin]].
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0,0,0], [2,2,2], (10,10,10));
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^3 -> R, x -> sin(pi*dot(x,x));
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  let r:Vh, r = sqrt(dot(cell_center(m),cell_center(m)));
+
+  write(gnuplot_raw_writer("raw_sin"), (name_output(r,"r"), name_output(fh, "f")));
+#+END_SRC
+
+#+NAME: writer-raw-sin-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/raw-sin.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 628,400
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/raw_sin.gnu)' lw 2 w lp
+#+END_SRC
+
+#+CAPTION: Example of produced gnuplot results from the ~gnuplot_raw_writer~.
+#+NAME: fig:writer-raw-sin
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-raw-sin-img
+
+******* ~gnuplot_raw_writer: string*R -> writer~
+
+This is the time series function in the case of the ~gnuplot_raw_writer~. It
+behaves similarly to [[gnuplot-1d-series]].
+
 ***** ~vtk~ writers
 
 For more complex post processing (including 3d), ~pugs~ can generate ~vtk~
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 7b404fcf670c5217721efba890728ce97e3e91a4..5ad47feff49b54cb5f90d2d3d9ddd68641204fbe 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -651,6 +651,27 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("cell_center",
+                            std::function(
+
+                              [](const std::shared_ptr<const MeshVariant>& mesh_v)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                return std::visit(
+                                  [&](auto&& mesh) {
+                                    using MeshType = mesh_type_t<decltype(mesh)>;
+                                    if constexpr (is_polygonal_mesh_v<MeshType>) {
+                                      return std::make_shared<DiscreteFunctionVariant>(
+                                        DiscreteFunctionP0(mesh_v,
+                                                           MeshDataManager::instance().getMeshData(*mesh).xj()));
+                                    } else {
+                                      throw NormalError("unexpected mesh type");
+                                    }
+                                  },
+                                  mesh_v->variant());
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("hyperelastic_dt", std::function(
 
                                                  [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index fc9d5703ed793f23c7f9768f4f2137585c3be318..4db29e203492c098a4807665ee50eb68b3be8287 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -9,6 +9,7 @@
 #include <mesh/Mesh.hpp>
 #include <output/GnuplotWriter.hpp>
 #include <output/GnuplotWriter1D.hpp>
+#include <output/GnuplotWriterRaw.hpp>
 #include <output/IWriter.hpp>
 #include <output/NamedDiscreteFunction.hpp>
 #include <output/NamedItemArrayVariant.hpp>
@@ -73,6 +74,23 @@ WriterModule::WriterModule()
 
                               ));
 
+  this->_addBuiltinFunction("gnuplot_raw_writer", std::function(
+
+                                                    [](const std::string& filename) -> std::shared_ptr<const IWriter> {
+                                                      return std::make_shared<GnuplotWriterRaw>(filename);
+                                                    }
+
+                                                    ));
+
+  this->_addBuiltinFunction("gnuplot_raw_writer",
+                            std::function(
+
+                              [](const std::string& filename, const double& period) -> std::shared_ptr<const IWriter> {
+                                return std::make_shared<GnuplotWriterRaw>(filename, period);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("name_output", std::function(
 
                                              [](std::shared_ptr<const DiscreteFunctionVariant> discrete_function,
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index 051e96480431aaa733059e2cf9c9518b29216812..fc975ba9b8234d8f70bd850a39bef58360e64108 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -4,6 +4,8 @@ add_library(
   PugsOutput
   GnuplotWriter.cpp
   GnuplotWriter1D.cpp
+  GnuplotWriterBase.cpp
+  GnuplotWriterRaw.cpp
   VTKWriter.cpp
   WriterBase.cpp)
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 991ec74537fffaaa2bef855f5c4af88f7acfa03f..83e2d8effaae57605d0bcf6f09f0a521c6524ff7 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -7,94 +7,14 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
+#include <utils/Demangle.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
 #include <utils/RevisionInfo.hpp>
 #include <utils/Stringify.hpp>
 
-#include <utils/Demangle.hpp>
-
 #include <fstream>
-#include <iomanip>
-
-std::string
-GnuplotWriter::_getDateAndVersionComment() const
-{
-  std::ostringstream os;
-
-  std::time_t now = std::time(nullptr);
-  os << "#  Generated by pugs: " << std::ctime(&now);
-  os << "#  version: " << RevisionInfo::version() << '\n';
-  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
-  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
-  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
-  os << '\n';
-
-  return os.str();
-}
-
-std::string
-GnuplotWriter::_getFilename() const
-{
-  std::ostringstream sout;
-  sout << m_base_filename;
-  if (m_period_manager.has_value()) {
-    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
-  }
-  sout << ".gnu";
-  return sout.str();
-}
-
-template <size_t Dimension>
-void
-GnuplotWriter::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_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;
-  for (const auto& i_named_item_data : output_named_item_data_set) {
-    const std::string name        = i_named_item_data.first;
-    const auto& item_data_variant = i_named_item_data.second;
-    std::visit(
-      [&](auto&& item_data) {
-        using ItemDataType = std::decay_t<decltype(item_data)>;
-        using DataType     = std::decay_t<typename ItemDataType::data_type>;
-        if constexpr (is_item_value_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            fout << ' ' << i++ << ':' << name;
-          } else if constexpr (is_tiny_vector_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else if constexpr (is_tiny_matrix_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
-              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
-                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
-              }
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else if constexpr (is_item_array_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else {
-          throw UnexpectedError("invalid ItemData type");
-        }
-      },
-      item_data_variant);
-  }
-  fout << "\n\n";
-}
 
 template <typename DataType>
 void
@@ -263,7 +183,7 @@ GnuplotWriter::_write(const MeshType& mesh,
     if (time.has_value()) {
       fout << "# time = " << *time << "\n\n";
     }
-    this->_writePreamble<MeshType::Dimension>(output_named_item_data_set, fout);
+    this->_writePreamble(MeshType::Dimension, output_named_item_data_set, true /*store coordinates*/, fout);
   }
 
   for (const auto& [name, item_data_variant] : output_named_item_data_set) {
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 3069b79da716d2640e8446822e23fcab37b98316..0c14204d930bc947832a484c617a9aa5965b04e6 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -1,7 +1,7 @@
 #ifndef GNUPLOT_WRITER_HPP
 #define GNUPLOT_WRITER_HPP
 
-#include <output/WriterBase.hpp>
+#include <output/GnuplotWriterBase.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -13,13 +13,9 @@ class MeshVariant;
 #include <optional>
 #include <string>
 
-class GnuplotWriter final : public WriterBase
+class GnuplotWriter final : public GnuplotWriterBase
 {
  private:
-  std::string _getDateAndVersionComment() const;
-
-  std::string _getFilename() const;
-
   template <typename DataType>
   void _writeCellData(const CellValue<DataType>& cell_value, CellId cell_id, std::ostream& fout) const;
 
@@ -32,9 +28,6 @@ class GnuplotWriter final : public WriterBase
   template <typename DataType>
   void _writeNodeData(const NodeArray<DataType>& node_value, NodeId cell_id, std::ostream& fout) const;
 
-  template <size_t Dimension>
-  void _writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const;
-
   template <typename ItemDataT>
   void _writeData(const ItemDataT& item_data,
                   [[maybe_unused]] CellId cell_id,
@@ -61,9 +54,11 @@ class GnuplotWriter final : public WriterBase
   void _writeMesh(const MeshVariant& mesh_v) const final;
 
  public:
-  GnuplotWriter(const std::string& base_filename) : WriterBase(base_filename) {}
+  GnuplotWriter(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
 
-  GnuplotWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
+  GnuplotWriter(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
+  {}
 
   ~GnuplotWriter() = default;
 };
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 84e8aa080c6802eb546a4d0183fb446201e31ac8..799a8d51d7dceafeea7f3555cbd7d73088d116c0 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -7,118 +7,13 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
+#include <utils/Demangle.hpp>
 #include <utils/Filesystem.hpp>
 #include <utils/Messenger.hpp>
 #include <utils/PugsTraits.hpp>
-#include <utils/RevisionInfo.hpp>
 #include <utils/Stringify.hpp>
 
-#include <utils/Demangle.hpp>
-
 #include <fstream>
-#include <iomanip>
-
-std::string
-GnuplotWriter1D::_getDateAndVersionComment() const
-{
-  std::ostringstream os;
-
-  std::time_t now = std::time(nullptr);
-  os << "#  Generated by pugs: " << std::ctime(&now);
-  os << "#  version: " << RevisionInfo::version() << '\n';
-  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
-  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
-  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
-  os << '\n';
-
-  return os.str();
-}
-
-std::string
-GnuplotWriter1D::_getFilename() const
-{
-  std::ostringstream sout;
-  sout << m_base_filename;
-  if (m_period_manager.has_value()) {
-    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
-  }
-  sout << ".gnu";
-  return sout.str();
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_cell_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::cell;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_face_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::face;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_edge_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::edge;
-}
-
-template <typename ItemDataType>
-bool
-GnuplotWriter1D::_is_node_data(const ItemDataType&) const
-{
-  return ItemDataType::item_t == ItemType::node;
-}
-
-void
-GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_data_set, std::ostream& fout) const
-{
-  fout << "# list of data\n";
-  fout << "# 1:x";
-  uint64_t i = 2;
-  for (const auto& i_named_item_data : output_named_item_data_set) {
-    const std::string name        = i_named_item_data.first;
-    const auto& item_data_variant = i_named_item_data.second;
-    std::visit(
-      [&](auto&& item_data) {
-        using ItemDataType = std::decay_t<decltype(item_data)>;
-        using DataType     = std::decay_t<typename ItemDataType::data_type>;
-        if constexpr (is_item_value_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            fout << ' ' << i++ << ':' << name;
-          } else if constexpr (is_tiny_vector_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else if constexpr (is_tiny_matrix_v<DataType>) {
-            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
-              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
-                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
-              }
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else if constexpr (is_item_array_v<ItemDataType>) {
-          if constexpr (std::is_arithmetic_v<DataType>) {
-            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
-              fout << ' ' << i++ << ':' << name << '[' << j << ']';
-            }
-          } else {
-            throw UnexpectedError("invalid data type");
-          }
-        } else {
-          throw UnexpectedError("invalid ItemData type");
-        }
-      },
-      item_data_variant);
-  }
-  fout << "\n\n";
-}
 
 template <typename DataType, ItemType item_type>
 size_t
@@ -339,7 +234,7 @@ GnuplotWriter1D::_write(const MeshType& mesh,
       fout << "# time = " << *time << "\n\n";
     }
 
-    _writePreamble(output_named_item_data_set, fout);
+    _writePreamble(MeshType::Dimension, output_named_item_data_set, true /*store coordinates*/, fout);
   }
 
   if (has_cell_data) {
@@ -355,7 +250,7 @@ GnuplotWriter1D::_writeMesh(const MeshVariant&) const
   std::ostringstream errorMsg;
   errorMsg << "gnuplot_1d_writer does not write meshes\n"
            << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
-           << "gnuplot_writer" << rang::fg::reset << "  instead";
+           << "gnuplot_writer" << rang::fg::reset << " instead";
   throw NormalError(errorMsg.str());
 }
 
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 66661918252856d8b6a069773c931c705a9b9330..32e4501873b175fb9af10123935d4666cd331404 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -1,7 +1,7 @@
 #ifndef GNUPLOT_WRITER_1D_HPP
 #define GNUPLOT_WRITER_1D_HPP
 
-#include <output/WriterBase.hpp>
+#include <output/GnuplotWriterBase.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -13,25 +13,9 @@ class MeshVariant;
 #include <optional>
 #include <string>
 
-class GnuplotWriter1D final : public WriterBase
+class GnuplotWriter1D final : public GnuplotWriterBase
 {
  private:
-  std::string _getDateAndVersionComment() const;
-
-  std::string _getFilename() const;
-
-  template <typename ItemDataType>
-  bool _is_cell_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_face_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_edge_data(const ItemDataType&) const;
-
-  template <typename ItemDataType>
-  bool _is_node_data(const ItemDataType&) const;
-
   template <typename DataType, ItemType item_type>
   size_t _itemDataNbRow(const ItemValue<DataType, item_type>&) const;
 
@@ -43,8 +27,6 @@ class GnuplotWriter1D final : public WriterBase
                        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 <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_value_set,
@@ -60,9 +42,10 @@ class GnuplotWriter1D final : public WriterBase
   void _writeMesh(const MeshVariant& mesh_v) const final;
 
  public:
-  GnuplotWriter1D(const std::string& base_filename) : WriterBase(base_filename) {}
+  GnuplotWriter1D(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
 
-  GnuplotWriter1D(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period)
+  GnuplotWriter1D(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
   {}
 
   ~GnuplotWriter1D() = default;
diff --git a/src/output/GnuplotWriterBase.cpp b/src/output/GnuplotWriterBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..110c03f77406c3e68e0ffb35b8528c5639948a94
--- /dev/null
+++ b/src/output/GnuplotWriterBase.cpp
@@ -0,0 +1,94 @@
+#include <output/GnuplotWriterBase.hpp>
+
+#include <output/OutputNamedItemValueSet.hpp>
+#include <utils/RevisionInfo.hpp>
+
+#include <ctime>
+#include <iomanip>
+#include <sstream>
+
+std::string
+GnuplotWriterBase::_getDateAndVersionComment() const
+{
+  std::ostringstream os;
+
+  std::time_t now = std::time(nullptr);
+  os << "#  Generated by pugs: " << std::ctime(&now);
+  os << "#  version: " << RevisionInfo::version() << '\n';
+  os << "#  tag:  " << RevisionInfo::gitTag() << '\n';
+  os << "#  HEAD: " << RevisionInfo::gitHead() << '\n';
+  os << "#  hash: " << RevisionInfo::gitHash() << " (" << ((RevisionInfo::gitIsClean()) ? "clean" : "dirty") << ")\n";
+  os << '\n';
+
+  return os.str();
+}
+
+std::string
+GnuplotWriterBase::_getFilename() const
+{
+  std::ostringstream sout;
+  sout << m_base_filename;
+  if (m_period_manager.has_value()) {
+    sout << '.' << std::setfill('0') << std::setw(4) << m_period_manager->nbSavedTimes();
+  }
+  sout << ".gnu";
+  return sout.str();
+}
+
+void
+GnuplotWriterBase::_writePreamble(const size_t& dimension,
+                                  const OutputNamedItemDataSet& output_named_item_data_set,
+                                  const bool& store_coordinates,
+                                  std::ostream& fout) const
+{
+  fout << "# list of data\n";
+  uint64_t i = 0;
+
+  fout << "#";
+  if (store_coordinates) {
+    fout << " 1:x";
+    if (dimension > 1) {
+      fout << " 2:y";
+    }
+    i = dimension;
+  }
+
+  for (const auto& i_named_item_data : output_named_item_data_set) {
+    const std::string name        = i_named_item_data.first;
+    const auto& item_data_variant = i_named_item_data.second;
+    std::visit(
+      [&](auto&& item_data) {
+        using ItemDataType = std::decay_t<decltype(item_data)>;
+        using DataType     = std::decay_t<typename ItemDataType::data_type>;
+        if constexpr (is_item_value_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            fout << ' ' << i++ << ':' << name;
+          } else if constexpr (is_tiny_vector_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.dimension(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else if constexpr (is_tiny_matrix_v<DataType>) {
+            for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
+              for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
+                fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+              }
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
+          }
+        } else if constexpr (is_item_array_v<ItemDataType>) {
+          if constexpr (std::is_arithmetic_v<DataType>) {
+            for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+              fout << ' ' << i++ << ':' << name << '[' << j << ']';
+            }
+          } else {
+            throw UnexpectedError("invalid data type");
+          }
+        } else {
+          throw UnexpectedError("invalid ItemData type");
+        }
+      },
+      item_data_variant);
+  }
+  fout << "\n\n";
+}
diff --git a/src/output/GnuplotWriterBase.hpp b/src/output/GnuplotWriterBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a99390facb8d2bd67cc323d676d305f26401dd6e
--- /dev/null
+++ b/src/output/GnuplotWriterBase.hpp
@@ -0,0 +1,58 @@
+#ifndef GNUPLOT_WRITER_BASE_HPP
+#define GNUPLOT_WRITER_BASE_HPP
+
+#include <output/WriterBase.hpp>
+
+#include <utils/PugsMacros.hpp>
+
+class GnuplotWriterBase : public WriterBase
+{
+ protected:
+  std::string _getDateAndVersionComment() const;
+
+  std::string _getFilename() const;
+
+  void _writePreamble(const size_t& dimension,
+                      const OutputNamedItemDataSet& output_named_item_data_set,
+                      const bool& store_coordinates,
+                      std::ostream& fout) const;
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_cell_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::cell;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_face_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::face;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_edge_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::edge;
+  }
+
+  template <typename ItemDataType>
+  PUGS_INLINE bool
+  _is_node_data(const ItemDataType&) const
+  {
+    return ItemDataType::item_t == ItemType::node;
+  }
+
+ public:
+  GnuplotWriterBase(const std::string& base_filename, const double& time_period)
+    : WriterBase(base_filename, time_period)
+  {}
+
+  GnuplotWriterBase(const std::string& base_filename) : WriterBase(base_filename) {}
+
+  virtual ~GnuplotWriterBase() = default;
+};
+
+#endif   // GNUPLOT_WRITER_BASE_HPP
diff --git a/src/output/GnuplotWriterRaw.cpp b/src/output/GnuplotWriterRaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac93ba0029f3d73583c98d62c85beb7901f3569f
--- /dev/null
+++ b/src/output/GnuplotWriterRaw.cpp
@@ -0,0 +1,269 @@
+#include <output/GnuplotWriterRaw.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <utils/Filesystem.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/PugsTraits.hpp>
+#include <utils/RevisionInfo.hpp>
+#include <utils/Stringify.hpp>
+
+#include <utils/Demangle.hpp>
+
+#include <fstream>
+#include <iomanip>
+
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriterRaw::_itemDataNbRow(const ItemValue<DataType, item_type>&) const
+{
+  if constexpr (std::is_arithmetic_v<DataType>) {
+    return 1;
+  } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
+    return DataType::Dimension;
+  } else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
+    return DataType{}.dimension();
+  } else {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+  }
+}
+
+template <typename DataType, ItemType item_type>
+size_t
+GnuplotWriterRaw::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array) const
+{
+  return item_array.sizeOfArrays();
+}
+
+template <MeshConcept MeshType, ItemType item_type>
+void
+GnuplotWriterRaw::_writeItemDatas(const MeshType& mesh,
+                                  const OutputNamedItemDataSet& output_named_item_data_set,
+                                  std::ostream& fout) const
+{
+  using ItemId = ItemIdT<item_type>;
+
+  const size_t number_of_columns = [&] {
+    size_t nb_columns = 0;
+    for (const auto& [name, item_data] : output_named_item_data_set) {
+      std::visit([&](auto&& value) { nb_columns += _itemDataNbRow(value); }, item_data);
+    }
+    return nb_columns;
+  }();
+
+  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) {
+        if (is_owned[item_id]) {
+          ++number_of_owned_items;
+        }
+      }
+
+      return number_of_owned_items;
+    } else {
+      return mesh.template numberOf<item_type>();
+    }
+  }();
+
+  Array<double> values{number_of_columns * number_of_owned_lines};
+
+  size_t column_number = 0;
+  for (const auto& [name, output_item_data] : output_named_item_data_set) {
+    std::visit(
+      [&](auto&& item_data) {
+        using ItemDataT = std::decay_t<decltype(item_data)>;
+        if constexpr (ItemDataT::item_t == item_type) {
+          if constexpr (is_item_value_v<ItemDataT>) {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  values[number_of_columns * index + column_number] = item_data[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_data[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{}.numberOfRows(); ++i) {
+                    for (size_t j = 0; j < DataT{}.numberOfColumns(); ++j) {
+                      values[k++] = item_data[item_id](i, j);
+                    }
+                  }
+                }
+                ++index;
+              }
+            }
+          } else {
+            using DataT  = std::decay_t<typename ItemDataT::data_type>;
+            size_t index = 0;
+            for (ItemId item_id = 0; item_id < item_data.numberOfItems(); ++item_id) {
+              if (is_owned[item_id]) {
+                if constexpr (std::is_arithmetic_v<DataT>) {
+                  const size_t k = number_of_columns * index + column_number;
+                  for (size_t j = 0; j < item_data.sizeOfArrays(); ++j) {
+                    values[k + j] = item_data[item_id][j];
+                  }
+                }
+                ++index;
+              }
+            }
+          }
+        }
+        column_number += _itemDataNbRow(item_data);
+      },
+      output_item_data);
+  }
+
+  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;
+    }
+
+    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';
+    }
+  }
+}
+
+template <MeshConcept MeshType>
+void
+GnuplotWriterRaw::_write(const MeshType& mesh,
+                         const OutputNamedItemDataSet& output_named_item_data_set,
+                         std::optional<double> time) const
+{
+  bool has_cell_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_cell_data |= std::visit([&](auto&& item_data) { return this->_is_cell_data(item_data); }, item_data_variant);
+  }
+
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    std::visit(
+      [&, name = name](auto&& item_data) {
+        if (this->_is_face_data(item_data)) {
+          std::ostringstream error_msg;
+          error_msg << "gnuplot_raw_writer does not support face data, cannot save variable \"" << rang::fgB::yellow
+                    << name << rang::fg::reset << '"';
+          throw NormalError(error_msg.str());
+        }
+      },
+      item_data_variant);
+  }
+
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    std::visit(
+      [&, name = name](auto&& item_data) {
+        if (this->_is_edge_data(item_data)) {
+          std::ostringstream error_msg;
+          error_msg << "gnuplot_1d_writer does not support edge data, cannot save variable \"" << rang::fgB::yellow
+                    << name << rang::fg::reset << '"';
+          throw NormalError(error_msg.str());
+        }
+      },
+      item_data_variant);
+  }
+
+  bool has_node_data = false;
+  for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+    has_node_data |= std::visit([&](auto&& item_data) { return this->_is_node_data(item_data); }, item_data_variant);
+  }
+
+  if (has_cell_data and has_node_data) {
+    throw NormalError("cannot store both node and cell data in the same gnuplot file");
+  }
+
+  std::ofstream fout;
+
+  if (parallel::rank() == 0) {
+    fout.open(_getFilename());
+    if (not fout) {
+      std::ostringstream error_msg;
+      error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilename() << rang::fg::reset << '"';
+      throw NormalError(error_msg.str());
+    }
+
+    fout.precision(15);
+    fout.setf(std::ios_base::scientific);
+    fout << _getDateAndVersionComment();
+
+    if (time.has_value()) {
+      fout << "# time = " << *time << "\n\n";
+    }
+
+    _writePreamble(MeshType::Dimension, output_named_item_data_set, false /*do not store coordinates*/, fout);
+  }
+
+  if (has_cell_data) {
+    this->_writeItemDatas<MeshType, ItemType::cell>(mesh, output_named_item_data_set, fout);
+  } else {   // has_node_value
+    this->_writeItemDatas<MeshType, ItemType::node>(mesh, output_named_item_data_set, fout);
+  }
+}
+
+void
+GnuplotWriterRaw::_writeMesh(const MeshVariant&) const
+{
+  std::ostringstream errorMsg;
+  errorMsg << "gnuplot_raw_writer does not write meshes\n"
+           << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
+           << "gnuplot_writer" << rang::fg::reset << " instead";
+  throw NormalError(errorMsg.str());
+}
+
+void
+GnuplotWriterRaw::_writeAtTime(const MeshVariant& mesh_v,
+                               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);
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      if constexpr (MeshType::Dimension == 1) {
+        this->_write(*p_mesh, output_named_item_data_set, time);
+      } else if constexpr (MeshType::Dimension == 2) {
+        std::ostringstream errorMsg;
+        errorMsg << "gnuplot_raw_writer is not available in dimension " << stringify(MeshType::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());
+      } else {
+        throw NormalError("gnuplot format is not available in dimension " + stringify(MeshType::Dimension));
+      }
+    },
+    mesh_v.variant());
+}
+
+void
+GnuplotWriterRaw::_write(const MeshVariant& mesh_v,
+                         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);
+
+  std::visit([&](auto&& p_mesh) { this->_write(*p_mesh, output_named_item_data_set, {}); }, mesh_v.variant());
+}
diff --git a/src/output/GnuplotWriterRaw.hpp b/src/output/GnuplotWriterRaw.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..481bb2a4bfc548c89ea684ec7a841df67bede57b
--- /dev/null
+++ b/src/output/GnuplotWriterRaw.hpp
@@ -0,0 +1,54 @@
+#ifndef GNUPLOT_WRITER_RAW_HPP
+#define GNUPLOT_WRITER_RAW_HPP
+
+#include <output/GnuplotWriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+
+class MeshVariant;
+
+#include <optional>
+#include <string>
+
+class GnuplotWriterRaw final : public GnuplotWriterBase
+{
+ private:
+  template <typename DataType, ItemType item_type>
+  size_t _itemDataNbRow(const ItemValue<DataType, item_type>&) const;
+
+  template <typename DataType, ItemType item_type>
+  size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
+
+  template <MeshConcept MeshType, ItemType item_type>
+  void _writeItemDatas(const MeshType& mesh,
+                       const OutputNamedItemDataSet& output_named_item_data_set,
+                       std::ostream& fout) const;
+
+  template <MeshConcept MeshType>
+  void _write(const MeshType& mesh,
+              const OutputNamedItemDataSet& output_named_item_value_set,
+              std::optional<double> time) const;
+
+  void _writeAtTime(const MeshVariant& mesh_v,
+                    const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+                    double time) const final;
+
+  void _write(const MeshVariant& mesh_v,
+              const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
+
+  void _writeMesh(const MeshVariant& mesh_v) const final;
+
+ public:
+  GnuplotWriterRaw(const std::string& base_filename) : GnuplotWriterBase(base_filename) {}
+
+  GnuplotWriterRaw(const std::string& base_filename, const double time_period)
+    : GnuplotWriterBase(base_filename, time_period)
+  {}
+
+  ~GnuplotWriterRaw() = default;
+};
+
+#endif   // GNUPLOT_WRITER_RAW_HPP