From 1579367a6d1eb4548beec45b93f3201f8299d410 Mon Sep 17 00:00:00 2001
From: Alexandre GANGLOFF <alexandre.gangloff@cea.fr>
Date: Fri, 15 Nov 2024 15:39:22 +0100
Subject: [PATCH] Create LocalDtHyperelasticSolver files
---
src/scheme/CMakeLists.txt | 1 +
.../Order2LocalDtHyperelasticSolver.cpp | 1605 +++++++++++++++++
.../Order2LocalDtHyperelasticSolver.hpp | 139 ++
3 files changed, 1745 insertions(+)
create mode 100644 src/scheme/Order2LocalDtHyperelasticSolver.cpp
create mode 100644 src/scheme/Order2LocalDtHyperelasticSolver.hpp
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index e63bff35e..aaaaafc39 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -16,6 +16,7 @@ add_library(
PolynomialReconstruction.cpp
Order2AcousticSolver.cpp
Order2HyperelasticSolver.cpp
+ Order2LocalDtHyperelasticSolver.cpp
)
target_link_libraries(
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
new file mode 100644
index 000000000..e809d74be
--- /dev/null
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -0,0 +1,1605 @@
+#include <scheme/Order2LocalDtHyperelasticSolver.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/StencilArray.hpp>
+#include <mesh/StencilManager.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/HyperelasticSolver.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+#include <variant>
+#include <vector>
+
+template <MeshConcept MeshType>
+class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver final
+ : public Order2LocalDtHyperelasticSolverHandler::IOrder2LocalDtHyperelasticSolver
+{
+ 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>;
+ using DiscreteTensorFunction = DiscreteFunctionP0<const Rdxd>;
+
+ class FixedBoundaryCondition;
+ class PressureBoundaryCondition;
+ class NormalStressBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class VelocityBoundaryCondition;
+ class CouplingBoundaryCondition;
+
+ using BoundaryCondition = std::variant<FixedBoundaryCondition,
+ PressureBoundaryCondition,
+ NormalStressBoundaryCondition,
+ SymmetryBoundaryCondition,
+ VelocityBoundaryCondition,
+ CouplingBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ NodeValuePerCell<const Rdxd>
+ _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) 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()};
+ const Rdxd I = identity;
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const size_t& nb_nodes = Ajr.numberOfSubValues(j);
+ const double& rhoaL_j = rhoaL[j];
+ const double& rhoaT_j = rhoaT[j];
+ for (size_t r = 0; r < nb_nodes; ++r) {
+ const Rdxd& M = tensorProduct(Cjr(j, r), njr(j, r));
+ Ajr(j, r) = rhoaL_j * M + rhoaT_j * (l2Norm(Cjr(j, r)) * I - M);
+ }
+ });
+
+ return Ajr;
+ }
+
+ NodeValuePerCell<const Rdxd>
+ _computeEucclhydAjr(const MeshType& mesh,
+ const DiscreteScalarFunction& rhoaL,
+ const DiscreteScalarFunction& rhoaT) 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; });
+ const Rdxd I = identity;
+ 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_aL = rhoaL[j];
+ const double& rho_aT = rhoaT[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]);
+ const Rdxd& M = tensorProduct(Nlr(l, rl), nlr(l, rl));
+ Ajr(j, R) += rho_aL * M + rho_aT * (l2Norm(Nlr(l, rl)) * I - M);
+ }
+ }
+ });
+
+ return Ajr;
+ }
+
+ NodeValuePerCell<const Rdxd>
+ _computeAjr(const SolverType& solver_type,
+ const MeshType& mesh,
+ const DiscreteScalarFunction& rhoaL,
+ const DiscreteScalarFunction& rhoaT) const
+ {
+ if constexpr (Dimension == 1) {
+ return _computeGlaceAjr(mesh, rhoaL, rhoaT);
+ } else {
+ switch (solver_type) {
+ case SolverType::Glace: {
+ return _computeGlaceAjr(mesh, rhoaL, rhoaT);
+ }
+ case SolverType::Eucclhyd: {
+ return _computeEucclhydAjr(mesh, rhoaL, rhoaT);
+ }
+ 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 DiscreteTensorFunction& sigma) 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] - sigma[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::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());
+
+ bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
+ } 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 if (dirichlet_bc_descriptor.name() == "normal-stress") {
+ const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension == 1) {
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ Array<const Rdxd> node_values =
+ InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
+ mesh_node_boundary.nodeList());
+
+ bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
+ } else {
+ MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ Array<const Rdxd> face_values =
+ InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::face>(normal_stress_id, mesh_data.xl(),
+ mesh_face_boundary.faceList());
+ bc_list.emplace_back(NormalStressBoundaryCondition{mesh_face_boundary, face_values});
+ }
+
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::coupling: {
+ const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
+ dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
+ bc_list.emplace_back(CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ coupling_bc_descriptor.label()));
+ 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 hyperelastic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
+ void _applyNormalStressBC(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 _applyCouplingBC(NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const;
+ void _applyCouplingBC(const MeshType& mesh,
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr,
+ const std::vector<NodeId>& map2) const;
+ void _applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& ur1,
+ NodeValue<Rd>& ur2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const;
+ void
+ _applyBoundaryConditions(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+ {
+ this->_applyPressureBC(bc_list, mesh, br);
+ this->_applyNormalStressBC(bc_list, mesh, br);
+ this->_applySymmetryBC(bc_list, Ar, br);
+ this->_applyVelocityBC(bc_list, Ar, br);
+ }
+
+ NodeValue<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;
+ }
+
+ NodeValuePerCell<Rd>
+ _computeFjr(const MeshType& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const NodeValue<const Rd>& ur,
+ const DiscreteVectorFunction& u,
+ const DiscreteTensorFunction& sigma) 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]]) + sigma[j] * Cjr(j, r);
+ }
+ });
+
+ return F;
+ }
+
+ std::tuple<std::vector<NodeId>, std::vector<NodeId>>
+ _computeMapping(std::shared_ptr<const MeshVariant>& i_mesh1,
+ std::shared_ptr<const MeshVariant>& i_mesh2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ {
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
+
+ const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+
+ std::vector<NodeId> map1;
+ std::vector<NodeId> map2;
+
+ NodeValue<Rd> xr1 = copy(mesh1.xr());
+ NodeValue<Rd> xr2 = copy(mesh2.xr());
+ for (const auto& boundary_condition1 : bc_list1) {
+ std::visit(
+ [&](auto&& bc1) {
+ using T1 = std::decay_t<decltype(bc1)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T1>) {
+ const auto& node_list1 = bc1.nodeList();
+ for (const auto& boundary_condition2 : bc_list2) {
+ std::visit(
+ [&](auto&& bc2) {
+ using T2 = std::decay_t<decltype(bc2)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T2>) {
+ if (bc1.label() == bc2.label()) {
+ const auto& node_list2 = bc2.nodeList();
+ for (size_t i = 0; i < node_list1.size(); i++) {
+ NodeId node_id1 = node_list1[i];
+ NodeId node_id2;
+ for (size_t j = 0; j < node_list2.size(); j++) {
+ node_id2 = node_list2[j];
+ double err = 0;
+ err = dot((xr1[node_id1] - xr2[node_id2]), (xr1[node_id1] - xr2[node_id2]));
+ if (sqrt(err) < 1e-14) {
+ map1.push_back(node_id1);
+ map2.push_back(node_id2);
+ }
+ }
+ }
+ }
+ }
+ },
+ boundary_condition2);
+ }
+ }
+ },
+ boundary_condition1);
+ }
+ return std::make_tuple(map1, map2);
+ }
+
+ public:
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+
+ const MeshType& mesh = *i_mesh->get<MeshType>();
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma);
+
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, 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, sigma);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+ }
+
+ 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>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
+ {
+ std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
+ std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
+ if (not i_mesh1) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+ if (not i_mesh2) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+ const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL1 = aL1_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT1 = aT1_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
+ const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL2 = aL2_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
+ NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
+
+ NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+ NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1);
+ NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+ NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+ this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
+
+ const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+ this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+ this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+
+ synchronize(Ar1);
+ synchronize(br1);
+ synchronize(Ar2);
+ synchronize(br2);
+
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+ NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
+ }
+
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ const std::vector<NodeId> map2) const
+ {
+ std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+
+ const MeshType& mesh = *i_mesh->get<MeshType>();
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma);
+
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+
+ synchronize(Ar);
+ synchronize(br);
+
+ NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+
+ this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+ }
+
+ 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>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
+ compute_fluxes2(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
+ {
+ std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
+ std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
+ if (not i_mesh1) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+ if (not i_mesh2) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+ const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL1 = aL1_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT1 = aT1_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
+ const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL2 = aL2_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
+ NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
+
+ NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+ NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1);
+ NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+ NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+
+ const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+ this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+ this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+
+ synchronize(Ar1);
+ synchronize(br1);
+ synchronize(Ar2);
+ synchronize(br2);
+
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
+ ur2,
+ Fjr2);
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ 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 DiscreteTensorFunction& CG,
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Fjr) const
+ {
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+ (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+ throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+ }
+
+ NodeValue<Rd> new_xr = copy(mesh.xr());
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * 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();
+ const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+ CellValue<double> new_rho = copy(rho.cellValues());
+ CellValue<Rd> new_u = copy(u.cellValues());
+ CellValue<double> new_E = copy(E.cellValues());
+ CellValue<Rdxd> new_CG = copy(CG.cellValues());
+
+ 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;
+ Rdxd gradv = zero;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ gradv += tensorProduct(ur[r], Cjr(j, R));
+ momentum_fluxes += Fjr(j, R);
+ energy_fluxes += dot(Fjr(j, R), ur[r]);
+ }
+ const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+ const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ const double dt_over_Vj = dt / Vj[j];
+ new_u[j] += dt_over_Mj * momentum_fluxes;
+ new_E[j] += dt_over_Mj * energy_fluxes;
+ new_CG[j] += dt_over_Vj * cauchy_green_fluxes;
+ new_CG[j] += transpose(new_CG[j]);
+ new_CG[j] *= 0.5;
+ });
+
+ 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::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ 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 DiscreteFunctionVariant>& CG,
+ const std::shared_ptr<const ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+ {
+ std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+
+ return this->apply_fluxes(dt, //
+ *i_mesh->get<MeshType>(), //
+ rho->get<DiscreteScalarFunction>(), //
+ u->get<DiscreteVectorFunction>(), //
+ E->get<DiscreteScalarFunction>(), //
+ CG->get<DiscreteTensorFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
+ Fjr->get<NodeValuePerCell<const Rd>>());
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ mesh_correction(const MeshType& mesh1,
+ const MeshType& mesh2,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ const DiscreteTensorFunction& CG,
+ std::vector<NodeId>& map1,
+ std::vector<NodeId>& map2) const
+ {
+ NodeValue<Rd> xr1 = copy(mesh1.xr());
+ NodeValue<Rd> new_xr2 = copy(mesh2.xr());
+
+ size_t n = map1.size();
+
+ for (size_t i = 0; i < n; i++) {
+ for (size_t j = 0; j < Dimension; j++) {
+ new_xr2[map2[i]][j] = xr1[map1[i]][j];
+ }
+ }
+
+ std::shared_ptr<const MeshType> new_mesh2 = std::make_shared<MeshType>(mesh2.shared_connectivity(), new_xr2);
+
+ CellValue<double> new_rho = copy(rho.cellValues());
+ CellValue<Rd> new_u = copy(u.cellValues());
+ CellValue<double> new_E = copy(E.cellValues());
+ CellValue<Rdxd> new_CG = copy(CG.cellValues());
+
+ return {std::make_shared<MeshVariant>(new_mesh2),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh2, new_u)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_E)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh2, new_CG))};
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ mesh_correction(const std::shared_ptr<const MeshVariant>& i_mesh1,
+ const std::shared_ptr<const MeshVariant>& i_mesh2,
+ 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>& CG,
+ std::vector<NodeId>& map1,
+ std::vector<NodeId>& map2) const
+ {
+ return this->mesh_correction(*i_mesh1->get<MeshType>(), *i_mesh2->get<MeshType>(),
+ rho->get<DiscreteScalarFunction>(), u->get<DiscreteVectorFunction>(),
+ E->get<DiscreteScalarFunction>(), CG->get<DiscreteTensorFunction>(), map1, map2);
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt1,
+ const size_t& q,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const double& mu,
+ const double& lambda) const
+ {
+ std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+
+ std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_CG2 = CG2;
+
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
+ compute_fluxes(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
+ bc_descriptor_list2, map1, map2);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
+
+ double dt2 = dt1 / q;
+ const double gamma = 1.4;
+ const TinyMatrix<Dimension> I = identity;
+ double sum_dt = 0;
+
+ size_t fluid = 0;
+ if ((mu == 0) and (lambda == 0)) {
+ fluid = 1;
+ }
+
+ for (size_t i = 0; i < q; i++) {
+ if (i == q - 1) {
+ dt2 = dt1 - sum_dt;
+ }
+
+ const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteTensorFunction& CG_d = new_CG2->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p = fluid * (gamma - 1) * rho_d * eps;
+ const DiscreteTensorFunction& sigma_d =
+ mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
+ const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d);
+
+ const DiscreteScalarFunction& aT_d = sqrt(mu / rho_d);
+
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 =
+ std::make_shared<const DiscreteFunctionVariant>(sigma_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 =
+ std::make_shared<const DiscreteFunctionVariant>(aL_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 =
+ std::make_shared<const DiscreteFunctionVariant>(aT_d);
+
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, new_rho2, new_aL2, new_aT2, new_u2, new_sigma2,
+ bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ apply_fluxes(dt2, new_rho2, new_u2, new_E2, new_CG2, sub_ur2, sub_Fjr2);
+
+ sum_dt += dt2;
+ }
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
+
+ return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const double& mu,
+ const double& lambda) const
+ {
+ std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+
+ std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_CG2 = CG2;
+
+ DiscreteScalarFunction aL_d = aL2->get<DiscreteScalarFunction>();
+ DiscreteScalarFunction aT_d = aT2->get<DiscreteScalarFunction>();
+ const std::shared_ptr<const DiscreteFunctionVariant>& c = std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+ double dt2 = 0.4 * hyperelastic_dt(c);
+
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
+ compute_fluxes2(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
+ bc_descriptor_list2, map1, map2);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
+ std::tie(new_m2,new_rho2,new_u2,new_E2,new_CG2) = apply_fluxes(dt2, rho2, u2, E2, CG2, ur2, Fjr2);
+
+ const double gamma = 1.4;
+ const TinyMatrix<Dimension> I = identity;
+ double sum_dt = dt2;
+
+ size_t fluid = 0;
+ if ((mu == 0) and (lambda == 0)) {
+ fluid = 1;
+ }
+
+ while(sum_dt < dt1) {
+
+ const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteTensorFunction& CG_d = new_CG2->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p = fluid * (gamma - 1) * rho_d * eps;
+ const DiscreteTensorFunction& sigma_d = mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
+ const DiscreteScalarFunction& new_aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d);
+ const DiscreteScalarFunction& new_aT_d = sqrt(mu / rho_d);
+
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 =
+ std::make_shared<const DiscreteFunctionVariant>(sigma_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_c =
+ std::make_shared<const DiscreteFunctionVariant>(new_aL_d + new_aT_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 =
+ std::make_shared<const DiscreteFunctionVariant>(new_aL_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 =
+ std::make_shared<const DiscreteFunctionVariant>(new_aT_d);
+
+ dt2 = 0.4 * hyperelastic_dt(new_c);
+ if(sum_dt + dt2 > dt1){
+ dt2 = dt1 - sum_dt;
+ }
+
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, new_rho2, new_aL2, new_aT2, new_u2, new_sigma2,
+ bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ apply_fluxes(dt2, new_rho2, new_u2, new_E2, new_CG2, sub_ur2, sub_Fjr2);
+
+ sum_dt += dt2;
+ }
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
+
+ return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
+ }
+
+ Order2LocalDtHyperelasticSolver() = default;
+ Order2LocalDtHyperelasticSolver(Order2LocalDtHyperelasticSolver&&) = default;
+ ~Order2LocalDtHyperelasticSolver() = default;
+};
+
+template <MeshConcept MeshType>
+void
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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>) {
+ MeshDataType& 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
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyNormalStressBC(
+ 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<NormalStressBoundaryCondition, T>) {
+ MeshDataType& 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
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
+ NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
+{
+ size_t n = map1.size();
+
+ for (size_t i = 0; i < n; i++) {
+ NodeId node_id1 = map1[i];
+ NodeId node_id2 = map2[i];
+
+ Ar1[node_id1] += Ar2[node_id2];
+ Ar2[node_id2] = Ar1[node_id1];
+
+ br1[node_id1] += br2[node_id2];
+ br2[node_id2] = br1[node_id1];
+ }
+}
+
+template <MeshConcept MeshType>
+void
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& ur1,
+ NodeValue<Rd>& ur2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
+{
+ size_t n = map1.size();
+ Rd lambda;
+
+ for(size_t i = 0; i<n; i++){
+
+ NodeId node_id1 = map1[i];
+ NodeId node_id2 = map2[i];
+
+ lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1]))*(ur2[node_id2]-ur1[node_id1]);
+
+ ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1])*lambda;
+ ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2])*lambda;
+ }
+}
+
+template <MeshConcept MeshType>
+void
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
+ const MeshType& mesh,
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr,
+ const std::vector<NodeId>& map2) const
+{
+ const size_t& n = map2.size();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ for (size_t i = 0; i < n; i++) {
+ const NodeId node_id = map2[i];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ 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];
+ Fjr(J, R) = CR_Fjr(J, R);
+ }
+ ur[node_id] = CR_ur[node_id];
+ }
+}
+
+template <MeshConcept MeshType>
+class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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 Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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 Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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 Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::NormalStressBoundaryCondition
+{
+ private:
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Array<const Rdxd> m_value_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Array<const Rdxd>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NormalStressBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const Rdxd>& value_list)
+ : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+ {}
+
+ ~NormalStressBoundaryCondition() = default;
+};
+
+template <>
+class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<Mesh<1>>::NormalStressBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const Array<const Rdxd> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const Rdxd>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NormalStressBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const Rdxd>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~NormalStressBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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 Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<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;
+};
+
+template <MeshConcept MeshType>
+class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver<MeshType>::CouplingBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const size_t m_label;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ size_t
+ label() const
+ {
+ return m_label;
+ }
+
+ CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const size_t label)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_label{label}
+ {
+ ;
+ }
+
+ ~CouplingBoundaryCondition() = default;
+};
+
+Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& i_mesh1,
+ const std::shared_ptr<const MeshVariant>& i_mesh2)
+{
+ if (not i_mesh1) {
+ throw NormalError("mesh1 not defined");
+ }
+
+ if (not i_mesh2) {
+ throw NormalError("mesh2 not defined");
+ }
+ std::visit(
+ [&](auto&& first_mesh, auto&& second_mesh) {
+ using FirstMeshType = mesh_type_t<decltype(first_mesh)>;
+ using SecondMeshType = mesh_type_t<decltype(second_mesh)>;
+
+ if constexpr (!std::is_same_v<FirstMeshType,
+ SecondMeshType>) { // <- ICI POUR LE TEST SUR LES TYPES DE MAILLAGES
+ throw NormalError("incompatible mesh types");
+ }
+ },
+ i_mesh1->variant(), i_mesh2->variant());
+
+ std::visit(
+ [&](auto&& mesh) {
+ using MeshType = mesh_type_t<decltype(mesh)>;
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ m_hyperelastic_solver = std::make_unique<Order2LocalDtHyperelasticSolver<MeshType>>();
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ i_mesh1->variant());
+}
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
new file mode 100644
index 000000000..306ed1b60
--- /dev/null
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -0,0 +1,139 @@
+#ifndef LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+#define LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
+class DiscreteFunctionVariant;
+
+class Order2LocalDtHyperelasticSolverHandler
+{
+ public:
+ enum class SolverType
+ {
+ Glace,
+ Eucclhyd
+ };
+
+ private:
+ struct IOrder2LocalDtHyperelasticSolver
+ {
+ virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(
+ const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ 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>,
+ 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 DiscreteFunctionVariant>& CG,
+ const std::shared_ptr<const ItemValueVariant>& ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) 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>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt1,
+ const size_t& q,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const double& mu,
+ const double& lambda) 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>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const double& mu,
+ const double& lambda) const = 0;
+
+ IOrder2LocalDtHyperelasticSolver() = default;
+ IOrder2LocalDtHyperelasticSolver(IOrder2LocalDtHyperelasticSolver&&) = default;
+ IOrder2LocalDtHyperelasticSolver& operator=(IOrder2LocalDtHyperelasticSolver&&) = default;
+
+ virtual ~IOrder2LocalDtHyperelasticSolver() = default;
+ };
+
+ template <MeshConcept MeshType>
+ class Order2LocalDtHyperelasticSolver;
+
+ std::unique_ptr<IOrder2LocalDtHyperelasticSolver> m_hyperelastic_solver;
+
+ public:
+ const IOrder2LocalDtHyperelasticSolver&
+ solver() const
+ {
+ return *m_hyperelastic_solver;
+ }
+
+ Order2LocalDtHyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh1,
+ const std::shared_ptr<const MeshVariant>& mesh2);
+};
+
+#endif // HYPERELASTIC_SOLVER_HPP
--
GitLab