From 62d439149d1acc072119cd37c9eccd980d5f8b00 Mon Sep 17 00:00:00 2001
From: Alexandre Gangloff <alexandre.gangloff@cea.fr>
Date: Thu, 6 Mar 2025 11:48:02 +0100
Subject: [PATCH] Fix bug on the limiter and add one state law
---
src/language/modules/LocalFSIModule.cpp | 37 +--
src/scheme/Order2HyperelasticSolver.cpp | 411 +++++++-----------------
src/scheme/Order2HyperelasticSolver.hpp | 11 +-
3 files changed, 143 insertions(+), 316 deletions(-)
diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 76d2bc03f..9915d4e8d 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -213,7 +213,7 @@ LocalFSIModule::LocalFSIModule()
));
- this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
+ this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -226,9 +226,10 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& p,
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 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 double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -237,12 +238,12 @@ LocalFSIModule::LocalFSIModule()
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
.apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
- sigma, p, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
+ sigma, p, bc_descriptor_list, lambda, mu, gamma, p_inf);
}
));
- this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
+ this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -255,11 +256,8 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& p,
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 std::shared_ptr<const DiscreteFunctionVariant>& lambda,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu,
const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -267,8 +265,8 @@ LocalFSIModule::LocalFSIModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
- .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
- sigma, p, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
+ .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 3, dt, rho, u, E, CG, aL, aT,
+ sigma, p, bc_descriptor_list, lambda, mu, lambda, mu);
}
));
@@ -286,13 +284,10 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& p,
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 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 double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -301,7 +296,7 @@ LocalFSIModule::LocalFSIModule()
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
.apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt, rho, u, E, CG, aL, aT,
- sigma, p, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
+ sigma, p, bc_descriptor_list, lambda, mu, gamma, p_inf);
}
));
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 6d0dd07db..97ef40f91 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -457,8 +457,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const NodeValuePerCell<const Rdxd>& Ajr,
DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
- DiscreteFunctionDPk<Dimension, const double> DPk_p,
- CellValue<double> lambda, NodeValuePerCell<double> coef) const
+ DiscreteFunctionDPk<Dimension, const double> DPk_p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -480,7 +479,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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];
- const Rdxd sigma_r = sqrt(lambda[J]+(1-lambda[J])*coef(J,R))*DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;
+ const Rdxd sigma_r = DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;
br += Ajr(J, R) * DPk_u[J](xr[r]) - sigma_r * Cjr(J, R);
}
@@ -497,8 +496,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const NodeValue<const Rd>& ur,
DiscreteFunctionDPk<Dimension, const Rd> DPk_uh,
DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
- DiscreteFunctionDPk<Dimension, const double> DPk_p,
- CellValue<double> lambda, NodeValuePerCell<double> coef) const
+ DiscreteFunctionDPk<Dimension, const double> DPk_p) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -515,7 +513,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const auto& cell_nodes = cell_to_node_matrix[j];
for(size_t r = 0; r < cell_nodes.size(); ++r) {
- Rdxd sigma_r = sqrt(lambda[j]+(1-lambda[j])*coef(j,r))*DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
+ Rdxd sigma_r = DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
}
});
@@ -600,7 +598,7 @@ void
CellValue<Rdxd> sym_grad_u{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId j){
- sym_grad_u[j] = grad_u[j] + transpose(grad_u[j]);
+ sym_grad_u[j] = (grad_u[j] + transpose(grad_u[j]));
});
CellValue<SmallArray<double>> eigunvalues{mesh.connectivity()};
@@ -617,48 +615,51 @@ void
n[cell_id] = zero;
t[cell_id] = zero;
l[cell_id] = zero;
- if(frobeniusNorm(sym_grad_u[cell_id]) > 1E-5){
- if constexpr(Dimension == 2){
- Rd nj;
- Rd tj;
- for(size_t i = 0; i < Dimension; ++i){
- nj[i] = eigunvectors[cell_id][0][i];
- tj[i] = eigunvectors[cell_id][1][i];
- }
- n[cell_id] = nj;
- t[cell_id] = tj;
- if(l2Norm(nj) != 0){
- n[cell_id] = (1/l2Norm(nj))*n[cell_id];
- }
- if(l2Norm(tj) != 0){
- t[cell_id] = (1/l2Norm(tj))*t[cell_id];
- }
- }
- else{
- if constexpr(Dimension == 3){
+ //if(frobeniusNorm(sym_grad_u[cell_id]) > 1E-5){
+ if constexpr(Dimension == 2){
Rd nj;
Rd tj;
- Rd lj;
for(size_t i = 0; i < Dimension; ++i){
nj[i] = eigunvectors[cell_id][0][i];
tj[i] = eigunvectors[cell_id][1][i];
- lj[i] = eigunvectors[cell_id][2][i];
}
n[cell_id] = nj;
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];
}
- if(l2Norm(lj) != 0){
- l[cell_id] = (1/l2Norm(lj))*l[cell_id];
+ } else {
+ if constexpr(Dimension == 3){
+ Rd nj;
+ Rd tj;
+ Rd lj;
+ for(size_t i = 0; i < Dimension; ++i){
+ nj[i] = eigunvectors[cell_id][0][i];
+ tj[i] = eigunvectors[cell_id][1][i];
+ lj[i] = eigunvectors[cell_id][2][i];
+ }
+ n[cell_id] = nj;
+ t[cell_id] = tj;
+ l[cell_id] = lj;
+ if(l2Norm(nj) != 0){
+ n[cell_id] = (1/l2Norm(nj))*n[cell_id];
+ }
+ if(l2Norm(tj) != 0){
+ t[cell_id] = (1/l2Norm(tj))*t[cell_id];
+ }
+ if(l2Norm(lj) != 0){
+ l[cell_id] = (1/l2Norm(lj))*l[cell_id];
+ }
+ } else {
+ Rd nj;
+ nj[0] = 1.;
+ n[cell_id] = nj;
}
}
- }
- }
+ //}
});
parallel_for(
@@ -757,175 +758,7 @@ void
});
}
- void
- _vector_limiter(const MeshType& mesh,
- const DiscreteFunctionP0<const Rd>& f,
- DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
- const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const
- {
- MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
- StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
- auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
- symmetry_boundary_descriptor_list);
-
- const auto xj = mesh_data.xj();
-
- CellValue<Rd> n{mesh.connectivity()};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto coefficients_p = DPk_ph.coefficients(cell_id);
- Rd grad_p = zero;
- for(size_t i = 1; i < coefficients_p.size(); ++i){
- grad_p[i-1] = coefficients_p[i];
- }
- const double norm_grad_p = l2Norm(grad_p);
- if(norm_grad_p == 0){
- n[cell_id] = zero;
- }
- else {
- n[cell_id] = (1/norm_grad_p) * grad_p;
- }
-
-
- });
-
- const CellValue<Rd> t{mesh.connectivity()};
- const CellValue<Rd> l{mesh.connectivity()};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- const Rd nj = n[cell_id];
- Rd a = zero;
- if(l2Norm(nj) != 0){
- if((nj[0] / l2Norm(nj) != 1) and (nj[0] / l2Norm(nj) != -1)){
- a[0] = 1.;
- }
- else {
- a[1] = 1.;
- }
- }
-
- Rd tj = a - dot(a,nj) * nj;
- const double& norm_tj = l2Norm(tj);
- if(norm_tj == 0){
- tj = zero;
- }
- else {
- tj = (1/norm_tj) * tj;
- }
- t[cell_id] = tj;
-
- Rd lj = zero;
- if(Dimension == 3){
- lj[0] = nj[1]*tj[2] - nj[2]*tj[1];
- lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
- lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
- const double& norm_lj = l2Norm(lj);
- if(norm_lj != 0){
- lj = (1/norm_lj) * lj;
- }
-
- }
- l[cell_id] = lj;
- });
-
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
- const double fn = dot(f[cell_id],n[cell_id]);
- const double ft = dot(f[cell_id],t[cell_id]);
- const double fl = dot(f[cell_id],l[cell_id]);
-
- double fn_min = fn;
- double fn_max = fn;
- double ft_min = ft;
- double ft_max = ft;
- double fl_min = fl;
- double fl_max = fl;
-
- const auto cell_stencil = stencil[cell_id];
- for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_stencil[i_cell]]);
- fn_min = std::min(fn_min,fn_k);
- fn_max = std::max(fn_max,fn_k);
-
- const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_stencil[i_cell]]);
- ft_min = std::min(ft_min,ft_k);
- ft_max = std::max(ft_max,ft_k);
-
- const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_stencil[i_cell]]);
- fl_min = std::min(fl_min,fl_k);
- fl_max = std::max(fl_max,fl_k);
- }
-
- double fn_bar_min = fn;
- double fn_bar_max = fn;
- double ft_bar_min = ft;
- double ft_bar_max = ft;
- double fl_bar_min = fl;
- double fl_bar_max = fl;
-
- for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- const CellId cell_k_id = cell_stencil[i_cell];
- const Rd f_xk = DPK_fh[cell_id](xj[cell_k_id]);
-
- const double fn_xk = dot(f_xk,n[cell_k_id]);
- fn_bar_min = std::min(fn_bar_min,fn_xk);
- fn_bar_max = std::max(fn_bar_max,fn_xk);
-
- const double ft_xk = dot(f_xk,t[cell_k_id]);
- ft_bar_min = std::min(ft_bar_min,ft_xk);
- ft_bar_max = std::max(ft_bar_max,ft_xk);
-
- const double fl_xk = dot(f_xk,l[cell_k_id]);
- fl_bar_min = std::min(fl_bar_min,fl_xk);
- fl_bar_max = std::max(fl_bar_max,fl_xk);
- }
-
- const double eps = 1E-14;
- double coef1_n = 1;
- if(std::abs(fn_bar_max - fn) > eps){
- coef1_n = (fn_max - fn) / (fn_bar_max - fn);
- }
- double coef2_n = 1;
- if(std::abs(fn_bar_min - fn) > eps){
- coef2_n = (fn_min - fn) / (fn_bar_min - fn);
- }
-
- double coef1_t = 1;
- if(std::abs(ft_bar_max - ft) > eps){
- coef1_t = (ft_max - ft) / (ft_bar_max - ft);
- }
- double coef2_t = 1;
- if(std::abs(ft_bar_min - ft) > eps){
- coef2_t = (ft_min - ft) / (ft_bar_min - ft);
- }
-
- double coef1_l = 1;
- if(std::abs(fl_bar_max - fl) > eps){
- coef1_l = (fl_max - fl) / (fl_bar_max - fl);
- }
- double coef2_l = 1;
- if(std::abs(fl_bar_min - fl) > eps){
- coef2_l = (fl_min - fl) / (fl_bar_min - fl);
- }
-
- const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
- const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
- const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
-
- auto coefficients = DPK_fh.coefficients(cell_id);
- coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id]
- + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id]
- + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
- for(size_t i = 1; i < coefficients.size(); ++i){
- coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id]
- + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id]
- + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
- }
- });
- }
-
- CellValue<double>
+ void
_matrix_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rdxd>& S,
DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) const
@@ -939,8 +772,8 @@ void
const auto xj = mesh_data.xj();
//A retirer + tard
- const auto xr = mesh.xr();
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ //const auto xr = mesh.xr();
+ //const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
//
CellValue<double> lambda{mesh.connectivity()};
@@ -961,22 +794,22 @@ void
double inv2_bar_min = inv2j;
double inv2_bar_max = inv2j;
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- for(size_t r = 0; r < cell_nodes.size(); ++r){
- const double inv2_xr = 0.5*trace(transpose(DPk_Sh[cell_id](xr[cell_nodes[r]]))*DPk_Sh[cell_id](xr[cell_nodes[r]]));
-
- inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
- inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
- }
-
- //for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- // const CellId cell_k_id = cell_stencil[i_cell];
- // const double inv2_xk = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
+ //const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ //for(size_t r = 0; r < cell_nodes.size(); ++r){
+ // const double inv2_xr = 0.5*trace(transpose(DPk_Sh[cell_id](xr[cell_nodes[r]]))*DPk_Sh[cell_id](xr[cell_nodes[r]]));
//
- // inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
- // inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
+ // inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
+ // inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
//}
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ const CellId cell_k_id = cell_stencil[i_cell];
+ const double inv2_xk = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
+
+ inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
+ inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
+ }
+
const double eps = 1E-14;
double coef1 = 1.;
double coef2 = 1.;
@@ -987,18 +820,17 @@ void
coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
}
- //const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
- const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
+ const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
+ //const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
lambda[cell_id] = lambda_inv2;
- //auto coefficients = DPk_Sh.coefficients(cell_id);
+ auto coefficients = DPk_Sh.coefficients(cell_id);
- //coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
- //for (size_t i = 1; i < coefficients.size(); ++i) {
- // coefficients[i] *= lambda_inv2;
- //}
+ coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda_inv2;
+ }
});
- return lambda;
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1009,20 +841,45 @@ void
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 DiscreteScalarFunction lambda,
+ const DiscreteScalarFunction mu,
+ const DiscreteScalarFunction gamma,
+ const DiscreteScalarFunction 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 DiscreteScalarFunction& p_d = (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);
+ (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
+ const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
+
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(p_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>,
+ const std::shared_ptr<const DiscreteFunctionVariant>>
+ _apply_NeoHook(const DiscreteScalarFunction rho,
+ const DiscreteVectorFunction u,
+ const DiscreteScalarFunction E,
+ const DiscreteTensorFunction CG,
+ const DiscreteScalarFunction lambda,
+ const DiscreteScalarFunction mu) const
+ {
+ const TinyMatrix<Dimension> I = identity;
+
+ const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
+ const DiscreteTensorFunction& sigma_d =
+ (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I);
+ const DiscreteScalarFunction& p_d = -1./3. * trace(sigma_d);
+ const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(p_d),
@@ -1038,23 +895,20 @@ void
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 DiscreteScalarFunction mu,
+ const DiscreteScalarFunction gamma,
+ const DiscreteScalarFunction 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_f - 1) * rho * eps - gamma_f * p_inf_f) + chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
+ (gamma - 1) * rho * eps - gamma * p_inf;
const DiscreteTensorFunction& sigma_d =
- chi_solid * 2 * mu / sqrt(det(CG)) * (CG - 1. / 3. * trace(CG) * I) - p_d * I;
+ 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);
+ (2 * mu) / rho + (gamma * (p_d + p_inf)) / rho);
+ const DiscreteScalarFunction& aT_d = sqrt(mu / rho);
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(p_d),
@@ -1071,28 +925,26 @@ void
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 std::shared_ptr<const DiscreteFunctionVariant>& lambda_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& mu_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_inf_v) 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;
+ const DiscreteScalarFunction& lambda_d = lambda_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& mu_d = mu_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& gamma_d = gamma_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& p_inf_d = p_inf_v->get<DiscreteScalarFunction>();
if (law == 1) {
- return _apply_NeoHook(rho_d, u_d, E_d, CG_d, chi_solid_d, lambda, mu, gamma_f, p_inf_f);
+ return _apply_NeoHook(rho_d, u_d, E_d, CG_d, lambda_d, mu_d, gamma_d, p_inf_d);
+ } else if (law == 2){
+ return _apply_NeoHook2(rho_d, u_d, E_d, CG_d, mu_d, gamma_d, p_inf_d);
} else {
- return _apply_NeoHook2(rho_d, u_d, E_d, CG_d, chi_solid_d, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
+ return _apply_NeoHook(rho_d, u_d, E_d, CG_d,lambda_d, mu_d);
}
}
@@ -1193,30 +1045,13 @@ void
const CellValue<const Rdxd> grad_u = this->_computeGradU(solver_type,rho_v,aL_v,aT_v,u_v,sigma_v,bc_descriptor_list);
this->_vector_limiter(mesh,u,u_lim,grad_u);
- CellValue<double> lambda = this->_matrix_limiter(mesh, S, S_lim);
+ this->_matrix_limiter(mesh, S, S_lim);
this->_scalar_limiter(mesh, p, p_lim);
- const NodeValue<const Rd>& xr = mesh.xr();
- const DiscreteScalarFunction& J2 = 0.5*trace(transpose(S)*S);
- NodeValuePerCell<double> coef{mesh.connectivity()};
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- 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){
- const double J2_p = 0.5*trace(transpose(DPk_Sh[j](xr[cell_nodes[r]]))*DPk_Sh[j](xr[cell_nodes[r]]));
- if(J2_p != 0.){
- coef(j,r) = J2[j]/J2_p;
- }else{
- coef(j,r) = 0;
- }
- }
- });
-
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_lim, S_lim, p_lim, lambda, coef);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -1225,7 +1060,7 @@ void
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
- NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -1622,23 +1457,23 @@ void
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
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
+ 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
{
auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, p, bc_descriptor_list);
auto [ur_o1, Fjr_o1] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
- std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
+ std::shared_ptr<const DiscreteFunctionVariant> lambda_app = shallowCopy(m_app, lambda);
+ std::shared_ptr<const DiscreteFunctionVariant> mu_app = shallowCopy(m_app, mu);
+ std::shared_ptr<const DiscreteFunctionVariant> gamma_app = shallowCopy(m_app, gamma);
+ std::shared_ptr<const DiscreteFunctionVariant> p_inf_app = shallowCopy(m_app, p_inf);
auto [sigma_app, p_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);
+ _apply_state_law(law,rho_app,u_app,E_app,CG_app,lambda_app,mu_app,gamma_app,p_inf_app);
auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app, p_app, bc_descriptor_list);
auto [ur_o1_app, Fjr_o1_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
diff --git a/src/scheme/Order2HyperelasticSolver.hpp b/src/scheme/Order2HyperelasticSolver.hpp
index 99f961604..f6b53ca1b 100644
--- a/src/scheme/Order2HyperelasticSolver.hpp
+++ b/src/scheme/Order2HyperelasticSolver.hpp
@@ -42,13 +42,10 @@ class Order2HyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
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 = 0;
+ 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 = 0;
IOrder2HyperelasticSolver() = default;
IOrder2HyperelasticSolver(IOrder2HyperelasticSolver&&) = default;
--
GitLab