Function to facilitate applying an external force (just a new wrapper…

…, before one could also apply `fext`) (#167)

include/FrictionQPotFEM/Generic2d.h

@@ -775,6 +775,47 @@ class System {
  return m_fext;
  }
 
+ /**
+ * @brief Modify the external force such that a shear stress is applied to the top boundary.
+ * @param sigxy The shear stress to be applied.
+ */
+ void applyShearStress(double sigxy)
+ {
+ double t = m_coor(m_coor.shape(0) - 1, 1);
+ double b = m_coor(0, 1);
+ auto top = xt::sort(
+ xt::flatten_indices(xt::argwhere(xt::isclose(xt::view(m_coor, xt::all(), 1), t))));
+ auto bot = xt::sort(
+ xt::flatten_indices(xt::argwhere(xt::isclose(xt::view(m_coor, xt::all(), 1), b))));
+ double h = m_coor(top(1), 0) - m_coor(top(0), 0);
+
+ for (size_t i = 0; i < top.size() - 1; i++) {
+ FRICTIONQPOTFEM_REQUIRE(xt::allclose(m_coor(top(i + 1), 0) - m_coor(top(i), 0), h));
+ }
+
+ for (size_t d = 0; d < 2; d++) {
+ auto dofs = xt::view(m_vector.dofs(), xt::keep(top), d);
+ if (d == 0) {
+ FRICTIONQPOTFEM_REQUIRE(!xt::any(xt::in1d(dofs, m_vector.iip())));
+ }
+ else {
+ FRICTIONQPOTFEM_REQUIRE(xt::all(xt::in1d(dofs, m_vector.iip())));
+ }
+ }
+
+ for (size_t d = 0; d < 2; d++) {
+ auto dofs = xt::view(m_vector.dofs(), xt::keep(bot), d);
+ FRICTIONQPOTFEM_REQUIRE(xt::all(xt::in1d(dofs, m_vector.iip())));
+ }
+
+ m_fext(top.front(), 0) = 0.5 * sigxy * h;
+ m_fext(top.back(), 0) = 0.5 * sigxy * h;
+
+ for (size_t i = 1; i < top.size() - 1; i++) {
+ m_fext(top(i), 0) = sigxy * h;
+ }
+ }
+
 public:
  /**
  Internal force.

python/main.cpp

@@ -208,6 +208,12 @@ PYBIND11_MODULE(_FrictionQPotFEM, mod)
  cls.def_property(
  "fext", &SUB::System::fext, &SUB::System::setFext<xt::pytensor<double, 2>>, "fext");
 
+ cls.def(
+ "applyShearStress",
+ &SUB::System::applyShearStress,
+ "Apply shear stress to the top boundary",
+ py::arg("sigxy"));
+
  cls.def_property_readonly("fint", &SUB::System::fint, "fint");
  cls.def_property_readonly("fmaterial", &SUB::System::fmaterial, "fmaterial");
  cls.def_property_readonly("fdamp", &SUB::System::fdamp, "fdamp");

tests/test_Generic2d.py

@@ -64,6 +64,71 @@ def test_version_compiler(self):
 
  self.assertEqual(parsed.day, datetime.date.today().day)
 
+ def test_applyStres(self):
+ """
+ Apply a stress to the top.
+ """
+
+ mesh = GooseFEM.Mesh.Quad4.Regular(3, 3)
+ coor = mesh.coor()
+ dofs = mesh.dofs()
+ dofs[mesh.nodesLeftOpenEdge(), ...] = dofs[mesh.nodesRightOpenEdge(), ...]
+ iip = np.concatenate(
+ (dofs[mesh.nodesBottomEdge(), :].ravel(), dofs[mesh.nodesTopEdge(), 1])
+ )
+
+ elastic = np.array([0, 1, 2, 6, 7, 8], dtype=np.uint64)
+ plastic = np.array([3, 4, 5], dtype=np.uint64)
+ epsy = np.cumsum(1000 * np.ones((plastic.size, 5)), axis=1)
+
+ system = FrictionQPotFEM.Generic2d.System(
+ coor=coor,
+ conn=mesh.conn(),
+ dofs=dofs,
+ iip=iip,
+ elastic_elem=elastic,
+ elastic_K=FrictionQPotFEM.moduli_toquad(np.ones(elastic.size)),
+ elastic_G=FrictionQPotFEM.moduli_toquad(np.ones(elastic.size)),
+ plastic_elem=plastic,
+ plastic_K=FrictionQPotFEM.moduli_toquad(np.ones(plastic.size)),
+ plastic_G=FrictionQPotFEM.moduli_toquad(np.ones(plastic.size)),
+ plastic_epsy=FrictionQPotFEM.epsy_initelastic_toquad(epsy),
+ dt=1,
+ rho=1,
+ alpha=1,
+ eta=0,
+ )
+
+ sigxy = 0.1
+ system.applyShearStress(sigxy)
+
+ vector = GooseFEM.VectorPartitioned(mesh.conn(), dofs, iip)
+ matrix = GooseFEM.MatrixPartitioned(mesh.conn(), dofs, iip)
+ solver = GooseFEM.MatrixPartitionedSolver()
+ elem = GooseFEM.Element.Quad4.Quadrature(vector.AsElement(coor))
+ mat = GMat.Elastic2d(np.ones([vector.nelem, elem.nip]), np.ones([vector.nelem, elem.nip]))
+
+ Ke = elem.Int_gradN_dot_tensor4_dot_gradNT_dV(mat.C)
+ matrix.assemble(Ke)
+
+ system.u = np.zeros_like(system.u)
+ solver.solve(matrix, system.fext, system.u)
+
+ ue = vector.AsElement(system.u)
+ elem.symGradN_vector(ue, mat.Eps)
+ system.refresh()
+ mat.refresh()
+
+ # stress must be uniform
+ self.assertTrue(np.allclose(mat.Sig[..., 0, 0], 0))
+ self.assertTrue(np.allclose(mat.Sig[..., 1, 1], 0))
+ self.assertTrue(np.allclose(mat.Sig[..., 0, 1], mat.Sig[0, 0, 0, 1]))
+ self.assertTrue(np.allclose(mat.Sig[..., 1, 0], mat.Sig[0, 0, 1, 0]))
+
+ dV = system.quad.AsTensor(2, system.quad.dV)
+ sigbar = np.average(system.Sig(), weights=dV, axis=(0, 1))
+ self.assertAlmostEqual(sigbar[0, 1], sigxy)
+
  def test_eventDrivenSimpleShear(self):
  """
  Simple test of event driven simple shear in a homogeneous system: