diff --git a/src/scheme/AcousticCompositeSolver.cpp b/src/scheme/AcousticCompositeSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c931ba8b8a9465c52801bda4f6d33fd09a4993a
--- /dev/null
+++ b/src/scheme/AcousticCompositeSolver.cpp
@@ -0,0 +1,1523 @@
+#include <scheme/AcousticCompositeSolver.hpp>
+
+#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/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+#include <variant>
+#include <vector>
+
+double
+composite_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+ const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+ return std::visit(
+ [&](auto&& p_mesh) -> double {
+ const auto& mesh = *p_mesh;
+
+ using MeshType = decltype(mesh);
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+ const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+ CellValue<double> local_dt{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
+
+ return min(local_dt);
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ c.meshVariant()->variant());
+}
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver final : public AcousticCompositeSolverHandler::IAcousticCompositeSolver
+{
+ private:
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ using MeshDataType = MeshData<MeshType>;
+
+ using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
+
+ class FixedBoundaryCondition;
+ class PressureBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class VelocityBoundaryCondition;
+ class ExternalFSIVelocityBoundaryCondition;
+
+ using BoundaryCondition = std::variant<FixedBoundaryCondition,
+ PressureBoundaryCondition,
+ SymmetryBoundaryCondition,
+ VelocityBoundaryCondition,
+ ExternalFSIVelocityBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ NodeValuePerCell<const Rdxd>
+ _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
+ NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const size_t& nb_nodes = Ajr.numberOfSubValues(j);
+ const double& rhoc_j = rhoc[j];
+ for (size_t r = 0; r < nb_nodes; ++r) {
+ Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r));
+ }
+ });
+
+ return Ajr;
+ }
+
+ EdgeValuePerCell<const Rdxd>
+ _computeGlaceAje(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
+ const EdgeValuePerCell<const Rd> nje = mesh_data.nje();
+
+ EdgeValuePerCell<Rdxd> Aje{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const size_t& nb_edges = Aje.numberOfSubValues(j);
+ const double& rhoc_j = rhoc[j];
+ for (size_t e = 0; e < nb_edges; ++e) {
+ Aje(j, e) = tensorProduct(rhoc_j * Cje(j, e), nje(j, e));
+ }
+ });
+
+ return Aje;
+ }
+
+ FaceValuePerCell<const Rdxd>
+ _computeGlaceAjf(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf();
+ const FaceValuePerCell<const Rd> njf = mesh_data.njf();
+
+ FaceValuePerCell<Rdxd> Ajf{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const size_t& nb_faces = Ajf.numberOfSubValues(j);
+ const double& rhoc_j = rhoc[j];
+ for (size_t f = 0; f < nb_faces; ++f) {
+ Ajf(j, f) = tensorProduct(rhoc_j * Cjf(j, f), njf(j, f));
+ }
+ });
+
+ return Ajf;
+ }
+
+
+
+ NodeValuePerCell<const Rdxd>
+ _computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+ const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+ parallel_for(
+ Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ const double& rho_c = rhoc[j];
+
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
+
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
+ }
+ }
+ return std::numeric_limits<size_t>::max();
+ };
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+ }
+ }
+ });
+
+ return Ajr;
+ }
+
+ EdgeValuePerCell<const Rdxd>
+ _computeEucclhydCompositeAje(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ /*MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const EdgeValuePerFace<const Rd> Nle = mesh_data.Nle();
+ const EdgeValuePerFace<const Rd> nle = mesh_data.nle();
+
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+*/
+ EdgeValuePerCell<Rdxd> Aje{mesh.connectivity()};
+
+ parallel_for(
+ Aje.numberOfValues(), PUGS_LAMBDA(size_t jr) { Aje[jr] = zero; });
+ /*
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ const double& rho_c = rhoc[j];
+
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_edge_matrix[l];
+
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
+ }
+ }
+ return std::numeric_limits<size_t>::max();
+ };
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Aje(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+ }
+ }
+ });
+ */
+ return Aje;
+ }
+
+ FaceValuePerCell<const Rdxd>
+ _computeEucclhydCompositeAjf(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ /*MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const FaceValue<const Rd> Nlf = mesh_data.Nlf();
+ const FaceValue<const Rd> nlf = mesh_data.nlf();
+
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+*/
+ FaceValuePerCell<Rdxd> Ajf{mesh.connectivity()};
+
+ parallel_for(
+ Ajf.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajf[jr] = zero; });
+ /*
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ const double& rho_c = rhoc[j];
+
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_edge_matrix[l];
+
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
+ }
+ }
+ return std::numeric_limits<size_t>::max();
+ };
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Aje(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+ }
+ }
+ });
+ */
+ return Ajf;
+ }
+
+
+
+ NodeValuePerCell<const Rdxd>
+ _computeAjr(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ if constexpr (Dimension == 1) {
+ return _computeGlaceAjr(mesh, rhoc);
+ } else {
+ switch (solver_type) {
+ case SolverType::GlaceComposite: {
+ return _computeGlaceAjr(mesh, rhoc);
+ }
+ case SolverType::EucclhydComposite: {
+ return _computeEucclhydAjr(mesh, rhoc);
+ }
+
+ default: {
+ throw UnexpectedError("invalid solver type");
+ }
+ }
+ }
+ }
+
+ EdgeValuePerCell<const Rdxd>
+ _computeAje(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ if constexpr (Dimension == 1) {
+ throw UnexpectedError("invalid solver type pour arete 1D..");
+ } else {
+ switch (solver_type) {
+ case SolverType::GlaceComposite: {
+ return _computeGlaceAje(mesh, rhoc);
+ }
+ case SolverType::EucclhydComposite: {
+ //return _computeEucclhydAje(mesh, rhoc);
+ }
+ default: {
+ throw UnexpectedError("invalid solver type");
+ }
+ }
+ }
+ }
+
+ FaceValuePerCell<const Rdxd>
+ _computeAjf(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
+ {
+ if constexpr (Dimension == 1) {
+ throw UnexpectedError("invalid solver type pour face 1D..");
+ } else {
+ switch (solver_type) {
+ case SolverType::GlaceComposite: {
+ return _computeGlaceAjf(mesh, rhoc);
+ }
+ case SolverType::EucclhydComposite: {
+ //return _computeEucclhydCompositeAjf(mesh, rhoc);
+ }
+ default: {
+ throw UnexpectedError("invalid solver type");
+ }
+ }
+ }
+ }
+
+
+ NodeValue<Rdxd>
+ _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();
+
+ NodeValue<Rdxd> Ar{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ Rdxd sum = zero;
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ sum += Ajr(J, R);
+ }
+ Ar[r] = sum;
+ });
+
+ return Ar;
+ }
+
+ NodeValue<Rd>
+ _computeBr(const Mesh<Dimension>& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rd> b{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ Rd br = zero;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R);
+ }
+
+ b[r] = br;
+ });
+
+ return b;
+ }
+
+ 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: {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fixed: {
+ bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::external: {
+ if constexpr (Dimension == 1) {
+ const ExternalBoundaryConditionDescriptor& extern_bc_descriptor =
+ dynamic_cast<const ExternalBoundaryConditionDescriptor&>(*bc_descriptor);
+ bc_list.emplace_back(
+ ExternalFSIVelocityBoundaryCondition(mesh, getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ extern_bc_descriptor.socket()));
+ } else {
+ throw NotImplementedError("external FSI velocity not implemented for dimension > 1");
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "velocity") {
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+
+ Array<const Rd> value_list =
+ InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
+ mesh.xr(),
+ mesh_node_boundary.nodeList());
+
+ bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
+ } else if (dirichlet_bc_descriptor.name() == "pressure") {
+ const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension == 1) {
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ Array<const double> node_values =
+ InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
+ mesh_node_boundary.nodeList());
+ PressureBoundaryCondition p(mesh_node_boundary, node_values);
+ bc_list.emplace_back(p);
+ } else {
+ MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ Array<const double> face_values =
+ InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
+ mesh_face_boundary.faceList());
+ bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
+ }
+
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ 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 acoustic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ 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 BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+ {
+ 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) const
+ {
+ NodeValue<Rd> u{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
+
+ return u;
+ }
+
+ FaceValue<const Rd>
+ _computeUf(const MeshType& mesh, //const FaceValue<Rdxd>& Af, const FaceValue<Rd>& bf,
+ const NodeValue<Rd>& Ur,
+ const EdgeValue<Rd>& Ue,
+ const DiscreteScalarFunction& Z
+ ) const
+ {
+ FaceValue<Rd> u{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfFaces(), PUGS_LAMBDA(FaceId f) { u[f] = zero;});//inverse(Ar[r]) * br[r]; });
+
+ return u;
+ }
+
+ EdgeValue<const Rd>
+ _computeUe(const MeshType& mesh,
+ // const EdgeValue<Rdxd>& Ae, const EdgeValue<Rd>& be,
+ const NodeValue<Rd>& Ur,
+ const DiscreteScalarFunction& Z
+ ) const
+ {
+ EdgeValue<Rd> u{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfEdges(), PUGS_LAMBDA(EdgeId e) { u[e] = zero; });//inverse(Ar[r]) * br[r]; });
+
+ return u;
+ }
+
+
+ NodeValuePerCell<Rd>
+ _computeFjr(const MeshType& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const NodeValue<const Rd>& ur,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ NodeValuePerCell<Rd> F{mesh.connectivity()};
+ parallel_for(
+ 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) {
+ F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
+ }
+ });
+
+ return F;
+ }
+
+ FaceValuePerCell<Rd>
+ _computeFjf(const MeshType& mesh,
+ const FaceValuePerCell<const Rdxd>& Ajf,
+ const FaceValue<const Rd>& uf,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf();
+
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ FaceValuePerCell<Rd> F{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ for (size_t f = 0; f < cell_faces.size(); ++f) {
+ F(j, f) = Ajf(j, f) * (u[j] - uf[cell_faces[f]]) + p[j] * Cjf(j, f);
+ }
+ });
+
+ return F;
+ }
+
+ EdgeValuePerCell<Rd>
+ _computeFje(const MeshType& mesh,
+ const EdgeValuePerCell<const Rdxd>& Aje,
+ const EdgeValue<const Rd>& ue,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& p) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
+
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+
+ EdgeValuePerCell<Rd> F{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_edges = cell_to_edge_matrix[j];
+
+ for (size_t e = 0; e < cell_edges.size(); ++e) {
+ F(j, e) = Aje(j, e) * (u[j] - ue[cell_edges[e]]) + p[j] * Cje(j, e);
+ }
+ });
+
+ return F;
+ }
+
+
+ public:
+
+ std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>
+ >
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, c_v, u_v, p_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+ const MeshType& mesh = *mesh_v->get<MeshType>();
+ //if constexpr (Dimension ==1)
+ //{
+ // throw NormalError("acoustic solver composite inutile dimension 1");
+ //}
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
+ EdgeValuePerCell<const Rdxd> Aje = this->_computeAje(solver_type, mesh, rho * c);
+ FaceValuePerCell<const Rdxd> Ajf = this->_computeAjf(solver_type, mesh, rho * c);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
+
+ 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);
+
+ NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+
+ auto edge_to_node_connectivity = mesh.connectivity().edgeToNodeMatrix();
+
+ EdgeValue<Rd> ue{mesh.connectivity()};
+ ue.fill(zero);
+
+ parallel_for(mesh.numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id)
+ {
+ const auto& edge_nodes = edge_to_node_connectivity[edge_id];
+ const NodeId& node0_id= edge_nodes[0];
+ const NodeId& node1_id= edge_nodes[1];
+ ue[edge_id]=.5*(ur[node0_id]+ur[node1_id]);
+ });
+
+ //EdgeValue<const Rd> ue = this->_computeUe(mesh, ur, rho * c);
+ EdgeValuePerCell<const Rd> Fje = this->_computeFje(mesh, Aje, ue, u, p);
+
+ auto face_to_node_connectivity = mesh.connectivity().faceToNodeMatrix();
+ FaceValue<Rd> uf{mesh.connectivity()};
+ uf.fill(zero);
+
+ parallel_for(mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id)
+ {
+ const auto& face_nodes = face_to_node_connectivity[face_id];
+ Rd V=zero;
+ for (size_t i_node=0;i_node<face_nodes.size();++i_node)
+ {
+ const NodeId& node= face_nodes[i_node];
+
+ V+=ur[node];
+ }
+ V*=1./face_nodes.size();
+ uf[face_id]=V;
+ });
+
+
+ //FaceValue<const Rd> uf = this->_computeUf(mesh, ur, ue, rho * c) ;
+ FaceValuePerCell<const Rd> Fjf = this->_computeFjf(mesh, Ajf, uf, u, p);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr),
+
+ std::make_shared<const ItemValueVariant>(ue),
+ std::make_shared<const SubItemValuePerItemVariant>(Fje),
+
+ std::make_shared<const ItemValueVariant>(uf),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjf)
+ );
+ }
+
+
+
+ //2d
+ /*
+ std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>
+ >
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, c_v, u_v, p_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+ const MeshType& mesh = *mesh_v->get<MeshType>();
+ //if constexpr (Dimension ==1)
+ //{
+ // throw NormalError("acoustic solver composite inutile dimension 1");
+ //}
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
+ EdgeValuePerCell<const Rdxd> Aje = this->_computeAje(solver_type, mesh, rho * c);
+ FaceValuePerCell<const Rdxd> Ajf = this->_computeAjf(solver_type, mesh, rho * c);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
+
+ 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);
+
+ NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+
+ FaceValue<Rd> uf{mesh.connectivity()};
+ uf.fill(zero);
+
+ auto face_to_node_connectivity = mesh.connectivity().faceToNodeMatrix();
+
+ parallel_for(mesh.numberOfFaces(), PUGS_LAMBDA(const FaceId face_id)
+ {
+ const auto& face_nodes = face_to_node_connectivity[face_id];
+ Rd V=zero;
+ for (size_t i_node=0;i_node<face_nodes.size();++i_node)
+ {
+ const NodeId& node= face_nodes[i_node];
+
+ V+=ur[node];
+ }
+ V*=1./face_nodes.size();
+ uf[face_id]=V;
+ });
+
+
+ //FaceValue<const Rd> uf = this->_computeUf(mesh, ur, ue, rho * c) ;
+ FaceValuePerCell<const Rd> Fjf = this->_computeFjf(mesh, Ajf, uf, u, p);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr),
+
+ std::make_shared<const ItemValueVariant>(uf),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjf)
+ );
+ }
+
+*/
+
+
+
+ // 3d
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply_fluxes(const double& dt,
+ const MeshType& mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Fjr,
+ const EdgeValue<const Rd>& ue,
+ const EdgeValuePerCell<const Rd>& Fje,
+ const FaceValue<const Rd>& uf,
+ const FaceValuePerCell<const Rd>& Fjf
+ ) const
+ {
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjr.connectivity_ptr()) or
+ (mesh.shared_connectivity() != ue.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fje.connectivity_ptr()) or
+ (mesh.shared_connectivity() != uf.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjf.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 * ur[r]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ CellValue<double> new_rho = copy(rho.cellValues());
+ CellValue<Rd> new_u = copy(u.cellValues());
+ CellValue<double> new_E = copy(E.cellValues());
+
+ const double teta_e=1.;
+ const double teta_f=0.;
+ const double teta_r=1.0-teta_e-teta_f;
+
+ double tetar=teta_r;
+ double tetaf=teta_e;
+ double tetae=0;
+
+ if (Dimension==3)
+ {
+ tetaf=teta_f;
+ tetae=teta_e;
+ }
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ momentum_fluxes +=tetar * Fjr(j, R);
+ energy_fluxes +=tetar * dot(Fjr(j, R), ur[r]);
+ }
+
+ const auto& cell_edges = cell_to_edge_matrix[j];
+
+ for (size_t R = 0; R < cell_edges.size(); ++R) {
+ const EdgeId r = cell_edges[R];
+ momentum_fluxes +=tetae * Fje(j, R);
+ energy_fluxes +=tetae * dot(Fje(j, R), ue[r]);
+ }
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ for (size_t R = 0; R < cell_faces.size(); ++R) {
+ const FaceId r = cell_faces[R];
+ momentum_fluxes += tetaf * Fjf(j, R);
+ energy_fluxes += tetaf * dot(Fjf(j, R), uf[r]);
+ }
+
+ const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ new_u[j] -= dt_over_Mj * momentum_fluxes;
+ new_E[j] -= dt_over_Mj * energy_fluxes;
+ });
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+ return {std::make_shared<MeshVariant>(new_mesh),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
+ }
+
+
+ // 2d ..
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply_fluxes(const double& dt,
+ const MeshType& mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Fjr,
+ //const EdgeValue<const Rd>& ue,
+ //const EdgeValuePerCell<const Rd>& Fje,
+ const FaceValue<const Rd>& uf,
+ const FaceValuePerCell<const Rd>& Fjf
+ ) const
+ {
+ throw NormalError(" Cas apply flux composite 2D spec not finished");
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ //const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjr.connectivity_ptr()) or
+ //(mesh.shared_connectivity() != ue.connectivity_ptr()) or
+ //(mesh.shared_connectivity() != Fje.connectivity_ptr()) or
+ (mesh.shared_connectivity() != uf.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjf.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 * ur[r]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+ CellValue<double> new_rho = copy(rho.cellValues());
+ CellValue<Rd> new_u = copy(u.cellValues());
+ CellValue<double> new_E = copy(E.cellValues());
+
+ const double teta_e=0.5;
+ const double teta_f=0.0;
+ const double teta_r=1.0-teta_e-teta_f;
+
+ double tetar=teta_r;
+ double tetaf=teta_e;
+ double tetae=0;
+
+ //if (Dimension==3)
+ //{
+ // tetaf=teta_f;
+ // tetae=teta_e;
+ //}
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ momentum_fluxes +=tetar * Fjr(j, R);
+ energy_fluxes +=tetar * dot(Fjr(j, R), ur[r]);
+ }
+
+ // const auto& cell_edges = cell_to_edge_matrix[j];
+
+ // for (size_t R = 0; R < cell_edges.size(); ++R) {
+ // const EdgeId r = cell_edges[R];
+ // momentum_fluxes +=tetae * Fje(j, R);
+ // energy_fluxes +=tetae * dot(Fje(j, R), ue[r]);
+ // }
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ for (size_t R = 0; R < cell_faces.size(); ++R) {
+ const FaceId r = cell_faces[R];
+ momentum_fluxes += tetaf * Fjf(j, R);
+ energy_fluxes += tetaf * dot(Fjf(j, R), uf[r]);
+ }
+
+ const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ new_u[j] -= dt_over_Mj * momentum_fluxes;
+ new_E[j] -= dt_over_Mj * energy_fluxes;
+ });
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+ return {std::make_shared<MeshVariant>(new_mesh),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
+ }
+
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply_fluxes(const double& dt,
+ 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 ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+ 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 mesh_v = getCommonMesh({rho_v, u_v, E_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 this->apply_fluxes(dt, //
+ *mesh_v->get<MeshType>(), //
+ rho_v->get<DiscreteScalarFunction>(), //
+ u_v->get<DiscreteVectorFunction>(), //
+ E_v->get<DiscreteScalarFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
+ Fjr->get<NodeValuePerCell<const Rd>>(),
+ ue->get<EdgeValue<const Rd>>(), //
+ Fje->get<EdgeValuePerCell<const Rd>>(),
+ uf->get<FaceValue<const Rd>>(), //
+ Fjf->get<FaceValuePerCell<const Rd>>()
+ );
+ }
+
+
+
+ // 2d
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply_fluxes(const double& dt,
+ 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 ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+ //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 mesh_v = getCommonMesh({rho_v, u_v, E_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 this->apply_fluxes(dt, //
+ *mesh_v->get<MeshType>(), //
+ rho_v->get<DiscreteScalarFunction>(), //
+ u_v->get<DiscreteVectorFunction>(), //
+ E_v->get<DiscreteScalarFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
+ Fjr->get<NodeValuePerCell<const Rd>>(),
+ //ue->get<EdgeValue<const Rd>>(), //
+ //Fje->get<EdgeValuePerCell<const Rd>>(),
+ uf->get<FaceValue<const Rd>>(), //
+ Fjf->get<FaceValuePerCell<const Rd>>()
+ );
+ }
+
+
+
+ // 3d
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt,
+ 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
+ {
+ auto [ur, Fjr, ue, Fje, uf, Fjf] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
+ return apply_fluxes(dt, rho, u, E, ur, Fjr, ue, Fje, uf, Fjf);
+ }
+
+ /* ci dessous marche pas .. faudrait passer template specialisation dim=2
+ //2d
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt,
+ 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
+ {
+ auto [ur, Fjr, uf, Fjf] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
+ return apply_fluxes(dt, rho, u, E, ur, Fjr, uf, Fjf);
+ }
+ */
+
+ AcousticCompositeSolver() = default;
+ AcousticCompositeSolver(AcousticCompositeSolver&&) = default;
+ ~AcousticCompositeSolver() = default;
+};
+
+template <MeshConcept MeshType>
+void
+AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& br) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ if constexpr (Dimension == 1) {
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ Assert(node_cell_list.size() == 1);
+
+ CellId node_cell_id = node_cell_list[0];
+ size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+ br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+ });
+ } else {
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ NodeId node_id = face_nodes[i_node];
+ br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+ const NodeId r = node_list[r_number];
+
+ Ar[r] = P * Ar[r] * P + nxn;
+ br[r] = P * br[r];
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+
+ Ar[node_id] = identity;
+ br[node_id] = value;
+ });
+ } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+
+ Ar[node_id] = identity;
+ br[node_id] = zero;
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
+ const DiscreteScalarFunction& p,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<ExternalFSIVelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList(p);
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+
+ Ar[node_id] = identity;
+ br[node_id] = value;
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::FixedBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+ ~FixedBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::PressureBoundaryCondition
+{
+ private:
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Array<const double> m_value_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
+ : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<Mesh<1>>::PressureBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const Array<const double> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
+{
+ private:
+ const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const Array<TinyVector<Dimension>> m_value_list;
+ const std::shared_ptr<const Socket> m_socket;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ Array<const TinyVector<Dimension>>
+ valueList(const DiscreteScalarFunction& p) const
+ {
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Parallelism is not supported");
+ }
+ if (m_value_list.size() != 1) {
+ throw NotImplementedError("Non connected boundaries are not supported");
+ }
+
+ const CellId cell_id = m_node_to_cell_matrix[m_mesh_node_boundary.nodeList()[0]][0];
+
+ write(*m_socket, p[cell_id]);
+ read(*m_socket, m_value_list[0]);
+ return m_value_list;
+ }
+
+ ExternalFSIVelocityBoundaryCondition(const Mesh<Dimension>& mesh,
+ const MeshNodeBoundary& mesh_node_boundary,
+ const std::shared_ptr<const Socket>& socket)
+ : m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
+ m_mesh_node_boundary{mesh_node_boundary},
+ m_value_list(mesh_node_boundary.nodeList().size()),
+ m_socket{socket}
+ {
+ if constexpr (Dimension > 1) {
+ throw NotImplementedError("external FSI velocity bc in dimension > 1");
+ }
+ }
+
+ ~ExternalFSIVelocityBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::VelocityBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ const Array<const TinyVector<Dimension>> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~VelocityBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class AcousticCompositeSolverHandler::AcousticCompositeSolver<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_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();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+AcousticCompositeSolverHandler::AcousticCompositeSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
+{
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ std::visit(
+ [&](auto&& mesh) {
+ using MeshType = mesh_type_t<decltype(mesh)>;
+ if constexpr ((is_polygonal_mesh_v<MeshType>) and (MeshType::Dimension>1))
+ {
+ m_acoustic_composite_solver = std::make_unique<AcousticCompositeSolver<MeshType>>();
+ } else {
+ if (MeshType::Dimension==1)
+ throw NormalError("Use classical AcousticSolver instead of Composite in 1D");
+ else
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/AcousticCompositeSolver.hpp b/src/scheme/AcousticCompositeSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..46fc548e9d26c77ef36c7e6efddec1f786359c55
--- /dev/null
+++ b/src/scheme/AcousticCompositeSolver.hpp
@@ -0,0 +1,99 @@
+#ifndef ACOUSTIC_COMPOSITE_SOLVER_HPP
+#define ACOUSTIC_COMPOSITE_SOLVER_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
+
+double composite_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c);
+
+class AcousticCompositeSolverHandler
+{
+ public:
+ enum class SolverType
+ {
+ GlaceComposite,
+ EucclhydComposite
+ };
+
+ private:
+ struct IAcousticCompositeSolver
+ {
+
+ virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>
+ >
+
+ compute_fluxes(
+ const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+
+ virtual std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply_fluxes(const double& dt,
+ 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 ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+ 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 = 0;
+
+
+ virtual std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt,
+ 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 = 0;
+
+ IAcousticCompositeSolver() = default;
+ IAcousticCompositeSolver(IAcousticCompositeSolver&&) = default;
+ IAcousticCompositeSolver& operator=(IAcousticCompositeSolver&&) = default;
+
+ virtual ~IAcousticCompositeSolver() = default;
+ };
+
+ template <MeshConcept MeshType>
+ class AcousticCompositeSolver;
+
+ std::unique_ptr<IAcousticCompositeSolver> m_acoustic_composite_solver;
+
+ public:
+ const IAcousticCompositeSolver&
+ solver() const
+ {
+ return *m_acoustic_composite_solver;
+ }
+
+ AcousticCompositeSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v);
+};
+
+#endif // ACOUSTIC_COMPOSITE_SOLVER_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index dadbe73bb24ad84e491248b89016baf85f0dbc31..8cdd43f9e51e3b36c0c93e4d4cbfcbdfd9325f71 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(
PugsScheme
AcousticSolver.cpp
+ AcousticCompositeSolver.cpp
HyperelasticSolver.cpp
DiscreteFunctionIntegrator.cpp
DiscreteFunctionInterpoler.cpp
diff --git a/tests/test_TestVolumes3D.cpp b/tests/test_TestVolumes3D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2bf95290526dbd3e6816f4cee2dd53a6640cb6bc
--- /dev/null
+++ b/tests/test_TestVolumes3D.cpp
@@ -0,0 +1,293 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <mesh/ItemArrayUtils.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+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
+
+ }
+ }
+}