Deprecated local energy barriers (#33)

include/FrictionQPotFEM/UniformSingleLayer2d.h

@@ -320,19 +320,6 @@ class HybridSystem : public System {
  xt::xtensor<double, 2> plastic_CurrentYieldRight() const override; // yield strain 'right'
  xt::xtensor<size_t, 2> plastic_CurrentIndex() const override; // current index in the landscape
 
- // Read the (tilted) energy landscape to local simple shear perturbation
- xt::xtensor<double, 2> plastic_ElementEnergyLandscapeForSimpleShear(
- const xt::xtensor<double, 1>& 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<double, 2> 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

include/FrictionQPotFEM/UniformSingleLayer2d.hpp

@@ -844,135 +844,6 @@ inline void HybridSystem::computeForceMaterial()
  xt::noalias(m_fmaterial) = m_felas + m_fplas;
 }
 
-inline xt::xtensor<double, 2> HybridSystem::plastic_ElementEnergyLandscapeForSimpleShear(
- const xt::xtensor<double, 1>& 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<double, 2> ret = xt::empty<double>({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<double>(m_nip);
- m_material_plas.setStrain(m_Eps_plas);
- return ret;
-}
-
-inline xt::xtensor<double, 2> 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<double>::infinity();
- xt::xtensor<double, 2> ret = inf * xt::ones<double>({m_N, size_t(2)});
-
- static constexpr size_t nip = 4;
- FRICTIONQPOTFEM_ASSERT(m_nip == nip);
- std::array<GM::Cusp*, nip> models;
- xt::xtensor<double, 1> trial_dgamma = xt::empty<double>({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<double>(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

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 "<FrictionQPotFEM.UniformSingleLayer2d.HybridSystem>";
  });

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<size_t, 1> iip = xt::empty<size_t>({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<size_t, 1>{0, 1, 2, 3, 5, 6, 7, 8},
- xt::xtensor<size_t, 1>{4});
-
- sys.setMassMatrix(xt::ones<double>({mesh.nelem()}));
- sys.setDampingMatrix(xt::ones<double>({mesh.nelem()}));
-
- sys.setElastic(
- 1.0 * xt::ones<double>({8}),
- 1.0 * xt::ones<double>({8}));
-
- sys.setPlastic(
- xt::xtensor<double, 1>{1.0},
- xt::xtensor<double, 1>{1.0},
- xt::xtensor<double, 2>{{1.0, 2.0, 3.0, 4.0}});
-
- sys.setDt(1.0);
-
- REQUIRE(sys.isHomogeneousElastic());
-
- xt::xtensor<double, 2> ret = {{1.0, 1.0}};
- REQUIRE(xt::allclose(ret, sys.plastic_ElementYieldBarrierForSimpleShear()));
-}
-
 SECTION("System::plastic_*")
 {
  GooseFEM::Mesh::Quad4::Regular mesh(1, 1);