FrictionQPotFEM 0.23.3
Loading...
Searching...
No Matches
MatrixPartitioned.h
Go to the documentation of this file.
1
11#ifndef GOOSEFEM_MATRIXPARTITIONED_H
12#define GOOSEFEM_MATRIXPARTITIONED_H
13
14#include "Matrix.h"
15#include "config.h"
16
17#include <Eigen/Eigen>
18#include <Eigen/Sparse>
19#include <Eigen/SparseCholesky>
20
21namespace GooseFEM {
22
23// forward declaration
24template <class>
25class MatrixPartitionedSolver;
26
34class MatrixPartitioned : public MatrixPartitionedBase<MatrixPartitioned> {
35private:
38
39protected:
40 Eigen::SparseMatrix<double> m_Auu;
41 Eigen::SparseMatrix<double> m_Aup;
42 Eigen::SparseMatrix<double> m_Apu;
43 Eigen::SparseMatrix<double> m_App;
44
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;
49
59
70
74 template <class>
76
77public:
78 MatrixPartitioned() = default;
79
91 {
92 m_conn = conn;
93 m_dofs = dofs;
94 m_iip = iip;
95 m_nelem = m_conn.shape(0);
96 m_nne = m_conn.shape(1);
97 m_nnode = m_dofs.shape(0);
98 m_ndim = m_dofs.shape(1);
99 m_ndof = xt::amax(m_dofs)() + 1;
100
102 GOOSEFEM_ASSERT(xt::amax(conn)() + 1 <= m_nnode);
103 GOOSEFEM_ASSERT(xt::amax(iip)() <= xt::amax(dofs)());
105
106 m_iiu = xt::setdiff1d(dofs, iip);
107 m_nnp = m_iip.size();
108 m_nnu = m_iiu.size();
109
111 m_part1d = Mesh::Reorder({m_iiu, m_iip}).apply(xt::eval(xt::arange<size_t>(m_ndof)));
112 m_Tuu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
113 m_Tup.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
114 m_Tpu.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
115 m_Tpp.reserve(m_nelem * m_nne * m_ndim * m_nne * m_ndim);
116 m_Auu.resize(m_nnu, m_nnu);
117 m_Aup.resize(m_nnu, m_nnp);
118 m_Apu.resize(m_nnp, m_nnu);
119 m_App.resize(m_nnp, m_nnp);
120 }
121
125 const Eigen::SparseMatrix<double>& data_uu() const
126 {
127 return m_Auu;
128 }
129
133 const Eigen::SparseMatrix<double>& data_up() const
134 {
135 return m_Aup;
136 }
137
141 const Eigen::SparseMatrix<double>& data_pu() const
142 {
143 return m_Apu;
144 }
145
149 const Eigen::SparseMatrix<double>& data_pp() const
150 {
151 return m_App;
152 }
153
154private:
155 template <class T>
156 void assemble_impl(const T& elemmat)
157 {
158 GOOSEFEM_ASSERT(xt::has_shape(elemmat, {m_nelem, m_nne * m_ndim, m_nne * m_ndim}));
159
160 m_Tuu.clear();
161 m_Tup.clear();
162 m_Tpu.clear();
163 m_Tpp.clear();
164
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) {
168
169 size_t di = m_part(m_conn(e, m), i);
170
171 for (size_t n = 0; n < m_nne; ++n) {
172 for (size_t j = 0; j < m_ndim; ++j) {
173
174 size_t dj = m_part(m_conn(e, n), j);
175
176 if (di < m_nnu && dj < m_nnu) {
177 m_Tuu.push_back(Eigen::Triplet<double>(
178 di, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
179 }
180 else if (di < m_nnu) {
181 m_Tup.push_back(Eigen::Triplet<double>(
182 di, dj - m_nnu, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
183 }
184 else if (dj < m_nnu) {
185 m_Tpu.push_back(Eigen::Triplet<double>(
186 di - m_nnu, dj, elemmat(e, m * m_ndim + i, n * m_ndim + j)));
187 }
188 else {
189 m_Tpp.push_back(Eigen::Triplet<double>(
190 di - m_nnu,
191 dj - m_nnu,
192 elemmat(e, m * m_ndim + i, n * m_ndim + j)));
193 }
194 }
195 }
196 }
197 }
198 }
199
200 m_Auu.setFromTriplets(m_Tuu.begin(), m_Tuu.end());
201 m_Aup.setFromTriplets(m_Tup.begin(), m_Tup.end());
202 m_Apu.setFromTriplets(m_Tpu.begin(), m_Tpu.end());
203 m_App.setFromTriplets(m_Tpp.begin(), m_Tpp.end());
204 m_changed = true;
205 }
206
207public:
215 void
218 const array_type::tensor<double, 2>& matrix)
219 {
220 GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
221 GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
222 GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
223 GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
224
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;
229
230 for (size_t i = 0; i < rows.size(); ++i) {
231 for (size_t j = 0; j < cols.size(); ++j) {
232 size_t di = m_part1d(rows(i));
233 size_t dj = m_part1d(cols(j));
234 double v = matrix(i, j);
235
236 if (di < m_nnu && dj < m_nnu) {
237 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
238 }
239 else if (di < m_nnu) {
240 Tup.push_back(Eigen::Triplet<double>(di, dj - m_nnu, v));
241 }
242 else if (dj < m_nnu) {
243 Tpu.push_back(Eigen::Triplet<double>(di - m_nnu, dj, v));
244 }
245 else {
246 Tpp.push_back(Eigen::Triplet<double>(di - m_nnu, dj - m_nnu, v));
247 }
248 }
249 }
250
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());
255 m_changed = true;
256 }
257
265 void
268 const array_type::tensor<double, 2>& matrix)
269 {
270 GOOSEFEM_ASSERT(rows.size() == matrix.shape(0));
271 GOOSEFEM_ASSERT(cols.size() == matrix.shape(1));
272 GOOSEFEM_ASSERT(xt::amax(cols)() < m_ndof);
273 GOOSEFEM_ASSERT(xt::amax(rows)() < m_ndof);
274
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;
279
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);
284
285 for (size_t i = 0; i < rows.size(); ++i) {
286 for (size_t j = 0; j < cols.size(); ++j) {
287 size_t di = m_part1d(rows(i));
288 size_t dj = m_part1d(cols(j));
289 double v = matrix(i, j);
290
291 if (di < m_nnu && dj < m_nnu) {
292 Tuu.push_back(Eigen::Triplet<double>(di, dj, v));
293 }
294 else if (di < m_nnu) {
295 Tup.push_back(Eigen::Triplet<double>(di, dj - m_nnu, v));
296 }
297 else if (dj < m_nnu) {
298 Tpu.push_back(Eigen::Triplet<double>(di - m_nnu, dj, v));
299 }
300 else {
301 Tpp.push_back(Eigen::Triplet<double>(di - m_nnu, dj - m_nnu, v));
302 }
303 }
304 }
305
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());
310 m_Auu += Auu;
311 m_Aup += Aup;
312 m_Apu += Apu;
313 m_App += App;
314 m_changed = true;
315 }
316
317private:
318 template <class T>
319 void todense_impl(T& ret) const
320 {
321 ret.fill(0.0);
322
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();
326 }
327 }
328
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();
332 }
333 }
334
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();
338 }
339 }
340
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();
344 }
345 }
346 }
347
348 template <class T>
349 void dot_nodevec_impl(const T& x, T& b) const
350 {
351 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
352 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
353
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;
358
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) {
362 if (m_part(m, i) < m_nnu) {
363 b(m, i) = B_u(m_part(m, i));
364 }
365 else {
366 b(m, i) = B_p(m_part(m, i) - m_nnu);
367 }
368 }
369 }
370 }
371
372 template <class T>
373 void dot_dofval_impl(const T& x, T& b) const
374 {
375 GOOSEFEM_ASSERT(b.size() == m_ndof);
376 GOOSEFEM_ASSERT(x.size() == m_ndof);
377
378 Eigen::VectorXd X_u = this->AsDofs_u(x);
379 Eigen::VectorXd X_p = this->AsDofs_p(x);
380
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;
383
384#pragma omp parallel for
385 for (size_t d = 0; d < m_nnu; ++d) {
386 b(m_iiu(d)) = B_u(d);
387 }
388
389#pragma omp parallel for
390 for (size_t d = 0; d < m_nnp; ++d) {
391 b(m_iip(d)) = B_p(d);
392 }
393 }
394
395 template <class T>
396 void reaction_nodevec_impl(const T& x, T& b) const
397 {
398 GOOSEFEM_ASSERT(xt::has_shape(x, {m_nnode, m_ndim}));
399 GOOSEFEM_ASSERT(xt::has_shape(b, {m_nnode, m_ndim}));
400
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;
404
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) {
408 if (m_part(m, i) >= m_nnu) {
409 b(m, i) = B_p(m_part(m, i) - m_nnu);
410 }
411 }
412 }
413 }
414
415 template <class T>
416 void reaction_dofval_impl(const T& x, T& b) const
417 {
418 GOOSEFEM_ASSERT(x.size() == m_ndof);
419 GOOSEFEM_ASSERT(b.size() == m_ndof);
420
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;
424
425#pragma omp parallel for
426 for (size_t d = 0; d < m_nnp; ++d) {
427 b(m_iip(d)) = B_p(d);
428 }
429 }
430
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
435 {
436 GOOSEFEM_ASSERT(x_u.size() == m_nnu);
437 GOOSEFEM_ASSERT(x_p.size() == m_nnp);
438 GOOSEFEM_ASSERT(b_p.size() == m_nnp);
439
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());
443 }
444
445private:
446 // Convert arrays (Eigen version of VectorPartitioned, which contains public functions)
447 Eigen::VectorXd AsDofs_u(const array_type::tensor<double, 1>& dofval) const
448 {
449 GOOSEFEM_ASSERT(dofval.size() == m_ndof);
450
451 Eigen::VectorXd dofval_u(m_nnu, 1);
452
453#pragma omp parallel for
454 for (size_t d = 0; d < m_nnu; ++d) {
455 dofval_u(d) = dofval(m_iiu(d));
456 }
457
458 return dofval_u;
459 }
460
461 Eigen::VectorXd AsDofs_u(const array_type::tensor<double, 2>& nodevec) const
462 {
463 GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
464
465 Eigen::VectorXd dofval_u = Eigen::VectorXd::Zero(m_nnu, 1);
466
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) {
470 if (m_part(m, i) < m_nnu) {
471 dofval_u(m_part(m, i)) = nodevec(m, i);
472 }
473 }
474 }
475
476 return dofval_u;
477 }
478
479 Eigen::VectorXd AsDofs_p(const array_type::tensor<double, 1>& dofval) const
480 {
481 GOOSEFEM_ASSERT(dofval.size() == m_ndof);
482
483 Eigen::VectorXd dofval_p(m_nnp, 1);
484
485#pragma omp parallel for
486 for (size_t d = 0; d < m_nnp; ++d) {
487 dofval_p(d) = dofval(m_iip(d));
488 }
489
490 return dofval_p;
491 }
492
493 Eigen::VectorXd AsDofs_p(const array_type::tensor<double, 2>& nodevec) const
494 {
495 GOOSEFEM_ASSERT(xt::has_shape(nodevec, {m_nnode, m_ndim}));
496
497 Eigen::VectorXd dofval_p = Eigen::VectorXd::Zero(m_nnp, 1);
498
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) {
502 if (m_part(m, i) >= m_nnu) {
503 dofval_p(m_part(m, i) - m_nnu) = nodevec(m, i);
504 }
505 }
506 }
507
508 return dofval_p;
509 }
510};
511
519template <class Solver = Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>>
521 : public MatrixSolverBase<MatrixPartitionedSolver<Solver>>,
522 public MatrixSolverPartitionedBase<MatrixPartitionedSolver<Solver>> {
523private:
526
527public:
528 MatrixPartitionedSolver() = default;
529
538 template <class M>
540 M& A,
543 {
544 GOOSEFEM_ASSERT(xt::has_shape(b_u, {A.nnu()}));
545 GOOSEFEM_ASSERT(xt::has_shape(x_p, {A.nnp()}));
546 array_type::tensor<double, 1> x_u = xt::empty_like(b_u);
547 this->solve_u_impl(A, b_u, x_p, x_u);
548 return x_u;
549 }
550
560 template <class M>
562 M& A,
566 {
567 GOOSEFEM_ASSERT(xt::has_shape(b_u, {A.nnu()}));
568 GOOSEFEM_ASSERT(xt::has_shape(x_p, {A.nnp()}));
569 GOOSEFEM_ASSERT(xt::has_shape(x_u, {A.nnu()}));
570 this->solve_u_impl(A, b_u, x_p, x_u);
571 }
572
573private:
574 template <class T>
575 void solve_nodevec_impl(MatrixPartitioned& A, const T& b, T& x)
576 {
577 this->factorize(A);
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));
581
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) {
585 if (A.m_part(m, i) < A.m_nnu) {
586 x(m, i) = X_u(A.m_part(m, i));
587 }
588 }
589 }
590 }
591
592 template <class T>
593 void solve_dofval_impl(MatrixPartitioned& A, const T& b, T& x)
594 {
595 this->factorize(A);
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));
599
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);
603 }
604 }
605
606 template <class T>
607 void solve_u_impl(MatrixPartitioned& A, const T& b_u, const T& x_p, T& x_u)
608 {
609 this->factorize(A);
610
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())));
615 }
616
617private:
618 Solver m_solver;
619 bool m_factor = true;
620
624 void factorize(MatrixPartitioned& A)
625 {
626 if (!A.m_changed && !m_factor) {
627 return;
628 }
629 m_solver.compute(A.m_Auu);
630 m_factor = false;
631 A.m_changed = false;
632 }
633};
634
635} // namespace GooseFEM
636
637#endif
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.
Definition: Matrix.h:415
array_type::tensor< size_t, 1 > m_iiu
See iiu()
Definition: Matrix.h:417
array_type::tensor< size_t, 1 > m_iip
See iip()
Definition: Matrix.h:418
const array_type::tensor< size_t, 1 > & iip() const
Prescribed DOFs.
Definition: Matrix.h:472
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.
Definition: Matrix.h:26
CRTP base class for a extra functions for a partitioned solver class.
Definition: Matrix.h:136
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
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