11#ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
12#define GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
61 m_iiu = xt::setdiff1d(
dofs,
iip);
65 m_Auu = xt::empty<double>({m_nnu});
66 m_App = xt::empty<double>({m_nnp});
67 m_inv_uu = xt::empty<double>({m_nnu});
72 void assemble_impl(
const T& elemmat)
80 for (
size_t e = 0; e <
m_nelem; ++e) {
81 for (
size_t m = 0; m <
m_nne; ++m) {
82 for (
size_t i = 0; i <
m_ndim; ++i) {
84 size_t d = m_part(
m_conn(e, m), i);
90 m_App(d - m_nnu) += elemmat(e, m *
m_ndim + i, m *
m_ndim + i);
100 void todense_impl(T& ret)
const
104#pragma omp parallel for
105 for (
size_t d = 0; d < m_nnu; ++d) {
106 ret(m_iiu(d), m_iiu(d)) = m_Auu(d);
109#pragma omp parallel for
110 for (
size_t d = 0; d < m_nnp; ++d) {
111 ret(m_iip(d), m_iip(d)) = m_App(d);
124#pragma omp parallel for
125 for (
size_t d = 0; d < m_nnu; ++d) {
126 m_Auu(d) = A(m_iiu(d));
129#pragma omp parallel for
130 for (
size_t d = 0; d < m_nnp; ++d) {
131 m_App(d) = A(m_iip(d));
145#pragma omp parallel for
146 for (
size_t d = 0; d < m_nnu; ++d) {
147 ret(m_iiu(d)) = m_Auu(d);
150#pragma omp parallel for
151 for (
size_t d = 0; d < m_nnp; ++d) {
152 ret(m_iip(d)) = m_App(d);
187 void dot_nodevec_impl(
const T& x, T& b)
const
192#pragma omp parallel for
193 for (
size_t m = 0; m <
m_nnode; ++m) {
194 for (
size_t i = 0; i <
m_ndim; ++i) {
196 size_t d = m_part(m, i);
199 b(m, i) = m_Auu(d) * x(m, i);
202 b(m, i) = m_App(d - m_nnu) * x(m, i);
209 void dot_dofval_impl(
const T& x, T& b)
const
214#pragma omp parallel for
215 for (
size_t d = 0; d < m_nnu; ++d) {
216 b(m_iiu(d)) = m_Auu(d) * x(m_iiu(d));
219#pragma omp parallel for
220 for (
size_t d = 0; d < m_nnp; ++d) {
221 b(m_iip(d)) = m_App(d) * x(m_iip(d));
232 array_type::tensor<double, 1>
236 this->
dot_u(x_u, x_p, b_u);
257#pragma omp parallel for
258 for (
size_t d = 0; d < m_nnu; ++d) {
259 b_u(d) = m_Auu(d) * x_u(d);
273 this->
dot_p(x_u, x_p, b_p);
294#pragma omp parallel for
295 for (
size_t d = 0; d < m_nnp; ++d) {
296 b_p(d) = m_App(d) * x_p(d);
302 void solve_nodevec_impl(
const T& b, T& x)
309#pragma omp parallel for
310 for (
size_t m = 0; m <
m_nnode; ++m) {
311 for (
size_t i = 0; i <
m_ndim; ++i) {
312 if (m_part(m, i) < m_nnu) {
313 x(m, i) = m_inv_uu(m_part(m, i)) * b(m, i);
320 void solve_dofval_impl(
const T& b, T& x)
327#pragma omp parallel for
328 for (
size_t d = 0; d < m_nnu; ++d) {
329 x(m_iiu(d)) = m_inv_uu(d) * b(m_iiu(d));
339 array_type::tensor<double, 1>
365#pragma omp parallel for
366 for (
size_t d = 0; d < m_nnu; ++d) {
367 x_u(d) = m_inv_uu(d) * b_u(d);
373 void reaction_nodevec_impl(
const T& x, T& b)
const
378#pragma omp parallel for
379 for (
size_t m = 0; m <
m_nnode; ++m) {
380 for (
size_t i = 0; i <
m_ndim; ++i) {
381 if (m_part(m, i) >= m_nnu) {
382 b(m, i) = m_App(m_part(m, i) - m_nnu) * x(m, i);
389 void reaction_dofval_impl(
const T& x, T& b)
const
394#pragma omp parallel for
395 for (
size_t d = 0; d < m_nnp; ++d) {
396 b(m_iip(d)) = m_App(d) * x(m_iip(d));
400 void reaction_p_impl(
401 const array_type::tensor<double, 1>& x_u,
402 const array_type::tensor<double, 1>& x_p,
403 array_type::tensor<double, 1>& b_p)
const
411#pragma omp parallel for
412 for (
size_t d = 0; d < m_nnp; ++d) {
413 b_p(d) = m_App(d) * x_p(d);
419 array_type::tensor<double, 1> m_Auu;
420 array_type::tensor<double, 1> m_App;
421 array_type::tensor<double, 1> m_inv_uu;
424 array_type::tensor<size_t, 2> m_part;
425 array_type::tensor<size_t, 1> m_iiu;
426 array_type::tensor<size_t, 1> m_iip;
439#pragma omp parallel for
440 for (
size_t d = 0; d < m_nnu; ++d) {
441 m_inv_uu(d) = 1.0 / m_Auu(d);
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 with tying.
Diagonal and partitioned matrix.
array_type::tensor< double, 1 > Solve_u(const array_type::tensor< double, 1 > &b_u, const array_type::tensor< double, 1 > &x_p)
MatrixDiagonalPartitioned(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.
void dot_u(const array_type::tensor< double, 1 > &x_u, const array_type::tensor< double, 1 > &x_p, array_type::tensor< double, 1 > &b_u) const
const array_type::tensor< double, 1 > & data_uu() const
Pointer to data.
array_type::tensor< double, 1 > Dot_p(const array_type::tensor< double, 1 > &x_u, const array_type::tensor< double, 1 > &x_p) const
array_type::tensor< double, 1 > data() const
Assemble to diagonal matrix (involves copies).
array_type::tensor< double, 1 > Todiagonal() const
Pointer to data.
array_type::tensor< double, 1 > Dot_u(const array_type::tensor< double, 1 > &x_u, const array_type::tensor< double, 1 > &x_p) const
void solve_u(const array_type::tensor< double, 1 > &b_u, const array_type::tensor< double, 1 > &x_p, array_type::tensor< double, 1 > &x_u)
void set(const array_type::tensor< double, 1 > &A)
Set all (diagonal) matrix components.
const array_type::tensor< double, 1 > & data_pp() const
Pointer to data.
void dot_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
CRTP base class for a partitioned matrix.
const array_type::tensor< size_t, 1 > & iip() const
Prescribed DOFs.
Reorder to lowest possible index, in specific order.
#define GOOSEFEM_ASSERT(expr)
All assertions are implementation as::
bool isDiagonal(const array_type::tensor< double, 3 > &elemmat)
Check that all of the matrices stored per elemmat (shape: [nelem, nne * ndim, nne * ndim]) are diagon...
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).