MatrixPartitionedTyings.h

Go to the documentation of this file.

 
#ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
#define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
 
#include "Matrix.h"
#include "config.h"
 
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
 
namespace GooseFEM {
 
// forward declaration
template <class>
class MatrixPartitionedTyingsSolver;
 
class MatrixPartitionedTyings : public MatrixPartitionedTyingsBase<MatrixPartitionedTyings> {
private:
    friend MatrixBase<MatrixPartitionedTyings>;
    friend MatrixPartitionedBase<MatrixPartitionedTyings>;
    friend MatrixPartitionedTyingsBase<MatrixPartitionedTyings>;
 
protected:
    Eigen::SparseMatrix<double> m_Auu; 
    Eigen::SparseMatrix<double> m_Aup; 
    Eigen::SparseMatrix<double> m_Apu; 
    Eigen::SparseMatrix<double> m_App; 
    Eigen::SparseMatrix<double> m_Aud; 
    Eigen::SparseMatrix<double> m_Apd; 
    Eigen::SparseMatrix<double> m_Adu; 
    Eigen::SparseMatrix<double> m_Adp; 
    Eigen::SparseMatrix<double> m_Add; 
    Eigen::SparseMatrix<double> m_ACuu; 
    Eigen::SparseMatrix<double> m_ACup; 
    Eigen::SparseMatrix<double> m_ACpu; 
    Eigen::SparseMatrix<double> m_ACpp; 
    std::vector<Eigen::Triplet<double>> m_Tuu; 
    std::vector<Eigen::Triplet<double>> m_Tup; 
    std::vector<Eigen::Triplet<double>> m_Tpu; 
    std::vector<Eigen::Triplet<double>> m_Tpp; 
    std::vector<Eigen::Triplet<double>> m_Tud; 
    std::vector<Eigen::Triplet<double>> m_Tpd; 
    std::vector<Eigen::Triplet<double>> m_Tdu; 
    std::vector<Eigen::Triplet<double>> m_Tdp; 
    std::vector<Eigen::Triplet<double>> m_Tdd; 
    Eigen::SparseMatrix<double> m_Cdu; 
    Eigen::SparseMatrix<double> m_Cdp; 
    Eigen::SparseMatrix<double> m_Cud; 
    Eigen::SparseMatrix<double> m_Cpd; 
 
    // grant access to solver class
    template <class>
    friend class MatrixPartitionedTyingsSolver;
 
public:
    MatrixPartitionedTyings() = default;
 
    MatrixPartitionedTyings(
        const array_type::tensor<size_t, 2>& conn,
        const array_type::tensor<size_t, 2>& dofs,
        const Eigen::SparseMatrix<double>& Cdu,
        const Eigen::SparseMatrix<double>& Cdp
    )
    {
        GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows());
 
        m_conn = conn;
        m_dofs = dofs;
        m_Cdu = Cdu;
        m_Cdp = Cdp;
        m_nnu = static_cast<size_t>(m_Cdu.cols());
        m_nnp = static_cast<size_t>(m_Cdp.cols());
        m_nnd = static_cast<size_t>(m_Cdp.rows());
        m_nni = m_nnu + m_nnp;
        m_ndof = m_nni + m_nnd;
        m_iiu = xt::arange<size_t>(m_nnu);
        m_iip = xt::arange<size_t>(m_nnu, m_nnu + m_nnp);
        m_iii = xt::arange<size_t>(m_nni);
        m_iid = xt::arange<size_t>(m_nni, m_nni + m_nnd);
        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_Cud = m_Cdu.transpose();
        m_Cpd = m_Cdp.transpose();
        m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tud.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tpd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tdu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tdp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Tdd.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
        m_Auu.resize(m_nnu, m_nnu);
        m_Aup.resize(m_nnu, m_nnp);
        m_Apu.resize(m_nnp, m_nnu);
        m_App.resize(m_nnp, m_nnp);
        m_Aud.resize(m_nnu, m_nnd);
        m_Apd.resize(m_nnp, m_nnd);
        m_Adu.resize(m_nnd, m_nnu);
        m_Adp.resize(m_nnd, m_nnp);
        m_Add.resize(m_nnd, m_nnd);
 
        GOOSEFEM_ASSERT(m_ndof <= m_nnode * m_ndim);
        GOOSEFEM_ASSERT(m_ndof == xt::amax(m_dofs)() + 1);
    }
 
    const Eigen::SparseMatrix<double>& data_uu() const
    {
        return m_Auu;
    }
 
    const Eigen::SparseMatrix<double>& data_up() const
    {
        return m_Aup;
    }
 
    const Eigen::SparseMatrix<double>& data_pu() const
    {
        return m_Apu;
    }
 
    const Eigen::SparseMatrix<double>& data_pp() const
    {
        return m_App;
    }
 
    const Eigen::SparseMatrix<double>& data_ud() const
    {
        return m_Aud;
    }
 
    const Eigen::SparseMatrix<double>& data_pd() const
    {
        return m_Apd;
    }
 
    const Eigen::SparseMatrix<double>& data_du() const
    {
        return m_Adu;
    }
 
    const Eigen::SparseMatrix<double>& data_dp() const
    {
        return m_Adp;
    }
 
    const Eigen::SparseMatrix<double>& data_dd() const
    {
        return m_Add;
    }
 
private:
    template <class T>
    void assemble_impl(const T& elemmat)
    {
        GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
 
        m_Tuu.clear();
        m_Tup.clear();
        m_Tpu.clear();
        m_Tpp.clear();
        m_Tud.clear();
        m_Tpd.clear();
        m_Tdu.clear();
        m_Tdp.clear();
        m_Tdd.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) {
 
                    size_t di = m_dofs(m_conn(e, m), i);
 
                    for (size_t n = 0; n < m_nne; ++n) {
                        for (size_t j = 0; j < m_ndim; ++j) {
 
                            size_t dj = m_dofs(m_conn(e, n), j);
 
                            if (di < m_nnu && dj < m_nnu) {
                                m_Tuu.push_back(Eigen::Triplet<double>(
                                    di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (di < m_nnu && dj < m_nni) {
                                m_Tup.push_back(Eigen::Triplet<double>(
                                    di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (di < m_nnu) {
                                m_Tud.push_back(Eigen::Triplet<double>(
                                    di, dj - m_nni, elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (di < m_nni && dj < m_nnu) {
                                m_Tpu.push_back(Eigen::Triplet<double>(
                                    di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (di < m_nni && dj < m_nni) {
                                m_Tpp.push_back(Eigen::Triplet<double>(
                                    di - m_nnu,
                                    dj - m_nnu,
                                    elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (di < m_nni) {
                                m_Tpd.push_back(Eigen::Triplet<double>(
                                    di - m_nnu,
                                    dj - m_nni,
                                    elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (dj < m_nnu) {
                                m_Tdu.push_back(Eigen::Triplet<double>(
                                    di - m_nni, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else if (dj < m_nni) {
                                m_Tdp.push_back(Eigen::Triplet<double>(
                                    di - m_nni,
                                    dj - m_nnu,
                                    elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                            else {
                                m_Tdd.push_back(Eigen::Triplet<double>(
                                    di - m_nni,
                                    dj - m_nni,
                                    elemmat(e, m * m_ndim + i, n * m_ndim + j)
                                ));
                            }
                        }
                    }
                }
            }
        }
 
        m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end());
        m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end());
        m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end());
        m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end());
        m_Aud.setFromTriplets(m_Tud.begin(), m_Tud.end());
        m_Apd.setFromTriplets(m_Tpd.begin(), m_Tpd.end());
        m_Adu.setFromTriplets(m_Tdu.begin(), m_Tdu.end());
        m_Adp.setFromTriplets(m_Tdp.begin(), m_Tdp.end());
        m_Add.setFromTriplets(m_Tdd.begin(), m_Tdd.end());
        m_changed = true;
    }
 
    // todo: test
    template <class T>
    void todense_impl(T& ret) const
    {
        ret.fill(0.0);
 
        for (int k = 0; k < m_Auu.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Auu, k); it; ++it) {
                ret(m_iiu(it.row()), m_iiu(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Aup.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Aup, k); it; ++it) {
                ret(m_iiu(it.row()), m_iip(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Apu.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Apu, k); it; ++it) {
                ret(m_iip(it.row()), m_iiu(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_App.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_App, k); it; ++it) {
                ret(m_iip(it.row()), m_iip(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Adu.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Adu, k); it; ++it) {
                ret(m_iid(it.row()), m_iiu(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Adp.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Adp, k); it; ++it) {
                ret(m_iid(it.row()), m_iip(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Aud.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Aud, k); it; ++it) {
                ret(m_iiu(it.row()), m_iid(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Apd.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Apd, k); it; ++it) {
                ret(m_iip(it.row()), m_iid(it.col())) = it.value();
            }
        }
 
        for (int k = 0; k < m_Add.outerSize(); ++k) {
            for (Eigen::SparseMatrix<double>::InnerIterator it(m_Add, k); it; ++it) {
                ret(m_iid(it.row()), m_iid(it.col())) = it.value();
            }
        }
    }
 
    template <class T>
    void dot_nodevec_impl(const T& x, T& b) const
    {
        UNUSED(x);
        UNUSED(b);
        throw std::runtime_error("Not yet implemented");
    }
 
    template <class T>
    void dot_dofval_impl(const T& x, T& b) const
    {
        UNUSED(x);
        UNUSED(b);
        throw std::runtime_error("Not yet implemented");
    }
 
    template <class T>
    void reaction_nodevec_impl(const T& x, T& b) const
    {
        UNUSED(x);
        UNUSED(b);
        throw std::runtime_error("Not yet implemented");
    }
 
    template <class T>
    void reaction_dofval_impl(const T& x, T& b) const
    {
        UNUSED(x);
        UNUSED(b);
        throw std::runtime_error("Not yet implemented");
    }
 
    void reaction_p_impl(
        const array_type::tensor<double, 1>& x_u,
        const array_type::tensor<double, 1>& x_p,
        array_type::tensor<double, 1>& b_p
    ) const
    {
        UNUSED(x_u);
        UNUSED(x_p);
        UNUSED(b_p);
        throw std::runtime_error("Not yet implemented");
    }
 
private:
    // Convert arrays (Eigen version of VectorPartitioned, which contains public functions)
    Eigen::VectorXd AsDofs_u(const array_type::tensor<double, 1>& dofval) const
    {
        GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
        Eigen::VectorXd dofval_u(m_nnu, 1);
 
#pragma omp parallel for
        for (size_t d = 0; d < m_nnu; ++d) {
            dofval_u(d) = dofval(m_iiu(d));
        }
 
        return dofval_u;
    }
 
    Eigen::VectorXd AsDofs_u(const array_type::tensor<double, 2>& nodevec) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
        Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(m_nnu, 1);
 
#pragma omp parallel for
        for (size_t m = 0; m < m_nnode; ++m) {
            for (size_t i = 0; i < m_ndim; ++i) {
                if (m_dofs(m, i) < m_nnu) {
                    dofval_u(m_dofs(m, i)) = nodevec(m, i);
                }
            }
        }
 
        return dofval_u;
    }
 
    Eigen::VectorXd AsDofs_p(const array_type::tensor<double, 1>& dofval) const
    {
        GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
        Eigen::VectorXd dofval_p(m_nnp, 1);
 
#pragma omp parallel for
        for (size_t d = 0; d < m_nnp; ++d) {
            dofval_p(d) = dofval(m_iip(d));
        }
 
        return dofval_p;
    }
 
    Eigen::VectorXd AsDofs_p(const array_type::tensor<double, 2>& nodevec) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
        Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(m_nnp, 1);
 
#pragma omp parallel for
        for (size_t m = 0; m < m_nnode; ++m) {
            for (size_t i = 0; i < m_ndim; ++i) {
                if (m_dofs(m, i) >= m_nnu && m_dofs(m, i) < m_nni) {
                    dofval_p(m_dofs(m, i) - m_nnu) = nodevec(m, i);
                }
            }
        }
 
        return dofval_p;
    }
 
    Eigen::VectorXd AsDofs_d(const array_type::tensor<double, 1>& dofval) const
    {
        GOOSEFEM_ASSERT(dofval.size() == m_ndof);
 
        Eigen::VectorXd dofval_d(m_nnd, 1);
 
#pragma omp parallel for
        for (size_t d = 0; d < m_nnd; ++d) {
            dofval_d(d) = dofval(m_iip(d));
        }
 
        return dofval_d;
    }
 
    Eigen::VectorXd AsDofs_d(const array_type::tensor<double, 2>& nodevec) const
    {
        GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
 
        Eigen::VectorXd dofval_d = Eigen::VectorXd::Zero(m_nnd, 1);
 
#pragma omp parallel for
        for (size_t m = 0; m < m_nnode; ++m) {
            for (size_t i = 0; i < m_ndim; ++i) {
                if (m_dofs(m, i) >= m_nni) {
                    dofval_d(m_dofs(m, i) - m_nni) = nodevec(m, i);
                }
            }
        }
 
        return dofval_d;
    }
};
 
template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
class MatrixPartitionedTyingsSolver
    : public MatrixSolverBase<MatrixPartitionedTyingsSolver<Solver>>,
      public MatrixSolverPartitionedBase<MatrixPartitionedTyingsSolver<Solver>> {
private:
    friend MatrixSolverBase<MatrixPartitionedTyingsSolver<Solver>>;
    friend MatrixSolverPartitionedBase<MatrixPartitionedTyingsSolver<Solver>>;
 
public:
    MatrixPartitionedTyingsSolver() = default;
 
private:
    template <class T>
    void solve_nodevec_impl(MatrixPartitionedTyings& A, const T& b, T& x)
    {
        this->factorize(A);
 
        Eigen::VectorXd B_u = A.AsDofs_u(b);
        Eigen::VectorXd B_d = A.AsDofs_d(b);
        Eigen::VectorXd X_p = A.AsDofs_p(x);
 
        B_u += A.m_Cud * B_d;
 
        Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.m_ACup * X_p));
        Eigen::VectorXd X_d = A.m_Cdu * X_u + A.m_Cdp * X_p;
 
#pragma omp parallel for
        for (size_t m = 0; m < A.m_nnode; ++m) {
            for (size_t i = 0; i < A.m_ndim; ++i) {
                if (A.m_dofs(m, i) < A.m_nnu) {
                    x(m, i) = X_u(A.m_dofs(m, i));
                }
                else if (A.m_dofs(m, i) >= A.m_nni) {
                    x(m, i) = X_d(A.m_dofs(m, i) - A.m_nni);
                }
            }
        }
    }
 
    template <class T>
    void solve_dofval_impl(MatrixPartitionedTyings& A, const T& b, T& x)
    {
        this->factorize(A);
 
        Eigen::VectorXd B_u = A.AsDofs_u(b);
        Eigen::VectorXd B_d = A.AsDofs_d(b);
        Eigen::VectorXd X_p = A.AsDofs_p(x);
 
        Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.m_ACup * X_p));
        Eigen::VectorXd X_d = A.m_Cdu * X_u + A.m_Cdp * X_p;
 
#pragma omp parallel for
        for (size_t d = 0; d < A.m_nnu; ++d) {
            x(A.m_iiu(d)) = X_u(d);
        }
 
#pragma omp parallel for
        for (size_t d = 0; d < A.m_nnd; ++d) {
            x(A.m_iid(d)) = X_d(d);
        }
    }
 
public:
    array_type::tensor<double, 1> Solve_u(
        MatrixPartitionedTyings& A,
        const array_type::tensor<double, 1>& b_u,
        const array_type::tensor<double, 1>& b_d,
        const array_type::tensor<double, 1>& x_p
    )
    {
        array_type::tensor<double, 1> x_u = xt::empty<double>({A.m_nnu});
        this->solve_u(A, b_u, b_d, x_p, x_u);
        return x_u;
    }
 
    void solve_u(
        MatrixPartitionedTyings& A,
        const array_type::tensor<double, 1>& b_u,
        const array_type::tensor<double, 1>& b_d,
        const array_type::tensor<double, 1>& x_p,
        array_type::tensor<double, 1>& x_u
    )
    {
        UNUSED(b_d);
        GOOSEFEM_ASSERT(b_u.size() == A.m_nnu);
        GOOSEFEM_ASSERT(b_d.size() == A.m_nnd);
        GOOSEFEM_ASSERT(x_p.size() == A.m_nnp);
        GOOSEFEM_ASSERT(x_u.size() == A.m_nnu);
 
        this->factorize(A);
 
        Eigen::Map<Eigen::VectorXd>(x_u.data(), x_u.size()).noalias() =
            m_solver.solve(Eigen::VectorXd(
                Eigen::Map<const Eigen::VectorXd>(b_u.data(), b_u.size()) -
                A.m_ACup * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size())
            ));
    }
 
private:
    Solver m_solver; 
    bool m_factor = true; 
 
    void factorize(MatrixPartitionedTyings& A)
    {
        if (!A.m_changed && !m_factor) {
            return;
        }
 
        A.m_ACuu = A.m_Auu + A.m_Aud * A.m_Cdu + A.m_Cud * A.m_Adu + A.m_Cud * A.m_Add * A.m_Cdu;
 
        A.m_ACup = A.m_Aup + A.m_Aud * A.m_Cdp + A.m_Cud * A.m_Adp + A.m_Cud * A.m_Add * A.m_Cdp;
 
        // A.m_ACpu = A.m_Apu + A.m_Apd * A.m_Cdu + A.m_Cpd * A.m_Adu
        //     + A.m_Cpd * A.m_Add * A.m_Cdu;
 
        // A.m_ACpp = A.m_App + A.m_Apd * A.m_Cdp + A.m_Cpd * A.m_Adp
        //     + A.m_Cpd * A.m_Add * A.m_Cdp;
 
        m_solver.compute(A.m_ACuu);
        m_factor = false;
        A.m_changed = false;
    }
};
 
} // namespace GooseFEM
 
#endif