diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 7597da325c343a0c8e3422c8f922688908c9ca13..05d090b86fd2ea8e0ad5f6f089dda57773a849d1 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -305,7 +305,7 @@ LocalFSIModule::LocalFSIModule()
));
- this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -326,7 +326,59 @@ LocalFSIModule::LocalFSIModule()
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, const size_t& q)
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1,
+ const size_t& q)
+ -> 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, 1, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_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 std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1)
-> std::tuple<
std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -341,14 +393,14 @@ LocalFSIModule::LocalFSIModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1,
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
- bc_descriptor_list2, mu, lambda);
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
}
));
- this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -369,7 +421,59 @@ LocalFSIModule::LocalFSIModule()
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)
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1,
+ const size_t& q)
+ -> 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, 2, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_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 std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1)
-> std::tuple<
std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -384,9 +488,9 @@ LocalFSIModule::LocalFSIModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt1, rho1, rho2, u1,
u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
- bc_descriptor_list2, mu, lambda);
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
}
));
@@ -610,7 +714,7 @@ LocalFSIModule::LocalFSIModule()
));
- this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
+ this->_addBuiltinFunction("local_dt_hyperelastic_glace_neohook_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -631,7 +735,59 @@ LocalFSIModule::LocalFSIModule()
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, const size_t& q)
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1,
+ const size_t& q)
+ -> 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, 1, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("local_dt_hyperelastic_glace_neohook_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 std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1)
-> std::tuple<
std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -646,14 +802,14 @@ LocalFSIModule::LocalFSIModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1,
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, 1, dt1, rho1, rho2, u1,
u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
- bc_descriptor_list2, mu, lambda);
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
}
));
- this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
+ this->_addBuiltinFunction("local_dt_hyperelastic_glace_neohook2_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -674,7 +830,59 @@ LocalFSIModule::LocalFSIModule()
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)
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1,
+ const size_t& q)
+ -> 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, 2, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("local_dt_hyperelastic_glace_neohook2_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 std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf2,
+ const double& dt1)
-> std::tuple<
std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -689,9 +897,9 @@ LocalFSIModule::LocalFSIModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1,
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, 2, dt1, rho1, rho2, u1,
u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
- bc_descriptor_list2, mu, lambda);
+ bc_descriptor_list2, lambda2, mu2, gamma2, p_inf2);
}
));
diff --git a/src/scheme/LocalDtHyperelasticSolver.cpp b/src/scheme/LocalDtHyperelasticSolver.cpp
index f5ec4f12ab7a5214a05b6d36863987ba58214db8..2c3cb3cedfcfacb9281315601e8d272fc7381653 100644
--- a/src/scheme/LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/LocalDtHyperelasticSolver.cpp
@@ -19,6 +19,7 @@
#include <scheme/HyperelasticSolver.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/LocalDtUtils.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
@@ -860,6 +861,132 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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>>
+ compute_sub_iterations(const SolverType& solver_type,
+ const size_t& law,
+ const double& dt,
+ const size_t& q,
+ 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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ std::vector<NodeId> map2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
+ {
+ std::shared_ptr sub_m = getCommonMesh({rho, u, E, CG});
+ std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_u = u;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_E = E;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_CG = CG;
+
+ std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = lambda;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_mu = mu;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_gamma = gamma;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf = p_inf;
+
+ double sub_dt = dt / q;
+ double sum_dt = 0.;
+ for(size_t i = 0; i < q; ++i){
+
+ if( i == q - 1){
+ sub_dt = dt - sum_dt;
+ }
+
+ const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);
+
+ const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
+ const std::shared_ptr<const DiscreteFunctionVariant>& c =
+ std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+ std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_ur2, sub_Fjr2);
+
+ sub_lambda = shallowCopy(sub_m,sub_lambda);
+ sub_mu = shallowCopy(sub_m,sub_mu);
+ sub_gamma = shallowCopy(sub_m,sub_gamma);
+ sub_p_inf = shallowCopy(sub_m,sub_p_inf);
+
+ sum_dt += sub_dt;
+
+ }
+
+ return {sub_m, sub_rho, sub_u, sub_E, sub_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>>
+ compute_sub_iterations(const SolverType& solver_type,
+ const size_t& law,
+ const double& CFL,
+ 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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ std::vector<NodeId> map2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf) const
+ {
+ std::shared_ptr sub_m = getCommonMesh({rho, u, E, CG});
+ std::shared_ptr<const DiscreteFunctionVariant> sub_rho = rho;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_u = u;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_E = E;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_CG = CG;
+
+ std::shared_ptr<const DiscreteFunctionVariant> sub_lambda = lambda;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_mu = mu;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_gamma = gamma;
+ std::shared_ptr<const DiscreteFunctionVariant> sub_p_inf = p_inf;
+
+ double sum_dt = 0;
+ while(sum_dt < dt){
+ const auto& [sigma, p, aL, aT] = apply_state_law<MeshType>(law,sub_rho,sub_u,sub_E,sub_CG,sub_lambda,sub_mu,sub_gamma,sub_p_inf);
+
+ const DiscreteScalarFunction& aL_d = aL->template get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT_d = aT->template get<DiscreteScalarFunction>();
+ const std::shared_ptr<const DiscreteFunctionVariant>& c =
+ std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+ double sub_dt = CFL * hyperelastic_dt(c);
+ if(sum_dt + sub_dt > dt){
+ sub_dt = dt - sum_dt;
+ }
+
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+
+ std::tie(sub_m, sub_rho, sub_u, sub_E, sub_CG) = apply_fluxes(sub_dt, sub_rho, sub_u, sub_E, sub_CG, sub_ur2, sub_Fjr2);
+
+ sub_lambda = shallowCopy(sub_m,sub_lambda);
+ sub_mu = shallowCopy(sub_m,sub_mu);
+ sub_gamma = shallowCopy(sub_m,sub_gamma);
+ sub_p_inf = shallowCopy(sub_m,sub_p_inf);
+
+ sum_dt += sub_dt;
+ }
+
+ return {sub_m, sub_rho, sub_u, sub_E, sub_CG};
+ }
+
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -871,6 +998,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt1,
const size_t& q,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -889,66 +1017,22 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& pinf2) 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;
+ std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
- auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+ auto [map1, map2] = _computeMapping(m1, 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;
- }
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations(solver_type,law,dt1,q,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,lambda2,mu2,gamma2,pinf2);
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);
@@ -967,6 +1051,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
+ const size_t& law,
const double& dt1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -984,74 +1069,23 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& pinf2) 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);
+ std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
- auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+ auto [map1, map2] = _computeMapping(m1, 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);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
- sum_dt += dt2;
- }
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations(solver_type,law,0.4,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,lambda2,mu2,gamma2,pinf2);
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);
diff --git a/src/scheme/LocalDtHyperelasticSolver.hpp b/src/scheme/LocalDtHyperelasticSolver.hpp
index 6f2a80a318871a503d68bdf1c500b4aaa36bde82..5fc0aa0f54dabd8a43e2494adc3990889acf2275 100644
--- a/src/scheme/LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/LocalDtHyperelasticSolver.hpp
@@ -51,49 +51,53 @@ class LocalDtHyperelasticSolverHandler
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;
+ 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 size_t& law,
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& pinf2) 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>>
+ 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 size_t& law,
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,
@@ -110,8 +114,10 @@ class LocalDtHyperelasticSolverHandler
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;
+ const std::shared_ptr<const DiscreteFunctionVariant>& lambda2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& pinf2) const = 0;
ILocalDtHyperelasticSolver() = default;
ILocalDtHyperelasticSolver(ILocalDtHyperelasticSolver&&) = default;
diff --git a/src/scheme/Order2Limiters.hpp b/src/scheme/Order2Limiters.hpp
index 49ef8dcdfd7d579a1c635ef1d21394b08eaeb537..277e02d6b7116247c622850fd2bf50ae815b5f2d 100644
--- a/src/scheme/Order2Limiters.hpp
+++ b/src/scheme/Order2Limiters.hpp
@@ -154,7 +154,7 @@ scalar_limiter(const MeshType& mesh,
t[cell_id] = tj;
l[cell_id] = lj;
if(l2Norm(nj) != 0){
- n[cell_id] = (1/l2Norm(nj))*n[cell_id];
+ n[cell_id] = (1/l2Norm(nj))*n[cell_id];
}
if(l2Norm(tj) != 0){
t[cell_id] = (1/l2Norm(tj))*t[cell_id];