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/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 38dc08b6e4bec82be8cab2666080563bd7cb9736..b31fe47d0eba87ff2e7f0823a530bf170e410d6c 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -212,6 +212,76 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
));
+ scheme_module._addBuiltinFunction("tensorProduct", std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<1> b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<2> b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> a,
+ const TinyVector<3>& b) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<1> a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<2> a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
+ scheme_module._addBuiltinFunction("tensorProduct",
+ std::function(
+
+ [](const TinyVector<3>& a, std::shared_ptr<const DiscreteFunctionVariant> b)
+ -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return tensorProduct(a, b);
+ }
+
+ ));
+
scheme_module._addBuiltinFunction("det", std::function(
[](std::shared_ptr<const DiscreteFunctionVariant> A)
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index afb9dbe0cf16f100e1367c9843325aa5bdec80ad..5f5199af2d850edc93dc961e00dbe4a4ead58116 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -71,6 +71,21 @@ MathModule::MathModule()
return dot(x, y);
}));
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<1> x, const TinyVector<1> y) -> TinyMatrix<1> {
+ return tensorProduct(x, y);
+ }));
+
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<2> x, const TinyVector<2> y) -> TinyMatrix<2> {
+ return tensorProduct(x, y);
+ }));
+
+ this->_addBuiltinFunction("tensorProduct",
+ std::function([](const TinyVector<3>& x, const TinyVector<3>& y) -> TinyMatrix<3> {
+ return tensorProduct(x, y);
+ }));
+
this->_addBuiltinFunction("det", std::function([](const TinyMatrix<1> A) -> double { return det(A); }));
this->_addBuiltinFunction("det", std::function([](const TinyMatrix<2>& A) -> double { return det(A); }));
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index f0c1b03df96cf1db239bb300cae2fbf797a17fd5..8481a1ed20958162866820c139983971e35f373b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -665,6 +665,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 837bc1252d64c5c851ac65561c049b79403ce980..6d84506501de277f1970c13a7c24cd587ed79f20 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -12,6 +12,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>
@@ -81,6 +82,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/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
index 715de1be2019070e517cb30639c199f721258cd2..961cc0e4092ba57c001866a52d7aed113e6b0f5d 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.cpp
@@ -315,6 +315,110 @@ template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<2>&
template std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<3>&,
const std::shared_ptr<const DiscreteFunctionVariant>&);
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& g_v)
+{
+ if (not hasSameMesh({f_v, g_v})) {
+ throw NormalError("operands are defined on different meshes");
+ }
+
+ return std::visit(
+ [&](auto&& f, auto&& g) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using TypeOfF = std::decay_t<decltype(f)>;
+ using TypeOfG = std::decay_t<decltype(g)>;
+ if constexpr (not std::is_same_v<TypeOfF, TypeOfG>) {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, g));
+ } else {
+ using DataType = std::decay_t<typename TypeOfF::data_type>;
+ if constexpr (is_discrete_function_P0_v<TypeOfF>) {
+ if constexpr (is_tiny_vector_v<DataType>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct(f, g));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ }
+ },
+ f_v->discreteFunction(), g_v->discreteFunction());
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>& f, const TinyVector<VectorDimension>& a)
+{
+ return std::visit(
+ [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
+ using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
+ if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(discrete_function0, a));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(f, a));
+ }
+ },
+ f->discreteFunction());
+}
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+tensorProduct(const TinyVector<VectorDimension>& a, const std::shared_ptr<const DiscreteFunctionVariant>& f)
+{
+ return std::visit(
+ [&](auto&& discrete_function0) -> std::shared_ptr<DiscreteFunctionVariant> {
+ using DiscreteFunction0Type = std::decay_t<decltype(discrete_function0)>;
+ if constexpr (is_discrete_function_P0_v<DiscreteFunction0Type>) {
+ using DataType = std::decay_t<typename DiscreteFunction0Type::data_type>;
+ if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (std::is_same_v<DataType, TinyVector<VectorDimension>>) {
+ return std::make_shared<DiscreteFunctionVariant>(tensorProduct<DataType>(a, discrete_function0));
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::invalidOperandType(f));
+ }
+ } else {
+ throw NormalError(EmbeddedDiscreteFunctionUtils::incompatibleOperandTypes(a, f));
+ }
+ },
+ f->discreteFunction());
+}
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<1>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<2>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<3>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<1>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<2>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(
+ const TinyVector<3>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
std::shared_ptr<const DiscreteFunctionVariant>
det(const std::shared_ptr<const DiscreteFunctionVariant>& A)
{
diff --git a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
index 66b784921da372cc7dfa2d616b94d6fbfe403ee3..bcda26113a4fac6d72e7b54edfe0f0e8d8d74b2d 100644
--- a/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedDiscreteFunctionMathFunctions.hpp
@@ -69,6 +69,17 @@ template <size_t VectorDimension>
std::shared_ptr<const DiscreteFunctionVariant> dot(const TinyVector<VectorDimension>&,
const std::shared_ptr<const DiscreteFunctionVariant>&);
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const std::shared_ptr<const DiscreteFunctionVariant>&,
+ const TinyVector<VectorDimension>&);
+
+template <size_t VectorDimension>
+std::shared_ptr<const DiscreteFunctionVariant> tensorProduct(const TinyVector<VectorDimension>&,
+ const std::shared_ptr<const DiscreteFunctionVariant>&);
+
std::shared_ptr<const DiscreteFunctionVariant> det(const std::shared_ptr<const DiscreteFunctionVariant>&);
std::shared_ptr<const DiscreteFunctionVariant> trace(const std::shared_ptr<const DiscreteFunctionVariant>&);
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 9e73b449adc517ba80379b413a606f4ef71b807a..4b249aa2e5f8bd11a6251f9edac463d69b8562df 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 7bfc671601b5a442619354d3c989572fa0d7cc2b..f1a1f030315c2991c814c8af2212a6d34bc609ab 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,
@@ -67,9 +60,11 @@ class GnuplotWriter final : public WriterBase
return Type::gnuplot;
}
- 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 f5d6dd7b2c0cd35696b21457194f72c3a4289c9a..10aafe1f5fbd478a8aa0d39b86ecd5f76b8aa64f 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
@@ -341,7 +236,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) {
@@ -357,7 +252,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 8b49e1935f5526603960e6995613cb3481eabd30..d105fdf05e52ef81bf37e76ecfa61eac8ed38edf 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,
@@ -66,9 +48,10 @@ class GnuplotWriter1D final : public WriterBase
return Type::gnuplot_1d;
}
- 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..d72b77162f65b2500d89313f3a0005ce9aaa8182
--- /dev/null
+++ b/src/output/GnuplotWriterRaw.hpp
@@ -0,0 +1,60 @@
+#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:
+ Type
+ type() const final
+ {
+ return Type::gnuplot_raw;
+ }
+
+ 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
diff --git a/src/output/IWriter.hpp b/src/output/IWriter.hpp
index 742d0cf349a6d64356cc6021ec0a380247c1d50b..85060e37e8abbc46c06fe70933c3a032721dd01c 100644
--- a/src/output/IWriter.hpp
+++ b/src/output/IWriter.hpp
@@ -15,6 +15,7 @@ class IWriter
{
gnuplot,
gnuplot_1d,
+ gnuplot_raw,
vtk
};
diff --git a/src/scheme/AxisBoundaryConditionDescriptor.hpp b/src/scheme/AxisBoundaryConditionDescriptor.hpp
index 07ab3be487f30af021fb51ceae84892091e3ae8b..1dbbf98997a380e05d2681e17962ab4181ce5d47 100644
--- a/src/scheme/AxisBoundaryConditionDescriptor.hpp
+++ b/src/scheme/AxisBoundaryConditionDescriptor.hpp
@@ -2,11 +2,11 @@
#define AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class AxisBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class AxisBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -16,23 +16,15 @@ class AxisBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return os;
}
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-
public:
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
return Type::axis;
}
- AxisBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
- : m_boundary_descriptor(boundary_descriptor)
+ AxisBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}
{
;
}
diff --git a/src/scheme/BoundaryConditionDescriptorBase.hpp b/src/scheme/BoundaryConditionDescriptorBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d05dfd6713e1f8026254567c9859dba8dc287785
--- /dev/null
+++ b/src/scheme/BoundaryConditionDescriptorBase.hpp
@@ -0,0 +1,34 @@
+#ifndef BOUNDARY_CONDITION_DESCRIPTOR_BASE_HPP
+#define BOUNDARY_CONDITION_DESCRIPTOR_BASE_HPP
+
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+class BoundaryConditionDescriptorBase : public IBoundaryConditionDescriptor
+{
+ protected:
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+ const std::shared_ptr<const IBoundaryDescriptor>&
+ boundaryDescriptor_shared() const final
+ {
+ return m_boundary_descriptor;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const final
+ {
+ return *m_boundary_descriptor;
+ }
+
+ BoundaryConditionDescriptorBase(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+ : m_boundary_descriptor{boundary_descriptor}
+ {}
+
+ BoundaryConditionDescriptorBase(const BoundaryConditionDescriptorBase&) = delete;
+ BoundaryConditionDescriptorBase(BoundaryConditionDescriptorBase&&) = delete;
+
+ virtual ~BoundaryConditionDescriptorBase() = default;
+};
+
+#endif // BOUNDARY_CONDITION_DESCRIPTOR_BASE_HPP
diff --git a/src/scheme/DirichletBoundaryConditionDescriptor.hpp b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
index 32e34081002989cb0d6164332b070c5ad712d1bb..2ddfa649cd151dec643984becae7ae5718ef000a 100644
--- a/src/scheme/DirichletBoundaryConditionDescriptor.hpp
+++ b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
@@ -3,11 +3,11 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class DirichletBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -19,7 +19,6 @@ class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
const std::string m_name;
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
const FunctionSymbolId m_rhs_symbol_id;
public:
@@ -35,12 +34,6 @@ class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return m_rhs_symbol_id;
}
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
@@ -48,9 +41,9 @@ class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
}
DirichletBoundaryConditionDescriptor(const std::string& name,
- std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor,
const FunctionSymbolId& rhs_symbol_id)
- : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+ : BoundaryConditionDescriptorBase{boundary_descriptor}, m_name{name}, m_rhs_symbol_id{rhs_symbol_id}
{
;
}
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index df43f52905837036d0ef843fbc5635f0d1274c6d..ea0ec56aeeebb28e4af381ae352a20c9abd43bd4 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -606,6 +606,51 @@ class DiscreteFunctionP0
return result;
}
+ template <typename LHSDataType>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(LHSDataType{}, DataType{}))>
+ tensorProduct(const LHSDataType& u, const DiscreteFunctionP0& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<LHSDataType>>);
+ Assert(v.m_cell_values.isBuilt());
+ DiscreteFunctionP0<decltype(tensorProduct(LHSDataType{}, DataType{}))> result{v.m_mesh_v};
+ parallel_for(
+ v.m_mesh_v->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u, v[cell_id]); });
+
+ return result;
+ }
+
+ template <typename RHSDataType>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(DataType{}, RHSDataType{}))>
+ tensorProduct(const DiscreteFunctionP0& u, const RHSDataType& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<RHSDataType>>);
+ Assert(u.m_cell_values.isBuilt());
+ DiscreteFunctionP0<decltype(tensorProduct(DataType{}, RHSDataType{}))> result{u.m_mesh_v};
+ parallel_for(
+ u.m_mesh_v->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u[cell_id], v); });
+
+ return result;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE friend DiscreteFunctionP0<decltype(tensorProduct(DataType{}, DataType2{}))>
+ tensorProduct(const DiscreteFunctionP0& u, const DiscreteFunctionP0<DataType2>& v)
+ {
+ static_assert(is_tiny_vector_v<std::decay_t<DataType>>);
+ static_assert(is_tiny_vector_v<std::decay_t<DataType2>>);
+ Assert(u.m_cell_values.isBuilt());
+ Assert(v.m_cell_values.isBuilt());
+ Assert(u.m_mesh_v->id() == v.m_mesh_v->id());
+ DiscreteFunctionP0<decltype(tensorProduct(DataType{}, DataType2{}))> result{u.m_mesh_v};
+ parallel_for(
+ u.m_mesh_v->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { result[cell_id] = tensorProduct(u[cell_id], v[cell_id]); });
+
+ return result;
+ }
+
PUGS_INLINE friend double
min(const DiscreteFunctionP0& f)
{
diff --git a/src/scheme/ExternalBoundaryConditionDescriptor.hpp b/src/scheme/ExternalBoundaryConditionDescriptor.hpp
index 05fa10eb13822de06bead5c17a674f6b7f94d24c..d015547f88862f7d16825f3888e1814a7702a7cd 100644
--- a/src/scheme/ExternalBoundaryConditionDescriptor.hpp
+++ b/src/scheme/ExternalBoundaryConditionDescriptor.hpp
@@ -2,12 +2,12 @@
#define EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <utils/Socket.hpp>
#include <memory>
-class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class ExternalBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -19,7 +19,6 @@ class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
const std::string m_name;
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
const std::shared_ptr<const Socket> m_socket;
public:
@@ -35,12 +34,6 @@ class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return m_socket;
}
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
@@ -48,9 +41,9 @@ class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
}
ExternalBoundaryConditionDescriptor(const std::string& name,
- std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor,
const std::shared_ptr<const Socket>& socket)
- : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_socket{socket}
+ : BoundaryConditionDescriptorBase{boundary_descriptor}, m_name{name}, m_socket{socket}
{}
ExternalBoundaryConditionDescriptor(const ExternalBoundaryConditionDescriptor&) = delete;
diff --git a/src/scheme/FixedBoundaryConditionDescriptor.hpp b/src/scheme/FixedBoundaryConditionDescriptor.hpp
index bc7ca9cca66f393a0458982fd0d169cff49b4a10..b8b5446248f040505d8a8f3653142076ae4d3d97 100644
--- a/src/scheme/FixedBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FixedBoundaryConditionDescriptor.hpp
@@ -3,11 +3,11 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class FixedBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class FixedBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -20,20 +20,14 @@ class FixedBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
public:
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
return Type::fixed;
}
- FixedBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
- : m_boundary_descriptor(boundary_descriptor)
+ FixedBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}
{}
FixedBoundaryConditionDescriptor(const FixedBoundaryConditionDescriptor&) = delete;
diff --git a/src/scheme/FourierBoundaryConditionDescriptor.hpp b/src/scheme/FourierBoundaryConditionDescriptor.hpp
index 327a98438e1dbab2e5b141f0856525fadb3eeeba..c04995634691554ce3174b7a6b23f0cf086ebed8 100644
--- a/src/scheme/FourierBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FourierBoundaryConditionDescriptor.hpp
@@ -3,11 +3,11 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class FourierBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -19,7 +19,6 @@ class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
const std::string m_name;
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
const FunctionSymbolId m_mass_symbol_id;
const FunctionSymbolId m_rhs_symbol_id;
@@ -42,12 +41,6 @@ class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return m_rhs_symbol_id;
}
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
@@ -55,11 +48,11 @@ class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
}
FourierBoundaryConditionDescriptor(const std::string& name,
- std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor,
const FunctionSymbolId& mass_symbol_id,
const FunctionSymbolId& rhs_symbol_id)
- : m_name{name},
- m_boundary_descriptor(boundary_descriptor),
+ : BoundaryConditionDescriptorBase{boundary_descriptor},
+ m_name{name},
m_mass_symbol_id{mass_symbol_id},
m_rhs_symbol_id{rhs_symbol_id}
{
diff --git a/src/scheme/FreeBoundaryConditionDescriptor.hpp b/src/scheme/FreeBoundaryConditionDescriptor.hpp
index 2d46f2a36a04cbaca0abc5b83c9ad1983e844841..8b6df26f7bbd155d610cf6bb5de0a47817a86830 100644
--- a/src/scheme/FreeBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FreeBoundaryConditionDescriptor.hpp
@@ -2,11 +2,11 @@
#define FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class FreeBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class FreeBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -16,15 +16,7 @@ class FreeBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return os;
}
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-
public:
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
@@ -32,7 +24,7 @@ class FreeBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
}
FreeBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
- : m_boundary_descriptor(boundary_descriptor)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}
{
;
}
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 1350aa4dba4f7f6739b7de6689c10554ceec6971..58e8b6076ff2c363049de6f107e6a949966d6c72 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -2,6 +2,7 @@
#define I_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <iostream>
+#include <memory>
class IBoundaryDescriptor;
@@ -32,7 +33,8 @@ class IBoundaryConditionDescriptor
return bcd._write(os);
}
- virtual const IBoundaryDescriptor& boundaryDescriptor() const = 0;
+ virtual const std::shared_ptr<const IBoundaryDescriptor>& boundaryDescriptor_shared() const = 0;
+ virtual const IBoundaryDescriptor& boundaryDescriptor() const = 0;
virtual Type type() const = 0;
diff --git a/src/scheme/InflowBoundaryConditionDescriptor.hpp b/src/scheme/InflowBoundaryConditionDescriptor.hpp
index e648c1593cb4316f0cc77a4f18946b0457d65b5e..0b460846983946468d4161672fbc1e9aeb3912be 100644
--- a/src/scheme/InflowBoundaryConditionDescriptor.hpp
+++ b/src/scheme/InflowBoundaryConditionDescriptor.hpp
@@ -3,11 +3,11 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class InflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class InflowBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -17,7 +17,6 @@ class InflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return os;
}
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
const FunctionSymbolId m_function_symbol_id;
public:
@@ -27,21 +26,15 @@ class InflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return m_function_symbol_id;
}
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
return Type::inflow;
}
- InflowBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ InflowBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor,
const FunctionSymbolId& function_symbol_id)
- : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
+ : BoundaryConditionDescriptorBase{boundary_descriptor}, m_function_symbol_id{function_symbol_id}
{
;
}
diff --git a/src/scheme/NeumannBoundaryConditionDescriptor.hpp b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
index d318d2842212e5b03816b2e63d6328127d15cf23..2d63df7451dba2a7a0cb222665757d1b6f692b09 100644
--- a/src/scheme/NeumannBoundaryConditionDescriptor.hpp
+++ b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
@@ -3,11 +3,11 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class NeumannBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -19,7 +19,6 @@ class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
const std::string m_name;
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
const FunctionSymbolId m_rhs_symbol_id;
public:
@@ -35,12 +34,6 @@ class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return m_rhs_symbol_id;
}
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
@@ -48,9 +41,9 @@ class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
}
NeumannBoundaryConditionDescriptor(const std::string& name,
- std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor,
const FunctionSymbolId& rhs_symbol_id)
- : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+ : BoundaryConditionDescriptorBase{boundary_descriptor}, m_name{name}, m_rhs_symbol_id{rhs_symbol_id}
{
;
}
diff --git a/src/scheme/OutflowBoundaryConditionDescriptor.hpp b/src/scheme/OutflowBoundaryConditionDescriptor.hpp
index 83ef7a775a81269a5d3930ed2e07d5634d99959a..bc999c501f1e4d539736e06a6ac59c76e39a170e 100644
--- a/src/scheme/OutflowBoundaryConditionDescriptor.hpp
+++ b/src/scheme/OutflowBoundaryConditionDescriptor.hpp
@@ -2,11 +2,11 @@
#define OUTFLOW_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class OutflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class OutflowBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -16,23 +16,15 @@ class OutflowBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return os;
}
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-
public:
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
return Type::outflow;
}
- OutflowBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
- : m_boundary_descriptor(boundary_descriptor)
+ OutflowBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}
{
;
}
diff --git a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
index 9364090dde18490a4803662c7eb770d9f6354b1c..1172bf7e40c5b2565dd82f216617298addda1e34 100644
--- a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
+++ b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
@@ -2,11 +2,11 @@
#define SYMMETRY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
#include <mesh/IBoundaryDescriptor.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
#include <memory>
-class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+class SymmetryBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
{
private:
std::ostream&
@@ -16,23 +16,15 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
return os;
}
- std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-
public:
- const IBoundaryDescriptor&
- boundaryDescriptor() const final
- {
- return *m_boundary_descriptor;
- }
-
Type
type() const final
{
return Type::symmetry;
}
- SymmetryBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
- : m_boundary_descriptor(boundary_descriptor)
+ SymmetryBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor)
+ : BoundaryConditionDescriptorBase{boundary_descriptor}
{
;
}
diff --git a/src/utils/Messenger.cpp b/src/utils/Messenger.cpp
index 05dfa43e5a22f5223fd6c1e8291302c5366f5015..14891477a2b804c4f0c8698fdfa51a8ea2470c8b 100644
--- a/src/utils/Messenger.cpp
+++ b/src/utils/Messenger.cpp
@@ -9,10 +9,10 @@ namespace parallel
Messenger* Messenger::m_instance = nullptr;
void
-Messenger::create(int& argc, char* argv[])
+Messenger::create(int& argc, char* argv[], bool parallel_output)
{
if (Messenger::m_instance == nullptr) {
- Messenger::m_instance = new Messenger(argc, argv);
+ Messenger::m_instance = new Messenger(argc, argv, parallel_output);
} else {
throw UnexpectedError("Messenger already created");
}
@@ -28,7 +28,7 @@ Messenger::destroy()
}
}
-Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
+Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[], bool parallel_output)
{
#ifdef PUGS_HAS_MPI
MPI_Init(&argc, &argv);
@@ -66,7 +66,7 @@ Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
return size;
}();
- if (m_rank != 0) {
+ if ((not parallel_output) and (m_rank != 0)) {
// LCOV_EXCL_START
std::cout.setstate(std::ios::badbit);
std::cerr.setstate(std::ios::badbit);
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index ac79ee8490882a1bac8a9fd74db58b6cd8dccdb7..a04c664d7dcbe36817d2748a6f27e0b8d8dc3481 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -85,7 +85,7 @@ class Messenger
};
static Messenger* m_instance;
- Messenger(int& argc, char* argv[]);
+ Messenger(int& argc, char* argv[], bool parallel_output);
#ifdef PUGS_HAS_MPI
MPI_Comm m_pugs_comm_world = MPI_COMM_WORLD;
@@ -406,7 +406,7 @@ class Messenger
}
public:
- static void create(int& argc, char* argv[]);
+ static void create(int& argc, char* argv[], bool parallel_output = false);
static void destroy();
PUGS_INLINE
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index f9b6c94985d4e273aa409c7e6aa50e27ef58d01d..b3d7d6c9a941f3b6e1c03e6d9cc0fb1bec24a1a7 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -86,9 +86,10 @@ setDefaultOMPEnvironment()
std::string
initialize(int& argc, char* argv[])
{
- bool enable_fpe = true;
- bool enable_signals = true;
- int nb_threads = -1;
+ bool enable_fpe = true;
+ bool enable_signals = true;
+ int nb_threads = -1;
+ bool parallel_output = false;
ParallelChecker::Mode pc_mode = ParallelChecker::Mode::automatic;
std::string pc_filename = ParallelChecker::instance().filename();
@@ -133,9 +134,11 @@ initialize(int& argc, char* argv[])
app.add_flag("-p,--pause-on-error", pause_on_error, "Pause for debugging on unexpected error [default: false]");
bool reproducible_sums = true;
- app.add_flag("--reproducible-sums,!--no-reproducible-sums", show_preamble,
+ app.add_flag("--reproducible-sums,!--no-reproducible-sums", reproducible_sums,
"Special treatment of array sums to ensure reproducibility [default: true]");
+ app.add_flag("--parallel-output", parallel_output, "All MPI processes output to console [default: false]");
+
std::map<std::string, ParallelChecker::Mode> pc_mode_map{{"auto", ParallelChecker::Mode::automatic},
{"write", ParallelChecker::Mode::write},
{"read", ParallelChecker::Mode::read}};
@@ -176,7 +179,8 @@ initialize(int& argc, char* argv[])
SignalManager::setPauseForDebug(pause_on_error);
ReproducibleSumManager::setReproducibleSums(reproducible_sums);
}
- parallel::Messenger::create(argc, argv);
+
+ parallel::Messenger::create(argc, argv, parallel_output);
PETScWrapper::initialize(argc, argv);
SLEPcWrapper::initialize(argc, argv);
diff --git a/src/utils/checkpointing/IWriterHFType.hpp b/src/utils/checkpointing/IWriterHFType.hpp
index 124cb2ee53f6bf477361cd381e98af366e626231..046998283e0c4a8567e1624741ddd920108413e7 100644
--- a/src/utils/checkpointing/IWriterHFType.hpp
+++ b/src/utils/checkpointing/IWriterHFType.hpp
@@ -8,7 +8,10 @@
HighFive::EnumType<IWriter::Type> PUGS_INLINE
create_enum_i_writer_type()
{
- return {{"gnuplot", IWriter::Type::gnuplot}, {"gnuplot_1d", IWriter::Type::gnuplot_1d}, {"vtk", IWriter::Type::vtk}};
+ return {{"gnuplot", IWriter::Type::gnuplot},
+ {"gnuplot_1d", IWriter::Type::gnuplot_1d},
+ {"gnuplot_raw", IWriter::Type::gnuplot_raw},
+ {"vtk", IWriter::Type::vtk}};
}
HIGHFIVE_REGISTER_TYPE(IWriter::Type, create_enum_i_writer_type)
diff --git a/src/utils/checkpointing/ReadIWriter.cpp b/src/utils/checkpointing/ReadIWriter.cpp
index 1f4dc1a55c4cefc0eaa9b8de05209b764448b3d0..d7037808859bd8916530e98ce1705214ef1c5972 100644
--- a/src/utils/checkpointing/ReadIWriter.cpp
+++ b/src/utils/checkpointing/ReadIWriter.cpp
@@ -4,6 +4,7 @@
#include <language/utils/EmbeddedData.hpp>
#include <output/GnuplotWriter.hpp>
#include <output/GnuplotWriter1D.hpp>
+#include <output/GnuplotWriterRaw.hpp>
#include <output/VTKWriter.hpp>
#include <utils/checkpointing/IWriterHFType.hpp>
@@ -30,6 +31,10 @@ readIWriter(const std::string& symbol_name, const HighFive::Group& symbol_table_
p_writer = std::make_shared<GnuplotWriter1D>(base_filename);
break;
}
+ case IWriter::Type::gnuplot_raw: {
+ p_writer = std::make_shared<GnuplotWriterRaw>(base_filename);
+ break;
+ }
case IWriter::Type::vtk: {
p_writer = std::make_shared<VTKWriter>(base_filename);
break;
diff --git a/src/utils/checkpointing/WriteIWriter.cpp b/src/utils/checkpointing/WriteIWriter.cpp
index f05c8e77c424856561ee9d579392390bd384650b..b04b6b6cd3c2cff62762285c1fbd64e466be1182 100644
--- a/src/utils/checkpointing/WriteIWriter.cpp
+++ b/src/utils/checkpointing/WriteIWriter.cpp
@@ -29,6 +29,7 @@ writeIWriter(const std::string& symbol_name,
switch (iwriter_p->type()) {
case IWriter::Type::gnuplot:
case IWriter::Type::gnuplot_1d:
+ case IWriter::Type::gnuplot_raw:
case IWriter::Type::vtk: {
const WriterBase& writer = dynamic_cast<const WriterBase&>(*iwriter_p);
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index 36b444e58f014f39f8d752b4be129a568a661fa8..014c7b123e182402002f97e44c8d3115274367d8 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -362,6 +362,40 @@ let s:R, s = dot(x,y);
CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
}
+ { // dot
+ tested_function_set.insert("tensorProduct:R^1*R^1");
+ std::string_view data = R"(
+import math;
+let x:R^1, x = -2;
+let y:R^1, y = 4;
+let s:R^1x1, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", TinyMatrix<1>{-2 * 4});
+ }
+
+ { // dot
+ tested_function_set.insert("tensorProduct:R^2*R^2");
+ std::string_view data = R"(
+import math;
+let x:R^2, x = [-2, 3];
+let y:R^2, y = [4, 3];
+let s:R^2x2, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", tensorProduct(TinyVector<2>{-2, 3}, TinyVector<2>{4, 3}));
+ }
+
+ { // dot
+ tested_function_set.insert("tensorProduct:R^3*R^3");
+ std::string_view data = R"(
+import math;
+let x:R^3, x = [-2, 3, 4];
+let y:R^3, y = [4, 3, 5];
+let s:R^3x3, s = tensorProduct(x,y);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s",
+ tensorProduct(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
+ }
+
{ // det
tested_function_set.insert("det:R^1x1");
std::string_view data = R"(
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 5b1133b4e6c60290b06aa90e5c68b64a1c1b4973..14c96b0db712daa220951f407e66eb1f774f126e 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -2729,6 +2729,53 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
}
+ SECTION("tensorProduct(uh,hv)")
+ {
+ DiscreteFunctionP0<TinyVector<2>> uh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ uh[cell_id] = TinyVector<2>{x + 1, 2 * x - 3};
+ });
+
+ DiscreteFunctionP0<TinyVector<2>> vh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ vh[cell_id] = TinyVector<2>{2.3 * x, 1 - x};
+ });
+
+ CHECK_STD_BINARY_MATH_FUNCTION(uh, vh, tensorProduct);
+ }
+
+ SECTION("tensorProduct(uh,v)")
+ {
+ DiscreteFunctionP0<TinyVector<2>> uh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ uh[cell_id] = TinyVector<2>{x + 1, 2 * x - 3};
+ });
+
+ const TinyVector<2> v{1, 2};
+
+ CHECK_STD_BINARY_MATH_FUNCTION_WITH_RHS_VALUE(uh, v, tensorProduct);
+ }
+
+ SECTION("tensorProduct(u,hv)")
+ {
+ const TinyVector<2> u{3, -2};
+
+ DiscreteFunctionP0<TinyVector<2>> vh{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ vh[cell_id] = TinyVector<2>{2.3 * x, 1 - x};
+ });
+
+ CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, tensorProduct);
+ }
+
SECTION("scalar sum")
{
const CellValue<const double> cell_value = positive_function.cellValues();
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
index 5843a8ae766549e0828551bae8d33ff8c595935d..41f3c0d2233b1884051f0eafd5df637ce37ed7e0 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions1D.cpp
@@ -502,6 +502,29 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -675,6 +698,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions1D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
index 11582b9820d9e2f6d910773b2084517a47069b37..919d4f8280e764d263129f70a849064867bb6f2a 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions2D.cpp
@@ -504,6 +504,30 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions2D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -677,6 +701,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions2D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp b/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
index edca0a4eb9c9e015a7e0ed31fd74ca068d4258c5..99994e9d0a902add718dee6149c3c040765a7106 100644
--- a/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
+++ b/tests/test_EmbeddedDiscreteFunctionMathFunctions3D.cpp
@@ -502,6 +502,30 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions3D", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("tensorProduct Vh*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R1_u, p_R1_v, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R2_u, p_R2_v, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VH2_TO_VH_FUNCTION_EVALUATION(p_R3_u, p_R3_v, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_other_mesh_R1_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, p_other_mesh_R2_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, p_other_mesh_R3_u),
+ "error: operands are defined on different meshes");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, p_R3_u),
+ "error: incompatible operand types Vh(P0:R^1) and Vh(P0:R^3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1x1_u, p_R1x1_u), "error: invalid operand type Vh(P0:R^1x1)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2x2_u, p_R2x2_u), "error: invalid operand type Vh(P0:R^2x2)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3x3_u, p_R3x3_u), "error: invalid operand type Vh(P0:R^3x3)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, p_Vector2_w), "error: invalid operand type Vh(P0Vector:R)");
+ }
+
SECTION("det Vh -> Vh")
{
CHECK_EMBEDDED_VH_TO_VH_FUNCTION_EVALUATION(p_R1x1_u, det, //
@@ -675,6 +699,48 @@ TEST_CASE("EmbeddedDiscreteFunctionVariantMathFunctions3D", "[scheme]")
REQUIRE_THROWS_WITH(dot((TinyVector<1>{-1}), p_R3_u), "error: incompatible operand types R^1 and Vh(P0:R^3)");
}
+ SECTION("tensorProduct Vh*Rd -> Vh")
+ {
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R1_u, (TinyVector<1>{3}), tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R2_u, (TinyVector<2>{-6, 2}), tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_VHxW_TO_VH_FUNCTION_EVALUATION(p_R3_u, (TinyVector<3>{-1, 5, 2}), tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_u, TinyVector<2>{1, 2}), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(p_Vector3_u, TinyVector<2>{1, 2}),
+ "error: incompatible operand types Vh(P0Vector:R) and R^2");
+
+ REQUIRE_THROWS_WITH(tensorProduct(p_R1_u, (TinyVector<2>{-6, 2})),
+ "error: incompatible operand types Vh(P0:R^1) and R^2");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R2_u, (TinyVector<3>{-1, 5, 2})),
+ "error: incompatible operand types Vh(P0:R^2) and R^3");
+ REQUIRE_THROWS_WITH(tensorProduct(p_R3_u, (TinyVector<1>{-1})),
+ "error: incompatible operand types Vh(P0:R^3) and R^1");
+ }
+
+ SECTION("tensorProduct Rd*Vh -> Vh")
+ {
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<1>{3}), p_R1_u, tensorProduct, //
+ DiscreteFunctionR1, DiscreteFunctionR1x1);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<2>{-6, 2}), p_R2_u, tensorProduct, //
+ DiscreteFunctionR2, DiscreteFunctionR2x2);
+ CHECK_EMBEDDED_WxVH_TO_VH_FUNCTION_EVALUATION((TinyVector<3>{-1, 5, 2}), p_R3_u, tensorProduct, //
+ DiscreteFunctionR3, DiscreteFunctionR3x3);
+
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(tensorProduct(TinyVector<2>{1, 2}, p_Vector3_u),
+ "error: incompatible operand types R^2 and Vh(P0Vector:R)");
+
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<2>{-6, 2}), p_R1_u),
+ "error: incompatible operand types R^2 and Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<3>{-1, 5, 2}), p_R2_u),
+ "error: incompatible operand types R^3 and Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(tensorProduct((TinyVector<1>{-1}), p_R3_u),
+ "error: incompatible operand types R^1 and Vh(P0:R^3)");
+ }
+
SECTION("sum_of_R* Vh -> R*")
{
REQUIRE(sum_of<double>(p_u) == sum(p_u->get<DiscreteFunctionR>().cellValues()));
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index f299d1e3678217d4fc6bcb4c08e417264b22a252..89a480e1ac4993b3cdbbd6702f5e2130a1872ec0 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -13,7 +13,7 @@ TEST_CASE("MathModule", "[language]")
MathModule math_module;
const auto& name_builtin_function = math_module.getNameBuiltinFunctionMap();
- REQUIRE(name_builtin_function.size() == 42);
+ REQUIRE(name_builtin_function.size() == 45);
SECTION("Z -> N")
{
@@ -458,6 +458,63 @@ TEST_CASE("MathModule", "[language]")
}
}
+ SECTION("(R^d, R^d) -> R^dxd")
+ {
+ SECTION("tensorProduct:R^1*R^1")
+ {
+ TinyVector<1> arg0{3};
+ TinyVector<1> arg1{2};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^1*R^1");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<1> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+ }
+
+ SECTION("tensorProduct:R^2*R^2")
+ {
+ TinyVector<2> arg0{3, 2};
+ TinyVector<2> arg1{-2, 5};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^2*R^2");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<2> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+ }
+
+ SECTION("tensorProduct:R^3*R^3")
+ {
+ TinyVector<3> arg0{3, 2, 4};
+ TinyVector<3> arg1{-2, 5, 2};
+
+ DataVariant arg0_variant = arg0;
+ DataVariant arg1_variant = arg1;
+
+ auto i_function = name_builtin_function.find("tensorProduct:R^3*R^3");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg0_variant, arg1_variant});
+
+ const TinyMatrix<3> result = tensorProduct(arg0, arg1);
+ REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+ }
+ }
+
SECTION("R^dxd -> double")
{
SECTION("det:R^1x1")
diff --git a/tests/test_checkpointing_IWriter.cpp b/tests/test_checkpointing_IWriter.cpp
index 33d00bd90660280159e90a3ae73b92cd76f04d95..f14156544e9896c3a59d0b8691fee8f78f475253 100644
--- a/tests/test_checkpointing_IWriter.cpp
+++ b/tests/test_checkpointing_IWriter.cpp
@@ -7,6 +7,7 @@
#include <language/utils/EmbeddedData.hpp>
#include <output/GnuplotWriter.hpp>
#include <output/GnuplotWriter1D.hpp>
+#include <output/GnuplotWriterRaw.hpp>
#include <output/NamedDiscreteFunction.hpp>
#include <output/VTKWriter.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
@@ -51,12 +52,18 @@ TEST_CASE("checkpointing_IWriter", "[utils/checkpointing]")
EmbeddedData{std::make_shared<DataHandler<const IWriter>>(p_gnuplot_writer)}, file,
useless_group, symbol_table_group);
- auto p_gnuplot_writer_1d = std::make_shared<const GnuplotWriter1D>("gnuplot_1d_basename", 1.12);
+ auto p_gnuplot_writer_1d = std::make_shared<const GnuplotWriter1D>("gnuplot_1d_basename", 1.21);
p_gnuplot_writer_1d->periodManager().value().setSaveTime(2);
checkpointing::writeIWriter("gnuplot_1d",
EmbeddedData{std::make_shared<DataHandler<const IWriter>>(p_gnuplot_writer_1d)}, file,
useless_group, symbol_table_group);
+ auto p_gnuplot_writer_raw = std::make_shared<const GnuplotWriterRaw>("gnuplot_raw_basename", 1.18);
+ p_gnuplot_writer_raw->periodManager().value().setSaveTime(2.3);
+ checkpointing::writeIWriter("gnuplot_raw",
+ EmbeddedData{std::make_shared<DataHandler<const IWriter>>(p_gnuplot_writer_raw)},
+ file, useless_group, symbol_table_group);
+
const std::string vtk_filename = path / "vtk_example";
auto mesh_v = MeshDataBaseForTests::get().cartesian1DMesh();
@@ -77,9 +84,10 @@ TEST_CASE("checkpointing_IWriter", "[utils/checkpointing]")
file.flush();
- EmbeddedData read_gnuplot_writer = checkpointing::readIWriter("gnuplot", symbol_table_group);
- EmbeddedData read_gnuplot_writer_1d = checkpointing::readIWriter("gnuplot_1d", symbol_table_group);
- EmbeddedData read_vtk_writer = checkpointing::readIWriter("vtk", symbol_table_group);
+ EmbeddedData read_gnuplot_writer = checkpointing::readIWriter("gnuplot", symbol_table_group);
+ EmbeddedData read_gnuplot_writer_1d = checkpointing::readIWriter("gnuplot_1d", symbol_table_group);
+ EmbeddedData read_gnuplot_writer_raw = checkpointing::readIWriter("gnuplot_raw", symbol_table_group);
+ EmbeddedData read_vtk_writer = checkpointing::readIWriter("vtk", symbol_table_group);
auto get_value = [](const EmbeddedData& embedded_data) -> const IWriter& {
return *dynamic_cast<const DataHandler<const IWriter>&>(embedded_data.get()).data_ptr();
@@ -100,10 +108,21 @@ TEST_CASE("checkpointing_IWriter", "[utils/checkpointing]")
REQUIRE(read_gp_writer_1d.type() == IWriter::Type::gnuplot_1d);
REQUIRE(read_gp_writer_1d.baseFilename() == "gnuplot_1d_basename");
REQUIRE(read_gp_writer_1d.periodManager().has_value());
- REQUIRE(read_gp_writer_1d.periodManager().value().timePeriod() == 1.12);
- REQUIRE(read_gp_writer_1d.periodManager().value().nextTime() == 3.12);
+ REQUIRE(read_gp_writer_1d.periodManager().value().timePeriod() == 1.21);
+ REQUIRE(read_gp_writer_1d.periodManager().value().nextTime() == Catch::Approx(3.21));
REQUIRE(not read_gp_writer_1d.signature().has_value());
+ REQUIRE_NOTHROW(get_value(read_gnuplot_writer_raw));
+ REQUIRE_NOTHROW(dynamic_cast<const GnuplotWriterRaw&>(get_value(read_gnuplot_writer_raw)));
+ const GnuplotWriterRaw& read_gp_writer_raw =
+ dynamic_cast<const GnuplotWriterRaw&>(get_value(read_gnuplot_writer_raw));
+ REQUIRE(read_gp_writer_raw.type() == IWriter::Type::gnuplot_raw);
+ REQUIRE(read_gp_writer_raw.baseFilename() == "gnuplot_raw_basename");
+ REQUIRE(read_gp_writer_raw.periodManager().has_value());
+ REQUIRE(read_gp_writer_raw.periodManager().value().timePeriod() == 1.18);
+ REQUIRE(read_gp_writer_raw.periodManager().value().nextTime() == Catch::Approx(3.48));
+ REQUIRE(not read_gp_writer_raw.signature().has_value());
+
REQUIRE_NOTHROW(get_value(read_vtk_writer));
REQUIRE_NOTHROW(dynamic_cast<const VTKWriter&>(get_value(read_vtk_writer)));
const VTKWriter& read_vtk = dynamic_cast<const VTKWriter&>(get_value(read_vtk_writer));
@@ -111,7 +130,7 @@ TEST_CASE("checkpointing_IWriter", "[utils/checkpointing]")
REQUIRE(read_vtk.baseFilename() == vtk_filename);
REQUIRE(read_vtk.periodManager().has_value());
REQUIRE(read_vtk.periodManager().value().timePeriod() == 0.02);
- REQUIRE(read_vtk.periodManager().value().nextTime() == 1.02);
+ REQUIRE(read_vtk.periodManager().value().nextTime() == Catch::Approx(1.02));
REQUIRE(read_vtk.signature().has_value());
REQUIRE(read_vtk.signature().value() == p_vtk_writer->signature().value());
}