Page MenuHomec4science
Files F61312101

MatrixPartitioned.h
No OneTemporary
Actions

File Metadata

Created
Sun, May 5, 21:14

MatrixPartitioned.h
View Options

/**
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 <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
namespace GooseFEM {
// forward declaration
template <class> 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<size_t, 2>& conn,
const xt::xtensor<size_t, 2>& dofs,
const xt::xtensor<size_t, 1>& iip);
/**
Number of unknown DOFs.
*/
size_t nnu() const;
/**
Number of prescribed DOFs.
*/
size_t nnp() const;
/**
Unknown DOFs [#nnu].
*/
xt::xtensor<size_t, 1> iiu() const;
/**
Prescribed DOFs [#nnp].
*/
xt::xtensor<size_t, 1> iip() const;
void assemble(const xt::xtensor<double, 3>& elemmat) override;
void set(
const xt::xtensor<size_t, 1>& rows,
const xt::xtensor<size_t, 1>& cols,
const xt::xtensor<double, 2>& matrix) override;
void add(
const xt::xtensor<size_t, 1>& rows,
const xt::xtensor<size_t, 1>& cols,
const xt::xtensor<double, 2>& matrix) override;
void todense(xt::xtensor<double, 2>& ret) const override;
void dot(const xt::xtensor<double, 2>& x, xt::xtensor<double, 2>& b) const override;
void dot(const xt::xtensor<double, 1>& x, xt::xtensor<double, 1>& 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
void reaction(
const xt::xtensor<double, 2>& x,
xt::xtensor<double, 2>& b) const; // modified with "b_p"
void reaction(
const xt::xtensor<double, 1>& x,
xt::xtensor<double, 1>& b) const; // modified with "b_p"
void reaction_p(
const xt::xtensor<double, 1>& x_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& b_p) const;
// Auto-allocation of the functions above
xt::xtensor<double, 2> Reaction(
const xt::xtensor<double, 2>& x, const xt::xtensor<double, 2>& b) const;
xt::xtensor<double, 1> Reaction(
const xt::xtensor<double, 1>& x, const xt::xtensor<double, 1>& b) const;
xt::xtensor<double, 1> Reaction_p(
const xt::xtensor<double, 1>& x_u, const xt::xtensor<double, 1>& x_p) const;
private:
// The matrix
Eigen::SparseMatrix<double> m_Auu;
Eigen::SparseMatrix<double> m_Aup;
Eigen::SparseMatrix<double> m_Apu;
Eigen::SparseMatrix<double> m_App;
// Matrix entries
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;
// Signal changes to data compare to the last inverse
bool m_changed = true;
// Bookkeeping
xt::xtensor<size_t, 2> m_part; // DOF-numbers per node, renumbered [nnode, ndim]
xt::xtensor<size_t, 1> m_iiu; // unknown DOFs [nnu]
xt::xtensor<size_t, 1> m_iip; // prescribed DOFs [nnp]
// Dimensions
size_t m_nnu; // number of unknown DOFs
size_t m_nnp; // number of prescribed DOFs
// grant access to solver class
template <class> friend class MatrixPartitionedSolver;
// Convert arrays (Eigen version of VectorPartitioned, which contains public functions)
Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 1>& dofval) const;
Eigen::VectorXd AsDofs_u(const xt::xtensor<double, 2>& nodevec) const;
Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 1>& dofval) const;
Eigen::VectorXd AsDofs_p(const xt::xtensor<double, 2>& nodevec) const;
};
template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
class MatrixPartitionedSolver {
public:
// Constructors
MatrixPartitionedSolver() = default;
// Solve:
// x_u = A_uu \ ( b_u - A_up * x_p )
void solve(
MatrixPartitioned& matrix,
const xt::xtensor<double, 2>& b,
xt::xtensor<double, 2>& x); // modified with "x_u"
void solve(
MatrixPartitioned& matrix,
const xt::xtensor<double, 1>& b,
xt::xtensor<double, 1>& x); // modified with "x_u"
void solve_u(
MatrixPartitioned& matrix,
const xt::xtensor<double, 1>& b_u,
const xt::xtensor<double, 1>& x_p,
xt::xtensor<double, 1>& x_u);
// Auto-allocation of the functions above
xt::xtensor<double, 2> Solve(
MatrixPartitioned& matrix,
const xt::xtensor<double, 2>& b,
const xt::xtensor<double, 2>& x);
xt::xtensor<double, 1> Solve(
MatrixPartitioned& matrix,
const xt::xtensor<double, 1>& b,
const xt::xtensor<double, 1>& x);
xt::xtensor<double, 1> Solve_u(
MatrixPartitioned& matrix,
const xt::xtensor<double, 1>& b_u,
const xt::xtensor<double, 1>& x_p);
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

Event Timeline

c4science · Help