12#ifndef GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
13#define GOOSEFEM_MATRIXPARTITIONEDTYINGS_H
19#include <Eigen/Sparse>
20#include <Eigen/SparseCholesky>
26class MatrixPartitionedTyingsSolver;
44 Eigen::SparseMatrix<double>
m_Auu;
45 Eigen::SparseMatrix<double>
m_Aup;
46 Eigen::SparseMatrix<double>
m_Apu;
47 Eigen::SparseMatrix<double>
m_App;
48 Eigen::SparseMatrix<double>
m_Aud;
49 Eigen::SparseMatrix<double>
m_Apd;
50 Eigen::SparseMatrix<double>
m_Adu;
51 Eigen::SparseMatrix<double>
m_Adp;
52 Eigen::SparseMatrix<double>
m_Add;
57 std::vector<Eigen::Triplet<double>>
m_Tuu;
58 std::vector<Eigen::Triplet<double>>
m_Tup;
59 std::vector<Eigen::Triplet<double>>
m_Tpu;
60 std::vector<Eigen::Triplet<double>>
m_Tpp;
61 std::vector<Eigen::Triplet<double>>
m_Tud;
62 std::vector<Eigen::Triplet<double>>
m_Tpd;
63 std::vector<Eigen::Triplet<double>>
m_Tdu;
64 std::vector<Eigen::Triplet<double>>
m_Tdp;
65 std::vector<Eigen::Triplet<double>>
m_Tdd;
66 Eigen::SparseMatrix<double>
m_Cdu;
67 Eigen::SparseMatrix<double>
m_Cdp;
68 Eigen::SparseMatrix<double>
m_Cud;
69 Eigen::SparseMatrix<double>
m_Cpd;
89 const Eigen::SparseMatrix<double>& Cdu,
90 const Eigen::SparseMatrix<double>& Cdp)
139 const Eigen::SparseMatrix<double>&
data_uu()
const
147 const Eigen::SparseMatrix<double>&
data_up()
const
155 const Eigen::SparseMatrix<double>&
data_pu()
const
163 const Eigen::SparseMatrix<double>&
data_pp()
const
171 const Eigen::SparseMatrix<double>&
data_ud()
const
179 const Eigen::SparseMatrix<double>&
data_pd()
const
187 const Eigen::SparseMatrix<double>&
data_du()
const
195 const Eigen::SparseMatrix<double>&
data_dp()
const
203 const Eigen::SparseMatrix<double>&
data_dd()
const
210 void assemble_impl(
const T& elemmat)
224 for (
size_t e = 0; e <
m_nelem; ++e) {
225 for (
size_t m = 0; m <
m_nne; ++m) {
226 for (
size_t i = 0; i <
m_ndim; ++i) {
230 for (
size_t n = 0; n <
m_nne; ++n) {
231 for (
size_t j = 0; j <
m_ndim; ++j) {
236 m_Tuu.push_back(Eigen::Triplet<double>(
240 m_Tup.push_back(Eigen::Triplet<double>(
243 else if (di <
m_nnu) {
244 m_Tud.push_back(Eigen::Triplet<double>(
248 m_Tpu.push_back(Eigen::Triplet<double>(
252 m_Tpp.push_back(Eigen::Triplet<double>(
257 else if (di <
m_nni) {
258 m_Tpd.push_back(Eigen::Triplet<double>(
263 else if (dj <
m_nnu) {
264 m_Tdu.push_back(Eigen::Triplet<double>(
267 else if (dj <
m_nni) {
268 m_Tdp.push_back(Eigen::Triplet<double>(
274 m_Tdd.push_back(Eigen::Triplet<double>(
299 void todense_impl(T& ret)
const
303 for (
int k = 0; k <
m_Auu.outerSize(); ++k) {
304 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Auu, k); it; ++it) {
305 ret(
m_iiu(it.row()),
m_iiu(it.col())) = it.value();
309 for (
int k = 0; k <
m_Aup.outerSize(); ++k) {
310 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Aup, k); it; ++it) {
311 ret(
m_iiu(it.row()),
m_iip(it.col())) = it.value();
315 for (
int k = 0; k <
m_Apu.outerSize(); ++k) {
316 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Apu, k); it; ++it) {
317 ret(
m_iip(it.row()),
m_iiu(it.col())) = it.value();
321 for (
int k = 0; k <
m_App.outerSize(); ++k) {
322 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_App, k); it; ++it) {
323 ret(
m_iip(it.row()),
m_iip(it.col())) = it.value();
327 for (
int k = 0; k <
m_Adu.outerSize(); ++k) {
328 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Adu, k); it; ++it) {
329 ret(
m_iid(it.row()),
m_iiu(it.col())) = it.value();
333 for (
int k = 0; k <
m_Adp.outerSize(); ++k) {
334 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Adp, k); it; ++it) {
335 ret(
m_iid(it.row()),
m_iip(it.col())) = it.value();
339 for (
int k = 0; k <
m_Aud.outerSize(); ++k) {
340 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Aud, k); it; ++it) {
341 ret(
m_iiu(it.row()),
m_iid(it.col())) = it.value();
345 for (
int k = 0; k <
m_Apd.outerSize(); ++k) {
346 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Apd, k); it; ++it) {
347 ret(
m_iip(it.row()),
m_iid(it.col())) = it.value();
351 for (
int k = 0; k <
m_Add.outerSize(); ++k) {
352 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Add, k); it; ++it) {
353 ret(
m_iid(it.row()),
m_iid(it.col())) = it.value();
359 void dot_nodevec_impl(
const T& x, T& b)
const
363 throw std::runtime_error(
"Not yet implemented");
367 void dot_dofval_impl(
const T& x, T& b)
const
371 throw std::runtime_error(
"Not yet implemented");
375 void reaction_nodevec_impl(
const T& x, T& b)
const
379 throw std::runtime_error(
"Not yet implemented");
383 void reaction_dofval_impl(
const T& x, T& b)
const
387 throw std::runtime_error(
"Not yet implemented");
390 void reaction_p_impl(
391 const array_type::tensor<double, 1>& x_u,
392 const array_type::tensor<double, 1>& x_p,
393 array_type::tensor<double, 1>& b_p)
const
398 throw std::runtime_error(
"Not yet implemented");
403 Eigen::VectorXd AsDofs_u(
const array_type::tensor<double, 1>& dofval)
const
407 Eigen::VectorXd dofval_u(
m_nnu, 1);
409#pragma omp parallel for
410 for (
size_t d = 0; d <
m_nnu; ++d) {
411 dofval_u(d) = dofval(
m_iiu(d));
417 Eigen::VectorXd AsDofs_u(
const array_type::tensor<double, 2>& nodevec)
const
421 Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(
m_nnu, 1);
423#pragma omp parallel for
424 for (
size_t m = 0; m <
m_nnode; ++m) {
425 for (
size_t i = 0; i <
m_ndim; ++i) {
427 dofval_u(
m_dofs(m, i)) = nodevec(m, i);
435 Eigen::VectorXd AsDofs_p(
const array_type::tensor<double, 1>& dofval)
const
439 Eigen::VectorXd dofval_p(
m_nnp, 1);
441#pragma omp parallel for
442 for (
size_t d = 0; d <
m_nnp; ++d) {
443 dofval_p(d) = dofval(
m_iip(d));
449 Eigen::VectorXd AsDofs_p(
const array_type::tensor<double, 2>& nodevec)
const
453 Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(
m_nnp, 1);
455#pragma omp parallel for
456 for (
size_t m = 0; m <
m_nnode; ++m) {
457 for (
size_t i = 0; i <
m_ndim; ++i) {
467 Eigen::VectorXd AsDofs_d(
const array_type::tensor<double, 1>& dofval)
const
471 Eigen::VectorXd dofval_d(
m_nnd, 1);
473#pragma omp parallel for
474 for (
size_t d = 0; d <
m_nnd; ++d) {
475 dofval_d(d) = dofval(
m_iip(d));
481 Eigen::VectorXd AsDofs_d(
const array_type::tensor<double, 2>& nodevec)
const
485 Eigen::VectorXd dofval_d = Eigen::VectorXd::Zero(
m_nnd, 1);
487#pragma omp parallel for
488 for (
size_t m = 0; m <
m_nnode; ++m) {
489 for (
size_t i = 0; i <
m_ndim; ++i) {
516template <
class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<
double>>>
533 Eigen::VectorXd B_u = A.AsDofs_u(b);
534 Eigen::VectorXd B_d = A.AsDofs_d(b);
535 Eigen::VectorXd X_p = A.AsDofs_p(x);
537 B_u += A.
m_Cud * B_d;
539 Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.
m_ACup * X_p));
540 Eigen::VectorXd X_d = A.
m_Cdu * X_u + A.
m_Cdp * X_p;
542#pragma omp parallel for
543 for (
size_t m = 0; m < A.
m_nnode; ++m) {
544 for (
size_t i = 0; i < A.
m_ndim; ++i) {
546 x(m, i) = X_u(A.
m_dofs(m, i));
560 Eigen::VectorXd B_u = A.AsDofs_u(b);
561 Eigen::VectorXd B_d = A.AsDofs_d(b);
562 Eigen::VectorXd X_p = A.AsDofs_p(x);
564 Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.
m_ACup * X_p));
565 Eigen::VectorXd X_d = A.
m_Cdu * X_u + A.
m_Cdp * X_p;
567#pragma omp parallel for
568 for (
size_t d = 0; d < A.
m_nnu; ++d) {
569 x(A.
m_iiu(d)) = X_u(d);
572#pragma omp parallel for
573 for (
size_t d = 0; d < A.
m_nnd; ++d) {
574 x(A.
m_iid(d)) = X_d(d);
597 this->
solve_u(A, b_u, b_d, x_p, x_u);
628 Eigen::Map<Eigen::VectorXd>(x_u.data(), x_u.size()).noalias() =
629 m_solver.solve(Eigen::VectorXd(
630 Eigen::Map<const Eigen::VectorXd>(b_u.data(), b_u.size()) -
631 A.
m_ACup * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size())));
636 bool m_factor =
true;
657 m_solver.compute(A.
m_ACuu);
CRTP base class for a matrix.
const array_type::tensor< size_t, 2 > & dofs() const
DOFs per node.
const array_type::tensor< size_t, 2 > & conn() const
Connectivity.
array_type::tensor< size_t, 2 > m_dofs
DOF-numbers per node [nnode, ndim].
array_type::tensor< size_t, 2 > m_conn
Connectivity [nelem, nne].
size_t m_nelem
See nelem().
bool m_changed
Signal changes to data.
size_t m_nnode
See nnode().
CRTP base class for a partitioned matrix.
array_type::tensor< size_t, 1 > m_iiu
See iiu()
array_type::tensor< size_t, 1 > m_iip
See iip()
CRTP base class for a partitioned matrix with tying.
array_type::tensor< size_t, 1 > m_iii
See iii()
array_type::tensor< size_t, 1 > m_iid
See iid()
Solver for MatrixPartitionedTyings().
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)
Same as Solve(MatrixPartitionedTyings&, const array_type::tensor<double, 2>&, const array_type::tenso...
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)
Same as Solve_u(MatrixPartitionedTyings&, const array_type::tensor<double, 1>&, const array_type::ten...
Sparse matrix from with dependent DOFs are eliminated, and the remaining (small) independent system i...
std::vector< Eigen::Triplet< double > > m_Tpp
Matrix entries.
Eigen::SparseMatrix< double > m_Auu
The matrix.
const Eigen::SparseMatrix< double > & data_pd() const
Pointer to data.
std::vector< Eigen::Triplet< double > > m_Tup
Matrix entries.
Eigen::SparseMatrix< double > m_Cdu
Tying matrix, see Tyings::Periodic::Cdu().
Eigen::SparseMatrix< double > m_App
The matrix.
Eigen::SparseMatrix< double > m_ACup
// The matrix for which the tyings have been applied.
Eigen::SparseMatrix< double > m_Aup
The matrix.
Eigen::SparseMatrix< double > m_Adp
The matrix.
Eigen::SparseMatrix< double > m_Apd
The matrix.
std::vector< Eigen::Triplet< double > > m_Tdu
Matrix entries.
std::vector< Eigen::Triplet< double > > m_Tdd
Matrix entries.
const Eigen::SparseMatrix< double > & data_uu() const
Pointer to data.
Eigen::SparseMatrix< double > m_Apu
The matrix.
const Eigen::SparseMatrix< double > & data_dd() const
Pointer to data.
const Eigen::SparseMatrix< double > & data_pu() const
Pointer to data.
const Eigen::SparseMatrix< double > & data_ud() const
Pointer to data.
Eigen::SparseMatrix< double > m_ACpp
// The matrix for which the tyings have been applied.
std::vector< Eigen::Triplet< double > > m_Tpd
Matrix entries.
Eigen::SparseMatrix< double > m_Add
The matrix.
Eigen::SparseMatrix< double > m_Aud
The matrix.
Eigen::SparseMatrix< double > m_ACpu
// The matrix for which the tyings have been applied.
const Eigen::SparseMatrix< double > & data_pp() const
Pointer to data.
const Eigen::SparseMatrix< double > & data_up() const
Pointer to data.
Eigen::SparseMatrix< double > m_Adu
The matrix.
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)
Constructor.
const Eigen::SparseMatrix< double > & data_du() const
Pointer to data.
std::vector< Eigen::Triplet< double > > m_Tud
Matrix entries.
Eigen::SparseMatrix< double > m_Cpd
Transpose of "m_Cdp".
Eigen::SparseMatrix< double > m_ACuu
// The matrix for which the tyings have been applied.
std::vector< Eigen::Triplet< double > > m_Tuu
Matrix entries.
const Eigen::SparseMatrix< double > & data_dp() const
Pointer to data.
Eigen::SparseMatrix< double > m_Cdp
Tying matrix, see Tyings::Periodic::Cdp().
std::vector< Eigen::Triplet< double > > m_Tpu
Matrix entries.
std::vector< Eigen::Triplet< double > > m_Tdp
Matrix entries.
Eigen::SparseMatrix< double > m_Cud
Transpose of "m_Cdu".
CRTP base class for a solver class.
CRTP base class for a extra functions for a partitioned solver class.
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.