DUNE PDELab (git)

matrixhelpers.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
4#define DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
5
6#include<utility>
7#include<vector>
8#include <unordered_map>
9#include <unordered_set>
10
11#include <dune/common/fmatrix.hh>
12#include <dune/istl/bcrsmatrix.hh>
13
14#include <dune/pdelab/ordering/orderingbase.hh>
15#include <dune/pdelab/backend/istl/tags.hh>
16
17namespace Dune {
18 namespace PDELab {
19
20#ifndef DOXYGEN
21
22 namespace ISTL {
23
24 template<typename RV, typename CV, typename block_type>
25 struct matrix_for_vectors;
26
27 template<typename B1, typename A1, typename B2, typename A2, typename block_type>
28 struct matrix_for_vectors<Dune::BlockVector<B1,A1>,Dune::BlockVector<B2,A2>,block_type>
29 {
30 typedef Dune::BCRSMatrix<block_type> type;
31 };
32
33 template<typename B1, int n1, typename B2, int n2, typename block_type>
34 struct matrix_for_vectors<Dune::FieldVector<B1,n1>,Dune::FieldVector<B2,n2>,block_type>
35 {
36 typedef Dune::FieldMatrix<block_type,n1,n2> type;
37 };
38
39 template<typename E, typename RV, typename CV, std::size_t blocklevel>
40 struct recursive_build_matrix_type
41 {
42 typedef typename matrix_for_vectors<RV,CV,typename recursive_build_matrix_type<E,typename RV::block_type,typename CV::block_type,blocklevel-1>::type>::type type;
43 };
44
45 template<typename E, typename RV, typename CV>
46 struct recursive_build_matrix_type<E,RV,CV,1>
47 {
48 typedef Dune::FieldMatrix<E,RV::dimension,CV::dimension> type;
49 };
50
51
52 template<typename E, typename RV, typename CV>
53 struct build_matrix_type
54 {
55
56#if DUNE_VERSION_LT(DUNE_ISTL,2,8)
57 static_assert(static_cast<int>(RV::blocklevel) == static_cast<int>(CV::blocklevel),"Both vectors must have identical blocking depth");
58
59 typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
60#else
61 static_assert(static_cast<int>(blockLevel<RV>()) == static_cast<int>(blockLevel<CV>()),"Both vectors must have identical blocking depth");
62
63 typedef typename recursive_build_matrix_type<E,RV,CV,blockLevel<RV>()>::type type;
64#endif
65
66 };
67
68 template<typename RowOrdering, typename ColOrdering, typename SubPattern_ = void>
69 class Pattern
70 : public std::vector<std::unordered_map<std::size_t,SubPattern_> >
71 {
72
73 public:
74
75 typedef SubPattern_ SubPattern;
76
77 template<typename RI, typename CI>
78 void add_link(const RI& ri, const CI& ci)
79 {
80 recursive_add_entry(ri.view(),ci.view());
81 }
82
83 template<typename RI, typename CI>
84 void recursive_add_entry(const RI& ri, const CI& ci)
85 {
86 this->resize(_row_ordering.blockCount());
87 std::pair<typename std::unordered_map<std::size_t,SubPattern>::iterator,bool> r = (*this)[ri.back()].insert(make_pair(ci.back(),SubPattern(_row_ordering.childOrdering(ri.back()),_col_ordering.childOrdering(ci.back()))));
88 r.first->second.recursive_add_entry(ri.back_popped(),ci.back_popped());
89 }
90
91 Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
92 : _row_ordering(row_ordering)
93 , _col_ordering(col_ordering)
94 {}
95
96 private:
97
98 const RowOrdering& _row_ordering;
99 const ColOrdering& _col_ordering;
100
101 };
102
103 template<typename RowOrdering, typename ColOrdering>
104 class Pattern<RowOrdering,ColOrdering,void>
105 : public std::vector<std::unordered_set<std::size_t> >
106 {
107
108 public:
109
110 typedef void SubPattern;
111
112 template<typename RI, typename CI>
113 void add_link(const RI& ri, const CI& ci)
114 {
115 recursive_add_entry(ri,ci);
116 }
117
118 template<typename RI, typename CI>
119 void recursive_add_entry(const RI& ri, const CI& ci)
120 {
121 this->resize(_row_ordering.blockCount());
122 (*this)[ri.back()].insert(ci.back());
123 }
124
125 Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
126 : _row_ordering(row_ordering)
127 , _col_ordering(col_ordering)
128 {}
129
130 private:
131
132 const RowOrdering& _row_ordering;
133 const ColOrdering& _col_ordering;
134
135 };
136
137#if DUNE_VERSION_LT(DUNE_ISTL,2,8)
138 template<typename M, int blocklevel = M::blocklevel>
139#else
140 template<typename M, int blocklevel = blockLevel<M>()>
141#endif
142 struct requires_pattern
143 {
144 static const bool value = requires_pattern<typename M::block_type,blocklevel-1>::value;
145 };
146
147 template<typename M>
148 struct requires_pattern<M,0>
149 {
150 static const bool value = false;
151 };
152
153 template<typename B, typename A, int blocklevel>
154 struct requires_pattern<Dune::BCRSMatrix<B,A>,blocklevel>
155 {
156 static const bool value = true;
157 };
158
159 template<typename M, typename RowOrdering, typename ColOrdering, bool pattern>
160 struct _build_pattern_type
161 {
162 typedef void type;
163 };
164
165 template<typename M, typename RowOrdering, typename ColOrdering>
166 struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
167 {
168 typedef Pattern<RowOrdering,ColOrdering,typename _build_pattern_type<typename M::block_type,RowOrdering,ColOrdering,requires_pattern<typename M::block_type>::value>::type> type;
169 };
170
171 template<typename M, typename GFSV, typename GFSU, typename Tag>
172 struct build_pattern_type
173 {
174
175 typedef OrderingBase<
176 typename GFSV::Ordering::Traits::DOFIndex,
177 typename GFSV::Ordering::Traits::ContainerIndex
178 > RowOrdering;
179
180 typedef OrderingBase<
181 typename GFSU::Ordering::Traits::DOFIndex,
182 typename GFSU::Ordering::Traits::ContainerIndex
183 > ColOrdering;
184
185 typedef typename _build_pattern_type<M,RowOrdering,ColOrdering,requires_pattern<M>::value>::type type;
186 };
187
188 template<typename M, typename GFSV, typename GFSU>
189 struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
190 {
191 typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
192 };
193
194
195 template<typename RI, typename CI, typename Block>
196 typename Block::field_type&
197 access_matrix_element(tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
198 {
199 assert(i == -1);
200 assert(j == -1);
201 return b[0][0];
202 }
203
204 template<typename RI, typename CI, typename Block>
205 typename Block::field_type&
206 access_matrix_element(tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
207 {
208 assert(i == 0);
209 assert(j == 0);
210 return b[ri[0]][ci[0]];
211 }
212
213 template<typename RI, typename CI, typename Block>
214 typename Block::field_type&
215 access_matrix_element(tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
216 {
217 assert(i == -1);
218 assert(j == 0);
219 return b[0][ci[0]];
220 }
221
222 template<typename RI, typename CI, typename Block>
223 typename Block::field_type&
224 access_matrix_element(tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
225 {
226 assert(i == 0);
227 assert(j == -1);
228 return b[ri[0]][0];
229 }
230
231 template<typename RI, typename CI, typename Block>
232 typename Block::field_type&
233 access_matrix_element(tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
234 {
235 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
236 }
237
238
239 template<typename RI, typename CI, typename Block>
240 const typename Block::field_type&
241 access_matrix_element(tags::field_matrix_1_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
242 {
243 assert(i == -1);
244 assert(j == -1);
245 return b[0][0];
246 }
247
248 template<typename RI, typename CI, typename Block>
249 const typename Block::field_type&
250 access_matrix_element(tags::field_matrix_n_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
251 {
252 assert(i == 0);
253 assert(j == 0);
254 return b[ri[0]][ci[0]];
255 }
256
257 template<typename RI, typename CI, typename Block>
258 const typename Block::field_type&
259 access_matrix_element(tags::field_matrix_1_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
260 {
261 assert(i == -1);
262 assert(j == 0);
263 return b[0][ci[0]];
264 }
265
266 template<typename RI, typename CI, typename Block>
267 const typename Block::field_type&
268 access_matrix_element(tags::field_matrix_n_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
269 {
270 assert(i == 0);
271 assert(j == -1);
272 return b[ri[0]][0];
273 }
274
275 template<typename RI, typename CI, typename Block>
276 const typename Block::field_type&
277 access_matrix_element(tags::bcrs_matrix, const Block& b, const RI& ri, const CI& ci, int i, int j)
278 {
279 return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
280 }
281
282
283
284 template<typename RI, typename Block>
285 void clear_matrix_row(tags::field_matrix_1_any, Block& b, const RI& ri, int i)
286 {
287 assert(i == -1);
288 b[0] = 0;
289 }
290
291 template<typename RI, typename Block>
292 void clear_matrix_row(tags::field_matrix_n_any, Block& b, const RI& ri, int i)
293 {
294 assert(i == 0);
295 b[ri[0]] = 0;
296 }
297
298 template<typename RI, typename Block>
299 void clear_matrix_row(tags::bcrs_matrix, Block& b, const RI& ri, int i)
300 {
301 typedef typename Block::ColIterator col_iterator_type;
302 const col_iterator_type end = b[ri[i]].end();
303 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
304 clear_matrix_row(container_tag(*cit),*cit,ri,i-1);
305 }
306
307
308 template<typename RI, typename Block>
309 void clear_matrix_row_block(tags::field_matrix_1_1, Block& b, const RI& ri, int i)
310 {
311 assert(i == -1);
312 b = 0;
313 }
314
315 template<typename RI, typename Block>
316 void clear_matrix_row_block(tags::field_matrix_1_any, Block& b, const RI& ri, int i)
317 {
318 DUNE_THROW(Dune::Exception,"Should never get here!");
319 }
320
321 template<typename RI, typename Block>
322 void clear_matrix_row_block(tags::field_matrix_n_any, Block& b, const RI& ri, int i)
323 {
324 assert(i == 0);
325 b = 0;
326 }
327
328 template<typename RI, typename Block>
329 void clear_matrix_row_block(tags::bcrs_matrix, Block& b, const RI& ri, int i)
330 {
331 typedef typename Block::ColIterator col_iterator_type;
332 const col_iterator_type end = b[ri[i]].end();
333 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
334 clear_matrix_row_block(container_tag(*cit),*cit,ri,i-1);
335 }
336
337
338
339 template<typename T, typename RI, typename CI, typename Block>
340 void write_matrix_element_if_exists(const T& v, tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
341 {
342 assert(i == -1);
343 assert(j == -1);
344 b[0][0] = v;
345 }
346
347 template<typename T, typename RI, typename CI, typename Block>
348 void write_matrix_element_if_exists(const T& v, tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
349 {
350 assert(i == 0);
351 assert(j == 0);
352 b[ri[0]][ci[0]] = v;
353 }
354
355 template<typename T, typename RI, typename CI, typename Block>
356 void write_matrix_element_if_exists(const T& v, tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
357 {
358 assert(i == -1);
359 assert(j == 0);
360 b[0][ci[0]] = v;
361 }
362
363 template<typename T, typename RI, typename CI, typename Block>
364 void write_matrix_element_if_exists(const T& v, tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
365 {
366 assert(i == 0);
367 assert(j == -1);
368 b[ri[0]][0] = v;
369 }
370
371 template<typename T, typename RI, typename CI, typename Block>
372 void write_matrix_element_if_exists(const T& v, tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
373 {
374 if (b.exists(ri[i],ci[j]))
375 write_matrix_element_if_exists(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
376 }
377
378
379
380
381 template<typename T, typename RI, typename CI, typename Block>
382 void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
383 {
384 assert(i == -1);
385 assert(j == -1);
386 b[0][0] = v;
387 }
388
389 template<typename T, typename RI, typename CI, typename Block>
390 void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
391 {
392 assert(i == 0);
393 assert(j == 0);
394 for (std::size_t i = 0; i < b.rows; ++i)
395 b[i][i] = v;
396 }
397
398 template<typename T, typename RI, typename CI, typename Block>
399 void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
400 {
401 DUNE_THROW(Dune::Exception,"Should never get here!");
402 }
403
404 template<typename T, typename RI, typename CI, typename Block>
405 void write_matrix_element_if_exists_to_block(const T& v, tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
406 {
407 DUNE_THROW(Dune::Exception,"Should never get here!");
408 }
409
410 template<typename T, typename RI, typename CI, typename Block>
411 void write_matrix_element_if_exists_to_block(const T& v, tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
412 {
413 if (b.exists(ri[i],ci[j]))
414 write_matrix_element_if_exists_to_block(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
415 }
416
417
418 template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
419 typename std::enable_if<
420 !std::is_same<typename Pattern::SubPattern,void>::value &&
421 requires_pattern<Container>::value
422 >::type
423 allocate_matrix(const OrderingV& ordering_v,
424 const OrderingU& ordering_u,
425 const Pattern& p,
426 Container& c)
427 {
428 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
429 c.setBuildMode(Container::random);
430
431 for (std::size_t i = 0; i < c.N(); ++i)
432 c.setrowsize(i,p[i].size());
433 c.endrowsizes();
434
435 for (std::size_t i = 0; i < c.N(); ++i)
436 for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
437 c.addindex(i,cit->first);
438 c.endindices();
439
440 for (std::size_t i = 0; i < c.N(); ++i)
441 for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
442 {
443 allocate_matrix(ordering_v.childOrdering(i),
444 ordering_u.childOrdering(cit->first),
445 cit->second,
446 c[i][cit->first]);
447 }
448 }
449
450 template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
451 typename std::enable_if<
452 !std::is_same<typename Pattern::SubPattern,void>::value &&
453 !requires_pattern<Container>::value
454 >::type
455 allocate_matrix(const OrderingV& ordering_v,
456 const OrderingU& ordering_u,
457 const Pattern& p,
458 Container& c)
459 {
460 for (std::size_t i = 0; i < c.N(); ++i)
461 for (typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
462 {
463 allocate_matrix(ordering_v.childOrdering(i),
464 ordering_u.childOrdering(cit->first),
465 cit->second,
466 c[i][cit->first]);
467 }
468 }
469
470 template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
471 typename std::enable_if<
472 std::is_same<typename Pattern::SubPattern,void>::value
473 >::type
474 allocate_matrix(const OrderingV& ordering_v,
475 const OrderingU& ordering_u,
476 const Pattern& p,
477 Container& c)
478 {
479 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
480 c.setBuildMode(Container::random);
481
482 for (std::size_t i = 0; i < c.N(); ++i)
483 c.setrowsize(i,p[i].size());
484 c.endrowsizes();
485
486 for (std::size_t i = 0; i < c.N(); ++i)
487 for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
488 c.addindex(i,*cit);
489 c.endindices();
490 }
491
492 } // namespace ISTL
493
494#endif // DOXYGEN
495
496 } // namespace PDELab
497} // namespace Dune
498
499#endif // DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
Dune::BlockVector
A vector of blocks with memory management.
Definition: bvector.hh:392
Dune::Exception
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:117
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:91
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)