From 7c3e7222921ca53db7b25bc0bca447ffa285edd5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Wed, 27 Jul 2022 08:53:03 +0200
Subject: [PATCH] Add handling of ItemValue and SubItemValuePerItem into the
language
- SubItemValuePerItem is associated to the pugs abstract type sub_item_value
- ItemValue is associated to the pugs abstract type item_value
This is done to enrich the data flow between functions but cannot
really be manipulated in the language.
However, one can interpolate to ItemData and can post process
them (only CellValue and NodeValue are supported by VTK and Gnuplot
writers).
- doc and tests have been updated
---
doc/userdoc.org | 184 +++++-
src/language/modules/MeshModule.cpp | 50 ++
src/language/modules/MeshModule.hpp | 14 +
src/language/modules/SchemeModule.cpp | 76 ++-
src/language/modules/WriterModule.cpp | 66 ++-
src/language/modules/WriterModule.hpp | 4 +-
src/language/utils/CMakeLists.txt | 1 +
.../ItemValueVariantFunctionInterpoler.cpp | 145 +++++
.../ItemValueVariantFunctionInterpoler.hpp | 38 ++
src/mesh/ItemValue.hpp | 7 +
src/mesh/ItemValueVariant.hpp | 113 ++++
src/mesh/MeshData.hpp | 35 +-
src/mesh/SubItemValuePerItemVariant.hpp | 163 ++++++
src/output/GnuplotWriter.cpp | 46 +-
src/output/GnuplotWriter.hpp | 7 +-
src/output/GnuplotWriter1D.cpp | 58 +-
src/output/GnuplotWriter1D.hpp | 13 +-
src/output/INamedDiscreteData.hpp | 27 +
src/output/IWriter.hpp | 26 +-
src/output/NamedDiscreteFunction.hpp | 12 +-
src/output/NamedItemValueVariant.hpp | 47 ++
src/output/OutputNamedItemValueSet.hpp | 36 ++
src/output/VTKWriter.cpp | 34 +-
src/output/VTKWriter.hpp | 7 +-
src/output/WriterBase.cpp | 346 ++++++++---
src/output/WriterBase.hpp | 44 +-
src/scheme/AcousticSolver.cpp | 216 ++++---
src/scheme/AcousticSolver.hpp | 38 +-
src/scheme/DiscreteFunctionInterpoler.cpp | 1 -
tests/CMakeLists.txt | 3 +
tests/test_ItemValueVariant.cpp | 101 ++++
...est_ItemValueVariantFunctionInterpoler.cpp | 547 ++++++++++++++++++
tests/test_SubItemValuePerItemVariant.cpp | 165 ++++++
33 files changed, 2406 insertions(+), 264 deletions(-)
create mode 100644 src/language/utils/ItemValueVariantFunctionInterpoler.cpp
create mode 100644 src/language/utils/ItemValueVariantFunctionInterpoler.hpp
create mode 100644 src/mesh/ItemValueVariant.hpp
create mode 100644 src/mesh/SubItemValuePerItemVariant.hpp
create mode 100644 src/output/INamedDiscreteData.hpp
create mode 100644 src/output/NamedItemValueVariant.hpp
create mode 100644 tests/test_ItemValueVariant.cpp
create mode 100644 tests/test_ItemValueVariantFunctionInterpoler.cpp
create mode 100644 tests/test_SubItemValuePerItemVariant.cpp
diff --git a/doc/userdoc.org b/doc/userdoc.org
index 12ce427db..a91d0d5b8 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2236,20 +2236,18 @@ follows the similar rules.
Let us give an example
#+NAME: functions-lifetime
#+BEGIN_SRC pugs :exports both :results output
- let f: R -> R, x -> 2.3*x+2;
{
- let f: R -> R, x -> 1.25*x+3;
- cout << "block 1: f(2) = " << f(2) << "\n";
{
let f: R -> R, x -> 3.25*x-0.3;
- cout << "block 2: f(2) = " << f(2) << "\n";
+ cout << "block 1.1: f(2) = " << f(2) << "\n";
}
+ let f: R -> R, x -> 1.25*x+3;
+ cout << "block 1: f(2) = " << f(2) << "\n";
}
for (let i:N, i = 1; i<=2; ++i) {
let f: R -> R, x -> i*x+2;
cout << "for(i=" << i << "): f(2) = " << f(2) << "\n";
}
- cout << "global: f(2) = " << f(2) << "\n";
#+END_SRC
Running this example produces
#+results: functions-lifetime
@@ -2627,6 +2625,41 @@ The information produced concerns
- the number of edges and their references,
- and the number of nodes and their references.
+
+***** ~item_type~
+
+This type is used to designate kinds of items (cell, face, edge or node).
+
+***** ~item_value~
+
+The type ~item_value~ is an abstract type use to designate arrays of
+values defined on the entities of a ~mesh~. Actually, these values are
+not associated to the mesh itself but on the *connectivity*. The values
+on the entities can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~
+and ~R^3x3~. Entities themselves can be cells, faces, edges or nodes.
+
+These variables are used to pass data from one function to another.
+
+#+BEGIN_warning
+By now, no mathematical operation is defined on ~item_value~ variables.
+#+END_warning
+
+Moreover, ~item_value~ variables can be post processed. Observe that
+edges or faces values cannot be post processed since neither ~VTK~ nor
+~Gnuplot~ can handle these data.
+
+***** ~sub_item_value~
+
+This abstract type is handle values defined on the sub items of the
+items of a ~mesh~. Similarly these values are attached to a *connectivity*
+and not to a mesh. Values can bed of type The values on the entities
+can be of type ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~ or ~R^3~. An example of
+~sub_item_value~ is the $\mathbf{C}_{jr}$ vectors which are defined at each
+node of each cell.
+
+These variables are used to pass data from one function to
+another. They cannot be post processed.
+
**** ~mesh~ provided functions
***** Boundary descriptor functions
@@ -2897,6 +2930,35 @@ The median dual mesh construction is not yet implemented in 3d and not
available in parallel
#+END_warning
+***** Item types
+The following functions are used to refer to ~item_type~
+- ~cell: void -> item_type~
+- ~face: void -> item_type~
+- ~edge: void -> item_type~
+- ~node: void -> item_type~
+
+***** ~interpolate: mesh*item_type*function -> item_value~
+
+This function is used to create an ~item_value~ by interpolating a user
+function to a given ~item_type~ using a mesh. The interpolation takes
+place at the center of the items.
+
+#+BEGIN_SRC pugs :exports both :results none
+ import mesh;
+
+ let m:mesh, m = cartesianMesh([0,0], [1,1], (10,10));
+ let f:R^2 -> R, x -> 2*x[0]-x[1];
+
+ let node_values:item_value, node_values = interpolate(m, node(), f);
+ let face_values:item_value, face_values = interpolate(m, face(), f);
+ let cell_values:item_value, cell_values = interpolate(m, cell(), f);
+#+END_SRC
+In this example, we set three arrays defined at all nodes, all the
+ faces and all the cells of the underlying connectivity of the mesh ~m~.
+- ~node_values~ interpolates ~f~ at the nodes of ~m~,
+- ~face_values~ interpolates ~f~ at the centers of the faces of ~m~,
+- ~cell_values~ interpolates ~f~ at the centers of the cells of ~m~.
+
***** ~transform: mesh*function -> mesh~
This function allows to compute a new mesh as the transformation of a
@@ -3426,6 +3488,8 @@ vector function.
Here we defined a $\vec{\mathbb{P}}_0(\mathbb{R})$ discrete function of
dimension 3.
+
+
****** ~integrate: mesh*quadrature*function -> Vh~ <<integrate-classic>>
This function integrates the given user function in each cell of a
@@ -4293,6 +4357,96 @@ of ~fh~.
#+ATTR_HTML: :width 300px;
#+RESULTS: writer-example2-2-img
+***** ~name_output: item_value*string -> output~
+
+This function give a name to an ~item_value~.
+
+#+BEGIN_SRC pugs :exports both :results none
+ import mesh;
+ import writer;
+
+ let m:mesh, m = cartesianMesh([0], [1], 100);
+
+ let x2:R^1 -> R, x -> x[0]*x[0];
+ let cell_values:item_value, cell_values = interpolate(m, node(), x2);
+
+ write(gnuplot_writer("writer-item-value-example1"), m, name_output(cell_values, "x^2"));
+#+END_SRC
+
+#+NAME: writer-item-value-example1-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-item-value-example1.png")
+ reset
+ set grid
+ set border
+ unset key
+ set xtics
+ set ytics
+ set square
+ set terminal png truecolor enhanced size 640,640
+ plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example1.gnu)' lw 2 w l
+#+END_SRC
+
+#+CAPTION: Simple illustration of the output.
+#+NAME: fig:writer-item-value-example1
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-item-value-example1-img
+
+
+#+BEGIN_warning
+One observes that we supplied a mesh to the writer function. The
+reason for that is that ~item_value~ refers a connectivity. It requires
+a ~mesh~ with the same connectivity to be written.
+
+However, is a discrete function (of type ~Vh~) build on the same
+connectivity is provided, there is no need to specify the ~mesh~.
+#+END_warning
+
+Let us illustrate it by an important second example.
+
+#+BEGIN_SRC pugs :exports both :results none
+ import mesh;
+ import scheme;
+ import writer;
+ import math;
+
+ let m:mesh, m = cartesianMesh([0], [1], 100);
+
+ let f:R^1 -> R, x -> 2*x[0]*x[0];
+ let g:R^1 -> R, x -> pow(sin(5*x[0]*x[0]),2);
+
+ let fh:Vh, fh = interpolate(m, P0(), f);
+ let cell_value:item_value, cell_value = interpolate(m, cell(), g);
+
+ let output_list:(output),
+ output_list = (name_output(fh, "f"),
+ name_output(cell_value, "g"));
+
+ write(gnuplot_1d_writer("writer-item-value-example2"), output_list);
+#+END_SRC
+
+Running this code produces the gnuplot file displayed in Figure
+[[fig:writer-item-value-example2]].
+
+#+NAME: writer-item-value-example2-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-item-value-example2.png")
+ reset
+ set grid
+ set border
+ set key
+ set xtics
+ set ytics
+ set square
+ set terminal png truecolor enhanced size 640,640
+ plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-item-value-example2.gnu)' u 1:3 lw 2 t "g" w l
+#+END_SRC
+
+#+CAPTION: The ~item_value~ is plotted on the same mesh as ~fh~.
+#+NAME: fig:writer-item-value-example2
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-item-value-example2-img
+
***** ~gnuplot~ writers <<gnuplot-writers>>
There are two ~gnuplot~ writers. One is dedicated to output of dimension
@@ -4652,6 +4806,14 @@ Trying to use variables defined on different meshes
gives the runtime error
#+results: cannot-use-diffrent-meshes-in-writer
+****** ~write: writer*mesh*(output) -> void~
+
+This variation permits to specify a ~mesh~. Its main utility is to allow
+to write data in the case where ~output~ list contains only ~item_value~
+data which are not related to a mesh but to a *connectivity*. One checks
+that the connectivity of the ~output~ data are the same as the
+connectivity of the ~mesh~.
+
****** ~write: writer*(output)*R -> void~
This ~write~ function is used to save time series data. The real
@@ -4680,6 +4842,13 @@ Trying to use a non time series ~writer~
gives the runtime error
#+results: cannot-use-non-time-series-writer
+****** ~write: writer*mesh*(output)*R -> void~
+
+This is a variation of the previous function where a mesh is
+specified. The function checks the validity of the ~output~ list:
+functions ~Vh~ must defined on the given ~mesh~ and ~item_value~ data mush
+share the same connectivity with the ~mesh~.
+
****** ~force_write: writer*(output)*R -> void~
One probably noticed that using the ~write~ function with a time series
@@ -4744,6 +4913,11 @@ times:
#+RESULTS: times-in-gp-1d-series-force
The last post processing time is now 1.
+****** ~force_write: writer*mesh*(output)*R -> void~
+
+This variation permits to specify a ~mesh~. Compatibility of output data
+with the ~mesh~ is checked.
+
****** ~write_mesh: writer*mesh -> void~ <<write-mesh>>.
This function just saves a ~mesh~ using a ~writer~. The ~gnuplot_1d_writer~
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 4fc75f4f2..9c5cf3aec 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -5,6 +5,7 @@
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/FunctionTable.hpp>
+#include <language/utils/ItemValueVariantFunctionInterpoler.hpp>
#include <language/utils/OStream.hpp>
#include <language/utils/OperatorRepository.hpp>
#include <language/utils/SymbolTable.hpp>
@@ -15,6 +16,7 @@
#include <mesh/GmshReader.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
#include <mesh/IZoneDescriptor.hpp>
+#include <mesh/ItemValueVariant.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshRelaxer.hpp>
#include <mesh/MeshTransformer.hpp>
@@ -22,6 +24,7 @@
#include <mesh/NamedZoneDescriptor.hpp>
#include <mesh/NumberedBoundaryDescriptor.hpp>
#include <mesh/NumberedZoneDescriptor.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
#include <utils/Exceptions.hpp>
#include <Kokkos_Core.hpp>
@@ -32,6 +35,43 @@ MeshModule::MeshModule()
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IZoneDescriptor>>);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const ItemType>>);
+
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const ItemValueVariant>>);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const SubItemValuePerItemVariant>>);
+
+ this->_addBuiltinFunction("cell", std::function(
+
+ []() -> std::shared_ptr<const ItemType> {
+ return std::make_shared<ItemType>(ItemType::cell);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("face", std::function(
+
+ []() -> std::shared_ptr<const ItemType> {
+ return std::make_shared<ItemType>(ItemType::face);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("edge", std::function(
+
+ []() -> std::shared_ptr<const ItemType> {
+ return std::make_shared<ItemType>(ItemType::edge);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("node", std::function(
+
+ []() -> std::shared_ptr<const ItemType> {
+ return std::make_shared<ItemType>(ItemType::node);
+ }
+
+ ));
+
this->_addBuiltinFunction("readGmsh", std::function(
[](const std::string& file_name) -> std::shared_ptr<const IMesh> {
@@ -74,6 +114,16 @@ MeshModule::MeshModule()
));
+ this->_addBuiltinFunction("interpolate",
+ std::function(
+
+ [](std::shared_ptr<const IMesh> mesh, std::shared_ptr<const ItemType> item_type,
+ const FunctionSymbolId& function_id) -> std::shared_ptr<const ItemValueVariant> {
+ return ItemValueVariantFunctionInterpoler{mesh, *item_type, function_id}.interpolate();
+ }
+
+ ));
+
this->_addBuiltinFunction("transform", std::function(
[](std::shared_ptr<const IMesh> p_mesh,
diff --git a/src/language/modules/MeshModule.hpp b/src/language/modules/MeshModule.hpp
index 904741258..bf9f0b714 100644
--- a/src/language/modules/MeshModule.hpp
+++ b/src/language/modules/MeshModule.hpp
@@ -20,6 +20,20 @@ template <>
inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IZoneDescriptor>> =
ASTNodeDataType::build<ASTNodeDataType::type_id_t>("zone");
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemType>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_type");
+
+class ItemValueVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const ItemValueVariant>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("item_value");
+
+class SubItemValuePerItemVariant;
+template <>
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const SubItemValuePerItemVariant>> =
+ ASTNodeDataType::build<ASTNodeDataType::type_id_t>("sub_item_value");
+
class MeshModule : public BuiltinModule
{
public:
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 01dcc66eb..e3c6f1b96 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -288,6 +288,24 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("glace_fluxes", std::function(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list)
+ -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>> {
+ return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
+ .solver()
+ .compute_fluxes(AcousticSolverHandler::SolverType::Glace, rho, c, u,
+ p, bc_descriptor_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("glace_solver",
std::function(
@@ -302,14 +320,29 @@ SchemeModule::SchemeModule()
-> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>> {
- return AcousticSolverHandler{AcousticSolverHandler::SolverType::Glace,
- rho,
- c,
- u,
- p,
- bc_descriptor_list}
+ return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
.solver()
- .apply(dt, rho, u, E);
+ .apply(AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
+ bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("eucclhyd_fluxes",
+ std::function(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list)
+ -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>> {
+ return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
+ .solver()
+ .compute_fluxes(AcousticSolverHandler::SolverType::Eucclhyd, rho, c, u, p,
+ bc_descriptor_list);
}
));
@@ -328,14 +361,29 @@ SchemeModule::SchemeModule()
-> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>> {
- return AcousticSolverHandler{AcousticSolverHandler::SolverType::Eucclhyd,
- rho,
- c,
- u,
- p,
- bc_descriptor_list}
+ return AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
+ .solver()
+ .apply(AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
+ bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("apply_acoustic_fluxes",
+ std::function(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& rho, //
+ const std::shared_ptr<const IDiscreteFunction>& u, //
+ const std::shared_ptr<const IDiscreteFunction>& E, //
+ const std::shared_ptr<const ItemValueVariant>& ur, //
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, //
+ const double& dt)
+ -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>> {
+ return AcousticSolverHandler{getCommonMesh({rho, u, E})} //
.solver()
- .apply(dt, rho, u, E);
+ .apply_fluxes(dt, rho, u, E, ur, Fjr);
}
));
diff --git a/src/language/modules/WriterModule.cpp b/src/language/modules/WriterModule.cpp
index 89d25001b..6523d56a9 100644
--- a/src/language/modules/WriterModule.cpp
+++ b/src/language/modules/WriterModule.cpp
@@ -4,11 +4,13 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/GmshReader.hpp>
+#include <mesh/ItemValueVariant.hpp>
#include <mesh/Mesh.hpp>
#include <output/GnuplotWriter.hpp>
#include <output/GnuplotWriter1D.hpp>
#include <output/IWriter.hpp>
#include <output/NamedDiscreteFunction.hpp>
+#include <output/NamedItemValueVariant.hpp>
#include <output/VTKWriter.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/IDiscreteFunction.hpp>
@@ -16,7 +18,7 @@
WriterModule::WriterModule()
{
- this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>>);
+ this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const INamedDiscreteData>>);
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IWriter>>);
@@ -71,15 +73,24 @@ WriterModule::WriterModule()
));
- this->_addBuiltinFunction("name_output",
- std::function(
+ this->_addBuiltinFunction("name_output", std::function(
- [](std::shared_ptr<const IDiscreteFunction> discrete_function,
- const std::string& name) -> std::shared_ptr<const NamedDiscreteFunction> {
- return std::make_shared<const NamedDiscreteFunction>(discrete_function, name);
- }
+ [](std::shared_ptr<const IDiscreteFunction> discrete_function,
+ const std::string& name) -> std::shared_ptr<const INamedDiscreteData> {
+ return std::make_shared<const NamedDiscreteFunction>(discrete_function,
+ name);
+ }
- ));
+ ));
+ this->_addBuiltinFunction("name_output", std::function(
+
+ [](std::shared_ptr<const ItemValueVariant> item_value_variant,
+ const std::string& name) -> std::shared_ptr<const INamedDiscreteData> {
+ return std::make_shared<const NamedItemValueVariant>(item_value_variant,
+ name);
+ }
+
+ ));
this->_addBuiltinFunction("write_mesh",
std::function(
@@ -93,7 +104,7 @@ WriterModule::WriterModule()
this->_addBuiltinFunction("write", std::function(
[](std::shared_ptr<const IWriter> writer,
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
named_discrete_function_list) -> void {
writer->write(named_discrete_function_list);
}
@@ -103,7 +114,7 @@ WriterModule::WriterModule()
this->_addBuiltinFunction("write", std::function(
[](std::shared_ptr<const IWriter> writer,
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
named_discrete_function_list,
const double& time) -> void {
writer->writeIfNeeded(named_discrete_function_list, time);
@@ -114,13 +125,46 @@ WriterModule::WriterModule()
this->_addBuiltinFunction("force_write", std::function(
[](std::shared_ptr<const IWriter> writer,
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>&
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
named_discrete_function_list,
const double& time) -> void {
writer->writeForced(named_discrete_function_list, time);
}
));
+
+ this->_addBuiltinFunction("write", std::function(
+
+ [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
+ named_discrete_function_list) -> void {
+ writer->writeOnMesh(mesh, named_discrete_function_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("write", std::function(
+
+ [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
+ named_discrete_function_list,
+ const double& time) -> void {
+ writer->writeOnMeshIfNeeded(mesh, named_discrete_function_list, time);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("force_write",
+ std::function(
+
+ [](std::shared_ptr<const IWriter> writer, std::shared_ptr<const IMesh> mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>&
+ named_discrete_function_list,
+ const double& time) -> void {
+ writer->writeOnMeshForced(mesh, named_discrete_function_list, time);
+ }
+
+ ));
}
void
diff --git a/src/language/modules/WriterModule.hpp b/src/language/modules/WriterModule.hpp
index adcd92530..f61eafb0e 100644
--- a/src/language/modules/WriterModule.hpp
+++ b/src/language/modules/WriterModule.hpp
@@ -6,13 +6,13 @@
#include <utils/PugsMacros.hpp>
class OutputNamedItemValueSet;
-class NamedDiscreteFunction;
+class INamedDiscreteData;
class IDiscreteFunction;
#include <string>
template <>
-inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const NamedDiscreteFunction>> =
+inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const INamedDiscreteData>> =
ASTNodeDataType::build<ASTNodeDataType::type_id_t>("output");
class IWriter;
diff --git a/src/language/utils/CMakeLists.txt b/src/language/utils/CMakeLists.txt
index b76655995..0fe1c798b 100644
--- a/src/language/utils/CMakeLists.txt
+++ b/src/language/utils/CMakeLists.txt
@@ -29,6 +29,7 @@ add_library(PugsLanguageUtils
FunctionSymbolId.cpp
IncDecOperatorRegisterForN.cpp
IncDecOperatorRegisterForZ.cpp
+ ItemValueVariantFunctionInterpoler.cpp
OFStream.cpp
OperatorRepository.cpp
UnaryOperatorRegisterForB.cpp
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
new file mode 100644
index 000000000..018dd7a98
--- /dev/null
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
@@ -0,0 +1,145 @@
+#include <language/utils/ItemValueVariantFunctionInterpoler.hpp>
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <memory>
+
+template <size_t Dimension, typename DataType, typename ValueType>
+std::shared_ptr<ItemValueVariant>
+ItemValueVariantFunctionInterpoler::_interpolate() const
+{
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(m_mesh);
+ using MeshDataType = MeshData<Dimension>;
+
+ switch (m_item_type) {
+ case ItemType::cell: {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+ return std::make_shared<ItemValueVariant>(
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(m_function_id,
+ mesh_data.xj()));
+ }
+ case ItemType::face: {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+ return std::make_shared<ItemValueVariant>(
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::face>(m_function_id,
+ mesh_data.xl()));
+ }
+ case ItemType::edge: {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+ return std::make_shared<ItemValueVariant>(
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::edge>(m_function_id,
+ mesh_data.xe()));
+ }
+ case ItemType::node: {
+ return std::make_shared<ItemValueVariant>(
+ InterpolateItemValue<DataType(TinyVector<Dimension>)>::template interpolate<ItemType::node>(m_function_id,
+ p_mesh->xr()));
+ }
+ default: {
+ throw UnexpectedError("invalid item type");
+ }
+ }
+}
+
+template <size_t Dimension>
+std::shared_ptr<ItemValueVariant>
+ItemValueVariantFunctionInterpoler::_interpolate() const
+{
+ const auto& function_descriptor = m_function_id.descriptor();
+ Assert(function_descriptor.domainMappingNode().children[1]->m_data_type == ASTNodeDataType::typename_t);
+
+ const ASTNodeDataType& data_type = function_descriptor.domainMappingNode().children[1]->m_data_type.contentType();
+
+ switch (data_type) {
+ case ASTNodeDataType::bool_t: {
+ return this->_interpolate<Dimension, bool, double>();
+ }
+ case ASTNodeDataType::unsigned_int_t: {
+ return this->_interpolate<Dimension, uint64_t, double>();
+ }
+ case ASTNodeDataType::int_t: {
+ return this->_interpolate<Dimension, int64_t, double>();
+ }
+ case ASTNodeDataType::double_t: {
+ return this->_interpolate<Dimension, double>();
+ }
+ case ASTNodeDataType::vector_t: {
+ switch (data_type.dimension()) {
+ case 1: {
+ return this->_interpolate<Dimension, TinyVector<1>>();
+ }
+ case 2: {
+ return this->_interpolate<Dimension, TinyVector<2>>();
+ }
+ case 3: {
+ return this->_interpolate<Dimension, TinyVector<3>>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+ case ASTNodeDataType::matrix_t: {
+ Assert(data_type.numberOfColumns() == data_type.numberOfRows(), "undefined matrix type");
+ switch (data_type.numberOfColumns()) {
+ case 1: {
+ return this->_interpolate<Dimension, TinyMatrix<1>>();
+ }
+ case 2: {
+ return this->_interpolate<Dimension, TinyMatrix<2>>();
+ }
+ case 3: {
+ return this->_interpolate<Dimension, TinyMatrix<3>>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid vector dimension " << rang::fgB::red << data_type.dimension() << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+ // LCOV_EXCL_START
+ default: {
+ std::ostringstream os;
+ os << "invalid interpolation value type: " << rang::fgB::red << dataTypeName(data_type) << rang::style::reset;
+
+ throw UnexpectedError(os.str());
+ }
+ // LCOV_EXCL_STOP
+ }
+}
+
+std::shared_ptr<ItemValueVariant>
+ItemValueVariantFunctionInterpoler::interpolate() const
+{
+ switch (m_mesh->dimension()) {
+ case 1: {
+ return this->_interpolate<1>();
+ }
+ case 2: {
+ return this->_interpolate<2>();
+ }
+ case 3: {
+ return this->_interpolate<3>();
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.hpp b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
new file mode 100644
index 000000000..4e03f72fc
--- /dev/null
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
@@ -0,0 +1,38 @@
+#ifndef ITEM_VALUE_VARIANT_FUNCTION_INTERPOLER_HPP
+#define ITEM_VALUE_VARIANT_FUNCTION_INTERPOLER_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <mesh/IZoneDescriptor.hpp>
+#include <mesh/ItemType.hpp>
+#include <mesh/ItemValueVariant.hpp>
+
+class ItemValueVariantFunctionInterpoler
+{
+ private:
+ std::shared_ptr<const IMesh> m_mesh;
+ const ItemType m_item_type;
+ const FunctionSymbolId m_function_id;
+
+ template <size_t Dimension, typename DataType, typename ValueType = DataType>
+ std::shared_ptr<ItemValueVariant> _interpolate() const;
+
+ template <size_t Dimension>
+ std::shared_ptr<ItemValueVariant> _interpolate() const;
+
+ public:
+ std::shared_ptr<ItemValueVariant> interpolate() const;
+
+ ItemValueVariantFunctionInterpoler(const std::shared_ptr<const IMesh>& mesh,
+ const ItemType& item_type,
+ const FunctionSymbolId& function_id)
+ : m_mesh{mesh}, m_item_type{item_type}, m_function_id{function_id}
+ {}
+
+ ItemValueVariantFunctionInterpoler(const ItemValueVariantFunctionInterpoler&) = delete;
+ ItemValueVariantFunctionInterpoler(ItemValueVariantFunctionInterpoler&&) = delete;
+
+ ~ItemValueVariantFunctionInterpoler() = default;
+};
+
+#endif // ITEM_VALUE_VARIANT_FUNCTION_INTERPOLER_HPP
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index 496c33242..da6bb88eb 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -43,6 +43,13 @@ class ItemValue
friend ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityWeakPtr>;
public:
+ // This is not the correct way to look at ItemValue, use with care
+ Array<const DataType>
+ arrayView() const
+ {
+ return m_values;
+ }
+
[[nodiscard]] friend PUGS_INLINE ItemValue<std::remove_const_t<DataType>, item_type, ConnectivityPtr>
copy(const ItemValue<DataType, item_type, ConnectivityPtr>& source)
{
diff --git a/src/mesh/ItemValueVariant.hpp b/src/mesh/ItemValueVariant.hpp
new file mode 100644
index 000000000..1bfff8f84
--- /dev/null
+++ b/src/mesh/ItemValueVariant.hpp
@@ -0,0 +1,113 @@
+#ifndef ITEM_VALUE_VARIANT_HPP
+#define ITEM_VALUE_VARIANT_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/ItemValue.hpp>
+#include <utils/Exceptions.hpp>
+
+class ItemValueVariant
+{
+ private:
+ using Variant = std::variant<NodeValue<const bool>,
+ NodeValue<const long int>,
+ NodeValue<const unsigned long int>,
+ NodeValue<const double>,
+ NodeValue<const TinyVector<1, double>>,
+ NodeValue<const TinyVector<2, double>>,
+ NodeValue<const TinyVector<3, double>>,
+ NodeValue<const TinyMatrix<1, 1, double>>,
+ NodeValue<const TinyMatrix<2, 2, double>>,
+ NodeValue<const TinyMatrix<3, 3, double>>,
+
+ EdgeValue<const bool>,
+ EdgeValue<const long int>,
+ EdgeValue<const unsigned long int>,
+ EdgeValue<const double>,
+ EdgeValue<const TinyVector<1, double>>,
+ EdgeValue<const TinyVector<2, double>>,
+ EdgeValue<const TinyVector<3, double>>,
+ EdgeValue<const TinyMatrix<1, 1, double>>,
+ EdgeValue<const TinyMatrix<2, 2, double>>,
+ EdgeValue<const TinyMatrix<3, 3, double>>,
+
+ FaceValue<const bool>,
+ FaceValue<const long int>,
+ FaceValue<const unsigned long int>,
+ FaceValue<const double>,
+ FaceValue<const TinyVector<1, double>>,
+ FaceValue<const TinyVector<2, double>>,
+ FaceValue<const TinyVector<3, double>>,
+ FaceValue<const TinyMatrix<1, 1, double>>,
+ FaceValue<const TinyMatrix<2, 2, double>>,
+ FaceValue<const TinyMatrix<3, 3, double>>,
+
+ CellValue<const bool>,
+ CellValue<const long int>,
+ CellValue<const unsigned long int>,
+ CellValue<const double>,
+ CellValue<const TinyVector<1, double>>,
+ CellValue<const TinyVector<2, double>>,
+ CellValue<const TinyVector<3, double>>,
+ CellValue<const TinyMatrix<1, 1, double>>,
+ CellValue<const TinyMatrix<2, 2, double>>,
+ CellValue<const TinyMatrix<3, 3, double>>>;
+
+ Variant m_item_value;
+
+ public:
+ PUGS_INLINE
+ const Variant&
+ itemValue() const
+ {
+ return m_item_value;
+ }
+
+ template <typename ItemValueT>
+ PUGS_INLINE auto
+ get() const
+ {
+ using DataType = typename ItemValueT::data_type;
+ constexpr ItemType item_type = ItemValueT::item_t;
+
+ if constexpr (std::is_same_v<ItemValueT, ItemValue<DataType, item_type>> or
+ std::is_same_v<ItemValueT, ItemValue<const DataType, item_type>> or
+ std::is_same_v<ItemValueT, WeakItemValue<DataType, item_type>> or
+ std::is_same_v<ItemValueT, WeakItemValue<const DataType, item_type>>) {
+ if (not std::holds_alternative<ItemValue<const DataType, item_type>>(this->m_item_value)) {
+ throw NormalError("invalid ItemValue type");
+ }
+ return std::get<ItemValue<const DataType, item_type>>(this->itemValue());
+ } else {
+ static_assert(std::is_same_v<ItemValueT, ItemValueT>, "invalid template argument");
+ }
+ }
+
+ template <typename DataType, ItemType item_type>
+ ItemValueVariant(const ItemValue<DataType, item_type>& item_value)
+ : m_item_value{ItemValue<const DataType, item_type>{item_value}}
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, bool> or //
+ std::is_same_v<std::remove_const_t<DataType>, long int> or //
+ std::is_same_v<std::remove_const_t<DataType>, unsigned long int> or //
+ std::is_same_v<std::remove_const_t<DataType>, double> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<1, 1, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<2, 2, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyMatrix<3, 3, double>>,
+ "ItemValue with this DataType is not allowed in variant");
+ }
+
+ ItemValueVariant& operator=(ItemValueVariant&&) = default;
+ ItemValueVariant& operator=(const ItemValueVariant&) = default;
+
+ ItemValueVariant(const ItemValueVariant&) = default;
+ ItemValueVariant(ItemValueVariant&&) = default;
+
+ ItemValueVariant() = delete;
+ ~ItemValueVariant() = default;
+};
+
+#endif // ITEM_VALUE_VARIANT_HPP
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index ada0f8fbd..9b0c89c8a 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -45,6 +45,7 @@ class MeshData : public IMeshData
FaceValue<const Rd> m_Nl;
FaceValue<const Rd> m_nl;
FaceValue<const double> m_ll;
+ EdgeValue<const Rd> m_xe;
PUGS_INLINE
void
@@ -136,7 +137,7 @@ class MeshData : public IMeshData
_computeFaceIsobarycenter()
{ // Computes vertices isobarycenter
if constexpr (Dimension == 1) {
- static_assert(Dimension != 1, "xl does not make sense in 1d");
+ m_xl = FaceValue<const Rd>{m_mesh.connectivity(), m_mesh.xr().arrayView()};
} else {
const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -155,6 +156,28 @@ class MeshData : public IMeshData
}
}
+ PUGS_INLINE
+ void
+ _computeEdgeCenter()
+ { // Computes vertices isobarycenter
+ if constexpr (Dimension < 3) {
+ m_xe = EdgeValue<const Rd>(m_mesh.connectivity(), this->xl().arrayView());
+ } else {
+ const NodeValue<const Rd>& xr = m_mesh.xr();
+
+ const auto& edge_to_node_matrix = m_mesh.connectivity().edgeToNodeMatrix();
+ EdgeValue<Rd> xe(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfEdges(), PUGS_LAMBDA(EdgeId edge_id) {
+ const auto& edge_nodes = edge_to_node_matrix[edge_id];
+ const NodeId node0 = edge_nodes[0];
+ const NodeId node1 = edge_nodes[1];
+ xe[edge_id] = 0.5 * (xr[node0] + xr[node1]);
+ });
+ m_xe = xe;
+ }
+ }
+
PUGS_INLINE
void
_computeCellIsoBarycenter()
@@ -614,6 +637,16 @@ class MeshData : public IMeshData
return m_xl;
}
+ PUGS_INLINE
+ EdgeValue<const Rd>
+ xe()
+ {
+ if (not m_xe.isBuilt()) {
+ this->_computeEdgeCenter();
+ }
+ return m_xe;
+ }
+
PUGS_INLINE
CellValue<const Rd>
cellIsoBarycenter()
diff --git a/src/mesh/SubItemValuePerItemVariant.hpp b/src/mesh/SubItemValuePerItemVariant.hpp
new file mode 100644
index 000000000..d711f400a
--- /dev/null
+++ b/src/mesh/SubItemValuePerItemVariant.hpp
@@ -0,0 +1,163 @@
+#ifndef SUB_ITEM_VALUE_PER_ITEM_VARIANT_HPP
+#define SUB_ITEM_VALUE_PER_ITEM_VARIANT_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <utils/Exceptions.hpp>
+
+class SubItemValuePerItemVariant
+{
+ private:
+ using Variant = std::variant<NodeValuePerEdge<const bool>,
+ NodeValuePerEdge<const long int>,
+ NodeValuePerEdge<const unsigned long int>,
+ NodeValuePerEdge<const double>,
+ NodeValuePerEdge<const TinyVector<1, double>>,
+ NodeValuePerEdge<const TinyVector<2, double>>,
+ NodeValuePerEdge<const TinyVector<3, double>>,
+
+ NodeValuePerFace<const bool>,
+ NodeValuePerFace<const long int>,
+ NodeValuePerFace<const unsigned long int>,
+ NodeValuePerFace<const double>,
+ NodeValuePerFace<const TinyVector<1, double>>,
+ NodeValuePerFace<const TinyVector<2, double>>,
+ NodeValuePerFace<const TinyVector<3, double>>,
+
+ NodeValuePerCell<const bool>,
+ NodeValuePerCell<const long int>,
+ NodeValuePerCell<const unsigned long int>,
+ NodeValuePerCell<const double>,
+ NodeValuePerCell<const TinyVector<1, double>>,
+ NodeValuePerCell<const TinyVector<2, double>>,
+ NodeValuePerCell<const TinyVector<3, double>>,
+
+ EdgeValuePerNode<const bool>,
+ EdgeValuePerNode<const long int>,
+ EdgeValuePerNode<const unsigned long int>,
+ EdgeValuePerNode<const double>,
+ EdgeValuePerNode<const TinyVector<1, double>>,
+ EdgeValuePerNode<const TinyVector<2, double>>,
+ EdgeValuePerNode<const TinyVector<3, double>>,
+
+ EdgeValuePerFace<const bool>,
+ EdgeValuePerFace<const long int>,
+ EdgeValuePerFace<const unsigned long int>,
+ EdgeValuePerFace<const double>,
+ EdgeValuePerFace<const TinyVector<1, double>>,
+ EdgeValuePerFace<const TinyVector<2, double>>,
+ EdgeValuePerFace<const TinyVector<3, double>>,
+
+ EdgeValuePerCell<const bool>,
+ EdgeValuePerCell<const long int>,
+ EdgeValuePerCell<const unsigned long int>,
+ EdgeValuePerCell<const double>,
+ EdgeValuePerCell<const TinyVector<1, double>>,
+ EdgeValuePerCell<const TinyVector<2, double>>,
+ EdgeValuePerCell<const TinyVector<3, double>>,
+
+ FaceValuePerNode<const bool>,
+ FaceValuePerNode<const long int>,
+ FaceValuePerNode<const unsigned long int>,
+ FaceValuePerNode<const double>,
+ FaceValuePerNode<const TinyVector<1, double>>,
+ FaceValuePerNode<const TinyVector<2, double>>,
+ FaceValuePerNode<const TinyVector<3, double>>,
+
+ FaceValuePerEdge<const bool>,
+ FaceValuePerEdge<const long int>,
+ FaceValuePerEdge<const unsigned long int>,
+ FaceValuePerEdge<const double>,
+ FaceValuePerEdge<const TinyVector<1, double>>,
+ FaceValuePerEdge<const TinyVector<2, double>>,
+ FaceValuePerEdge<const TinyVector<3, double>>,
+
+ FaceValuePerCell<const bool>,
+ FaceValuePerCell<const long int>,
+ FaceValuePerCell<const unsigned long int>,
+ FaceValuePerCell<const double>,
+ FaceValuePerCell<const TinyVector<1, double>>,
+ FaceValuePerCell<const TinyVector<2, double>>,
+ FaceValuePerCell<const TinyVector<3, double>>,
+
+ CellValuePerNode<const bool>,
+ CellValuePerNode<const long int>,
+ CellValuePerNode<const unsigned long int>,
+ CellValuePerNode<const double>,
+ CellValuePerNode<const TinyVector<1, double>>,
+ CellValuePerNode<const TinyVector<2, double>>,
+ CellValuePerNode<const TinyVector<3, double>>,
+
+ CellValuePerEdge<const bool>,
+ CellValuePerEdge<const long int>,
+ CellValuePerEdge<const unsigned long int>,
+ CellValuePerEdge<const double>,
+ CellValuePerEdge<const TinyVector<1, double>>,
+ CellValuePerEdge<const TinyVector<2, double>>,
+ CellValuePerEdge<const TinyVector<3, double>>,
+
+ CellValuePerFace<const bool>,
+ CellValuePerFace<const long int>,
+ CellValuePerFace<const unsigned long int>,
+ CellValuePerFace<const double>,
+ CellValuePerFace<const TinyVector<1, double>>,
+ CellValuePerFace<const TinyVector<2, double>>,
+ CellValuePerFace<const TinyVector<3, double>>>;
+
+ Variant m_sub_item_value_per_item;
+
+ public:
+ PUGS_INLINE
+ const Variant&
+ itemValue() const
+ {
+ return m_sub_item_value_per_item;
+ }
+
+ template <typename SubItemValuePerItemT>
+ PUGS_INLINE auto
+ get() const
+ {
+ using DataType = typename SubItemValuePerItemT::data_type;
+ using ItemOfItemTypeT = typename SubItemValuePerItemT::ItemOfItemType;
+
+ if constexpr (std::is_same_v<SubItemValuePerItemT, SubItemValuePerItem<DataType, ItemOfItemTypeT>> or
+ std::is_same_v<SubItemValuePerItemT, SubItemValuePerItem<const DataType, ItemOfItemTypeT>> or
+ std::is_same_v<SubItemValuePerItemT, WeakSubItemValuePerItem<DataType, ItemOfItemTypeT>> or
+ std::is_same_v<SubItemValuePerItemT, WeakSubItemValuePerItem<const DataType, ItemOfItemTypeT>>) {
+ if (not std::holds_alternative<SubItemValuePerItem<const DataType, ItemOfItemTypeT>>(
+ this->m_sub_item_value_per_item)) {
+ throw NormalError("invalid SubItemValuePerItem type");
+ }
+ return std::get<SubItemValuePerItem<const DataType, ItemOfItemTypeT>>(this->m_sub_item_value_per_item);
+ } else {
+ static_assert(std::is_same_v<SubItemValuePerItemT, SubItemValuePerItemT>, "invalid template argument");
+ }
+ }
+
+ template <typename DataType, typename ItemOfItemTypeT>
+ SubItemValuePerItemVariant(const SubItemValuePerItem<DataType, ItemOfItemTypeT>& sub_item_value_per_item)
+ : m_sub_item_value_per_item{SubItemValuePerItem<const DataType, ItemOfItemTypeT>{sub_item_value_per_item}}
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, bool> or //
+ std::is_same_v<std::remove_const_t<DataType>, long int> or //
+ std::is_same_v<std::remove_const_t<DataType>, unsigned long int> or //
+ std::is_same_v<std::remove_const_t<DataType>, double> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<1, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<2, double>> or //
+ std::is_same_v<std::remove_const_t<DataType>, TinyVector<3, double>>,
+ "SubItemValuePerItem with this DataType is not allowed in variant");
+ }
+
+ SubItemValuePerItemVariant& operator=(SubItemValuePerItemVariant&&) = default;
+ SubItemValuePerItemVariant& operator=(const SubItemValuePerItemVariant&) = default;
+
+ SubItemValuePerItemVariant(const SubItemValuePerItemVariant&) = default;
+ SubItemValuePerItemVariant(SubItemValuePerItemVariant&&) = default;
+
+ SubItemValuePerItemVariant() = delete;
+ ~SubItemValuePerItemVariant() = default;
+};
+
+#endif // SUB_ITEM_VALUE_PER_ITEM_VARIANT_HPP
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index ae0ba69af..ef976be76 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -255,6 +255,34 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
this->_writePreamble<MeshType::Dimension>(output_named_item_data_set, fout);
}
+ for (const auto& [name, item_data_variant] : output_named_item_data_set) {
+ std::visit(
+ [&, name = name](auto&& item_data) {
+ using ItemDataType = std::decay_t<decltype(item_data)>;
+ if (ItemDataType::item_t == ItemType::face) {
+ std::ostringstream error_msg;
+ error_msg << "gnuplot_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) {
+ using ItemDataType = std::decay_t<decltype(item_data)>;
+ if (ItemDataType::item_t == ItemType::edge) {
+ std::ostringstream error_msg;
+ error_msg << "gnuplot_writer does not support edge data, cannot save variable \"" << rang::fgB::yellow << name
+ << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+ },
+ item_data_variant);
+ }
+
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
if (i_rank == parallel::rank()) {
std::ofstream fout(_getFilename(), std::ios_base::app);
@@ -268,13 +296,11 @@ GnuplotWriter::_write(const std::shared_ptr<const MeshType>& mesh,
}
void
-GnuplotWriter::_writeAtTime(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
- double time) const
+GnuplotWriter::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
@@ -312,12 +338,10 @@ GnuplotWriter::_writeMesh(const std::shared_ptr<const IMesh>& p_mesh) const
}
void
-GnuplotWriter::_write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+GnuplotWriter::_write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 112b46062..7d4c0c8a4 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -50,11 +50,12 @@ class GnuplotWriter final : public WriterBase
const OutputNamedItemDataSet& output_named_item_data_set,
std::optional<double> time) const;
- void _writeAtTime(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const final;
- void _write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const final;
+ void _write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 73e5a4e64..e98ebdd17 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -50,6 +50,20 @@ 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
@@ -269,6 +283,32 @@ GnuplotWriter1D::_write(const std::shared_ptr<const MeshType>& mesh,
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_1d_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);
@@ -311,13 +351,11 @@ GnuplotWriter1D::_writeMesh(const std::shared_ptr<const IMesh>&) const
}
void
-GnuplotWriter1D::_writeAtTime(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
- double time) const
+GnuplotWriter1D::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
@@ -338,12 +376,10 @@ GnuplotWriter1D::_writeAtTime(
}
void
-GnuplotWriter1D::_write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+GnuplotWriter1D::_write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index fad94c136..93f814046 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -22,6 +22,12 @@ class GnuplotWriter1D final : public WriterBase
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;
@@ -43,11 +49,12 @@ class GnuplotWriter1D final : public WriterBase
const OutputNamedItemDataSet& output_named_item_value_set,
std::optional<double> time) const;
- void _writeAtTime(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const final;
- void _write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const final;
+ void _write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
diff --git a/src/output/INamedDiscreteData.hpp b/src/output/INamedDiscreteData.hpp
new file mode 100644
index 000000000..5d22e82c2
--- /dev/null
+++ b/src/output/INamedDiscreteData.hpp
@@ -0,0 +1,27 @@
+#ifndef I_NAMED_DISCRETE_DATA_HPP
+#define I_NAMED_DISCRETE_DATA_HPP
+
+#include <string>
+
+class INamedDiscreteData
+{
+ public:
+ enum class Type
+ {
+ item_value,
+ discrete_function
+ };
+
+ virtual Type type() const = 0;
+
+ virtual const std::string& name() const = 0;
+
+ INamedDiscreteData(const INamedDiscreteData&) = default;
+ INamedDiscreteData(INamedDiscreteData&&) = default;
+
+ INamedDiscreteData() = default;
+
+ virtual ~INamedDiscreteData() = default;
+};
+
+#endif // I_NAMED_DISCRETE_DATA_HPP
diff --git a/src/output/IWriter.hpp b/src/output/IWriter.hpp
index dab629604..a277bdea3 100644
--- a/src/output/IWriter.hpp
+++ b/src/output/IWriter.hpp
@@ -1,7 +1,7 @@
#ifndef I_WRITER_HPP
#define I_WRITER_HPP
-#include <output/NamedDiscreteFunction.hpp>
+#include <output/INamedDiscreteData.hpp>
#include <memory>
#include <vector>
@@ -13,17 +13,27 @@ class IWriter
public:
virtual void writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
- virtual void write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const = 0;
+ virtual void write(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const = 0;
- virtual void writeIfNeeded(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
- double time) const = 0;
+ virtual void writeIfNeeded(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const = 0;
+
+ virtual void writeForced(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const = 0;
- virtual void writeForced(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ virtual void writeOnMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const = 0;
+
+ virtual void writeOnMeshIfNeeded(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const = 0;
+ virtual void writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const = 0;
+
IWriter() = default;
virtual ~IWriter() = default;
};
diff --git a/src/output/NamedDiscreteFunction.hpp b/src/output/NamedDiscreteFunction.hpp
index 1b0850040..d56d44369 100644
--- a/src/output/NamedDiscreteFunction.hpp
+++ b/src/output/NamedDiscreteFunction.hpp
@@ -1,20 +1,28 @@
#ifndef NAMED_DISCRETE_FUNCTION_HPP
#define NAMED_DISCRETE_FUNCTION_HPP
+#include <output/INamedDiscreteData.hpp>
+
#include <memory>
#include <string>
class IDiscreteFunction;
-class NamedDiscreteFunction
+class NamedDiscreteFunction final : public INamedDiscreteData
{
private:
std::shared_ptr<const IDiscreteFunction> m_discrete_function;
std::string m_name;
public:
+ Type
+ type() const final
+ {
+ return INamedDiscreteData::Type::discrete_function;
+ }
+
const std::string&
- name() const
+ name() const final
{
return m_name;
}
diff --git a/src/output/NamedItemValueVariant.hpp b/src/output/NamedItemValueVariant.hpp
new file mode 100644
index 000000000..0415582b7
--- /dev/null
+++ b/src/output/NamedItemValueVariant.hpp
@@ -0,0 +1,47 @@
+#ifndef NAMED_ITEM_VALUE_VARIANT_HPP
+#define NAMED_ITEM_VALUE_VARIANT_HPP
+
+#include <output/INamedDiscreteData.hpp>
+
+#include <memory>
+#include <string>
+
+class ItemValueVariant;
+
+class NamedItemValueVariant final : public INamedDiscreteData
+{
+ private:
+ std::shared_ptr<const ItemValueVariant> m_item_value_variant;
+ std::string m_name;
+
+ public:
+ Type
+ type() const final
+ {
+ return INamedDiscreteData::Type::item_value;
+ }
+
+ const std::string&
+ name() const final
+ {
+ return m_name;
+ }
+
+ const std::shared_ptr<const ItemValueVariant>
+ itemValueVariant() const
+ {
+ return m_item_value_variant;
+ }
+
+ NamedItemValueVariant(const std::shared_ptr<const ItemValueVariant>& item_value_variant, const std::string& name)
+ : m_item_value_variant{item_value_variant}, m_name{name}
+ {}
+
+ NamedItemValueVariant(const NamedItemValueVariant&) = default;
+ NamedItemValueVariant(NamedItemValueVariant&&) = default;
+
+ NamedItemValueVariant() = default;
+ ~NamedItemValueVariant() = default;
+};
+
+#endif // NAMED_ITEM_VALUE_VARIANT_HPP
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index 80c3d2052..aa2c5e2ed 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -72,6 +72,30 @@ class OutputNamedItemDataSet
NodeValue<const TinyMatrix<2, 2, double>>,
NodeValue<const TinyMatrix<3, 3, double>>,
+ EdgeValue<const bool>,
+ EdgeValue<const int>,
+ EdgeValue<const long int>,
+ EdgeValue<const unsigned long int>,
+ EdgeValue<const double>,
+ EdgeValue<const TinyVector<1, double>>,
+ EdgeValue<const TinyVector<2, double>>,
+ EdgeValue<const TinyVector<3, double>>,
+ EdgeValue<const TinyMatrix<1, 1, double>>,
+ EdgeValue<const TinyMatrix<2, 2, double>>,
+ EdgeValue<const TinyMatrix<3, 3, double>>,
+
+ FaceValue<const bool>,
+ FaceValue<const int>,
+ FaceValue<const long int>,
+ FaceValue<const unsigned long int>,
+ FaceValue<const double>,
+ FaceValue<const TinyVector<1, double>>,
+ FaceValue<const TinyVector<2, double>>,
+ FaceValue<const TinyVector<3, double>>,
+ FaceValue<const TinyMatrix<1, 1, double>>,
+ FaceValue<const TinyMatrix<2, 2, double>>,
+ FaceValue<const TinyMatrix<3, 3, double>>,
+
CellValue<const bool>,
CellValue<const int>,
CellValue<const long int>,
@@ -90,6 +114,18 @@ class OutputNamedItemDataSet
NodeArray<const unsigned long int>,
NodeArray<const double>,
+ EdgeArray<const bool>,
+ EdgeArray<const int>,
+ EdgeArray<const long int>,
+ EdgeArray<const unsigned long int>,
+ EdgeArray<const double>,
+
+ FaceArray<const bool>,
+ FaceArray<const int>,
+ FaceArray<const long int>,
+ FaceArray<const unsigned long int>,
+ FaceArray<const double>,
+
CellArray<const bool>,
CellArray<const int>,
CellArray<const long int>,
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 277c0d5cc..342a384e7 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -389,6 +389,26 @@ VTKWriter::_write(const std::shared_ptr<const MeshType>& mesh,
}
fout << "</PCells>\n";
+ for (const auto& [name, item_value_variant] : output_named_item_data_set) {
+ std::visit(
+ [&, name = name](auto&& item_value) {
+ using ItemValueType = std::decay_t<decltype(item_value)>;
+ if constexpr (ItemValueType::item_t == ItemType::edge) {
+ std::ostringstream error_msg;
+ error_msg << "VTK format does not support edge data, cannot save variable \"" << rang::fgB::yellow << name
+ << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+ if constexpr (ItemValueType::item_t == ItemType::face) {
+ std::ostringstream error_msg;
+ error_msg << "VTK format does not support face data, cannot save variable \"" << rang::fgB::yellow << name
+ << rang::fg::reset << '"';
+ throw NormalError(error_msg.str());
+ }
+ },
+ item_value_variant);
+ }
+
fout << "<PPointData>\n";
for (const auto& [name, item_value_variant] : output_named_item_data_set) {
std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
@@ -662,12 +682,11 @@ VTKWriter::_writeMesh(const std::shared_ptr<const IMesh>& mesh) const
}
void
-VTKWriter::_writeAtTime(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+VTKWriter::_writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
@@ -689,11 +708,10 @@ VTKWriter::_writeAtTime(const std::vector<std::shared_ptr<const NamedDiscreteFun
}
void
-VTKWriter::_write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+VTKWriter::_write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
- std::shared_ptr mesh = this->_getMesh(named_discrete_function_list);
-
- OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_function_list);
+ OutputNamedItemDataSet output_named_item_data_set = this->_getOutputNamedItemDataSet(named_discrete_data_list);
switch (mesh->dimension()) {
case 1: {
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index c0adbd6c4..1b6fc35f6 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -100,11 +100,12 @@ class VTKWriter final : public WriterBase
const OutputNamedItemDataSet& output_named_item_data_set,
std::optional<double> time) const;
- void _writeAtTime(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const final;
- void _write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const final;
+ void _write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 2fe733d54..0ce785cd9 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -1,6 +1,9 @@
#include <output/WriterBase.hpp>
#include <mesh/IMesh.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <output/NamedDiscreteFunction.hpp>
+#include <output/NamedItemValueVariant.hpp>
#include <output/OutputNamedItemValueSet.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
@@ -10,9 +13,9 @@
template <typename DiscreteFunctionType>
void
-WriterBase::_register(const std::string& name,
- const DiscreteFunctionType& discrete_function,
- OutputNamedItemDataSet& named_item_data_set)
+WriterBase::_registerDiscreteFunction(const std::string& name,
+ const DiscreteFunctionType& discrete_function,
+ OutputNamedItemDataSet& named_item_data_set)
{
if constexpr (DiscreteFunctionType::handled_data_type == IDiscreteFunction::HandledItemDataType::value) {
named_item_data_set.add(NamedItemData{name, discrete_function.cellValues()});
@@ -23,30 +26,30 @@ WriterBase::_register(const std::string& name,
template <size_t Dimension, template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
void
-WriterBase::_register(const std::string& name,
- const IDiscreteFunction& i_discrete_function,
- OutputNamedItemDataSet& named_item_data_set)
+WriterBase::_registerDiscreteFunction(const std::string& name,
+ const IDiscreteFunction& i_discrete_function,
+ OutputNamedItemDataSet& named_item_data_set)
{
const ASTNodeDataType& data_type = i_discrete_function.dataType();
switch (data_type) {
case ASTNodeDataType::bool_t: {
- _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, bool>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name, dynamic_cast<const DiscreteFunctionType<Dimension, bool>&>(i_discrete_function),
+ named_item_data_set);
break;
}
case ASTNodeDataType::unsigned_int_t: {
- _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, uint64_t>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name, dynamic_cast<const DiscreteFunctionType<Dimension, uint64_t>&>(i_discrete_function),
+ named_item_data_set);
break;
}
case ASTNodeDataType::int_t: {
- _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, int64_t>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name, dynamic_cast<const DiscreteFunctionType<Dimension, int64_t>&>(i_discrete_function),
+ named_item_data_set);
break;
}
case ASTNodeDataType::double_t: {
- _register(name, dynamic_cast<const DiscreteFunctionType<Dimension, double>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name, dynamic_cast<const DiscreteFunctionType<Dimension, double>&>(i_discrete_function),
+ named_item_data_set);
break;
}
case ASTNodeDataType::vector_t: {
@@ -56,21 +59,24 @@ WriterBase::_register(const std::string& name,
} else {
switch (data_type.dimension()) {
case 1: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<1, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<1, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
case 2: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<2, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<2, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
case 3: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<3, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyVector<3, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
default: {
@@ -88,21 +94,24 @@ WriterBase::_register(const std::string& name,
Assert(data_type.numberOfRows() == data_type.numberOfColumns(), "invalid matrix dimensions");
switch (data_type.numberOfRows()) {
case 1: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, 1, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, 1, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
case 2: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, 2, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, 2, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
case 3: {
- _register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, 3, double>>&>(i_discrete_function),
- named_item_data_set);
+ _registerDiscreteFunction(name,
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, 3, double>>&>(
+ i_discrete_function),
+ named_item_data_set);
break;
}
default: {
@@ -120,21 +129,22 @@ WriterBase::_register(const std::string& name,
template <template <size_t Dimension, typename DataType> typename DiscreteFunctionType>
void
-WriterBase::_register(const NamedDiscreteFunction& named_discrete_function, OutputNamedItemDataSet& named_item_data_set)
+WriterBase::_registerDiscreteFunction(const NamedDiscreteFunction& named_discrete_function,
+ OutputNamedItemDataSet& named_item_data_set)
{
const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
const std::string& name = named_discrete_function.name();
switch (i_discrete_function.mesh()->dimension()) {
case 1: {
- _register<1, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
+ _registerDiscreteFunction<1, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
break;
}
case 2: {
- _register<2, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
+ _registerDiscreteFunction<2, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
break;
}
case 3: {
- _register<3, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
+ _registerDiscreteFunction<3, DiscreteFunctionType>(name, i_discrete_function, named_item_data_set);
break;
}
default: {
@@ -143,54 +153,207 @@ WriterBase::_register(const NamedDiscreteFunction& named_discrete_function, Outp
}
}
+void
+WriterBase::_checkConnectivity(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
+{
+ Assert(named_discrete_data_list.size() > 0);
+
+ std::shared_ptr<const IConnectivity> connectivity = [&]() -> std::shared_ptr<const IConnectivity> {
+ switch (mesh->dimension()) {
+ case 1: {
+ return dynamic_cast<const Mesh<Connectivity<1>>&>(*mesh).shared_connectivity();
+ }
+ case 2: {
+ return dynamic_cast<const Mesh<Connectivity<2>>&>(*mesh).shared_connectivity();
+ }
+ case 3: {
+ return dynamic_cast<const Mesh<Connectivity<3>>&>(*mesh).shared_connectivity();
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
+ }();
+
+ for (size_t i = 0; i < named_discrete_data_list.size(); ++i) {
+ const auto& named_discrete_data = named_discrete_data_list[i];
+
+ if (named_discrete_data->type() == INamedDiscreteData::Type::item_value) {
+ const NamedItemValueVariant& named_item_value_variant =
+ dynamic_cast<const NamedItemValueVariant&>(*named_discrete_data);
+
+ std::visit(
+ [&](auto&& item_value) {
+ if (item_value.connectivity_ptr() != connectivity) {
+ std::ostringstream error_msg;
+ error_msg << "The variable " << rang::fgB::yellow << named_item_value_variant.name() << rang::fg::reset
+ << " is not defined on the provided same connectivity as the mesh\n";
+ throw NormalError(error_msg.str());
+ }
+ },
+ named_item_value_variant.itemValueVariant()->itemValue());
+ }
+ }
+}
+
+void
+WriterBase::_checkMesh(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
+{
+ Assert(named_discrete_data_list.size() > 0);
+
+ for (size_t i = 0; i < named_discrete_data_list.size(); ++i) {
+ const auto& named_discrete_data = named_discrete_data_list[i];
+
+ if (named_discrete_data->type() == INamedDiscreteData::Type::discrete_function) {
+ const NamedDiscreteFunction& named_discrete_function =
+ dynamic_cast<const NamedDiscreteFunction&>(*named_discrete_data);
+
+ if (mesh != named_discrete_function.discreteFunction()->mesh()) {
+ std::ostringstream error_msg;
+ error_msg << "The variable " << rang::fgB::yellow << named_discrete_function.name() << rang::fg::reset
+ << " is not defined on the provided mesh\n";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+}
+
std::shared_ptr<const IMesh>
-WriterBase::_getMesh(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+WriterBase::_getMesh(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
- Assert(named_discrete_function_list.size() > 0);
+ Assert(named_discrete_data_list.size() > 0);
+
+ std::map<std::shared_ptr<const IMesh>, std::string> mesh_set;
+ std::map<std::shared_ptr<const IConnectivity>, std::string> connectivity_set;
+
+ for (size_t i = 0; i < named_discrete_data_list.size(); ++i) {
+ const auto& named_discrete_data = named_discrete_data_list[i];
+
+ switch (named_discrete_data->type()) {
+ case INamedDiscreteData::Type::discrete_function: {
+ const NamedDiscreteFunction& named_discrete_function =
+ dynamic_cast<const NamedDiscreteFunction&>(*named_discrete_data);
- std::shared_ptr mesh = named_discrete_function_list[0]->discreteFunction()->mesh();
- for (size_t i = 1; i < named_discrete_function_list.size(); ++i) {
- if (mesh != named_discrete_function_list[i]->discreteFunction()->mesh()) {
+ std::shared_ptr mesh = named_discrete_function.discreteFunction()->mesh();
+ mesh_set[mesh] = named_discrete_function.name();
+
+ switch (mesh->dimension()) {
+ case 1: {
+ connectivity_set[dynamic_cast<const Mesh<Connectivity<1>>&>(*mesh).shared_connectivity()] =
+ named_discrete_function.name();
+ break;
+ }
+ case 2: {
+ connectivity_set[dynamic_cast<const Mesh<Connectivity<2>>&>(*mesh).shared_connectivity()] =
+ named_discrete_function.name();
+ break;
+ }
+ case 3: {
+ connectivity_set[dynamic_cast<const Mesh<Connectivity<3>>&>(*mesh).shared_connectivity()] =
+ named_discrete_function.name();
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
+ break;
+ }
+ case INamedDiscreteData::Type::item_value: {
+ const NamedItemValueVariant& named_item_value_variant =
+ dynamic_cast<const NamedItemValueVariant&>(*named_discrete_data);
+
+ std::visit([&](
+ auto&&
+ item_value) { connectivity_set[item_value.connectivity_ptr()] = named_item_value_variant.name(); },
+ named_item_value_variant.itemValueVariant()->itemValue());
+ }
+ }
+ }
+
+ if (mesh_set.size() != 1) {
+ if (mesh_set.size() == 0) {
+ throw NormalError("cannot find any mesh associated to output quantities");
+ } else {
std::ostringstream error_msg;
- error_msg << "discrete functions must be defined on the same mesh!\n"
- << rang::fgB::yellow << "note:" << rang::style::reset << " cannot write " << rang::fgB::blue
- << named_discrete_function_list[0]->name() << rang::style::reset << " and " << rang::fgB::blue
- << named_discrete_function_list[i]->name() << rang::style::reset << " in the same file.";
+ error_msg << "cannot save data using different " << rang::fgB::red << "meshes" << rang::fg::reset
+ << " in the same file!\n";
+ error_msg << rang::fgB::yellow << "note:" << rang::fg::reset
+ << "the following variables are defined on different meshes:";
+ for (const auto& [mesh, name] : mesh_set) {
+ error_msg << "\n- " << name;
+ }
throw NormalError(error_msg.str());
}
}
- return mesh;
+ if (connectivity_set.size() > 1) {
+ std::ostringstream error_msg;
+ error_msg << "cannot save data using different " << rang::fgB::red << "connectivities" << rang::fg::reset
+ << " in the same file!\n";
+ error_msg << rang::fgB::yellow << "note:" << rang::fg::reset
+ << "the following variables are defined on different connectivities:";
+ for (const auto& [connectivity, name] : connectivity_set) {
+ error_msg << "\n- " << name;
+ }
+ throw NormalError(error_msg.str());
+ }
+
+ return mesh_set.begin()->first;
}
OutputNamedItemDataSet
WriterBase::_getOutputNamedItemDataSet(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
OutputNamedItemDataSet named_item_data_set;
- for (auto& named_discrete_function : named_discrete_function_list) {
- const IDiscreteFunction& i_discrete_function = *named_discrete_function->discreteFunction();
+ for (auto& named_discrete_data : named_discrete_data_list) {
+ switch (named_discrete_data->type()) {
+ case INamedDiscreteData::Type::discrete_function: {
+ const NamedDiscreteFunction& named_discrete_function =
+ dynamic_cast<const NamedDiscreteFunction&>(*named_discrete_data);
- switch (i_discrete_function.descriptor().type()) {
- case DiscreteFunctionType::P0: {
- WriterBase::_register<DiscreteFunctionP0>(*named_discrete_function, named_item_data_set);
+ const IDiscreteFunction& i_discrete_function = *named_discrete_function.discreteFunction();
+
+ switch (i_discrete_function.descriptor().type()) {
+ case DiscreteFunctionType::P0: {
+ WriterBase::_registerDiscreteFunction<DiscreteFunctionP0>(named_discrete_function, named_item_data_set);
+ break;
+ }
+ case DiscreteFunctionType::P0Vector: {
+ WriterBase::_registerDiscreteFunction<DiscreteFunctionP0Vector>(named_discrete_function, named_item_data_set);
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << "the type of discrete function of " << rang::fgB::blue << named_discrete_data->name()
+ << rang::style::reset << " is not supported";
+ throw NormalError(error_msg.str());
+ }
+ }
break;
}
- case DiscreteFunctionType::P0Vector: {
- WriterBase::_register<DiscreteFunctionP0Vector>(*named_discrete_function, named_item_data_set);
+ case INamedDiscreteData::Type::item_value: {
+ const NamedItemValueVariant& named_item_value_variant =
+ dynamic_cast<const NamedItemValueVariant&>(*named_discrete_data);
+
+ const std::string& name = named_item_value_variant.name();
+
+ const ItemValueVariant& item_value_variant = *named_item_value_variant.itemValueVariant();
+
+ std::visit([&](auto&& item_value) { named_item_data_set.add(NamedItemData{name, item_value}); },
+ item_value_variant.itemValue());
break;
}
default: {
- std::ostringstream error_msg;
- error_msg << "the type of discrete function of " << rang::fgB::blue << named_discrete_function->name()
- << rang::style::reset << " is not supported";
- throw NormalError(error_msg.str());
+ throw UnexpectedError("invalid discrete data type");
}
}
}
-
return named_item_data_set;
}
@@ -201,7 +364,18 @@ WriterBase::WriterBase(const std::string& base_filename, const double& time_peri
WriterBase::WriterBase(const std::string& base_filename) : m_base_filename{base_filename} {}
void
-WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+WriterBase::write(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
+{
+ if (m_period_manager.has_value()) {
+ throw NormalError("this writer requires time value");
+ } else {
+ std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
+ this->_write(mesh, named_discrete_data_list);
+ }
+}
+
+void
+WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const
{
if (m_period_manager.has_value()) {
@@ -210,7 +384,8 @@ WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteF
return; // output already performed
if (time >= m_period_manager->nextTime()) {
- this->_writeAtTime(named_discrete_function_list, time);
+ std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
+ this->_writeAtTime(mesh, named_discrete_data_list, time);
m_period_manager->setSaveTime(time);
}
} else {
@@ -219,14 +394,14 @@ WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteF
}
void
-WriterBase::writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+WriterBase::writeForced(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const
{
if (m_period_manager.has_value()) {
if (time == m_period_manager->getLastTime())
return; // output already performed
-
- this->_writeAtTime(named_discrete_function_list, time);
+ std::shared_ptr<const IMesh> mesh = _getMesh(named_discrete_data_list);
+ this->_writeAtTime(mesh, named_discrete_data_list, time);
m_period_manager->setSaveTime(time);
} else {
throw NormalError("this writer does not allow time value");
@@ -234,12 +409,49 @@ WriterBase::writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFun
}
void
-WriterBase::write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const
+WriterBase::writeOnMesh(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const
{
if (m_period_manager.has_value()) {
throw NormalError("this writer requires time value");
} else {
- this->_write(named_discrete_function_list);
+ _checkMesh(mesh, named_discrete_data_list);
+ _checkConnectivity(mesh, named_discrete_data_list);
+ this->_write(mesh, named_discrete_data_list);
+ }
+}
+
+void
+WriterBase::writeOnMeshIfNeeded(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const
+{
+ if (m_period_manager.has_value()) {
+ if (time == m_period_manager->getLastTime())
+ return; // output already performed
+ _checkMesh(mesh, named_discrete_data_list);
+ _checkConnectivity(mesh, named_discrete_data_list);
+ this->_writeAtTime(mesh, named_discrete_data_list, time);
+ m_period_manager->setSaveTime(time);
+ } else {
+ throw NormalError("this writer does not allow time value");
+ }
+}
+
+void
+WriterBase::writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const
+{
+ if (m_period_manager.has_value()) {
+ if (time == m_period_manager->getLastTime())
+ return; // output already performed
+ _checkMesh(mesh, named_discrete_data_list);
+ _checkConnectivity(mesh, named_discrete_data_list);
+ this->_writeAtTime(mesh, named_discrete_data_list, time);
+ m_period_manager->setSaveTime(time);
+ } else {
+ throw NormalError("this writer does not allow time value");
}
}
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
index 2f29120b5..6f6a8e143 100644
--- a/src/output/WriterBase.hpp
+++ b/src/output/WriterBase.hpp
@@ -10,6 +10,8 @@
class IMesh;
class OutputNamedItemDataSet;
+class IDiscreteFunction;
+class NamedDiscreteFunction;
class WriterBase : public IWriter
{
@@ -86,39 +88,55 @@ class WriterBase : public IWriter
private:
template <typename DiscreteFunctionType>
- static void _register(const std::string& name, const DiscreteFunctionType&, OutputNamedItemDataSet&);
+ static void _registerDiscreteFunction(const std::string& name, const DiscreteFunctionType&, OutputNamedItemDataSet&);
template <size_t Dimension, template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
- static void _register(const std::string& name, const IDiscreteFunction&, OutputNamedItemDataSet&);
+ static void _registerDiscreteFunction(const std::string& name, const IDiscreteFunction&, OutputNamedItemDataSet&);
template <template <size_t DimensionT, typename DataTypeT> typename DiscreteFunctionType>
- static void _register(const NamedDiscreteFunction&, OutputNamedItemDataSet&);
+ static void _registerDiscreteFunction(const NamedDiscreteFunction&, OutputNamedItemDataSet&);
protected:
+ void _checkConnectivity(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const;
+
+ void _checkMesh(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const;
+
std::shared_ptr<const IMesh> _getMesh(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const;
OutputNamedItemDataSet _getOutputNamedItemDataSet(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const;
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const;
- virtual void _writeAtTime(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
- double time) const = 0;
+ virtual void _writeAtTime(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const = 0;
- virtual void _write(
- const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const = 0;
+ virtual void _write(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const = 0;
virtual void _writeMesh(const std::shared_ptr<const IMesh>& mesh) const = 0;
public:
- void write(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list) const final;
+ void write(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
- void writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ void writeIfNeeded(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const final;
- void writeForced(const std::vector<std::shared_ptr<const NamedDiscreteFunction>>& named_discrete_function_list,
+ void writeForced(const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
double time) const final;
+ void writeOnMesh(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list) const final;
+
+ void writeOnMeshIfNeeded(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const final;
+ void writeOnMeshForced(const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const INamedDiscreteData>>& named_discrete_data_list,
+ double time) const final;
+
void writeMesh(const std::shared_ptr<const IMesh>& mesh) const final;
WriterBase() = delete;
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 1ef88c67a..589b1d380 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -2,9 +2,11 @@
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -87,14 +89,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
using BoundaryConditionList = std::vector<BoundaryCondition>;
- const SolverType m_solver_type;
- BoundaryConditionList m_boundary_condition_list;
-
- NodeValue<const Rd> m_ur;
- NodeValuePerCell<const Rd> m_Fjr;
-
NodeValuePerCell<const Rdxd>
- _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
+ _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -116,7 +112,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
NodeValuePerCell<const Rdxd>
- _computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
+ _computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -164,12 +160,12 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
NodeValuePerCell<const Rdxd>
- _computeAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
+ _computeAjr(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
{
if constexpr (Dimension == 1) {
return _computeGlaceAjr(mesh, rhoc);
} else {
- switch (m_solver_type) {
+ switch (solver_type) {
case SolverType::Glace: {
return _computeGlaceAjr(mesh, rhoc);
}
@@ -184,7 +180,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
NodeValue<Rdxd>
- _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr)
+ _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
{
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
@@ -212,7 +208,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
_computeBr(const Mesh<Connectivity<Dimension>>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const DiscreteVectorFunction& u,
- const DiscreteScalarFunction& p)
+ const DiscreteScalarFunction& p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -243,7 +239,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
BoundaryConditionList
_getBCList(const MeshType& mesh,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
BoundaryConditionList bc_list;
@@ -327,21 +323,27 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
return bc_list;
}
- void _applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br);
- void _applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
- void _applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
- void _applyExternalVelocityBC(const DiscreteScalarFunction& p, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+ void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
+ void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+ void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+ void _applyExternalVelocityBC(const BoundaryConditionList& bc_list,
+ const DiscreteScalarFunction& p,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const;
void
- _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+ _applyBoundaryConditions(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
- this->_applyPressureBC(mesh, br);
- this->_applySymmetryBC(Ar, br);
- this->_applyVelocityBC(Ar, br);
+ this->_applyPressureBC(bc_list, mesh, br);
+ this->_applySymmetryBC(bc_list, Ar, br);
+ this->_applyVelocityBC(bc_list, Ar, br);
}
NodeValue<const Rd>
- _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br)
+ _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
NodeValue<Rd> u{mesh.connectivity()};
parallel_for(
@@ -355,7 +357,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
const NodeValuePerCell<const Rdxd>& Ajr,
const NodeValue<const Rd>& ur,
const DiscreteVectorFunction& u,
- const DiscreteScalarFunction& p)
+ const DiscreteScalarFunction& p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -376,46 +378,71 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
return F;
}
- AcousticSolver(SolverType solver_type,
- const MeshType& mesh,
- const DiscreteScalarFunction& rho,
- const DiscreteScalarFunction& c,
- const DiscreteVectorFunction& u,
- const DiscreteScalarFunction& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
- : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(mesh, bc_descriptor_list)}
+ public:
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const IDiscreteFunction>& i_rho,
+ const std::shared_ptr<const IDiscreteFunction>& i_c,
+ const std::shared_ptr<const IDiscreteFunction>& i_u,
+ const std::shared_ptr<const IDiscreteFunction>& i_p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
- NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rho * c);
+ std::shared_ptr i_mesh = getCommonMesh({i_rho, i_c, i_u, i_p});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({i_rho, i_c, i_u, i_p}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+ const DiscreteScalarFunction& rho = dynamic_cast<const DiscreteScalarFunction&>(*i_rho);
+ const DiscreteScalarFunction& c = dynamic_cast<const DiscreteScalarFunction&>(*i_c);
+ const DiscreteVectorFunction& u = dynamic_cast<const DiscreteVectorFunction&>(*i_u);
+ const DiscreteScalarFunction& p = dynamic_cast<const DiscreteScalarFunction&>(*i_p);
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
- this->_applyBoundaryConditions(mesh, Ar, br);
- this->_applyExternalVelocityBC(p, Ar, br);
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+ this->_applyExternalVelocityBC(bc_list, p, Ar, br);
synchronize(Ar);
synchronize(br);
- m_ur = this->_computeUr(mesh, Ar, br);
- m_Fjr = this->_computeFjr(mesh, Ajr, m_ur, u, p);
+ NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
- public:
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
- apply(const double& dt,
- const MeshType& mesh,
- const DiscreteFunctionP0<Dimension, double>& rho,
- const DiscreteFunctionP0<Dimension, Rd>& u,
- const DiscreteFunctionP0<Dimension, double>& E) const
+ apply_fluxes(const double& dt,
+ const MeshType& mesh,
+ const DiscreteFunctionP0<Dimension, double>& rho,
+ const DiscreteFunctionP0<Dimension, Rd>& u,
+ const DiscreteFunctionP0<Dimension, double>& E,
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Fjr) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+ throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+ }
+
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
- mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
@@ -433,8 +460,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
double energy_fluxes = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
- momentum_fluxes += m_Fjr(j, R);
- energy_fluxes += dot(m_Fjr(j, R), m_ur[r]);
+ momentum_fluxes += Fjr(j, R);
+ energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
new_u[j] -= dt_over_Mj * momentum_fluxes;
@@ -455,10 +482,12 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>>
- apply(const double& dt,
- const std::shared_ptr<const IDiscreteFunction>& rho,
- const std::shared_ptr<const IDiscreteFunction>& u,
- const std::shared_ptr<const IDiscreteFunction>& E) const
+ apply_fluxes(const double& dt,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
{
std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
if (not i_mesh) {
@@ -469,26 +498,31 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
throw NormalError("acoustic solver expects P0 functions");
}
- return this->apply(dt, dynamic_cast<const MeshType&>(*i_mesh), dynamic_cast<const DiscreteScalarFunction&>(*rho),
- dynamic_cast<const DiscreteVectorFunction&>(*u),
- dynamic_cast<const DiscreteScalarFunction&>(*E));
- }
-
- AcousticSolver(SolverType solver_type,
- const std::shared_ptr<const IMesh>& mesh,
- const std::shared_ptr<const IDiscreteFunction>& rho,
- const std::shared_ptr<const IDiscreteFunction>& c,
- const std::shared_ptr<const IDiscreteFunction>& u,
- const std::shared_ptr<const IDiscreteFunction>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
- : AcousticSolver{solver_type,
- dynamic_cast<const MeshType&>(*mesh),
- dynamic_cast<const DiscreteScalarFunction&>(*rho),
- dynamic_cast<const DiscreteScalarFunction&>(*c),
- dynamic_cast<const DiscreteVectorFunction&>(*u),
- dynamic_cast<const DiscreteScalarFunction&>(*p),
- bc_descriptor_list}
- {}
+ return this->apply_fluxes(dt, //
+ dynamic_cast<const MeshType&>(*i_mesh), //
+ dynamic_cast<const DiscreteScalarFunction&>(*rho), //
+ dynamic_cast<const DiscreteVectorFunction&>(*u), //
+ dynamic_cast<const DiscreteScalarFunction&>(*E), //
+ ur->get<NodeValue<const Rd>>(), //
+ Fjr->get<NodeValuePerCell<const Rd>>());
+ }
+
+ std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>>
+ apply(const SolverType& solver_type,
+ const double& dt,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
+ return apply_fluxes(dt, rho, u, E, ur, Fjr);
+ }
AcousticSolver() = default;
AcousticSolver(AcousticSolver&&) = default;
@@ -497,9 +531,11 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
template <size_t Dimension>
void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br)
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& br) const
{
- for (const auto& boundary_condition : m_boundary_condition_list) {
+ for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
@@ -560,9 +596,11 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshTyp
template <size_t Dimension>
void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
- for (const auto& boundary_condition : m_boundary_condition_list) {
+ for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
@@ -589,9 +627,11 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(NodeValue<Rdx
template <size_t Dimension>
void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
- for (const auto& boundary_condition : m_boundary_condition_list) {
+ for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
@@ -624,11 +664,12 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdx
template <size_t Dimension>
void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const DiscreteScalarFunction& p,
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
+ const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
- NodeValue<Rd>& br)
+ NodeValue<Rd>& br) const
{
- for (const auto& boundary_condition : m_boundary_condition_list) {
+ for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
@@ -840,34 +881,23 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryConditio
~SymmetryBoundaryCondition() = default;
};
-AcousticSolverHandler::AcousticSolverHandler(
- SolverType solver_type,
- const std::shared_ptr<const IDiscreteFunction>& rho,
- const std::shared_ptr<const IDiscreteFunction>& c,
- const std::shared_ptr<const IDiscreteFunction>& u,
- const std::shared_ptr<const IDiscreteFunction>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
{
- std::shared_ptr i_mesh = getCommonMesh({rho, c, u, p});
if (not i_mesh) {
throw NormalError("discrete functions are not defined on the same mesh");
}
- if (not checkDiscretizationType({rho, c, u, p}, DiscreteFunctionType::P0)) {
- throw NormalError("acoustic solver expects P0 functions");
- }
-
switch (i_mesh->dimension()) {
case 1: {
- m_acoustic_solver = std::make_unique<AcousticSolver<1>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
+ m_acoustic_solver = std::make_unique<AcousticSolver<1>>();
break;
}
case 2: {
- m_acoustic_solver = std::make_unique<AcousticSolver<2>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
+ m_acoustic_solver = std::make_unique<AcousticSolver<2>>();
break;
}
case 3: {
- m_acoustic_solver = std::make_unique<AcousticSolver<3>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
+ m_acoustic_solver = std::make_unique<AcousticSolver<3>>();
break;
}
default: {
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 2ed970757..489b8e795 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -8,6 +8,8 @@
class IDiscreteFunction;
class IBoundaryConditionDescriptor;
class IMesh;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
@@ -23,14 +25,39 @@ class AcousticSolverHandler
private:
struct IAcousticSolver
{
+ virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(
+ const SolverType& solver_type,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+ virtual std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>>
+ apply_fluxes(const double& dt,
+ const std::shared_ptr<const IDiscreteFunction>& rho,
+ const std::shared_ptr<const IDiscreteFunction>& u,
+ const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
+
virtual std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>,
std::shared_ptr<const IDiscreteFunction>>
- apply(const double& dt,
+ apply(const SolverType& solver_type,
+ const double& dt,
const std::shared_ptr<const IDiscreteFunction>& rho,
const std::shared_ptr<const IDiscreteFunction>& u,
- const std::shared_ptr<const IDiscreteFunction>& E) const = 0;
+ const std::shared_ptr<const IDiscreteFunction>& E,
+ const std::shared_ptr<const IDiscreteFunction>& c,
+ const std::shared_ptr<const IDiscreteFunction>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
IAcousticSolver() = default;
IAcousticSolver(IAcousticSolver&&) = default;
@@ -51,12 +78,7 @@ class AcousticSolverHandler
return *m_acoustic_solver;
}
- AcousticSolverHandler(SolverType solver_type,
- const std::shared_ptr<const IDiscreteFunction>& rho,
- const std::shared_ptr<const IDiscreteFunction>& c,
- const std::shared_ptr<const IDiscreteFunction>& u,
- const std::shared_ptr<const IDiscreteFunction>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+ AcousticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
};
#endif // ACOUSTIC_SOLVER_HPP
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index daf7646ee..e19073f1e 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -185,7 +185,6 @@ DiscreteFunctionInterpoler::_interpolate() const
std::shared_ptr<IDiscreteFunction>
DiscreteFunctionInterpoler::interpolate() const
{
- std::shared_ptr<IDiscreteFunction> discrete_function;
switch (m_mesh->dimension()) {
case 1: {
return this->_interpolate<1>();
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 196938400..3f1743981 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -175,6 +175,8 @@ add_executable (mpi_unit_tests
test_ItemArrayUtils.cpp
test_ItemValue.cpp
test_ItemValueUtils.cpp
+ test_ItemValueVariant.cpp
+ test_ItemValueVariantFunctionInterpoler.cpp
test_MeshEdgeBoundary.cpp
test_MeshFaceBoundary.cpp
test_MeshFlatEdgeBoundary.cpp
@@ -189,6 +191,7 @@ add_executable (mpi_unit_tests
test_Partitioner.cpp
test_RandomEngine.cpp
test_SubItemValuePerItem.cpp
+ test_SubItemValuePerItemVariant.cpp
test_SubItemValuePerItemUtils.cpp
test_SubItemArrayPerItem.cpp
test_SubItemArrayPerItemUtils.cpp
diff --git a/tests/test_ItemValueVariant.cpp b/tests/test_ItemValueVariant.cpp
new file mode 100644
index 000000000..4d8fca46d
--- /dev/null
+++ b/tests/test_ItemValueVariant.cpp
@@ -0,0 +1,101 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("ItemValueVariant", "[mesh]")
+{
+ std::shared_ptr mesh = MeshDataBaseForTests::get().hybrid2DMesh();
+
+ const Connectivity<2>& connectivity = *mesh->shared_connectivity();
+
+ using R1 = TinyVector<1>;
+ using R2 = TinyVector<2>;
+ using R3 = TinyVector<3>;
+
+ using R1x1 = TinyMatrix<1>;
+ using R2x2 = TinyMatrix<2>;
+ using R3x3 = TinyMatrix<3>;
+
+ SECTION("NodeValue<double>")
+ {
+ NodeValue<double> node_value{connectivity};
+ ItemValueVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<NodeValue<const double>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValue<int64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<uint64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R1x1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R2x2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R3x3>>(), "error: invalid ItemValue type");
+
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const double>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const double>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const double>>(), "error: invalid ItemValue type");
+ }
+
+ SECTION("EdgeValue<R3>")
+ {
+ EdgeValue<R3> node_value{connectivity};
+ ItemValueVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<int64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<uint64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<double>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R2>>(), "error: invalid ItemValue type");
+ REQUIRE_NOTHROW(v.get<EdgeValue<const R3>>());
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R1x1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R2x2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R3x3>>(), "error: invalid ItemValue type");
+
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R3>>(), "error: invalid ItemValue type");
+ }
+
+ SECTION("FaceValue<R3x3>")
+ {
+ FaceValue<R3x3> node_value{connectivity};
+ ItemValueVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<FaceValue<int64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<uint64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<double>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R1x1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const R2x2>>(), "error: invalid ItemValue type");
+ REQUIRE_NOTHROW(v.get<FaceValue<const R3x3>>());
+
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const R3x3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const R3x3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R3x3>>(), "error: invalid ItemValue type");
+ }
+
+ SECTION("CellValue<int64_t>")
+ {
+ CellValue<int64_t> node_value{connectivity};
+ ItemValueVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<CellValue<const int64_t>>());
+ REQUIRE_THROWS_WITH(v.get<CellValue<const uint64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const double>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R3>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R1x1>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R2x2>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<CellValue<const R3x3>>(), "error: invalid ItemValue type");
+
+ REQUIRE_THROWS_WITH(v.get<NodeValue<const int64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValue<const int64_t>>(), "error: invalid ItemValue type");
+ REQUIRE_THROWS_WITH(v.get<FaceValue<const int64_t>>(), "error: invalid ItemValue type");
+ }
+}
diff --git a/tests/test_ItemValueVariantFunctionInterpoler.cpp b/tests/test_ItemValueVariantFunctionInterpoler.cpp
new file mode 100644
index 000000000..8689b1086
--- /dev/null
+++ b/tests/test_ItemValueVariantFunctionInterpoler.cpp
@@ -0,0 +1,547 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <language/ast/ASTBuilder.hpp>
+#include <language/ast/ASTModulesImporter.hpp>
+#include <language/ast/ASTNodeDataTypeBuilder.hpp>
+#include <language/ast/ASTNodeExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionEvaluationExpressionBuilder.hpp>
+#include <language/ast/ASTNodeFunctionExpressionBuilder.hpp>
+#include <language/ast/ASTNodeTypeCleaner.hpp>
+#include <language/ast/ASTSymbolTableBuilder.hpp>
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/SymbolTable.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <language/utils/ItemValueVariantFunctionInterpoler.hpp>
+
+#include <pegtl/string_input.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("ItemValueVariantFunctionInterpoler", "[scheme]")
+{
+ auto same_item_value = [](auto f, auto g) -> bool {
+ using ItemIdType = typename decltype(f)::index_type;
+ for (ItemIdType item_id = 0; item_id < f.numberOfItems(); ++item_id) {
+ if (f[item_id] != g[item_id]) {
+ return false;
+ }
+ }
+
+ return true;
+ };
+
+ SECTION("1D")
+ {
+ constexpr size_t Dimension = 1;
+
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ auto xj = MeshDataManager::instance().getMeshData(*mesh_1d).xj();
+ auto xr = mesh_1d->xr();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_1d: R^1 -> B, x -> (exp(2 * x[0]) + 3 > 4);
+let N_scalar_non_linear_1d: R^1 -> N, x -> floor(3 * x[0] * x[0] + 2);
+let Z_scalar_non_linear_1d: R^1 -> Z, x -> floor(exp(2 * x[0]) - 1);
+let R_scalar_non_linear_1d: R^1 -> R, x -> 2 * exp(x[0]) + 3;
+let R1_non_linear_1d: R^1 -> R^1, x -> 2 * exp(x[0]);
+let R2_non_linear_1d: R^1 -> R^2, x -> [2 * exp(x[0]), -3*x[0]];
+let R3_non_linear_1d: R^1 -> R^3, x -> [2 * exp(x[0]) + 3, x[0] - 2, 3];
+let R1x1_non_linear_1d: R^1 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[0]) + 3);
+let R2x2_non_linear_1d: R^1 -> R^2x2, x -> [[2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0])], [3, x[0] * x[0]]];
+let R3x3_non_linear_1d: R^1 -> R^3x3, x -> [[2 * exp(x[0]) * sin(x[0]) + 3, sin(x[0] - 2 * x[0]), 3], [x[0] * x[0], -4*x[0], 2*x[0]+1], [3, -6*x[0], exp(x[0])]];
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("B_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<bool> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::exp(2 * x[0]) + 3 > 4;
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_data_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_data_variant->get<CellValue<bool>>()));
+ }
+
+ SECTION("N_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ NodeValue<uint64_t> node_value{mesh_1d->connectivity()};
+ parallel_for(
+ node_value.numberOfItems(), PUGS_LAMBDA(const NodeId node_id) {
+ const TinyVector<Dimension>& x = xr[node_id];
+ node_value[node_id] = std::floor(3 * x[0] * x[0] + 2);
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::node, function_symbol_id);
+ std::shared_ptr item_data_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(node_value, item_data_variant->get<NodeValue<uint64_t>>()));
+ }
+
+ SECTION("Z_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<int64_t> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = std::floor(std::exp(2 * x[0]) - 1);
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<int64_t>>()));
+ }
+
+ SECTION("R_scalar_non_linear_1d")
+ {
+ auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<double> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = 2 * std::exp(x[0]) + 3;
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<double>>()));
+ }
+
+ SECTION("R1_non_linear_1d")
+ {
+ using DataType = TinyVector<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0])};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+
+ SECTION("R2_non_linear_1d")
+ {
+ using DataType = TinyVector<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]), -3 * x[0]};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+
+ SECTION("R3_non_linear_1d")
+ {
+ using DataType = TinyVector<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) + 3, x[0] - 2, 3};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+
+ SECTION("R1x1_non_linear_1d")
+ {
+ using DataType = TinyMatrix<1>;
+
+ auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+
+ SECTION("R2x2_non_linear_1d")
+ {
+ using DataType = TinyMatrix<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] =
+ DataType{2 * std::exp(x[0]) * std::sin(x[0]) + 3, std::sin(x[0] - 2 * x[0]), 3, x[0] * x[0]};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+
+ SECTION("R3x3_non_linear_1d")
+ {
+ using DataType = TinyMatrix<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_1d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ CellValue<DataType> cell_value{mesh_1d->connectivity()};
+ parallel_for(
+ cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension>& x = xj[cell_id];
+ cell_value[cell_id] = DataType{2 * exp(x[0]) * std::sin(x[0]) + 3,
+ std::sin(x[0] - 2 * x[0]),
+ 3,
+ x[0] * x[0],
+ -4 * x[0],
+ 2 * x[0] + 1,
+ 3,
+ -6 * x[0],
+ std::exp(x[0])};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_1d, ItemType::cell, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(cell_value, item_value_variant->get<CellValue<DataType>>()));
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ constexpr size_t Dimension = 2;
+
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ auto xl = MeshDataManager::instance().getMeshData(*mesh_2d).xl();
+
+ std::string_view data = R"(
+import math;
+let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
+let R2_non_linear_2d: R^2 -> R^2, x -> [2 * exp(x[0]), -3*x[1]];
+let R3x3_non_linear_2d: R^2 -> R^3x3, x -> [[2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3], [x[1] * x[0], -4*x[1], 2*x[0]+1], [3, -6*x[0], exp(x[1])]];
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("B_scalar_non_linear_2d")
+ {
+ using DataType = bool;
+
+ auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ FaceValue<DataType> face_value{mesh_2d->connectivity()};
+ parallel_for(
+ face_value.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+ const TinyVector<Dimension>& x = xl[face_id];
+ face_value[face_id] = std::exp(2 * x[0]) < 2 * x[1];
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(face_value, item_value_variant->get<FaceValue<DataType>>()));
+ }
+
+ SECTION("R2_non_linear_2d")
+ {
+ using DataType = TinyVector<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ FaceValue<DataType> face_value{mesh_2d->connectivity()};
+ parallel_for(
+ face_value.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+ const TinyVector<Dimension>& x = xl[face_id];
+ face_value[face_id] = DataType{2 * std::exp(x[0]), -3 * x[1]};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(face_value, item_value_variant->get<FaceValue<DataType>>()));
+ }
+
+ SECTION("R3x3_non_linear_2d")
+ {
+ using DataType = TinyMatrix<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ FaceValue<DataType> face_value{mesh_2d->connectivity()};
+ parallel_for(
+ face_value.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) {
+ const TinyVector<Dimension>& x = xl[face_id];
+
+ face_value[face_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+ std::sin(x[1] - 2 * x[0]),
+ 3,
+ x[1] * x[0],
+ -4 * x[1],
+ 2 * x[0] + 1,
+ 3,
+ -6 * x[0],
+ std::exp(x[1])};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_2d, ItemType::face, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(face_value, item_value_variant->get<FaceValue<DataType>>()));
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ constexpr size_t Dimension = 3;
+
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ auto xe = MeshDataManager::instance().getMeshData(*mesh_3d).xe();
+
+ std::string_view data = R"(
+import math;
+let R_scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]+x[2]) + 3 * x[1];
+let R3_non_linear_3d: R^3 -> R^3, x -> [2 * exp(x[0]) + 3, x[1] - 2, 3 * x[2]];
+let R2x2_non_linear_3d: R^3 -> R^2x2, x -> [[2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0])], [3, x[1] * x[0] - x[2]]];
+)";
+ TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+ auto ast = ASTBuilder::build(input);
+
+ ASTModulesImporter{*ast};
+ ASTNodeTypeCleaner<language::import_instruction>{*ast};
+
+ ASTSymbolTableBuilder{*ast};
+ ASTNodeDataTypeBuilder{*ast};
+
+ ASTNodeTypeCleaner<language::var_declaration>{*ast};
+ ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+ ASTNodeExpressionBuilder{*ast};
+
+ std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+ TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+ position.byte = data.size(); // ensure that variables are declared at this point
+
+ SECTION("R_scalar_non_linear_3d")
+ {
+ auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ EdgeValue<double> edge_value{mesh_3d->connectivity()};
+ parallel_for(
+ edge_value.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+ const TinyVector<Dimension>& x = xe[edge_id];
+ edge_value[edge_id] = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(edge_value, item_value_variant->get<EdgeValue<double>>()));
+ }
+
+ SECTION("R3_non_linear_3d")
+ {
+ using DataType = TinyVector<3>;
+
+ auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ EdgeValue<DataType> edge_value{mesh_3d->connectivity()};
+ parallel_for(
+ edge_value.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+ const TinyVector<Dimension>& x = xe[edge_id];
+ edge_value[edge_id] = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(edge_value, item_value_variant->get<EdgeValue<DataType>>()));
+ }
+
+ SECTION("R2x2_non_linear_3d")
+ {
+ using DataType = TinyMatrix<2>;
+
+ auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
+ REQUIRE(found);
+ REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+ FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+ EdgeValue<DataType> edge_value{mesh_3d->connectivity()};
+ parallel_for(
+ edge_value.numberOfItems(), PUGS_LAMBDA(const EdgeId edge_id) {
+ const TinyVector<Dimension>& x = xe[edge_id];
+ edge_value[edge_id] =
+ DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
+ });
+
+ ItemValueVariantFunctionInterpoler interpoler(mesh_3d, ItemType::edge, function_symbol_id);
+ std::shared_ptr item_value_variant = interpoler.interpolate();
+
+ REQUIRE(same_item_value(edge_value, item_value_variant->get<EdgeValue<DataType>>()));
+ }
+ }
+ }
+ }
+}
diff --git a/tests/test_SubItemValuePerItemVariant.cpp b/tests/test_SubItemValuePerItemVariant.cpp
new file mode 100644
index 000000000..d37e9ab40
--- /dev/null
+++ b/tests/test_SubItemValuePerItemVariant.cpp
@@ -0,0 +1,165 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SubItemValuePerItemVariant", "[mesh]")
+{
+ std::shared_ptr mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+
+ const Connectivity<3>& connectivity = *mesh->shared_connectivity();
+
+ using R1 = TinyVector<1>;
+ using R2 = TinyVector<2>;
+ using R3 = TinyVector<3>;
+
+ SECTION("NodeValuePerCell<double>")
+ {
+ NodeValuePerCell<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<NodeValuePerCell<double>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("NodeValuePerFace<R1>")
+ {
+ NodeValuePerFace<R1> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<const double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<const int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<const uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<NodeValuePerFace<const R1>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<const R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<const R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("NodeValuePerEdge<int64_t>")
+ {
+ NodeValuePerEdge<int64_t> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<NodeValuePerEdge<int64_t>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("EdgeValuePerCell<R2>")
+ {
+ EdgeValuePerCell<R2> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerCell<double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerCell<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerCell<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerCell<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<EdgeValuePerCell<R2>>());
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerCell<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("EdgeValuePerFace<R1>")
+ {
+ EdgeValuePerFace<R1> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerFace<const double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerFace<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerFace<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<EdgeValuePerFace<R1>>());
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerFace<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerFace<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("EdgeValuePerNode<double>")
+ {
+ EdgeValuePerNode<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<EdgeValuePerNode<double>>());
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerNode<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerNode<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerNode<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerNode<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<EdgeValuePerNode<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("FaceValuePerCell<uint64_t>")
+ {
+ FaceValuePerCell<uint64_t> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerCell<double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerCell<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<FaceValuePerCell<uint64_t>>());
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerCell<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerCell<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerCell<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("FaceValuePerEdge<double>")
+ {
+ FaceValuePerEdge<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<FaceValuePerEdge<double>>());
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerEdge<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerEdge<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerEdge<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerEdge<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerEdge<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("FaceValuePerNode<double>")
+ {
+ FaceValuePerNode<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<FaceValuePerNode<double>>());
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerNode<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerNode<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerNode<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerNode<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<FaceValuePerNode<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("NodeValuePerCell<double>")
+ {
+ NodeValuePerCell<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<NodeValuePerCell<double>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerCell<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+
+ SECTION("NodeValuePerFace<R3>")
+ {
+ NodeValuePerFace<R3> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<double>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerFace<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_NOTHROW(v.get<NodeValuePerFace<R3>>());
+ }
+
+ SECTION("NodeValuePerEdge<double>")
+ {
+ NodeValuePerEdge<double> node_value{connectivity};
+ SubItemValuePerItemVariant v(node_value);
+ REQUIRE_NOTHROW(v.get<NodeValuePerEdge<double>>());
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<int64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<uint64_t>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R1>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R2>>(), "error: invalid SubItemValuePerItem type");
+ REQUIRE_THROWS_WITH(v.get<NodeValuePerEdge<R3>>(), "error: invalid SubItemValuePerItem type");
+ }
+}
--
GitLab