Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77493665
Matrix.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Aug 14, 19:41
Size
23 KB
Mime Type
text/x-c++
Expires
Fri, Aug 16, 19:41 (2 d)
Engine
blob
Format
Raw Data
Handle
19878806
Attached To
rGOOSEFEM GooseFEM
Matrix.h
View Options
/**
* @file
* @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 <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
namespace GooseFEM {
// forward declaration
template <class>
class MatrixSolver;
/**
* CRTP base class for a solver class.
*/
template <class D>
class MatrixSolverBase {
public:
/**
* Underlying type.
*/
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:
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b nodevec [nelem, ndim].
* @param x (overwritten) nodevec [nelem, ndim].
*/
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);
}
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b dofval [ndof].
* @param x (overwritten) dofval [ndof].
*/
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);
}
};
/**
* CRTP base class for a solver class.
*/
template <class D>
class MatrixSolverSingleBase {
public:
/**
* Underlying type.
*/
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:
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b nodevec [nelem, ndim].
* @return x nodevec [nelem, ndim].
*/
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;
}
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b dofval [ndof].
* @return x dofval [ndof].
*/
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;
}
};
/**
* CRTP base class for a extra functions for a partitioned solver class.
*/
template <class D>
class MatrixSolverPartitionedBase {
public:
/**
* Underlying type.
*/
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:
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b nodevec [nelem, ndim].
* @param x nodevec [nelem, ndim].
* @return x nodevec [nelem, ndim].
*/
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;
}
/**
* Solve \f$ x = A^{-1} b \f$.
*
* @param A GooseFEM (sparse) matrix, see e.g. GooseFEM::Matrix().
* @param b dofval [ndof].
* @param x nodevec [nelem, ndim].
* @return x dofval [ndof].
*/
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;
}
};
/**
* CRTP base class for a matrix.
*/
template <class D>
class MatrixBase {
protected:
array_type::tensor<size_t, 2> m_conn; ///< Connectivity [#nelem, #nne].
array_type::tensor<size_t, 2> 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.
public:
/**
* Underlying type.
*/
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:
/**
* Number of elements.
* @return Unsigned integer.
*/
size_t nelem() const
{
return derived_cast().m_nelem;
}
/**
* Number of nodes per element.
* @return Unsigned integer.
*/
size_t nne() const
{
return derived_cast().m_nne;
}
/**
* Number of nodes.
* @return Unsigned integer.
*/
size_t nnode() const
{
return derived_cast().m_nnode;
}
/**
* Number of dimensions.
* @return Unsigned integer.
*/
size_t ndim() const
{
return derived_cast().m_ndim;
}
/**
* Number of DOFs.
* @return Unsigned integer.
*/
size_t ndof() const
{
return derived_cast().m_ndof;
}
/**
* DOFs per node.
* @return [#nnode, #ndim].
*/
const array_type::tensor<size_t, 2>& dofs() const
{
return derived_cast().m_dofs;
}
/**
* Connectivity.
* @return [#nelem, #nne].
*/
const array_type::tensor<size_t, 2>& conn() const
{
return derived_cast().m_conn;
}
/**
* Shape of "dofval".
* @return [#ndof].
*/
std::array<size_t, 1> shape_dofval() const
{
return std::array<size_t, 1>{derived_cast().m_ndof};
}
/**
* Shape of "nodevec".
* @return [#nnode, #ndim].
*/
std::array<size_t, 2> shape_nodevec() const
{
return std::array<size_t, 2>{derived_cast().m_nnode, derived_cast().m_ndim};
}
/**
* Shape of "elemmat".
* @return [#nelem, #nne * #ndim, #nne * #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
};
}
/**
* Assemble from "elemmat".
* @param elemmat [#nelem, #nne * #ndim, #nne * #ndim].
*/
template <class T>
void assemble(const T& elemmat)
{
GOOSEFEM_ASSERT(xt::has_shape(elemmat, this->shape_elemmat()));
derived_cast().assemble_impl(elemmat);
}
/**
* Copy as dense matrix.
* @return [#ndof, #ndof].
*/
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;
}
/**
* Copy to dense matrix.
* @param ret overwritten [#ndof, #ndof].
*/
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);
}
/**
* Dot-product \f$ b_i = A_{ij} x_j \f$.
*
* @param x nodevec [#nelem, #ndim].
* @return b nodevec [#nelem, #ndim].
*/
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;
}
/**
* Dot-product \f$ b_i = A_{ij} x_j \f$.
*
* @param x dofval [#ndof].
* @return b dofval [#ndof].
*/
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;
}
/**
* Dot-product \f$ b_i = A_{ij} x_j \f$.
*
* @param x nodevec [#nelem, #ndim].
* @param b (overwritten) nodevec [#nelem, #ndim].
*/
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);
}
/**
* Dot-product \f$ b_i = A_{ij} x_j \f$.
*
* @param x dofval [#ndof].
* @param b (overwritten) dofval [#ndof].
*/
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);
}
};
/**
* CRTP base class for a partitioned matrix.
*/
template <class D>
class MatrixPartitionedBase : public MatrixBase<D> {
protected:
array_type::tensor<size_t, 1> m_iiu; ///< See iiu()
array_type::tensor<size_t, 1> m_iip; ///< See iip()
size_t m_nnu; ///< See #nnu
size_t m_nnp; ///< See #nnp
public:
/**
* Underlying type.
*/
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:
/**
* Number of unknown DOFs.
* @return Unsigned integer.
*/
size_t nnu() const
{
return derived_cast().m_nnu;
}
/**
* Number of prescribed DOFs.
* @return Unsigned integer.
*/
size_t nnp() const
{
return derived_cast().m_nnp;
}
/**
* Unknown DOFs.
* @return [#nnu].
*/
const array_type::tensor<size_t, 1>& iiu() const
{
return derived_cast().m_iiu;
}
/**
* Prescribed DOFs.
* @return [#nnp].
*/
const array_type::tensor<size_t, 1>& iip() const
{
return derived_cast().m_iip;
}
/**
* Right-hand-size for corresponding to the prescribed DOFs:
*
* \f$ b_p = A_{pu} * x_u + 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.
*/
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;
}
/**
* Same as Reaction(const array_type::tensor<double, 2>&, const array_type::tensor<double, 2>&),
* but of "dofval" input and output.
*
* @param x "dofval" [#ndof].
* @param b "dofval" [#ndof].
* @return Copy of `b` with \f$ b_p \f$ overwritten.
*/
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;
}
/**
* Same as Reaction(const array_type::tensor<double, 2>&, const array_type::tensor<double, 2>&),
* but inserting in-place.
*
* @param x "nodevec" [#nnode, #ndim].
* @param b "nodevec" [#nnode, #ndim], \f$ b_p \f$ overwritten.
*/
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);
}
/**
* Same as Reaction(const array_type::tensor<double, 1>&, const array_type::tensor<double, 1>&),
* but inserting in-place.
*
* @param x "dofval" [#ndof].
* @param b "dofval" [#ndof], \f$ b_p \f$ overwritten.
*/
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);
}
/**
* Same as Reaction(const array_type::tensor<double, 1>&, const array_type::tensor<double, 1>&),
* but with partitioned input and output.
*
* @param x_u unknown "dofval" [#nnu].
* @param x_p prescribed "dofval" [#nnp].
* @return b_p prescribed "dofval" [#nnp].
*/
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;
}
/**
* Same as Reaction_p(const array_type::tensor<double, 1>&, const array_type::tensor<double,
* 1>&), 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 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);
}
};
/**
* CRTP base class for a partitioned matrix with tying.
*/
template <class D>
class MatrixPartitionedTyingsBase : public MatrixPartitionedBase<D> {
protected:
array_type::tensor<size_t, 1> m_iii; ///< See iii()
array_type::tensor<size_t, 1> m_iid; ///< See iid()
size_t m_nni; ///< See #nni
size_t m_nnd; ///< See #nnd
public:
/**
* Underlying type.
*/
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:
/**
* Number of independent DOFs.
* @return Unsigned integer.
*/
size_t nni() const
{
return derived_cast().m_nni;
}
/**
* Number of dependent DOFs.
* @return Unsigned integer.
*/
size_t nnd() const
{
return derived_cast().m_nnd;
}
/**
* Independent DOFs.
* @return [#nnu].
*/
const array_type::tensor<size_t, 1>& iii() const
{
return derived_cast().m_iii;
}
/**
* Dependent DOFs.
* @return [#nnp].
*/
const array_type::tensor<size_t, 1>& iid() const
{
return derived_cast().m_iid;
}
};
/**
* Sparse matrix.
*
* See GooseFEM::Vector() for bookkeeping definitions.
*/
class Matrix : public MatrixBase<Matrix> {
private:
friend MatrixBase<Matrix>;
private:
bool m_changed = true; ///< Signal changes to data.
Eigen::SparseMatrix<double> m_A; ///< The matrix.
std::vector<Eigen::Triplet<double>> m_T; ///< Matrix entries.
/**
* Class to solve the system (allowing single factorisation for multiple right-hand-sides).
*/
template <class>
friend class MatrixSolver;
public:
Matrix() = default;
/**
* Constructor.
*
* @param conn connectivity [#nelem, #nne].
* @param dofs DOFs per node [#nnode, #ndim].
*/
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);
}
/**
* Pointer to data.
*/
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:
/**
* Overwrite matrix.
*
* @param rows Row numbers [m].
* @param cols Column numbers [n].
* @param matrix Data entries `matrix(i, j)` for `rows(i), cols(j)` [m, n].
*/
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;
}
/**
* Add matrix.
*
* @param rows Row numbers [m].
* @param cols Column numbers [n].
* @param matrix Data entries `matrix(i, j)` for `rows(i), cols(j)` [m, n].
*/
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:
/**
* Convert "nodevec" to "dofval" (overwrite entries that occur more than once).
*
* @param nodevec input [#nnode, #ndim]
* @return dofval output [#ndof]
*/
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;
}
/**
* Convert "dofval" to "nodevec" (overwrite entries that occur more than once).
*
* @param dofval input [#ndof]
* @param nodevec output [#nnode, #ndim]
*/
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));
}
}
}
};
/**
* Solve \f$ x = A^{-1} b \f$, for `A` of the GooseFEM::Matrix() class.
* You can solve for multiple right-hand-sides using one factorisation.
*/
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; ///< Solver.
bool m_factor = true; ///< Signal to force factorization.
/**
* Compute inverse (evaluated by "solve").
*/
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
Event Timeline
Log In to Comment