Generic2d.h

Go to the documentation of this file.

 
#ifndef FRICTIONQPOTFEM_GENERIC2D_H
#define FRICTIONQPOTFEM_GENERIC2D_H
 
#include <algorithm>
#include <string>
#include <tuple>
 
#include "config.h"
 
#include <xtensor/xnorm.hpp>
#include <xtensor/xset_operation.hpp>
#include <xtensor/xshape.hpp>
#include <xtensor/xtensor.hpp>
 
#include <GooseFEM/GooseFEM.h>
#include <GooseFEM/Matrix.h>
#include <GooseFEM/MatrixPartitioned.h>
 
#include <GMatElastoPlasticQPot/Cartesian2d.h>
#include <GMatElastoPlasticQPot/version.h>
 
#include <GMatTensor/Cartesian2d.h>
#include <GMatTensor/version.h>
 
namespace FrictionQPotFEM {
 
namespace detail {
 
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
 
inline array_type::tensor<double, 3> epsy_initelastic_toquad(
    const array_type::tensor<double, 2>& arg,
    size_t nip = GooseFEM::Element::Quad4::Gauss::nip())
{
    array_type::tensor<double, 3> ret = xt::empty<double>({arg.shape(0), nip, arg.shape(1) + 1});
 
    for (size_t e = 0; e < arg.shape(0); ++e) {
        for (size_t q = 0; q < nip; ++q) {
            std::copy(&arg(e, 0), &arg(e, 0) + arg.shape(1), &ret(e, q, 1));
            ret(e, q, 0) = -arg(e, 0);
        }
    }
 
    return ret;
}
 
inline array_type::tensor<double, 2> moduli_toquad(
    const array_type::tensor<double, 1>& arg,
    size_t nip = GooseFEM::Element::Quad4::Gauss::nip())
{
    array_type::tensor<double, 2> ret = xt::empty<double>({arg.shape(0), nip});
 
    for (size_t e = 0; e < arg.shape(0); ++e) {
        for (size_t q = 0; q < nip; ++q) {
            ret(e, q) = arg(e);
        }
    }
 
    return ret;
}
 
template <class T>
inline typename T::value_type getuniform(const T& arg)
{
    if (xt::allclose(arg.flat(0), arg)) {
        return arg.flat(0);
    }
 
    throw std::runtime_error("Values not uniform");
}
 
namespace Generic2d {
 
namespace detail {
 
bool is_same(double a, double b)
{
    return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b &&
           std::nextafter(a, std::numeric_limits<double>::max()) >= b;
}
 
} // namespace detail
 
inline std::vector<std::string> version_dependencies()
{
    return GMatTensor::version_dependencies();
}
 
inline std::vector<std::string> version_compiler()
{
    return GMatTensor::version_compiler();
}
 
class System {
 
public:
    System() = default;
 
public:
    virtual ~System(){};
 
public:
    template <class C, class E, class L, class M, class Y>
    System(
        const C& coor,
        const E& conn,
        const E& dofs,
        const L& iip,
        const L& elastic_elem,
        const M& elastic_K,
        const M& elastic_G,
        const L& plastic_elem,
        const M& plastic_K,
        const M& plastic_G,
        const Y& plastic_epsy,
        double dt,
        double rho,
        double alpha,
        double eta)
    {
        this->initSystem(
            coor,
            conn,
            dofs,
            iip,
            elastic_elem,
            elastic_K,
            elastic_G,
            plastic_elem,
            plastic_K,
            plastic_G,
            plastic_epsy,
            dt,
            rho,
            alpha,
            eta);
    }
 
protected:
    template <class C, class E, class L, class M, class Y>
    void initSystem(
        const C& coor,
        const E& conn,
        const E& dofs,
        const L& iip,
        const L& elastic_elem,
        const M& elastic_K,
        const M& elastic_G,
        const L& plastic_elem,
        const M& plastic_K,
        const M& plastic_G,
        const Y& plastic_epsy,
        double dt,
        double rho,
        double alpha,
        double eta)
    {
        FRICTIONQPOTFEM_ASSERT(coor.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(coor.shape(1) == 2);
        FRICTIONQPOTFEM_ASSERT(conn.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(dofs.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(dofs.shape(1) == 2);
        FRICTIONQPOTFEM_ASSERT(iip.dimension() == 1);
        FRICTIONQPOTFEM_ASSERT(GooseFEM::is_unique(iip));
        FRICTIONQPOTFEM_ASSERT(elastic_elem.dimension() == 1);
        FRICTIONQPOTFEM_ASSERT(elastic_K.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(elastic_K.shape(0) == elastic_elem.size());
        FRICTIONQPOTFEM_ASSERT(elastic_G.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(elastic_G.shape(0) == elastic_elem.size());
        FRICTIONQPOTFEM_ASSERT(plastic_elem.dimension() == 1);
        FRICTIONQPOTFEM_ASSERT(plastic_K.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(plastic_K.shape(0) == plastic_elem.size());
        FRICTIONQPOTFEM_ASSERT(plastic_G.dimension() == 2);
        FRICTIONQPOTFEM_ASSERT(plastic_G.shape(0) == plastic_elem.size());
        FRICTIONQPOTFEM_ASSERT(plastic_epsy.dimension() == 3);
        FRICTIONQPOTFEM_ASSERT(plastic_epsy.shape(0) == plastic_elem.size());
 
        m_inc = 0;
        m_full_outdated = true;
 
        m_coor = coor;
        m_elem_elas = elastic_elem;
        m_elem_plas = plastic_elem;
 
        array_type::tensor<size_t, 2> conn_elas = xt::view(conn, xt::keep(m_elem_elas), xt::all());
        array_type::tensor<size_t, 2> conn_plas = xt::view(conn, xt::keep(m_elem_plas), xt::all());
 
        m_nnode = m_coor.shape(0);
        m_ndim = m_coor.shape(1);
        m_nelem = conn.shape(0);
        m_nne = conn.shape(1);
 
        m_nelem_elas = m_elem_elas.size();
        m_nelem_plas = m_elem_plas.size();
        m_N = m_nelem_plas;
 
#ifdef FRICTIONQPOTFEM_ENABLE_ASSERT
        // check that "elem_plastic" and "elem_plastic" together span all elements
        array_type::tensor<size_t, 1> elem = xt::concatenate(xt::xtuple(m_elem_elas, m_elem_plas));
        FRICTIONQPOTFEM_ASSERT(xt::sort(elem) == xt::arange<size_t>(m_nelem));
        // check that all "iip" or in "dofs"
        FRICTIONQPOTFEM_ASSERT(xt::all(xt::isin(iip, dofs)));
#endif
 
        m_vector = GooseFEM::VectorPartitioned(conn, dofs, iip);
        m_vector_elas = GooseFEM::VectorPartitioned(conn_elas, dofs, iip);
        m_vector_plas = GooseFEM::VectorPartitioned(conn_plas, dofs, iip);
 
        m_quad = GooseFEM::Element::Quad4::Quadrature(m_vector.AsElement(m_coor));
        m_quad_elas = GooseFEM::Element::Quad4::Quadrature(m_vector_elas.AsElement(m_coor));
        m_quad_plas = GooseFEM::Element::Quad4::Quadrature(m_vector_plas.AsElement(m_coor));
        m_nip = m_quad.nip();
 
        FRICTIONQPOTFEM_ASSERT(elastic_K.shape(1) == m_nip);
        FRICTIONQPOTFEM_ASSERT(elastic_G.shape(1) == m_nip);
        FRICTIONQPOTFEM_ASSERT(plastic_K.shape(1) == m_nip);
        FRICTIONQPOTFEM_ASSERT(plastic_G.shape(1) == m_nip);
        FRICTIONQPOTFEM_ASSERT(plastic_epsy.shape(1) == m_nip);
 
        m_u = m_vector.allocate_nodevec(0.0);
        m_v = m_vector.allocate_nodevec(0.0);
        m_a = m_vector.allocate_nodevec(0.0);
        m_v_n = m_vector.allocate_nodevec(0.0);
        m_a_n = m_vector.allocate_nodevec(0.0);
 
        m_ue = m_vector.allocate_elemvec(0.0);
        m_fe = m_vector.allocate_elemvec(0.0);
        m_ue_elas = m_vector_elas.allocate_elemvec(0.0);
        m_fe_elas = m_vector_elas.allocate_elemvec(0.0);
        m_ue_plas = m_vector_plas.allocate_elemvec(0.0);
        m_fe_plas = m_vector_plas.allocate_elemvec(0.0);
 
        m_fmaterial = m_vector.allocate_nodevec(0.0);
        m_felas = m_vector.allocate_nodevec(0.0);
        m_fplas = m_vector.allocate_nodevec(0.0);
        m_fdamp = m_vector.allocate_nodevec(0.0);
        m_fvisco = m_vector.allocate_nodevec(0.0);
        m_ftmp = m_vector.allocate_nodevec(0.0);
        m_fint = m_vector.allocate_nodevec(0.0);
        m_fext = m_vector.allocate_nodevec(0.0);
        m_fres = m_vector.allocate_nodevec(0.0);
 
        m_Eps = m_quad.allocate_qtensor<2>(0.0);
        m_Sig = m_quad.allocate_qtensor<2>(0.0);
        m_Epsdot_plas = m_quad_plas.allocate_qtensor<2>(0.0);
 
        m_M = GooseFEM::MatrixDiagonalPartitioned(conn, dofs, iip);
        m_D = GooseFEM::MatrixDiagonal(conn, dofs);
 
        m_K_elas = GooseFEM::Matrix(conn_elas, dofs);
 
        // allocated to strain-free state, matching #m_u
        m_elas = GMatElastoPlasticQPot::Cartesian2d::Elastic<2>(elastic_K, elastic_G);
        m_plas = GMatElastoPlasticQPot::Cartesian2d::Cusp<2>(plastic_K, plastic_G, plastic_epsy);
 
        m_K_elas.assemble(m_quad_elas.Int_gradN_dot_tensor4_dot_gradNT_dV(m_elas.C()));
 
        this->setDt(dt);
        this->setRho(rho);
        this->setAlpha(alpha);
        this->setEta(eta);
    }
public:
    virtual size_t N() const
    {
        return m_nelem_plas;
    }
 
public:
    virtual std::string type() const
    {
        return "FrictionQPotFEM.Generic2d.System";
    }
 
public:
    double rho() const
    {
        return m_rho;
    }
 
public:
    void setRho(double rho)
    {
        return this->setMassMatrix(xt::eval(rho * xt::ones<double>({m_nelem})));
    }
 
public:
    template <class T>
    void setMassMatrix(const T& val_elem)
    {
        FRICTIONQPOTFEM_ASSERT(val_elem.size() == m_nelem);
 
        if (xt::allclose(val_elem, val_elem(0))) {
            m_rho = val_elem(0);
        }
        else {
            m_rho = 0.0;
        }
 
        GooseFEM::Element::Quad4::Quadrature nodalQuad(
            m_vector.AsElement(m_coor),
            GooseFEM::Element::Quad4::Nodal::xi(),
            GooseFEM::Element::Quad4::Nodal::w());
 
        array_type::tensor<double, 2> val_quad = xt::empty<double>({m_nelem, nodalQuad.nip()});
        for (size_t q = 0; q < nodalQuad.nip(); ++q) {
            xt::view(val_quad, xt::all(), q) = val_elem;
        }
 
        m_M.assemble(nodalQuad.Int_N_scalar_NT_dV(val_quad));
    }
 
public:
    void setEta(double eta)
    {
        if (detail::is_same(eta, 0.0)) {
            m_set_visco = false;
            m_eta = 0.0;
        }
        else {
            m_set_visco = true;
            m_eta = eta;
        }
    }
 
public:
    double eta() const
    {
        return m_eta;
    }
 
public:
    void setAlpha(double alpha)
    {
        if (detail::is_same(alpha, 0.0)) {
            m_set_D = false;
            m_alpha = 0.0;
            return;
        }
 
        return this->setDampingMatrix(xt::eval(alpha * xt::ones<double>({m_nelem})));
    }
 
public:
    double alpha() const
    {
        return m_alpha;
    }
 
public:
    template <class T>
    void setDampingMatrix(const T& val_elem)
    {
        FRICTIONQPOTFEM_ASSERT(val_elem.size() == m_nelem);
        m_set_D = true;
 
        if (xt::allclose(val_elem, val_elem(0))) {
            m_alpha = val_elem(0);
            if (detail::is_same(m_alpha, 0.0)) {
                m_set_D = false;
            }
        }
        else {
            m_alpha = 0.0;
        }
 
        GooseFEM::Element::Quad4::Quadrature nodalQuad(
            m_vector.AsElement(m_coor),
            GooseFEM::Element::Quad4::Nodal::xi(),
            GooseFEM::Element::Quad4::Nodal::w());
 
        array_type::tensor<double, 2> val_quad = xt::empty<double>({m_nelem, nodalQuad.nip()});
        for (size_t q = 0; q < nodalQuad.nip(); ++q) {
            xt::view(val_quad, xt::all(), q) = val_elem;
        }
 
        m_D.assemble(nodalQuad.Int_N_scalar_NT_dV(val_quad));
    }
 
public:
    bool isHomogeneousElastic() const
    {
        auto k = this->K();
        auto g = this->G();
 
        return xt::allclose(k, k.data()[0]) && xt::allclose(g, g.data()[0]);
    }
 
public:
    double dt() const
    {
        return m_dt;
    }
 
    void setDt(double dt)
    {
        m_dt = dt;
    }
 
public:
    template <class T>
    void setU(const T& u)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(u, {m_nnode, m_ndim}));
        xt::noalias(m_u) = u;
        m_full_outdated = true;
        this->updated_u();
    }
 
public:
    template <class T>
    void setV(const T& v)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(v, {m_nnode, m_ndim}));
        xt::noalias(m_v) = v;
        this->updated_v();
    }
 
public:
    template <class T>
    void setA(const T& a)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(a, {m_nnode, m_ndim}));
        xt::noalias(m_a) = a;
    }
 
public:
    virtual void updated_u()
    {
        this->computeForceMaterial();
    }
 
    void updated_v()
    {
        if (m_set_D) {
            m_D.dot(m_v, m_fdamp);
        }
 
        // m_ue_plas, m_fe_plas, m_ftmp are temporaries that can be reused
        if (m_set_visco) {
            m_vector_plas.asElement(m_v, m_ue_plas);
            m_quad_plas.symGradN_vector(m_ue_plas, m_Epsdot_plas);
            m_quad_plas.int_gradN_dot_tensor2_dV(m_Epsdot_plas, m_fe_plas);
            if (!m_set_D) {
                m_vector_plas.assembleNode(m_fe_plas, m_fdamp);
                m_fdamp *= m_eta;
            }
            else {
                m_vector_plas.assembleNode(m_fe_plas, m_ftmp);
                m_ftmp *= m_eta;
                m_fdamp += m_ftmp;
            }
        }
    }
 
public:
    template <class T>
    void setFext(const T& fext)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(fext, {m_nnode, m_ndim}));
        xt::noalias(m_fext) = fext;
    }
 
public:
    const auto& elastic_elem() const
    {
        return m_elem_elas;
    }
 
public:
    const auto& plastic_elem() const
    {
        return m_elem_plas;
    }
 
public:
    const auto& conn() const
    {
        return m_vector.conn();
    }
 
public:
    const auto& coor() const
    {
        return m_coor;
    }
 
public:
    const auto& dofs() const
    {
        return m_vector.dofs();
    }
 
public:
    const auto& u() const
    {
        return m_u;
    }
 
public:
    const auto& v() const
    {
        return m_v;
    }
 
public:
    const auto& a() const
    {
        return m_a;
    }
 
public:
    auto& mass() const
    {
        return m_M;
    }
 
public:
    auto& damping() const
    {
        return m_D;
    }
 
public:
    const auto& fext()
    {
        this->computeInternalExternalResidualForce();
        return m_fext;
    }
 
    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:
    const auto& fint()
    {
        this->computeInternalExternalResidualForce();
        return m_fint;
    }
 
public:
    const auto& fmaterial() const
    {
        return m_fmaterial;
    }
 
public:
    const auto& fdamp() const
    {
        return m_fdamp;
    }
 
public:
    double residual() const
    {
        double r_fres = xt::norm_l2(m_fres)();
        double r_fext = xt::norm_l2(m_fext)();
        if (r_fext != 0.0) {
            return r_fres / r_fext;
        }
        return r_fres;
    }
 
public:
    void quench()
    {
        m_v.fill(0.0);
        m_a.fill(0.0);
        this->updated_v();
    }
 
public:
    double t() const
    {
        return m_inc * m_dt;
    }
 
    void setT(double arg)
    {
        this->setInc(static_cast<size_t>(std::round(arg / m_dt)));
        FRICTIONQPOTFEM_REQUIRE(xt::allclose(this->t(), arg));
    }
 
    size_t inc() const
    {
        return m_inc;
    }
 
    virtual void setInc(size_t arg)
    {
        m_inc = arg;
    }
 
public:
    const auto& dV() const
    {
        return m_quad.dV();
    }
 
public:
    const GooseFEM::VectorPartitioned& vector() const
    {
        return m_vector;
    }
 
public:
    const GooseFEM::Element::Quad4::Quadrature& quad() const
    {
        return m_quad;
    }
 
public:
    GMatElastoPlasticQPot::Cartesian2d::Elastic<2>& elastic()
    {
        return m_elas;
    }
 
public:
    GMatElastoPlasticQPot::Cartesian2d::Cusp<2>& plastic()
    {
        return m_plas;
    }
 
public:
    array_type::tensor<double, 2> K() const
    {
        array_type::tensor<double, 2> ret = xt::empty<double>({m_nelem, m_nip});
 
        const auto& ret_elas = m_elas.K();
        const auto& ret_plas = m_plas.K();
        size_t n = xt::strides(ret_elas, 0);
        FRICTIONQPOTFEM_ASSERT(n == m_nip);
 
        for (size_t e = 0; e < m_nelem_elas; ++e) {
            std::copy(
                &ret_elas(e, xt::missing),
                &ret_elas(e, xt::missing) + n,
                &ret(m_elem_elas(e), xt::missing));
        }
 
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            std::copy(
                &ret_plas(e, xt::missing),
                &ret_plas(e, xt::missing) + n,
                &ret(m_elem_plas(e), xt::missing));
        }
 
        return ret;
    }
 
public:
    array_type::tensor<double, 2> G() const
    {
        array_type::tensor<double, 2> ret = xt::empty<double>({m_nelem, m_nip});
 
        const auto& ret_elas = m_elas.G();
        const auto& ret_plas = m_plas.G();
        size_t n = xt::strides(ret_elas, 0);
        FRICTIONQPOTFEM_ASSERT(n == m_nip);
 
        for (size_t e = 0; e < m_nelem_elas; ++e) {
            std::copy(
                &ret_elas(e, xt::missing),
                &ret_elas(e, xt::missing) + n,
                &ret(m_elem_elas(e), xt::missing));
        }
 
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            std::copy(
                &ret_plas(e, xt::missing),
                &ret_plas(e, xt::missing) + n,
                &ret(m_elem_plas(e), xt::missing));
        }
 
        return ret;
    }
 
public:
    const array_type::tensor<double, 4>& Sig()
    {
        this->computeEpsSig();
        return m_Sig;
    }
 
public:
    const array_type::tensor<double, 4>& Eps()
    {
        this->computeEpsSig();
        return m_Eps;
    }
 
public:
    array_type::tensor<double, 4> Epsdot() const
    {
        return m_quad.SymGradN_vector(m_vector.AsElement(m_v));
    }
 
public:
    array_type::tensor<double, 4> Epsddot() const
    {
        return m_quad.SymGradN_vector(m_vector.AsElement(m_a));
    }
 
public:
    GooseFEM::MatrixPartitioned stiffness() const
    {
        auto ret = m_quad.allocate_qtensor<4>(0.0);
        const auto& ret_plas = m_plas.C();
        const auto& ret_elas = m_elas.C();
        size_t n = xt::strides(ret_elas, 0);
        FRICTIONQPOTFEM_ASSERT(n == m_nip * 16);
 
        for (size_t e = 0; e < m_nelem_elas; ++e) {
            std::copy(
                &ret_elas(e, xt::missing),
                &ret_elas(e, xt::missing) + n,
                &ret(m_elem_elas(e), xt::missing));
        }
 
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            std::copy(
                &ret_plas(e, xt::missing),
                &ret_plas(e, xt::missing) + n,
                &ret(m_elem_plas(e), xt::missing));
        }
 
        GooseFEM::MatrixPartitioned K(m_vector.conn(), m_vector.dofs(), m_vector.iip());
        K.assemble(m_quad.Int_gradN_dot_tensor4_dot_gradNT_dV(ret));
        return K;
    }
 
public:
    const array_type::tensor<double, 4>& plastic_Epsdot()
    {
        // m_ue_plas, m_fe_plas are temporaries that can be reused
        if (!m_set_visco) {
            m_vector_plas.asElement(m_v, m_ue_plas);
            m_quad_plas.symGradN_vector(m_ue_plas, m_Epsdot_plas);
        }
 
        return m_Epsdot_plas;
    }
 
protected:
    void computeEpsSig()
    {
        if (!m_full_outdated) {
            return;
        }
 
        auto& Eps_elas = m_elas.Eps();
        auto& Sig_elas = m_elas.Sig();
        const auto& Eps_plas = m_plas.Eps();
        const auto& Sig_plas = m_plas.Sig();
 
        m_vector_elas.asElement(m_u, m_ue_elas);
        m_quad_elas.symGradN_vector(m_ue_elas, Eps_elas);
        m_elas.refresh();
 
        size_t n = xt::strides(Eps_elas, 0);
        FRICTIONQPOTFEM_ASSERT(n == m_nip * 4);
 
        for (size_t e = 0; e < m_nelem_elas; ++e) {
            std::copy(
                &Eps_elas(e, xt::missing),
                &Eps_elas(e, xt::missing) + n,
                &m_Eps(m_elem_elas(e), xt::missing));
            std::copy(
                &Sig_elas(e, xt::missing),
                &Sig_elas(e, xt::missing) + n,
                &m_Sig(m_elem_elas(e), xt::missing));
        }
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            std::copy(
                &Eps_plas(e, xt::missing),
                &Eps_plas(e, xt::missing) + n,
                &m_Eps(m_elem_plas(e), xt::missing));
            std::copy(
                &Sig_plas(e, xt::missing),
                &Sig_plas(e, xt::missing) + n,
                &m_Sig(m_elem_plas(e), xt::missing));
        }
 
        m_full_outdated = false;
    }
 
    void computeForceMaterial()
    {
        m_full_outdated = true;
 
        m_vector_plas.asElement(m_u, m_ue_plas);
        m_quad_plas.symGradN_vector(m_ue_plas, m_plas.Eps());
        m_plas.refresh();
        m_quad_plas.int_gradN_dot_tensor2_dV(m_plas.Sig(), m_fe_plas);
        m_vector_plas.assembleNode(m_fe_plas, m_fplas);
 
        m_K_elas.dot(m_u, m_felas);
 
        xt::noalias(m_fmaterial) = m_felas + m_fplas;
    }
 
protected:
    virtual void computeInternalExternalResidualForce()
    {
        xt::noalias(m_fint) = m_fmaterial + m_fdamp;
        m_vector.copy_p(m_fint, m_fext);
        xt::noalias(m_fres) = m_fext - m_fint;
    }
 
public:
    void refresh()
    {
        this->updated_u();
        this->updated_v();
        this->computeInternalExternalResidualForce();
    }
 
public:
    template <class T>
    double eventDriven_setDeltaU(const T& delta_u, bool autoscale = true)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(delta_u, m_u.shape()));
        m_pert_u = delta_u;
 
        m_vector_plas.asElement(delta_u, m_ue_plas);
        auto Eps = m_quad_plas.SymGradN_vector(m_ue_plas);
        m_pert_Epsd_plastic = GMatTensor::Cartesian2d::Deviatoric(Eps);
 
        if (!autoscale) {
            return 1.0;
        }
 
        auto deps = xt::amax(GMatElastoPlasticQPot::Cartesian2d::Epsd(m_pert_Epsd_plastic))();
        auto d = xt::amin(xt::diff(m_plas.epsy(), 1))();
        double c = 0.1 * d / deps;
 
        m_pert_u *= c;
        m_pert_Epsd_plastic *= c;
 
        return c;
    }
 
    const auto& eventDriven_deltaU() const
    {
        return m_pert_u;
    }
 
    virtual double eventDrivenStep(
        double deps,
        bool kick,
        int direction = 1,
        bool yield_element = false,
        bool iterative = false)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(m_pert_u, m_u.shape()));
        FRICTIONQPOTFEM_ASSERT(direction == 1 || direction == -1);
 
        auto eps = GMatElastoPlasticQPot::Cartesian2d::Epsd(m_plas.Eps());
        auto Epsd_plastic = GMatTensor::Cartesian2d::Deviatoric(m_plas.Eps());
        auto idx = m_plas.i();
        const auto& epsy_l = m_plas.epsy_left();
        const auto& epsy_r = m_plas.epsy_right();
 
        FRICTIONQPOTFEM_WIP(iterative || direction > 0 || !xt::any(xt::equal(idx, 0)));
 
        array_type::tensor<double, 2> target;
 
        if (direction > 0 && kick) { // direction > 0 && kick
            target = epsy_r + 0.5 * deps;
        }
        else if (direction > 0) { // direction > 0 && !kick
            target = epsy_r - 0.5 * deps;
        }
        else if (kick) { // direction < 0 && kick
            target = epsy_l - 0.5 * deps;
        }
        else { // direction < 0 && !kick
            target = epsy_l + 0.5 * deps;
        }
 
        if (!kick) {
            auto d = target - eps;
            if (direction > 0 && xt::any(d < 0.0)) {
                return 0.0;
            }
            if (direction < 0 && xt::any(d > 0.0)) {
                return 0.0;
            }
        }
 
        auto scale = xt::empty_like(target);
 
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            for (size_t q = 0; q < m_nip; ++q) {
 
                double e_t = Epsd_plastic(e, q, 0, 0);
                double g_t = Epsd_plastic(e, q, 0, 1);
                double e_d = m_pert_Epsd_plastic(e, q, 0, 0);
                double g_d = m_pert_Epsd_plastic(e, q, 0, 1);
                double epsd_target = target(e, q);
 
                double a = e_d * e_d + g_d * g_d;
                double b = 2.0 * (e_t * e_d + g_t * g_d);
                double c = e_t * e_t + g_t * g_t - epsd_target * epsd_target;
                double D = std::sqrt(b * b - 4.0 * a * c);
 
                FRICTIONQPOTFEM_REQUIRE(b >= 0.0 || iterative);
                scale(e, q) = (-b + D) / (2.0 * a);
            }
        }
 
        if (kick) {
            size_t nyield = m_plas.epsy().shape(m_plas.epsy().dimension() - 1);
            size_t bound = direction > 0 ? nyield - 2 : 0;
            if (xt::all(xt::equal(idx, bound))) {
                return 0.0;
            }
            scale = xt::where(xt::equal(idx, bound), std::numeric_limits<double>::max(), scale);
        }
 
        double ret;
        size_t e;
        size_t q;
 
        if (!iterative) {
 
            auto index = xt::unravel_index(xt::argmin(xt::abs(scale))(), scale.shape());
            e = index[0];
            q = index[1];
 
            if (kick && yield_element) {
                q = xt::argmax(xt::view(xt::abs(scale), e, xt::all()))();
            }
 
            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));
                }
            }
        }
 
        else {
 
            double dir = static_cast<double>(direction);
            auto steps = xt::sort(xt::ravel(xt::eval(xt::abs(scale))));
            auto u_n = m_u;
            xt::xtensor<long, 2> jdx = xt::cast<long>(m_plas.i());
            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>(m_plas.i());
                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;
                    }
                }
            }
 
            // 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 step = 0; step < 1100; ++step) {
 
                ret = 0.5 * (right + left);
                this->setU(u_n + dir * ret * m_pert_u);
                auto jdx_new = xt::cast<long>(m_plas.i());
                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(step < 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>(m_plas.i());
                auto S = xt::abs(jdx_new - jdx);
 
                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>(m_plas.i());
                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));
        }
 
        this->setU(m_u + ret * m_pert_u);
 
        const auto& idx_new = m_plas.i();
        FRICTIONQPOTFEM_REQUIRE(kick || xt::all(xt::equal(idx, idx_new)));
        FRICTIONQPOTFEM_REQUIRE(!kick || xt::any(xt::not_equal(idx, idx_new)));
 
        if (!iterative) {
            auto eps_new = GMatElastoPlasticQPot::Cartesian2d::Epsd(m_plas.Eps())(e, q);
            auto eps_target = target(e, q);
            FRICTIONQPOTFEM_REQUIRE(xt::allclose(eps_new, eps_target));
        }
 
        return ret;
    }
 
    void timeStep()
    {
        // history
 
        m_inc++;
        xt::noalias(m_v_n) = m_v;
        xt::noalias(m_a_n) = m_a;
 
        // new displacement
 
        xt::noalias(m_u) = m_u + m_dt * m_v + 0.5 * std::pow(m_dt, 2.0) * m_a;
        this->updated_u();
 
        // estimate new velocity, update corresponding force
 
        xt::noalias(m_v) = m_v_n + m_dt * m_a_n;
        this->updated_v();
 
        // compute residual force & solve
 
        this->computeInternalExternalResidualForce();
        m_M.solve(m_fres, m_a);
 
        // re-estimate new velocity, update corresponding force
 
        xt::noalias(m_v) = m_v_n + 0.5 * m_dt * (m_a_n + m_a);
        this->updated_v();
 
        // compute residual force & solve
 
        this->computeInternalExternalResidualForce();
        m_M.solve(m_fres, m_a);
 
        // new velocity, update corresponding force
 
        xt::noalias(m_v) = m_v_n + 0.5 * m_dt * (m_a_n + m_a);
        this->updated_v();
 
        // compute residual force & solve
 
        this->computeInternalExternalResidualForce();
        m_M.solve(m_fres, m_a);
    }
 
    long timeSteps(size_t n, size_t nmargin = 0)
    {
        FRICTIONQPOTFEM_REQUIRE(n + 1 < std::numeric_limits<long>::max());
 
        size_t nyield = m_plas.epsy().shape(2);
        size_t nmax = nyield - nmargin;
        long step;
 
        FRICTIONQPOTFEM_ASSERT(nmargin < nyield);
        FRICTIONQPOTFEM_REQUIRE(xt::all(m_plas.i() <= nmax));
 
        for (step = 1; step < static_cast<long>(n + 1); ++step) {
 
            this->timeStep();
 
            if (nmargin > 0) {
                if (xt::any(m_plas.i() > nmax)) {
                    return -step;
                }
            }
        }
 
        return step;
    }
 
    size_t timeStepsUntilEvent(
        size_t nmargin = 0,
        double tol = 1e-5,
        size_t niter_tol = 20,
        size_t max_iter = 1e7)
    {
        FRICTIONQPOTFEM_ASSERT(tol < 1.0);
        FRICTIONQPOTFEM_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
 
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
        size_t nyield = m_plas.epsy().shape(2);
        size_t nmax = nyield - nmargin;
        auto i_n = m_plas.i();
        size_t step;
 
        FRICTIONQPOTFEM_ASSERT(nmargin < nyield);
        FRICTIONQPOTFEM_REQUIRE(xt::all(m_plas.i() <= nmax));
 
        for (step = 1; step < max_iter + 1; ++step) {
 
            this->timeStep();
 
            if (xt::any(xt::not_equal(m_plas.i(), i_n))) {
                return step;
            }
 
            residuals.roll_insert(this->residual());
 
            if ((residuals.descending() && residuals.all_less(tol)) || residuals.all_less(tol2)) {
                this->quench();
                return 0;
            }
        }
 
        return step;
    }
 
    std::tuple<
        array_type::tensor<double, 1>,
        array_type::tensor<double, 2>,
        array_type::tensor<double, 2>>
    minimise_highfrequency(
        const array_type::tensor<size_t, 1>& nodes,
        const array_type::tensor<size_t, 1>& top,
        size_t t_step = 1,
        size_t nmargin = 0,
        double tol = 1e-5,
        size_t niter_tol = 20,
        size_t max_iter = 1e7)
    {
        FRICTIONQPOTFEM_ASSERT(tol < 1.0);
        FRICTIONQPOTFEM_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
 
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
 
        size_t nyield = m_plas.epsy().shape(2);
        size_t nmax = nyield - nmargin;
 
        FRICTIONQPOTFEM_ASSERT(nmargin < nyield);
        FRICTIONQPOTFEM_REQUIRE(xt::all(m_plas.i() <= nmax));
        FRICTIONQPOTFEM_REQUIRE(nmargin == 0 || xt::all(m_plas.i() >= nmargin));
 
        std::vector<double> fext;
        std::vector<std::vector<double>> ux(nodes.size());
        std::vector<std::vector<double>> uy(nodes.size());
 
        for (size_t step = 1; step < max_iter + 1; ++step) {
 
            this->timeStep();
            residuals.roll_insert(this->residual());
 
            if (step % t_step == 0) {
 
                double f = 0.0;
                for (auto& n : top) {
                    f += m_fext(n, 0);
                }
                fext.push_back(f / static_cast<double>(top.size() - 1)); // account for periodicity
 
                for (size_t i = 0; i < nodes.size(); ++i) {
                    size_t n = nodes(i);
                    ux[i].push_back(m_u(n, 0));
                    uy[i].push_back(m_u(n, 1));
                }
            }
 
            if (nmargin > 0) {
                if (xt::any(m_plas.i() > nmax) || xt::any(m_plas.i() < nmargin)) {
                    throw std::runtime_error("Out of bounds");
                }
            }
 
            if ((residuals.descending() && residuals.all_less(tol)) || residuals.all_less(tol2)) {
                this->quench();
 
                size_t n = fext.size();
                array_type::tensor<double, 1> ret_fext = xt::empty<double>({n});
                std::copy(fext.begin(), fext.end(), ret_fext.begin());
                fext.resize(0);
 
                array_type::tensor<double, 2> ret_ux = xt::empty<double>({nodes.size(), n});
                for (size_t i = 0; i < nodes.size(); ++i) {
                    std::copy(ux[i].begin(), ux[i].end(), ret_ux.begin() + i * n);
                }
                ux.resize(0);
 
                array_type::tensor<double, 2> ret_uy = xt::empty<double>({nodes.size(), n});
                for (size_t i = 0; i < nodes.size(); ++i) {
                    std::copy(uy[i].begin(), uy[i].end(), ret_uy.begin() + i * n);
                }
                uy.resize(0);
 
                return std::make_tuple(ret_fext, ret_ux, ret_uy);
            }
        }
 
        throw std::runtime_error("No convergence found");
    }
 
    template <class T>
    long flowSteps(size_t n, const T& v, size_t nmargin = 0)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(v, m_u.shape()));
        FRICTIONQPOTFEM_REQUIRE(n + 1 < std::numeric_limits<long>::max());
 
        size_t nyield = m_plas.epsy().shape(2);
        size_t nmax = nyield - nmargin;
        long step;
 
        FRICTIONQPOTFEM_ASSERT(nmargin < nyield);
        FRICTIONQPOTFEM_REQUIRE(xt::all(m_plas.i() <= nmax));
 
        for (step = 1; step < static_cast<long>(n + 1); ++step) {
 
            m_u += v * m_dt;
            this->timeStep();
 
            if (nmargin > 0) {
                if (xt::any(m_plas.i() > nmax)) {
                    return -step;
                }
            }
        }
 
        return step;
    }
 
    long minimise(
        size_t nmargin = 0,
        double tol = 1e-5,
        size_t niter_tol = 20,
        size_t max_iter = 1e7,
        bool time_activity = false,
        bool max_iter_is_error = true)
    {
        FRICTIONQPOTFEM_ASSERT(tol < 1.0);
        FRICTIONQPOTFEM_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
 
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
 
        size_t nyield = m_plas.epsy().shape(2);
        size_t nmax = nyield - nmargin;
        array_type::tensor<long, 2> i_n = xt::cast<long>(m_plas.i());
        long s = 0;
        long s_n = 0;
        bool init = true;
        long step;
 
        FRICTIONQPOTFEM_ASSERT(nmargin < nyield);
        FRICTIONQPOTFEM_REQUIRE(xt::all(m_plas.i() <= nmax));
 
        for (step = 1; step < static_cast<long>(max_iter + 1); ++step) {
 
            this->timeStep();
            residuals.roll_insert(this->residual());
 
            if (nmargin > 0) {
                if (xt::any(m_plas.i() > nmax)) {
                    return -step;
                }
            }
 
            if (time_activity) {
                array_type::tensor<long, 2> i = xt::cast<long>(m_plas.i());
                s = xt::sum(xt::abs(i - i_n))();
                if (s != s_n) {
                    if (init) {
                        init = false;
                        m_qs_inc_first = m_inc;
                    }
                    m_qs_inc_last = m_inc;
                }
                s_n = s;
            }
 
            if ((residuals.descending() && residuals.all_less(tol)) || residuals.all_less(tol2)) {
                this->quench();
                return 0;
            }
        }
 
        if (max_iter_is_error) {
            throw std::runtime_error("No convergence found");
        }
 
        return step;
    }
 
    size_t quasistaticActivityFirst() const
    {
        return m_qs_inc_first;
    }
 
    size_t quasistaticActivityLast() const
    {
        return m_qs_inc_last;
    }
 
    template <class T>
    size_t minimise_truncate(
        const T& idx_n,
        size_t A_truncate = 0,
        size_t S_truncate = 0,
        double tol = 1e-5,
        size_t niter_tol = 20,
        size_t max_iter = 1e7)
    {
        FRICTIONQPOTFEM_REQUIRE(xt::has_shape(idx_n, std::array<size_t, 1>{m_N}));
        FRICTIONQPOTFEM_REQUIRE(S_truncate == 0);
        FRICTIONQPOTFEM_REQUIRE(A_truncate > 0);
        FRICTIONQPOTFEM_ASSERT(tol < 1.0);
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
 
        for (size_t step = 1; step < max_iter + 1; ++step) {
 
            this->timeStep();
            residuals.roll_insert(this->residual());
 
            if ((residuals.descending() && residuals.all_less(tol)) || residuals.all_less(tol2)) {
                this->quench();
                return step;
            }
 
            array_type::tensor<size_t, 1> idx = xt::view(m_plas.i(), xt::all(), 0);
 
            if (static_cast<size_t>(xt::sum(xt::not_equal(idx_n, idx))()) >= A_truncate) {
                return 0;
            }
        }
 
        std::cout << "residuals = " << xt::adapt(residuals.data()) << std::endl;
        bool converged = false;
        FRICTIONQPOTFEM_REQUIRE(converged == true);
        return 0; // irrelevant, the code never goes here
    }
 
    size_t minimise_truncate(
        size_t A_truncate = 0,
        size_t S_truncate = 0,
        double tol = 1e-5,
        size_t niter_tol = 20,
        size_t max_iter = 1e7)
    {
        array_type::tensor<size_t, 1> idx_n = xt::view(m_plas.i(), xt::all(), 0);
        return this->minimise_truncate(idx_n, A_truncate, S_truncate, tol, niter_tol, max_iter);
    }
 
    array_type::tensor<double, 2> affineSimpleShear(double delta_gamma) const
    {
        array_type::tensor<double, 2> ret = xt::zeros_like(m_u);
 
        for (size_t n = 0; n < m_nnode; ++n) {
            ret(n, 0) += 2.0 * delta_gamma * (m_coor(n, 1) - m_coor(0, 1));
        }
 
        return ret;
    }
 
    array_type::tensor<double, 2> affineSimpleShearCentered(double delta_gamma) const
    {
        array_type::tensor<double, 2> ret = xt::zeros_like(m_u);
        size_t ll = m_vector.conn()(m_elem_plas(0), 0);
        size_t ul = m_vector.conn()(m_elem_plas(0), 3);
        double y0 = (m_coor(ul, 1) + m_coor(ll, 1)) / 2.0;
 
        for (size_t n = 0; n < m_nnode; ++n) {
            ret(n, 0) += 2.0 * delta_gamma * (m_coor(n, 1) - y0);
        }
 
        return ret;
    }
 
protected:
    array_type::tensor<double, 2> m_coor; 
    size_t m_N; 
    size_t m_nelem; 
    size_t m_nelem_elas; 
    size_t m_nelem_plas; 
    size_t m_nne; 
    size_t m_nnode; 
    size_t m_ndim; 
    size_t m_nip; 
    array_type::tensor<size_t, 1> m_elem_elas; 
    array_type::tensor<size_t, 1> m_elem_plas; 
    GooseFEM::Element::Quad4::Quadrature m_quad; 
    GooseFEM::Element::Quad4::Quadrature m_quad_elas; 
    GooseFEM::Element::Quad4::Quadrature m_quad_plas; 
    GooseFEM::VectorPartitioned m_vector; 
    GooseFEM::VectorPartitioned m_vector_elas; 
    GooseFEM::VectorPartitioned m_vector_plas; 
    GooseFEM::MatrixDiagonalPartitioned m_M; 
    GooseFEM::MatrixDiagonal m_D; 
    GMatElastoPlasticQPot::Cartesian2d::Elastic<2> m_elas; 
    GMatElastoPlasticQPot::Cartesian2d::Cusp<2> m_plas; 
 
    array_type::tensor<double, 2> m_u;
 
    array_type::tensor<double, 2> m_v;
 
    array_type::tensor<double, 2> m_a; 
    array_type::tensor<double, 2> m_v_n; 
    array_type::tensor<double, 2> m_a_n; 
    array_type::tensor<double, 3> m_ue; 
    array_type::tensor<double, 3> m_fe; 
    array_type::tensor<double, 3> m_ue_elas; 
    array_type::tensor<double, 3> m_fe_elas; 
    array_type::tensor<double, 3> m_ue_plas; 
    array_type::tensor<double, 3> m_fe_plas; 
    array_type::tensor<double, 2> m_fmaterial; 
    array_type::tensor<double, 2> m_felas; 
    array_type::tensor<double, 2> m_fplas; 
    array_type::tensor<double, 2> m_fdamp; 
    array_type::tensor<double, 2> m_fvisco; 
    array_type::tensor<double, 2> m_ftmp; 
    array_type::tensor<double, 2> m_fint; 
    array_type::tensor<double, 2> m_fext; 
    array_type::tensor<double, 2> m_fres; 
    array_type::tensor<double, 4> m_Eps; 
    array_type::tensor<double, 4> m_Sig; 
    array_type::tensor<double, 4> m_Epsdot_plas; 
    GooseFEM::Matrix m_K_elas; 
    size_t m_qs_inc_first = 0; 
    size_t m_qs_inc_last = 0; 
    size_t m_inc; 
    double m_dt; 
    double m_eta; 
    double m_rho; 
    double m_alpha; 
    bool m_set_D; 
    bool m_set_visco; 
    bool m_full_outdated; 
    array_type::tensor<double, 2> m_pert_u; 
    array_type::tensor<double, 4> m_pert_Epsd_plastic; 
};
 
} // namespace Generic2d
} // namespace FrictionQPotFEM
 
#endif