Matrix.h

Go to the documentation of this file.

 
#ifndef GOOSEFEM_MATRIX_H
#define GOOSEFEM_MATRIX_H
 
#include "config.h"
 
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
 
namespace GooseFEM {
 
// forward declaration
template <class>
class MatrixSolver;
 
template <class D>
class MatrixSolverBase {
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    template <class M>
    void solve(M& A, const array_type::tensor<double, 2>& b, array_type::tensor<double, 2>& x)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_nodevec()));
        GOOSEFEM_ASSERT(xt::has_shape(x, A.shape_nodevec()));
        derived_cast().solve_nodevec_impl(A, b, x);
    }
 
    template <class M>
    void solve(M& A, const array_type::tensor<double, 1>& b, array_type::tensor<double, 1>& x)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_dofval()));
        GOOSEFEM_ASSERT(xt::has_shape(x, A.shape_dofval()));
        derived_cast().solve_dofval_impl(A, b, x);
    }
};
 
template <class D>
class MatrixSolverSingleBase {
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    template <class M>
    array_type::tensor<double, 2> Solve(M& A, const array_type::tensor<double, 2>& b)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_nodevec()));
        array_type::tensor<double, 2> x = xt::empty_like(b);
        derived_cast().solve_nodevec_impl(A, b, x);
        return x;
    }
 
    template <class M>
    array_type::tensor<double, 1> Solve(M& A, const array_type::tensor<double, 1>& b)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_dofval()));
        array_type::tensor<double, 1> x = xt::empty_like(b);
        derived_cast().solve_dofval_impl(A, b, x);
        return x;
    }
};
 
template <class D>
class MatrixSolverPartitionedBase {
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    template <class M>
    array_type::tensor<double, 2>
    Solve(M& A, const array_type::tensor<double, 2>& b, const array_type::tensor<double, 2>& x)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_nodevec()));
        GOOSEFEM_ASSERT(xt::has_shape(x, A.shape_nodevec()));
        array_type::tensor<double, 2> ret = x;
        derived_cast().solve_nodevec_impl(A, b, ret);
        return ret;
    }
 
    template <class M>
    array_type::tensor<double, 1>
    Solve(M& A, const array_type::tensor<double, 1>& b, const array_type::tensor<double, 1>& x)
    {
        GOOSEFEM_ASSERT(xt::has_shape(b, A.shape_dofval()));
        GOOSEFEM_ASSERT(xt::has_shape(x, A.shape_dofval()));
        array_type::tensor<double, 1> ret = x;
        derived_cast().solve_dofval_impl(A, b, ret);
        return ret;
    }
};
 
template <class D>
class MatrixBase {
protected:
    array_type::tensor<size_t, 2> m_conn; 
    array_type::tensor<size_t, 2> m_dofs; 
 
    size_t m_nelem; 
    size_t m_nne; 
    size_t m_nnode; 
    size_t m_ndim; 
    size_t m_ndof; 
 
    bool m_changed = true; 
 
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    size_t nelem() const
    {
        return derived_cast().m_nelem;
    }
 
    size_t nne() const
    {
        return derived_cast().m_nne;
    }
 
    size_t nnode() const
    {
        return derived_cast().m_nnode;
    }
 
    size_t ndim() const
    {
        return derived_cast().m_ndim;
    }
 
    size_t ndof() const
    {
        return derived_cast().m_ndof;
    }
 
    const array_type::tensor<size_t, 2>& dofs() const
    {
        return derived_cast().m_dofs;
    }
 
    const array_type::tensor<size_t, 2>& conn() const
    {
        return derived_cast().m_conn;
    }
 
    std::array<size_t, 1> shape_dofval() const
    {
        return std::array<size_t, 1>{derived_cast().m_ndof};
    }
 
    std::array<size_t, 2> shape_nodevec() const
    {
        return std::array<size_t, 2>{derived_cast().m_nnode, derived_cast().m_ndim};
    }
 
    std::array<size_t, 3> shape_elemmat() const
    {
        return std::array<size_t, 3>{
            derived_cast().m_nelem,
            derived_cast().m_nne * derived_cast().m_ndim,
            derived_cast().m_nne * derived_cast().m_ndim
        };
    }
 
    template <class T>
    void assemble(const T& elemmat)
    {
        GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat()));
        derived_cast().assemble_impl(elemmat);
    }
 
    array_type::tensor<double, 2> Todense() const
    {
        size_t ndof = derived_cast().m_ndof;
        array_type::tensor<double, 2> ret = xt::empty<double>({ndof, ndof});
        derived_cast().todense_impl(ret);
        return ret;
    }
 
    template <class T>
    void todense(T& ret) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(ret, {derived_cast().m_ndof, derived_cast().m_ndof}));
        derived_cast().todense_impl(ret);
    }
 
    array_type::tensor<double, 2> Dot(const array_type::tensor<double, 2>& x) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_nodevec()));
        array_type::tensor<double, 2> b = xt::empty_like(x);
        derived_cast().dot_nodevec_impl(x, b);
        return b;
    }
 
    array_type::tensor<double, 1> Dot(const array_type::tensor<double, 1>& x) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_dofval()));
        array_type::tensor<double, 1> b = xt::empty_like(x);
        derived_cast().dot_dofval_impl(x, b);
        return b;
    }
 
    void dot(const array_type::tensor<double, 2>& x, array_type::tensor<double, 2>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_nodevec()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_nodevec()));
        derived_cast().dot_nodevec_impl(x, b);
    }
 
    void dot(const array_type::tensor<double, 1>& x, array_type::tensor<double, 1>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_dofval()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_dofval()));
        derived_cast().dot_dofval_impl(x, b);
    }
};
 
template <class D>
class MatrixPartitionedBase : public MatrixBase<D> {
protected:
    array_type::tensor<size_t, 1> m_iiu; 
    array_type::tensor<size_t, 1> m_iip; 
 
    size_t m_nnu; 
    size_t m_nnp; 
 
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    size_t nnu() const
    {
        return derived_cast().m_nnu;
    }
 
    size_t nnp() const
    {
        return derived_cast().m_nnp;
    }
 
    const array_type::tensor<size_t, 1>& iiu() const
    {
        return derived_cast().m_iiu;
    }
 
    const array_type::tensor<size_t, 1>& iip() const
    {
        return derived_cast().m_iip;
    }
 
    array_type::tensor<double, 2>
    Reaction(const array_type::tensor<double, 2>& x, const array_type::tensor<double, 2>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_nodevec()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_nodevec()));
        array_type::tensor<double, 2> ret = b;
        derived_cast().reaction_nodevec_impl(x, ret);
        return ret;
    }
 
    array_type::tensor<double, 1>
    Reaction(const array_type::tensor<double, 1>& x, const array_type::tensor<double, 1>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_dofval()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_dofval()));
        array_type::tensor<double, 1> ret = b;
        derived_cast().reaction_dofval_impl(x, ret);
        return ret;
    }
 
    void reaction(const array_type::tensor<double, 2>& x, array_type::tensor<double, 2>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_nodevec()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_nodevec()));
        derived_cast().reaction_nodevec_impl(x, b);
    }
 
    void reaction(const array_type::tensor<double, 1>& x, array_type::tensor<double, 1>& b) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(x, this->shape_dofval()));
        GOOSEFEM_ASSERT(xt::has_shape(b, this->shape_dofval()));
        derived_cast().reaction_dofval_impl(x, b);
    }
 
    array_type::tensor<double, 1> Reaction_p(
        const array_type::tensor<double, 1>& x_u,
        const array_type::tensor<double, 1>& x_p
    ) const
    {
        array_type::tensor<double, 1> b_p = xt::empty<double>({m_nnp});
        derived_cast().reaction_p_impl(x_u, x_p, b_p);
        return b_p;
    }
 
    void reaction_p(
        const array_type::tensor<double, 1>& x_u,
        const array_type::tensor<double, 1>& x_p,
        array_type::tensor<double, 1>& b_p
    ) const
    {
        derived_cast().reaction_p_impl(x_u, x_p, b_p);
    }
};
 
template <class D>
class MatrixPartitionedTyingsBase : public MatrixPartitionedBase<D> {
protected:
    array_type::tensor<size_t, 1> m_iii; 
    array_type::tensor<size_t, 1> m_iid; 
 
    size_t m_nni; 
    size_t m_nnd; 
 
public:
    using derived_type = D;
 
private:
    auto derived_cast() -> derived_type&
    {
        return *static_cast<derived_type*>(this);
    }
 
    auto derived_cast() const -> const derived_type&
    {
        return *static_cast<const derived_type*>(this);
    }
 
public:
    size_t nni() const
    {
        return derived_cast().m_nni;
    }
 
    size_t nnd() const
    {
        return derived_cast().m_nnd;
    }
 
    const array_type::tensor<size_t, 1>& iii() const
    {
        return derived_cast().m_iii;
    }
 
    const array_type::tensor<size_t, 1>& iid() const
    {
        return derived_cast().m_iid;
    }
};
 
class Matrix : public MatrixBase<Matrix> {
private:
    friend MatrixBase<Matrix>;
 
private:
    bool m_changed = true; 
 
    Eigen::SparseMatrix<double> m_A; 
 
    std::vector<Eigen::Triplet<double>> m_T; 
 
    template <class>
    friend class MatrixSolver;
 
public:
    Matrix() = default;
 
    Matrix(const array_type::tensor<size_t, 2>& conn, const array_type::tensor<size_t, 2>& dofs)
    {
        m_conn = conn;
        m_dofs = dofs;
        m_nelem = m_conn.shape(0);
        m_nne = m_conn.shape(1);
        m_nnode = m_dofs.shape(0);
        m_ndim = m_dofs.shape(1);
        m_ndof = xt::amax(m_dofs)() + 1;
        m_T.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_A.resize(m_ndof, m_ndof);
 
        GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
        GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
    }
 
    const Eigen::SparseMatrix<double>& data() const
    {
        return m_A;
    }
 
private:
    template <class T>
    void assemble_impl(const T& elemmat)
    {
        m_T.clear();
 
        for (size_t e = 0; e < m_nelem; ++e) {
            for (size_t m = 0; m < m_nne; ++m) {
                for (size_t i = 0; i < m_ndim; ++i) {
                    for (size_t n = 0; n < m_nne; ++n) {
                        for (size_t j = 0; j < m_ndim; ++j) {
                            m_T.push_back(Eigen::Triplet<double>(
                                m_dofs(m_conn(e, m), i),
                                m_dofs(m_conn(e, n), j),
                                elemmat(e, m * m_ndim + i, n * m_ndim + j)
                            ));
                        }
                    }
                }
            }
        }
 
        m_A.setFromTriplets(m_T.begin(), m_T.end());
        m_changed = true;
    }
 
public:
    void
    set(const array_type::tensor<size_t, 1>& rows,
        const array_type::tensor<size_t, 1>& cols,
        const array_type::tensor<double, 2>& matrix)
    {
        GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
        GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
        GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
        GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
 
        std::vector<Eigen::Triplet<double>> T;
 
        for (size_t i = 0; i < rows.size(); ++i) {
            for (size_t j = 0; j < cols.size(); ++j) {
                T.push_back(Eigen::Triplet<double>(rows(i), cols(j), matrix(i, j)));
            }
        }
 
        m_A.setFromTriplets(T.begin(), T.end());
        m_changed = true;
    }
 
    void
    add(const array_type::tensor<size_t, 1>& rows,
        const array_type::tensor<size_t, 1>& cols,
        const array_type::tensor<double, 2>& matrix)
    {
        GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
        GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
        GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
        GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
 
        std::vector<Eigen::Triplet<double>> T;
 
        Eigen::SparseMatrix<double> A(m_ndof, m_ndof);
 
        for (size_t i = 0; i < rows.size(); ++i) {
            for (size_t j = 0; j < cols.size(); ++j) {
                T.push_back(Eigen::Triplet<double>(rows(i), cols(j), matrix(i, j)));
            }
        }
 
        A.setFromTriplets(T.begin(), T.end());
        m_A += A;
        m_changed = true;
    }
 
private:
    template <class T>
    void todense_impl(T& ret) const
    {
        ret.fill(0.0);
 
        for (int k = 0; k < m_A.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_A, k); it; ++it) {
                ret(it.row(), it.col()) = it.value();
            }
        }
    }
 
    template <class T>
    void dot_nodevec_impl(const T& x, T& b) const
    {
        this->Eigen_asNode_dofval_nodevec(m_A * this->Eigen_AsDofs_nodevec(x), b);
    }
 
    template <class T>
    void dot_dofval_impl(const T& x, T& b) const
    {
        Eigen::Map<Eigen::VectorXd>(b.data(), b.size()).noalias() =
            m_A * Eigen::Map<const Eigen::VectorXd>(x.data(), x.size());
    }
 
private:
    template <class T>
    Eigen::VectorXd Eigen_AsDofs_nodevec(const T& nodevec) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
        Eigen::VectorXd dofval = Eigen::VectorXd::Zero(m_ndof, 1);
 
#pragma omp parallel for
        for (size_t m = 0; m < m_nnode; ++m) {
            for (size_t i = 0; i < m_ndim; ++i) {
                dofval(m_dofs(m, i)) = nodevec(m, i);
            }
        }
 
        return dofval;
    }
 
    template <class T>
    void Eigen_asNode_dofval_nodevec(const Eigen::VectorXd& dofval, T& nodevec) const
    {
        GOOSEFEM_ASSERT(static_cast<size_t>(dofval.size()) == m_ndof);
        GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
#pragma omp parallel for
        for (size_t m = 0; m < m_nnode; ++m) {
            for (size_t i = 0; i < m_ndim; ++i) {
                nodevec(m, i) = dofval(m_dofs(m, i));
            }
        }
    }
};
 
template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
class MatrixSolver : public MatrixSolverBase<MatrixSolver<Solver>>,
                     public MatrixSolverSingleBase<MatrixSolver<Solver>> {
private:
    friend MatrixSolverBase<MatrixSolver<Solver>>;
    friend MatrixSolverSingleBase<MatrixSolver<Solver>>;
 
public:
    MatrixSolver() = default;
 
private:
    template <class T>
    void solve_nodevec_impl(Matrix& A, const T& b, T& x)
    {
        this->factorize(A);
        Eigen::VectorXd X = m_solver.solve(A.Eigen_AsDofs_nodevec(b));
        A.Eigen_asNode_dofval_nodevec(X, x);
    }
 
    template <class T>
    void solve_dofval_impl(Matrix& A, const T& b, T& x)
    {
        this->factorize(A);
        Eigen::Map<Eigen::VectorXd>(x.data(), x.size()).noalias() =
            m_solver.solve(Eigen::Map<const Eigen::VectorXd>(b.data(), A.m_ndof));
    }
 
private:
    Solver m_solver; 
    bool m_factor = true; 
 
    void factorize(Matrix& A)
    {
        if (!A.m_changed && !m_factor) {
            return;
        }
        m_solver.compute(A.m_A);
        m_factor = false;
        A.m_changed = false;
    }
};
 
} // namespace GooseFEM
 
#endif