detail.h

Go to the documentation of this file.

 
#ifndef FRICTIONQPOTSPRINGBLOCK_DETAIL_H
#define FRICTIONQPOTSPRINGBLOCK_DETAIL_H
 
#include <string>
 
#include <xtensor/xadapt.hpp>
#include <xtensor/xnorm.hpp>
#include <xtensor/xtensor.hpp>
 
#include "config.h"
 
#include <GooseFEM/Iterate.h>
#include <GooseFEM/version.h>
#include <prrng.h>
 
#include <GMatTensor/version.h>
 
namespace FrictionQPotSpringBlock {
 
namespace detail {
 
inline prrng::distribution string_to_distribution(const std::string& str)
{
    if (str == "random") {
        return prrng::distribution::random;
    }
 
    if (str == "delta") {
        return prrng::distribution::delta;
    }
 
    if (str == "exponential") {
        return prrng::distribution::exponential;
    }
 
    if (str == "power") {
        return prrng::distribution::power;
    }
 
    if (str == "gamma") {
        return prrng::distribution::gamma;
    }
 
    if (str == "pareto") {
        return prrng::distribution::pareto;
    }
 
    if (str == "weibull") {
        return prrng::distribution::weibull;
    }
 
    if (str == "normal") {
        return prrng::distribution::normal;
    }
 
    throw std::runtime_error("Unknown distribution: " + str);
}
 
template <class U, class Y, class I>
inline bool check_disorder(const U& u, const Y& yield, const I& i)
{
    if (i.dimension() != 1 || yield.dimension() != 2) {
        return false;
    }
 
    if (u.size() != i.size() || u.size() != yield.shape(0)) {
        return false;
    }
 
    using index_type = typename I::value_type;
    index_type n = static_cast<ptrdiff_t>(yield.shape(1));
 
    for (size_t p = 0; p < yield.shape(0); ++p) {
        if (i(p) <= 0 || i(p) >= n - 2) {
            return false;
        }
 
        if (u.flat(p) < yield(p, i(p)) || u.flat(p) > yield(p, i(p) + 1)) {
            return false;
        }
 
        if (!std::is_sorted(&yield(p, 0), &yield(p, 0) + yield.shape(1))) {
            return false;
        }
    }
 
    return true;
}
 
template <class Generator>
class Cuspy {
protected:
    using stype = typename Generator::size_type; 
    stype m_N; 
    double m_mu; 
    Generator* m_chunk; 
 
public:
    Cuspy() = default;
 
    Cuspy(double mu, Generator* chunk) : m_mu(mu), m_chunk(chunk)
    {
        m_N = m_chunk->generators().size();
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG((stype)u.size() == m_N && (stype)f.size() == m_N);
 
        m_chunk->align(u);
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const double, 2> yield = xt::adapt(
            m_chunk->data().data(),
            m_chunk->data().size(),
            xt::no_ownership(),
            std::array<stype, 2>{m_N, m_chunk->chunk_size()}
        );
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const ptrdiff_t, 1> i = xt::adapt(
            m_chunk->chunk_index_at_align().data(),
            m_chunk->chunk_index_at_align().size(),
            xt::no_ownership(),
            std::array<stype, 1>{m_N}
        );
 
        FRICTIONQPOTSPRINGBLOCK_DEBUG(check_disorder(u, yield, i));
 
        for (stype p = 0; p < m_N; ++p) {
            const auto* l = &yield(p, i(p));
            f.flat(p) = 0.5 * (*(l) + *(l + 1)) - u.flat(p);
        }
 
        f *= m_mu;
    }
 
    template <class T>
    double maxUniformDisplacement(const T& u, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(direction == 1 || direction == -1);
        m_chunk->align(u);
 
        if (direction > 0) {
            return xt::amin(m_chunk->template right_of_align<T>() - u)();
        }
 
        return xt::amin(u - m_chunk->template left_of_align<T>())();
    }
 
    template <class T>
    void trigger(T& u, size_t p, double eps, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG((stype)p < m_N);
        m_chunk->align(u);
 
        if (direction > 0) {
            u.flat(p) = m_chunk->template right_of_align<T>().flat(p) + 0.5 * eps;
        }
        else {
            u.flat(p) = m_chunk->template left_of_align<T>().flat(p) - 0.5 * eps;
        }
    }
};
 
template <class Generator>
class SemiSmooth {
protected:
    using stype = typename Generator::size_type; 
    stype m_N; 
    double m_mu; 
    double m_kappa; 
    Generator* m_chunk; 
 
public:
    SemiSmooth() = default;
 
    SemiSmooth(double mu, double kappa, Generator* chunk) : m_mu(mu), m_kappa(kappa), m_chunk(chunk)
    {
        m_N = m_chunk->generators().size();
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG((stype)u.size() == m_N && (stype)f.size() == m_N);
 
        m_chunk->align(u);
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const double, 2> yield = xt::adapt(
            m_chunk->data().data(),
            m_chunk->data().size(),
            xt::no_ownership(),
            std::array<stype, 2>{m_N, m_chunk->chunk_size()}
        );
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const ptrdiff_t, 1> i = xt::adapt(
            m_chunk->chunk_index_at_align().data(),
            m_chunk->chunk_index_at_align().size(),
            xt::no_ownership(),
            std::array<stype, 1>{m_N}
        );
 
        FRICTIONQPOTSPRINGBLOCK_DEBUG(check_disorder(u, yield, i));
 
        for (stype p = 0; p < m_N; ++p) {
            auto* y = &yield(p, i(p));
            double xi = 0.5 * (*(y) + *(y + 1));
            double u_r = (m_mu * xi + m_kappa * *(y + 1)) / (m_mu + m_kappa);
            double u_l = (m_mu * xi + m_kappa * *(y)) / (m_mu + m_kappa);
            double up = u.flat(p);
            if (up < u_l) {
                f.flat(p) = m_kappa * (up - *(y));
            }
            else if (up <= u_r) {
                f.flat(p) = m_mu * (0.5 * (*(y) + *(y + 1)) - up);
            }
            else {
                f.flat(p) = m_kappa * (up - *(y + 1));
            }
        }
    }
 
    template <class T>
    double maxUniformDisplacement(const T& u, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(direction == 1 || direction == -1);
        m_chunk->align(u);
 
        bool positive = direction > 0;
        std::vector<double> du;
        du.reserve(m_N);
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const double, 2> yield = xt::adapt(
            m_chunk->data().data(),
            m_chunk->data().size(),
            xt::no_ownership(),
            std::array<stype, 2>{m_N, m_chunk->chunk_size()}
        );
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const ptrdiff_t, 1> i = xt::adapt(
            m_chunk->chunk_index_at_align().data(),
            m_chunk->chunk_index_at_align().size(),
            xt::no_ownership(),
            std::array<stype, 1>{m_N}
        );
 
        FRICTIONQPOTSPRINGBLOCK_DEBUG(check_disorder(u, yield, i));
 
        for (stype p = 0; p < m_N; ++p) {
            auto* y = &yield(p, i(p));
            double xi = 0.5 * (*(y) + *(y + 1));
            double u_r = (m_mu * xi + m_kappa * *(y + 1)) / (m_mu + m_kappa);
            double u_l = (m_mu * xi + m_kappa * *(y)) / (m_mu + m_kappa);
            double up = u.flat(p);
 
            if (up < u_l) {
                return 0.0;
            }
            else if (up <= u_r) {
                if (positive) {
                    du.push_back(u_r - up);
                }
                else {
                    du.push_back(up - u_l);
                }
            }
            else {
                return 0.0;
            }
        }
 
        return *std::min_element(du.begin(), du.end());
    }
 
    template <class T>
    void trigger(T& u, size_t p, double eps, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(p < m_N);
 
        if (direction > 0) {
            u.flat(p) = m_chunk->template right_of_align<T>().flat(p) + 0.5 * eps;
        }
        else {
            u.flat(p) = m_chunk->template left_of_align<T>().flat(p) - 0.5 * eps;
        }
    }
};
 
template <class Generator>
class Smooth {
protected:
    using stype = typename Generator::size_type; 
    stype m_N; 
    double m_mu; 
    Generator* m_chunk; 
 
public:
    Smooth() = default;
 
    Smooth(double mu, Generator* chunk) : m_mu(mu), m_chunk(chunk)
    {
        m_N = m_chunk->generators().size();
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG((stype)u.size() == m_N && (stype)f.size() == m_N);
 
        m_chunk->align(u);
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const double, 2> yield = xt::adapt(
            m_chunk->data().data(),
            m_chunk->data().size(),
            xt::no_ownership(),
            std::array<stype, 2>{m_N, m_chunk->chunk_size()}
        );
 
        // this does not allocate data, but only creates a view
        xt::xtensor_pointer<const ptrdiff_t, 1> i = xt::adapt(
            m_chunk->chunk_index_at_align().data(),
            m_chunk->chunk_index_at_align().size(),
            xt::no_ownership(),
            std::array<stype, 1>{m_N}
        );
 
        for (stype p = 0; p < m_N; ++p) {
            auto* y = &yield(p, i(p));
            double up = u.flat(p);
            double umin = 0.5 * (*(y) + *(y + 1));
            double dy = 0.5 * (*(y + 1) - *(y));
            f.flat(p) = -m_mu * dy / M_PI * std::sin(M_PI * (up - umin) / dy);
        }
    }
 
    template <class T>
    double maxUniformDisplacement(const T& u, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(direction == 1 || direction == -1);
        (void)(u);
        (void)(direction);
        throw std::runtime_error("Operation not possible.");
        return 0.0;
    }
 
    template <class T>
    void trigger(T& u, size_t p, double eps, int direction) const
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(p < m_N);
 
        if (direction > 0) {
            u.flat(p) = m_chunk->template right_of_align<T>().flat(p) + 0.5 * eps;
        }
        else {
            u.flat(p) = m_chunk->template left_of_align<T>().flat(p) - 0.5 * eps;
        }
    }
};
 
class Laplace1d {
protected:
    size_type m_N; 
    double m_k; 
 
public:
    Laplace1d() = default;
 
    Laplace1d(double k, size_type N) : m_N(N), m_k(k)
    {
    }
 
    double k() const
    {
        return m_k;
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(u.dimension() == 1 && f.dimension() == 1);
        FRICTIONQPOTSPRINGBLOCK_DEBUG((size_type)u.size() == m_N && (size_type)f.size() == m_N);
 
        for (size_type p = 1; p < m_N - 1; ++p) {
            f(p) = u(p - 1) - 2 * u(p) + u(p + 1);
        }
        f.front() = u.back() - 2 * u.front() + u(1);
        f.back() = u(m_N - 2) - 2 * u.back() + u.front();
 
        f *= m_k;
    }
};
 
class Laplace2d {
protected:
    size_type m_rows; 
    size_type m_cols; 
    double m_k; 
 
public:
    Laplace2d() = default;
 
    Laplace2d(double k, size_type rows, size_type cols) : m_rows(rows), m_cols(cols), m_k(k)
    {
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(u.dimension() == 2 && f.dimension() == 2);
        FRICTIONQPOTSPRINGBLOCK_DEBUG(xt::has_shape(u, f.shape()));
        FRICTIONQPOTSPRINGBLOCK_DEBUG(
            (size_type)u.shape(0) == m_rows && (size_type)u.shape(1) == m_cols
        );
 
        std::array<size_type, 2> edge_rows = {0, m_rows - 1};
        std::array<size_type, 2> edge_cols = {0, m_cols - 1};
 
        // first and last row, all columns
        for (size_type r = 0; r < 2; ++r) {
            for (size_type j = 0; j < m_cols; ++j) {
                size_type i = edge_rows[r];
                f(i, j) = u.periodic(i - 1, j) + u.periodic(i + 1, j) + u.periodic(i, j - 1) +
                          u.periodic(i, j + 1) - 4 * u(i, j);
            }
        }
 
        // first and last column, 'interior' rows
        for (size_type i = 1; i < m_rows - 1; ++i) {
            for (size_type c = 0; c < 2; ++c) {
                size_type j = edge_cols[c];
                f(i, j) = u.periodic(i - 1, j) + u.periodic(i + 1, j) + u.periodic(i, j - 1) +
                          u.periodic(i, j + 1) - 4 * u(i, j);
            }
        }
 
        // interior
        for (size_type i = 1; i < m_rows - 1; ++i) {
            for (size_type j = 1; j < m_cols - 1; ++j) {
                f(i, j) = u(i - 1, j) + u(i + 1, j) + u(i, j - 1) + u(i, j + 1) - 4 * u(i, j);
            }
        }
 
        f *= m_k;
    }
};
 
class QuarticGradient1d {
protected:
    size_type m_N; 
    double m_k2; 
    double m_k4; 
    double m_k4_4; 
 
public:
    QuarticGradient1d() = default;
 
    QuarticGradient1d(double k2, double k4, size_type N)
        : m_N(N), m_k2(k2), m_k4(k4), m_k4_4(0.25 * k4)
    {
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        for (size_type p = 1; p < m_N - 1; ++p) {
            double du = u(p + 1) - u(p - 1);
            f(p) = (u(p - 1) - 2 * u(p) + u(p + 1)) * (m_k2 + m_k4_4 * du * du);
        }
 
        double duf = u(1) - u.back();
        f.front() = (u.back() - 2 * u.front() + u(1)) * (m_k2 + m_k4_4 * duf * duf);
 
        double dub = u.front() - u(m_N - 2);
        f.back() = (u(m_N - 2) - 2 * u.back() + u.front()) * (m_k2 + m_k4_4 * dub * dub);
    }
};
 
class QuarticGradient2d {
protected:
    size_type m_rows; 
    size_type m_cols; 
    double m_k2; 
    double m_k4; 
 
public:
    QuarticGradient2d() = default;
 
    QuarticGradient2d(double k2, double k4, size_type rows, size_type cols)
        : m_rows(rows), m_cols(cols), m_k2(k2), m_k4(k4)
    {
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(u.dimension() == 2 && f.dimension() == 2);
        FRICTIONQPOTSPRINGBLOCK_DEBUG(xt::has_shape(u, f.shape()));
        FRICTIONQPOTSPRINGBLOCK_DEBUG(
            (size_type)u.shape(0) == m_rows && (size_type)u.shape(1) == m_cols
        );
 
        double mk4_3 = m_k4 / 3.0;
        double mk4_23 = 2.0 * mk4_3;
        std::array<size_type, 2> edge_rows = {0, m_rows - 1};
        std::array<size_type, 2> edge_cols = {0, m_cols - 1};
 
        // first and last row, all columns
        for (size_type r = 0; r < 2; ++r) {
            for (size_type j = 0; j < m_cols; ++j) {
                size_type i = edge_rows[r];
                double l = u.periodic(i + 1, j) + u.periodic(i - 1, j) + u.periodic(i, j + 1) +
                           u.periodic(i, j - 1) - 4 * u(i, j);
                double dudx = 0.5 * (u.periodic(i + 1, j) - u.periodic(i - 1, j));
                double dudy = 0.5 * (u.periodic(i, j + 1) - u.periodic(i, j - 1));
                double d2udxdy = 0.25 * (u.periodic(i + 1, j + 1) - u.periodic(i + 1, j - 1) -
                                         u.periodic(i - 1, j + 1) + u.periodic(i - 1, j - 1));
                double d2udx2 = u.periodic(i + 1, j) - 2 * u(i, j) + u.periodic(i - 1, j);
                double d2udy2 = u.periodic(i, j + 1) - 2 * u(i, j) + u.periodic(i, j - 1);
 
                f(i, j) =
                    l * (m_k2 + mk4_3) + mk4_23 * (dudx * dudx * d2udx2 + dudy * dudy * d2udy2 +
                                                   2.0 * dudx * dudy * d2udxdy);
            }
        }
 
        // first and last column, 'interior' rows
        for (size_type i = 1; i < m_rows - 1; ++i) {
            for (size_type c = 0; c < 2; ++c) {
                size_type j = edge_cols[c];
                double l = u.periodic(i + 1, j) + u.periodic(i - 1, j) + u.periodic(i, j + 1) +
                           u.periodic(i, j - 1) - 4 * u(i, j);
                double dudx = 0.5 * (u.periodic(i + 1, j) - u.periodic(i - 1, j));
                double dudy = 0.5 * (u.periodic(i, j + 1) - u.periodic(i, j - 1));
                double d2udxdy = 0.25 * (u.periodic(i + 1, j + 1) - u.periodic(i + 1, j - 1) -
                                         u.periodic(i - 1, j + 1) + u.periodic(i - 1, j - 1));
                double d2udx2 = u.periodic(i + 1, j) - 2 * u(i, j) + u.periodic(i - 1, j);
                double d2udy2 = u.periodic(i, j + 1) - 2 * u(i, j) + u.periodic(i, j - 1);
 
                f(i, j) =
                    l * (m_k2 + mk4_3) + mk4_23 * (dudx * dudx * d2udx2 + dudy * dudy * d2udy2 +
                                                   2.0 * dudx * dudy * d2udxdy);
            }
        }
 
        // interior
        for (size_type i = 1; i < m_rows - 1; ++i) {
            for (size_type j = 1; j < m_cols - 1; ++j) {
                double l = u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1) - 4 * u(i, j);
                double dudx = 0.5 * (u(i + 1, j) - u(i - 1, j));
                double dudy = 0.5 * (u(i, j + 1) - u(i, j - 1));
                double d2udxdy =
                    0.25 * (u(i + 1, j + 1) - u(i + 1, j - 1) - u(i - 1, j + 1) + u(i - 1, j - 1));
                double d2udx2 = u(i + 1, j) - 2 * u(i, j) + u(i - 1, j);
                double d2udy2 = u(i, j + 1) - 2 * u(i, j) + u(i, j - 1);
 
                f(i, j) =
                    l * (m_k2 + mk4_3) + mk4_23 * (dudx * dudx * d2udx2 + dudy * dudy * d2udy2 +
                                                   2.0 * dudx * dudy * d2udxdy);
            }
        }
    }
};
 
class Quartic1d {
protected:
    size_type m_N; 
    double m_a1; 
    double m_a2; 
 
public:
    Quartic1d() = default;
 
    Quartic1d(double a1, double a2, size_type N) : m_N(N), m_a1(a1), m_a2(a2)
    {
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        FRICTIONQPOTSPRINGBLOCK_DEBUG(u.dimension() == 1 && f.dimension() == 1);
        FRICTIONQPOTSPRINGBLOCK_DEBUG((size_type)u.size() == m_N && (size_type)f.size() == m_N);
 
        for (size_type p = 1; p < m_N - 1; ++p) {
            double dup = u(p + 1) - u(p);
            double dun = u(p - 1) - u(p);
            f(p) = m_a1 * (u(p - 1) - 2 * u(p) + u(p + 1)) +
                   m_a2 * (dup * dup * dup + dun * dun * dun);
        }
 
        {
            double dup = u(1) - u.front();
            double dun = u.back() - u.front();
            f.front() = m_a1 * (u.back() - 2 * u.front() + u(1)) +
                        m_a2 * (dup * dup * dup + dun * dun * dun);
        }
 
        {
            double dup = u.front() - u.back();
            double dun = u(m_N - 2) - u.back();
            f.back() = m_a1 * (u(m_N - 2) - 2 * u.back() + u.front()) +
                       m_a2 * (dup * dup * dup + dun * dun * dun);
        }
    }
};
 
class LongRange1d {
protected:
    size_type m_N; 
    double m_k; 
    double m_alpha; 
    ptrdiff_t m_n; 
    ptrdiff_t m_m; 
    array_type::tensor<double, 1> m_prefactor; 
 
public:
    LongRange1d() = default;
 
    LongRange1d(double k, double alpha, size_type N) : m_N(N), m_k(k), m_alpha(alpha)
    {
        m_n = static_cast<ptrdiff_t>(m_N);
        m_m = (m_n - m_n % 2) / 2;
 
        m_prefactor = xt::empty<double>({m_N});
 
        for (ptrdiff_t d = 0; d < m_n; ++d) {
            if (d == 0) {
                m_prefactor(0) = 0.0;
            }
            else {
                m_prefactor(d) = m_k / std::pow(d, m_alpha + 1.0);
            }
        }
    }
 
    template <class T>
    void force(const T& u, T& f)
    {
        for (ptrdiff_t p = 0; p < m_n; ++p) {
            double fp = 0.0;
            double up = u(p);
            for (ptrdiff_t i = 0; i < m_n; ++i) {
                if (i == p) {
                    continue;
                }
                ptrdiff_t d = std::abs(i - p);
                if (d > m_m) {
                    d = m_n - d;
                }
                fp += (u(i) - up) * m_prefactor(d);
            }
            f(p) = fp;
        }
    }
};
 
template <size_t rank>
class RandomNormalForcing {
protected:
    std::array<size_type, rank> m_shape; 
    size_type m_N; 
    array_type::tensor<double, rank> m_f_thermal; 
    double m_mean; 
    double m_stddev; 
    array_type::tensor<ptrdiff_t, rank> m_next; 
    array_type::tensor<ptrdiff_t, rank> m_dinc; 
    prrng::pcg32 m_rng; 
 
public:
    RandomNormalForcing() = default;
 
    template <class T, class S>
    RandomNormalForcing(
        const S& shape,
        double mean,
        double stddev,
        uint64_t seed,
        const T& dinc_init,
        const T& dinc
    )
    {
        std::copy(shape.cbegin(), shape.cend(), m_shape.begin());
        m_N = std::accumulate(m_shape.cbegin(), m_shape.cend(), 1, std::multiplies<size_type>{});
        m_f_thermal = xt::zeros<double>(m_shape);
        m_rng.seed(seed);
        m_mean = mean;
        m_stddev = stddev;
        m_next = dinc_init;
        m_dinc = dinc;
    }
 
    template <class T, class S>
    void force(const T& u, T& f, S inc)
    {
        (void)(u);
 
        for (size_type p = 0; p < m_N; ++p) {
            if (inc >= m_next(p)) {
                m_f_thermal.flat(p) = m_rng.normal(m_mean, m_stddev);
                m_next.flat(p) += m_dinc.flat(p);
            }
        }
 
        std::copy(m_f_thermal.begin(), m_f_thermal.end(), f.begin());
    }
 
    uint64_t state() const
    {
        return m_rng.state();
    }
 
    void set_state(uint64_t state)
    {
        m_rng.restore(state);
    }
 
    const auto& f_thermal() const
    {
        return m_f_thermal;
    }
 
    void set_f_thermal(const array_type::tensor<double, rank>& f_thermal)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(f_thermal, m_f_thermal.shape()));
        m_f_thermal = f_thermal;
    }
 
    const auto& next()
    {
        return m_next;
    }
 
    void set_next(const array_type::tensor<ptrdiff_t, rank>& next)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(next, m_next.shape()));
        m_next = next;
    }
};
 
class Overdamped {
public:
    Overdamped()
    {
    }
};
 
class None {
public:
    None()
    {
    }
};
 
template <
    size_t rank,
    class Potential,
    class Generator,
    class Interactions = void,
    class External = void,
    class Minimisation = void>
class System {
protected:
    size_type m_N; 
    array_type::tensor<double, rank> m_f; 
    array_type::tensor<double, rank> m_f_thermal; 
    array_type::tensor<double, rank> m_f_potential; 
    array_type::tensor<double, rank> m_f_interactions; 
    array_type::tensor<double, rank> m_f_frame; 
    array_type::tensor<double, rank> m_f_damping; 
    array_type::tensor<double, rank> m_u; 
    array_type::tensor<double, rank> m_v; 
    array_type::tensor<double, rank> m_a; 
    array_type::tensor<double, rank> m_v_n; 
    array_type::tensor<double, rank> m_a_n; 
    ptrdiff_t m_inc = 0; 
    ptrdiff_t m_qs_inc_first = 0; 
    ptrdiff_t m_qs_inc_last = 0; 
    double m_dt; 
    double m_eta; 
    double m_m; 
    double m_inv_m; 
    double m_mu; 
    double m_k_frame; 
    double m_u_frame = 0.0; 
    Potential* m_potential; 
    Generator* m_chunk; 
    Interactions* m_interactions; 
    External* m_external; 
 
protected:
    void initSystem(
        double m,
        double eta,
        double k_frame,
        double mu,
        double dt,
        Potential* potential,
        Generator* chunk,
        Interactions* interactions = nullptr,
        External* external = nullptr
    )
    {
        m_N = static_cast<size_type>(chunk->generators().size());
        m_m = m;
        m_inv_m = 1.0 / m;
        m_eta = eta;
        m_mu = mu;
        m_k_frame = k_frame;
        m_u_frame = 0.0;
        m_dt = dt;
 
        m_f = xt::zeros<double>(chunk->generators().shape());
        m_f_potential = xt::zeros<double>(chunk->generators().shape());
        m_f_interactions = xt::zeros<double>(chunk->generators().shape());
        m_f_frame = xt::zeros<double>(chunk->generators().shape());
        m_f_damping = xt::zeros<double>(chunk->generators().shape());
 
        if constexpr (!std::is_same<External, void>::value) {
            m_f_thermal = xt::zeros<double>(chunk->generators().shape());
        }
 
        m_u = xt::zeros<double>(chunk->generators().shape());
        m_v = xt::zeros<double>(chunk->generators().shape());
        m_a = xt::zeros<double>(chunk->generators().shape());
        m_v_n = xt::zeros<double>(chunk->generators().shape());
        m_a_n = xt::zeros<double>(chunk->generators().shape());
 
        m_potential = potential;
        m_chunk = chunk;
        m_interactions = interactions;
        m_external = external;
 
        this->refresh();
    }
 
public:
    const auto& chunk() const
    {
        return m_chunk;
    }
 
    auto size() const
    {
        return m_chunk->generators().size();
    }
 
    const auto& shape() const
    {
        return m_chunk->generators().shape();
    }
 
    auto dt() const
    {
        return m_dt;
    }
 
    auto mu() const
    {
        return m_mu;
    }
 
    auto eta() const
    {
        return m_eta;
    }
 
    auto m() const
    {
        return m_m;
    }
 
    auto k_frame() const
    {
        return m_k_frame;
    }
 
    const auto& external() const
    {
        return m_external;
    }
 
    void set_t(double arg)
    {
        m_inc = static_cast<decltype(m_inc)>(std::round(arg / m_dt));
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::allclose(this->t(), arg));
    }
 
    void set_inc(ptrdiff_t arg)
    {
        m_inc = arg;
        m_qs_inc_first = arg;
        m_qs_inc_last = arg;
        this->updated_inc();
    }
 
    void set_u_frame(double arg)
    {
        m_u_frame = arg;
        this->computeForceFrame();
        this->computeForce();
    }
 
    double u_frame() const
    {
        return m_u_frame;
    }
 
    void set_u(const array_type::tensor<double, rank>& arg)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(arg, m_u.shape()));
        xt::noalias(m_u) = arg;
        this->updated_u();
    }
 
    void set_v(const array_type::tensor<double, rank>& arg)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(arg, m_v.shape()));
        xt::noalias(m_v) = arg;
        this->updated_v();
    }
 
    void set_a(const array_type::tensor<double, rank>& arg)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(arg, m_a.shape()));
        xt::noalias(m_a) = arg;
    }
 
    void refresh()
    {
        this->updated_u();
        this->updated_v();
        this->updated_inc();
    }
 
protected:
    void computeForce()
    {
        if constexpr (std::is_same<External, void>::value) {
            xt::noalias(m_f) = m_f_frame + m_f_potential + m_f_interactions + m_f_damping;
        }
        else {
            xt::noalias(m_f) =
                m_f_frame + m_f_potential + m_f_interactions + m_f_damping + m_f_thermal;
        }
    }
 
    void computeForcePotential()
    {
        m_potential->force(m_u, m_f_potential);
    }
 
    void computeForceInteractions()
    {
        if constexpr (!std::is_same<Interactions, void>::value) {
            m_interactions->force(m_u, m_f_interactions);
        }
    }
 
    void computeForceFrame()
    {
        xt::noalias(m_f_frame) = m_k_frame * (m_u_frame - m_u);
    }
 
    void computeForceDamping()
    {
        xt::noalias(m_f_damping) = -m_eta * m_v;
    }
 
    void updated_inc()
    {
        if constexpr (!std::is_same<External, void>::value) {
            m_external->force(m_u, m_f_thermal, m_inc);
        }
        this->computeForce();
    }
 
    void updated_u()
    {
        this->computeForcePotential();
        this->computeForceInteractions();
        this->computeForceFrame();
        this->computeForce();
    }
 
    void updated_v()
    {
        this->computeForceDamping();
        this->computeForce();
    }
 
public:
    const auto& u() const
    {
        return m_u;
    }
 
    const auto& v() const
    {
        return m_v;
    }
 
    const auto& a() const
    {
        return m_a;
    }
 
    const auto& f() const
    {
        return m_f;
    }
 
    const auto& f_potential() const
    {
        return m_f_potential;
    }
 
    const auto& f_frame() const
    {
        return m_f_frame;
    }
 
    const auto& f_interactions() const
    {
        return m_f_interactions;
    }
 
    const auto& f_damping() const
    {
        return m_f_damping;
    }
 
    auto t() const
    {
        return static_cast<double>(m_inc) * m_dt;
    }
 
    auto inc() const
    {
        return m_inc;
    }
 
    double temperature() const
    {
        return 0.5 * m_m * xt::norm_sq(m_v)() / static_cast<double>(m_N);
    }
 
    double residual() const
    {
        double r_fres = xt::norm_l2(m_f)();
        double r_fext = xt::norm_l2(m_f_frame)();
        if (r_fext != 0.0) {
            return r_fres / r_fext;
        }
        return r_fres;
    }
 
    void quench()
    {
        m_v.fill(0.0);
        m_a.fill(0.0);
        this->updated_v();
    }
 
    void timeStep()
    {
        m_inc++;
        if constexpr (!std::is_same<External, void>::value) {
            this->updated_inc();
        }
 
        xt::noalias(m_v_n) = m_v;
        xt::noalias(m_a_n) = m_a;
 
        xt::noalias(m_u) = m_u + m_dt * m_v + 0.5 * m_dt * m_dt * m_a;
        this->updated_u();
 
        xt::noalias(m_v) = m_v_n + m_dt * m_a_n;
        this->updated_v();
 
        xt::noalias(m_a) = m_f * m_inv_m;
 
        xt::noalias(m_v) = m_v_n + 0.5 * m_dt * (m_a_n + m_a);
        this->updated_v();
 
        xt::noalias(m_a) = m_f * m_inv_m;
 
        xt::noalias(m_v) = m_v_n + 0.5 * m_dt * (m_a_n + m_a);
        this->updated_v();
 
        xt::noalias(m_a) = m_f * m_inv_m;
 
        if (xt::any(xt::isnan(m_u))) {
            throw std::runtime_error("NaN entries found");
        }
    }
 
    void timeSteps(size_t n)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(n + 1 < std::numeric_limits<long>::max());
        for (size_t step = 0; step < n; ++step) {
            this->timeStep();
        }
    }
 
    size_t timeStepsUntilEvent(double tol = 1e-5, size_t niter_tol = 10, size_t max_iter = 1e9)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(tol < 1.0);
        FRICTIONQPOTSPRINGBLOCK_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
 
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
        auto i_n = m_chunk->index_at_align();
        size_t step;
 
        for (step = 1; step < max_iter + 1; ++step) {
 
            this->timeStep();
 
            if (xt::any(xt::not_equal(m_chunk->index_at_align(), 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;
    }
 
    void flowSteps(size_t n, double v_frame)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(n + 1 < std::numeric_limits<long>::max());
 
        for (size_t step = 0; step < n; ++step) {
            m_u_frame += v_frame * m_dt;
            this->timeStep();
        }
    }
 
    size_t minimise(
        double tol = 1e-5,
        size_t niter_tol = 10,
        size_t max_iter = 1e9,
        bool time_activity = false,
        bool max_iter_is_error = true
    )
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(tol < 1.0);
        FRICTIONQPOTSPRINGBLOCK_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
 
        size_t step;
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList res(niter_tol);
 
        if constexpr (std::is_same<Minimisation, None>::value) {
            throw std::runtime_error("Minimisation not implementated");
        }
        else if constexpr (std::is_same<Minimisation, Overdamped>::value) {
            static_assert(std::is_same<Interactions, Laplace1d>::value);
            FRICTIONQPOTSPRINGBLOCK_ASSERT(!time_activity);
            (void)(time_activity);
 
            double uneigh;
            double u;
            double umin;
            ptrdiff_t i;
            ptrdiff_t j;
            m_qs_inc_first = m_inc; // unused
            m_qs_inc_last = m_inc; // unused
            double k = m_interactions->k();
 
            for (step = 1; step < max_iter + 1; ++step) {
 
                // "misuse" unused variable
                xt::noalias(m_v_n) = m_u;
 
                for (size_type p = 0; p < m_N; ++p) {
 
                    if (p == 0) {
                        uneigh = m_v_n.back() + m_v_n(1);
                    }
                    else if (p == m_N - 1) {
                        uneigh = m_v_n(m_N - 2) + m_v_n.front();
                    }
                    else {
                        uneigh = m_v_n(p - 1) + m_v_n(p + 1);
                    }
 
                    i = m_chunk->chunk_index_at_align()(p);
                    auto* y = &m_chunk->data()(p, 0);
 
                    while (true) {
                        umin = 0.5 * (*(y + i) + *(y + i + 1));
                        u = (k * uneigh + m_k_frame * m_u_frame + m_mu * umin) /
                            (2 * k + m_k_frame + m_mu);
                        m_chunk->align(p, u);
                        j = m_chunk->chunk_index_at_align()(p);
                        if (j == i) {
                            break;
                        }
                        i = j;
                    }
                    m_u(p) = u;
                    m_f_frame(p) = m_k_frame * (m_u_frame - u);
                    m_f_potential(p) = m_mu * (umin - u);
                }
 
                this->computeForceInteractions();
                xt::noalias(m_f) = m_f_potential + m_f_interactions + m_f_frame;
                res.roll_insert(this->residual());
 
                if ((res.descending() && res.all_less(tol)) || res.all_less(tol2)) {
                    this->quench(); // no dynamics are run: make sure that the user is not confused
                    return 0;
                }
            }
        }
        else {
            array_type::tensor<ptrdiff_t, rank> i_n;
            long s = 0;
            long s_n = 0;
            bool init = true;
 
            if (time_activity) {
                i_n = m_chunk->index_at_align();
            }
 
            for (step = 1; step < max_iter + 1; ++step) {
                this->timeStep();
                res.roll_insert(this->residual());
 
                if (time_activity) {
                    s = xt::sum(xt::abs(m_chunk->index_at_align() - 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 ((res.descending() && res.all_less(tol)) || res.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;
    }
 
    size_t minimise_truncate(
        array_type::tensor<ptrdiff_t, rank> i_n,
        size_t A_truncate = 0,
        size_t S_truncate = 0,
        double tol = 1e-5,
        size_t niter_tol = 10,
        size_t max_iter = 1e9,
        bool time_activity = true,
        bool max_iter_is_error = true
    )
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(tol < 1.0);
        FRICTIONQPOTSPRINGBLOCK_ASSERT(max_iter + 1 < std::numeric_limits<long>::max());
        FRICTIONQPOTSPRINGBLOCK_ASSERT(xt::has_shape(i_n, m_u.shape()));
        FRICTIONQPOTSPRINGBLOCK_ASSERT(time_activity);
        (void)(time_activity);
 
        double tol2 = tol * tol;
        GooseFEM::Iterate::StopList residuals(niter_tol);
 
        long s = 0;
        long s_n = 0;
        bool init = true;
        size_t step;
 
        for (step = 1; step < max_iter + 1; ++step) {
 
            this->timeStep();
            residuals.roll_insert(this->residual());
 
            size_t a = (size_t)xt::sum(xt::not_equal(m_chunk->index_at_align(), i_n))();
            s = xt::sum(xt::abs(m_chunk->index_at_align() - 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 (A_truncate > 0 && a >= A_truncate) {
                return step;
            }
 
            if (S_truncate > 0 && s >= (long)S_truncate) {
                return step;
            }
        }
 
        if (max_iter_is_error) {
            throw std::runtime_error("No convergence found");
        }
 
        return step;
    }
 
    double maxUniformDisplacement(int direction)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(direction == 1 || direction == -1);
        return m_potential->maxUniformDisplacement(m_u, direction);
    };
 
    double eventDrivenStep(double eps, bool kick, int direction = 1)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT(direction == 1 || direction == -1);
 
        if (direction > 0 && !kick) {
            double du = this->maxUniformDisplacement(direction);
            if (du < 0.5 * eps) {
                return 0.0;
            }
            return this->advanceUniformly(du - 0.5 * eps, false);
        }
 
        if (direction > 0 && kick) {
            return this->advanceUniformly(eps, false);
        }
 
        // direction < 0
 
        if (!kick) {
            double du = this->maxUniformDisplacement(direction);
            if (du < 0.5 * eps) {
                return 0.0;
            }
            return this->advanceUniformly(0.5 * eps - du, false);
        }
 
        return this->advanceUniformly(-eps, false);
    }
 
    void trigger(size_t p, double eps, int direction = 1)
    {
        FRICTIONQPOTSPRINGBLOCK_ASSERT((size_type)p < m_N);
        return m_potential->trigger(m_u, p, eps, direction);
        this->updated_u();
    };
 
    void advanceToFixedForce(double f_frame, bool allow_plastic = false)
    {
        auto i_n = m_chunk->index_at_align();
        this->advanceUniformly((f_frame - xt::mean(m_f_frame)()) / m_mu, false);
        FRICTIONQPOTSPRINGBLOCK_REQUIRE(
            allow_plastic || xt::all(xt::equal(m_chunk->index_at_align(), i_n))
        );
    }
 
protected:
    double advanceUniformly(double du, bool input_is_frame = true)
    {
        double du_particles;
        double du_frame;
 
        if (input_is_frame) {
            du_frame = du;
            du_particles = du * m_k_frame / (m_k_frame + m_mu);
        }
        else {
            du_particles = du;
            du_frame = du * (m_k_frame + m_mu) / m_k_frame;
        }
 
        m_u += du_particles;
        m_u_frame += du_frame;
        this->updated_u();
 
        if (input_is_frame) {
            return du_particles;
        }
 
        return du_frame;
    }
};
 
} // namespace detail
} // namespace FrictionQPotSpringBlock
 
#endif