From 2f44457d48929e6be075c5ede3f12eaf40374f36 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 25 May 2023 15:29:06 +0200
Subject: [PATCH] add perfect plasticity into hypoelastic solver
---
src/language/modules/SchemeModule.cpp | 17 +++++++++++------
src/scheme/HypoelasticSolver.cpp | 14 +++++++++++---
src/scheme/HypoelasticSolver.hpp | 2 ++
tests/test_TinyMatrix.cpp | 2 +-
4 files changed, 25 insertions(+), 10 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index dd5010f70..9237f5db7 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -611,6 +611,7 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -618,10 +619,11 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return HypoelasticSolverHandler{getCommonMesh({rho, u, E, aL, aT, sigmad, p, mu})}
+ return HypoelasticSolverHandler{
+ getCommonMesh({rho, u, E, aL, aT, sigmad, p, mu, yieldy})}
.solver()
.apply(HypoelasticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, sigmad, aL, aT,
- p, mu, bc_descriptor_list);
+ p, mu, yieldy, bc_descriptor_list);
}
));
@@ -658,6 +660,7 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -665,10 +668,11 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return HypoelasticSolverHandler{getCommonMesh({rho, u, E, aL, aT, sigmad, mu, p})}
+ return HypoelasticSolverHandler{
+ getCommonMesh({rho, u, E, aL, aT, sigmad, mu, yieldy, p})}
.solver()
.apply(HypoelasticSolverHandler::SolverType::Glace, dt, rho, u, E, sigmad, aL, aT, p,
- mu, bc_descriptor_list);
+ mu, yieldy, bc_descriptor_list);
}
));
@@ -681,6 +685,7 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& E, //
const std::shared_ptr<const DiscreteFunctionVariant>& sigmad, //
const std::shared_ptr<const DiscreteFunctionVariant>& mu, //
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy, //
const std::shared_ptr<const ItemValueVariant>& ur, //
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, //
const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
@@ -688,9 +693,9 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return HypoelasticSolverHandler{getCommonMesh({rho, u, E, sigmad, mu})} //
+ return HypoelasticSolverHandler{getCommonMesh({rho, u, E, sigmad, mu, yieldy})} //
.solver()
- .apply_fluxes(dt, rho, u, E, sigmad, mu, ur, Fjr);
+ .apply_fluxes(dt, rho, u, E, sigmad, mu, yieldy, ur, Fjr);
}
));
diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp
index 0bea644c9..f41632e76 100644
--- a/src/scheme/HypoelasticSolver.cpp
+++ b/src/scheme/HypoelasticSolver.cpp
@@ -452,6 +452,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
const DiscreteScalarFunction& E,
const DiscreteTensorFunction& sigmad,
const DiscreteScalarFunction& mu,
+ const DiscreteScalarFunction& yieldy,
const NodeValue<const Rd>& ur,
const NodeValuePerCell<const Rd>& Fjr) const
{
@@ -475,6 +476,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_sigmad = copy(sigmad.cellValues());
+ auto chi = [&](const double a) -> double { return (a < 0 ? 0 : 1); };
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -489,17 +491,20 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
- // TODO
const Rdxd I = identity;
+ const double norms = l2Norm(sigmad[j]);
+ const double yieldf = norms - std::sqrt(2 / 3) * yieldy[j];
+ const Rdxd N = 1. / norms * sigmad[j];
const Rdxd D = 0.5 * (gradv + transpose(gradv));
const Rdxd Dd = D - 1 / 3 * trace(D) * I;
+ const Rdxd Dp = chi(yieldf) * scalarProduct(N, D) * N;
const Rdxd gradva = 0.5 * (gradv - transpose(gradv));
const Rdxd jaumann = sigmad[j] * gradva - gradva * sigmad[j];
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
- new_sigmad[j] += dt_over_Vj * (2 * mu[j] * Dd - jaumann);
+ new_sigmad[j] += dt_over_Vj * (2 * mu[j] * (Dd - Dp) - jaumann);
});
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
@@ -524,6 +529,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
{
@@ -543,6 +549,7 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
E->get<DiscreteScalarFunction>(), //
sigmad->get<DiscreteTensorFunction>(), //
mu->get<DiscreteScalarFunction>(), //
+ yieldy->get<DiscreteScalarFunction>(), //
ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}
@@ -562,10 +569,11 @@ class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolv
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigmad, p, bc_descriptor_list);
- return apply_fluxes(dt, rho, u, E, sigmad, mu, ur, Fjr);
+ return apply_fluxes(dt, rho, u, E, sigmad, mu, yieldy, ur, Fjr);
}
HypoelasticSolver() = default;
diff --git a/src/scheme/HypoelasticSolver.hpp b/src/scheme/HypoelasticSolver.hpp
index f0239743f..f67346f00 100644
--- a/src/scheme/HypoelasticSolver.hpp
+++ b/src/scheme/HypoelasticSolver.hpp
@@ -48,6 +48,7 @@ class HypoelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& sigmad,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
@@ -66,6 +67,7 @@ class HypoelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::shared_ptr<const DiscreteFunctionVariant>& mu,
+ const std::shared_ptr<const DiscreteFunctionVariant>& yieldy,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
IHypoelasticSolver() = default;
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 0af7e5931..ef7e06364 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -303,7 +303,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
TinyMatrix<3, 4, int> B(0, 0, 0, -1, 1, -1, 1, -1, 1, -1, 1, -1);
// TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
- REQUIRE(norm(B) == 3);
+ REQUIRE(l2Norm(B) == 3);
}
#ifndef NDEBUG
SECTION("output with signaling NaN")
--
GitLab