diff --git a/include/GooseFEM/Matrix.h b/include/GooseFEM/Matrix.h index a1cd6b5..12831ee 100644 --- a/include/GooseFEM/Matrix.h +++ b/include/GooseFEM/Matrix.h @@ -1,239 +1,239 @@ /** Sparse matrix. \file Matrix.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIX_H #define GOOSEFEM_MATRIX_H #include "config.h" #include #include #include namespace GooseFEM { // forward declaration template class MatrixSolver; /** Sparse matrix. See Vector() for bookkeeping definitions. */ class Matrix { public: Matrix() = default; /** Constructor. \param conn connectivity [#nelem, #nne]. \param dofs DOFs per node [#nnode, #ndim]. */ Matrix(const xt::xtensor& conn, const xt::xtensor& dofs); /** \return Number of elements. */ size_t nelem() const; /** \return Number of nodes per element. */ size_t nne() const; /** \return Number of nodes. */ size_t nnode() const; /** \return Number of dimensions. */ size_t ndim() const; /** \return Number of DOFs. */ size_t ndof() const; /** \return DOFs per node [#nnode, #ndim] */ xt::xtensor dofs() const; /** Assemble from matrices stored per element. \param elemmat [#nelem, #nne * #ndim, #nne * #ndim]. */ virtual void assemble(const xt::xtensor& elemmat); /** Overwrite matrix with dense matrix. \param rows Row numbers in the matrix [n]. \param cols Column numbers in the matrix [n]. \param matrix Data entries on (rows, cols) [n]. */ virtual void set( const xt::xtensor& rows, const xt::xtensor& cols, const xt::xtensor& matrix); /** Add a dense matrix to the matrix. \param rows Row numbers in the matrix [n]. \param cols Column numbers in the matrix [n]. \param matrix Data entries on (rows, cols) [n]. */ virtual void add( const xt::xtensor& rows, const xt::xtensor& cols, const xt::xtensor& matrix); /** \return Matrix as dense matrix [#ndof, #ndof]. */ xt::xtensor Todense() const; /** Covert matrix to dense matrix. \param ret overwritten [#ndof, #ndof]. */ virtual void todense(xt::xtensor& ret) const; /** Dot-product \f$ b_i = A_{ij} x_j \f$. \param x nodevec [#nelem, #ndim]. \return b nodevec overwritten [#nelem, #ndim]. */ xt::xtensor Dot(const xt::xtensor& x) const; /** Same as Dot(const xt::xtensor&, xt::xtensor& b) but writing to preallocated data. \param x nodevec [#nelem, #ndim]. \param b nodevec overwritten [#nelem, #ndim]. */ virtual void dot(const xt::xtensor& x, xt::xtensor& b) const; /** Dot-product \f$ b_i = A_{ij} x_j \f$. \param x dofval [#ndof]. \return b dofval overwritten [#ndof]. */ xt::xtensor Dot(const xt::xtensor& x) const; /** Same as Dot(const xt::xtensor&, xt::xtensor& b) but writing to preallocated data. \param x dofval [#ndof]. \param b dofval overwritten [#ndof]. */ virtual void dot(const xt::xtensor& x, xt::xtensor& b) const; protected: xt::xtensor m_conn; ///< Connectivity [#nelem, #nne]. xt::xtensor m_dofs; ///< DOF-numbers per node [#nnode, #ndim]. size_t m_nelem; ///< See nelem(). size_t m_nne; ///< See nne(). size_t m_nnode; ///< See nnode(). size_t m_ndim; ///< See ndim(). size_t m_ndof; ///< See ndof(). bool m_changed = true; ///< Signal changes to data. private: Eigen::SparseMatrix m_A; ///< The matrix. std::vector> m_T; ///< Matrix entries. template friend class MatrixSolver; ///< Grant access to solver class /** Convert "nodevec" to "dofval" (overwrite entries that occur more than once). \param nodevec input [#nnode, #ndim] \return dofval output [#ndof] */ Eigen::VectorXd AsDofs(const xt::xtensor& nodevec) const; /** Convert "dofval" to "nodevec" (overwrite entries that occur more than once). \param dofval input [#ndof] \param nodevec output [#nnode, #ndim] */ void asNode(const Eigen::VectorXd& dofval, xt::xtensor& nodevec) const; }; /** Solver for Matrix(). The idea is that this solver class can be used to solve for multiple right-hand-sides using one factorisation. */ template >> class MatrixSolver { public: MatrixSolver() = default; /** - Solve \f$ x_u = A_{uu}^{-1} (b_u - A_{up} * x_p) \f$. + Solve \f$ x = A^{-1} b \f$. \param A sparse matrix, see Matrix(). \param b nodevec [nelem, ndim]. \return x nodevec [nelem, ndim]. */ xt::xtensor Solve(Matrix& A, const xt::xtensor& b); /** Same as Solve(Matrix&, const xt::xtensor&) but writing to preallocated data. \param A sparse matrix, see Matrix(). \param b nodevec [nelem, ndim]. \param x nodevec overwritten [nelem, ndim]. */ void solve(Matrix& A, const xt::xtensor& b, xt::xtensor& x); /** - Solve \f$ x_u = A_{uu}^{-1} (b_u - A_{up} * x_p) \f$. + Solve \f$ x = A^{-1} b \f$. \param A sparse matrix, see Matrix(). \param b dofval [ndof]. \return x dofval [ndof]. */ xt::xtensor Solve(Matrix& A, const xt::xtensor& b); /** Same as Solve(Matrix&, const xt::xtensor&) but writing to preallocated data. \param A sparse matrix, see Matrix(). \param b dofval [ndof]. \param x dofval overwritten [ndof]. */ void solve(Matrix& A, const xt::xtensor& b, xt::xtensor& x); private: Solver m_solver; ///< Solver. bool m_factor = true; ///< Signal to force factorization. void factorize(Matrix& A); ///< Compute inverse (evaluated by "solve"). }; } // namespace GooseFEM #include "Matrix.hpp" #endif diff --git a/include/GooseFEM/MatrixPartitioned.h b/include/GooseFEM/MatrixPartitioned.h index 46ea09d..7baf514 100644 --- a/include/GooseFEM/MatrixPartitioned.h +++ b/include/GooseFEM/MatrixPartitioned.h @@ -1,194 +1,294 @@ /** 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: - // b_p = A_pu * x_u + A_pp * x_p = A_pp * x_p + /** + 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; // modified with "b_p" + 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; // modified with "b_p" + 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; + xt::xtensor& b_p) const; - // Auto-allocation of the functions above +private: - xt::xtensor Reaction( - const xt::xtensor& x, const xt::xtensor& b) const; - xt::xtensor Reaction( - const xt::xtensor& x, const xt::xtensor& b) const; - xt::xtensor Reaction_p( - const xt::xtensor& x_u, const xt::xtensor& x_p) const; + 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 -private: - // The matrix - Eigen::SparseMatrix m_Auu; - Eigen::SparseMatrix m_Aup; - Eigen::SparseMatrix m_Apu; - Eigen::SparseMatrix m_App; - - // Matrix entries - std::vector> m_Tuu; - std::vector> m_Tup; - std::vector> m_Tpu; - std::vector> m_Tpp; - - // Signal changes to data compare to the last inverse - bool m_changed = true; - - // Bookkeeping - xt::xtensor m_part; // DOF-numbers per node, renumbered [nnode, ndim] - xt::xtensor m_iiu; // unknown DOFs [nnu] - xt::xtensor m_iip; // prescribed DOFs [nnp] - - // Dimensions - size_t m_nnu; // number of unknown DOFs - size_t m_nnp; // number of prescribed DOFs + /** + 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: - // Constructors - MatrixPartitionedSolver() = default; - // Solve: - // x_u = A_uu \ ( b_u - A_up * x_p ) - void solve( - MatrixPartitioned& matrix, - const xt::xtensor& b, - xt::xtensor& x); // modified with "x_u" - - void solve( - MatrixPartitioned& matrix, - const xt::xtensor& b, - xt::xtensor& x); // modified with "x_u" + MatrixPartitionedSolver() = default; - void solve_u( - MatrixPartitioned& matrix, - const xt::xtensor& b_u, - const xt::xtensor& x_p, - xt::xtensor& x_u); + /** + Solve \f$ x_u = A_{uu}^{-1} (b_u - A_{up} * x_p) \f$. - // Auto-allocation of the functions above + \param A sparse matrix, see Matrix(). + \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& matrix, + 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 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 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& matrix, + 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 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$. + + \param A sparse matrix, see Matrix(). + \param b_u unknown dofval [nnu]. + \param x_p prescribed dofval [nnp] + \return x_u unknown dofval [nnu]. + */ xt::xtensor Solve_u( - MatrixPartitioned& matrix, + 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 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") + 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 3e4a6cb..3d5fb10 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.h +++ b/include/GooseFEM/MatrixPartitionedTyings.h @@ -1,182 +1,175 @@ /** 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; -class MatrixPartitionedTyings { +class MatrixPartitionedTyings : public Matrix { public: - // Constructors + MatrixPartitionedTyings() = default; MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp); // Dimensions - size_t nelem() const; // number of elements - size_t nne() const; // number of nodes per element - size_t nnode() const; // number of nodes - size_t ndim() const; // number of dimensions - size_t ndof() const; // number of DOFs 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 // DOF lists - xt::xtensor dofs() const; // DOFs 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 // Assemble from matrices stored per element [nelem, nne*ndim, nne*ndim] - void assemble(const xt::xtensor& elemmat); + 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; - // Signal changes to data - bool m_changed = true; - // Bookkeeping - xt::xtensor m_conn; // connectivity [nelem, nne ] - xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] 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_nelem; // number of elements - size_t m_nne; // number of nodes per element - size_t m_nnode; // number of nodes - size_t m_ndim; // number of dimensions - size_t m_ndof; // number of DOFs 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; // grant access to solver class template friend class MatrixPartitionedTyingsSolver; // 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; }; 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 void solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, xt::xtensor& x); // updates x_u and x_d void solve_u( MatrixPartitionedTyings& matrix, const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p, xt::xtensor& x_u); // Auto-allocation of the functions above xt::xtensor Solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, const xt::xtensor& x); xt::xtensor Solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, const xt::xtensor& x); xt::xtensor Solve_u( MatrixPartitionedTyings& matrix, const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p); private: 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 diff --git a/include/GooseFEM/MatrixPartitionedTyings.hpp b/include/GooseFEM/MatrixPartitionedTyings.hpp index 49fb8e4..a285f95 100644 --- a/include/GooseFEM/MatrixPartitionedTyings.hpp +++ b/include/GooseFEM/MatrixPartitionedTyings.hpp @@ -1,457 +1,430 @@ /** Implementation of MatrixPartitionedTyings.h \file MatrixPartitionedTyings.hpp \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP #define GOOSEFEM_MATRIXPARTITIONEDTYINGS_HPP #include "MatrixPartitionedTyings.h" namespace GooseFEM { inline MatrixPartitionedTyings::MatrixPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp) - : m_conn(conn), m_dofs(dofs), m_Cdu(Cdu), m_Cdp(Cdp) { GOOSEFEM_ASSERT(Cdu.rows() == Cdp.rows()); + m_conn = conn; + m_dofs = dofs; + m_Cdu = Cdu; + m_Cdp = Cdp; m_nnu = static_cast(m_Cdu.cols()); m_nnp = static_cast(m_Cdp.cols()); m_nnd = static_cast(m_Cdp.rows()); m_nni = m_nnu + m_nnp; m_ndof = m_nni + m_nnd; m_iiu = xt::arange(m_nnu); m_iip = xt::arange(m_nnu, m_nnu + m_nnp); m_iid = xt::arange(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); } -inline size_t MatrixPartitionedTyings::nelem() const -{ - return m_nelem; -} - -inline size_t MatrixPartitionedTyings::nne() const -{ - return m_nne; -} - -inline size_t MatrixPartitionedTyings::nnode() const -{ - return m_nnode; -} - -inline size_t MatrixPartitionedTyings::ndim() const -{ - return m_ndim; -} - -inline size_t MatrixPartitionedTyings::ndof() const -{ - return m_ndof; -} - inline size_t MatrixPartitionedTyings::nnu() const { return m_nnu; } inline size_t MatrixPartitionedTyings::nnp() const { return m_nnp; } inline size_t MatrixPartitionedTyings::nni() const { return m_nni; } inline size_t MatrixPartitionedTyings::nnd() const { return m_nnd; } -inline xt::xtensor MatrixPartitionedTyings::dofs() const -{ - return m_dofs; -} - inline xt::xtensor MatrixPartitionedTyings::iiu() const { return m_iiu; } inline xt::xtensor MatrixPartitionedTyings::iip() const { return m_iip; } inline xt::xtensor MatrixPartitionedTyings::iii() const { return xt::arange(m_nni); } inline xt::xtensor MatrixPartitionedTyings::iid() const { return m_iid; } inline void MatrixPartitionedTyings::assemble(const xt::xtensor& 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( 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( 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( 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( 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( 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( 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( 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( di - m_nni, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j))); } else { m_Tdd.push_back(Eigen::Triplet( 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_u(const xt::xtensor& 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_u(const xt::xtensor& 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_p(const xt::xtensor& 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_p(const xt::xtensor& 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_d(const xt::xtensor& 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; } inline Eigen::VectorXd MatrixPartitionedTyings::AsDofs_d(const xt::xtensor& 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 inline void MatrixPartitionedTyingsSolver::factorize(MatrixPartitionedTyings& matrix) { if (!matrix.m_changed && !m_factor) { return; } matrix.m_ACuu = matrix.m_Auu + matrix.m_Aud * matrix.m_Cdu + matrix.m_Cud * matrix.m_Adu + matrix.m_Cud * matrix.m_Add * matrix.m_Cdu; matrix.m_ACup = matrix.m_Aup + matrix.m_Aud * matrix.m_Cdp + matrix.m_Cud * matrix.m_Adp + matrix.m_Cud * matrix.m_Add * matrix.m_Cdp; // matrix.m_ACpu = matrix.m_Apu + matrix.m_Apd * matrix.m_Cdu + matrix.m_Cpd * matrix.m_Adu // + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdu; // matrix.m_ACpp = matrix.m_App + matrix.m_Apd * matrix.m_Cdp + matrix.m_Cpd * matrix.m_Adp // + matrix.m_Cpd * matrix.m_Add * matrix.m_Cdp; m_solver.compute(matrix.m_ACuu); m_factor = false; matrix.m_changed = false; } template inline void MatrixPartitionedTyingsSolver::solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(xt::has_shape(b, {matrix.m_nnode, matrix.m_ndim})); GOOSEFEM_ASSERT(xt::has_shape(x, {matrix.m_nnode, matrix.m_ndim})); this->factorize(matrix); Eigen::VectorXd B_u = matrix.AsDofs_u(b); Eigen::VectorXd B_d = matrix.AsDofs_d(b); Eigen::VectorXd X_p = matrix.AsDofs_p(x); B_u += matrix.m_Cud * B_d; Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p)); Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p; #pragma omp parallel for for (size_t m = 0; m < matrix.m_nnode; ++m) { for (size_t i = 0; i < matrix.m_ndim; ++i) { if (matrix.m_dofs(m, i) < matrix.m_nnu) { x(m, i) = X_u(matrix.m_dofs(m, i)); } else if (matrix.m_dofs(m, i) >= matrix.m_nni) { x(m, i) = X_d(matrix.m_dofs(m, i) - matrix.m_nni); } } } } template inline void MatrixPartitionedTyingsSolver::solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, xt::xtensor& x) { GOOSEFEM_ASSERT(b.size() == matrix.m_ndof); GOOSEFEM_ASSERT(x.size() == matrix.m_ndof); this->factorize(matrix); Eigen::VectorXd B_u = matrix.AsDofs_u(b); Eigen::VectorXd B_d = matrix.AsDofs_d(b); Eigen::VectorXd X_p = matrix.AsDofs_p(x); Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - matrix.m_ACup * X_p)); Eigen::VectorXd X_d = matrix.m_Cdu * X_u + matrix.m_Cdp * X_p; #pragma omp parallel for for (size_t d = 0; d < matrix.m_nnu; ++d) { x(matrix.m_iiu(d)) = X_u(d); } #pragma omp parallel for for (size_t d = 0; d < matrix.m_nnd; ++d) { x(matrix.m_iid(d)) = X_d(d); } } template inline void MatrixPartitionedTyingsSolver::solve_u( MatrixPartitionedTyings& matrix, const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p, xt::xtensor& x_u) { UNUSED(b_d); GOOSEFEM_ASSERT(b_u.size() == matrix.m_nnu); GOOSEFEM_ASSERT(b_d.size() == matrix.m_nnd); GOOSEFEM_ASSERT(x_p.size() == matrix.m_nnp); GOOSEFEM_ASSERT(x_u.size() == matrix.m_nnu); this->factorize(matrix); Eigen::Map(x_u.data(), x_u.size()).noalias() = m_solver.solve(Eigen::VectorXd( Eigen::Map(b_u.data(), b_u.size()) - matrix.m_ACup * Eigen::Map(x_p.data(), x_p.size()))); } template inline xt::xtensor MatrixPartitionedTyingsSolver::Solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, const xt::xtensor& x) { xt::xtensor ret = x; this->solve(matrix, b, ret); return ret; } template inline xt::xtensor MatrixPartitionedTyingsSolver::Solve( MatrixPartitionedTyings& matrix, const xt::xtensor& b, const xt::xtensor& x) { xt::xtensor ret = x; this->solve(matrix, b, ret); return ret; } template inline xt::xtensor MatrixPartitionedTyingsSolver::Solve_u( MatrixPartitionedTyings& matrix, const xt::xtensor& b_u, const xt::xtensor& b_d, const xt::xtensor& x_p) { xt::xtensor x_u = xt::empty({matrix.m_nnu}); this->solve_u(matrix, b_u, b_d, x_p, x_u); return x_u; } } // namespace GooseFEM #endif diff --git a/include/GooseFEM/VectorPartitionedTyings.h b/include/GooseFEM/VectorPartitionedTyings.h index 4a485a4..4156462 100644 --- a/include/GooseFEM/VectorPartitionedTyings.h +++ b/include/GooseFEM/VectorPartitionedTyings.h @@ -1,159 +1,160 @@ /** Methods to switch between storage types based on a mesh and DOFs that are partitioned in: - unknown DOFs - prescribed DOFs - dependent DOFs \file VectorPartitionedTyings.h \copyright Copyright 2017. Tom de Geus. All rights reserved. \license This project is released under the GNU Public License (GPLv3). */ #ifndef GOOSEFEM_VECTORPARTITIONEDTYINGS_H #define GOOSEFEM_VECTORPARTITIONEDTYINGS_H #include "config.h" +#include "Vector.h" #include #include namespace GooseFEM { /** Class to switch between storage types. In particular: - "nodevec": nodal vectors [#nnode, #ndim]. - "elemvec": nodal vectors stored per element [nelem, #nne, #ndim]. - "dofval": DOF values [#ndof]. - "dofval_u": DOF values (Unknown), `== dofval[iiu]`, [#nnu]. - "dofval_p": DOF values (Prescribed), `== dofval[iiu]`, [#nnp]. */ class VectorPartitionedTyings : public Vector { public: VectorPartitionedTyings() = 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(). \param Cdi See Tyings::Periodic::Cdi(). */ VectorPartitionedTyings( const xt::xtensor& conn, const xt::xtensor& dofs, const Eigen::SparseMatrix& Cdu, const Eigen::SparseMatrix& Cdp, const Eigen::SparseMatrix& Cdi); /** \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; /** Independent unknown DOFs. \return List of DOF numbers. */ xt::xtensor iiu() const; /** Independent prescribed DOFs. \return List of DOF numbers. */ xt::xtensor iip() const; /** Copy (part of) "dofval" to another "dofval": dofval_dest[iip()] = dofval_src[iip()]. \param dofval_src DOF values, iip() updated, [#ndof]. \param dofval_dest DOF values, iip() updated, [#ndof]. */ void copy_p( const xt::xtensor& dofval_src, xt::xtensor& dofval_dest) const; /** Convert to "dofval" (overwrite entries that occur more than once). Only the independent DOFs are retained. \param nodevec nodal vectors [#nnode, #ndim]. \return dofval[iii()] [#nni]. */ xt::xtensor AsDofs_i(const xt::xtensor& nodevec) const; /** Same as InterpQuad_vector(), but writing to preallocated return. \param nodevec [#nnode, #ndim]. \param dofval_i [#nni]. \param apply_tyings If `true` the dependent DOFs are eliminated. */ void asDofs_i( const xt::xtensor& nodevec, xt::xtensor& dofval_i, bool apply_tyings = true) const; private: 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_Cdi; ///< Tying matrix, see Tyings::Periodic::Cdi(). Eigen::SparseMatrix m_Cud; ///< Transpose of "m_Cdu". Eigen::SparseMatrix m_Cpd; ///< Transpose of "m_Cdp". Eigen::SparseMatrix m_Cid; ///< Transpose of "m_Cdi". /** Convert to "dofval" (overwrite entries that occur more than once). Only the dependent DOFs are retained. \param nodevec nodal vectors [#nnode, #ndim]. \return dofval[iid()] [#nnd]. */ Eigen::VectorXd Eigen_asDofs_d(const xt::xtensor& nodevec) const; }; } // namespace GooseFEM #include "VectorPartitionedTyings.hpp" #endif