dune-pdelab  2.7-git
matrixhelpers.hh
Go to the documentation of this file.
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 
16 
17 namespace 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  {
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 
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<
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<
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<
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
p
const P & p
Definition: constraints.hh:148
orderingbase.hh
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
value
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
tags.hh
Dune::PDELab::ISTL::container_tag
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:234