11#ifndef GOOSEFEM_MATRIXPARTITIONED_H
12#define GOOSEFEM_MATRIXPARTITIONED_H
18#include <Eigen/Sparse>
19#include <Eigen/SparseCholesky>
25class MatrixPartitionedSolver;
40 Eigen::SparseMatrix<double>
m_Auu;
41 Eigen::SparseMatrix<double>
m_Aup;
42 Eigen::SparseMatrix<double>
m_Apu;
43 Eigen::SparseMatrix<double>
m_App;
45 std::vector<Eigen::Triplet<double>>
m_Tuu;
46 std::vector<Eigen::Triplet<double>>
m_Tup;
47 std::vector<Eigen::Triplet<double>>
m_Tpu;
48 std::vector<Eigen::Triplet<double>>
m_Tpp;
125 const Eigen::SparseMatrix<double>&
data_uu()
const
133 const Eigen::SparseMatrix<double>&
data_up()
const
141 const Eigen::SparseMatrix<double>&
data_pu()
const
149 const Eigen::SparseMatrix<double>&
data_pp()
const
156 void assemble_impl(
const T& elemmat)
165 for (
size_t e = 0; e <
m_nelem; ++e) {
166 for (
size_t m = 0; m <
m_nne; ++m) {
167 for (
size_t i = 0; i <
m_ndim; ++i) {
171 for (
size_t n = 0; n <
m_nne; ++n) {
172 for (
size_t j = 0; j <
m_ndim; ++j) {
177 m_Tuu.push_back(Eigen::Triplet<double>(
180 else if (di <
m_nnu) {
181 m_Tup.push_back(Eigen::Triplet<double>(
184 else if (dj <
m_nnu) {
185 m_Tpu.push_back(Eigen::Triplet<double>(
189 m_Tpp.push_back(Eigen::Triplet<double>(
225 std::vector<Eigen::Triplet<double>> Tuu;
226 std::vector<Eigen::Triplet<double>> Tup;
227 std::vector<Eigen::Triplet<double>> Tpu;
228 std::vector<Eigen::Triplet<double>> Tpp;
230 for (
size_t i = 0; i < rows.size(); ++i) {
231 for (
size_t j = 0; j < cols.size(); ++j) {
234 double v = matrix(i, j);
237 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
239 else if (di <
m_nnu) {
240 Tup.push_back(Eigen::Triplet<double>(di, dj -
m_nnu, v));
242 else if (dj <
m_nnu) {
243 Tpu.push_back(Eigen::Triplet<double>(di -
m_nnu, dj, v));
246 Tpp.push_back(Eigen::Triplet<double>(di -
m_nnu, dj -
m_nnu, v));
251 m_Auu.setFromTriplets(Tuu.begin(), Tuu.end());
252 m_Aup.setFromTriplets(Tup.begin(), Tup.end());
253 m_Apu.setFromTriplets(Tpu.begin(), Tpu.end());
254 m_App.setFromTriplets(Tpp.begin(), Tpp.end());
275 std::vector<Eigen::Triplet<double>> Tuu;
276 std::vector<Eigen::Triplet<double>> Tup;
277 std::vector<Eigen::Triplet<double>> Tpu;
278 std::vector<Eigen::Triplet<double>> Tpp;
280 Eigen::SparseMatrix<double> Auu(
m_nnu,
m_nnu);
281 Eigen::SparseMatrix<double> Aup(
m_nnu,
m_nnp);
282 Eigen::SparseMatrix<double> Apu(
m_nnp,
m_nnu);
283 Eigen::SparseMatrix<double> App(
m_nnp,
m_nnp);
285 for (
size_t i = 0; i < rows.size(); ++i) {
286 for (
size_t j = 0; j < cols.size(); ++j) {
289 double v = matrix(i, j);
292 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
294 else if (di <
m_nnu) {
295 Tup.push_back(Eigen::Triplet<double>(di, dj -
m_nnu, v));
297 else if (dj <
m_nnu) {
298 Tpu.push_back(Eigen::Triplet<double>(di -
m_nnu, dj, v));
301 Tpp.push_back(Eigen::Triplet<double>(di -
m_nnu, dj -
m_nnu, v));
306 Auu.setFromTriplets(Tuu.begin(), Tuu.end());
307 Aup.setFromTriplets(Tup.begin(), Tup.end());
308 Apu.setFromTriplets(Tpu.begin(), Tpu.end());
309 App.setFromTriplets(Tpp.begin(), Tpp.end());
319 void todense_impl(T& ret)
const
323 for (
int k = 0; k <
m_Auu.outerSize(); ++k) {
324 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Auu, k); it; ++it) {
325 ret(
m_iiu(it.row()),
m_iiu(it.col())) = it.value();
329 for (
int k = 0; k <
m_Aup.outerSize(); ++k) {
330 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Aup, k); it; ++it) {
331 ret(
m_iiu(it.row()),
m_iip(it.col())) = it.value();
335 for (
int k = 0; k <
m_Apu.outerSize(); ++k) {
336 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_Apu, k); it; ++it) {
337 ret(
m_iip(it.row()),
m_iiu(it.col())) = it.value();
341 for (
int k = 0; k <
m_App.outerSize(); ++k) {
342 for (Eigen::SparseMatrix<double>::InnerIterator it(
m_App, k); it; ++it) {
343 ret(
m_iip(it.row()),
m_iip(it.col())) = it.value();
349 void dot_nodevec_impl(
const T& x, T& b)
const
354 Eigen::VectorXd X_u = this->AsDofs_u(x);
355 Eigen::VectorXd X_p = this->AsDofs_p(x);
356 Eigen::VectorXd B_u =
m_Auu * X_u +
m_Aup * X_p;
357 Eigen::VectorXd B_p =
m_Apu * X_u +
m_App * X_p;
359#pragma omp parallel for
360 for (
size_t m = 0; m <
m_nnode; ++m) {
361 for (
size_t i = 0; i <
m_ndim; ++i) {
363 b(m, i) = B_u(
m_part(m, i));
373 void dot_dofval_impl(
const T& x, T& b)
const
378 Eigen::VectorXd X_u = this->AsDofs_u(x);
379 Eigen::VectorXd X_p = this->AsDofs_p(x);
381 Eigen::VectorXd B_u =
m_Auu * X_u +
m_Aup * X_p;
382 Eigen::VectorXd B_p =
m_Apu * X_u +
m_App * X_p;
384#pragma omp parallel for
385 for (
size_t d = 0; d <
m_nnu; ++d) {
386 b(
m_iiu(d)) = B_u(d);
389#pragma omp parallel for
390 for (
size_t d = 0; d <
m_nnp; ++d) {
391 b(
m_iip(d)) = B_p(d);
396 void reaction_nodevec_impl(
const T& x, T& b)
const
401 Eigen::VectorXd X_u = this->AsDofs_u(x);
402 Eigen::VectorXd X_p = this->AsDofs_p(x);
403 Eigen::VectorXd B_p =
m_Apu * X_u +
m_App * X_p;
405#pragma omp parallel for
406 for (
size_t m = 0; m <
m_nnode; ++m) {
407 for (
size_t i = 0; i <
m_ndim; ++i) {
416 void reaction_dofval_impl(
const T& x, T& b)
const
421 Eigen::VectorXd X_u = this->AsDofs_u(x);
422 Eigen::VectorXd X_p = this->AsDofs_p(x);
423 Eigen::VectorXd B_p =
m_Apu * X_u +
m_App * X_p;
425#pragma omp parallel for
426 for (
size_t d = 0; d <
m_nnp; ++d) {
427 b(
m_iip(d)) = B_p(d);
431 void reaction_p_impl(
432 const array_type::tensor<double, 1>& x_u,
433 const array_type::tensor<double, 1>& x_p,
434 array_type::tensor<double, 1>& b_p)
const
440 Eigen::Map<Eigen::VectorXd>(b_p.data(), b_p.size()).noalias() =
441 m_Apu * Eigen::Map<const Eigen::VectorXd>(x_u.data(), x_u.size()) +
442 m_App * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size());
447 Eigen::VectorXd AsDofs_u(
const array_type::tensor<double, 1>& dofval)
const
451 Eigen::VectorXd dofval_u(
m_nnu, 1);
453#pragma omp parallel for
454 for (
size_t d = 0; d <
m_nnu; ++d) {
455 dofval_u(d) = dofval(
m_iiu(d));
461 Eigen::VectorXd AsDofs_u(
const array_type::tensor<double, 2>& nodevec)
const
465 Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(
m_nnu, 1);
467#pragma omp parallel for
468 for (
size_t m = 0; m <
m_nnode; ++m) {
469 for (
size_t i = 0; i <
m_ndim; ++i) {
471 dofval_u(
m_part(m, i)) = nodevec(m, i);
479 Eigen::VectorXd AsDofs_p(
const array_type::tensor<double, 1>& dofval)
const
483 Eigen::VectorXd dofval_p(
m_nnp, 1);
485#pragma omp parallel for
486 for (
size_t d = 0; d <
m_nnp; ++d) {
487 dofval_p(d) = dofval(
m_iip(d));
493 Eigen::VectorXd AsDofs_p(
const array_type::tensor<double, 2>& nodevec)
const
497 Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(
m_nnp, 1);
499#pragma omp parallel for
500 for (
size_t m = 0; m <
m_nnode; ++m) {
501 for (
size_t i = 0; i <
m_ndim; ++i) {
519template <
class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<
double>>>
547 this->solve_u_impl(A, b_u, x_p, x_u);
570 this->solve_u_impl(A, b_u, x_p, x_u);
578 Eigen::VectorXd B_u = A.AsDofs_u(b);
579 Eigen::VectorXd X_p = A.AsDofs_p(x);
580 Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.
m_Aup * X_p));
582#pragma omp parallel for
583 for (
size_t m = 0; m < A.
m_nnode; ++m) {
584 for (
size_t i = 0; i < A.
m_ndim; ++i) {
586 x(m, i) = X_u(A.
m_part(m, i));
593 void solve_dofval_impl(MatrixPartitioned& A,
const T& b, T& x)
596 Eigen::VectorXd B_u = A.AsDofs_u(b);
597 Eigen::VectorXd X_p = A.AsDofs_p(x);
598 Eigen::VectorXd X_u = m_solver.solve(Eigen::VectorXd(B_u - A.m_Aup * X_p));
600#pragma omp parallel for
601 for (
size_t d = 0; d < A.m_nnu; ++d) {
602 x(A.m_iiu(d)) = X_u(d);
607 void solve_u_impl(MatrixPartitioned& A,
const T& b_u,
const T& x_p, T& x_u)
611 Eigen::Map<Eigen::VectorXd>(x_u.data(), x_u.size()).noalias() =
612 m_solver.solve(Eigen::VectorXd(
613 Eigen::Map<const Eigen::VectorXd>(b_u.data(), b_u.size()) -
614 A.m_Aup * Eigen::Map<const Eigen::VectorXd>(x_p.data(), x_p.size())));
619 bool m_factor =
true;
624 void factorize(MatrixPartitioned& A)
626 if (!A.m_changed && !m_factor) {
629 m_solver.compute(A.m_Auu);
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()
const array_type::tensor< size_t, 1 > & iip() const
Prescribed DOFs.
Solve for A of the MatrixPartitioned() class.
array_type::tensor< double, 1 > Solve_u(M &A, const array_type::tensor< double, 1 > &b_u, const array_type::tensor< double, 1 > &x_p)
Solve .
void solve_u(M &A, const array_type::tensor< double, 1 > &b_u, const array_type::tensor< double, 1 > &x_p, array_type::tensor< double, 1 > &x_u)
Same as Solve .
Sparse matrix partitioned in an unknown and a prescribed part.
std::vector< Eigen::Triplet< double > > m_Tuu
Matrix entries.
const Eigen::SparseMatrix< double > & data_uu() const
Pointer to data.
const Eigen::SparseMatrix< double > & data_pu() const
Pointer to data.
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)
Add matrix.
const Eigen::SparseMatrix< double > & data_up() const
Pointer to data.
array_type::tensor< size_t, 1 > m_part1d
Map real DOF to DOF in partitioned system.
Eigen::SparseMatrix< double > m_App
The matrix.
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)
Overwrite matrix.
Eigen::SparseMatrix< double > m_Aup
The matrix.
std::vector< Eigen::Triplet< double > > m_Tup
Matrix entries.
array_type::tensor< size_t, 2 > m_part
Renumbered DOFs per node, such that:
Eigen::SparseMatrix< double > m_Auu
The matrix.
std::vector< Eigen::Triplet< double > > m_Tpu
Matrix entries.
const Eigen::SparseMatrix< double > & data_pp() const
Pointer to data.
Eigen::SparseMatrix< double > m_Apu
The matrix.
MatrixPartitioned(const array_type::tensor< size_t, 2 > &conn, const array_type::tensor< size_t, 2 > &dofs, const array_type::tensor< size_t, 1 > &iip)
Constructor.
std::vector< Eigen::Triplet< double > > m_Tpp
Matrix entries.
CRTP base class for a solver class.
CRTP base class for a extra functions for a partitioned solver class.
Reorder to lowest possible index, in specific order.
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Toolbox to perform finite element computations.
bool is_unique(const T &arg)
Returns true is a list is unique (has not duplicate items).