From 56905a81175c3ccf880447ca886b3ee39363f330 Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Fri, 27 Nov 2020 11:41:58 +0100 Subject: [PATCH 1/5] Getting energy landscape --- .../FrictionQPotFEM/UniformSingleLayer2d.h | 6 + .../FrictionQPotFEM/UniformSingleLayer2d.hpp | 127 ++++++++++++++++++ python/main.cpp | 32 ++++- 3 files changed, 164 insertions(+), 1 deletion(-) diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index 3158c2d65..1afb55004 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.h +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.h @@ -311,6 +311,12 @@ 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); + protected: // mesh parameters diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp index 84b403e74..cf1d68a20 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp @@ -794,6 +794,133 @@ inline void HybridSystem::computeForceMaterial() xt::noalias(m_fmaterial) = m_felas + m_fplas; } +inline xt::xtensor HybridSystem::plastic_ElementEnergyLandscapeForSimpleShear( + const xt::xtensor& dgamma, + bool tilted) +{ + auto Eps0 = m_Eps_plas; + auto dgamma_inc = xt::diff(dgamma); + auto dV_plas = xt::view(m_quad.dV(), xt::keep(m_elem_plas), xt::all()); + FRICTIONQPOTFEM_ASSERT(xt::all(dgamma_inc >= 0.0)); + + xt::xtensor ret = xt::empty({m_nelem_plas, dgamma.size()}); + xt::view(ret, xt::all(), 0) = xt::sum(m_material_plas.Energy() * dV_plas, 1); + + for (size_t i = 0; i < dgamma_inc.size(); ++i) { + xt::view(m_Eps_plas, xt::all(), xt::all(), 0, 1) += 0.5 * dgamma_inc(i); + xt::view(m_Eps_plas, xt::all(), xt::all(), 1, 0) += 0.5 * dgamma_inc(i); + m_material_plas.setStrain(m_Eps_plas); + xt::view(ret, xt::all(), i + 1) = xt::sum(m_material_plas.Energy() * dV_plas, 1); + } + + if (tilted) { + for (size_t e = 0; e < m_nelem_plas; ++e) { + auto elem = xt::view(m_conn, m_elem_plas(e), xt::all()); + double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); + 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()) -= (0.5 * h * f * dgamma); + } + } + + // normalise to energy density + FRICTIONQPOTFEM_ASSERT(xt::allclose(dV_plas, dV_plas(0, 0))); + ret /= (dV_plas(0, 0) * static_cast(m_nip)); + + // restore strain: internally strain is always assumed to be match the current displacement field + m_Eps_plas = Eps0; + m_material_plas.setStrain(m_Eps_plas); + + return ret; +} + +inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) +{ + static constexpr size_t nip = 4; + FRICTIONQPOTFEM_ASSERT(m_nip == nip); + + // TODO: assert on h + auto elem = xt::view(m_conn, m_elem_plas(0), xt::all()); + double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); + + // TODO: assert on V + auto dV_plas = xt::view(m_quad.dV(), xt::keep(m_elem_plas), xt::all()); + double dV = dV_plas(0, 0); + + double inf = std::numeric_limits::infinity(); + xt::xtensor ret = inf * xt::ones({m_N, size_t(2)}); + + std::array models; + xt::xtensor dgammas = 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; + + for (size_t i = 0; i < 1000; ++i) { + + 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; + } + dgammas(q) = -Eps(0, 1) + std::sqrt(std::pow(eps_new, 2.0) + std::pow(Eps(0, 1), 2.0) - std::pow(eps, 2.0)); + } + + double dgamma = xt::amin(dgammas)(); + 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); + double E = - 0.5 * h * f * dgamma; + delta_gamma += dgamma; + + std::cout << delta_gamma << ", "; + + for (size_t q = 0; q < m_nip; ++q) { + auto Eps = models[q]->Strain(); + Eps(0, 1) += dgamma + 1e-12; // TODO: as input parameter + Eps(1, 0) += dgamma + 1e-12; + models[q]->setStrain(Eps); + E += models[q]->energy() * dV; + } + + if (E < E_n) { + ret(e, 0) = delta_gamma; + ret(e, 1) = (E_n - E0) / (dV * static_cast(m_nip)); + break; + } + + E_n = E; + } + + std::cout << std::endl; + } + + m_material_plas.setStrain(m_Eps_plas); + + return ret; +} + +// inline std::tuple, xt::xtensor> +// HybridSystem::plastic_ElementEnergyLandscapeForSimpleShearEventDriven( +// size_t plastic_element, +// bool tilted) +// { +// double dv = quad.dV()(m_elem_plas(plastic), 0); + + + +// } + } // namespace UniformSingleLayer2d } // namespace FrictionQPotFEM diff --git a/python/main.cpp b/python/main.cpp index 6aa796e18..3dc4eb1f4 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,19 @@ 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) + .def("__repr__", [](const SM::HybridSystem&) { return ""; }); From 9e2f05c1da1ff01768bca1157273d015a7fa2d03 Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Fri, 27 Nov 2020 14:38:31 +0100 Subject: [PATCH 2/5] Working? --- .../FrictionQPotFEM/UniformSingleLayer2d.h | 3 +- .../FrictionQPotFEM/UniformSingleLayer2d.hpp | 102 +++++++++++------- 2 files changed, 66 insertions(+), 39 deletions(-) diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index 1afb55004..7c2f07f11 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.h +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.h @@ -315,7 +315,8 @@ class HybridSystem : public System { const xt::xtensor& dgamma, // simple shear perturbation bool tilted = true); // tilt based on current element internal force - xt::xtensor plastic_ElementEnergyBarrierForSimpleShear(bool tilted = true); + // xt::xtensor plastic_ElementEnergyBarrierForSimpleShear(bool tilted = true); + std::tuple, xt::xtensor> plastic_ElementEnergyBarrierForSimpleShear(bool tilted = true); protected: diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp index cf1d68a20..ff32b2b36 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp @@ -795,46 +795,56 @@ inline void HybridSystem::computeForceMaterial() } inline xt::xtensor HybridSystem::plastic_ElementEnergyLandscapeForSimpleShear( - const xt::xtensor& dgamma, + const xt::xtensor& Delta_gamma, bool tilted) { - auto Eps0 = m_Eps_plas; - auto dgamma_inc = xt::diff(dgamma); auto dV_plas = xt::view(m_quad.dV(), xt::keep(m_elem_plas), xt::all()); - FRICTIONQPOTFEM_ASSERT(xt::all(dgamma_inc >= 0.0)); + double dV = dV_plas(0, 0); + FRICTIONQPOTFEM_ASSERT(xt::allclose(dV_plas, 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_nelem_plas, dgamma.size()}); - xt::view(ret, xt::all(), 0) = xt::sum(m_material_plas.Energy() * dV_plas, 1); + xt::xtensor ret = xt::empty({m_nelem_plas, Delta_gamma.size()}); + xt::view(ret, xt::all(), 0) = xt::sum(m_material_plas.Energy() * dV, 1); - for (size_t i = 0; i < dgamma_inc.size(); ++i) { - xt::view(m_Eps_plas, xt::all(), xt::all(), 0, 1) += 0.5 * dgamma_inc(i); - xt::view(m_Eps_plas, xt::all(), xt::all(), 1, 0) += 0.5 * dgamma_inc(i); - m_material_plas.setStrain(m_Eps_plas); - xt::view(ret, xt::all(), i + 1) = xt::sum(m_material_plas.Energy() * dV_plas, 1); + for (size_t i = 0; i < dgamma.size(); ++i) { + xt::view(Eps, xt::all(), xt::all(), 0, 1) += 0.5 * dgamma(i); + xt::view(Eps, xt::all(), xt::all(), 1, 0) += 0.5 * 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_nelem_plas; ++e) { auto elem = xt::view(m_conn, m_elem_plas(e), xt::all()); double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); - 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()) -= (0.5 * h * f * dgamma); + 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()) -= (0.5 * h * f * Delta_gamma); } } - // normalise to energy density - FRICTIONQPOTFEM_ASSERT(xt::allclose(dV_plas, dV_plas(0, 0))); - ret /= (dV_plas(0, 0) * static_cast(m_nip)); + ret /= dV * static_cast(m_nip); - // restore strain: internally strain is always assumed to be match the current displacement field - m_Eps_plas = Eps0; m_material_plas.setStrain(m_Eps_plas); return ret; } -inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) +// inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) +inline std::tuple, xt::xtensor> HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) { + 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)); + + // TODO: input + size_t niter = 20; + static constexpr size_t nip = 4; FRICTIONQPOTFEM_ASSERT(m_nip == nip); @@ -842,15 +852,13 @@ inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpl auto elem = xt::view(m_conn, m_elem_plas(0), xt::all()); double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); - // TODO: assert on V - auto dV_plas = xt::view(m_quad.dV(), xt::keep(m_elem_plas), xt::all()); - double dV = dV_plas(0, 0); - double inf = std::numeric_limits::infinity(); xt::xtensor ret = inf * xt::ones({m_N, size_t(2)}); + xt::xtensor out_E = inf * xt::ones({m_N, niter + size_t(1)}); + xt::xtensor out_gamma = inf * xt::ones({m_N, niter + size_t(1)}); std::array models; - xt::xtensor dgammas = xt::empty({nip}); + xt::xtensor trial_dgamma = xt::empty({nip}); for (size_t e = 0; e < m_N; ++e) { @@ -864,8 +872,13 @@ inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpl double E_n = E0; - for (size_t i = 0; i < 1000; ++i) { + out_E(e, 0) = E0 / (dV * static_cast(m_nip)); + out_gamma(e, 0) = 0.0; + + for (size_t i = 0; i < niter; ++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(); @@ -875,39 +888,52 @@ inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpl if (eps < epsp) { eps_new = epsp; } - dgammas(q) = -Eps(0, 1) + std::sqrt(std::pow(eps_new, 2.0) + std::pow(Eps(0, 1), 2.0) - std::pow(eps, 2.0)); + eps_new += 1e-12; // TODO: figure out a way around this + 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)); } - double dgamma = xt::amin(dgammas)(); - 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); - double E = - 0.5 * h * f * dgamma; + // increment strain for integration points with the same value + double E = 0.0; + double dgamma = xt::amin(trial_dgamma)(); delta_gamma += dgamma; - std::cout << delta_gamma << ", "; - for (size_t q = 0; q < m_nip; ++q) { auto Eps = models[q]->Strain(); - Eps(0, 1) += dgamma + 1e-12; // TODO: as input parameter - Eps(1, 0) += dgamma + 1e-12; + 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; + } + + // energy barrier found: store the last known configuration, this will be the maximum if (E < E_n) { - ret(e, 0) = delta_gamma; + ret(e, 0) = delta_gamma - dgamma; ret(e, 1) = (E_n - E0) / (dV * static_cast(m_nip)); - break; + // break; // TODO: uncomment } + out_E(e, i + 1) = E / (dV * static_cast(m_nip)); + out_gamma(e, i + 1) = delta_gamma; + + // store 'history' E_n = E; } - - std::cout << std::endl; } m_material_plas.setStrain(m_Eps_plas); - return ret; + // return ret; + return std::make_tuple(out_E, out_gamma); } // inline std::tuple, xt::xtensor> From 86667b23ac29729663824146b34d30186cebdc84 Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Fri, 27 Nov 2020 16:17:00 +0100 Subject: [PATCH 3/5] Improving implementation --- .../FrictionQPotFEM/UniformSingleLayer2d.h | 4 + .../FrictionQPotFEM/UniformSingleLayer2d.hpp | 101 ++++++++++-------- 2 files changed, 63 insertions(+), 42 deletions(-) diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index 7c2f07f11..2b1a62e70 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 diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp index ff32b2b36..4127d4eab 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; @@ -798,72 +816,60 @@ inline xt::xtensor HybridSystem::plastic_ElementEnergyLandscapeForSim const xt::xtensor& Delta_gamma, bool tilted) { - 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)); + 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_nelem_plas, Delta_gamma.size()}); + 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) += 0.5 * dgamma(i); - xt::view(Eps, xt::all(), xt::all(), 1, 0) += 0.5 * dgamma(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_nelem_plas; ++e) { - auto elem = xt::view(m_conn, m_elem_plas(e), xt::all()); - double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); - 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()) -= (0.5 * h * f * Delta_gamma); + 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) inline std::tuple, xt::xtensor> HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) { - 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)); + size_t max_iter = 20; + double pert = 1e-12; // increment to cusps + small perturbation, to advance event driven - // TODO: input - size_t niter = 20; - - static constexpr size_t nip = 4; - FRICTIONQPOTFEM_ASSERT(m_nip == nip); - - // TODO: assert on h - auto elem = xt::view(m_conn, m_elem_plas(0), xt::all()); - double h = m_coor(elem(3), 1) - m_coor(elem(0), 1); + 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)}); - xt::xtensor out_E = inf * xt::ones({m_N, niter + size_t(1)}); - xt::xtensor out_gamma = inf * xt::ones({m_N, niter + size_t(1)}); + xt::xtensor out_E = inf * xt::ones({m_N, max_iter + size_t(1)}); + xt::xtensor out_gamma = inf * xt::ones({m_N, max_iter + size_t(1)}); + + 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; + double Delta_gamma = 0.0; for (size_t q = 0; q < m_nip; ++q) { models[q] = m_material_plas.refCusp({e, q}); @@ -871,11 +877,12 @@ inline std::tuple, xt::xtensor> HybridSystem:: } double E_n = E0; + bool negative_slope = false; out_E(e, 0) = E0 / (dV * static_cast(m_nip)); out_gamma(e, 0) = 0.0; - for (size_t i = 0; i < niter; ++i) { + 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 @@ -888,7 +895,7 @@ inline std::tuple, xt::xtensor> HybridSystem:: if (eps < epsp) { eps_new = epsp; } - eps_new += 1e-12; // TODO: figure out a way around this + 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)); @@ -897,7 +904,7 @@ inline std::tuple, xt::xtensor> HybridSystem:: // increment strain for integration points with the same value double E = 0.0; double dgamma = xt::amin(trial_dgamma)(); - delta_gamma += dgamma; + Delta_gamma += dgamma; for (size_t q = 0; q < m_nip; ++q) { auto Eps = models[q]->Strain(); @@ -908,22 +915,32 @@ inline std::tuple, xt::xtensor> HybridSystem:: } 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; + 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; + } + } + if (negative_slope) { + if (E > E_n) { + negative_slope = false; + // E0 = E_n; // TODO: check + } } // energy barrier found: store the last known configuration, this will be the maximum - if (E < E_n) { - ret(e, 0) = delta_gamma - dgamma; + if (E < E_n || !negative_slope) { + ret(e, 0) = Delta_gamma - dgamma; ret(e, 1) = (E_n - E0) / (dV * static_cast(m_nip)); // break; // TODO: uncomment } out_E(e, i + 1) = E / (dV * static_cast(m_nip)); - out_gamma(e, i + 1) = delta_gamma; + out_gamma(e, i + 1) = Delta_gamma; // store 'history' E_n = E; From 422f1d1a6ca1f2868dec320ffa877fcade70358c Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Fri, 27 Nov 2020 18:09:05 +0100 Subject: [PATCH 4/5] Finishing implementation --- .../FrictionQPotFEM/UniformSingleLayer2d.h | 6 ++- .../FrictionQPotFEM/UniformSingleLayer2d.hpp | 41 ++++--------------- python/main.cpp | 4 +- 3 files changed, 16 insertions(+), 35 deletions(-) diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index 2b1a62e70..951d53593 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.h +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.h @@ -319,8 +319,10 @@ class HybridSystem : public System { const xt::xtensor& dgamma, // simple shear perturbation bool tilted = true); // tilt based on current element internal force - // xt::xtensor plastic_ElementEnergyBarrierForSimpleShear(bool tilted = true); - std::tuple, xt::xtensor> plastic_ElementEnergyBarrierForSimpleShear(bool tilted = true); + 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: diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp index 4127d4eab..943e9c06f 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp @@ -846,21 +846,17 @@ inline xt::xtensor HybridSystem::plastic_ElementEnergyLandscapeForSim return ret; } -// inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) -inline std::tuple, xt::xtensor> HybridSystem::plastic_ElementEnergyBarrierForSimpleShear(bool tilted) +inline xt::xtensor HybridSystem::plastic_ElementEnergyBarrierForSimpleShear( + bool tilted, + size_t max_iter, + double pert) { - size_t max_iter = 20; - double pert = 1e-12; // increment to cusps + small perturbation, to advance event driven - 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)}); - xt::xtensor out_E = inf * xt::ones({m_N, max_iter + size_t(1)}); - xt::xtensor out_gamma = inf * xt::ones({m_N, max_iter + size_t(1)}); - static constexpr size_t nip = 4; FRICTIONQPOTFEM_ASSERT(m_nip == nip); std::array models; @@ -879,9 +875,6 @@ inline std::tuple, xt::xtensor> HybridSystem:: double E_n = E0; bool negative_slope = false; - out_E(e, 0) = E0 / (dV * static_cast(m_nip)); - out_gamma(e, 0) = 0.0; - for (size_t i = 0; i < max_iter; ++i) { // for each integration point: compute the increment in strain to reach @@ -921,27 +914,23 @@ inline std::tuple, xt::xtensor> HybridSystem:: } if (i == 0) { - if (E < E_n) { + if (E <= E_n) { negative_slope = true; } } - if (negative_slope) { + else if (negative_slope) { if (E > E_n) { negative_slope = false; - // E0 = E_n; // TODO: check } } // energy barrier found: store the last known configuration, this will be the maximum - if (E < E_n || !negative_slope) { + if (E < E_n && !negative_slope) { ret(e, 0) = Delta_gamma - dgamma; ret(e, 1) = (E_n - E0) / (dV * static_cast(m_nip)); - // break; // TODO: uncomment + break; } - out_E(e, i + 1) = E / (dV * static_cast(m_nip)); - out_gamma(e, i + 1) = Delta_gamma; - // store 'history' E_n = E; } @@ -949,21 +938,9 @@ inline std::tuple, xt::xtensor> HybridSystem:: m_material_plas.setStrain(m_Eps_plas); - // return ret; - return std::make_tuple(out_E, out_gamma); + return ret; } -// inline std::tuple, xt::xtensor> -// HybridSystem::plastic_ElementEnergyLandscapeForSimpleShearEventDriven( -// size_t plastic_element, -// bool tilted) -// { -// double dv = quad.dV()(m_elem_plas(plastic), 0); - - - -// } - } // namespace UniformSingleLayer2d } // namespace FrictionQPotFEM diff --git a/python/main.cpp b/python/main.cpp index 3dc4eb1f4..d2f7766e5 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -121,7 +121,9 @@ PYBIND11_MODULE(FrictionQPotFEM, m) "plastic_ElementEnergyBarrierForSimpleShear", &SM::HybridSystem::plastic_ElementEnergyBarrierForSimpleShear, "plastic_ElementEnergyBarrierForSimpleShear", - py::arg("titled") = true) + py::arg("titled") = true, + py::arg("max_iter") = 100, + py::arg("perturbation") = 1e-12) .def("__repr__", [](const SM::HybridSystem&) { return ""; From 0897c5ef7efcf7bd72fc724228f17904887d4b8b Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Fri, 27 Nov 2020 18:09:46 +0100 Subject: [PATCH 5/5] Ignoring Python build --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) 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