diff --git a/src/main.cpp b/src/main.cpp
index 2aee459b7de4503a5f61dfab71922e3cacb21716..a0ab842a3ca1b0961ce64fddfeb1fb04d68cf714 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -120,7 +120,9 @@ int main(int argc, char *argv[])
VTKWriter vtk_writer("mesh", 0.01);
while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -132,7 +134,9 @@ int main(int argc, char *argv[])
t += dt;
++iteration;
}
- vtk_writer.write(mesh, t, true); // forces last output
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
@@ -227,7 +231,9 @@ int main(int argc, char *argv[])
VTKWriter vtk_writer("mesh", 0.01);
while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -239,7 +245,9 @@ int main(int argc, char *argv[])
t += dt;
++iteration;
}
- vtk_writer.write(mesh, t, true); // forces last output
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
@@ -323,7 +331,9 @@ int main(int argc, char *argv[])
VTKWriter vtk_writer("mesh", 0.01);
while((t<tmax) and (iteration<itermax)) {
- vtk_writer.write(mesh, t);
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -334,7 +344,9 @@ int main(int argc, char *argv[])
t += dt;
++iteration;
}
- vtk_writer.write(mesh, t, true); // forces last output
+ vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
+ NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 60195176262fc8b62d2be7cdbaf55c01af5573d5..bf7e86ebca7d279a03586a6ad183b24887457b23 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -62,6 +62,16 @@ class ItemToItemMatrix
const ConnectivityMatrix& m_connectivity_matrix;
public:
+ auto numberOfEntries() const
+ {
+ return m_connectivity_matrix.numEntries();
+ }
+
+ const auto& entries() const
+ {
+ return m_connectivity_matrix.entries();
+ }
+
PASTIS_INLINE
auto operator[](const SourceItemId& source_id) const
{
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4f056098840f338f6afd869acad83f2f567d7d2f
--- /dev/null
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -0,0 +1,121 @@
+#ifndef OUTPUT_NAMED_ITEM_VALUE_SET_HPP
+#define OUTPUT_NAMED_ITEM_VALUE_SET_HPP
+
+#include <ItemValue.hpp>
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+
+#include <variant>
+#include <map>
+#include <string>
+
+template <typename DataType,
+ ItemType item_type>
+class NamedItemValue
+{
+ private:
+ std::string m_name;
+ ItemValue<const DataType,item_type> m_item_value;
+
+ public:
+ constexpr const std::string& name() const
+ {
+ return m_name;
+ }
+
+ constexpr const ItemValue<const DataType,item_type>& itemValue() const
+ {
+ return m_item_value;
+ }
+
+ NamedItemValue& operator=(const NamedItemValue&) = default;
+ NamedItemValue& operator=(NamedItemValue&&) = default;
+
+ NamedItemValue(const std::string& name, const ItemValue<DataType,item_type>& item_value)
+ : m_name(name),
+ m_item_value(item_value)
+ {
+ ;
+ }
+
+ NamedItemValue(const std::string& name, const ItemValue<const DataType,item_type>& item_value)
+ : m_name(name),
+ m_item_value(item_value)
+ {
+ ;
+ }
+
+ NamedItemValue(const NamedItemValue&) = default;
+ NamedItemValue(NamedItemValue&&) = default;
+ ~NamedItemValue() = default;
+};
+
+class OutputNamedItemValueSet
+{
+ public:
+ using ItemValueVariant = std::variant<NodeValue<const int>,
+ NodeValue<const long int>,
+ NodeValue<const double>,
+ NodeValue<const TinyVector<1,double>>,
+ NodeValue<const TinyVector<2,double>>,
+ NodeValue<const TinyVector<3,double>>,
+
+ CellValue<const int>,
+ CellValue<const long int>,
+ CellValue<const double>,
+ CellValue<const TinyVector<1,double>>,
+ CellValue<const TinyVector<2,double>>,
+ CellValue<const TinyVector<3,double>>
+ >;
+
+ private:
+ std::map<std::string, ItemValueVariant> m_name_itemvariant_map;
+
+ template <typename DataType,
+ ItemType item_type>
+ PASTIS_FORCEINLINE
+ constexpr void _doInsert(const NamedItemValue<DataType, item_type>& named_itemvalue)
+ {
+ if (m_name_itemvariant_map.find(named_itemvalue.name()) == m_name_itemvariant_map.end()) {
+ m_name_itemvariant_map[named_itemvalue.name()] = named_itemvalue.itemValue();
+ }
+ }
+
+ template <typename DataType,
+ ItemType item_type,
+ typename... Args>
+ PASTIS_FORCEINLINE
+ constexpr void _unpackVariadicInput(const NamedItemValue<DataType, item_type>& named_itemvalue,
+ Args&&... args)
+ {
+ _doInsert(named_itemvalue);
+ if constexpr (sizeof...(args) > 0) {
+ this->_unpackVariadicInput(std::forward<Args>(args)...);
+ }
+ }
+
+ public:
+ auto begin() const
+ {
+ return m_name_itemvariant_map.begin();
+ }
+
+ auto end() const
+ {
+ return m_name_itemvariant_map.end();
+ }
+
+ template <typename ...DataType,
+ ItemType ...item_type>
+ OutputNamedItemValueSet(NamedItemValue<DataType, item_type>... named_itemvalue)
+ {
+ _unpackVariadicInput(named_itemvalue...);
+ }
+
+ OutputNamedItemValueSet(const OutputNamedItemValueSet&) = default;
+ OutputNamedItemValueSet() = default;
+ ~OutputNamedItemValueSet() = default;
+};
+
+
+#endif // OUTPUT_NAMED_ITEM_VALUE_SET_HPP
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 276bad33d5624f018b9d512743f4f7695a6213e5..ffed42738c967dde7b799cd7607eece281ce1ed5 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -9,6 +9,7 @@
#include <IConnectivity.hpp>
#include <ItemValue.hpp>
+#include <OutputNamedItemValueSet.hpp>
class VTKWriter
{
@@ -18,9 +19,117 @@ class VTKWriter
double m_last_time;
const double m_time_period;
+ template <typename DataType> struct VTKType {};
+
+ template <typename DataType>
+ void _write_array(std::ofstream& os,
+ const std::string& name,
+ const Array<DataType>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\">\n";
+ for (typename Array<DataType>::index_type i=0;
+ i<item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <size_t N,
+ typename DataType>
+ void _write_array(std::ofstream& os,
+ const std::string& name,
+ const Array<TinyVector<N, DataType>>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N << "\">\n";
+ for (typename Array<DataType>::index_type i=0;
+ i<item_value.size(); ++i) {
+ for (size_t j=0; j<N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <typename DataType>
+ void _write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const DataType>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\">\n";
+ for (NodeId i=0; i<item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <size_t N,
+ typename DataType>
+ void _write_node_value(std::ofstream& os,
+ const std::string& name,
+ const NodeValue<const TinyVector<N, DataType>>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N << "\">\n";
+ for (NodeId i=0; i<item_value.size(); ++i) {
+ for (size_t j=0; j<N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <typename DataType>
+ void _write_node_value(std::ofstream&,
+ const std::string&,
+ const CellValue<const DataType>&) {}
+
+ template <typename DataType>
+ void _write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const DataType>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\">\n";
+ for (CellId i=0; i<item_value.size(); ++i) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i] << ' ';
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <size_t N,
+ typename DataType>
+ void _write_cell_value(std::ofstream& os,
+ const std::string& name,
+ const CellValue<const TinyVector<N, DataType>>& item_value)
+ {
+ os << "<DataArray type=\"" << VTKType<DataType>::name
+ << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N << "\">\n";
+ for (CellId i=0; i<item_value.size(); ++i) {
+ for (size_t j=0; j<N; ++j) {
+ // The following '+' enforces integer output for char types
+ os << +item_value[i][j] << ' ';
+ }
+ }
+ os << "\n</DataArray>\n";
+ }
+
+ template <typename DataType>
+ void _write_cell_value(std::ofstream&,
+ const std::string&,
+ const NodeValue<const DataType>&) {}
+
public:
template <typename MeshType>
void write(const MeshType& mesh,
+ const OutputNamedItemValueSet& output_named_item_value_set,
const double& time,
const bool& forced_output = false)
{
@@ -38,96 +147,98 @@ class VTKWriter
fout << "<UnstructuredGrid>\n";
fout << "<Piece NumberOfPoints=\""<< mesh.numberOfNodes()
<< "\" NumberOfCells=\"" << mesh.numberOfCells() << "\">\n";
-
+ fout << "<CellData>\n";
+ for(const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name=name](auto&& item_value) {
+ return this->_write_cell_value(fout, name, item_value);
+ }, item_value_variant);
+ }
+ fout << "</CellData>\n";
+ fout << "<PointData>\n";
+ for(const auto& [name, item_value_variant] : output_named_item_value_set) {
+ std::visit([&, name=name](auto&& item_value) {
+ return this->_write_node_value(fout, name, item_value);
+ }, item_value_variant);
+ }
+ fout << "</PointData>\n";
fout << "<Points>\n";
- fout << "<DataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n";
- using Rd = TinyVector<MeshType::dimension>;
- const NodeValue<const Rd>& xr = mesh.xr();
- if constexpr(MeshType::dimension == 1) {
- for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
- for (unsigned short i=0; i<1; ++i) {
- fout << xr[r][i] << ' ';
- }
- fout << "0 0 "; // VTK requires 3 components
- }
- } else if constexpr (MeshType::dimension == 2) {
+ {
+ using Rd = TinyVector<MeshType::dimension>;
+ const NodeValue<const Rd>& xr = mesh.xr();
+ Array<TinyVector<3>> positions(mesh.numberOfNodes());
for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
- for (unsigned short i=0; i<2; ++i) {
- fout << xr[r][i] << ' ';
+ for (unsigned short i=0; i<MeshType::dimension; ++i) {
+ positions[r][i] = xr[r][i];
}
- fout << "0 "; // VTK requires 3 components
- }
- } else {
- for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
- for (unsigned short i=0; i<3; ++i) {
- fout << xr[r][i] << ' ';
+ for (unsigned short i=MeshType::dimension; i<3; ++i) {
+ positions[r][i] = 0;
}
}
+ _write_array(fout, "Positions", positions);
}
- fout << '\n';
- fout << "</DataArray>\n";
fout << "</Points>\n";
fout << "<Cells>\n";
- fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
-
- const auto& cell_to_node_matrix
- = mesh.connectivity().cellToNodeMatrix();
+ {
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().cellToNodeMatrix();
+ Array<int> connectivity(cell_to_node_matrix.numberOfEntries());
- for (CellId j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- for (unsigned short r=0; r<cell_nodes.size(); ++r) {
- fout << cell_nodes[r] << ' ';
+ for (size_t i=0; i<cell_to_node_matrix.numberOfEntries(); ++i) {
+ connectivity[i] = cell_to_node_matrix.entries()[i];
}
+ _write_array(fout, "connectivity", connectivity);
}
- fout << '\n';
- fout << "</DataArray>\n";
- fout << "<DataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
{
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().cellToNodeMatrix();
+ Array<unsigned int> offsets(mesh.numberOfCells());
unsigned int offset=0;
for (CellId j=0; j<mesh.numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix[j];
offset += cell_nodes.size();
- fout << offset << ' ';
+ offsets[j]=offset;
}
+ _write_array(fout, "offsets", offsets);
}
- fout << '\n';
- fout << "</DataArray>\n";
-
- fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
- for (CellId j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- switch (cell_nodes.size()) {
- case 2: {
- fout << "3 ";
- break;
- }
- case 3: {
- fout << "5 ";
- break;
- }
- case 4: {
- if (mesh.meshDimension() == 3) {
- fout << "10 ";
- } else {
- fout << "9 ";
+
+ {
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().cellToNodeMatrix();
+ Array<int8_t> types(mesh.numberOfCells());
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ switch (cell_nodes.size()) {
+ case 2: {
+ types[j]=3;
+ break;
+ }
+ case 3: {
+ types[j]=5;
+ break;
+ }
+ case 4: {
+ if (mesh.meshDimension() == 3) {
+ types[j]=10;
+ } else {
+ types[j]=9;
+ }
+ break;
+ }
+ case 8: {
+ types[j]=12;
+ break;
+ }
+ default: {
+ types[j]=7;
+ break;
}
- break;
- }
- case 8: {
- fout << "12 ";
- break;
- }
- default: {
- fout << "7 ";
- break;
}
}
+ _write_array(fout, "types", types);
}
- fout << '\n';
- fout << "</DataArray>\n";
fout << "</Cells>\n";
fout << "</Piece>\n";
@@ -136,6 +247,7 @@ class VTKWriter
m_file_number++;
}
+
VTKWriter(const std::string& base_filename,
const double time_period)
: m_base_filename(base_filename),
@@ -147,4 +259,55 @@ class VTKWriter
~VTKWriter() = default;
};
+template <>
+struct VTKWriter::VTKType<int8_t>
+{
+ inline const static std::string name{"Int8"};
+};
+
+template <> struct VTKWriter::VTKType<uint8_t>
+{
+ inline const static std::string name{"UInt8"};
+};
+
+template <> struct VTKWriter::VTKType<int16_t>
+{
+ inline const static std::string name{"Int16"};
+};
+
+template <> struct VTKWriter::VTKType<uint16_t>
+{
+ inline const static std::string name{"UInt16"};
+};
+
+template <> struct VTKWriter::VTKType<int32_t>
+{
+ inline const static std::string name{"Int32"};
+};
+
+template <> struct VTKWriter::VTKType<uint32_t>
+{
+ inline const static std::string name{"UInt32"};
+};
+
+template <> struct VTKWriter::VTKType<int64_t>
+{
+ inline const static std::string name{"Int64"};
+};
+
+template <> struct VTKWriter::VTKType<uint64_t>
+{
+ inline const static std::string name{"UInt64"};
+};
+
+template <> struct VTKWriter::VTKType<float>
+{
+ inline const static std::string name{"Float32"};
+};
+
+template <> struct VTKWriter::VTKType<double>
+{
+ inline const static std::string name{"Float64"};
+};
+
#endif // VTK_WRITER_HPP