FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
MatrixDiagonalPartitioned.h
Go to the documentation of this file.
1
11#ifndef GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
12#define GOOSEFEM_MATRIXDIAGONALPARTITIONED_H
13
14#include "MatrixDiagonal.h"
15#include "Mesh.h"
16#include "config.h"
17
18namespace GooseFEM {
19
25class MatrixDiagonalPartitioned : public MatrixPartitionedBase<MatrixDiagonalPartitioned>,
26 public MatrixDiagonalBase<MatrixDiagonalPartitioned> {
27private:
31
32public:
33 MatrixDiagonalPartitioned() = default;
34
46 {
47 m_conn = conn;
48 m_dofs = dofs;
49 m_nelem = m_conn.shape(0);
50 m_nne = m_conn.shape(1);
51 m_nnode = m_dofs.shape(0);
52 m_ndim = m_dofs.shape(1);
53 m_ndof = xt::amax(m_dofs)() + 1;
54
56 GOOSEFEM_ASSERT(xt::amax(conn)() + 1 <= m_nnode);
58 GOOSEFEM_ASSERT(xt::amax(iip)() <= xt::amax(dofs)());
59
60 m_iip = iip;
61 m_iiu = xt::setdiff1d(dofs, iip);
62 m_nnp = m_iip.size();
63 m_nnu = m_iiu.size();
64 m_part = Mesh::Reorder({m_iiu, m_iip}).apply(m_dofs);
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});
68 }
69
70private:
71 template <class T>
72 void assemble_impl(const T& elemmat)
73 {
74 GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
76
77 m_Auu.fill(0.0);
78 m_App.fill(0.0);
79
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) {
83
84 size_t d = m_part(m_conn(e, m), i);
85
86 if (d < m_nnu) {
87 m_Auu(d) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
88 }
89 else {
90 m_App(d - m_nnu) += elemmat(e, m * m_ndim + i, m * m_ndim + i);
91 }
92 }
93 }
94 }
95
96 m_changed = true;
97 }
98
99 template <class T>
100 void todense_impl(T& ret) const
101 {
102 ret.fill(0.0);
103
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);
107 }
108
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);
112 }
113 }
114
115public:
121 {
122 GOOSEFEM_ASSERT(A.size() == m_ndof);
123
124#pragma omp parallel for
125 for (size_t d = 0; d < m_nnu; ++d) {
126 m_Auu(d) = A(m_iiu(d));
127 }
128
129#pragma omp parallel for
130 for (size_t d = 0; d < m_nnp; ++d) {
131 m_App(d) = A(m_iip(d));
132 }
133
134 m_changed = true;
135 }
136
142 {
143 array_type::tensor<double, 1> ret = xt::zeros<double>({m_ndof});
144
145#pragma omp parallel for
146 for (size_t d = 0; d < m_nnu; ++d) {
147 ret(m_iiu(d)) = m_Auu(d);
148 }
149
150#pragma omp parallel for
151 for (size_t d = 0; d < m_nnp; ++d) {
152 ret(m_iip(d)) = m_App(d);
153 }
154
155 return ret;
156 }
157
163 {
164 return m_Auu;
165 }
166
172 {
173 return m_App;
174 }
175
181 {
182 return this->data();
183 }
184
185private:
186 template <class T>
187 void dot_nodevec_impl(const T& x, T& b) const
188 {
189 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
190 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
191
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) {
195
196 size_t d = m_part(m, i);
197
198 if (d < m_nnu) {
199 b(m, i) = m_Auu(d) * x(m, i);
200 }
201 else {
202 b(m, i) = m_App(d - m_nnu) * x(m, i);
203 }
204 }
205 }
206 }
207
208 template <class T>
209 void dot_dofval_impl(const T& x, T& b) const
210 {
211 GOOSEFEM_ASSERT(x.size() == m_ndof);
212 GOOSEFEM_ASSERT(b.size() == m_ndof);
213
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));
217 }
218
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));
222 }
223 }
224
225public:
232 array_type::tensor<double, 1>
234 {
235 array_type::tensor<double, 1> b_u = xt::empty<double>({m_nnu});
236 this->dot_u(x_u, x_p, b_u);
237 return b_u;
238 }
239
246 void dot_u(
250 {
251 UNUSED(x_p);
252
253 GOOSEFEM_ASSERT(x_u.size() == m_nnu);
254 GOOSEFEM_ASSERT(x_p.size() == m_nnp);
255 GOOSEFEM_ASSERT(b_u.size() == m_nnu);
256
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);
260 }
261 }
262
271 {
272 array_type::tensor<double, 1> b_p = xt::empty<double>({m_nnp});
273 this->dot_p(x_u, x_p, b_p);
274 return b_p;
275 }
276
283 void dot_p(
287 {
288 UNUSED(x_u);
289
290 GOOSEFEM_ASSERT(x_u.size() == m_nnu);
291 GOOSEFEM_ASSERT(x_p.size() == m_nnp);
292 GOOSEFEM_ASSERT(b_p.size() == m_nnp);
293
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);
297 }
298 }
299
300private:
301 template <class T>
302 void solve_nodevec_impl(const T& b, T& x)
303 {
304 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
305 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
306
307 this->factorize();
308
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);
314 }
315 }
316 }
317 }
318
319 template <class T>
320 void solve_dofval_impl(const T& b, T& x)
321 {
322 GOOSEFEM_ASSERT(b.size() == m_ndof);
323 GOOSEFEM_ASSERT(x.size() == m_ndof);
324
325 this->factorize();
326
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));
330 }
331 }
332
333public:
339 array_type::tensor<double, 1>
341 {
342 array_type::tensor<double, 1> x_u = xt::empty<double>({m_nnu});
343 this->solve_u(b_u, x_p, x_u);
344 return x_u;
345 }
346
356 {
357 UNUSED(x_p);
358
359 GOOSEFEM_ASSERT(b_u.size() == m_nnu);
360 GOOSEFEM_ASSERT(x_p.size() == m_nnp);
361 GOOSEFEM_ASSERT(x_u.size() == m_nnu);
362
363 this->factorize();
364
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);
368 }
369 }
370
371private:
372 template <class T>
373 void reaction_nodevec_impl(const T& x, T& b) const
374 {
375 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
376 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
377
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);
383 }
384 }
385 }
386 }
387
388 template <class T>
389 void reaction_dofval_impl(const T& x, T& b) const
390 {
391 GOOSEFEM_ASSERT(x.size() == m_ndof);
392 GOOSEFEM_ASSERT(b.size() == m_ndof);
393
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));
397 }
398 }
399
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
404 {
405 UNUSED(x_u);
406
407 GOOSEFEM_ASSERT(x_u.size() == m_nnu);
408 GOOSEFEM_ASSERT(x_p.size() == m_nnp);
409 GOOSEFEM_ASSERT(b_p.size() == m_nnp);
410
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);
414 }
415 }
416
417private:
418 // The diagonal matrix, and its inverse (re-used to solve different RHS)
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;
422
423 // Bookkeeping
424 array_type::tensor<size_t, 2> m_part; // DOF-numbers per node, renumbered [nnode, ndim]
425 array_type::tensor<size_t, 1> m_iiu; // DOF-numbers that are unknown [nnu]
426 array_type::tensor<size_t, 1> m_iip; // DOF-numbers that are prescribed [nnp]
427
428 // Dimensions
429 size_t m_nnu; // number of unknown DOFs
430 size_t m_nnp; // number of prescribed DOFs
431
432 // Compute inverse (automatically evaluated by "solve")
433 void factorize()
434 {
435 if (!m_changed) {
436 return;
437 }
438
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);
442 }
443
444 m_changed = false;
445 }
446};
447
448} // namespace GooseFEM
449
450#endif
Diagonal matrix.
Generic mesh operations.
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
size_t m_nelem
See nelem().
Definition: Matrix.h:203
size_t m_nne
See nne().
Definition: Matrix.h:204
bool m_changed
Signal changes to data.
Definition: Matrix.h:209
size_t m_nnode
See nnode().
Definition: Matrix.h:205
size_t m_ndim
See ndim().
Definition: Matrix.h:206
size_t m_ndof
See ndof().
Definition: Matrix.h:207
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.
Definition: Matrix.h:415
const array_type::tensor< size_t, 1 > & iip() const
Prescribed DOFs.
Definition: Matrix.h:472
Reorder to lowest possible index, in specific order.
Definition: Mesh.h:2387
#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
bool is_unique(const T &arg)
Returns true is a list is unique (has not duplicate items).
Definition: assertions.h:20