diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h index 7baf514..18a2fc7 100644 --- a/include/GooseFEM/MatrixPartitioned.h +++ b/include/GooseFEM/MatrixPartitioned.h @@ -1,294 +1,295 @@ /** Sparse matrix that is partitioned in: - unknown DOFs - prescribed DOFs \file MatrixPartitioned.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIXPARTITIONED_H #define GOOSEFEM_MATRIXPARTITIONED_H #include "config.h" #include "Matrix.h" #include #include #include namespace GooseFEM { // forward declaration template class MatrixPartitionedSolver; /** Sparse matrix partitioned in an unknown and a prescribed part. In particular: \f$ \begin{bmatrix} A_{uu} & A_{up} \\ A_{pu} & A_{pp} \end{bmatrix} \f$ See VectorPartitioned() for bookkeeping definitions. */ class MatrixPartitioned : public Matrix { public: MatrixPartitioned() = default; /** Constructor. \param conn connectivity [#nelem, #nne]. \param dofs DOFs per node [#nnode, #ndim]. \param iip prescribed DOFs [#nnp]. */ MatrixPartitioned( const xt::xtensor& conn, const xt::xtensor& dofs, const xt::xtensor& iip); /** Number of unknown DOFs. */ size_t nnu() const; /** Number of prescribed DOFs. */ size_t nnp() const; /** Unknown DOFs [#nnu]. */ xt::xtensor iiu() const; /** Prescribed DOFs [#nnp]. */ xt::xtensor iip() const; void assemble(const xt::xtensor& elemmat) override; void set( const xt::xtensor& rows, const xt::xtensor& cols, const xt::xtensor& matrix) override; void add( const xt::xtensor& rows, const xt::xtensor& cols, const xt::xtensor& matrix) override; void todense(xt::xtensor& ret) const override; void dot(const xt::xtensor& x, xt::xtensor& b) const override; void dot(const xt::xtensor& x, xt::xtensor& b) const override; /** Get right-hand-size for corresponding to the prescribed DOFs. \f$ b_p = A_{pu} * x_u + A_{pp} * x_p = A_{pp} * x_p \f$ and assemble them to the appropriate places in "nodevec". \param x "nodevec" [#nnode, #ndim]. \param b "nodevec" [#nnode, #ndim]. \return Copy of `b` with \f$ b_p \f$ overwritten. */ xt::xtensor Reaction( const xt::xtensor& x, const xt::xtensor& b) const; /** Same as Reaction(const xt::xtensor&, const xt::xtensor&), but inserting in-place. \param x "nodevec" [#nnode, #ndim]. \param b "nodevec" [#nnode, #ndim], \f$ b_p \f$ overwritten. */ void reaction( const xt::xtensor& x, xt::xtensor& b) const; /** Same as Reaction(const xt::xtensor&, const xt::xtensor&), but of "dofval" input and output. \param x "dofval" [#ndof]. \param b "dofval" [#ndof]. \return Copy of `b` with \f$ b_p \f$ overwritten. */ xt::xtensor Reaction( const xt::xtensor& x, const xt::xtensor& b) const; /** Same as Reaction(const xt::xtensor&, const xt::xtensor&), but inserting in-place. \param x "dofval" [#ndof]. \param b "dofval" [#ndof], \f$ b_p \f$ overwritten. */ void reaction( const xt::xtensor& x, xt::xtensor& b) const; /** Same as Reaction(const xt::xtensor&, const xt::xtensor&), but with partitioned input and output. \param x_u unknown "dofval" [#nnu]. \param x_p prescribed "dofval" [#nnp]. \return b_p prescribed "dofval" [#nnp]. */ xt::xtensor Reaction_p( const xt::xtensor& x_u, const xt::xtensor& x_p) const; /** Same as Reaction_p(const xt::xtensor&, const xt::xtensor&), but writing to preallocated output. \param x_u unknown "dofval" [#nnu]. \param x_p prescribed "dofval" [#nnp]. \param b_p (overwritten) prescribed "dofval" [#nnp]. */ void reaction_p( const xt::xtensor& x_u, const xt::xtensor& x_p, xt::xtensor& b_p) const; private: Eigen::SparseMatrix m_Auu; ///< The matrix. Eigen::SparseMatrix m_Aup; ///< The matrix. Eigen::SparseMatrix m_Apu; ///< The matrix. Eigen::SparseMatrix m_App; ///< The matrix. std::vector> m_Tuu; ///< Matrix entries. std::vector> m_Tup; ///< Matrix entries. std::vector> m_Tpu; ///< Matrix entries. std::vector> m_Tpp; ///< Matrix entries. xt::xtensor m_iiu; ///< See iiu() xt::xtensor m_iip; ///< See iip() size_t m_nnu; ///< See #nnu size_t m_nnp; ///< See #nnp /** Renumbered DOFs per node, such that iiu = arange(nnu) iip = nnu + arange(nnp) making is much simpler to slice. */ xt::xtensor m_part; // grant access to solver class template friend class MatrixPartitionedSolver; private: // Convert arrays (Eigen version of VectorPartitioned, which contains public functions) Eigen::VectorXd AsDofs_u(const xt::xtensor& dofval) const; Eigen::VectorXd AsDofs_u(const xt::xtensor& nodevec) const; Eigen::VectorXd AsDofs_p(const xt::xtensor& dofval) const; Eigen::VectorXd AsDofs_p(const xt::xtensor& nodevec) const; }; /** Solver for MatrixPartitioned(). The idea is that this solver class can be used to solve for multiple right-hand-sides using one factorisation. */ template >> class MatrixPartitionedSolver { public: MatrixPartitionedSolver() = default; /** Solve \f$ x_u = A_{uu}^{-1} (b_u - A_{up} * x_p) \f$. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b nodevec [nelem, ndim]. \param x nodevec [nelem, ndim], used to read \f$ x_p \f$. \return x nodevec [nelem, ndim], \f$ x_u \f$ filled, \f$ x_p \f$ copied. */ xt::xtensor Solve( MatrixPartitioned& A, const xt::xtensor& b, const xt::xtensor& x); /** Same as Solve(MatrixPartitioned&, const xt::xtensor&, const xt::xtensor&), but filling \f$ x_u \f$ in place. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b nodevec [nelem, ndim]. \param x nodevec [nelem, ndim], \f$ x_p \f$ read, \f$ x_u \f$ filled. */ void solve(MatrixPartitioned& A, const xt::xtensor& b, xt::xtensor& x); /** Same as Solve(MatrixPartitioned&, const xt::xtensor&, const xt::xtensor&), but for "dofval" input and output. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b dofval [ndof]. \param x dofval [ndof], used to read \f$ x_p \f$. \return x dofval [ndof], \f$ x_u \f$ filled, \f$ x_p \f$ copied. */ xt::xtensor Solve( MatrixPartitioned& A, const xt::xtensor& b, const xt::xtensor& x); /** Same as Solve(MatrixPartitioned&, const xt::xtensor&, const xt::xtensor&), but filling \f$ x_u \f$ in place. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b dofval [ndof]. \param x dofval [ndof], \f$ x_p \f$ read, \f$ x_u \f$ filled. */ void solve(MatrixPartitioned& A, const xt::xtensor& b, xt::xtensor& x); /** - Solve \f$ x_u = A_{uu}^{-1} (b_u - A_{up} * x_p) \f$. + Same as Solve(MatrixPartitioned&, const xt::xtensor&, const xt::xtensor&), + but with partitioned input and output. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b_u unknown dofval [nnu]. \param x_p prescribed dofval [nnp] \return x_u unknown dofval [nnu]. */ xt::xtensor Solve_u( MatrixPartitioned& A, const xt::xtensor& b_u, const xt::xtensor& x_p); /** Same as Solve_u(MatrixPartitioned&, const xt::xtensor&, const xt::xtensor&), but writing to pre-allocated output. - \param A sparse matrix, see Matrix(). + \param A sparse matrix, see MatrixPartitioned(). \param b_u unknown dofval [nnu]. \param x_p prescribed dofval [nnp] \param x_u (overwritten) unknown dofval [nnu]. */ void solve_u( MatrixPartitioned& A, const xt::xtensor& b_u, const xt::xtensor& x_p, xt::xtensor& x_u); private: Solver m_solver; ///< solver bool m_factor = true; ///< signal to force factorization void factorize(MatrixPartitioned& matrix); ///< compute inverse (evaluated by "solve") }; } // namespace GooseFEM #include "MatrixPartitioned.hpp" #endif diff --git a/include/GooseFEM/MatrixPartitionedTyings.h b/include/GooseFEM/MatrixPartitionedTyings.h index 3d5fb10..b282f48 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.h +++ b/include/GooseFEM/MatrixPartitionedTyings.h @@ -1,175 +1,286 @@ /** Sparse matrix that is partitioned in: - unknown DOFs - prescribed DOFs - tied DOFs \file MatrixPartitionedTyings.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H #include "config.h" #include "Matrix.h" #include #include #include namespace GooseFEM { // forward declaration template class MatrixPartitionedTyingsSolver; +/** +Sparse matrix from with dependent DOFs are eliminated, +and the remaining (small) independent system is partitioned in an unknown and a prescribed part. +In particular: + +\f$ A_{ii} = \begin{bmatrix} A_{uu} & A_{up} \\ A_{pu} & A_{pp} \end{bmatrix} \f$ + +See VectorPartitionedTyings() for bookkeeping definitions. +*/ class MatrixPartitionedTyings : public Matrix { public: MatrixPartitionedTyings() = default; + /** + Constructor. + + \param conn connectivity [#nelem, #nne]. + \param dofs DOFs per node [#nnode, #ndim]. + \param Cdu See Tyings::Periodic::Cdu(). + \param Cdp See Tyings::Periodic::Cdp(). + */ MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp); - // Dimensions - size_t nnu() const; // number of independent, unknown DOFs - size_t nnp() const; // number of independent, prescribed DOFs - size_t nni() const; // number of independent DOFs - size_t nnd() const; // number of dependent DOFs + /** + \return Number of dependent DOFs. + */ + size_t nnd() const; + + /** + \return Number of independent DOFs. + */ + size_t nni() const; + + /** + \return Number of independent unknown DOFs. + */ + size_t nnu() const; + + /** + \return Number of independent prescribed DOFs. + */ + size_t nnp() const; + + /** + Dependent DOFs. + + \return List of DOF numbers. + */ + xt::xtensor iid() const; + + /** + Independent DOFs. + + \return List of DOF numbers. + */ + xt::xtensor iii() const; - // DOF lists - xt::xtensor iiu() const; // independent, unknown DOFs - xt::xtensor iip() const; // independent, prescribed DOFs - xt::xtensor iii() const; // independent DOFs - xt::xtensor iid() const; // dependent DOFs + /** + Independent unknown DOFs. + + \return List of DOF numbers. + */ + xt::xtensor iiu() const; + + /** + Independent prescribed DOFs. + + \return List of DOF numbers. + */ + xt::xtensor iip() const; - // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] void assemble(const xt::xtensor& elemmat) override; private: + using Matrix::set; using Matrix::add; using Matrix::Todense; using Matrix::todense; using Matrix::Dot; using Matrix::dot; private: - // The matrix - Eigen::SparseMatrix m_Auu; - Eigen::SparseMatrix m_Aup; - Eigen::SparseMatrix m_Apu; - Eigen::SparseMatrix m_App; - Eigen::SparseMatrix m_Aud; - Eigen::SparseMatrix m_Apd; - Eigen::SparseMatrix m_Adu; - Eigen::SparseMatrix m_Adp; - Eigen::SparseMatrix m_Add; - - // The matrix for which the tyings have been applied - Eigen::SparseMatrix m_ACuu; - Eigen::SparseMatrix m_ACup; - Eigen::SparseMatrix m_ACpu; - Eigen::SparseMatrix m_ACpp; - - // Matrix entries - std::vector> m_Tuu; - std::vector> m_Tup; - std::vector> m_Tpu; - std::vector> m_Tpp; - std::vector> m_Tud; - std::vector> m_Tpd; - std::vector> m_Tdu; - std::vector> m_Tdp; - std::vector> m_Tdd; - - // Bookkeeping - xt::xtensor m_iiu; // unknown DOFs [nnu] - xt::xtensor m_iip; // prescribed DOFs [nnp] - xt::xtensor m_iid; // dependent DOFs [nnd] - - // Dimensions - size_t m_nnu; // number of independent, unknown DOFs - size_t m_nnp; // number of independent, prescribed DOFs - size_t m_nni; // number of independent DOFs - size_t m_nnd; // number of dependent DOFs - - // Tyings - Eigen::SparseMatrix m_Cdu; - Eigen::SparseMatrix m_Cdp; - Eigen::SparseMatrix m_Cud; - Eigen::SparseMatrix m_Cpd; + + Eigen::SparseMatrix m_Auu; ///< The matrix. + Eigen::SparseMatrix m_Aup; ///< The matrix. + Eigen::SparseMatrix m_Apu; ///< The matrix. + Eigen::SparseMatrix m_App; ///< The matrix. + Eigen::SparseMatrix m_Aud; ///< The matrix. + Eigen::SparseMatrix m_Apd; ///< The matrix. + Eigen::SparseMatrix m_Adu; ///< The matrix. + Eigen::SparseMatrix m_Adp; ///< The matrix. + Eigen::SparseMatrix m_Add; ///< The matrix. + Eigen::SparseMatrix m_ACuu; ///< // The matrix for which the tyings have been applied. + Eigen::SparseMatrix m_ACup; ///< // The matrix for which the tyings have been applied. + Eigen::SparseMatrix m_ACpu; ///< // The matrix for which the tyings have been applied. + Eigen::SparseMatrix m_ACpp; ///< // The matrix for which the tyings have been applied. + std::vector> m_Tuu; ///< Matrix entries. + std::vector> m_Tup; ///< Matrix entries. + std::vector> m_Tpu; ///< Matrix entries. + std::vector> m_Tpp; ///< Matrix entries. + std::vector> m_Tud; ///< Matrix entries. + std::vector> m_Tpd; ///< Matrix entries. + std::vector> m_Tdu; ///< Matrix entries. + std::vector> m_Tdp; ///< Matrix entries. + std::vector> m_Tdd; ///< Matrix entries. + xt::xtensor m_iiu; ///< See iiu() + xt::xtensor m_iip; ///< See iip() + xt::xtensor m_iid; ///< See iid() + size_t m_nnu; ///< See #nnu + size_t m_nnp; ///< See #nnp + size_t m_nni; ///< See #nni + size_t m_nnd; ///< See #nnd + Eigen::SparseMatrix m_Cdu; ///< Tying matrix, see Tyings::Periodic::Cdu(). + Eigen::SparseMatrix m_Cdp; ///< Tying matrix, see Tyings::Periodic::Cdp(). + Eigen::SparseMatrix m_Cud; ///< Transpose of "m_Cdu". + Eigen::SparseMatrix m_Cpd; ///< Transpose of "m_Cdp". // grant access to solver class template friend class MatrixPartitionedTyingsSolver; +private: + // Convert arrays (Eigen version of VectorPartitioned, which contains public functions) Eigen::VectorXd AsDofs_u(const xt::xtensor& dofval) const; Eigen::VectorXd AsDofs_u(const xt::xtensor& nodevec) const; Eigen::VectorXd AsDofs_p(const xt::xtensor& dofval) const; Eigen::VectorXd AsDofs_p(const xt::xtensor& nodevec) const; Eigen::VectorXd AsDofs_d(const xt::xtensor& dofval) const; Eigen::VectorXd AsDofs_d(const xt::xtensor& nodevec) const; }; +/** +Solver for MatrixPartitionedTyings(). +The idea is that this solver class can be used to solve for multiple right-hand-sides +using one factorisation. +*/ template >> class MatrixPartitionedTyingsSolver { public: - // Constructors + MatrixPartitionedTyingsSolver() = default; - // Solve: - // A' = A_ii + K_id * C_di + C_di^T * K_di + C_di^T * K_dd * C_di - // b' = b_i + C_di^T * b_d - // x_u = A'_uu \ ( b'_u - A'_up * x_p ) - // x_i = [x_u, x_p] - // x_d = C_di * x_i - void solve( - MatrixPartitionedTyings& matrix, - const xt::xtensor& b, - xt::xtensor& x); // updates x_u and x_d + /** + Solve as follows. - void solve( - MatrixPartitionedTyings& matrix, - const xt::xtensor& b, - xt::xtensor& x); // updates x_u and x_d + \f$ A' = A_{ii} + A_{id} * C_{di} + C_{di}^T * A_{di} + C_{di}^T * A_{dd} * C_{di} \f$ - void solve_u( - MatrixPartitionedTyings& matrix, - const xt::xtensor& b_u, - const xt::xtensor& b_d, - const xt::xtensor& x_p, - xt::xtensor& x_u); + \f$ b' = b_i + C_{di}^T * b_d \f$ + + \f$ x_u = A'_{uu} \ ( b'_u - A'_{up} * x_p ) \f$ - // Auto-allocation of the functions above + \f$ x_i = \begin{bmatrix} x_u \\ x_p \end{bmatrix} \f$ + + \f$ x_d = C_{di} * x_i \f$ + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b nodevec [nelem, ndim]. + \param x nodevec [nelem, ndim], used to read \f$ x_p \f$. + \return x nodevec [nelem, ndim], \f$ x_u \f$ filled, \f$ x_p \f$ copied. + */ xt::xtensor Solve( - MatrixPartitionedTyings& matrix, + MatrixPartitionedTyings& A, const xt::xtensor& b, const xt::xtensor& x); + /** + Same as + Solve(MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&), + but filling \f$ x_u \f$ and \f$ x_d \f$ in place. + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b nodevec [nelem, ndim]. + \param x nodevec [nelem, ndim], \f$ x_p \f$ read, \f$ x_u \f$ and \f$ x_d \f$ filled. + */ + void solve( + MatrixPartitionedTyings& A, + const xt::xtensor& b, + xt::xtensor& x); + + /** + Same as + Solve(MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&), + but for "dofval" input and output. + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b dofval [ndof]. + \param x dofval [ndof], used to read \f$ x_p \f$. + \return x dofval [ndof], \f$ x_u \f$ and \f$ x_d \f$ filled, \f$ x_p \f$ copied. + */ xt::xtensor Solve( - MatrixPartitionedTyings& matrix, + MatrixPartitionedTyings& A, const xt::xtensor& b, const xt::xtensor& x); + /** + Same as + Solve(MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&), + but filling \f$ x_u \f$ and \f$ x_d \f$ in place. + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b dofval [ndof]. + \param x dofval [ndof], \f$ x_p \f$ read, \f$ x_u \f$ and \f$ x_d \f$ filled. + */ + void solve( + MatrixPartitionedTyings& A, + const xt::xtensor& b, + xt::xtensor& x); + + /** + Same as + Solve(MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&), + but with partitioned input and output. + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b_u unknown dofval [nnu]. + \param b_d dependent dofval [nnd]. + \param x_p prescribed dofval [nnp] + \return x_u unknown dofval [nnu]. + */ xt::xtensor Solve_u( - MatrixPartitionedTyings& matrix, + MatrixPartitionedTyings& A, const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p); + /** + Same as + Solve_u(MatrixPartitionedTyings&, const xt::xtensor&, const xt::xtensor&, const xt::xtensor&), + but writing to pre-allocated output. + + \param A sparse matrix, see MatrixPartitionedTyings(). + \param b_u unknown dofval [nnu]. + \param b_d dependent dofval [nnd]. + \param x_p prescribed dofval [nnp] + \param x_u (overwritten) unknown dofval [nnu]. + */ + void solve_u( + MatrixPartitionedTyings& A, + const xt::xtensor& b_u, + const xt::xtensor& b_d, + const xt::xtensor& x_p, + xt::xtensor& x_u); + private: - Solver m_solver; // solver - bool m_factor = true; // signal to force factorization - void factorize(MatrixPartitionedTyings& matrix); // compute inverse (evaluated by "solve") + Solver m_solver; ///< solver + bool m_factor = true; ///< signal to force factorization + void factorize(MatrixPartitionedTyings& matrix); ///< compute inverse (evaluated by "solve") }; } // namespace GooseFEM #include "MatrixPartitionedTyings.hpp" #endif