diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 1f51634718f8f5bda9b510ddb83ba78bf0e8ec71..cb85a06b4090113896b7ce60e5a146e6b6921eb3 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -815,6 +815,49 @@ SchemeModule::SchemeModule() )); + this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver", + std::function( + + [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1, + const std::shared_ptr<const DiscreteFunctionVariant>& u1, + const std::shared_ptr<const DiscreteFunctionVariant>& E1, + const std::shared_ptr<const DiscreteFunctionVariant>& CG1, + const std::shared_ptr<const DiscreteFunctionVariant>& aL1, + const std::shared_ptr<const DiscreteFunctionVariant>& aT1, + const std::shared_ptr<const DiscreteFunctionVariant>& sigma1, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list1, + const std::shared_ptr<const DiscreteFunctionVariant>& rho2, + const std::shared_ptr<const DiscreteFunctionVariant>& u2, + const std::shared_ptr<const DiscreteFunctionVariant>& E2, + const std::shared_ptr<const DiscreteFunctionVariant>& CG2, + const std::shared_ptr<const DiscreteFunctionVariant>& aL2, + const std::shared_ptr<const DiscreteFunctionVariant>& aT2, + const std::shared_ptr<const DiscreteFunctionVariant>& sigma2, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list2, + const double& mu, const double& lambda, const double& dt1) + -> 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>> { + return LocalDtHyperelasticSolverHandler{getCommonMesh( + {rho1, u1, E1, CG1, aL1, aT1, sigma1}), + getCommonMesh( + {rho2, u2, E2, CG2, aL2, aT2, sigma2})} + .solver() + .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, + u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1, + bc_descriptor_list2, mu, lambda); + } + + )); + this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver", std::function( @@ -858,6 +901,49 @@ SchemeModule::SchemeModule() )); + this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver", + std::function( + + [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1, + const std::shared_ptr<const DiscreteFunctionVariant>& u1, + const std::shared_ptr<const DiscreteFunctionVariant>& E1, + const std::shared_ptr<const DiscreteFunctionVariant>& CG1, + const std::shared_ptr<const DiscreteFunctionVariant>& aL1, + const std::shared_ptr<const DiscreteFunctionVariant>& aT1, + const std::shared_ptr<const DiscreteFunctionVariant>& sigma1, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list1, + const std::shared_ptr<const DiscreteFunctionVariant>& rho2, + const std::shared_ptr<const DiscreteFunctionVariant>& u2, + const std::shared_ptr<const DiscreteFunctionVariant>& E2, + const std::shared_ptr<const DiscreteFunctionVariant>& CG2, + const std::shared_ptr<const DiscreteFunctionVariant>& aL2, + const std::shared_ptr<const DiscreteFunctionVariant>& aT2, + const std::shared_ptr<const DiscreteFunctionVariant>& sigma2, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list2, + const double& mu, const double& lambda, const double& dt1) + -> 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>> { + return LocalDtHyperelasticSolverHandler{getCommonMesh( + {rho1, u1, E1, CG1, aL1, aT1, sigma1}), + getCommonMesh( + {rho2, u2, E2, CG2, aL2, aT2, sigma2})} + .solver() + .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1, + u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1, + bc_descriptor_list2, mu, lambda); + } + + )); + this->_addBuiltinFunction("apply_hyperelastic_fluxes", std::function( diff --git a/src/scheme/LocalDtHyperelasticSolver.cpp b/src/scheme/LocalDtHyperelasticSolver.cpp index 0bfe1e06c868ce28554d8cd4e9e2128dc5edaf4f..a69082fba1164e05d7a35bc85fd7cd9d4190fe61 100644 --- a/src/scheme/LocalDtHyperelasticSolver.cpp +++ b/src/scheme/LocalDtHyperelasticSolver.cpp @@ -16,6 +16,7 @@ #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/SymmetryBoundaryConditionDescriptor.hpp> @@ -328,6 +329,12 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final 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, @@ -608,6 +615,91 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final 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>, @@ -864,6 +956,110 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final 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; + aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d); + 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>(aL_d + aT_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); + + 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}; + } + LocalDtHyperelasticSolver() = default; LocalDtHyperelasticSolver(LocalDtHyperelasticSolver&&) = default; ~LocalDtHyperelasticSolver() = default; @@ -1095,6 +1291,30 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCou } } +template <MeshConcept MeshType> +void +LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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 LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC( diff --git a/src/scheme/LocalDtHyperelasticSolver.hpp b/src/scheme/LocalDtHyperelasticSolver.hpp index 731e0810a4d11a233df73fc3b2ee24af429e9699..6f2a80a318871a503d68bdf1c500b4aaa36bde82 100644 --- a/src/scheme/LocalDtHyperelasticSolver.hpp +++ b/src/scheme/LocalDtHyperelasticSolver.hpp @@ -82,6 +82,37 @@ class LocalDtHyperelasticSolverHandler 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; + ILocalDtHyperelasticSolver() = default; ILocalDtHyperelasticSolver(ILocalDtHyperelasticSolver&&) = default; ILocalDtHyperelasticSolver& operator=(ILocalDtHyperelasticSolver&&) = default;