From 6b2d417ae691374ae1121ac1c5e156f154cda746 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 9 Mar 2021 22:00:13 +0100
Subject: [PATCH] Add `gnuplot_1d_writer`

The aim of this variant is to provide a more classic version of the
output: one value per cell or per node.

It can only be used for 1d data.

It remains to make it work in parallel.
---
 src/language/modules/WriterModule.cpp |  11 +
 src/output/CMakeLists.txt             |   1 +
 src/output/GnuplotWriter.cpp          |  16 +-
 src/output/GnuplotWriter1D.cpp        | 312 ++++++++++++++++++++++++++
 src/output/GnuplotWriter1D.hpp        |  75 +++++++
 5 files changed, 405 insertions(+), 10 deletions(-)
 create mode 100644 src/output/GnuplotWriter1D.cpp
 create mode 100644 src/output/GnuplotWriter1D.hpp

diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index 611952939..9f2fd98f1 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -6,6 +6,7 @@
 #include <mesh/GmshReader.hpp>
 #include <mesh/Mesh.hpp>
 #include <output/GnuplotWriter.hpp>
+#include <output/GnuplotWriter1D.hpp>
 #include <output/IWriter.hpp>
 #include <output/NamedDiscreteFunction.hpp>
 #include <output/VTKWriter.hpp>
@@ -39,6 +40,16 @@ WriterModule::WriterModule()
 
                               ));
 
+  this->_addBuiltinFunction("gnuplot_1d_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<GnuplotWriter1D>(filename, period);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("name_output",
                             std::make_shared<BuiltinFunctionEmbedder<
                               std::shared_ptr<const NamedDiscreteFunction>(std::shared_ptr<const IDiscreteFunction>,
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index a3e679965..1f3380903 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -3,6 +3,7 @@
 add_library(
   PugsOutput
   GnuplotWriter.cpp
+  GnuplotWriter1D.cpp
   VTKWriter.cpp
   WriterBase.cpp)
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index f503c996c..4b1bf84b8 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -177,7 +177,7 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
                       const OutputNamedItemValueSet& output_named_item_value_set,
                       double time) const
 {
-  if (parallel::rank() == 0) {
+  if (parallel::size() == 1) {
     std::ofstream fout{_getFilename()};
     fout.precision(15);
     fout.setf(std::ios_base::scientific);
@@ -186,17 +186,13 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
     fout << "\n# time = " << time << "\n\n";
 
     this->_writePreamble<MeshType::Dimension>(output_named_item_value_set, fout);
-  }
 
-  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);
+    fout.precision(15);
+    fout.setf(std::ios_base::scientific);
 
-      this->_writeValues(mesh, output_named_item_value_set, fout);
-    }
-    parallel::barrier();
+    this->_writeValues(mesh, output_named_item_value_set, fout);
+  } else {
+    throw NotImplementedError("not implemented in parallel");
   }
 }
 
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
new file mode 100644
index 000000000..bafef0856
--- /dev/null
+++ b/src/output/GnuplotWriter1D.cpp
@@ -0,0 +1,312 @@
+#include <output/GnuplotWriter1D.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
+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;
+  sout << '.' << std::setfill('0') << std::setw(4) << m_saved_times.size() << ".gnu";
+  return sout.str();
+}
+
+template <typename DataType>
+bool
+GnuplotWriter1D::_is_cell_value(const CellValue<const DataType>&) const
+{
+  return true;
+}
+
+template <typename DataType>
+bool
+GnuplotWriter1D::_is_cell_value(const NodeValue<const DataType>&) const
+{
+  return false;
+}
+
+template <typename DataType>
+bool
+GnuplotWriter1D::_is_node_value(const CellValue<const DataType>&) const
+{
+  return false;
+}
+
+template <typename DataType>
+bool
+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;
+  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
+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
+{
+  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
+GnuplotWriter1D::_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
+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];
+    }
+  } 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
+GnuplotWriter1D::_writeNodeValue(const CellValue<ValueType>&, NodeId, std::ostream&) const
+{}
+
+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';
+    }
+  } 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
+GnuplotWriter1D::_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
+GnuplotWriter1D::writeMesh(const std::shared_ptr<const IMesh>&) const
+{
+  throw UnexpectedError("This function should not be called");
+}
+
+void
+GnuplotWriter1D::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: {
+    std::ostringstream errorMsg;
+    errorMsg << "gnuplot_1d_writer is not available in dimension " << std::to_string(mesh->dimension()) << '\n'
+             << rang::style::bold << "note:" << rang::style::reset << " one can use " << rang::fgB::blue
+             << "gnuplot_writer" << rang::style::reset << " in dimension 2";
+    throw NormalError(errorMsg.str());
+  }
+  default: {
+    throw NormalError("gnuplot format is not available in dimension " + std::to_string(mesh->dimension()));
+  }
+  }
+}
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
new file mode 100644
index 000000000..a299b14ee
--- /dev/null
+++ b/src/output/GnuplotWriter1D.hpp
@@ -0,0 +1,75 @@
+#ifndef GNUPLOT_WRITER_1D_HPP
+#define GNUPLOT_WRITER_1D_HPP
+
+#include <output/WriterBase.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <output/OutputNamedItemValueSet.hpp>
+
+class IMesh;
+
+#include <string>
+
+class GnuplotWriter1D : 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;
+
+  GnuplotWriter1D(const std::string& base_filename, const double time_period) : WriterBase(base_filename, time_period)
+  {}
+
+  ~GnuplotWriter1D() = default;
+};
+
+#endif   // GNUPLOT_WRITER_1D_HPP
-- 
GitLab