9#ifndef GOOSEFEM_MATRIXDIAGONAL_H
10#define GOOSEFEM_MATRIXDIAGONAL_H
56 derived_cast().solve_nodevec_impl(b, x);
73 derived_cast().solve_dofval_impl(b, x);
89 derived_cast().solve_nodevec_impl(b, x);
104 derived_cast().solve_dofval_impl(b, x);
132 template <
class C,
class D>
142 m_A = xt::empty<double>({
m_ndof});
143 m_inv = xt::empty<double>({
m_ndof});
151 void assemble_impl(
const T& elemmat)
158 for (
size_t e = 0; e <
m_nelem; ++e) {
159 for (
size_t m = 0; m <
m_nne; ++m) {
160 for (
size_t i = 0; i <
m_ndim; ++i) {
170 void todense_impl(T& ret)
const
174#pragma omp parallel for
175 for (
size_t d = 0; d <
m_ndof; ++d) {
188 std::copy(A.begin(), A.end(), m_A.begin());
212 void dot_nodevec_impl(
const T& x, T& b)
const
217#pragma omp parallel for
218 for (
size_t m = 0; m <
m_nnode; ++m) {
219 for (
size_t i = 0; i <
m_ndim; ++i) {
220 b(m, i) = m_A(
m_dofs(m, i)) * x(m, i);
226 void dot_dofval_impl(
const T& x, T& b)
const
231 xt::noalias(b) = m_A * x;
236 void solve_nodevec_impl(
const T& b, T& x)
243#pragma omp parallel for
244 for (
size_t m = 0; m <
m_nnode; ++m) {
245 for (
size_t i = 0; i <
m_ndim; ++i) {
246 x(m, i) = m_inv(
m_dofs(m, i)) * b(m, i);
252 void solve_dofval_impl(
const T& b, T& x)
257 xt::noalias(x) = m_inv * b;
261 array_type::tensor<double, 1> m_A;
262 array_type::tensor<double, 1> m_inv;
273#pragma omp parallel for
274 for (
size_t d = 0; d <
m_ndof; ++d) {
275 m_inv(d) = 1.0 / m_A(d);
Convenience methods for integration point data.
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.
void solve(const array_type::tensor< double, 2 > &b, array_type::tensor< double, 2 > &x)
Solve .
void solve(const array_type::tensor< double, 1 > &b, array_type::tensor< double, 1 > &x)
Solve .
array_type::tensor< double, 2 > Solve(const array_type::tensor< double, 2 > &b)
Solve .
array_type::tensor< double, 1 > Solve(const array_type::tensor< double, 1 > &b)
Solve .
D derived_type
Underlying type.
const array_type::tensor< double, 1 > & Todiagonal() const
Copy as diagonal matrix.
void set(const array_type::tensor< double, 1 > &A)
Set all (diagonal) matrix components.
const array_type::tensor< double, 1 > & data() const
Underlying matrix.
MatrixDiagonal(const C &conn, const D &dofs)
Constructor.
#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.