diff --git a/.gitignore b/.gitignore index 3fd142de5..fc0b6ef93 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ # Build directory build +# Python build +*.egg-info +dist + # Prerequisites *.d diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index 3158c2d65..951d53593 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.h +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.h @@ -90,6 +90,10 @@ class System { auto coor() const; // nodal coordinates auto dofs() const; // DOFs per node + // Basic properties of the layer + double plastic_h() const; // element height + double plastic_dV() const; // integration point volume + // Extract nodal quantities. auto u() const; // displacements auto fmaterial() const; // material resistance @@ -311,6 +315,15 @@ class HybridSystem : public System { xt::xtensor plastic_CurrentYieldRight() const override; // yield strain 'right' xt::xtensor plastic_CurrentIndex() const override; // current index in the landscape + xt::xtensor plastic_ElementEnergyLandscapeForSimpleShear( + const xt::xtensor& dgamma, // simple shear perturbation + bool tilted = true); // tilt based on current element internal force + + xt::xtensor plastic_ElementEnergyBarrierForSimpleShear( + bool tilted = true, // tilt based on current element internal force + size_t max_iter = 100, // maximum number of iterations to find a barrier + double perturbation = 1e-12); // small perturbation, to advance event driven + protected: // mesh parameters diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp index 84b403e74..943e9c06f 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp @@ -263,6 +263,24 @@ inline auto System::u() const return m_u; } +inline double System::plastic_h() const +{ + auto bot = xt::view(m_conn, xt::keep(m_elem_plas), 0); + auto top = xt::view(m_conn, xt::keep(m_elem_plas), 3); + auto h_plas = xt::view(m_coor, xt::keep(top), 1) - xt::view(m_coor, xt::keep(bot), 1); + double h = h_plas(0); + FRICTIONQPOTFEM_ASSERT(xt::allclose(h_plas, h)); + return h; +} + +inline double System::plastic_dV() const +{ + auto dV_plas = xt::view(m_quad.dV(), xt::keep(m_elem_plas), xt::all()); + double dV = dV_plas(0, 0); + FRICTIONQPOTFEM_ASSERT(xt::allclose(dV_plas, dV)); + return dV; +} + inline auto System::fmaterial() const { return m_fmaterial; @@ -794,6 +812,135 @@ inline void HybridSystem::computeForceMaterial() xt::noalias(m_fmaterial) = m_felas + m_fplas; } +inline xt::xtensor HybridSystem::plastic_ElementEnergyLandscapeForSimpleShear( + const xt::xtensor& Delta_gamma, + bool tilted) +{ + double h = this->plastic_h(); + double dV = this->plastic_dV(); + + auto Eps = m_Eps_plas; + auto dgamma = xt::diff(Delta_gamma); + FRICTIONQPOTFEM_ASSERT(xt::all(dgamma >= 0.0)); + + xt::xtensor ret = xt::empty({m_N, Delta_gamma.size()}); + xt::view(ret, xt::all(), 0) = xt::sum(m_material_plas.Energy() * dV, 1); + + for (size_t i = 0; i < dgamma.size(); ++i) { + xt::view(Eps, xt::all(), xt::all(), 0, 1) += dgamma(i); + xt::view(Eps, xt::all(), xt::all(), 1, 0) += dgamma(i); + m_material_plas.setStrain(Eps); + xt::view(ret, xt::all(), i + 1) = xt::sum(m_material_plas.Energy() * dV, 1); + } + + if (tilted) { + for (size_t e = 0; e < m_N; ++e) { + double f = m_fe_plas(e, 2, 0) + m_fe_plas(e, 3, 0) + - m_fe_plas(e, 0, 0) - m_fe_plas(e, 1, 0); + xt::view(ret, e, xt::all()) -= h * f * Delta_gamma; + } + } + + ret /= dV * static_cast(m_nip); + m_material_plas.setStrain(m_Eps_plas); + return ret; +} + +inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear( + bool tilted, + size_t max_iter, + double pert) +{ + double h = this->plastic_h(); + double dV = this->plastic_dV(); + + double inf = std::numeric_limits::infinity(); + xt::xtensor ret = inf * xt::ones({m_N, size_t(2)}); + + static constexpr size_t nip = 4; + FRICTIONQPOTFEM_ASSERT(m_nip == nip); + std::array models; + xt::xtensor trial_dgamma = xt::empty({nip}); + + for (size_t e = 0; e < m_N; ++e) { + + double E0 = 0.0; + double Delta_gamma = 0.0; + + for (size_t q = 0; q < m_nip; ++q) { + models[q] = m_material_plas.refCusp({e, q}); + E0 += models[q]->energy() * dV; + } + + double E_n = E0; + bool negative_slope = false; + + for (size_t i = 0; i < max_iter; ++i) { + + // for each integration point: compute the increment in strain to reach + // the next minimum or the next cusp + for (size_t q = 0; q < m_nip; ++q) { + auto Eps = models[q]->Strain(); + double epsy_r = models[q]->currentYieldRight(); + double epsp = models[q]->epsp(); + double eps = GM::Epsd(Eps)(); + double eps_new = epsy_r; + if (eps < epsp) { + eps_new = epsp; + } + eps_new += pert; + trial_dgamma(q) = + - Eps(0, 1) + + std::sqrt(std::pow(eps_new, 2.0) + std::pow(Eps(0, 1), 2.0) - std::pow(eps, 2.0)); + } + + // increment strain for integration points with the same value + double E = 0.0; + double dgamma = xt::amin(trial_dgamma)(); + Delta_gamma += dgamma; + + for (size_t q = 0; q < m_nip; ++q) { + auto Eps = models[q]->Strain(); + Eps(0, 1) += dgamma; + Eps(1, 0) += dgamma; + models[q]->setStrain(Eps); + E += models[q]->energy() * dV; + } + + if (tilted) { + double f = m_fe_plas(e, 2, 0) + m_fe_plas(e, 3, 0) + - m_fe_plas(e, 0, 0) - m_fe_plas(e, 1, 0); + E -= h * f * Delta_gamma; + } + + if (i == 0) { + if (E <= E_n) { + negative_slope = true; + } + } + else if (negative_slope) { + if (E > E_n) { + negative_slope = false; + } + } + + // energy barrier found: store the last known configuration, this will be the maximum + if (E < E_n && !negative_slope) { + ret(e, 0) = Delta_gamma - dgamma; + ret(e, 1) = (E_n - E0) / (dV * static_cast(m_nip)); + break; + } + + // store 'history' + E_n = E; + } + } + + m_material_plas.setStrain(m_Eps_plas); + + return ret; +} + } // namespace UniformSingleLayer2d } // namespace FrictionQPotFEM diff --git a/python/main.cpp b/python/main.cpp index 6aa796e18..d2f7766e5 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -28,7 +28,7 @@ PYBIND11_MODULE(FrictionQPotFEM, m) namespace SM = FrictionQPotFEM::UniformSingleLayer2d; - sm.def("versionInfo", &SM::versionInfo, "Return version information."); + // sm.def("versionInfo", &SM::versionInfo, "Return version information."); py::class_(sm, "HybridSystem") @@ -48,6 +48,18 @@ PYBIND11_MODULE(FrictionQPotFEM, m) py::arg("elem_elastic"), py::arg("elem_plastic")) + .def( + "setMassMatrix", + &SM::HybridSystem::setMassMatrix, + "setMassMatrix", + py::arg("rho_elem")) + + .def( + "setDampingMatrix", + &SM::HybridSystem::setDampingMatrix, + "setDampingMatrix", + py::arg("alpha_elem")) + .def( "setElastic", &SM::HybridSystem::setElastic, @@ -83,6 +95,11 @@ PYBIND11_MODULE(FrictionQPotFEM, m) .def("material_plastic", &SM::HybridSystem::material_plastic, "material_plastic") .def("Eps", &SM::HybridSystem::Eps, "Eps") .def("Sig", &SM::HybridSystem::Sig, "Sig") + .def("plastic_Eps", &SM::HybridSystem::plastic_Eps, "plastic_Eps") + .def("plastic_Sig", &SM::HybridSystem::plastic_Sig, "plastic_Sig") + .def("plastic_CurrentYieldLeft", &SM::HybridSystem::plastic_CurrentYieldLeft, "plastic_CurrentYieldLeft") + .def("plastic_CurrentYieldRight", &SM::HybridSystem::plastic_CurrentYieldRight, "plastic_CurrentYieldRight") + .def("plastic_CurrentIndex", &SM::HybridSystem::plastic_CurrentIndex, "plastic_CurrentIndex") .def("timeStep", &SM::HybridSystem::timeStep, "timeStep") .def( @@ -93,6 +110,21 @@ PYBIND11_MODULE(FrictionQPotFEM, m) py::arg("niter_tol") = 20, py::arg("max_iter") = 1000000) + .def( + "plastic_ElementEnergyLandscapeForSimpleShear", + &SM::HybridSystem::plastic_ElementEnergyLandscapeForSimpleShear, + "plastic_ElementEnergyLandscapeForSimpleShear", + py::arg("dgamma"), + py::arg("titled") = true) + + .def( + "plastic_ElementEnergyBarrierForSimpleShear", + &SM::HybridSystem::plastic_ElementEnergyBarrierForSimpleShear, + "plastic_ElementEnergyBarrierForSimpleShear", + py::arg("titled") = true, + py::arg("max_iter") = 100, + py::arg("perturbation") = 1e-12) + .def("__repr__", [](const SM::HybridSystem&) { return ""; });