Adding addAffineSimpleShearCentered (#41)

include/FrictionQPotFEM/UniformSingleLayer2d.h

@@ -153,10 +153,16 @@ class System {
  auto plastic_signOfPerturbation(const xt::xtensor<double, 2>& delta_u);
 
  // Add affine simple shear (may be negative to subtract affine simple shear).
+ // The displacement of the bottom boundary is zero, while it is maximal for the top boundary.
  // The input is the strain value, the shear component deformation gradient is twice that.
  // Return deformation gradient that is applied to the system.
  double addAffineSimpleShear(double delta_gamma);
 
+ // Similar to "addAffineSimpleShear" with the difference that the displacement is zero
+ // exactly in the middle, while the displacement at the bottom and the top boundary is maximal
+ // (with a negative displacement for the bottom boundary).
+ double addAffineSimpleShearCentered(double delta_gamma);
+
  // Add event driven simple shear step.
  // Return deformation gradient that is applied to the system.
  double addSimpleShearEventDriven(

include/FrictionQPotFEM/UniformSingleLayer2d.hpp

@@ -558,6 +558,21 @@ inline double System::addAffineSimpleShear(double delta_gamma)
  return delta_gamma * 2.0;
 }
 
+inline double System::addAffineSimpleShearCentered(double delta_gamma)
+{
+ size_t ll = m_conn(m_elem_plas(0), 0);
+ size_t ul = m_conn(m_elem_plas(0), 3);
+ double y0 = (m_coor(ul, 1) + m_coor(ll, 1)) / 2.0;
+ auto u_new = this->u();
+
+ for (size_t n = 0; n < m_nnode; ++n) {
+ u_new(n, 0) += 2.0 * delta_gamma * (m_coor(n, 1) - y0);
+ }
+ this->setU(u_new);
+
+ return delta_gamma * 2.0;
+}
+
 inline double System::addSimpleShearEventDriven(
  double deps_kick, bool kick, double direction, bool dry_run)
 {

test/UniformSingleLayer2d.cpp

@@ -75,6 +75,48 @@ TEST_CASE("FrictionQPotFEM::UniformSingleLayer2d", "UniformSingleLayer2d.h")
  }
  }
 
+ SECTION("System::addAffineSimpleShearCentered")
+ {
+ GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
+
+ FrictionQPotFEM::UniformSingleLayer2d::System sys(
+ mesh.coor(),
+ mesh.conn(),
+ mesh.dofs(),
+ xt::arange<size_t>(mesh.nnode() * mesh.ndim()),
+ xt::xtensor<size_t, 1>{0, 1, 2, 6, 7, 8},
+ xt::xtensor<size_t, 1>{3, 4, 5});
+
+ sys.setMassMatrix(xt::ones<double>({mesh.nelem()}));
+ sys.setDampingMatrix(xt::ones<double>({mesh.nelem()}));
+
+ sys.setElastic(
+ xt::ones<double>({6}),
+ xt::ones<double>({6}));
+
+ sys.setPlastic(
+ xt::xtensor<double, 1>{1.0, 1.0, 1.0},
+ xt::xtensor<double, 1>{1.0, 1.0, 1.0},
+ xt::xtensor<double, 2>{{1.0, 2.0, 3.0, 4.0}, {1.0, 2.0, 3.0, 4.0}, {1.0, 2.0, 3.0, 4.0}});
+
+ sys.setDt(1.0);
+
+ auto plastic = sys.plastic();
+ auto conn = sys.conn();
+ auto bot = xt::view(conn, xt::keep(plastic), 0); // missing last node, but ok for test
+ auto top = xt::view(conn, xt::keep(plastic), 3);
+
+ for (size_t i = 0; i < 10; ++i) {
+ double delta_gamma = 0.01;
+ double gamma = delta_gamma * static_cast<double>(i + 1);
+ sys.addAffineSimpleShearCentered(delta_gamma);
+ auto u = sys.u();
+ auto du = xt::eval(xt::view(u, xt::keep(top), 1) + xt::view(u, xt::keep(bot), 1));
+ REQUIRE(xt::allclose(xt::view(sys.Eps(), xt::all(), xt::all(), 0, 1), gamma));
+ REQUIRE(xt::allclose(du, 0.0));
+ }
+ }
+
  SECTION("System::addSimpleShearEventDriven")
  {
  GooseFEM::Mesh::Quad4::Regular mesh(1, 1);