[BREAKING CHANGE] eventDriven: adding iterative search (replaces "fal…

…lback" that is immediately deprecated)

include/FrictionQPotFEM/Generic2d.h

@@ -499,11 +499,10 @@ class System {
  but of that element the displacement update is the maximum of the element, such that all
  integration points of the element are forced to yield.
 
- \param fallback
- If `true` an elastic step (`kick = false`) is skipped if the outcome of applying the
- perturbation is too ambiguous. This can happen if the perturbation is not analytical.
- This way the next step, where `kick = true`, changes the state slightly to ensure the
- integrity of the event driven protocol.
+ \param iterative
+ If `true` the step is iteratively searched.
+ This is more costly by recommended, if the perturbation is non-analytical
+ (and contains numerical errors).
 
  \return
  Factor with which the displacement perturbation, see eventDriven_deltaU(), is scaled.
@@ -513,7 +512,7 @@ class System {
  bool kick,
  int direction = 1,
  bool yield_element = false,
- bool fallback = false);
+ bool iterative = false);
 
  /**
  Make a time-step: apply velocity-Verlet integration.

include/FrictionQPotFEM/Generic2d.hpp

@@ -594,7 +594,7 @@ inline double System::eventDrivenStep(
  bool kick,
  int direction,
  bool yield_element,
- bool fallback)
+ bool iterative)
 {
  FRICTIONQPOTFEM_ASSERT(xt::has_shape(m_pert_u, m_u.shape()));
  FRICTIONQPOTFEM_ASSERT(direction == 1 || direction == -1);
@@ -604,7 +604,7 @@ inline double System::eventDrivenStep(
  auto epsy_l = this->plastic_CurrentYieldLeft();
  auto epsy_r = this->plastic_CurrentYieldRight();
 
- FRICTIONQPOTFEM_WIP(direction > 0 || !xt::any(xt::equal(idx, 0)));
+ FRICTIONQPOTFEM_WIP(iterative || direction > 0 || !xt::any(xt::equal(idx, 0)));
 
  xt::xtensor<double, 2> target;
 
@@ -643,43 +643,166 @@ inline double System::eventDrivenStep(
  }
  }
 
- auto index = xt::unravel_index(xt::argmin(xt::abs(scale))(), scale.shape());
- size_t e = index[0];
- size_t q = index[1];
+ double ret;
+ size_t e;
+ size_t q;
 
- if (kick && yield_element) {
- q = xt::argmax(xt::view(xt::abs(scale), e, xt::all()))();
- }
+ if (!iterative) {
 
- double c = scale(e, q);
+ auto index = xt::unravel_index(xt::argmin(xt::abs(scale))(), scale.shape());
+ e = index[0];
+ q = index[1];
 
- if ((direction > 0 && c < 0) || (direction < 0 && c > 0)) {
- if (!kick) {
- return 0.0;
+ if (kick && yield_element) {
+ q = xt::argmax(xt::view(xt::abs(scale), e, xt::all()))();
  }
- else {
- FRICTIONQPOTFEM_REQUIRE((direction > 0 && c > 0) || (direction < 0 && c < 0));
+
+ ret = scale(e, q);
+
+ if ((direction > 0 && ret < 0) || (direction < 0 && ret > 0)) {
+ if (!kick) {
+ return 0.0;
+ }
+ else {
+ FRICTIONQPOTFEM_REQUIRE((direction > 0 && ret > 0) || (direction < 0 && ret < 0));
+ }
  }
  }
 
- this->setU(m_u + c * m_pert_u);
+ else {
+
+ double dir = static_cast<double>(direction);
+ auto steps = xt::sort(xt::ravel(xt::eval(xt::abs(scale))));
+ auto u_n = this->u();
+ auto jdx = xt::cast<long>(this->plastic_CurrentIndex());
+ size_t i;
+ long nip = static_cast<long>(m_nip);
+
+ // find first step that:
+ // if (!kick || (kick && !yield_element)): is plastic for the first time
+ // if (kick && yield_element): yields the element for the first time
+
+ for (i = 0; i < steps.size(); ++i) {
+ this->setU(u_n + dir * steps(i) * m_pert_u);
+ auto jdx_new = xt::cast<long>(this->plastic_CurrentIndex());
+ auto S = xt::abs(jdx_new - jdx);
+ if (!yield_element || !kick) {
+ if (xt::any(S > 0)) {
+ break;
+ }
+ }
+ else {
+ if (xt::any(xt::equal(xt::sum(S > 0, 1), nip))) {
+ break;
+ }
+ }
+ }
 
- auto idx_new = this->plastic_CurrentIndex();
- auto eps_new = GMatElastoPlasticQPot::Cartesian2d::Epsd(this->plastic_Eps())(e, q);
- auto eps_target = target(e, q);
+ // iterate to acceptable step
+
+ double right = steps(i);
+ double left = 0.0;
+ ret = right;
+
+ if (i > 0) {
+ left = steps(i - 1);
+ }
+
+ // iterate to actual step
+
+ for (size_t iiter = 0; iiter < 1100; ++iiter) {
+
+ ret = 0.5 * (right + left);
+ this->setU(u_n + dir * ret * m_pert_u);
+ auto jdx_new = xt::cast<long>(this->plastic_CurrentIndex());
+ auto S = xt::abs(jdx_new - jdx);
+
+ if (!kick) {
+ if (xt::any(S > 0)) {
+ right = ret;
+ }
+ else {
+ left = ret;
+ }
+ }
+ else if (yield_element) {
+ if (xt::any(xt::equal(xt::sum(S > 0, 1), nip))) {
+ right = ret;
+ }
+ else {
+ left = ret;
+ }
+ }
+ else {
+ if (xt::any(S > 0)) {
+ right = ret;
+ }
+ else {
+ left = ret;
+ }
+ }
+
+ if ((right - left) / left < 1e-5) {
+ break;
+ }
+ FRICTIONQPOTFEM_REQUIRE(iiter < 1000);
+ }
+
+ // final assertion: make sure that "left" and "right" are still bounds
+
+ {
+ this->setU(u_n + dir * left * m_pert_u);
+ auto jdx_new = xt::cast<long>(this->plastic_CurrentIndex());
+ auto S = xt::abs(jdx_new - jdx);
 
- if (fallback) {
- if (!kick && xt::any(xt::not_equal(idx, idx_new))) {
- this->setU(m_u - c * m_pert_u);
- return 0.0;
+ FRICTIONQPOTFEM_REQUIRE(kick || xt::all(xt::equal(S, 0)));
+ if (kick && yield_element) {
+ FRICTIONQPOTFEM_REQUIRE(!xt::any(xt::equal(xt::sum(S > 0, 1), nip)));
+ }
+ else if (kick) {
+ FRICTIONQPOTFEM_REQUIRE(xt::all(xt::equal(S, 0)));
+ }
  }
+ {
+ this->setU(u_n + dir * right * m_pert_u);
+ auto jdx_new = xt::cast<long>(this->plastic_CurrentIndex());
+ auto S = xt::abs(jdx_new - jdx);
+ FRICTIONQPOTFEM_REQUIRE(!xt::all(xt::equal(S, 0)));
+
+ FRICTIONQPOTFEM_REQUIRE(kick || !xt::all(xt::equal(S, 0)));
+ if (kick && yield_element) {
+ FRICTIONQPOTFEM_REQUIRE(xt::any(xt::equal(xt::sum(S > 0, 1), nip)));
+ }
+ else if (kick) {
+ FRICTIONQPOTFEM_REQUIRE(xt::any(S > 0));
+ }
+ }
+
+ // "output"
+
+ if (!kick) {
+ ret = dir * left;
+ }
+ else {
+ ret = dir * right;
+ }
+ this->setU(u_n);
+ FRICTIONQPOTFEM_REQUIRE((direction > 0 && ret >= 0) || (direction < 0 && ret <= 0));
  }
 
- FRICTIONQPOTFEM_REQUIRE(xt::allclose(eps_new, eps_target));
+ this->setU(m_u + ret * m_pert_u);
+
+ auto idx_new = this->plastic_CurrentIndex();
  FRICTIONQPOTFEM_REQUIRE(kick || xt::all(xt::equal(idx, idx_new)));
  FRICTIONQPOTFEM_REQUIRE(!kick || xt::any(xt::not_equal(idx, idx_new)));
 
- return c;
+ if (!iterative) {
+ auto eps_new = GMatElastoPlasticQPot::Cartesian2d::Epsd(this->plastic_Eps())(e, q);
+ auto eps_target = target(e, q);
+ FRICTIONQPOTFEM_REQUIRE(xt::allclose(eps_new, eps_target));
+ }
+
+ return ret;
 }
 
 inline void System::timeStep()

python/main.cpp

@@ -230,7 +230,7 @@ PYBIND11_MODULE(_FrictionQPotFEM, mod)
  py::arg("kick"),
  py::arg("direction") = +1,
  py::arg("yield_element") = false,
- py::arg("fallback") = false);
+ py::arg("iterative") = false);
 
  cls.def("timeStep", &SUB::System::timeStep, "timeStep");
  cls.def("timeSteps", &SUB::System::timeSteps, "timeSteps", py::arg("n"));

tests/Generic2d.py

@@ -132,31 +132,46 @@ def test_eventDrivenSimpleShear_random(self):
  kicks[1::2] = True
 
  for inc, kick in enumerate(kicks):
- idx_n = system.plastic_CurrentIndex().astype(int)
- system.eventDrivenStep(deps, kick)
- idx = system.plastic_CurrentIndex().astype(int)
+ idx_n = system.plastic_CurrentIndex()
+ u_n = system.u()
 
- if kick and inc == 0:
- self.assertTrue(np.sum(idx - idx_n) > 1)
- elif kick:
- self.assertTrue(np.sum(idx - idx_n) == 4)
+ system.eventDrivenStep(deps, kick, iterative=True)
+ idx = system.plastic_CurrentIndex()
+ if kick:
+ self.assertTrue(not np.all(idx == idx_n))
+ else:
+ self.assertTrue(np.all(idx == idx_n))
+
+ system.setU(u_n)
+ system.eventDrivenStep(deps, kick)
+ idx = system.plastic_CurrentIndex()
+ if kick:
+ self.assertTrue(not np.all(idx == idx_n))
  else:
  self.assertTrue(np.all(idx == idx_n))
 
  for kick in kicks:
- idx_n = system.plastic_CurrentIndex().astype(int)
+ idx_n = system.plastic_CurrentIndex()
+ u_n = system.u()
+
+ system.setU(u_n)
+ system.eventDrivenStep(deps, kick, -1, iterative=True)
+ idx = system.plastic_CurrentIndex()
+ if kick:
+ self.assertTrue(not np.all(idx == idx_n))
+ else:
+ self.assertTrue(np.all(idx == idx_n))
 
  if np.any(idx_n == 0):
  with self.assertRaises(IndexError):
  system.eventDrivenStep(deps, kick, -1)
  break
 
+ system.setU(u_n)
  system.eventDrivenStep(deps, kick, -1)
-
- idx = system.plastic_CurrentIndex().astype(int)
-
+ idx = system.plastic_CurrentIndex()
  if kick:
- self.assertTrue(np.sum(idx - idx_n) == -4)
+ self.assertTrue(not np.all(idx == idx_n))
  else:
  self.assertTrue(np.all(idx == idx_n))