Working around xtensor bug

include/FrictionQPotFEM/Generic2d.h

@@ -29,6 +29,40 @@
 
 namespace FrictionQPotFEM {
 
+namespace detail {
+
+/**
+Drop-in replacement for xt::average(a, w, {0, 1});
+There appears to be a library bug.
+*/
+array_type::tensor<double, 2> myaverage(
+    const array_type::tensor<double, 4>& a,
+    const array_type::tensor<double, 4>& w,
+    const std::array<size_t, 2>& axis)
+{
+    FRICTIONQPOTFEM_ASSERT(axis[0] == 0);
+    FRICTIONQPOTFEM_ASSERT(axis[1] == 1);
+    FRICTIONQPOTFEM_ASSERT(xt::has_shape(a, w.shape()));
+
+    array_type::tensor<double, 2> ret = xt::zeros<double>({a.shape(2), a.shape(3)});
+    array_type::tensor<double, 2> norm = xt::zeros<double>({a.shape(2), a.shape(3)});
+
+    for (size_t i = 0; i < a.shape(0); ++i) {
+        for (size_t j = 0; j < a.shape(1); ++j) {
+            for (size_t k = 0; k < a.shape(2); ++k) {
+                for (size_t l = 0; l < a.shape(3); ++l) {
+                    ret(k, l) += a(i, j, k, l) * w(i, j, k, l);
+                    norm(k, l) += w(i, j, k, l);
+                }
+            }
+        }
+    }
+
+    return ret / norm;
+}
+
+} // namespace detail
+
 /**
 Convert array of yield strains stored per element [nelem, n]:
 

include/FrictionQPotFEM/UniformSingleLayer2d.h

@@ -166,12 +166,11 @@ class System : public Generic2d::System {
     {
         FRICTIONQPOTFEM_ASSERT(this->isHomogeneousElastic());
 
-        auto u_new = this->u();
         auto dV = m_quad.AsTensor<2>(this->dV());
         double G = m_plas.G().flat(0);
 
-        array_type::tensor<double, 2> Epsbar = xt::average(this->Eps(), dV, {0, 1});
-        array_type::tensor<double, 2> Sigbar = xt::average(this->Sig(), dV, {0, 1});
+        array_type::tensor<double, 2> Epsbar = detail::myaverage(this->Eps(), dV, {0, 1});
+        array_type::tensor<double, 2> Sigbar = detail::myaverage(this->Sig(), dV, {0, 1});
         array_type::tensor<double, 2> Epsd = GMatTensor::Cartesian2d::Deviatoric(Epsbar);
         double epsxx = Epsd(0, 0);
         double epsxy = Epsd(0, 1);
@@ -194,16 +193,15 @@ class System : public Generic2d::System {
         }
 
         for (size_t n = 0; n < m_nnode; ++n) {
-            u_new(n, 0) += direction * dgamma * (m_coor(n, 1) - m_coor(0, 1));
+            m_u(n, 0) += direction * dgamma * (m_coor(n, 1) - m_coor(0, 1));
         }
-        this->setU(u_new);
+        this->updated_u();
 
         // sanity check
         // ------------
 
-        Sigbar = xt::average(this->Sig(), dV, {0, 1});
+        Sigbar = detail::myaverage(this->Sig(), dV, {0, 1});
         sig = GMatElastoPlasticQPot::Cartesian2d::Sigd(Sigbar)();
-
         FRICTIONQPOTFEM_REQUIRE(std::abs(target_stress - sig) / sig < 1e-4);
 
         return direction * dgamma;

tests/test_UniformSingleLayer2d.py

@@ -61,6 +61,56 @@ def test_eventDrivenSimpleShear(self):
                 self.assertTrue(np.allclose(GMat.Epsd(system.plastic.Eps), eps_expect))
                 self.assertTrue(system.residual() < 1e-5)
 
+    def test_addSimpleShearToFixedStress(self):
+
+        mesh = GooseFEM.Mesh.Quad4.Regular(3, 3)
+        coor = mesh.coor()
+        dofs = mesh.dofs()
+        dofs[mesh.nodesLeftOpenEdge(), ...] = dofs[mesh.nodesRightOpenEdge(), ...]
+
+        elastic = np.array([0, 1, 2, 6, 7, 8], dtype=np.uint64)
+        plastic = np.array([3, 4, 5], dtype=np.uint64)
+        epsy = 1e-3 * np.cumsum(np.ones((plastic.size, 5)), axis=1)
+
+        system = FrictionQPotFEM.UniformSingleLayer2d.System(
+            coor=coor,
+            conn=mesh.conn(),
+            dofs=dofs,
+            iip=dofs[np.concatenate((mesh.nodesBottomEdge(), mesh.nodesTopEdge())), :].ravel(),
+            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=0.01,
+            rho=1,
+            alpha=0.1,
+            eta=0,
+        )
+
+        system.initEventDrivenSimpleShear()
+        dV = system.quad.AsTensor(2, system.quad.dV)
+
+        for step in range(4):
+
+            u_n = np.copy(system.u)
+            sig_n = GMat.Sigd(np.average(system.Sig(), weights=dV, axis=(0, 1)))
+
+            system.eventDrivenStep(1e-4, step % 2 == 0, direction=1)
+            system.minimise()
+
+            sig = GMat.Sigd(np.average(system.Sig(), weights=dV, axis=(0, 1)))
+
+        target = sig_n + 0.5 * (sig - sig_n)
+        system.u = u_n
+        system.addSimpleShearToFixedStress(target)
+
+        self.assertTrue(
+            np.isclose(target, GMat.Sigd(np.average(system.Sig(), weights=dV, axis=(0, 1))))
+        )
+
 
 if __name__ == "__main__":