From caab12af3dcb28bda5e81a91bf5bbb55ff60d6d2 Mon Sep 17 00:00:00 2001
From: Alexandre GANGLOFF <alexandre.gangloff@cea.fr>
Date: Wed, 20 Nov 2024 17:22:12 +0100
Subject: [PATCH] Add order 2 extension in apply on LocalDt
---
src/language/modules/SchemeModule.cpp | 44 ++++
.../Order2LocalDtHyperelasticSolver.cpp | 243 ++++++++++++++----
.../Order2LocalDtHyperelasticSolver.hpp | 4 +-
3 files changed, 240 insertions(+), 51 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 07271feaa..1d6583a48 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -49,6 +49,7 @@
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/Order2AcousticSolver.hpp>
#include <scheme/Order2HyperelasticSolver.hpp>
+#include <scheme/Order2LocalDtHyperelasticSolver.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VariableBCDescriptor.hpp>
@@ -960,6 +961,49 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver_order2",
+ 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 Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+ {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+ getCommonMesh(
+ {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+ .solver()
+ .apply(Order2LocalDtHyperelasticSolverHandler::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(
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index 80f34ce82..d74835e0b 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -865,7 +865,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
{
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");
+ throw NormalError("discrete functions are not defined on the same mesh ici1");
}
if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
@@ -942,14 +942,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
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");
+ throw NormalError("discrete functions are not defined on the same mesh ici2");
}
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");
+ throw NormalError("discrete functions are not defined on the same mesh ici2 mesh2");
}
if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
@@ -1119,7 +1119,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
{
std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
if (not i_mesh) {
- throw NormalError("discrete functions are not defined on the same mesh");
+ throw NormalError("discrete functions are not defined on the same mesh ici3");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
@@ -1224,7 +1224,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
- throw NormalError("discrete functions are not defined on the same mesh");
+ throw NormalError("discrete functions are not defined on the same mesh ici4");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
@@ -1301,6 +1301,105 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
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_classic_midpoint_method(const SolverType& solver_type,
+ const double& dt1,
+ 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,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ std::vector<NodeId>& map2,
+ size_t fluid,
+ const double& lambda,
+ const double& mu) const
+ {
+
+ std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma});
+ std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
+ std::shared_ptr<const DiscreteFunctionVariant> u_app = u;
+ std::shared_ptr<const DiscreteFunctionVariant> E_app = E;
+ std::shared_ptr<const DiscreteFunctionVariant> CG_app = CG;
+
+ std::shared_ptr new_m = getCommonMesh({rho, aL, aT, u, sigma});
+ std::shared_ptr<const DiscreteFunctionVariant> new_rho = rho;
+ std::shared_ptr<const DiscreteFunctionVariant> new_u = u;
+ std::shared_ptr<const DiscreteFunctionVariant> new_E = E;
+ std::shared_ptr<const DiscreteFunctionVariant> new_CG = CG;
+ std::shared_ptr<const DiscreteFunctionVariant> new_sigma = sigma;
+ std::shared_ptr<const DiscreteFunctionVariant> new_aL = aL;
+ std::shared_ptr<const DiscreteFunctionVariant> new_aT = aT;
+
+ DiscreteScalarFunction new_aL_d = aL->get<DiscreteScalarFunction>();
+ DiscreteScalarFunction new_aT_d = aT->get<DiscreteScalarFunction>();
+
+ const TinyMatrix<Dimension> I = identity;
+ const double& gamma = 1.4;
+ double sum_dt = 0;
+
+ while(sum_dt < dt1){
+
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_c = std::make_shared<const DiscreteFunctionVariant>(new_aL_d + new_aT_d);
+
+ double dt2 = 0.4 * hyperelastic_dt(new_c);
+ if(sum_dt + dt2 > dt1){
+ dt2 = dt1 - sum_dt;
+ }
+
+ auto [ur, Fjr] = compute_fluxes(solver_type,new_rho,new_aL,new_aT,new_u,new_sigma,bc_descriptor_list,CR_ur,CR_Fjr,map2);
+ std::tie(m_app, rho_app, u_app, E_app, CG_app) = apply_fluxes(dt2/2,new_rho,new_u,new_E,new_CG,ur,Fjr);
+
+ 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 = 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_d * I;
+ const DiscreteScalarFunction& aL_d_app = sqrt((lambda + 2 * mu) / rho_d + gamma * p_d / rho_d);
+ const DiscreteScalarFunction& aT_d_app = sqrt(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_app);
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT_app =
+ std::make_shared<const DiscreteFunctionVariant>(aT_d_app);
+
+ auto [ur_app, Fjr_app] = compute_fluxes(solver_type,rho_app,aL_app,aT_app,u_app,sigma_app,bc_descriptor_list,CR_ur,CR_Fjr,map2);
+ std::tie(new_m,new_rho,new_u,new_E,new_CG) = order2_apply_fluxes(dt2,new_rho,new_u,new_E,new_CG,rho_app,CG_app,ur_app,Fjr_app);
+
+ const DiscreteScalarFunction& new_rho_d = new_rho->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& new_u_d = new_u->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& new_E_d = new_E->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& new_eps = new_E_d - 0.5 * dot(new_u_d,new_u_d);
+ const DiscreteTensorFunction& new_CG_d = CG_app->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& new_p_d = fluid*(gamma-1) * new_rho_d * new_eps;
+ const DiscreteTensorFunction& new_sigma_d = (mu / sqrt(det(new_CG_d)) * (new_CG_d - I) + lambda / sqrt(det(new_CG_d)) * log(sqrt(det(new_CG_d))) * I) - new_p_d * I;
+ new_aL_d = sqrt((lambda + 2 * mu) / new_rho_d + gamma * new_p_d / new_rho_d);
+ new_aT_d = sqrt(mu / new_rho_d);
+
+ new_sigma = std::make_shared<const DiscreteFunctionVariant>(new_sigma_d);
+ new_aL = std::make_shared<const DiscreteFunctionVariant>(new_aL_d);
+ new_aT = std::make_shared<const DiscreteFunctionVariant>(new_aT_d);
+
+ sum_dt += dt2;
+
+ }
+
+ return {new_m, new_rho, new_u, new_E, new_CG};
+
+ }
+
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1428,76 +1527,122 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const double& mu,
const double& lambda) const
{
- std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr 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 m1_app = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+ std::shared_ptr<const DiscreteFunctionVariant> rho1_app = rho1;
+ std::shared_ptr<const DiscreteFunctionVariant> u1_app = u1;
+ std::shared_ptr<const DiscreteFunctionVariant> E1_app = E1;
+ std::shared_ptr<const DiscreteFunctionVariant> CG1_app = CG1;
+
+ std::shared_ptr m2_app = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr<const DiscreteFunctionVariant> rho2_app = rho2;
+ std::shared_ptr<const DiscreteFunctionVariant> u2_app = u2;
+ std::shared_ptr<const DiscreteFunctionVariant> E2_app = E2;
+ std::shared_ptr<const DiscreteFunctionVariant> CG2_app = CG2;
+ std::shared_ptr<const DiscreteFunctionVariant> aL2_app = aL2;
+ std::shared_ptr<const DiscreteFunctionVariant> aT2_app = aT2;
+ std::shared_ptr<const DiscreteFunctionVariant> sigma2_app = sigma2;
+
+
+ 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); //Fluxes t^n
+ std::tie(m1_app,rho1_app,u1_app,E1_app,CG1_app) = apply_fluxes(dt1/2, rho1, u1, E1, CG1, ur1, Fjr1); //Fluid t^n+(dt*0.5)
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);
+ if(dt2 > dt1/2){
+ dt2 = dt1/2;
+ }
- auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+ std::tie(m2_app,rho2_app,u2_app,E2_app,CG2_app) = apply_fluxes(dt2, rho2, u2, E2, CG2, ur2, Fjr2); //Solid t^n + dt2
- 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 double gamma = 1.4; //à modif
const TinyMatrix<Dimension> I = identity;
double sum_dt = dt2;
size_t fluid = 0;
- if ((mu == 0) and (lambda == 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;
+ while(sum_dt < dt1/2){ //Partie a revoir avec loi de comportement diff
+ const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2_d = u2_app->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E2_d = E2_app->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d, u2_d);
+ const DiscreteTensorFunction& CG2_d = CG2_app->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
+ const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
+ const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
+ const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
+
+ sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
+ aL2_app = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
+ aT2_app = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_app =
+ std::make_shared<const DiscreteFunctionVariant>(aL2_d_app + aT2_d_app);
+
+ dt2 = 0.4 * hyperelastic_dt(c_app);
+ if(sum_dt + dt2 > dt1/2){
+ dt2 = dt1/2 - 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);
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, rho2_app, aL2_app, aT2_app, u2_app, sigma2_app,
+ 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);
+ std::tie(m2_app, rho2_app, u2_app, E2_app, CG2_app) =
+ apply_fluxes(dt2, rho2_app, u2_app, E2_app, CG2_app, 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);
+ const DiscreteScalarFunction& rho1_d = rho1_app->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u1_d = u1_app->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E1_d = E1_app->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& eps1 = E1_d - 0.5 * dot(u1_d,u1_d);
+ const DiscreteTensorFunction& CG1_d = CG1_app->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p1_d = (1-fluid)*(gamma-1) * rho1_d * eps1;
+ const DiscreteTensorFunction& sigma1_d = fluid * (mu / sqrt(det(CG1_d)) * (CG1_d - I) + lambda / sqrt(det(CG1_d)) * log(sqrt(det(CG1_d))) * I) - p1_d * I;
+ const DiscreteScalarFunction& aL1_d = sqrt(fluid * (lambda + 2 * mu) / rho1_d + gamma * p1_d / rho1_d);
+ const DiscreteScalarFunction& aT1_d = sqrt(fluid * mu / rho1_d);
+
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_app = std::make_shared<const DiscreteFunctionVariant>(sigma1_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1_app = std::make_shared<const DiscreteFunctionVariant>(aL1_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1_app = std::make_shared<const DiscreteFunctionVariant>(aT1_d);
+
+ const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2_d = u2_app->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& E2_d = E2_app->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d, u2_d);
+ const DiscreteTensorFunction& CG2_d = CG2_app->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p2 = fluid * (gamma - 1) * rho2_d * eps2;
+ const DiscreteTensorFunction& sigma2_d = mu / sqrt(det(CG2_d)) * (CG2_d - I) + lambda / sqrt(det(CG2_d)) * log(sqrt(det(CG2_d))) * I - p2 * I;
+ const DiscreteScalarFunction& aL2_d_app = sqrt((lambda + 2 * mu) / rho2_d + gamma * p2 / rho2_d);
+ const DiscreteScalarFunction& aT2_d_app = sqrt(mu / rho2_d);
+
+ sigma2_app = std::make_shared<const DiscreteFunctionVariant>(sigma2_d);
+ aL2_app = std::make_shared<const DiscreteFunctionVariant>(aL2_d_app);
+ aT2_app = std::make_shared<const DiscreteFunctionVariant>(aT2_d_app);
+
+ auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes2(solver_type,rho1_app, aL1_app, aT1_app, u1_app, sigma1_app, bc_descriptor_list1,
+ rho2_app, aL2_app, aT2_app, u2_app, sigma2_app, bc_descriptor_list2, map1, map2); //Fluxes t^n + dt*0.5
+
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = order2_apply_fluxes(dt1, rho1, u1, E1, CG1, rho1_app, CG1_app, ur1_app, Fjr1_app);
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
+ CR_ur_app,CR_Fjr_app,map2,fluid,lambda,mu);
+
+ 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;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index 306ed1b60..58fe0b852 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -1,5 +1,5 @@
-#ifndef LOCAL_DT_HYPERELASTIC_SOLVER_HPP
-#define LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+#ifndef ORDER2_LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+#define ORDER2_LOCAL_DT_HYPERELASTIC_SOLVER_HPP
#include <mesh/MeshTraits.hpp>
--
GitLab