Exploring energy landscape to local simple shear perturbation (#24)

.gitignore

@@ -1,6 +1,10 @@
 # Build directory
 build
 
+# Python build
+*.egg-info
+dist
+
 # Prerequisites
 *.d
 

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<double, 2> plastic_CurrentYieldRight() const override; // yield strain 'right'
  xt::xtensor<size_t, 2> plastic_CurrentIndex() const override; // current index in the landscape
 
+ 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
+
+ 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

@@ -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<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;
+}
+
 } // namespace UniformSingleLayer2d
 } // namespace FrictionQPotFEM
 

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