diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 98d8d86cf98c26faf0975ea03c65ae3a6abe9fa5..3ed298e6e6ad59e51a123171f9ba5ddd5f0d20c9 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -223,6 +223,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
@@ -236,7 +237,7 @@ LocalFSIModule::LocalFSIModule()
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
.apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
- sigma, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
+ sigma, p, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
}
));
@@ -251,6 +252,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
@@ -266,7 +268,7 @@ LocalFSIModule::LocalFSIModule()
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
.apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
- sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
+ sigma, p, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
}
));
@@ -281,6 +283,7 @@ LocalFSIModule::LocalFSIModule()
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
@@ -298,7 +301,7 @@ LocalFSIModule::LocalFSIModule()
return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
.solver()
.apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt, rho, u, E, CG, aL, aT,
- sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
+ sigma, p, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
}
));
diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp
index c8f2d9e7bde77263421e5999584c16195a9e1afb..0c2c166ef77412a986d594f65c25373df544a1cf 100644
--- a/src/scheme/Order2AcousticSolver.cpp
+++ b/src/scheme/Order2AcousticSolver.cpp
@@ -410,7 +410,6 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
_vector_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rd>& f,
DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
- //const DiscreteFunctionP0<const double>& p,
const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const
{
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
@@ -421,30 +420,25 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
const auto xj = mesh_data.xj();
- //Calcul du gradiant
CellValue<Rd> n{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto coefficients_p = DPk_ph.coefficients(cell_id);
Rd grad_p = zero;
for(size_t i = 1; i < coefficients_p.size(); ++i){
- grad_p[i] = coefficients_p[i];
+ grad_p[i-1] = coefficients_p[i];
}
const double norm_grad_p = l2Norm(grad_p);
if(norm_grad_p == 0){
n[cell_id] = zero;
- }
- else{
- for(size_t d = 0; d < Dimension; ++d){
- n[cell_id][d] = grad_p[d] / norm_grad_p;
- }
+ }
+ else {
+ n[cell_id] = (1/norm_grad_p) * grad_p;
}
});
-
- //Construction des vecteurs orthogonaux
const CellValue<Rd> t{mesh.connectivity()};
const CellValue<Rd> l{mesh.connectivity()};
parallel_for(
@@ -466,9 +460,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
tj = zero;
}
else {
- for(size_t d = 0; d < Dimension; ++d){
- tj[d] = tj[d] / norm_tj;
- }
+ tj = (1/norm_tj) * tj;
}
t[cell_id] = tj;
@@ -478,14 +470,11 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
const double& norm_lj = l2Norm(lj);
- for(size_t d = 0; d < Dimension; ++d){
- lj[d] = lj[d] / norm_lj;
- }
+ lj = (1/norm_lj) * lj;
}
l[cell_id] = lj;
});
- //Limiteurs
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
const double fn = dot(f[cell_id],n[cell_id]);
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index ec3ec68c8914b6b359bf484a278d102e1ecf726a..6875c5005faaf096a2b61a0625a29a9ed389057f 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -335,214 +335,368 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
return u;
}
+ 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;
+ }
+
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;
+ }
+
+ 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;
+ }
+
+ 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_sigma) const
+ 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 auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
- //const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const TinyMatrix<Dimension> I = identity;
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValuePerCell<Rd> F{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);
- //
- // 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];
- // F(J, R) = -Ajr(J, R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J, R);
- // }
- // });
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) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + DPk_sigma[j](xr[cell_nodes[r]])*Cjr(j,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;
}
- void
- _tensor_limiter(const MeshType& mesh,
- const DiscreteFunctionP0<const Rdxd>& f,
- DiscreteFunctionDPk<Dimension, Rdxd>& 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;
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();
- parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- const Rdxd fj = f[cell_id];
+ CellValue<double> alph_J2{mesh.connectivity()};
- Rdxd f_min = fj;
- Rdxd f_max = fj;
+ parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+ const double fj = f[cell_id];
- const auto cell_stencil = stencil[cell_id];
- for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
- for (size_t row = 0; row < Dimension; ++row) {
- for (size_t col = 0; col < Dimension; ++col) {
- f_min(row, col) = std::min(f_min(row, col), f[cell_stencil[i_cell]](row, col));
- f_max(row, col) = std::max(f_max(row, col), f[cell_stencil[i_cell]](row, col));
- }
- }
- }
+ double f_min = fj;
+ double f_max = fj;
- Rdxd f_bar_min = fj;
- Rdxd f_bar_max = fj;
+ 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]]);
+ }
- for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
- const CellId cell_k_id = cell_stencil[i_cell];
- const Rdxd f_xk = DPk_fh[cell_id](xj[cell_k_id]);
+ double f_bar_min = fj;
+ double f_bar_max = fj;
- for (size_t row = 0; row < Dimension; ++row) {
- for (size_t col = 0; col < Dimension; ++col) {
- f_bar_min(row, col) = std::min(f_bar_min(row, col), f_xk(row, col));
- f_bar_max(row, col) = std::max(f_bar_max(row, col), f_xk(row, col));
- }
- }
- }
+ 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]);
- const double eps = 1E-14;
- Rdxd coef1;
- for (size_t row = 0; row < Dimension; ++row) {
- for (size_t col = 0; col < Dimension; ++col) {
- coef1(row, col) = 1;
- if (std::abs(f_bar_max(row, col) - fj(row, col)) > eps) {
- coef1(row, col) = (f_max(row, col) - fj(row, col)) / (f_bar_max(row, col) - fj(row, col));
- }
- }
- }
+ f_bar_min = std::min(f_bar_min, f_xk);
+ f_bar_max = std::max(f_bar_max, f_xk);
+ }
- Rdxd coef2;
- for (size_t row = 0; row < Dimension; ++row) {
- for (size_t col = 0; col < Dimension; ++col) {
- coef2(row, col) = 1;
- if (std::abs(f_bar_min(row, col) - fj(row, col)) > eps) {
- coef2(row, col) = (f_min(row, col) - fj(row, col)) / ((f_bar_min(row, col) - fj(row, col)));
- }
- }
- }
+ const double eps = 1E-14;
+ double coef1 = 1;
+ if (std::abs(f_bar_max - fj) > eps) {
+ coef1 = (f_max - fj) / ((f_bar_max - fj));
+ }
- double min_coef1 = coef1(0, 0);
- double min_coef2 = coef2(0, 0);
- for (size_t row = 0; row < Dimension; ++row) {
- for (size_t col = 0; col < Dimension; ++col) {
- min_coef1 = std::min(min_coef1, coef1(row, col));
- min_coef2 = std::min(min_coef2, coef2(row, col));
- }
- }
+ 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);
+ auto coefficients = DPk_fh.coefficients(cell_id);
- coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
- for (size_t i = 1; i < coefficients.size(); ++i) {
- coefficients[i] *= lambda;
- }
- });
+ coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ });
}
- void
+ void
_vector_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rd>& f,
- DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const
+ DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
+ const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const
{
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
symmetry_boundary_descriptor_list);
-
+
const auto xj = mesh_data.xj();
+ CellValue<Rd> n{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- const Rd fj = f[cell_id];
+ auto coefficients_p = DPk_ph.coefficients(cell_id);
+ Rd grad_p = zero;
+ for(size_t i = 1; i < coefficients_p.size(); ++i){
+ grad_p[i-1] = coefficients_p[i];
+ }
+ const double norm_grad_p = l2Norm(grad_p);
+ if(norm_grad_p == 0){
+ n[cell_id] = zero;
+ }
+ else {
+ n[cell_id] = (1/norm_grad_p) * grad_p;
+ }
- Rd f_min = fj;
- Rd 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]);
+ const CellValue<Rd> t{mesh.connectivity()};
+ const CellValue<Rd> l{mesh.connectivity()};
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const Rd nj = n[cell_id];
+ Rd a = zero;
+ if(l2Norm(nj) != 0){
+ if((nj[0] / l2Norm(nj) != 1) and (nj[0] / l2Norm(nj) != -1)){
+ a[0] = 1.;
}
+ else {
+ a[1] = 1.;
+ }
+ }
+
+ Rd tj = a - dot(a,nj) * nj;
+ const double& norm_tj = l2Norm(tj);
+ if(norm_tj == 0){
+ tj = zero;
+ }
+ else {
+ tj = (1/norm_tj) * tj;
+ }
+ t[cell_id] = tj;
+
+ Rd lj = zero;
+ if(Dimension == 3){
+ lj[0] = nj[1]*tj[2] - nj[2]*tj[1];
+ lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
+ lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
+ const double& norm_lj = l2Norm(lj);
+ if(norm_lj != 0){
+ lj = (1/norm_lj) * lj;
+ }
+
+ }
+ l[cell_id] = lj;
+ });
+
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+ const double fn = dot(f[cell_id],n[cell_id]);
+ const double ft = dot(f[cell_id],t[cell_id]);
+ const double fl = dot(f[cell_id],l[cell_id]);
+
+ double fn_min = fn;
+ double fn_max = fn;
+ double ft_min = ft;
+ double ft_max = ft;
+ double fl_min = fl;
+ double fl_max = fl;
+
+ const auto cell_stencil = stencil[cell_id];
+ for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+ const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_stencil[i_cell]]);
+ fn_min = std::min(fn_min,fn_k);
+ fn_max = std::max(fn_max,fn_k);
+
+ const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_stencil[i_cell]]);
+ ft_min = std::min(ft_min,ft_k);
+ ft_max = std::max(ft_max,ft_k);
+
+ const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_stencil[i_cell]]);
+ fl_min = std::min(fl_min,fl_k);
+ fl_max = std::max(fl_max,fl_k);
}
- Rd f_bar_min = fj;
- Rd f_bar_max = fj;
+ double fn_bar_min = fn;
+ double fn_bar_max = fn;
+ double ft_bar_min = ft;
+ double ft_bar_max = ft;
+ double fl_bar_min = fl;
+ double fl_bar_max = fl;
- for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ 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 Rd 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]);
- }
+ const double fn_xk = dot(f_xk,n[cell_k_id]);
+ fn_bar_min = std::min(fn_bar_min,fn_xk);
+ fn_bar_max = std::max(fn_bar_max,fn_xk);
+
+ const double ft_xk = dot(f_xk,t[cell_k_id]);
+ ft_bar_min = std::min(ft_bar_min,ft_xk);
+ ft_bar_max = std::max(ft_bar_max,ft_xk);
+
+ const double fl_xk = dot(f_xk,l[cell_k_id]);
+ fl_bar_min = std::min(fl_bar_min,fl_xk);
+ fl_bar_max = std::max(fl_bar_max,fl_xk);
}
const double eps = 1E-14;
- 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]));
- }
+ 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);
}
-
- 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_t = 1;
+ if(std::abs(ft_bar_max - ft) > eps){
+ coef1_t = (ft_max - ft) / (ft_bar_max - ft);
}
-
- 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_t = 1;
+ if(std::abs(ft_bar_min - ft) > eps){
+ coef2_t = (ft_min - ft) / (ft_bar_min - ft);
}
- const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
-
- auto coefficients = DPk_fh.coefficients(cell_id);
+ 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);
+ }
- coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
- for (size_t i = 1; i < coefficients.size(); ++i) {
- coefficients[i] *= lambda;
+ const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
+ const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
+ const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
+
+ auto coefficients = DPK_fh.coefficients(cell_id);
+ coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id]
+ + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id]
+ + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
+ for(size_t i = 1; i < coefficients.size(); ++i){
+ coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id]
+ + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id]
+ + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
}
- });
+ });
}
- CellValue<double>
- _scalar_limiter_coefficient(const MeshType& mesh,
- const DiscreteFunctionP0<const double>& f,
- const DiscreteFunctionDPk<Dimension, const double>& DPk_fh) const
+ 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;
@@ -552,49 +706,71 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const auto xj = mesh_data.xj();
- CellValue<double> lambda{mesh.connectivity()};
+ //A retirer + tard
+ const auto xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ //
+ 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 fj = f[cell_id];
-
- double f_min = fj;
- double f_max = fj;
+ mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+ const double inv2j = inv2[cell_id];
+
+ double inv2_min = inv2j;
+ double inv2_max = inv2j;
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]]);
+ 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]]);
}
- double f_bar_min = fj;
- double f_bar_max = fj;
+ double inv2_bar_min = inv2j;
+ double inv2_bar_max = inv2j;
- 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]);
+ 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]]));
- f_bar_min = std::min(f_bar_min, f_xk);
- f_bar_max = std::max(f_bar_max, f_xk);
+ 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;
- if (std::abs(f_bar_max - fj) > eps) {
- coef1 = (f_max - fj) / ((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 double inv2_xk = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
+//
+ // inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
+ // inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
+ //}
+
+ const double eps = 1E-14;
+ double coef1 = 1.;
+ double coef2 = 1.;
+ 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);
}
- double coef2 = 1.;
- if (std::abs(f_bar_min - fj) > eps) {
- coef2 = (f_min - fj) / ((f_bar_min - fj));
- }
+ //const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
+ const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
+ lambda[cell_id] = lambda_inv2;
- lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
- });
+ //auto coefficients = DPk_Sh.coefficients(cell_id);
+
+ //coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
+ //for (size_t i = 1; i < coefficients.size(); ++i) {
+ // coefficients[i] *= lambda_inv2;
+ //}
+ });
return lambda;
- }
+ }
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,
@@ -616,11 +792,14 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
- return {std::make_shared<DiscreteFunctionVariant>(sigma_d), std::make_shared<DiscreteFunctionVariant>(aL_d),
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(p_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
std::make_shared<DiscreteFunctionVariant>(aT_d)};
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
+ const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>>
_apply_NeoHook2(const DiscreteScalarFunction rho,
@@ -638,18 +817,21 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
const DiscreteScalarFunction& p_d =
- (1 - chi_solid) * ((gamma_f - 1) * rho * eps) + chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
+ (1 - chi_solid) * ((gamma_f - 1) * rho * eps - gamma_f * p_inf_f) + chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
const DiscreteTensorFunction& sigma_d =
chi_solid * 2 * mu / sqrt(det(CG)) * (CG - 1. / 3. * trace(CG) * I) - p_d * I;
const DiscreteScalarFunction& aL_d = sqrt(
chi_solid * (2 * mu) / rho + ((1 - chi_solid) * gamma_f * p_d + chi_solid * (gamma_s * (p_d + p_inf_s))) / rho);
const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
- return {std::make_shared<DiscreteFunctionVariant>(sigma_d), std::make_shared<DiscreteFunctionVariant>(aL_d),
+ return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
+ std::make_shared<DiscreteFunctionVariant>(p_d),
+ std::make_shared<DiscreteFunctionVariant>(aL_d),
std::make_shared<DiscreteFunctionVariant>(aT_d)};
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
+ const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>>
_apply_state_law(const size_t law,
@@ -695,7 +877,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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");
- } // - p* identity; // - chi_solid*P_zero_air_repo*identity;
+ }
if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
@@ -708,6 +890,55 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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);
+
+ NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+
+ 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>>
+ 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::shared_ptr<const DiscreteFunctionVariant>& p_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, p_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ } // - p* identity; // - chi_solid*P_zero_air_repo*identity;
+
+ if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v, p_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>();
+ 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;
for (auto&& bc_descriptor : bc_descriptor_list) {
@@ -719,102 +950,40 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
symmetry_boundary_descriptor_list);
- // 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(mesh_v, 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);
- // }
- // }
- //
- // NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ;
- // const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- //
- // parallel_for(
- // mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
- // const auto& cell_nodes = cell_to_node_matrix[j];
- // for(size_t 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;
- // }
- // }
- // }
- // }
- //);
- //
- // 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];
- //
- // for(size_t R = 0; R < cell_nodes.size(); ++R){
- // const NodeId r = cell_nodes[R];
- //
- // for(size_t l = 0; l < sigma_coef.size(); ++l){
- // const size_t& i = row[l];
- // const size_t& k = col[l];
- //
- // auto coefficients = sigma_coef[l].coefficients(j);
- //
- // 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];
- // }
- //
- // sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
- // if(i != k){
- // sigma_lim(j,R)(i,k) *= 2;
- // }
- // }
- //
- // sigma_lim(j,R) += transpose(sigma_lim(j,R));
- // sigma_lim(j,R) *= 0.5;
- // }
- //
- // }
- //);
- auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, sigma_v);
- auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
- auto DPk_sigmah = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
+ 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> sigma_lim = copy(DPk_sigmah);
- this->_vector_limiter(mesh, u, u_lim);
- this->_tensor_limiter(mesh, sigma, sigma_lim);
+ DiscreteFunctionDPk<Dimension, Rdxd> S_lim = copy(DPk_Sh);
+ DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
+
+ 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);
+
+ const NodeValue<const Rd>& xr = mesh.xr();
+ const DiscreteScalarFunction& J2 = 0.5*trace(transpose(S)*S);
+ NodeValuePerCell<double> coef{mesh.connectivity()};
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for(size_t r = 0; r < cell_nodes.size(); ++r){
+ const double J2_p = 0.5*trace(transpose(DPk_Sh[j](xr[cell_nodes[r]]))*DPk_Sh[j](xr[cell_nodes[r]]));
+ if(J2_p != 0.){
+ coef(j,r) = J2[j]/J2_p;
+ }else{
+ coef(j,r) = 0;
+ }
+ }
+ });
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
- NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, 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);
@@ -823,7 +992,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
- NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
+ NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -840,21 +1009,24 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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();
@@ -864,29 +1036,105 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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(
@@ -910,7 +1158,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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 mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
@@ -927,8 +1177,11 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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>,
@@ -944,21 +1197,24 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const DiscreteTensorFunction& CG,
const DiscreteScalarFunction& rho_app,
const DiscreteTensorFunction& CG_app,
- 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();
@@ -968,6 +1224,13 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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];
@@ -988,8 +1251,75 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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(
@@ -1015,7 +1345,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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 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) {
@@ -1035,7 +1367,9 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
rho_app->get<DiscreteScalarFunction>(), //
CG_app->get<DiscreteTensorFunction>(), //
ur->get<NodeValue<const Rd>>(), //
- Fjr->get<NodeValuePerCell<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>,
@@ -1053,6 +1387,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
const double& lambda,
@@ -1062,22 +1397,26 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
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;
+ //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);
- auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
- // auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, Fjr);
+ auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
- // std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
+ std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
- // auto [sigma_app, aL_app, aT_app] =
- // _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
+ auto [sigma_app, p_app, aL_app, aT_app] =
+ _apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
- // auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
- return apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
+ 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/Order2HyperelasticSolver.hpp b/src/scheme/Order2HyperelasticSolver.hpp
index 0051bd0d7966323dae9d69a03ac9d4b2193478b5..99f9616045268e9775be368a34b7a2f3ca81e14b 100644
--- a/src/scheme/Order2HyperelasticSolver.hpp
+++ b/src/scheme/Order2HyperelasticSolver.hpp
@@ -25,30 +25,6 @@ class Order2HyperelasticSolverHandler
private:
struct IOrder2HyperelasticSolver
{
- 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>,
@@ -64,6 +40,7 @@ class Order2HyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
const double& lambda,