From 1bc2762644cc09c96e91a56ebf7388dc928ad345 Mon Sep 17 00:00:00 2001
From: Alexandre GANGLOFF <alexandre.gangloff@cea.fr>
Date: Thu, 28 Nov 2024 10:30:42 +0100
Subject: [PATCH] Add choice for the contructive law for hyper mono method
---
src/language/modules/SchemeModule.cpp | 50 +++++++++--
src/scheme/Order2HyperelasticSolver.cpp | 110 +++++++++++++++++++-----
src/scheme/Order2HyperelasticSolver.hpp | 9 +-
3 files changed, 141 insertions(+), 28 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 1d6583a48..33c1f7b50 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -823,7 +823,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("hyperelastic_eucclhyd_solver_order2",
+ this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -836,6 +836,8 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+ const double& lambda,
+ const double& mu,
const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -843,13 +845,13 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
- .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, CG, aL, aT,
- sigma, bc_descriptor_list, chi_solid);
+ .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
+ sigma, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
}
));
- this->_addBuiltinFunction("hyperelastic_glace_solver_order2",
+ this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -862,6 +864,10 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf,
const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -869,8 +875,40 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
- .apply(Order2HyperelasticSolverHandler::SolverType::Glace, dt, rho, u, E, CG, aL, aT,
- sigma, bc_descriptor_list, chi_solid);
+ .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
+ sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook2_solver_order2",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s,
+ const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
+ .solver()
+ .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt, rho, u, E, CG, aL, aT,
+ sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
}
));
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 9c2bd2591..c249ab60f 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -491,6 +491,86 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
return lambda;
}
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_NeoHook(const DiscreteScalarFunction rho,
+ const DiscreteVectorFunction u,
+ const DiscreteScalarFunction E,
+ const DiscreteTensorFunction CG,
+ const DiscreteScalarFunction chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) const
+ {
+ const TinyMatrix<Dimension> I = identity;
+
+ const DiscreteScalarFunction& eps = E - 0.5 * dot(u,u);
+ const DiscreteScalarFunction& p_d = (1-chi_solid)*(gamma-1) * rho * eps - p_inf*gamma;
+ const DiscreteTensorFunction& sigma_d = chi_solid * (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
+ const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
+
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
+ std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ }
+
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_NeoHook2(const DiscreteScalarFunction rho,
+ const DiscreteVectorFunction u,
+ const DiscreteScalarFunction E,
+ const DiscreteTensorFunction CG,
+ const DiscreteScalarFunction chi_solid,
+ const double& mu,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const
+ {
+ const TinyMatrix<Dimension> I = identity;
+
+ const DiscreteScalarFunction& eps = E - 0.5 * dot(u,u);
+ const DiscreteScalarFunction& p_d = (1-chi_solid)*(gamma_f-1) * rho * eps - p_inf_f*gamma_f + chi_solid * (gamma_s - 1) * rho * eps - gamma_s * p_inf_s;
+ const DiscreteTensorFunction& sigma_d = chi_solid * 2*mu/sqrt(det(CG))*(CG-1./3. * trace(CG)*I) - p_d * I;
+ const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (2 * mu) / rho + ((1-chi_solid)*gamma_f * p_d + chi_solid * (gamma_s*(p_d+p_inf_s))) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
+
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
+ std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ }
+
+ std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_state_law(const size_t law,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid_v,
+ const double& lambda,
+ const double& mu,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const
+ {
+ const DiscreteScalarFunction& chi_solid_d = chi_solid_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& rho_d = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u_d = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E_d = E_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& CG_d = CG_v->get<DiscreteTensorFunction>();
+
+ //const std::shared_ptr<const DiscreteFunctionVariant> sigma;
+ //const std::shared_ptr<const DiscreteFunctionVariant> aL;
+ //const std::shared_ptr<const DiscreteFunctionVariant> aT;
+
+ if(law == 1){
+ return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,mu,gamma_f,p_inf_f);
+ } else {
+ return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid_d,lambda,gamma_f,p_inf_f,gamma_s,p_inf_s);
+ }
+ }
+
public:
std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes(const SolverType& solver_type,
@@ -851,6 +931,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
@@ -860,7 +941,13 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid) const
+ const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const
{
std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma});
std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
@@ -871,28 +958,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final : public O
auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
std::tie(m_app,rho_app,u_app,E_app,CG_app) = apply_fluxes(dt/2, rho, u, E, CG, ur, Fjr);
- const double& gamma = 1.4;
- const double& youngs = 7000;
- const double& nu = 0;
- const double& mu = youngs/(2*(1-nu));
- const double& lambda = youngs * nu / ((1+nu) * (1 - 2 * nu));
- const TinyMatrix<Dimension> I = identity;
std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
- const DiscreteScalarFunction& chi_solid_d = chi_solid_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& rho_d = rho_app->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u_d = u_app->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& E_d = E_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteTensorFunction& CG_d = CG_app->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p_d = (1-chi_solid_d)*(gamma-1) * rho_d * eps;
- const DiscreteTensorFunction& sigma_d = chi_solid_d * (mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I) - p_d * I;
- const DiscreteScalarFunction& aL_d = sqrt(chi_solid_d * (lambda + 2 * mu) / rho_d + gamma * p_d / rho_d);
- const DiscreteScalarFunction& aT_d = sqrt(chi_solid_d * mu / rho_d);
-
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma_app = std::make_shared<const DiscreteFunctionVariant>(sigma_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& aL_app = std::make_shared<const DiscreteFunctionVariant>(aL_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& aT_app = std::make_shared<const DiscreteFunctionVariant>(aT_d);
+ auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
return order2_apply_fluxes(dt, rho, u, E, CG, rho_app, CG_app, ur_app, Fjr_app);
diff --git a/src/scheme/Order2HyperelasticSolver.hpp b/src/scheme/Order2HyperelasticSolver.hpp
index 960959df5..0051bd0d7 100644
--- a/src/scheme/Order2HyperelasticSolver.hpp
+++ b/src/scheme/Order2HyperelasticSolver.hpp
@@ -55,6 +55,7 @@ class Order2HyperelasticSolverHandler
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
@@ -64,7 +65,13 @@ class Order2HyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid) const = 0;
+ const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma_f,
+ const double& p_inf_f,
+ const double& gamma_s,
+ const double& p_inf_s) const = 0;
IOrder2HyperelasticSolver() = default;
IOrder2HyperelasticSolver(IOrder2HyperelasticSolver&&) = default;
--
GitLab