UniformMultiLayerIndividualDrive2d.h

Go to the documentation of this file.

 
#ifndef FRICTIONQPOTFEM_UNIFORMMULTILAYERINDIVIDUALDRIVE2D_H
#define FRICTIONQPOTFEM_UNIFORMMULTILAYERINDIVIDUALDRIVE2D_H
 
#include "Generic2d.h"
#include "config.h"
 
namespace FrictionQPotFEM {
 
namespace UniformMultiLayerIndividualDrive2d {
 
inline std::vector<std::string> version_dependencies()
{
    return Generic2d::version_dependencies();
}
 
inline std::vector<std::string> version_compiler()
{
    return Generic2d::version_compiler();
}
 
class System : public Generic2d::System {
 
public:
    System() = default;
 
    virtual ~System(){};
 
    System(
        const array_type::tensor<double, 2>& coor,
        const array_type::tensor<size_t, 2>& conn,
        const array_type::tensor<size_t, 2>& dofs,
        const array_type::tensor<size_t, 1>& iip,
        const std::vector<array_type::tensor<size_t, 1>>& elem,
        const std::vector<array_type::tensor<size_t, 1>>& node,
        const array_type::tensor<bool, 1>& layer_is_plastic,
        const array_type::tensor<double, 2>& elastic_K,
        const array_type::tensor<double, 2>& elastic_G,
        const array_type::tensor<double, 2>& plastic_K,
        const array_type::tensor<double, 2>& plastic_G,
        const array_type::tensor<double, 3>& plastic_epsy,
        double dt,
        double rho,
        double alpha,
        double eta,
        const array_type::tensor<bool, 2>& drive_is_active,
        double k_drive)
    {
        this->init(
            coor,
            conn,
            dofs,
            iip,
            elem,
            node,
            layer_is_plastic,
            elastic_K,
            elastic_G,
            plastic_K,
            plastic_G,
            plastic_epsy,
            dt,
            rho,
            alpha,
            eta,
            drive_is_active,
            k_drive);
    }
 
protected:
    void init(
        const array_type::tensor<double, 2>& coor,
        const array_type::tensor<size_t, 2>& conn,
        const array_type::tensor<size_t, 2>& dofs,
        const array_type::tensor<size_t, 1>& iip,
        const std::vector<array_type::tensor<size_t, 1>>& elem,
        const std::vector<array_type::tensor<size_t, 1>>& node,
        const array_type::tensor<bool, 1>& layer_is_plastic,
        const array_type::tensor<double, 2>& elastic_K,
        const array_type::tensor<double, 2>& elastic_G,
        const array_type::tensor<double, 2>& plastic_K,
        const array_type::tensor<double, 2>& plastic_G,
        const array_type::tensor<double, 3>& plastic_epsy,
        double dt,
        double rho,
        double alpha,
        double eta,
        const array_type::tensor<bool, 2>& drive_is_active,
        double k_drive)
    {
        FRICTIONQPOTFEM_ASSERT(layer_is_plastic.size() == elem.size());
        FRICTIONQPOTFEM_ASSERT(layer_is_plastic.size() == node.size());
 
        using size_type = typename decltype(m_layerdrive_active)::size_type;
 
        m_n_layer = node.size();
        m_layer_node = node;
        m_layer_elem = elem;
        m_layer_is_plastic = layer_is_plastic;
 
        std::array<size_type, 1> shape1 = {static_cast<size_type>(m_n_layer)};
        std::array<size_type, 2> shape2 = {static_cast<size_type>(m_n_layer), 2};
        m_layerdrive_active.resize(shape2);
        m_layerdrive_targetubar.resize(shape2);
        m_layer_ubar.resize(shape2);
        m_layer_dV1.resize(shape2);
        m_slice_index.resize(shape1);
 
        m_layerdrive_targetubar.fill(0.0);
        m_layerdrive_active.fill(false);
 
        size_t n_elem_plas = 0;
        size_t n_elem_elas = 0;
        size_t n_layer_plas = 0;
        size_t n_layer_elas = 0;
 
        for (size_t i = 0; i < elem.size(); ++i) {
            if (m_layer_is_plastic(i)) {
                n_elem_plas += elem[i].size();
                n_layer_plas++;
            }
            else {
                n_elem_elas += elem[i].size();
                n_layer_elas++;
            }
        }
 
        array_type::tensor<size_t, 1> plas = xt::empty<size_t>({n_elem_plas});
        array_type::tensor<size_t, 1> elas = xt::empty<size_t>({n_elem_elas});
 
        std::array<size_type, 1> size_plas = {static_cast<size_type>(n_layer_plas) + 1};
        std::array<size_type, 1> size_elas = {static_cast<size_type>(n_layer_elas) + 1};
        m_slice_plas.resize(size_plas);
        m_slice_elas.resize(size_elas);
        m_slice_plas(0) = 0;
        m_slice_elas(0) = 0;
        m_N = 0;
 
        n_elem_plas = 0;
        n_elem_elas = 0;
        n_layer_plas = 0;
        n_layer_elas = 0;
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            if (m_layer_is_plastic(i)) {
                size_t l = m_slice_plas(n_layer_plas);
                size_t u = n_elem_plas + elem[i].size();
 
                m_slice_index(i) = n_layer_plas;
                m_slice_plas(n_layer_plas + 1) = u;
 
                xt::view(plas, xt::range(l, u)) = elem[i];
 
                n_elem_plas += elem[i].size();
                n_layer_plas++;
 
                FRICTIONQPOTFEM_REQUIRE(m_N == elem[i].size() || m_N == 0);
                m_N = elem[i].size();
            }
            else {
                size_t l = m_slice_elas(n_layer_elas);
                size_t u = n_elem_elas + elem[i].size();
 
                m_slice_index(i) = n_layer_elas;
                m_slice_elas(n_layer_elas + 1) = u;
 
                xt::view(elas, xt::range(l, u)) = elem[i];
 
                n_elem_elas += elem[i].size();
                n_layer_elas++;
            }
        }
 
        this->initSystem(
            coor,
            conn,
            dofs,
            iip,
            elas,
            elastic_K,
            elastic_G,
            plas,
            plastic_K,
            plastic_G,
            plastic_epsy,
            dt,
            rho,
            alpha,
            eta);
 
        m_fdrive = m_vector.allocate_nodevec(0.0);
        m_ud = m_vector.allocate_dofval(0.0);
        m_uq = m_quad.allocate_qtensor<1>(0.0);
        m_dV = m_quad.dV();
 
        size_t nip = m_quad.nip();
        m_layer_dV1.fill(0.0);
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (auto& e : m_layer_elem[i]) {
                for (size_t q = 0; q < nip; ++q) {
                    for (size_t d = 0; d < 2; ++d) {
                        m_layer_dV1(i, d) += m_dV(e, q);
                    }
                }
            }
        }
 
// sanity check nodes per layer
#ifdef FRICTIONQPOTFEM_ENABLE_ASSERT
        for (size_t i = 0; i < m_n_layer; ++i) {
            auto e = this->layerElements(i);
            auto n = xt::unique(xt::view(m_vector.conn(), xt::keep(e)));
            FRICTIONQPOTFEM_ASSERT(xt::all(xt::equal(xt::sort(n), xt::sort(node[i]))));
        }
#endif
 
// sanity check elements per layer + slice of elas/plas element sets
#ifdef FRICTIONQPOTFEM_ENABLE_ASSERT
        for (size_t i = 0; i < m_n_layer; ++i) {
            array_type::tensor<size_t, 1> e;
            size_t j = m_slice_index(i);
            if (m_layer_is_plastic(i)) {
                e = xt::view(m_elem_plas, xt::range(m_slice_plas(j), m_slice_plas(j + 1)));
            }
            else {
                e = xt::view(m_elem_elas, xt::range(m_slice_elas(j), m_slice_elas(j + 1)));
            }
            FRICTIONQPOTFEM_ASSERT(xt::all(xt::equal(xt::sort(e), xt::sort(elem[i]))));
        }
#endif
 
        this->layerSetDriveStiffness(k_drive);
        this->layerSetTargetActive(drive_is_active);
    }
public:
    size_t N() const override
    {
        return m_N;
    }
 
    std::string type() const override
    {
        return "FrictionQPotFEM.UniformMultiLayerIndividualDrive2d.System";
    }
 
    size_t nlayer() const
    {
        return m_n_layer;
    }
 
    array_type::tensor<size_t, 1> layerElements(size_t i) const
    {
        FRICTIONQPOTFEM_ASSERT(i < m_n_layer);
        return m_layer_elem[i];
    }
 
    array_type::tensor<size_t, 1> layerNodes(size_t i) const
    {
        FRICTIONQPOTFEM_ASSERT(i < m_n_layer);
        return m_layer_node[i];
    }
 
    const array_type::tensor<bool, 1>& layerIsPlastic() const
    {
        return m_layer_is_plastic;
    }
 
    void layerSetDriveStiffness(double k, bool symmetric = true)
    {
        m_drive_spring_symmetric = symmetric;
        m_drive_k = k;
        this->computeLayerUbarActive();
        this->computeForceFromTargetUbar();
    }
 
    template <class S, class T>
    double initEventDriven(const S& ubar, const T& active)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(ubar, m_layerdrive_targetubar.shape()));
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(active, m_layerdrive_active.shape()));
 
        // backup system
 
        auto u0 = m_u;
        auto v0 = m_v;
        auto a0 = m_a;
        auto active0 = m_layerdrive_active;
        auto ubar0 = m_layerdrive_targetubar;
        auto inc0 = m_inc;
 
        array_type::tensor<double, 3> epsy0 = m_plas.epsy();
        array_type::tensor<double, 3> rigid = xt::empty<double>({m_nelem_plas, m_nip, size_t(2)});
 
        double infty = std::numeric_limits<double>::max();
 
        for (size_t e = 0; e < m_nelem_plas; ++e) {
            for (size_t q = 0; q < m_nip; ++q) {
                rigid(e, q, 0) = -infty;
                rigid(e, q, 1) = infty;
            }
        }
 
        m_plas.set_epsy(rigid);
 
        // perturbation
 
        m_u.fill(0.0);
        m_v.fill(0.0);
        m_a.fill(0.0);
        this->updated_u();
        FRICTIONQPOTFEM_ASSERT(xt::all(xt::equal(m_plas.i(), 0)));
        this->layerSetTargetActive(active);
        this->layerSetTargetUbar(ubar);
        this->minimise(0);
        FRICTIONQPOTFEM_ASSERT(xt::all(xt::equal(m_plas.i(), 0)));
        auto up = m_u;
        m_u.fill(0.0);
 
        // restore yield strains
 
        m_plas.set_epsy(epsy0);
 
        // store result
 
        auto c = this->eventDriven_setDeltaU(up);
        m_pert_layerdrive_active = active;
        m_pert_layerdrive_targetubar = c * ubar;
 
        // restore system
 
        this->setU(u0);
        this->setV(v0);
        this->setA(a0);
        this->layerSetTargetActive(active0);
        this->layerSetTargetUbar(ubar0);
        this->setInc(inc0);
 
        return c;
    }
 
    template <class S, class T, class U>
    double initEventDriven(const S& ubar, const T& active, const U& delta_u)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(ubar, m_layerdrive_targetubar.shape()));
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(active, m_layerdrive_active.shape()));
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(delta_u, m_u.shape()));
        double c = this->eventDriven_setDeltaU(delta_u);
        m_pert_layerdrive_active = active;
        m_pert_layerdrive_targetubar = c * ubar;
        return c;
    }
 
    const array_type::tensor<double, 2>& eventDriven_deltaUbar() const
    {
        return m_pert_layerdrive_targetubar;
    }
 
    const array_type::tensor<bool, 2>& eventDriven_targetActive() const
    {
        return m_pert_layerdrive_active;
    }
 
    double eventDrivenStep(
        double deps,
        bool kick,
        int direction = +1,
        bool yield_element = false,
        bool fallback = false) override
    {
        double c =
            Generic2d::System::eventDrivenStep(deps, kick, direction, yield_element, fallback);
        this->layerSetTargetUbar(m_layerdrive_targetubar + c * m_pert_layerdrive_targetubar);
        return c;
    }
 
    template <class T>
    void layerSetTargetActive(const T& active)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(active, m_layerdrive_active.shape()));
        m_layerdrive_active = active;
        this->computeLayerUbarActive();
        this->computeForceFromTargetUbar();
    }
 
    array_type::tensor<double, 2> layerUbar()
    {
        // Recompute needed because computeLayerUbarActive() only computes the average
        // on layers with an active spring.
        // This function, in contrast, returns the average on all layers.
 
        m_layer_ubar.fill(0.0);
        size_t nip = m_quad.nip();
 
        m_vector.asElement(m_u, m_ue);
        m_quad.interpQuad_vector(m_ue, m_uq);
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (auto& e : m_layer_elem[i]) {
                for (size_t d = 0; d < 2; ++d) {
                    for (size_t q = 0; q < nip; ++q) {
                        m_layer_ubar(i, d) += m_uq(e, q, d) * m_dV(e, q);
                    }
                }
            }
        }
 
        m_layer_ubar /= m_layer_dV1;
 
        return m_layer_ubar;
    }
 
    const array_type::tensor<double, 2>& layerTargetUbar() const
    {
        return m_layerdrive_targetubar;
    }
 
    const array_type::tensor<bool, 2>& layerTargetActive() const
    {
        return m_layerdrive_active;
    }
 
    template <class T>
    void layerSetTargetUbar(const T& ubar)
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(ubar, m_layerdrive_targetubar.shape()));
        m_layerdrive_targetubar = ubar;
        this->computeForceFromTargetUbar(); // average displacement and other forces do not change
    }
 
    template <class S, class T>
    void layerSetUbar(const S& ubar, const T& prescribe)
    {
        auto current = this->layerUbar();
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (size_t d = 0; d < 2; ++d) {
                if (prescribe(i, d)) {
                    double du = ubar(i, d) - current(i, d);
                    for (auto& n : m_layer_node[i]) {
                        m_u(n, d) += du;
                    }
                }
            }
        }
 
        this->computeLayerUbarActive();
        this->computeForceFromTargetUbar();
        this->computeForceMaterial();
    }
 
    template <class T>
    array_type::tensor<double, 2>
    layerTargetUbar_affineSimpleShear(double delta_gamma, const T& height) const
    {
        FRICTIONQPOTFEM_ASSERT(xt::has_shape(height, {m_n_layer}));
 
        auto ret = xt::zeros_like(m_layerdrive_targetubar);
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            ret(i, 0) += 2.0 * delta_gamma * height(i);
        }
 
        return ret;
    }
 
    const array_type::tensor<double, 2>& fdrive() const
    {
        return m_fdrive;
    }
 
    array_type::tensor<double, 2> layerFdrive() const
    {
        array_type::tensor<double, 2> ret = xt::zeros<double>({m_n_layer, size_t(2)});
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (size_t d = 0; d < 2; ++d) {
                if (m_layerdrive_active(i, d)) {
                    ret(i, d) = m_drive_k * (m_layerdrive_targetubar(i, d) - m_layer_ubar(i, d));
                }
            }
        }
 
        return ret;
    }
 
protected:
    void computeLayerUbarActive()
    {
        m_layer_ubar.fill(0.0);
        size_t nip = m_quad.nip();
 
        m_vector.asElement(m_u, m_ue);
        m_quad.interpQuad_vector(m_ue, m_uq);
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (size_t d = 0; d < 2; ++d) {
                if (m_layerdrive_active(i, d)) {
                    for (auto& e : m_layer_elem[i]) {
                        for (size_t q = 0; q < nip; ++q) {
                            m_layer_ubar(i, d) += m_uq(e, q, d) * m_dV(e, q);
                        }
                    }
                }
            }
        }
 
        m_layer_ubar /= m_layer_dV1;
    }
 
    void computeForceFromTargetUbar()
    {
        m_uq.fill(0.0); // pre-allocated value that an be freely used
        size_t nip = m_quad.nip();
 
        for (size_t i = 0; i < m_n_layer; ++i) {
            for (size_t d = 0; d < 2; ++d) {
                if (m_layerdrive_active(i, d)) {
                    double f = m_drive_k * (m_layer_ubar(i, d) - m_layerdrive_targetubar(i, d));
                    if (m_drive_spring_symmetric || f < 0) { // buckle under compression
                        for (auto& e : m_layer_elem[i]) {
                            for (size_t q = 0; q < nip; ++q) {
                                m_uq(e, q, d) = f;
                            }
                        }
                    }
                }
            }
        }
 
        m_quad.int_N_vector_dV(m_uq, m_ue);
        m_vector.assembleDofs(m_ue, m_ud);
        m_vector.asNode(m_ud, m_fdrive);
    }
 
    void updated_u() override
    {
        this->computeLayerUbarActive();
        this->computeForceFromTargetUbar();
        this->computeForceMaterial();
    }
 
    void computeInternalExternalResidualForce() override
    {
        xt::noalias(m_fint) = m_fdrive + m_fmaterial + m_fdamp;
        m_vector.copy_p(m_fint, m_fext);
        xt::noalias(m_fres) = m_fext - m_fint;
    }
 
protected:
    size_t m_N; 
    size_t m_n_layer; 
    std::vector<array_type::tensor<size_t, 1>> m_layer_node; 
    std::vector<array_type::tensor<size_t, 1>> m_layer_elem; 
    array_type::tensor<bool, 1> m_layer_is_plastic; 
    array_type::tensor<size_t, 1>
        m_slice_index; 
    array_type::tensor<size_t, 1> m_slice_plas; 
    array_type::tensor<size_t, 1> m_slice_elas; 
 
    bool m_drive_spring_symmetric; 
    double m_drive_k; 
    array_type::tensor<bool, 2> m_layerdrive_active; 
    array_type::tensor<double, 2>
        m_layerdrive_targetubar; 
    array_type::tensor<double, 2> m_layer_ubar; 
    array_type::tensor<double, 2> m_fdrive; 
    array_type::tensor<double, 2>
        m_layer_dV1; 
    array_type::tensor<double, 2> m_dV; 
    array_type::tensor<double, 3> m_uq; 
    array_type::tensor<double, 1> m_ud; 
 
    array_type::tensor<bool, 2> m_pert_layerdrive_active; 
    array_type::tensor<double, 2>
        m_pert_layerdrive_targetubar; 
};
 
} // namespace UniformMultiLayerIndividualDrive2d
} // namespace FrictionQPotFEM
 
#endif