diff --git a/.gitignore b/.gitignore
index cb98282bfb904d03ac706cc8740ee1ad922c8851..066e7213900207034976978af9c91ee0dca2317d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,7 +1,7 @@
*.o
*.a
*~
-build/
+build*/
CMakeFiles/
CMakeCache.txt
.\#*
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index f60cedd06cbea496f19ca991c47dba984fa0611e..43b6e5c3ddf1844d2d3390215c03ca92e8c9e5cd 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -18,8 +18,8 @@
#include <mesh/MeshRandomizer.hpp>
#include <mesh/MeshSmoother.hpp>
#include <mesh/MeshTraits.hpp>
-#include <scheme/AcousticSolver.hpp>
#include <scheme/AcousticCompositeSolver.hpp>
+#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
@@ -40,8 +40,10 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/InflowBoundaryConditionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolver.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VariableBCDescriptor.hpp>
#include <utils/Socket.hpp>
@@ -394,6 +396,18 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("inflow_list",
+ std::function(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const std::vector<FunctionSymbolId>& function_id_list)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<InflowListBoundaryConditionDescriptor>(boundary,
+ function_id_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("outflow", std::function(
[](std::shared_ptr<const IBoundaryDescriptor> boundary)
@@ -420,29 +434,29 @@ SchemeModule::SchemeModule()
}
));
-
- this->_addBuiltinFunction("composite_glace_fluxes", std::function(
- [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list)
- -> std::tuple<std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>,
- std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>,
- std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>
- > {
- return AcousticCompositeSolverHandler{getCommonMesh({rho, c, u, p})}
- .solver()
- .compute_fluxes(AcousticCompositeSolverHandler::SolverType::GlaceComposite, rho, c, u,
- p, bc_descriptor_list);
- }
+ this->_addBuiltinFunction("composite_glace_fluxes",
+ std::function(
- ));
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list)
+ -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>,
+ std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>,
+ std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>> {
+ return AcousticCompositeSolverHandler{getCommonMesh({rho, c, u, p})}
+ .solver()
+ .compute_fluxes(AcousticCompositeSolverHandler::SolverType::GlaceComposite, rho, c, u,
+ p, bc_descriptor_list);
+ }
+
+ ));
this->_addBuiltinFunction("glace_solver",
std::function(
@@ -466,8 +480,7 @@ SchemeModule::SchemeModule()
));
-
- this->_addBuiltinFunction("composite_glace_solver",
+ this->_addBuiltinFunction("composite_glace_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -483,14 +496,30 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return AcousticCompositeSolverHandler{getCommonMesh({rho, u, E, c, p})}
.solver()
- .apply(AcousticCompositeSolverHandler::SolverType::GlaceComposite, dt, rho, u, E, c, p,
- bc_descriptor_list);
+ .apply(AcousticCompositeSolverHandler::SolverType::GlaceComposite, dt, rho, u, E, c,
+ p, bc_descriptor_list);
}
));
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version1",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver(rho, u, E, c, p, bc_descriptor_list, dt);
+ }
+
+ ));
-
this->_addBuiltinFunction("eucclhyd_fluxes",
std::function(
@@ -510,7 +539,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("composite_eucclhyd_fluxes",
+ this->_addBuiltinFunction("composite_eucclhyd_fluxes",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -520,22 +549,19 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list)
-> std::tuple<std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>,
- std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>,
- std::shared_ptr<const ItemValueVariant>,
- std::shared_ptr<const SubItemValuePerItemVariant>
- > {
+ std::shared_ptr<const SubItemValuePerItemVariant>,
+ std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>,
+ std::shared_ptr<const ItemValueVariant>,
+ std::shared_ptr<const SubItemValuePerItemVariant>> {
return AcousticCompositeSolverHandler{getCommonMesh({rho, c, u, p})}
.solver()
- .compute_fluxes(AcousticCompositeSolverHandler::SolverType::EucclhydComposite, rho, c, u, p,
- bc_descriptor_list);
+ .compute_fluxes(AcousticCompositeSolverHandler::SolverType::EucclhydComposite, rho, c,
+ u, p, bc_descriptor_list);
}
));
-
-
this->_addBuiltinFunction("eucclhyd_solver",
std::function(
@@ -574,13 +600,12 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return AcousticCompositeSolverHandler{getCommonMesh({rho, u, E, c, p})}
.solver()
- .apply(AcousticCompositeSolverHandler::SolverType::EucclhydComposite, dt, rho, u, E, c, p,
- bc_descriptor_list);
+ .apply(AcousticCompositeSolverHandler::SolverType::EucclhydComposite, dt, rho, u, E,
+ c, p, bc_descriptor_list);
}
));
-
this->_addBuiltinFunction("apply_acoustic_fluxes",
std::function(
@@ -599,7 +624,7 @@ SchemeModule::SchemeModule()
}
));
- this->_addBuiltinFunction("apply_acoustic_composite_fluxes",
+ this->_addBuiltinFunction("apply_acoustic_composite_fluxes",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho, //
@@ -607,17 +632,17 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& E, //
const std::shared_ptr<const ItemValueVariant>& ur, //
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, //
- const std::shared_ptr<const ItemValueVariant>& ue, //
+ const std::shared_ptr<const ItemValueVariant>& ue, //
const std::shared_ptr<const SubItemValuePerItemVariant>& Fje, //
- const std::shared_ptr<const ItemValueVariant>& uf, //
- const std::shared_ptr<const SubItemValuePerItemVariant>& Fjf, //
+ const std::shared_ptr<const ItemValueVariant>& uf, //
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjf, //
const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return AcousticCompositeSolverHandler{getCommonMesh({rho, u, E})} //
- .solver()
- .apply_fluxes(dt, rho, u, E, ur, Fjr, ue, Fje, uf, Fjf);
+ .solver()
+ .apply_fluxes(dt, rho, u, E, ur, Fjr, ue, Fje, uf, Fjf);
}
));
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 38f8e7f29697c476dca64f0a8006f447aef1e24c..a23d778a6a41df49f06af3c00fa95272c79852d7 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -31,27 +31,27 @@ class MeshData<Mesh<Dimension>>
const MeshType& m_mesh;
NodeValuePerFace<const Rd> m_Nlr;
NodeValuePerFace<const Rd> m_nlr;
-
+
NodeValuePerCell<const Rd> m_Cjr;
EdgeValuePerCell<const Rd> m_Cje;
FaceValuePerCell<const Rd> m_Cjf;
-
+
NodeValuePerCell<const double> m_ljr;
NodeValuePerCell<const Rd> m_njr;
EdgeValuePerCell<const double> m_lje;
EdgeValuePerCell<const Rd> m_nje;
FaceValuePerCell<const double> m_ljf;
FaceValuePerCell<const Rd> m_njf;
-
+
FaceValue<const Rd> m_xl;
CellValue<const Rd> m_cell_centroid;
CellValue<const Rd> m_cell_iso_barycenter;
CellValue<const double> m_Vj;
-
+
CellValue<const double> m_sum_r_ljr;
CellValue<const double> m_sum_e_lje;
CellValue<const double> m_sum_f_ljf;
-
+
FaceValue<const Rd> m_Nl;
FaceValue<const Rd> m_nl;
FaceValue<const double> m_ll;
@@ -143,8 +143,7 @@ class MeshData<Mesh<Dimension>>
m_nlr = nlr;
}
}
-
-
+
PUGS_INLINE
void
_computeFaceIsobarycenter()
@@ -505,8 +504,6 @@ class MeshData<Mesh<Dimension>>
}
}
-
-
PUGS_INLINE
void
_computeCjf()
@@ -528,9 +525,9 @@ class MeshData<Mesh<Dimension>>
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
for (size_t R = 0; R < cell_nodes.size(); ++R) {
- int Rp = (R + 1) % cell_nodes.size();
+ int Rp = (R + 1) % cell_nodes.size();
const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
- Cjf(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
+ Cjf(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
}
});
m_Cjf = Cjf;
@@ -545,16 +542,15 @@ class MeshData<Mesh<Dimension>>
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-
- const auto& cell_faces = cell_to_face_matrix[j];
+ const auto& cell_faces = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_faces.size(); ++l) {
- const FaceId& F = cell_faces[l];
+ const FaceId& F = cell_faces[l];
- if (cell_face_is_reversed(j,l))
- Cjf(j,l)=-Nl[F];
- else
- Cjf(j,l)=Nl[F];
+ if (cell_face_is_reversed(j, l))
+ Cjf(j, l) = -Nl[F];
+ else
+ Cjf(j, l) = Nl[F];
}
});
@@ -562,15 +558,14 @@ class MeshData<Mesh<Dimension>>
}
}
-
- PUGS_INLINE
+ PUGS_INLINE
void
- _compute_ljf()
+ _compute_ljf()
{
auto Cjf = this->Cjf();
-
+
if constexpr (Dimension == 1) {
- FaceValuePerCell<double> ljf(m_mesh.connectivity());
+ FaceValuePerCell<double> ljf(m_mesh.connectivity());
parallel_for(
ljf.numberOfValues(), PUGS_LAMBDA(size_t jf) { ljf[jf] = 1; });
m_ljf = ljf;
@@ -583,8 +578,7 @@ class MeshData<Mesh<Dimension>>
}
}
-
- PUGS_INLINE
+ PUGS_INLINE
void
_compute_njf()
{
@@ -602,7 +596,6 @@ class MeshData<Mesh<Dimension>>
}
}
-
PUGS_INLINE
void
_computeCje()
@@ -624,17 +617,17 @@ class MeshData<Mesh<Dimension>>
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
for (size_t R = 0; R < cell_nodes.size(); ++R) {
- int Rp = (R + 1) % cell_nodes.size();
+ int Rp = (R + 1) % cell_nodes.size();
const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
- Cje(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
+ Cje(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
}
});
m_Cje = Cje;
} else if (Dimension == 3) {
auto Nlr = this->Nlr();
- const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
- //const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ // const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
const auto& cell_to_edge_matrix = m_mesh.connectivity().cellToEdgeMatrix();
const auto& edge_to_face_matrix = m_mesh.connectivity().edgeToFaceMatrix();
@@ -646,63 +639,61 @@ class MeshData<Mesh<Dimension>>
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_faces = cell_to_face_matrix[j];
- const auto& cell_edges = cell_to_edge_matrix[j];
-
+ const auto& cell_faces = cell_to_face_matrix[j];
+ const auto& cell_edges = cell_to_edge_matrix[j];
+
for (size_t i_edge = 0; i_edge < cell_edges.size(); ++i_edge) {
- const EdgeId edge_id = cell_edges[i_edge];
- auto edges_faces = edge_to_face_matrix[edge_id];
- auto edges_nodes = edge_to_node_matrix[edge_id];
-
- const NodeId& node0_id= edges_nodes[0];
- const NodeId& node1_id= edges_nodes[1];
-
- int nbFacEdge=0;
-
- // On ne garde que si les 2 faces qui sont dans la maille
- //bool findTwoFaces=false;
- for (size_t i_face = 0; i_face < edges_faces.size(); ++i_face) {
- const FaceId face_id = edges_faces[i_face];
-
- for (size_t i_face_cell = 0; i_face_cell < cell_faces.size(); ++i_face_cell) {
- const FaceId face_id_cell = cell_faces[i_face_cell];
-
- if (face_id_cell==face_id) {
-
- double sign = (cell_face_is_reversed(j, i_face_cell)) ? -1 : 1;
-
- const auto& face_nodes = face_to_node_matrix[face_id_cell];
-
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- if (node0_id==face_nodes[rl] or node1_id==face_nodes[rl])
- Cje(j,i_edge)+=.5*sign*Nlr(face_id_cell,rl);
- }
- //findTwoFaces=true;
- ++nbFacEdge;
- break;
- }
- }
- }
-
- if (nbFacEdge!=2)
- {
- throw NormalError("Probleme finding edge mesh ") ;
- }
- }
+ const EdgeId edge_id = cell_edges[i_edge];
+ auto edges_faces = edge_to_face_matrix[edge_id];
+ auto edges_nodes = edge_to_node_matrix[edge_id];
+
+ const NodeId& node0_id = edges_nodes[0];
+ const NodeId& node1_id = edges_nodes[1];
+
+ int nbFacEdge = 0;
+
+ // On ne garde que si les 2 faces qui sont dans la maille
+ // bool findTwoFaces=false;
+ for (size_t i_face = 0; i_face < edges_faces.size(); ++i_face) {
+ const FaceId face_id = edges_faces[i_face];
+
+ for (size_t i_face_cell = 0; i_face_cell < cell_faces.size(); ++i_face_cell) {
+ const FaceId face_id_cell = cell_faces[i_face_cell];
+
+ if (face_id_cell == face_id) {
+ double sign = (cell_face_is_reversed(j, i_face_cell)) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_matrix[face_id_cell];
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ if (node0_id == face_nodes[rl] or node1_id == face_nodes[rl])
+ Cje(j, i_edge) += .5 * sign * Nlr(face_id_cell, rl);
+ }
+ // findTwoFaces=true;
+ ++nbFacEdge;
+ break;
+ }
+ }
+ }
+
+ if (nbFacEdge != 2) {
+ throw NormalError("Probleme finding edge mesh ");
+ }
+ }
});
-
+
m_Cje = Cje;
}
- }
+ }
- PUGS_INLINE
+ PUGS_INLINE
void
- _compute_lje()
+ _compute_lje()
{
auto Cje = this->Cje();
-
+
if constexpr (Dimension == 1) {
- EdgeValuePerCell<double> lje(m_mesh.connectivity());
+ EdgeValuePerCell<double> lje(m_mesh.connectivity());
parallel_for(
lje.numberOfValues(), PUGS_LAMBDA(size_t je) { lje[je] = 1; });
m_lje = lje;
@@ -715,8 +706,7 @@ class MeshData<Mesh<Dimension>>
}
}
-
- PUGS_INLINE
+ PUGS_INLINE
void
_compute_nje()
{
@@ -734,9 +724,6 @@ class MeshData<Mesh<Dimension>>
}
}
-
-
-
PUGS_INLINE
void
_computeSumOverRljr()
@@ -764,7 +751,7 @@ class MeshData<Mesh<Dimension>>
m_sum_r_ljr = sum_r_ljr;
}
-
+
PUGS_INLINE
void
_computeSumOverElje()
@@ -821,8 +808,6 @@ class MeshData<Mesh<Dimension>>
m_sum_f_ljf = sum_f_ljf;
}
-
-
void
_checkCellVolume()
{
@@ -930,7 +915,6 @@ class MeshData<Mesh<Dimension>>
return m_njr;
}
-
PUGS_INLINE
EdgeValuePerCell<const Rd>
Cje()
@@ -991,8 +975,6 @@ class MeshData<Mesh<Dimension>>
return m_njf;
}
-
-
PUGS_INLINE
FaceValue<const Rd>
xl()
@@ -1053,7 +1035,7 @@ class MeshData<Mesh<Dimension>>
}
return m_sum_r_ljr;
}
-
+
PUGS_INLINE
CellValue<const double>
sumOverELje()
@@ -1063,7 +1045,7 @@ class MeshData<Mesh<Dimension>>
}
return m_sum_e_lje;
}
-
+
PUGS_INLINE
CellValue<const double>
sumOverFLjf()
@@ -1074,7 +1056,6 @@ class MeshData<Mesh<Dimension>>
return m_sum_f_ljf;
}
-
private:
// MeshData **must** be constructed through MeshDataManager
friend class MeshDataManager;
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 8cdd43f9e51e3b36c0c93e4d4cbfcbdfd9325f71..e285dbde1250fae2ebe063291db6bba6cc8bc2ec 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -11,6 +11,7 @@ add_library(
DiscreteFunctionVectorIntegrator.cpp
DiscreteFunctionVectorInterpoler.cpp
FluxingAdvectionSolver.cpp
+ RusanovEulerianCompositeSolver.cpp
)
target_link_libraries(
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 1350aa4dba4f7f6739b7de6689c10554ceec6971..e1169aed79c785ffe7c768ee4762b50e122dabb2 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -17,6 +17,7 @@ class IBoundaryConditionDescriptor
fixed,
free,
inflow,
+ inflow_list,
neumann,
outflow,
symmetry
diff --git a/src/scheme/InflowListBoundaryConditionDescriptor.hpp b/src/scheme/InflowListBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3effc4315d886930708de629713f28a3010f6365
--- /dev/null
+++ b/src/scheme/InflowListBoundaryConditionDescriptor.hpp
@@ -0,0 +1,55 @@
+#ifndef INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+
+class InflowListBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "inflow_list(" << ',' << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+ const std::vector<FunctionSymbolId> m_function_symbol_id_list;
+
+ public:
+ const std::vector<FunctionSymbolId>&
+ functionSymbolIdList() const
+ {
+ return m_function_symbol_id_list;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const final
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::inflow_list;
+ }
+
+ InflowListBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const std::vector<FunctionSymbolId>& function_symbol_id_list)
+ : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id_list{function_symbol_id_list}
+ {
+ ;
+ }
+
+ InflowListBoundaryConditionDescriptor(const InflowListBoundaryConditionDescriptor&) = delete;
+ InflowListBoundaryConditionDescriptor(InflowListBoundaryConditionDescriptor&&) = delete;
+
+ ~InflowListBoundaryConditionDescriptor() = default;
+};
+
+#endif // INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..633c324cd3bc0bcd1f0b2e9482e78dec5aed1c96
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -0,0 +1,499 @@
+#include <scheme/RusanovEulerianCompositeSolver.hpp>
+
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshEdgeBoundary.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatEdgeBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class RusanovEulerianCompositeSolver
+{
+ private:
+ using MeshType = MeshTypeT;
+
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ class SymmetryBoundaryCondition;
+ class InflowListBoundaryCondition;
+ class OutflowBoundaryCondition;
+
+ using BoundaryCondition =
+ std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList
+ _getBCList(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ BoundaryConditionList bc_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::inflow_list: {
+ const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
+ dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
+ std::ostringstream error_msg;
+ error_msg << "invalid number of functions for inflow boundary "
+ << inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
+ << inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
+ throw NormalError(error_msg.str());
+ }
+
+ if constexpr (Dimension == 2) {
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> node_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ mesh.xr(), node_boundary.nodeList());
+
+ auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
+
+ auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> face_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xl, face_boundary.faceList());
+
+ bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, face_boundary, node_values, face_values));
+ } else {
+ static_assert(Dimension == 3);
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> node_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ mesh.xr(), node_boundary.nodeList());
+
+ auto xe = MeshDataManager::instance().getMeshData(mesh).xe();
+
+ auto edge_boundary = getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> edge_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::edge>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xe, edge_boundary.edgeList());
+
+ auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
+
+ auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> face_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xl, face_boundary.faceList());
+
+ bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, edge_boundary, face_boundary, node_values,
+ edge_values, face_values));
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::outflow: {
+ std::cout << "outflow not implemented yet\n";
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for Rusanov Eulerian Composite solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ public:
+ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ solve(const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteFunctionP0<const double>& rho_n,
+ const DiscreteFunctionP0<const Rd>& u_n,
+ const DiscreteFunctionP0<const double>& E_n,
+ const DiscreteFunctionP0<const double>& c_n,
+ const DiscreteFunctionP0<const double>& p_n,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt) const
+ {
+ auto rho = copy(rho_n);
+ auto u = copy(u_n);
+ auto E = copy(E_n);
+ auto c = copy(c_n);
+ auto p = copy(p_n);
+
+ auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list);
+
+ return std::make_tuple(std::make_shared<const DiscreteFunctionVariant>(rho),
+ std::make_shared<const DiscreteFunctionVariant>(u),
+ std::make_shared<const DiscreteFunctionVariant>(E));
+ }
+
+ RusanovEulerianCompositeSolver() = default;
+ ~RusanovEulerianCompositeSolver() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver<MeshType>::SymmetryBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<2>>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
+ const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<3>>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary;
+ const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_flat_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_flat_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
+ const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
+ const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary),
+ m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
+ m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<2>>::InflowListBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Table<const double> m_node_array_list;
+ const Table<const double> m_face_array_list;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Table<const double>&
+ nodeArrayList() const
+ {
+ return m_node_array_list;
+ }
+
+ const Table<const double>&
+ faceArrayList() const
+ {
+ return m_face_array_list;
+ }
+
+ InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshFaceBoundary& mesh_face_boundary,
+ const Table<const double>& node_array_list,
+ const Table<const double>& face_array_list)
+ : m_mesh_node_boundary(mesh_node_boundary),
+ m_mesh_face_boundary(mesh_face_boundary),
+ m_node_array_list(node_array_list),
+ m_face_array_list(face_array_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<3>>::InflowListBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshEdgeBoundary m_mesh_edge_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Table<const double> m_node_array_list;
+ const Table<const double> m_edge_array_list;
+ const Table<const double> m_face_array_list;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Table<const double>&
+ nodeArrayList() const
+ {
+ return m_node_array_list;
+ }
+
+ const Table<const double>&
+ edgeArrayList() const
+ {
+ return m_edge_array_list;
+ }
+
+ const Table<const double>&
+ faceArrayList() const
+ {
+ return m_face_array_list;
+ }
+
+ InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshEdgeBoundary& mesh_edge_boundary,
+ const MeshFaceBoundary& mesh_face_boundary,
+ const Table<const double>& node_array_list,
+ const Table<const double>& edge_array_list,
+ const Table<const double>& face_array_list)
+ : m_mesh_node_boundary(mesh_node_boundary),
+ m_mesh_edge_boundary(mesh_edge_boundary),
+ m_mesh_face_boundary(mesh_face_boundary),
+ m_node_array_list(node_array_list),
+ m_edge_array_list(edge_array_list),
+ m_face_array_list(face_array_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver<MeshType>::OutflowBoundaryCondition
+{
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver(
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt)
+{
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v, c_v, p_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ return std::visit(
+ PUGS_LAMBDA(auto&& p_mesh)
+ ->std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ if constexpr (Dimension == 1) {
+ throw NormalError("RusanovEulerianCompositeSolver is not available in 1D");
+ } else {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ return RusanovEulerianCompositeSolver<MeshType>{}.solve(p_mesh,
+ rho_v->get<DiscreteFunctionP0<const double>>(),
+ u_v->get<DiscreteFunctionP0<const Rd>>(),
+ E_v->get<DiscreteFunctionP0<const double>>(),
+ c_v->get<DiscreteFunctionP0<const double>>(),
+ p_v->get<DiscreteFunctionP0<const double>>(),
+ bc_descriptor_list, dt);
+ } else {
+ throw NormalError("RusanovEulerianCompositeSolver is only defined on polygonal meshes");
+ }
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolver.hpp b/src/scheme/RusanovEulerianCompositeSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..69dfb6be8c4b762dad42414feebb20d2eec0b337
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver.hpp
@@ -0,0 +1,23 @@
+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver(
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt);
+
+#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP
diff --git a/tests/test_TestVolumes3D.cpp b/tests/test_TestVolumes3D.cpp
index 2bf95290526dbd3e6816f4cee2dd53a6640cb6bc..7006dc7ebeb2a17ba25615e57304af00b9b8c0df 100644
--- a/tests/test_TestVolumes3D.cpp
+++ b/tests/test_TestVolumes3D.cpp
@@ -1,8 +1,8 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-#include <mesh/ItemArrayUtils.hpp>
#include <mesh/DualMeshManager.hpp>
+#include <mesh/ItemArrayUtils.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
@@ -16,278 +16,263 @@
TEST_CASE("TestVolumes3D", "[scheme]")
{
SECTION("3D")
- {
- using R3 = TinyVector<3>;
-
- auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
- /*
- auto f = [](const R3& X) -> double {
- const double x = X[0];
- const double y = X[1];
- const double z = X[2];
- return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
- };
- */
- std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
- mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
- mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(hybrid_mesh)));
-
- for (const auto& mesh_info : mesh_list) {
- auto mesh_name = mesh_info.first;
- auto mesh = mesh_info.second->get<Mesh<3>>();
- std::clog << " name " << mesh_name << "\n";
- auto& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
- auto Cjr = mesh_data.Cjr();
- auto Cje = mesh_data.Cje();
- auto Cjf = mesh_data.Cjf();
-
- auto xr = mesh->xr();
- auto Vj = mesh_data.Vj();
-
- auto xl = mesh_data.xl(); // Les centres des faces
- auto xe = mesh_data.xe(); // Les centres des aretes
- auto Nl = mesh_data.Nl(); // Normale aux faces .. au "centre" des faces
- auto nl = mesh_data.nl(); // Normale unite aux faces .. au "centre" des faces
-
- NodeValuePerFace<const R3> Nlr=mesh_data.Nlr();
-
- auto cell_to_node_connectivity = mesh->connectivity().cellToNodeMatrix();
- auto cell_to_face_connectivity = mesh->connectivity().cellToFaceMatrix();
-
- const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
-
- auto cell_to_edge_connectivity = mesh->connectivity().cellToEdgeMatrix();
- auto edge_to_node_connectivity = mesh->connectivity().edgeToNodeMatrix();
- auto edge_to_face_connectivity = mesh->connectivity().edgeToFaceMatrix();
- auto face_to_node_connectivity = mesh->connectivity().faceToNodeMatrix();
-
- FaceValuePerCell<R3> face_normal_per_cell(mesh->connectivity());
- face_normal_per_cell.fill(zero);
- EdgeValuePerCell<R3> edge_normal_per_cell(mesh->connectivity());
- edge_normal_per_cell.fill(zero);
- CellValue<double> tabErrVolCjr(mesh->connectivity());
- CellValue<double> tabErrVolCjf(mesh->connectivity());
- CellValue<double> tabErrVolCje(mesh->connectivity());
- tabErrVolCjr.fill(0.);
- tabErrVolCjf.fill(0.);
- tabErrVolCje.fill(0.);
-
- CellValue<double> tabErrStokesCjr(mesh->connectivity());
- CellValue<double> tabErrStokesCjf(mesh->connectivity());
- CellValue<double> tabErrStokesCje(mesh->connectivity());
- tabErrStokesCjr.fill(0.);
- tabErrStokesCjf.fill(0.);
- tabErrStokesCje.fill(0.);
-
- NodeValue<double> tabErrConservation_r_Cjr(mesh->connectivity());
- FaceValue<double> tabErrConservation_f_Cjf(mesh->connectivity());
- EdgeValue<double> tabErrConservation_e_Cje(mesh->connectivity());
- tabErrConservation_r_Cjr.fill(0.);
- tabErrConservation_f_Cjf.fill(0.);
- tabErrConservation_e_Cje.fill(0.);
-
- parallel_for(mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-
- // ********************************
- // Travail sur les nodes
- // ********************************
-
- auto cell_nodes = cell_to_node_connectivity[cell_id];
- double SumCjrXr=0;
- R3 SumCjr=R3{0,0,0};
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- const NodeId node_id = cell_nodes[i_node];
- SumCjrXr+=dot(Cjr(cell_id,i_node),xr[node_id]);
- SumCjr+=Cjr(cell_id,i_node);
- }
- SumCjrXr/=3.;
- tabErrVolCjr[cell_id]=fabs(Vj[cell_id]-SumCjrXr);
- tabErrStokesCjr[cell_id]=l2Norm(SumCjr);
-
- // ********************************
- // Travail sur les faces
- // ********************************
-
- auto cell_faces = cell_to_face_connectivity[cell_id];
- const auto& face_is_reversed = cell_face_is_reversed.itemArray(cell_id);
-
- double SumCjrXf=0;
- R3 SumCjf=R3{0,0,0};
- double SumCjrXf_mesh=0;
- R3 SumCjf_mesh=R3{0,0,0};
- for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
- const FaceId face_id = cell_faces[i_face];
- SumCjrXf_mesh+=dot(Cjf(cell_id,i_face),xl[face_id]);
- SumCjf_mesh+=Cjf(cell_id,i_face);
-
- if (face_is_reversed[i_face])
- {
- face_normal_per_cell[cell_id][i_face]=-Nl[face_id];
- SumCjrXf-=dot(Nl[face_id],xl[face_id]);
- SumCjf-=Nl[face_id];
- }
- else
- {
- face_normal_per_cell[cell_id][i_face]=Nl[face_id];
- SumCjrXf+=dot(Nl[face_id],xl[face_id]);
- SumCjf+=Nl[face_id];
- }
- }
- SumCjrXf/=3.;
- SumCjrXf_mesh/=3.;
- tabErrVolCjf[cell_id]=fabs(Vj[cell_id]-SumCjrXf);
- tabErrStokesCjf[cell_id]=l2Norm(SumCjf);
-
- // ********************************
- // Travail sur les edges
- // ********************************
- //auto cell_faces = cell_to_face_connectivity[cell_id];
- //const auto& face_is_reversed = cell_face_is_reversed.itemArray(cell_id);
- //const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
-
- auto cell_edges = cell_to_edge_connectivity[cell_id];
- double SumCjrXe=0;
- R3 SumCje=zero;
- double SumCjrXe_mesh=0;
- R3 SumCje_mesh=zero;
-
- for (size_t i_edge = 0; i_edge < cell_edges.size(); ++i_edge) {
- const EdgeId edge_id = cell_edges[i_edge];
-
- SumCjrXe_mesh+=dot(Cje(cell_id,i_edge),xe[edge_id]);
- SumCje_mesh+=Cje(cell_id,i_edge);
-
- auto edges_faces = edge_to_face_connectivity[edge_id];
- //auto cell_faces = cell_to_face_connectivity[edge_id];
- auto edges_nodes = edge_to_node_connectivity[edge_id];
-
- const NodeId& node0_id= edges_nodes[0];
- const NodeId& node1_id= edges_nodes[1];
-
- int nbFacEdge=0;
- // On ne garde que si les 2 faces qui sont dans la maille
- for (size_t i_face = 0; i_face < edges_faces.size(); ++i_face) {
- const FaceId face_id = edges_faces[i_face];
-
- for (size_t i_face_cell = 0; i_face_cell < cell_faces.size(); ++i_face_cell) {
- const FaceId face_id_cell = cell_faces[i_face_cell];
-
- // La face autour de l arete appartient à la maille courante : On peut faire le calcul avec la face
- if (face_id_cell==face_id)
- {
-
- double sign = (cell_face_is_reversed(cell_id, i_face_cell)) ? -1 : 1;
-
- const auto& face_nodes = face_to_node_connectivity[face_id_cell];
-
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- if (node0_id==face_nodes[rl] or node1_id==face_nodes[rl])
- edge_normal_per_cell(cell_id,i_edge)+=.5*sign*Nlr(face_id_cell,rl);
- }
-
- ++nbFacEdge;
- break;
- }
- }
- }
-
- if (nbFacEdge!=2)
- {
- exit(1);
- }
-
- R3 Cje_l=edge_normal_per_cell[cell_id][i_edge];
- //R3 Cje_l=Cje(cell_id,i_edge);//edge_normal_per_cell[cell_id][i_edge];
-
- SumCjrXe+=dot(Cje_l,xe[edge_id]);
- SumCje+=Cje_l;
- }
- SumCjrXe/=3.;
- SumCjrXe_mesh/=3.;
-
- tabErrVolCje[cell_id]=fabs(Vj[cell_id]-SumCjrXe);
- tabErrStokesCje[cell_id]=l2Norm(SumCje);
-
- });
-
- auto node_to_cell_connectivity = mesh->connectivity().nodeToCellMatrix();
- auto face_to_cell_connectivity = mesh->connectivity().faceToCellMatrix();
- auto edge_to_cell_connectivity = mesh->connectivity().edgeToCellMatrix();
-
- auto is_boundary_node = mesh->connectivity().isBoundaryNode();
- auto is_boundary_face = mesh->connectivity().isBoundaryFace();
- auto is_boundary_edge = mesh->connectivity().isBoundaryEdge();
-
- parallel_for(mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id)
- {
- R3 SumCjr=R3{0,0,0};
- if (not is_boundary_node[node_id])
- {
-
- for (size_t i_node = 0; i_node < node_to_cell_connectivity[node_id].size(); ++i_node) {
- const CellId cell_id= node_to_cell_connectivity[node_id][i_node];
-
-
- for (size_t i_node_cell = 0; i_node_cell < cell_to_node_connectivity[cell_id].size(); ++i_node_cell) {
- const NodeId node_cell_id = cell_to_node_connectivity[cell_id][i_node_cell];
- if (node_id==node_cell_id)
- {
- SumCjr+=Cjr(cell_id,i_node_cell);
- break;
- }
- }
- }
- }
- tabErrConservation_r_Cjr[node_id]=l2Norm(SumCjr);
- });
-
- parallel_for(mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
- R3 SumCjf=zero;
- if (not is_boundary_face[face_id])
- {
- for (size_t i_face = 0; i_face < face_to_cell_connectivity[face_id].size(); ++i_face) {
- const CellId cell_id= face_to_cell_connectivity[face_id][i_face];
-
- for (size_t i_face_cell = 0; i_face_cell < cell_to_face_connectivity[cell_id].size(); ++i_face_cell) {
- const FaceId face_cell_id = cell_to_face_connectivity[cell_id][i_face_cell];
- if (face_id==face_cell_id)
- {
- SumCjf+=face_normal_per_cell[cell_id][i_face_cell];
- break;
- }
- }
- }
- }
- tabErrConservation_f_Cjf[face_id]=l2Norm(SumCjf)+1e-13;//R3(1e-13,1e-56,1e-34};
- });
-
- parallel_for(mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
- R3 SumCje=zero;
-
- if (not is_boundary_edge[edge_id])
- {
- for (size_t i_edge = 0; i_edge < edge_to_cell_connectivity[edge_id].size(); ++i_edge) {
- const CellId cell_id= edge_to_cell_connectivity[edge_id][i_edge];
-
- for (size_t i_edge_cell = 0; i_edge_cell < cell_to_edge_connectivity[cell_id].size(); ++i_edge_cell) {
- const EdgeId edge_cell_id = cell_to_edge_connectivity[cell_id][i_edge_cell];
- if (edge_id==edge_cell_id)
- {
- SumCje+=edge_normal_per_cell[cell_id][i_edge_cell];
- //SumCje+=Cje(cell_id,i_edge_cell);//edge_normal_per_cell[cell_id][i_edge_cell];
- break;
- }
- }
- }
- }
- tabErrConservation_e_Cje[edge_id]=l2Norm(SumCje);
- });
- //std::clog << " Max Cje et tab " << max(dif_edge_normal_per_cell) << "\n" ;
- std::clog << " Max Error (Stokes/Vol/Cons) Cjr " << max(tabErrStokesCjr) << "/" << max(tabErrVolCjr) << "/" << max(tabErrConservation_r_Cjr) << " Max Error (Stokes/Vol/Cons) Cjf " << max(tabErrStokesCjf) << "/" << max(tabErrVolCjf) << "/" << max(tabErrConservation_f_Cjf) << " Max Error (Stokes/Vol/Cons) Cje "<< max(tabErrStokesCje) << "/" << max(tabErrVolCje) << "/"<< max(tabErrConservation_e_Cje) << "\n";
-
- //1.00868e-16/1.73472e-18/2.399e-17 Max Error (Stokes/Vol/Cons) Cjf 9.10193e-17/1.82146e-17/2.22507e-308 Max Error (Stokes/Vol/Cons) Cje
-
- }
+ {
+ using R3 = TinyVector<3>;
+
+ auto hybrid_mesh = MeshDataBaseForTests::get().hybrid3DMesh();
+ /*
+ auto f = [](const R3& X) -> double {
+ const double x = X[0];
+ const double y = X[1];
+ const double z = X[2];
+ return x * x + 2 * x * y + 3 * y * y + 2 * z * z - z + 1;
+ };
+ */
+ std::vector<std::pair<std::string, decltype(hybrid_mesh)>> mesh_list;
+ mesh_list.push_back(std::make_pair("hybrid mesh", hybrid_mesh));
+ mesh_list.push_back(std::make_pair("diamond mesh", DualMeshManager::instance().getDiamondDualMesh(hybrid_mesh)));
+
+ for (const auto& mesh_info : mesh_list) {
+ auto mesh_name = mesh_info.first;
+ auto mesh = mesh_info.second->get<Mesh<3>>();
+
+ auto& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ auto Cjr = mesh_data.Cjr();
+ auto Cje = mesh_data.Cje();
+ auto Cjf = mesh_data.Cjf();
+
+ auto xr = mesh->xr();
+ auto Vj = mesh_data.Vj();
+
+ auto xl = mesh_data.xl(); // Les centres des faces
+ auto xe = mesh_data.xe(); // Les centres des aretes
+ auto Nl = mesh_data.Nl(); // Normale aux faces .. au "centre" des faces
+ auto nl = mesh_data.nl(); // Normale unite aux faces .. au "centre" des faces
+
+ NodeValuePerFace<const R3> Nlr = mesh_data.Nlr();
+
+ auto cell_to_node_connectivity = mesh->connectivity().cellToNodeMatrix();
+ auto cell_to_face_connectivity = mesh->connectivity().cellToFaceMatrix();
+
+ const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+ auto cell_to_edge_connectivity = mesh->connectivity().cellToEdgeMatrix();
+ auto edge_to_node_connectivity = mesh->connectivity().edgeToNodeMatrix();
+ auto edge_to_face_connectivity = mesh->connectivity().edgeToFaceMatrix();
+ auto face_to_node_connectivity = mesh->connectivity().faceToNodeMatrix();
+
+ FaceValuePerCell<R3> face_normal_per_cell(mesh->connectivity());
+ face_normal_per_cell.fill(zero);
+ EdgeValuePerCell<R3> edge_normal_per_cell(mesh->connectivity());
+ edge_normal_per_cell.fill(zero);
+ CellValue<double> tabErrVolCjr(mesh->connectivity());
+ CellValue<double> tabErrVolCjf(mesh->connectivity());
+ CellValue<double> tabErrVolCje(mesh->connectivity());
+ tabErrVolCjr.fill(0.);
+ tabErrVolCjf.fill(0.);
+ tabErrVolCje.fill(0.);
+
+ CellValue<double> tabErrStokesCjr(mesh->connectivity());
+ CellValue<double> tabErrStokesCjf(mesh->connectivity());
+ CellValue<double> tabErrStokesCje(mesh->connectivity());
+ tabErrStokesCjr.fill(0.);
+ tabErrStokesCjf.fill(0.);
+ tabErrStokesCje.fill(0.);
+
+ NodeValue<double> tabErrConservation_r_Cjr(mesh->connectivity());
+ FaceValue<double> tabErrConservation_f_Cjf(mesh->connectivity());
+ EdgeValue<double> tabErrConservation_e_Cje(mesh->connectivity());
+ tabErrConservation_r_Cjr.fill(0.);
+ tabErrConservation_f_Cjf.fill(0.);
+ tabErrConservation_e_Cje.fill(0.);
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ // ********************************
+ // Travail sur les nodes
+ // ********************************
+
+ auto cell_nodes = cell_to_node_connectivity[cell_id];
+ double SumCjrXr = 0;
+ R3 SumCjr = R3{0, 0, 0};
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ SumCjrXr += dot(Cjr(cell_id, i_node), xr[node_id]);
+ SumCjr += Cjr(cell_id, i_node);
+ }
+ SumCjrXr /= 3.;
+ tabErrVolCjr[cell_id] = fabs(Vj[cell_id] - SumCjrXr);
+ tabErrStokesCjr[cell_id] = l2Norm(SumCjr);
+
+ // ********************************
+ // Travail sur les faces
+ // ********************************
+
+ auto cell_faces = cell_to_face_connectivity[cell_id];
+ const auto& face_is_reversed = cell_face_is_reversed.itemArray(cell_id);
+
+ double SumCjrXf = 0;
+ R3 SumCjf = R3{0, 0, 0};
+ double SumCjrXf_mesh = 0;
+ R3 SumCjf_mesh = R3{0, 0, 0};
+ for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+ const FaceId face_id = cell_faces[i_face];
+ SumCjrXf_mesh += dot(Cjf(cell_id, i_face), xl[face_id]);
+ SumCjf_mesh += Cjf(cell_id, i_face);
+
+ if (face_is_reversed[i_face]) {
+ face_normal_per_cell[cell_id][i_face] = -Nl[face_id];
+ SumCjrXf -= dot(Nl[face_id], xl[face_id]);
+ SumCjf -= Nl[face_id];
+ } else {
+ face_normal_per_cell[cell_id][i_face] = Nl[face_id];
+ SumCjrXf += dot(Nl[face_id], xl[face_id]);
+ SumCjf += Nl[face_id];
+ }
+ }
+ SumCjrXf /= 3.;
+ SumCjrXf_mesh /= 3.;
+ tabErrVolCjf[cell_id] = fabs(Vj[cell_id] - SumCjrXf);
+ tabErrStokesCjf[cell_id] = l2Norm(SumCjf);
+
+ // ********************************
+ // Travail sur les edges
+ // ********************************
+ // auto cell_faces = cell_to_face_connectivity[cell_id];
+ // const auto& face_is_reversed = cell_face_is_reversed.itemArray(cell_id);
+ // const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+ auto cell_edges = cell_to_edge_connectivity[cell_id];
+ double SumCjrXe = 0;
+ R3 SumCje = zero;
+ double SumCjrXe_mesh = 0;
+ R3 SumCje_mesh = zero;
+
+ for (size_t i_edge = 0; i_edge < cell_edges.size(); ++i_edge) {
+ const EdgeId edge_id = cell_edges[i_edge];
+
+ SumCjrXe_mesh += dot(Cje(cell_id, i_edge), xe[edge_id]);
+ SumCje_mesh += Cje(cell_id, i_edge);
+
+ auto edges_faces = edge_to_face_connectivity[edge_id];
+ // auto cell_faces = cell_to_face_connectivity[edge_id];
+ auto edges_nodes = edge_to_node_connectivity[edge_id];
+
+ const NodeId& node0_id = edges_nodes[0];
+ const NodeId& node1_id = edges_nodes[1];
+
+ int nbFacEdge = 0;
+ // On ne garde que si les 2 faces qui sont dans la maille
+ for (size_t i_face = 0; i_face < edges_faces.size(); ++i_face) {
+ const FaceId face_id = edges_faces[i_face];
+
+ for (size_t i_face_cell = 0; i_face_cell < cell_faces.size(); ++i_face_cell) {
+ const FaceId face_id_cell = cell_faces[i_face_cell];
+
+ // La face autour de l arete appartient à la maille courante : On peut faire le calcul avec la face
+ if (face_id_cell == face_id) {
+ double sign = (cell_face_is_reversed(cell_id, i_face_cell)) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_connectivity[face_id_cell];
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ if (node0_id == face_nodes[rl] or node1_id == face_nodes[rl])
+ edge_normal_per_cell(cell_id, i_edge) += .5 * sign * Nlr(face_id_cell, rl);
+ }
+
+ ++nbFacEdge;
+ break;
+ }
+ }
+ }
+
+ if (nbFacEdge != 2) {
+ throw UnexpectedError("invaid number of edges");
+ }
+
+ R3 Cje_l = edge_normal_per_cell[cell_id][i_edge];
+ // R3 Cje_l=Cje(cell_id,i_edge);//edge_normal_per_cell[cell_id][i_edge];
+
+ SumCjrXe += dot(Cje_l, xe[edge_id]);
+ SumCje += Cje_l;
+ }
+ SumCjrXe /= 3.;
+ SumCjrXe_mesh /= 3.;
+
+ tabErrVolCje[cell_id] = fabs(Vj[cell_id] - SumCjrXe);
+ tabErrStokesCje[cell_id] = l2Norm(SumCje);
+ });
+
+ std::clog << " Max Error Vol et Stokes " << max(tabErrVolCje) << " / " << max(tabErrStokesCje);
+ auto node_to_cell_connectivity = mesh->connectivity().nodeToCellMatrix();
+ auto face_to_cell_connectivity = mesh->connectivity().faceToCellMatrix();
+ auto edge_to_cell_connectivity = mesh->connectivity().edgeToCellMatrix();
+
+ auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+ auto is_boundary_face = mesh->connectivity().isBoundaryFace();
+ auto is_boundary_edge = mesh->connectivity().isBoundaryEdge();
+
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ R3 SumCjr = R3{0, 0, 0};
+ if (not is_boundary_node[node_id]) {
+ for (size_t i_node = 0; i_node < node_to_cell_connectivity[node_id].size(); ++i_node) {
+ const CellId cell_id = node_to_cell_connectivity[node_id][i_node];
+
+ for (size_t i_node_cell = 0; i_node_cell < cell_to_node_connectivity[cell_id].size(); ++i_node_cell) {
+ const NodeId node_cell_id = cell_to_node_connectivity[cell_id][i_node_cell];
+ if (node_id == node_cell_id) {
+ SumCjr += Cjr(cell_id, i_node_cell);
+ break;
+ }
+ }
+ }
+ }
+ tabErrConservation_r_Cjr[node_id] = l2Norm(SumCjr);
+ });
+
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(const FaceId face_id) {
+ R3 SumCjf = zero;
+ if (not is_boundary_face[face_id]) {
+ for (size_t i_face = 0; i_face < face_to_cell_connectivity[face_id].size(); ++i_face) {
+ const CellId cell_id = face_to_cell_connectivity[face_id][i_face];
+
+ for (size_t i_face_cell = 0; i_face_cell < cell_to_face_connectivity[cell_id].size(); ++i_face_cell) {
+ const FaceId face_cell_id = cell_to_face_connectivity[cell_id][i_face_cell];
+ if (face_id == face_cell_id) {
+ SumCjf += face_normal_per_cell[cell_id][i_face_cell];
+ break;
+ }
+ }
+ }
+ }
+ tabErrConservation_f_Cjf[face_id] = l2Norm(SumCjf);
+ });
+
+ parallel_for(
+ mesh->numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+ R3 SumCje = zero;
+
+ if (not is_boundary_edge[edge_id]) {
+ for (size_t i_edge = 0; i_edge < edge_to_cell_connectivity[edge_id].size(); ++i_edge) {
+ const CellId cell_id = edge_to_cell_connectivity[edge_id][i_edge];
+
+ for (size_t i_edge_cell = 0; i_edge_cell < cell_to_edge_connectivity[cell_id].size(); ++i_edge_cell) {
+ const EdgeId edge_cell_id = cell_to_edge_connectivity[cell_id][i_edge_cell];
+ if (edge_id == edge_cell_id) {
+ SumCje += edge_normal_per_cell[cell_id][i_edge_cell];
+ // SumCje+=Cje(cell_id,i_edge_cell);//edge_normal_per_cell[cell_id][i_edge_cell];
+ break;
+ }
+ }
+ }
+ }
+ tabErrConservation_e_Cje[edge_id] = l2Norm(SumCje);
+ });
+
+ REQUIRE(sum(tabErrConservation_e_Cje) == Catch::Approx(0).margin(1E-13));
}
+ }
}