diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
index a2a456f0ebc91476b2dbb541e9fdc6eecc77c83c..edce9e025f0feb4cd39acaadcb537451c9251065 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -1,7 +1,6 @@
 #include <language/algorithms/AcousticSolverAlgorithm.hpp>
 
 #include <language/utils/InterpolateItemValue.hpp>
-#include <output/VTKWriter.hpp>
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryDescriptor.hpp>
@@ -178,18 +177,9 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
       mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
   }
 
-  VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
-
   while ((t < tmax) and (iteration < itermax)) {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
-    vtk_writer.write(mesh,
-                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
-                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
-                     t);
-
     double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
     if (t + dt > tmax) {
       dt = tmax - t;
@@ -210,16 +200,6 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
   std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
             << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
             << rang::style::reset << '\n';
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-    vtk_writer.write(mesh,
-                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
-                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
-                     t, true);   // forces last output
-  }
 }
 
 template AcousticSolverAlgorithm<1>::AcousticSolverAlgorithm(
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index c447ed188ce794aeea9ddec00775254d405e211b..611952939906e0afce5ea4682428ec1eefd36eaa 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -5,6 +5,7 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
+#include <output/GnuplotWriter.hpp>
 #include <output/IWriter.hpp>
 #include <output/NamedDiscreteFunction.hpp>
 #include <output/VTKWriter.hpp>
@@ -12,113 +13,6 @@
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
-template <size_t Dimension, typename DataType>
-void
-WriterModule::registerDiscreteFunctionP0(const std::string& name,
-                                         const IDiscreteFunction& i_discrete_function,
-                                         OutputNamedItemValueSet& named_item_value_set)
-{
-  const DiscreteFunctionP0<Dimension, DataType>& discrete_function =
-    dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(i_discrete_function);
-  named_item_value_set.add(NamedItemValue{name, discrete_function.cellValues()});
-}
-
-template <size_t Dimension>
-void
-WriterModule::registerDiscreteFunctionP0(const std::string& name,
-                                         const IDiscreteFunction& i_discrete_function,
-                                         OutputNamedItemValueSet& named_item_value_set)
-{
-  const ASTNodeDataType& data_type = i_discrete_function.dataType();
-  switch (data_type) {
-  case ASTNodeDataType::bool_t: {
-    registerDiscreteFunctionP0<Dimension, bool>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case ASTNodeDataType::unsigned_int_t: {
-    registerDiscreteFunctionP0<Dimension, uint64_t>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case ASTNodeDataType::int_t: {
-    registerDiscreteFunctionP0<Dimension, int64_t>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case ASTNodeDataType::double_t: {
-    registerDiscreteFunctionP0<Dimension, double>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case ASTNodeDataType::vector_t: {
-    switch (data_type.dimension()) {
-    case 1: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<1, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 2: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<2, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 3: {
-      registerDiscreteFunctionP0<Dimension, TinyVector<3, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    default: {
-      throw UnexpectedError("invalid vector dimension");
-    }
-    }
-    break;
-  }
-  case ASTNodeDataType::matrix_t: {
-    Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
-    switch (data_type.nbRows()) {
-    case 1: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<1, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 2: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<2, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    case 3: {
-      registerDiscreteFunctionP0<Dimension, TinyMatrix<3, double>>(name, i_discrete_function, named_item_value_set);
-      break;
-    }
-    default: {
-      throw UnexpectedError("invalid matrix dimension");
-    }
-    }
-    break;
-  }
-  default: {
-    throw UnexpectedError("invalid data type " + dataTypeName(data_type));
-  }
-  }
-}
-
-void
-WriterModule::registerDiscreteFunctionP0(const NamedDiscreteFunction& named_discrete_function,
-                                         OutputNamedItemValueSet& named_item_value_set)
-{
-  const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
-  const std::string& name                      = named_discrete_function.name();
-  switch (i_discrete_function.mesh()->dimension()) {
-  case 1: {
-    registerDiscreteFunctionP0<1>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case 2: {
-    registerDiscreteFunctionP0<2>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  case 3: {
-    registerDiscreteFunctionP0<3>(name, i_discrete_function, named_item_value_set);
-    break;
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
 WriterModule::WriterModule()
 {
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>>);
@@ -135,6 +29,16 @@ WriterModule::WriterModule()
 
                               ));
 
+  this->_addBuiltinFunction("gnuplot_writer",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IWriter>(const std::string&,
+                                                                                                    const double&)>>(
+
+                              [](const std::string& filename, const double& period) {
+                                return std::make_shared<GnuplotWriter>(filename, period);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("name_output",
                             std::make_shared<BuiltinFunctionEmbedder<
                               std::shared_ptr<const NamedDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
@@ -152,42 +56,7 @@ WriterModule::WriterModule()
                                                                           std::shared_ptr<const IMesh>)>>(
 
                               [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> p_mesh) -> void {
-                                static double number = 0;
-
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  const std::shared_ptr mesh =
-                                    std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(p_mesh);
-
-                                  writer->write(mesh,
-                                                OutputNamedItemValueSet{
-                                                  NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
-                                                number, true);
-                                  break;
-                                }
-                                case 2: {
-                                  const std::shared_ptr mesh =
-                                    std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(p_mesh);
-
-                                  writer->write(mesh,
-                                                OutputNamedItemValueSet{
-                                                  NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
-                                                number, true);
-                                  break;
-                                }
-                                case 3: {
-                                  const std::shared_ptr mesh =
-                                    std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(p_mesh);
-
-                                  writer->write(mesh,
-                                                OutputNamedItemValueSet{
-                                                  NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
-                                                number, true);
-                                  break;
-                                }
-                                }
-
-                                number++;
+                                writer->writeMesh(p_mesh);
                               }
 
                               ));
@@ -201,70 +70,21 @@ WriterModule::WriterModule()
                                  const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
                                    named_discrete_function_list,
                                  const double& time) -> void {
-                                Assert(named_discrete_function_list.size() > 0);
-
-                                std::shared_ptr p_mesh = named_discrete_function_list[0]->discreteFunction()->mesh();
-                                for (size_t i = 1; i < named_discrete_function_list.size(); ++i) {
-                                  if (p_mesh != named_discrete_function_list[i]->discreteFunction()->mesh()) {
-                                    std::ostringstream error_msg;
-                                    error_msg << "discrete functions must be defined on the same mesh!\n"
-                                              << rang::fgB::yellow << "note:" << rang::style::reset << " cannot write "
-                                              << rang::fgB::blue << named_discrete_function_list[0]->name()
-                                              << rang::style::reset << " and " << rang::fgB::blue
-                                              << named_discrete_function_list[i]->name() << rang::style::reset
-                                              << " in the same file.";
-                                    throw NormalError(error_msg.str());
-                                  }
-                                }
-
-                                OutputNamedItemValueSet named_item_value_set;
-
-                                for (auto& named_discrete_function : named_discrete_function_list) {
-                                  const IDiscreteFunction& i_discrete_function =
-                                    *named_discrete_function->discreteFunction();
-
-                                  switch (i_discrete_function.descriptor().type()) {
-                                  case DiscreteFunctionType::P0: {
-                                    WriterModule::registerDiscreteFunctionP0(*named_discrete_function,
-                                                                             named_item_value_set);
-                                    break;
-                                  }
-                                  default: {
-                                    std::ostringstream error_msg;
-                                    error_msg << "the type of discrete function of " << rang::fgB::blue
-                                              << named_discrete_function->name() << rang::style::reset
-                                              << " is not supported";
-                                    throw NormalError(error_msg.str());
-                                  }
-                                  }
-                                }
-
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  using MeshType = Mesh<Connectivity<1>>;
-                                  const std::shared_ptr<const MeshType> mesh =
-                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+                                writer->writeIfNeeded(named_discrete_function_list, time);
+                              }
 
-                                  writer->write(mesh, named_item_value_set, time, true);
-                                  break;
-                                }
-                                case 2: {
-                                  using MeshType = Mesh<Connectivity<2>>;
-                                  const std::shared_ptr<const MeshType> mesh =
-                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+                              ));
 
-                                  writer->write(mesh, named_item_value_set, time, true);
-                                  break;
-                                }
-                                case 3: {
-                                  using MeshType = Mesh<Connectivity<3>>;
-                                  const std::shared_ptr<const MeshType> mesh =
-                                    std::dynamic_pointer_cast<const MeshType>(p_mesh);
+  this->_addBuiltinFunction("force_write",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IWriter>,
+                                   const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&, const double&)>>(
 
-                                  writer->write(mesh, named_item_value_set, time, true);
-                                  break;
-                                }
-                                }
+                              [](std::shared_ptr<const IWriter> writer,
+                                 const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+                                   named_discrete_function_list,
+                                 const double& time) -> void {
+                                writer->writeForced(named_discrete_function_list, time);
                               }
 
                               ));
diff --git a/src/language/modules/WriterModule.hpp b/src/language/modules/WriterModule.hpp
index 08f17851e15153baa4b8cdcfb9dcec25fa7c984e..14a961d7a812b34aac1c5010ddda674749fdbfb8 100644
--- a/src/language/modules/WriterModule.hpp
+++ b/src/language/modules/WriterModule.hpp
@@ -22,15 +22,6 @@ inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IWriter>> =
 
 class WriterModule : public BuiltinModule
 {
- private:
-  template <size_t Dimension, typename DataType>
-  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
-
-  template <size_t Dimension>
-  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
-
-  static void registerDiscreteFunctionP0(const NamedDiscreteFunction&, OutputNamedItemValueSet&);
-
  public:
   std::string_view
   name() const final
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index e5b860bc0bf2781a16d7453537e6657600d55250..a3e6799657d2ed422fee5ab7018939ea080cf4c2 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -2,7 +2,9 @@
 
 add_library(
   PugsOutput
-  VTKWriter.cpp)
+  GnuplotWriter.cpp
+  VTKWriter.cpp
+  WriterBase.cpp)
 
 # ------------------- Installation --------------------
 # temporary version workaround
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1fcfd4dc0165b4eec1a15ae366903bee82f4879
--- /dev/null
+++ b/src/output/GnuplotWriter.cpp
@@ -0,0 +1,332 @@
+#include <output/GnuplotWriter.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <utils/Messenger.hpp>
+#include <utils/PugsTraits.hpp>
+#include <utils/RevisionInfo.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;
+  sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".gnu";
+  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
+{
+  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_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;
+    std::visit(
+      [&](auto&& item_value) {
+        using ItemValueType = std::decay_t<decltype(item_value)>;
+        using DataType      = std::decay_t<typename ItemValueType::data_type>;
+        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{}.nbRows(); ++j) {
+            for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+              fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
+            }
+          }
+        } else {
+          throw UnexpectedError("invalid data type");
+        }
+      },
+      item_value_variant);
+  }
+  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
+{
+  const auto& value = cell_value[cell_id];
+  if constexpr (std::is_arithmetic_v<DataType>) {
+    fout << ' ' << value;
+  } else if constexpr (is_tiny_vector_v<std::decay_t<DataType>>) {
+    for (size_t i = 0; i < value.dimension(); ++i) {
+      fout << ' ' << value[i];
+    }
+  } 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);
+      }
+    }
+  } else {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<DataType>());
+  }
+}
+
+template <typename MeshType>
+void
+GnuplotWriter::_writeCellValues(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);
+
+    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) {
+    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);
+
+      for (auto [name, item_value] : output_named_item_value_set) {
+        std::visit([&](auto&& cell_value) { _writeCellValue(cell_value, cell_id, os); }, 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';
+    }
+  } else {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+}
+
+template <typename ValueType>
+void
+GnuplotWriter::_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];
+    }
+  } 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 {
+    throw UnexpectedError("invalid data type for cell value output: " + demangle<ValueType>());
+  }
+}
+
+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);
+  }
+
+  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);
+  }
+
+  if (has_cell_value and has_node_value) {
+    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();
+
+  fout << "\n# time = " << time << "\n\n";
+
+  _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);
+  }
+}
+
+void
+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: {
+    _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);
+    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);
+    break;
+  }
+  default: {
+    throw NormalError("gnuplot format is not available in dimension " + std::to_string(p_mesh->dimension()));
+  }
+  }
+}
+
+void
+GnuplotWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                     double time) const
+{
+  std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
+
+  OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+
+  switch (mesh->dimension()) {
+  case 1: {
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
+    break;
+  }
+  case 2: {
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
+    break;
+  }
+  default: {
+    throw NormalError("gnuplot format is not available in dimension " + std::to_string(mesh->dimension()));
+  }
+  }
+}
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0006d5ccfde75c50f0ca8d73c0479fa2d6d6f00e
--- /dev/null
+++ b/src/output/GnuplotWriter.hpp
@@ -0,0 +1,74 @@
+#ifndef GNUPLOT_WRITER_HPP
+#define GNUPLOT_WRITER_HPP
+
+#include <output/WriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+
+class IMesh;
+
+#include <string>
+
+class GnuplotWriter : public WriterBase
+{
+ private:
+  std::string _getDateAndVersionComment() const;
+
+  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 MeshType>
+  void _write(const std::shared_ptr<const MeshType>& mesh,
+              const OutputNamedItemValueSet& output_named_item_value_set,
+              double time) const;
+
+ public:
+  void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+
+  void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+             double time) const final;
+
+  GnuplotWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
+
+  ~GnuplotWriter() = default;
+};
+
+#endif   // GNUPLOT_WRITER_HPP
diff --git a/src/output/IWriter.hpp b/src/output/IWriter.hpp
index 98df0858ee121db43119ba085e783ae757c3022d..02317ad98b055567133f9b969476fedbd2a5c0e1 100644
--- a/src/output/IWriter.hpp
+++ b/src/output/IWriter.hpp
@@ -1,16 +1,25 @@
 #ifndef I_WRITER_HPP
 #define I_WRITER_HPP
 
-#include <output/OutputNamedItemValueSet.hpp>
+#include <output/NamedDiscreteFunction.hpp>
+
+#include <memory>
+#include <vector>
+
 class IMesh;
 
 class IWriter
 {
  public:
-  virtual void write(const std::shared_ptr<const IMesh>& mesh,
-                     const OutputNamedItemValueSet& output_named_item_value_set,
-                     double time,
-                     bool forced_output = false) const = 0;
+  virtual void writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
+
+  virtual void writeIfNeeded(
+    const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+    double time) const = 0;
+
+  virtual void writeForced(
+    const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+    double time) const = 0;
 
   IWriter()          = default;
   virtual ~IWriter() = default;
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index a0b4921abba97cbe2ff6e08be2358114d91d3e84..48f81dd376bd96c36820a535b2ff2fa20a4e831b 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -16,7 +16,7 @@ VTKWriter::_getFilenamePVTU() const
 {
   std::ostringstream sout;
   sout << m_base_filename;
-  sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".pvtu";
+  sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".pvtu";
   return sout.str();
 }
 
@@ -44,7 +44,7 @@ VTKWriter::_getFilenameVTU(int rank_number) const
   if (parallel::size() > 1) {
     sout << '-' << std::setfill('0') << std::setw(4) << rank_number;
   }
-  sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".vtu";
+  sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".vtu";
   return sout.str();
 }
 
@@ -273,22 +273,13 @@ template <typename MeshType>
 void
 VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
                   const OutputNamedItemValueSet& given_output_named_item_value_set,
-                  double time,
-                  bool forced_output) const
+                  double time) const
 {
   OutputNamedItemValueSet output_named_item_value_set{given_output_named_item_value_set};
   // Adding basic mesh information
   output_named_item_value_set.add(NamedItemValue{"cell_number", mesh->connectivity().cellNumber()});
   output_named_item_value_set.add(NamedItemValue{"node_number", mesh->connectivity().nodeNumber()});
 
-  if (time == m_last_time)
-    return;   // output already performed
-  if ((time - m_last_time >= m_time_period) or forced_output) {
-    m_last_time = time;
-  } else {
-    return;
-  }
-
   if (parallel::rank() == 0) {   // write PVTK file
     std::ofstream fout(_getFilenamePVTU());
     fout << "<?xml version=\"1.0\"?>\n";
@@ -519,29 +510,72 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << "</UnstructuredGrid>\n";
     fout << "</VTKFile>\n";
   }
-  m_file_number++;
+
+  if (parallel::rank() == 0) {   // write PVD file
+    std::ofstream fout(m_base_filename + ".pvd");
+
+    fout << "<?xml version=\"1.0\"?>\n";
+    fout << "<VTKFile type=\"Collection\" version=\"0.1\">\n";
+
+    fout << "<Collection>\n";
+    for (size_t i_time = 0; i_time < m_saved_times.size(); ++i_time) {
+      std::ostringstream sout;
+      sout << m_base_filename;
+      sout << '.' << std::setfill('0') << std::setw(4) << i_time << ".pvtu";
+
+      fout << "<DataSet timestep=\"" << m_saved_times[i_time] << "\" "
+           << "file=\"" << sout.str() << "\"/>\n";
+    }
+    fout << "<DataSet timestep=\"" << time << "\" group=\"\" part=\"0\" file=\"" << _getFilenamePVTU() << "\"/>\n";
+
+    fout << "</Collection>\n";
+    fout << "</VTKFile>\n";
+  }
+}
+
+void
+VTKWriter::writeMesh(const std::shared_ptr<const IMesh>& mesh) const
+{
+  OutputNamedItemValueSet output_named_item_value_set;
+
+  switch (mesh->dimension()) {
+  case 1: {
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, 0);
+    break;
+  }
+  case 2: {
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, 0);
+    break;
+  }
+  case 3: {
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, 0);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
 }
 
 void
-VTKWriter::write(const std::shared_ptr<const IMesh>& mesh,
-                 const OutputNamedItemValueSet& output_named_item_value_set,
-                 double time,
-                 bool forced_output) const
+VTKWriter::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                 double time) const
 {
+  std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
+
+  OutputNamedItemValueSet output_named_item_value_set = this->_getOutputNamedItemValueSet(named_discrete_function_list);
+
   switch (mesh->dimension()) {
   case 1: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time,
-                 forced_output);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), output_named_item_value_set, time);
     break;
   }
   case 2: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time,
-                 forced_output);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), output_named_item_value_set, time);
     break;
   }
   case 3: {
-    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, time,
-                 forced_output);
+    this->_write(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), output_named_item_value_set, time);
     break;
   }
   default: {
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 52a4bd9d94918d01cad6fccbced121688883baad..704ba63f990ead2b6b8e8f64cafdda627c2be64b 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -1,7 +1,7 @@
 #ifndef VTK_WRITER_HPP
 #define VTK_WRITER_HPP
 
-#include <output/IWriter.hpp>
+#include <output/WriterBase.hpp>
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
@@ -11,14 +11,9 @@ class IMesh;
 
 #include <string>
 
-class VTKWriter : public IWriter
+class VTKWriter : public WriterBase
 {
  private:
-  const std::string m_base_filename;
-  mutable unsigned int m_file_number;
-  mutable double m_last_time;
-  const double m_time_period;
-
   std::string _getDateAndVersionComment() const;
 
   std::string _getFilenamePVTU() const;
@@ -101,21 +96,15 @@ class VTKWriter : public IWriter
   template <typename MeshType>
   void _write(const std::shared_ptr<const MeshType>& mesh,
               const OutputNamedItemValueSet& output_named_item_value_set,
-              double time,
-              bool forced_output = false) const;
+              double time) const;
 
  public:
-  void write(const std::shared_ptr<const IMesh>& mesh,
-             const OutputNamedItemValueSet& output_named_item_value_set,
-             double time,
-             bool forced_output = false) const final;
-
-  VTKWriter(const std::string& base_filename, const double time_period)
-    : m_base_filename(base_filename),
-      m_file_number(0),
-      m_last_time(-std::numeric_limits<double>::max()),
-      m_time_period(time_period)
-  {}
+  void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
+
+  void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+             double time) const final;
+
+  VTKWriter(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period) {}
 
   ~VTKWriter() = default;
 };
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0359cf7b3d33c8aa9fb33f8b62c671425e5958ba
--- /dev/null
+++ b/src/output/WriterBase.cpp
@@ -0,0 +1,201 @@
+#include <output/WriterBase.hpp>
+
+#include <mesh/IMesh.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <utils/Exceptions.hpp>
+
+template <size_t Dimension, typename DataType>
+void
+WriterBase::registerDiscreteFunctionP0(const std::string& name,
+                                       const IDiscreteFunction& i_discrete_function,
+                                       OutputNamedItemValueSet& named_item_value_set)
+{
+  const DiscreteFunctionP0<Dimension, DataType>& discrete_function =
+    dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(i_discrete_function);
+  named_item_value_set.add(NamedItemValue{name, discrete_function.cellValues()});
+}
+
+template <size_t Dimension>
+void
+WriterBase::registerDiscreteFunctionP0(const std::string& name,
+                                       const IDiscreteFunction& i_discrete_function,
+                                       OutputNamedItemValueSet& named_item_value_set)
+{
+  const ASTNodeDataType& data_type = i_discrete_function.dataType();
+  switch (data_type) {
+  case ASTNodeDataType::bool_t: {
+    registerDiscreteFunctionP0<Dimension, bool>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case ASTNodeDataType::unsigned_int_t: {
+    registerDiscreteFunctionP0<Dimension, uint64_t>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case ASTNodeDataType::int_t: {
+    registerDiscreteFunctionP0<Dimension, int64_t>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case ASTNodeDataType::double_t: {
+    registerDiscreteFunctionP0<Dimension, double>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case ASTNodeDataType::vector_t: {
+    switch (data_type.dimension()) {
+    case 1: {
+      registerDiscreteFunctionP0<Dimension, TinyVector<1, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    case 2: {
+      registerDiscreteFunctionP0<Dimension, TinyVector<2, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    case 3: {
+      registerDiscreteFunctionP0<Dimension, TinyVector<3, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    default: {
+      throw UnexpectedError("invalid vector dimension");
+    }
+    }
+    break;
+  }
+  case ASTNodeDataType::matrix_t: {
+    Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
+    switch (data_type.nbRows()) {
+    case 1: {
+      registerDiscreteFunctionP0<Dimension, TinyMatrix<1, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    case 2: {
+      registerDiscreteFunctionP0<Dimension, TinyMatrix<2, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    case 3: {
+      registerDiscreteFunctionP0<Dimension, TinyMatrix<3, double>>(name, i_discrete_function, named_item_value_set);
+      break;
+    }
+    default: {
+      throw UnexpectedError("invalid matrix dimension");
+    }
+    }
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid data type " + dataTypeName(data_type));
+  }
+  }
+}
+
+void
+WriterBase::registerDiscreteFunctionP0(const NamedDiscreteFunction& named_discrete_function,
+                                       OutputNamedItemValueSet& named_item_value_set)
+{
+  const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
+  const std::string& name                      = named_discrete_function.name();
+  switch (i_discrete_function.mesh()->dimension()) {
+  case 1: {
+    registerDiscreteFunctionP0<1>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case 2: {
+    registerDiscreteFunctionP0<2>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  case 3: {
+    registerDiscreteFunctionP0<3>(name, i_discrete_function, named_item_value_set);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+WriterBase::_getMesh(
+  const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+{
+  Assert(named_discrete_function_list.size() > 0);
+
+  std::shared_ptr mesh = named_discrete_function_list[0]->discreteFunction()->mesh();
+  for (size_t i = 1; i < named_discrete_function_list.size(); ++i) {
+    if (mesh != named_discrete_function_list[i]->discreteFunction()->mesh()) {
+      std::ostringstream error_msg;
+      error_msg << "discrete functions must be defined on the same mesh!\n"
+                << rang::fgB::yellow << "note:" << rang::style::reset << " cannot write " << rang::fgB::blue
+                << named_discrete_function_list[0]->name() << rang::style::reset << " and " << rang::fgB::blue
+                << named_discrete_function_list[i]->name() << rang::style::reset << " in the same file.";
+      throw NormalError(error_msg.str());
+    }
+  }
+
+  return mesh;
+}
+
+OutputNamedItemValueSet
+WriterBase::_getOutputNamedItemValueSet(
+  const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+{
+  OutputNamedItemValueSet named_item_value_set;
+
+  for (auto& named_discrete_function : named_discrete_function_list) {
+    const IDiscreteFunction& i_discrete_function = *named_discrete_function->discreteFunction();
+
+    switch (i_discrete_function.descriptor().type()) {
+    case DiscreteFunctionType::P0: {
+      WriterBase::registerDiscreteFunctionP0(*named_discrete_function, named_item_value_set);
+      break;
+    }
+    default: {
+      std::ostringstream error_msg;
+      error_msg << "the type of discrete function of " << rang::fgB::blue << named_discrete_function->name()
+                << rang::style::reset << " is not supported";
+      throw NormalError(error_msg.str());
+    }
+    }
+  }
+
+  return named_item_value_set;
+}
+
+double
+WriterBase::getLastTime() const
+{
+  if (m_saved_times.size() > 0) {
+    return m_saved_times[m_saved_times.size() - 1];
+  } else {
+    return -std::numeric_limits<double>::max();
+  }
+}
+
+WriterBase::WriterBase(const std::string& base_filename, const double& time_period)
+  : m_base_filename{base_filename}, m_time_period{time_period}
+{}
+
+void
+WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                          double time) const
+{
+  const double last_time = getLastTime();
+  if (time == last_time)
+    return;   // output already performed
+
+  if (time >= last_time + m_time_period) {
+    this->write(named_discrete_function_list, time);
+    m_saved_times.push_back(time);
+  }
+}
+
+void
+WriterBase::writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                        double time) const
+{
+  if (time == getLastTime())
+    return;   // output already performed
+
+  this->write(named_discrete_function_list, time);
+  m_saved_times.push_back(time);
+}
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e9a08ad123361178327ddb984e5a3da5bdfa909
--- /dev/null
+++ b/src/output/WriterBase.hpp
@@ -0,0 +1,54 @@
+#ifndef WRITER_BASE_HPP
+#define WRITER_BASE_HPP
+
+#include <output/IWriter.hpp>
+
+#include <string>
+
+class IMesh;
+class OutputNamedItemValueSet;
+
+class WriterBase : public IWriter
+{
+ protected:
+  const std::string m_base_filename;
+  const double m_time_period;
+
+  mutable std::vector<double> m_saved_times;
+
+ private:
+  template <size_t Dimension, typename DataType>
+  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+
+  template <size_t Dimension>
+  static void registerDiscreteFunctionP0(const std::string& name, const IDiscreteFunction&, OutputNamedItemValueSet&);
+
+  static void registerDiscreteFunctionP0(const NamedDiscreteFunction&, OutputNamedItemValueSet&);
+
+ protected:
+  std::shared_ptr<const IMesh> _getMesh(
+    const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+
+  OutputNamedItemValueSet _getOutputNamedItemValueSet(
+    const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+
+  virtual void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                     double time) const = 0;
+
+ public:
+  double getLastTime() const;
+
+  void writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                     double time) const final;
+
+  void writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+                   double time) const final;
+
+  WriterBase() = delete;
+
+  WriterBase(const std::string& base_filename, const double& time_period);
+
+  virtual ~WriterBase() = default;
+};
+
+#endif   // WRITER_BASE_HPP