From 128e7df7f98cb4f3b1d596682abc8fbe3aa6890f Mon Sep 17 00:00:00 2001 From: Tom de Geus Date: Sat, 12 Dec 2020 14:39:33 +0100 Subject: [PATCH] Deprecated local energy barriers (#33) --- .../FrictionQPotFEM/UniformSingleLayer2d.h | 13 -- .../FrictionQPotFEM/UniformSingleLayer2d.hpp | 129 ------------------ python/main.cpp | 22 --- test/UniformSingleLayer2d.cpp | 42 ------ 4 files changed, 206 deletions(-) diff --git a/include/FrictionQPotFEM/UniformSingleLayer2d.h b/include/FrictionQPotFEM/UniformSingleLayer2d.h index bb2fe3fe3..d20f8eeb4 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.h +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.h @@ -320,19 +320,6 @@ class HybridSystem : public System { xt::xtensor plastic_CurrentYieldRight() const override; // yield strain 'right' xt::xtensor plastic_CurrentIndex() const override; // current index in the landscape - // Read the (tilted) energy landscape to local simple shear perturbation - xt::xtensor plastic_ElementEnergyLandscapeForSimpleShear( - const xt::xtensor& dgamma, // simple shear perturbation - bool tilted = true); // tilt based on current element internal force - - // Read the (first) energy barrier in the tilted energy landscape in - // "plastic_ElementEnergyLandscapeForSimpleShear". - // Returns array of shape (N, 2) with columns (delta_E, delta_epsxy). - 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 84019ba80..80fa86ee5 100644 --- a/include/FrictionQPotFEM/UniformSingleLayer2d.hpp +++ b/include/FrictionQPotFEM/UniformSingleLayer2d.hpp @@ -844,135 +844,6 @@ 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; -} - inline LocalTriggerFineLayer::LocalTriggerFineLayer(const System& sys) { // Copy / allocate local variables diff --git a/python/main.cpp b/python/main.cpp index 0874ddbd9..abbc4901d 100644 --- a/python/main.cpp +++ b/python/main.cpp @@ -180,28 +180,6 @@ 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( - "plastic_ElementYieldBarrierForSimpleShear", - &SM::HybridSystem::plastic_ElementYieldBarrierForSimpleShear, - "plastic_ElementYieldBarrierForSimpleShear", - py::arg("deps_kick") = 0.0, - py::arg("iquad") = 0) - .def("__repr__", [](const SM::HybridSystem&) { return ""; }); diff --git a/test/UniformSingleLayer2d.cpp b/test/UniformSingleLayer2d.cpp index eaed729b1..14f7173a3 100644 --- a/test/UniformSingleLayer2d.cpp +++ b/test/UniformSingleLayer2d.cpp @@ -186,48 +186,6 @@ SECTION("System::triggerElementWithLocalSimpleShear") REQUIRE(xt::allclose(eps_p, 1.0 + delta_eps / 2.0)); } -SECTION("System::plastic_ElementYieldBarrierForSimpleShear") -{ - GooseFEM::Mesh::Quad4::Regular mesh(3, 3); - - auto dofs = mesh.dofs(); - auto top = mesh.nodesTopEdge(); - auto bottom = mesh.nodesBottomEdge(); - size_t nfix = top.size(); - xt::xtensor iip = xt::empty({2 * mesh.ndim() * nfix}); - xt::view(iip, xt::range(0 * nfix, 1 * nfix)) = xt::view(dofs, xt::keep(bottom), 0); - xt::view(iip, xt::range(1 * nfix, 2 * nfix)) = xt::view(dofs, xt::keep(bottom), 1); - xt::view(iip, xt::range(2 * nfix, 3 * nfix)) = xt::view(dofs, xt::keep(top), 0); - xt::view(iip, xt::range(3 * nfix, 4 * nfix)) = xt::view(dofs, xt::keep(top), 1); - - FrictionQPotFEM::UniformSingleLayer2d::System sys( - mesh.coor(), - mesh.conn(), - dofs, - iip, - xt::xtensor{0, 1, 2, 3, 5, 6, 7, 8}, - xt::xtensor{4}); - - sys.setMassMatrix(xt::ones({mesh.nelem()})); - sys.setDampingMatrix(xt::ones({mesh.nelem()})); - - sys.setElastic( - 1.0 * xt::ones({8}), - 1.0 * xt::ones({8})); - - sys.setPlastic( - xt::xtensor{1.0}, - xt::xtensor{1.0}, - xt::xtensor{{1.0, 2.0, 3.0, 4.0}}); - - sys.setDt(1.0); - - REQUIRE(sys.isHomogeneousElastic()); - - xt::xtensor ret = {{1.0, 1.0}}; - REQUIRE(xt::allclose(ret, sys.plastic_ElementYieldBarrierForSimpleShear())); -} - SECTION("System::plastic_*") { GooseFEM::Mesh::Quad4::Regular mesh(1, 1);