diff --git a/doc/userdoc.org b/doc/userdoc.org
index 12ce427db3bdd792588d6e32a12c9b0af980c02e..a91d0d5b850289bb953ab1d172d03bbbaea3e919 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 4fc75f4f2b61058b339c748f822720107c92dd1a..9c5cf3aec1e11d16f652357db53f68db54f6ed79 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 9047412586f24200789ebea289e9b8c4f0fbd27d..bf9f0b714b327a1216f7ec3e54fc98a3ec65d70f 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 01dcc66eb23da733c9b5d5922d8bdb55156d73ca..e3c6f1b96bf5ae660f7e3648ab70858edd4c20a3 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 89d25001b6f70ceaec5c8465718cbd01da32a928..6523d56a93d4fe0b4f953fb6065cd4a15be7330c 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 adcd92530298a99f088743e7ae6583bc09faa93c..f61eafb0e4847c248e9def011bf269b3b9a4022e 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 b76655995fd813239304c75c8189a4836c20e6d5..0fe1c798ba77d98b4453ded4d4c3334c0d2199e5 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 0000000000000000000000000000000000000000..018dd7a9837997e4faf595898773f7c847f47f20
--- /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 0000000000000000000000000000000000000000..4e03f72fc41ec74cf628e4ced7d331b6354d6f6a
--- /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 496c33242317e49f19842441e9c0fbec0f16ef70..da6bb88eb5a9d92dda08ee7ed535aae104639d95 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 0000000000000000000000000000000000000000..1bfff8f84248f29b0cab089e8e246a8d59c369cc
--- /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 ada0f8fbdbcf0b8d2075ebbb8ad05cafa3a63a2e..9b0c89c8a6df4c4aab44b4b2f1aadc0615265fe8 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 0000000000000000000000000000000000000000..d711f400a869b4d2a9431e9425ab7b5e53e09b5b
--- /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 ae0ba69aff214a7d0f4a9957113bdac52939315d..ef976be76e21b8d28ecc307a58a77161ffa2804c 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 112b460625c9df4a436d55863daa989986611243..7d4c0c8a430df6dec0eb21a8246b848dcc41cc39 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 73e5a4e64a0e85eda1a5571503d28c8b221f2cdc..e98ebdd17780b385783eff189538c0ae90f6dcc6 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 fad94c13687e3a8ff5c1989401bb39f9bf92f505..93f81404659e8ef9aa869ad59b9a78885381a70a 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 0000000000000000000000000000000000000000..5d22e82c2e9d599cb695970ec9ba7049efdf08bf
--- /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 dab629604c0b51f2be9233436af42fe7dee7134b..a277bdea3c7d33051ddf228b3ce95b052655b454 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 1b0850040a1fd944c5e3d3ba361f7d764bece3dd..d56d44369bb2e5492c078a35a0cf0a01e5524d42 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 0000000000000000000000000000000000000000..0415582b728fb614019e5452f73fc7d22e371001
--- /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 80c3d205200a7604840f48282b4b8cff8c137dc8..aa2c5e2ed67140d47c7d4ddf1f57caecc2f28360 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 277c0d5ccd3a7c80c622a4a8f54a20daf9ff13d1..342a384e7f0ea4989b53fd9fd39ede1d0973171c 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 c0adbd6c4e86cff510924018f80b5fb27229001a..1b6fc35f6ee02dc7a441f9dc956165af32c33276 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 2fe733d540e89c7f7e1bdea6e801e221e6a96eb8..0ce785cd9e6bd3f3985466d6f2f90e38fc49b802 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 2f29120b596f481bdd437cb42e7fd2fc665cf911..6f6a8e143b5b181bed10d76223f809851cda4229 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 1ef88c67af2968ba891874da558e206a3dd86964..589b1d380d6ba1fbac10a61290552fe85ab1d193 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 2ed97075717f8984afb22799648b4b26fbc9327b..489b8e795d843f9677daf562214bd0266f803baa 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 daf7646ee266c7c874c17b6b7403388fc33c3ce6..e19073f1e93bd70d649f66f7ca0bd2f49ef8a102 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 196938400cab3735d0633b1d1ba41efc00358f4b..3f174398186c2b5fd0bcc56aef017f40936f9704 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 0000000000000000000000000000000000000000..4d8fca46dac3d6ccd0c3e88a5d6fa5fc9c065fc4
--- /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 0000000000000000000000000000000000000000..8689b1086a816b10552b4b8c6f6ae78bcced7670
--- /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 0000000000000000000000000000000000000000..d37e9ab403421c7f706583537b4ca724f25f4774
--- /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");
+ }
+}