diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 3ed298e6e6ad59e51a123171f9ba5ddd5f0d20c9..76d2bc03fb856acf075d69e91fc6c15d459d50ec 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -402,6 +402,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -411,6 +412,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2,
const double& lambda1, const double& mu1,
@@ -431,13 +433,13 @@ LocalFSIModule::LocalFSIModule()
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
.apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
- u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, p1, p2, bc_descriptor_list1,
bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
}
));
- this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+ this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook2_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -447,6 +449,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
@@ -456,6 +459,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2,
const double& lambda1, const double& mu1,
@@ -476,8 +480,8 @@ LocalFSIModule::LocalFSIModule()
getCommonMesh(
{rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
- u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt1, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, p1, p2, bc_descriptor_list1,
bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, p_inf1, gamma2, p_inf2);
}
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index 96f1db1a0872e31e481c651f0a90e7533d1dfdf4..6d0dd07db4451722c3b3518910d9c40c91360122 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -326,72 +326,6 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
this->_applyVelocityBC(bc_list, Ar, br);
}
- CellValue<Rdxd>
- _computeEigenVectorsMatrix(const MeshType& mesh,
- const CellValue<Rdxd>& grad_v,
- const double& eps) const
- {
- CellValue<Rdxd> A = copy(grad_v);
-
- CellValue<Rdxd> V{mesh.connectivity()};
- const TinyMatrix<Dimension> I = identity;
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- V[j] = I;
-
- while(true){
- double max_val = 0;
- size_t p = 0;
- size_t q = 1;
- for(size_t i = 0; i < Dimension; ++i){
- for(size_t k = i + 1; k < Dimension; ++k){
- if(abs(A[j](i,k)) > max_val){
- max_val = abs(A[j](i,k));
- p = i;
- q = k;
- }
- }
- }
- if(max_val < eps){
- break;
- }
-
- const double& theta = 0.5*(A[j](q,q) - A[j](p,p)) / A[j](p, q);
- double t;
- if(theta >= 0){
- t = 1/(theta + sqrt(theta*theta + 1));
- }else{
- t = -1/(-theta + sqrt(theta*theta + 1));
- }
- const double& c = 1/sqrt(t*t + 1);
- const double& s = t * c;
-
- for(size_t i = 0; i < Dimension; ++i){
- if((i != p) and (i!=q)){
- const double& temp_p = c * A[j](i,p) - s * A[j](i,q);
- const double& temp_q = q * A[j](i,p) + c * A[j](i,q);
- A[j](i,p) = temp_p;
- A[j](i,q) = temp_q;
- A[j](p,i) = temp_p;
- A[j](q,i) = temp_q;
- }
- }
- A[j](p,p) = A[j](p,p) - t * A[j](p,q);
- A[j](q,q) = A[j](q,q) + t * A[j](p,q);
- A[j](p,q) = 0.;
- A[j](q,p) = 0.;
-
- for(size_t i = 0; i < Dimension; ++i){
- V[j](i,p) = c * V[j](i,p) - s * V[j](i,q);
- V[j](i,q) = s * V[j](i,p) + c * V[j](i,q);
- }
- }
-
- });
-
- return V;
- }
-
CellValue<Rdxd>
_computeGradU(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
@@ -674,8 +608,6 @@ void
for(CellId j = 0; j < mesh.numberOfCells(); ++j){
EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
}
-
- //CellValue<Rdxd> V = _computeEigenVectorsMatrix(mesh,sym_grad_u,1E-4);
CellValue<Rd> n{mesh.connectivity()};
CellValue<Rd> t{mesh.connectivity()};
@@ -689,10 +621,7 @@ void
if constexpr(Dimension == 2){
Rd nj;
Rd tj;
- //Rdxd Vj = V[cell_id];
for(size_t i = 0; i < Dimension; ++i){
- // nj[i] = Vj(i,0);
- // tj[i] = Vj(i,1);
nj[i] = eigunvectors[cell_id][0][i];
tj[i] = eigunvectors[cell_id][1][i];
}
@@ -710,11 +639,7 @@ void
Rd nj;
Rd tj;
Rd lj;
- // Rdxd Vj = V[cell_id];
for(size_t i = 0; i < Dimension; ++i){
- // nj[i] = Vj(i,0);
- // tj[i] = Vj(i,1);
- // lj[i] = Vj(i,2);
nj[i] = eigunvectors[cell_id][0][i];
tj[i] = eigunvectors[cell_id][1][i];
lj[i] = eigunvectors[cell_id][2][i];
@@ -1268,7 +1193,6 @@ 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);
- //this->_vector_limiter(mesh, u, u_lim, DPk_ph);
CellValue<double> lambda = this->_matrix_limiter(mesh, S, S_lim);
this->_scalar_limiter(mesh, p, p_lim);
@@ -1706,12 +1630,6 @@ void
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;
- //std::shared_ptr<const DiscreteFunctionVariant> u_app = u;
- //std::shared_ptr<const DiscreteFunctionVariant> E_app = E;
- //std::shared_ptr<const DiscreteFunctionVariant> CG_app = CG;
-
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);
@@ -1725,7 +1643,6 @@ void
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);
return order2_apply_fluxes(dt, rho, u, E, CG, rho_app, CG_app, ur_app, ur_o1_app, Fjr_app, Fjr_o1_app);
- //return apply_fluxes(dt, rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
}
Order2HyperelasticSolver() = default;
diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
index a80ae83fed4234da7fd65f7ef9248406d496fbf5..a98d525a12cbf7f0458f72155bd8b8ae807ba510 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp
@@ -1,3 +1,4 @@
+#include <algebra/EigenvalueSolver.hpp>
#include <scheme/Order2LocalDtHyperelasticSolver.hpp>
#include <language/utils/InterpolateItemValue.hpp>
@@ -221,6 +222,78 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return b;
}
+ NodeValue<Rd>
+ _computeBr(const Mesh<Dimension>& mesh,
+ 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
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+ const auto& xr = mesh.xr();
+ const TinyMatrix<Dimension> I = identity;
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rd> b{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ Rd br = zero;
+ 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;
+
+ br += Ajr(J, R) * DPk_u[J](xr[r]) - sigma_r * Cjr(J, R);
+ }
+
+ b[r] = br;
+ });
+
+ return b;
+ }
+
+ NodeValue<Rd>
+ _computeBr(const Mesh<Dimension>& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const DiscreteVectorFunction& u,
+ const DiscreteTensorFunction& sigma) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rd> b{mesh.connectivity()};
+
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ Rd br = zero;
+ 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];
+ br += Ajr(J, R) * u[J] - sigma[J] * Cjr(J, R);
+ }
+
+ b[r] = br;
+ });
+
+ return b;
+ }
+
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
@@ -354,6 +427,63 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
this->_applyVelocityBC(bc_list, Ar, br);
}
+ CellValue<Rdxd>
+ _computeGradU(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
+
+ const MeshType& mesh = *mesh_v->get<MeshType>();
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+
+ 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, sigma);
+
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+
+ synchronize(Ar);
+ synchronize(br);
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+ const NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+
+ CellValue<Rdxd> grad_u{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ auto cell_nodes = cell_to_node_matrix[j];
+
+ Rdxd gradv = zero;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ NodeId r = cell_nodes[R];
+ gradv += tensorProduct(ur[r], Cjr(j, R));
+ }
+ grad_u[j] = gradv;
+ });
+
+ return grad_u;
+ }
+
NodeValue<Rd>
_computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
@@ -395,6 +525,64 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return F;
}
+ NodeValuePerCell<Rd>
+ _computeFjr(const MeshType& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ 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
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValue<const Rd>& xr = mesh.xr();
+ const TinyMatrix<Dimension> I = identity;
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ NodeValuePerCell<Rd> F{mesh.connectivity()};
+
+ 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) {
+ 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;
+ F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
+ }
+ });
+
+ return F;
+ }
+
+ NodeValuePerCell<Rd>
+ _computeFjr(const MeshType& mesh,
+ const NodeValuePerCell<const Rdxd>& Ajr,
+ const NodeValue<const Rd>& ur,
+ const DiscreteVectorFunction& u,
+ const DiscreteTensorFunction& sigma) const
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+ NodeValuePerCell<Rd> F{mesh.connectivity()};
+ 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) {
+ F(j, r) = -Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + sigma[j] * Cjr(j, r);
+ }
+ });
+
+ return F;
+ }
+
std::tuple<std::vector<NodeId>, std::vector<NodeId>>
_computeMapping(std::shared_ptr<const MeshVariant>& i_mesh1,
std::shared_ptr<const MeshVariant>& i_mesh2,
@@ -450,10 +638,10 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return std::make_tuple(map1, map2);
}
- void
- _vector_limiter(const MeshType& mesh,
- const DiscreteFunctionP0<const Rd>& f,
- DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const
+ void
+ _scalar_limiter(const MeshType& mesh,
+ const DiscreteFunctionP0<const double>& f,
+ DiscreteFunctionDPk<Dimension, double>& DPk_fh) const
{
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
@@ -462,58 +650,44 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const auto xj = mesh_data.xj();
+ CellValue<double> alph_J2{mesh.connectivity()};
+
parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
- const Rd fj = f[cell_id];
+ const double fj = f[cell_id];
- Rd f_min = fj;
- Rd f_max = fj;
+ double f_min = fj;
+ double f_max = fj;
const auto cell_stencil = stencil[cell_id];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- for(size_t dim = 0; dim < Dimension; ++dim){
- f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]);
- f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]);
- }
+ f_min = std::min(f_min, f[cell_stencil[i_cell]]);
+ f_max = std::max(f_max, f[cell_stencil[i_cell]]);
}
- Rd f_bar_min = fj;
- Rd f_bar_max = fj;
+ double f_bar_min = fj;
+ double f_bar_max = fj;
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 f_xk = DPk_fh[cell_id](xj[cell_k_id]);
- for(size_t dim = 0; dim < Dimension; ++dim){
- f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
- f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
- }
+ f_bar_min = std::min(f_bar_min, f_xk);
+ f_bar_max = std::max(f_bar_max, f_xk);
}
const double eps = 1E-14;
- Rd coef1;
- for(size_t dim = 0; dim < Dimension; ++dim){
- coef1[dim] = 1;
- if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
- coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
- }
- }
-
- Rd coef2;
- for(size_t dim = 0; dim < Dimension; ++dim){
- coef2[dim] = 1;
- if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
- coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
- }
+ double coef1 = 1;
+ if (std::abs(f_bar_max - fj) > eps) {
+ coef1 = (f_max - fj) / ((f_bar_max - fj));
}
- double min_coef1 = coef1[0];
- double min_coef2 = coef2[0];
- for(size_t dim = 1; dim < Dimension; ++dim){
- min_coef1 = std::min(min_coef1,coef1[dim]);
- min_coef2 = std::min(min_coef2,coef2[dim]);
+ double coef2 = 1.;
+ if (std::abs(f_bar_min - fj) > eps) {
+ coef2 = (f_min - fj) / ((f_bar_min - fj));
}
- const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+ const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+ alph_J2[cell_id] = lambda;
auto coefficients = DPk_fh.coefficients(cell_id);
@@ -524,155 +698,244 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
});
}
- CellValue<double>
- _scalar_limiter_coefficient(const MeshType& mesh,
- const DiscreteFunctionP0<const double>& f,
- const DiscreteFunctionDPk<Dimension, const double>& DPk_fh) const
+ void
+ _vector_limiter(const MeshType& mesh,
+ const DiscreteFunctionP0<const Rd>& f,
+ DiscreteFunctionDPk<Dimension, Rd> DPk_fh,
+ const CellValue<const Rdxd>& grad_u) 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);
-
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+ symmetry_boundary_descriptor_list);
+
const auto xj = mesh_data.xj();
- CellValue<double> lambda{mesh.connectivity()};
-
- parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
- const double fj = f[cell_id];
+ 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]);
+ });
- double f_min = fj;
- double f_max = fj;
+ CellValue<SmallArray<double>> eigunvalues{mesh.connectivity()};
+ CellValue<std::vector<SmallVector<double>>> eigunvectors{mesh.connectivity()};
+ for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+ EigenvalueSolver{}.computeForSymmetricMatrix(sym_grad_u[j],eigunvalues[j],eigunvectors[j]);
+ }
- const auto cell_stencil = stencil[cell_id];
- for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- f_min = std::min(f_min, f[cell_stencil[i_cell]]);
- f_max = std::max(f_max, f[cell_stencil[i_cell]]);
- }
+ CellValue<Rd> n{mesh.connectivity()};
+ CellValue<Rd> t{mesh.connectivity()};
+ CellValue<Rd> l{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ 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){
+ 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];
+ }
+ }
+ }
+ }
+ });
- double f_bar_min = fj;
- double f_bar_max = fj;
+ 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_id]);
+ 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_id]);
+ 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_id]);
+ fl_min = std::min(fl_min,fl_k);
+ fl_max = std::max(fl_max,fl_k);
+ }
- for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
- const CellId cell_k_id = cell_stencil[i_cell];
- const double f_xk = DPk_fh[cell_id](xj[cell_k_id]);
+ 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;
- f_bar_min = std::min(f_bar_min, f_xk);
- f_bar_max = std::max(f_bar_max, f_xk);
- }
+ 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 eps = 1E-14;
- double coef1 = 1;
- if (std::abs(f_bar_max - fj) > eps) {
- coef1 = (f_max - fj) / ((f_bar_max - fj));
- }
+ const double fn_xk = dot(f_xk,n[cell_id]);
+ fn_bar_min = std::min(fn_bar_min,fn_xk);
+ fn_bar_max = std::max(fn_bar_max,fn_xk);
- double coef2 = 1.;
- if (std::abs(f_bar_min - fj) > eps) {
- coef2 = (f_min - fj) / ((f_bar_min - fj));
- }
+ const double ft_xk = dot(f_xk,t[cell_id]);
+ ft_bar_min = std::min(ft_bar_min,ft_xk);
+ ft_bar_max = std::max(ft_bar_max,ft_xk);
- lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
- });
+ const double fl_xk = dot(f_xk,l[cell_id]);
+ fl_bar_min = std::min(fl_bar_min,fl_xk);
+ fl_bar_max = std::max(fl_bar_max,fl_xk);
+ }
- return lambda;
- }
+ 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);
+ }
- NodeValuePerCell<Rdxd>
- _compute_tensor_reconstruction(const std::shared_ptr<const MeshVariant>& v_mesh,
- DiscreteTensorFunction sigma,
- PolynomialReconstructionDescriptor reconstruction_descriptor) const
- {
- const MeshType& mesh = *v_mesh->get<MeshType>();
-
- CellValue<double> limiter_j{mesh.connectivity()};
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
- limiter_j[j] = 1;
- }
- );
-
- const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
- std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
- std::vector<size_t> row;
- std::vector<size_t> col;
- for(size_t i = 0; i < Dimension; ++i){
- for(size_t k = i; k < Dimension; ++k){
-
- CellValue<double> coef{mesh.connectivity()};
- for(CellId j = 0; j <mesh.numberOfCells(); ++j){
- coef[j] = sigma_j[j](i,k);
- }
-
- const auto& coef_scalar_function = DiscreteScalarFunction(v_mesh, coef);
- auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
- auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
-
- const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef);
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
- });
-
- sigma_coef.push_back(copy(DPk_coef));
- row.push_back(i);
- col.push_back(k);
- }
+ 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);
}
- NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ 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];
+ }
+ });
+ }
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
- const auto& cell_nodes = cell_to_node_matrix[j];
- for(size_t i = 0; i < Dimension; ++i){
- for(size_t k = 0; k < Dimension; ++k){
- for(size_t R = 0; R < cell_nodes.size(); ++R){
- sigma_lim(j,R)(i,k) = 0;
- }
- }
- }
- }
- );
+ CellValue<double>
+ _matrix_limiter(const MeshType& mesh,
+ const DiscreteFunctionP0<const Rdxd>& S,
+ DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) 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 NodeValue<const Rd>& xr = copy(mesh.xr());
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
+ const auto xj = mesh_data.xj();
- for(size_t R = 0; R < cell_nodes.size(); ++R){
- const NodeId r = cell_nodes[R];
+ //A retirer + tard
+ const auto xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ //
- for(size_t l = 0; l < sigma_coef.size(); ++l){
- const size_t& i = row[l];
- const size_t& k = col[l];
+ CellValue<double> lambda{mesh.connectivity()};
+ const DiscreteScalarFunction& inv2 = 0.5*trace(transpose(S)*S);
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+ const double inv2j = inv2[cell_id];
+
+ double inv2_min = inv2j;
+ double inv2_max = inv2j;
- auto coefficients = sigma_coef[l].coefficients(j);
+ const auto cell_stencil = stencil[cell_id];
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ inv2_min = std::min(inv2_min, inv2[cell_stencil[i_cell]]);
+ inv2_max = std::max(inv2_max, inv2[cell_stencil[i_cell]]);
+ }
- coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
- for (size_t indice = 1; indice < coefficients.size(); ++indice) {
- coefficients[indice] *= limiter_j[j];
- }
+ double inv2_bar_min = inv2j;
+ double inv2_bar_max = inv2j;
- sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
- if(i != k){
- sigma_lim(j,R)(i,k) *= 2;
- }
- }
+ 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]]));
- sigma_lim(j,R) += transpose(sigma_lim(j,R));
- sigma_lim(j,R) *= 0.5;
- }
+ inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
+ inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
+ }
- }
- );
+ const double eps = 1E-14;
+ double coef1 = 1.;
+ double coef2 = 1.;
+ if(std::abs(inv2_bar_max - inv2j) > eps){
+ coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j);
+ }
+ if(std::abs(inv2_bar_min - inv2j) > eps){
+ coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
+ }
- return sigma_lim;
- }
+ const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
+ lambda[cell_id] = lambda_inv2;
+ });
+ return lambda;
+ }
- std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ 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,
@@ -693,10 +956,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(aL_d),
- std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ std::make_shared<DiscreteFunctionVariant>(aT_d),
+ std::make_shared<DiscreteFunctionVariant>(p_d)};
}
- std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ 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_NeoHook2(const DiscreteScalarFunction rho,
const DiscreteVectorFunction u,
const DiscreteScalarFunction E,
@@ -716,10 +983,14 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(aL_d),
- std::make_shared<DiscreteFunctionVariant>(aT_d)};
+ std::make_shared<DiscreteFunctionVariant>(aT_d),
+ std::make_shared<DiscreteFunctionVariant>(p_d)};
}
- std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>, const std::shared_ptr<const DiscreteFunctionVariant>>
+ 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_state_law(const size_t law,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
@@ -739,70 +1010,187 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
if(law == 1){
return _apply_NeoHook(rho_d,u_d,E_d,CG_d,chi_solid,lambda,mu,gamma,p_inf);
} else {
- return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid,lambda,gamma,p_inf);
+ return _apply_NeoHook2(rho_d,u_d,E_d,CG_d,chi_solid,mu,gamma,p_inf);
}
}
+ NodeValuePerCell<double>
+ _compute_inv2_coefficient(const MeshType& mesh,
+ const DiscreteTensorFunction& S,
+ const DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) const
+ {
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const NodeValue<const Rd>& xr = mesh.xr();
+
+ const DiscreteScalarFunction& J2 = 0.5 * trace(transpose(S)*S);
+ NodeValuePerCell<double> coef{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const 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;
+ }
+ }
+ });
+
+ return coef;
+ }
+
public:
- std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
compute_fluxes(const SolverType& solver_type,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
{
- std::shared_ptr i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
- if (not i_mesh) {
+ std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
+ std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
+ if (not i_mesh1) {
throw NormalError("discrete functions are not defined on the same mesh");
}
- if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+ 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");
+ }
- const MeshType& mesh = *i_mesh->get<MeshType>();
- const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
- const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+ if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
- std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+ const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL1 = aL1_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT1 = aT1_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p1 = p1_v->get<DiscreteScalarFunction>();
- for(auto&& bc_descriptor : bc_descriptor_list){
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
+ const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL2 = aL2_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p2 = p2_v->get<DiscreteScalarFunction>();
+
+ const TinyMatrix<Dimension> I = identity;
+ const DiscreteTensorFunction& S1 = sigma1 + p1*I;
+ const DiscreteTensorFunction& S2 = sigma2 + p2*I;
+
+ const std::shared_ptr<const DiscreteFunctionVariant>& S1_v = std::make_shared<DiscreteFunctionVariant>(S1);
+ const std::shared_ptr<const DiscreteFunctionVariant>& S2_v = std::make_shared<DiscreteFunctionVariant>(S2);
+
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+
+ for(auto&& bc_descriptor : bc_descriptor_list1){
if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
- symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+ for(auto&& bc_descriptor : bc_descriptor_list2){
+ if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
+ symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
}
}
- PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list);
+ PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
+ symmetry_boundary_descriptor_list1);
+ PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
+ symmetry_boundary_descriptor_list2);
- NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
- auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v, S1_v, p1_v);
+ auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto DPk_Sh1 = reconstruction1[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+ auto DPk_ph1 = reconstruction1[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
- DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
- this->_vector_limiter(mesh,u,u_lim);
+ DiscreteFunctionDPk<Dimension, Rd> u1_lim = copy(DPk_uh1);
+ DiscreteFunctionDPk<Dimension, Rdxd> S1_lim = copy(DPk_Sh1);
+ DiscreteFunctionDPk<Dimension, double> p1_lim = copy(DPk_ph1);
- NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
+ const CellValue<const Rdxd> grad_u1 = this->_computeGradU(solver_type,rho1_v,aL1_v,aT1_v,u1_v,sigma1_v,bc_descriptor_list1);
+ this->_vector_limiter(mesh1,u1,u1_lim,grad_u1);
+ CellValue<double> lambda1 = this->_matrix_limiter(mesh1, S1, S1_lim);
+ this->_scalar_limiter(mesh1, p1, p1_lim);
- NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
- NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
+ auto copy_DPk_Sh1 = copy(DPk_Sh1);
+ const NodeValuePerCell<double> coef1 = _compute_inv2_coefficient(mesh1, S1, copy_DPk_Sh1);
- const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
- this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+ auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v, S2_v, p2_v);
+ auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto DPk_Sh2 = reconstruction2[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+ auto DPk_ph2 = reconstruction2[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
- synchronize(Ar);
- synchronize(br);
+ DiscreteFunctionDPk<Dimension, Rd> u2_lim = copy(DPk_uh2);
+ DiscreteFunctionDPk<Dimension, Rdxd> S2_lim = copy(DPk_Sh2);
+ DiscreteFunctionDPk<Dimension, double> p2_lim = copy(DPk_ph2);
- NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
- NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+ const CellValue<const Rdxd> grad_u2 = this->_computeGradU(solver_type,rho2_v,aL2_v,aT2_v,u2_v,sigma2_v,bc_descriptor_list2);
+ this->_vector_limiter(mesh2,u2,u2_lim,grad_u2);
+ CellValue<double> lambda2 = this->_matrix_limiter(mesh2, S2, S2_lim);
+ this->_scalar_limiter(mesh2, p2, p2_lim);
- return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+ auto copy_DPk_Sh2 = copy(DPk_Sh2);
+ const NodeValuePerCell<double> coef2 = _compute_inv2_coefficient(mesh2, S2, copy_DPk_Sh2);
+
+ NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
+ NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
+
+ NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+ NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1_lim, S1_lim, p1_lim, lambda1, coef1);
+ NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+ NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2_lim, S2_lim, p2_lim, lambda2, coef2);
+
+ const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+ this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+ this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+
+ synchronize(Ar1);
+ synchronize(br1);
+ synchronize(Ar2);
+ synchronize(br2);
+
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+
+ this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, S1_lim, p1_lim, lambda1, coef1);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, S2_lim, p2_lim, lambda2, coef2);
+
+ NodeValue<Rd> CR_ur = copy(ur2);
+ NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
}
std::tuple<const std::shared_ptr<const ItemValueVariant>,
@@ -872,32 +1260,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
}
}
- PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list1);
- PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list2);
-
- NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
- NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
-
- auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
- auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
- auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
- auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
-
- DiscreteFunctionDPk<Dimension, Rd> u1_lim = copy(DPk_uh1);
- this->_vector_limiter(mesh1,u1,u1_lim);
- DiscreteFunctionDPk<Dimension, Rd> u2_lim = copy(DPk_uh2);
- this->_vector_limiter(mesh2,u2,u2_lim);
-
NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
- NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim);
+ NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1);
NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
- NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
- this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
+ NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2);
const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
@@ -910,11 +1279,15 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
synchronize(br2);
NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
- NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
- NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
- NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
+
+ this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+
+ NodeValue<Rd> CR_ur = copy(ur2);
+ NodeValuePerCell<Rd> CR_Fjr = copy(Fjr2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
@@ -929,6 +1302,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
NodeValue<Rd> CR_ur,
NodeValuePerCell<Rd> CR_Fjr,
@@ -949,6 +1323,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
+ const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
+
+ const TinyMatrix<Dimension> I = identity;
+ const DiscreteTensorFunction& S = sigma + p*I;
+ const std::shared_ptr<const DiscreteFunctionVariant>& S_v = std::make_shared<DiscreteFunctionVariant>(S);
std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
@@ -961,148 +1340,89 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
symmetry_boundary_descriptor_list);
- NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor);
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v);
- auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, S_v, p_v);
+ auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto DPk_Sh = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+ auto DPk_ph = reconstruction[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
+ DiscreteFunctionDPk<Dimension, Rdxd> S_lim = copy(DPk_Sh);
+ DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
+
+ 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->_scalar_limiter(mesh, p, p_lim);
- DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
- this->_vector_limiter(mesh,u,u_lim);
+ auto copy_DPk_Sh = copy(DPk_Sh);
+ const NodeValuePerCell<double> coef = _compute_inv2_coefficient(mesh, S, copy_DPk_Sh);
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, sigma_lim);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim, lambda, coef);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
- synchronize(Ar);
- synchronize(br);
-
- NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
- NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
-
- this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
-
- return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr));
- }
-
- std::tuple<const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- NodeValue<Rd>,
- NodeValuePerCell<Rd>>
- compute_fluxes2(const SolverType& solver_type,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const
- {
- std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
- std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
- if (not i_mesh1) {
- throw NormalError("discrete functions are not defined on the same mesh 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 ici2 mesh2");
- }
-
- if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
- throw NormalError("acoustic solver expects P0 functions");
- }
+ synchronize(Ar);
+ synchronize(br);
- const MeshType& mesh1 = *i_mesh1->get<MeshType>();
- const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& aL1 = aL1_v->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& aT1 = aT1_v->get<DiscreteScalarFunction>();
- const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+ NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
- const MeshType& mesh2 = *i_mesh2->get<MeshType>();
- const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
- const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
- const DiscreteScalarFunction& aL2 = aL2_v->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>();
- const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+ this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
- std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1;
- std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2;
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+ }
- for(auto&& bc_descriptor : bc_descriptor_list1){
- if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
- symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared());
- }
- }
- for(auto&& bc_descriptor : bc_descriptor_list2){
- if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){
- symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared());
- }
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ const std::vector<NodeId> map2) const
+ {
+ 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");
}
- PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list1);
- PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list2);
-
- NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1);
- NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2);
+ if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("hyperelastic solver expects P0 functions");
+ }
- auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v);
- auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
- auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v);
- auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ const MeshType& mesh = *i_mesh->get<MeshType>();
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
+ const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
- DiscreteFunctionDPk<Dimension, Rd> u1_lim = copy(DPk_uh1);
- this->_vector_limiter(mesh1,u1,u1_lim);
- DiscreteFunctionDPk<Dimension, Rd> u2_lim = copy(DPk_uh2);
- this->_vector_limiter(mesh2,u2,u2_lim);
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
- NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
- NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma);
- NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
- NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim);
- NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
- NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim);
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
- const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
- const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
- this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
- this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+ synchronize(Ar);
+ synchronize(br);
- synchronize(Ar1);
- synchronize(br1);
- synchronize(Ar2);
- synchronize(br2);
+ NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
- NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
- NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
- NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim);
- NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim);
+ this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
- return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
- std::make_shared<const ItemValueVariant>(ur2),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
- ur2,
- Fjr2);
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1116,21 +1436,24 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteTensorFunction& CG,
- const NodeValue<const Rd>& ur,
- const NodeValuePerCell<const Rd>& Fjr) const
+ const NodeValue<const Rd>& ur_o2,
+ const NodeValue<const Rd>& ur_o1,
+ const NodeValuePerCell<const Rd>& Fjr_o2,
+ const NodeValuePerCell<const Rd>& Fjr_o1) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ NodeValue<Rd> ur = copy(ur_o2);
+ NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
+
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
- NodeValue<Rd> new_xr = copy(mesh.xr());
- parallel_for(
- mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
- std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+ 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);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
@@ -1140,29 +1463,105 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
+ CellValue<double> new_rho_stencil = copy(rho.cellValues());
+ CellValue<Rd> new_u_stencil = copy(u.cellValues());
+ CellValue<double> new_E_stencil = copy(E.cellValues());
+ CellValue<Rdxd> new_CG_stencil = copy(CG.cellValues());
+
+ CellValue<int> cell_flag{mesh.connectivity()};
+
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
+ auto cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
- const NodeId r = cell_nodes[R];
+ NodeId r = cell_nodes[R];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
- const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
- const double dt_over_Mj = dt / (rho[j] * Vj[j]);
- const double dt_over_Vj = dt / Vj[j];
+ Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+ double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_CG[j] += dt_over_Vj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
+
+ double new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
+ if(new_eps < 0.){
+ cell_flag[j] = 1;
+ }else{
+ cell_flag[j] = 0;
+ }
});
+ for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+ if(cell_flag[j] == 1){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ Rdxd gradv = zero;
+ std::vector<NodeId> NodeList;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ NodeId r = cell_nodes[R];
+ NodeList.push_back(r);
+ Fjr(j,R) = Fjr_o1(j,R);
+ ur[r] = ur_o1[r];
+ gradv += tensorProduct(ur[r], Cjr(j, R));
+ momentum_fluxes += Fjr(j, R);
+ energy_fluxes += dot(Fjr(j, R), ur[r]);
+ }
+ Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+ double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ double dt_over_Vj = dt / Vj[j];
+ new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
+ new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
+ new_CG[j] = new_CG_stencil[j] + dt_over_Vj * cauchy_green_fluxes;
+ new_CG[j] += transpose(new_CG[j]);
+ new_CG[j] *= 0.5;
+
+ auto cell_stencil = stencil[j];
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ CellId J = cell_stencil[i_cell];
+ auto i_cell_nodes = cell_to_node_matrix[J];
+ momentum_fluxes = zero;
+ energy_fluxes = 0;
+ gradv = zero;
+ for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
+ NodeId i_node = i_cell_nodes[R];
+ for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
+ if(NodeList[node_test] == i_node){
+ Fjr(J,R) = Fjr_o1(J,R);
+ }
+ }
+ gradv += tensorProduct(ur[i_node], Cjr(J, R));
+ momentum_fluxes += Fjr(J, R);
+ energy_fluxes += dot(Fjr(J, R), ur[i_node]);
+ }
+ cauchy_green_fluxes = gradv * CG[J] + CG[J] * transpose(gradv);
+ dt_over_Mj = dt / (rho[J] * Vj[J]);
+ dt_over_Vj = dt / Vj[J];
+ new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
+ new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
+ new_CG[J] = new_CG_stencil[J] + dt_over_Vj * cauchy_green_fluxes;
+ new_CG[J] += transpose(new_CG[J]);
+ new_CG[J] *= 0.5;
+ }
+ }
+ }
+
+ NodeValue<Rd> new_xr = copy(mesh.xr());
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
@@ -1186,11 +1585,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& CG,
const std::shared_ptr<const ItemValueVariant>& ur,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+ const std::shared_ptr<const ItemValueVariant>& ur_o1,
+ const std::shared_ptr<const SubItemValuePerItemVariant> Fjr,
+ const std::shared_ptr<const SubItemValuePerItemVariant> Fjr_o1) const
{
- std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
- if (not i_mesh) {
- throw NormalError("discrete functions are not defined on the same mesh ici3");
+ std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
@@ -1198,13 +1599,16 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
}
return this->apply_fluxes(dt, //
- *i_mesh->get<MeshType>(), //
+ *mesh_v->get<MeshType>(), //
rho->get<DiscreteScalarFunction>(), //
u->get<DiscreteVectorFunction>(), //
E->get<DiscreteScalarFunction>(), //
CG->get<DiscreteTensorFunction>(), //
- ur->get<NodeValue<const Rd>>(), //
- Fjr->get<NodeValuePerCell<const Rd>>());
+ ur->get<NodeValue<const Rd>>(),
+ ur_o1->get<NodeValue<const Rd>>(), //
+ Fjr->get<NodeValuePerCell<const Rd>>(),
+ Fjr_o1->get<NodeValuePerCell<const Rd>>()
+ );
}
std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1213,28 +1617,31 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(const double& dt,
- const MeshType& mesh,
- const DiscreteScalarFunction& rho,
- const DiscreteVectorFunction& u,
- const DiscreteScalarFunction& E,
- const DiscreteTensorFunction& CG,
- const DiscreteScalarFunction& rho_app,
- const DiscreteTensorFunction& CG_app,
- const NodeValue<const Rd>& ur,
- const NodeValuePerCell<const Rd>& Fjr) const
+ const MeshType& mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ const DiscreteTensorFunction& CG,
+ const DiscreteScalarFunction& rho_app,
+ const DiscreteTensorFunction& CG_app,
+ const NodeValue<const Rd>& ur_o2,
+ const NodeValue<const Rd>& ur_o1,
+ const NodeValuePerCell<const Rd>& Fjr_o2,
+ const NodeValuePerCell<const Rd>& Fjr_o1) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ NodeValue<Rd> ur = copy(ur_o2);
+ NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
+
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
- NodeValue<Rd> new_xr = copy(mesh.xr());
- parallel_for(
- mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
- std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+ 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);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
@@ -1244,6 +1651,13 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
+ CellValue<double> new_rho_stencil = copy(rho.cellValues());
+ CellValue<Rd> new_u_stencil = copy(u.cellValues());
+ CellValue<double> new_E_stencil = copy(E.cellValues());
+ CellValue<Rdxd> new_CG_stencil = copy(CG.cellValues());
+
+ CellValue<int> cell_flag{mesh.connectivity()};
+
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
@@ -1264,8 +1678,75 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
+
+ const double& new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
+ if(new_eps < 0.){
+ cell_flag[j] = 1;
+ }else{
+ cell_flag[j] = 0;
+ }
});
+ for(CellId j = 0; j < mesh.numberOfCells(); ++j){
+ if(cell_flag[j] == 1){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ Rdxd gradv = zero;
+ std::vector<NodeId> NodeList;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ NodeId r = cell_nodes[R];
+ NodeList.push_back(r);
+ Fjr(j,R) = Fjr_o1(j,R);
+ ur[r] = ur_o1[r];
+ gradv += tensorProduct(ur[r], Cjr(j, R));
+ momentum_fluxes += Fjr(j, R);
+ energy_fluxes += dot(Fjr(j, R), ur[r]);
+ }
+ Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
+ double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
+ new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
+ new_CG[j] = new_CG_stencil[j] + dt_over_Mj * cauchy_green_fluxes;
+ new_CG[j] += transpose(new_CG[j]);
+ new_CG[j] *= 0.5;
+
+ auto cell_stencil = stencil[j];
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ CellId J = cell_stencil[i_cell];
+ auto i_cell_nodes = cell_to_node_matrix[J];
+ momentum_fluxes = zero;
+ energy_fluxes = 0;
+ gradv = zero;
+ for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
+ NodeId i_node = i_cell_nodes[R];
+ for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
+ if(NodeList[node_test] == i_node){
+ Fjr(J,R) = Fjr_o1(J,R);
+ }
+ }
+ gradv += tensorProduct(ur[i_node], Cjr(J, R));
+ momentum_fluxes += Fjr(J, R);
+ energy_fluxes += dot(Fjr(J, R), ur[i_node]);
+ }
+ cauchy_green_fluxes = rho_app[J] * (gradv * CG_app[J] + CG_app[J] * transpose(gradv));
+ dt_over_Mj = dt / (rho[J] * Vj[J]);
+ new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
+ new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
+ new_CG[J] = new_CG_stencil[J] + dt_over_Mj * cauchy_green_fluxes;
+ new_CG[J] += transpose(new_CG[J]);
+ new_CG[J] *= 0.5;
+ }
+ }
+ }
+
+ NodeValue<Rd> new_xr = copy(mesh.xr());
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
@@ -1284,34 +1765,38 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(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::shared_ptr<const DiscreteFunctionVariant>& rho_app,
- const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
- const std::shared_ptr<const ItemValueVariant>& ur,
- const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+ 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>& rho_app,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
+ const std::shared_ptr<const ItemValueVariant>& ur,
+ const std::shared_ptr<const ItemValueVariant>& ur_o1,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_o1) const
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
- throw NormalError("discrete functions are not defined on the same mesh ici4");
+ throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
- return this->order2_apply_fluxes(dt, //
- *mesh_v->get<MeshType>(), //
- rho->get<DiscreteScalarFunction>(), //
- u->get<DiscreteVectorFunction>(), //
- E->get<DiscreteScalarFunction>(), //
- CG->get<DiscreteTensorFunction>(), //
- rho_app->get<DiscreteScalarFunction>(), //
- CG_app->get<DiscreteTensorFunction>(), //
- ur->get<NodeValue<const Rd>>(), //
- Fjr->get<NodeValuePerCell<const Rd>>());
+ return this->order2_apply_fluxes(dt, //
+ *mesh_v->get<MeshType>(), //
+ rho->get<DiscreteScalarFunction>(), //
+ u->get<DiscreteVectorFunction>(), //
+ E->get<DiscreteScalarFunction>(), //
+ CG->get<DiscreteTensorFunction>(), //
+ rho_app->get<DiscreteScalarFunction>(), //
+ CG_app->get<DiscreteTensorFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
+ ur_o1->get<NodeValue<const Rd>>(),
+ Fjr->get<NodeValuePerCell<const Rd>>(),
+ Fjr_o1->get<NodeValuePerCell<const Rd>>());
}
std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1373,173 +1858,155 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
}
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 size_t& law,
- 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 chi_solid,
- const double& lambda,
- const double& mu,
- const double& gamma,
- const double& p_inf) 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;
-
- double sum_dt = 0;
-
- while(sum_dt < dt1){
-
- DiscreteScalarFunction new_aL_d = new_aL->get<DiscreteScalarFunction>();
- DiscreteScalarFunction new_aT_d = new_aT->get<DiscreteScalarFunction>();
- 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 auto [sigma_app, aL_app, aT_app] = _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid,lambda,mu,gamma,p_inf);
-
- 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);
-
- std::tie(new_sigma,new_aL,new_aT) = _apply_state_law(law,new_rho,new_u,new_E,new_CG,chi_solid,lambda,mu,gamma,p_inf);
-
- 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>,
- 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
+ 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 size_t& chi_solid,
+ const double& lambda,
+ const double& mu,
+ const double& gamma,
+ const double& p_inf) 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;
-
- auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+ std::shared_ptr sub_m = getCommonMesh({rho, u});
+ 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;
- 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;
+ double sum_dt = 0;
+ while(sum_dt < dt){
+ const auto& [sigma, aL, aT, p] = _apply_state_law(law,sub_rho,sub_u,sub_E,sub_CG,chi_solid,lambda,mu,gamma,p_inf);
- size_t fluid = 0;
- if ((mu == 0) and (lambda == 0)) {
- fluid = 1;
- }
+ 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);
- for (size_t i = 0; i < q; i++) {
- if (i == q - 1) {
- dt2 = dt1 - sum_dt;
+ double sub_dt = CFL * hyperelastic_dt(c);
+ if(sum_dt + sub_dt > dt){
+ sub_dt = dt - 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);
+ auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, p, bc_descriptor_list, CR_ur, CR_Fjr, map2);
+ auto [sub_ur2_o1, sub_Fjr2_o1] = compute_fluxes(solver_type, sub_rho, aL, aT, sub_u, sigma, bc_descriptor_list, CR_ur, CR_Fjr, map2);
- const DiscreteScalarFunction& aT_d = sqrt(mu / rho_d);
+ 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_ur2_o1, sub_Fjr2, sub_Fjr2_o1);
- 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;
+ sum_dt += sub_dt;
}
- 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};
+ 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>,
+ // 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
+ //{
+ // 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;
+//
+ // auto [map1, map2] = _computeMapping(m1, new_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;
+ // }
+//
+ // 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};
+ //}
+
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1567,6 +2034,8 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
const double& lambda_f,
@@ -1581,79 +2050,51 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
std::shared_ptr m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
- 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 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;
+ std::shared_ptr<const DiscreteFunctionVariant> p2_app = p2;
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;
- }
-
- 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_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,p1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,p2,bc_descriptor_list2,map1,map2);
+ auto [ur1_o1, Fjr1_o1, ur2_o1, Fjr2_o1, CR_ur_o1, CR_Fjr_o1] = compute_fluxes(solver_type,rho1,aL1,aT1,u1,sigma1,bc_descriptor_list1,rho2,aL2,aT2,u2,sigma2,bc_descriptor_list2,map1,map2);
- double sum_dt = dt2;
+ //auto [m1_app, rho1_app, u1_app, E1_app, CG1_app] = apply_fluxes(dt1/2., rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
size_t chi_solid = 0;
- if((mu_s!=0) and (lambda_s!=0)){
+ if((mu_s!=0) or (lambda_s!=0)){
chi_solid = 1;
}
const size_t chi_fluid = 1-chi_solid;
- while(sum_dt < dt1/2){ //Partie a revoir avec loi de comportement diff
- std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
-
- DiscreteScalarFunction aL2_d_app = aL2_app->get<DiscreteScalarFunction>();
- DiscreteScalarFunction aT2_d_app = aT2_app->get<DiscreteScalarFunction>();
- 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, rho2_app, aL2_app, aT2_app, u2_app, sigma2_app,
- bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+ auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = compute_sub_iterations(solver_type,law,0.04,dt1,rho2,u2,E2,CG2,bc_descriptor_list2,CR_ur,CR_Fjr,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
- 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;
- }
+ //const auto& [sigma1_app, aL1_app, aT1_app, p1_app] = _apply_state_law(law,rho1_app,u1_app,E1_app,CG1_app,chi_fluid,lambda_f,mu_f,gamma_f,p_inf_f);
+ //std::tie(sigma2_app, aL2_app, aT2_app, p2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
- const auto& [sigma1_app, aL1_app, aT1_app] = _apply_state_law(law,rho1_app,u1_app,E1_app,CG1_app,chi_fluid,lambda_f,mu_f,gamma_f,p_inf_f);
- std::tie(sigma2_app, aL2_app, aT2_app) = _apply_state_law(law,rho2_app,u2_app,E2_app,CG2_app,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ //auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CR_ur_app, CR_Fjr_app] = compute_fluxes(solver_type,rho1_app, aL1_app, aT1_app, u1_app, sigma1_app, p1_app, bc_descriptor_list1,
+ // rho2_app, aL2_app, aT2_app, u2_app, sigma2_app, p2_app, bc_descriptor_list2, map1, map2);
+ //auto [ur1_app_o1, Fjr1_app_o1, ur2_app_o1, Fjr2_app_o1, CR_ur_app_o1, CR_Fjr_app_o1] = compute_fluxes(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);
- 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, ur1_app_o1, Fjr1_app, Fjr1_app_o1);
+ //auto [new_m2, new_rho2, new_u2, new_E2, new_CG2] = _compute_classic_midpoint_method(solver_type,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,p2,bc_descriptor_list2,
+ // CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
- 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,law,dt1,rho2,u2,E2,CG2,aL2,aT2,sigma2,bc_descriptor_list2,
- CR_ur_app,CR_Fjr_app,map2,chi_solid,lambda_s,mu_s,gamma_s,p_inf_s);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, ur1_o1, Fjr1, Fjr1_o1);
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/Order2LocalDtHyperelasticSolver.hpp b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
index da09ae6138414ed0a8adcf6999b4789d9549be36..b0682e12a0875449ed3ac4dbf762a56a8ca9288f 100644
--- a/src/scheme/Order2LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/Order2LocalDtHyperelasticSolver.hpp
@@ -26,61 +26,37 @@ class Order2LocalDtHyperelasticSolverHandler
private:
struct IOrder2LocalDtHyperelasticSolver
{
- virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>>
- compute_fluxes(
- const SolverType& solver_type,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& aL,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) 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>>
- apply_fluxes(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::shared_ptr<const ItemValueVariant>& ur,
- 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;
+ //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;
virtual std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -109,6 +85,8 @@ class Order2LocalDtHyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
const double& lambda_f,