FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
MatrixDiagonal.h
Go to the documentation of this file.
1
9#ifndef GOOSEFEM_MATRIXDIAGONAL_H
10#define GOOSEFEM_MATRIXDIAGONAL_H
11
12#include "Element.h"
13#include "Matrix.h"
14#include "config.h"
15
16namespace GooseFEM {
17
21template <class D>
23public:
27 using derived_type = D;
28
29private:
30 auto derived_cast() -> derived_type&
31 {
32 return *static_cast<derived_type*>(this);
33 }
34
35 auto derived_cast() const -> const derived_type&
36 {
37 return *static_cast<const derived_type*>(this);
38 }
39
40public:
54 {
55 array_type::tensor<double, 2> x = xt::empty_like(b);
56 derived_cast().solve_nodevec_impl(b, x);
57 return x;
58 }
59
71 {
72 array_type::tensor<double, 1> x = xt::empty_like(b);
73 derived_cast().solve_dofval_impl(b, x);
74 return x;
75 }
76
88 {
89 derived_cast().solve_nodevec_impl(b, x);
90 }
91
103 {
104 derived_cast().solve_dofval_impl(b, x);
105 }
106};
107
115class MatrixDiagonal : public MatrixBase<MatrixDiagonal>,
116 public MatrixDiagonalBase<MatrixDiagonal> {
117private:
120
121public:
122 MatrixDiagonal() = default;
123
132 template <class C, class D>
133 MatrixDiagonal(const C& conn, const D& dofs)
134 {
135 m_conn = conn;
136 m_dofs = dofs;
137 m_nelem = m_conn.shape(0);
138 m_nne = m_conn.shape(1);
139 m_nnode = m_dofs.shape(0);
140 m_ndim = m_dofs.shape(1);
141 m_ndof = xt::amax(m_dofs)() + 1;
142 m_A = xt::empty<double>({m_ndof});
143 m_inv = xt::empty<double>({m_ndof});
144
145 GOOSEFEM_ASSERT(xt::amax(m_conn)() + 1 <= m_nnode);
147 }
148
149private:
150 template <class T>
151 void assemble_impl(const T& elemmat)
152 {
153 GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
155
156 m_A.fill(0.0);
157
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) {
161 m_A(m_dofs(m_conn(e, m), i)) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
162 }
163 }
164 }
165
166 m_changed = true;
167 }
168
169 template <class T>
170 void todense_impl(T& ret) const
171 {
172 ret.fill(0.0);
173
174#pragma omp parallel for
175 for (size_t d = 0; d < m_ndof; ++d) {
176 ret(d, d) = m_A(d);
177 }
178 }
179
180public:
186 {
187 GOOSEFEM_ASSERT(A.size() == m_ndof);
188 std::copy(A.begin(), A.end(), m_A.begin());
189 m_changed = true;
190 }
191
196 [[deprecated]] const array_type::tensor<double, 1>& Todiagonal() const
197 {
198 return m_A;
199 }
200
206 {
207 return m_A;
208 }
209
210private:
211 template <class T>
212 void dot_nodevec_impl(const T& x, T& b) const
213 {
214 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
215 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
216
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);
221 }
222 }
223 }
224
225 template <class T>
226 void dot_dofval_impl(const T& x, T& b) const
227 {
228 GOOSEFEM_ASSERT(x.size() == m_ndof);
229 GOOSEFEM_ASSERT(b.size() == m_ndof);
230
231 xt::noalias(b) = m_A * x;
232 }
233
234private:
235 template <class T>
236 void solve_nodevec_impl(const T& b, T& x)
237 {
238 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
239 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
240
241 this->factorize();
242
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);
247 }
248 }
249 }
250
251 template <class T>
252 void solve_dofval_impl(const T& b, T& x)
253 {
254 GOOSEFEM_ASSERT(b.size() == m_ndof);
255 GOOSEFEM_ASSERT(x.size() == m_ndof);
256 this->factorize();
257 xt::noalias(x) = m_inv * b;
258 }
259
260private:
261 array_type::tensor<double, 1> m_A;
262 array_type::tensor<double, 1> m_inv;
263
267 void factorize()
268 {
269 if (!m_changed) {
270 return;
271 }
272
273#pragma omp parallel for
274 for (size_t d = 0; d < m_ndof; ++d) {
275 m_inv(d) = 1.0 / m_A(d);
276 }
277
278 m_changed = false;
279 }
280};
281
282} // namespace GooseFEM
283
284#endif
Convenience methods for integration point data.
CRTP base class for a matrix.
Definition: Matrix.h:198
const array_type::tensor< size_t, 2 > & dofs() const
DOFs per node.
Definition: Matrix.h:278
const array_type::tensor< size_t, 2 > & conn() const
Connectivity.
Definition: Matrix.h:287
array_type::tensor< size_t, 2 > m_dofs
DOF-numbers per node [#nnode, #ndim].
Definition: Matrix.h:201
array_type::tensor< size_t, 2 > m_conn
Connectivity [#nelem, #nne].
Definition: Matrix.h:200
bool m_changed
Signal changes to data.
Definition: Matrix.h:209
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::
Definition: config.h:96
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...
Definition: Element.h:120
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:176
Toolbox to perform finite element computations.
Definition: Allocate.h:14