dune-pdelab  2.4-dev
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  static_assert(static_cast<int>(RV::blocklevel) == static_cast<int>(CV::blocklevel),"Both vectors must have identical blocking depth");
57 
58  typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
59 
60  };
61 
62  template<typename RowOrdering, typename ColOrdering, typename SubPattern_ = void>
63  class Pattern
64  : public std::vector<std::unordered_map<std::size_t,SubPattern_> >
65  {
66 
67  public:
68 
69  typedef SubPattern_ SubPattern;
70 
71  template<typename RI, typename CI>
72  void add_link(const RI& ri, const CI& ci)
73  {
74  recursive_add_entry(ri.view(),ci.view());
75  }
76 
77  template<typename RI, typename CI>
78  void recursive_add_entry(const RI& ri, const CI& ci)
79  {
80  this->resize(_row_ordering.blockCount());
81  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()))));
82  r.first->second.recursive_add_entry(ri.back_popped(),ci.back_popped());
83  }
84 
85  Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
86  : _row_ordering(row_ordering)
87  , _col_ordering(col_ordering)
88  {}
89 
90  private:
91 
92  const RowOrdering& _row_ordering;
93  const ColOrdering& _col_ordering;
94 
95  };
96 
97  template<typename RowOrdering, typename ColOrdering>
98  class Pattern<RowOrdering,ColOrdering,void>
99  : public std::vector<std::unordered_set<std::size_t> >
100  {
101 
102  public:
103 
104  typedef void SubPattern;
105 
106  template<typename RI, typename CI>
107  void add_link(const RI& ri, const CI& ci)
108  {
109  recursive_add_entry(ri,ci);
110  }
111 
112  template<typename RI, typename CI>
113  void recursive_add_entry(const RI& ri, const CI& ci)
114  {
115  this->resize(_row_ordering.blockCount());
116  (*this)[ri.back()].insert(ci.back());
117  }
118 
119  Pattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
120  : _row_ordering(row_ordering)
121  , _col_ordering(col_ordering)
122  {}
123 
124  private:
125 
126  const RowOrdering& _row_ordering;
127  const ColOrdering& _col_ordering;
128 
129  };
130 
131  template<typename M, int blocklevel = M::blocklevel>
132  struct requires_pattern
133  {
134  static const bool value = requires_pattern<typename M::block_type,blocklevel-1>::value;
135  };
136 
137  template<typename M>
138  struct requires_pattern<M,0>
139  {
140  static const bool value = false;
141  };
142 
143  template<typename B, typename A, int blocklevel>
144  struct requires_pattern<Dune::BCRSMatrix<B,A>,blocklevel>
145  {
146  static const bool value = true;
147  };
148 
149  template<typename M, typename RowOrdering, typename ColOrdering, bool pattern>
150  struct _build_pattern_type
151  {
152  typedef void type;
153  };
154 
155  template<typename M, typename RowOrdering, typename ColOrdering>
156  struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
157  {
159  };
160 
161  template<typename M, typename GFSV, typename GFSU, typename Tag>
162  struct build_pattern_type
163  {
164 
165  typedef OrderingBase<
166  typename GFSV::Ordering::Traits::DOFIndex,
167  typename GFSV::Ordering::Traits::ContainerIndex
168  > RowOrdering;
169 
170  typedef OrderingBase<
171  typename GFSU::Ordering::Traits::DOFIndex,
172  typename GFSU::Ordering::Traits::ContainerIndex
173  > ColOrdering;
174 
176  };
177 
178  template<typename M, typename GFSV, typename GFSU>
179  struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
180  {
181  typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
182  };
183 
184 
185  template<typename RI, typename CI, typename Block>
186  typename Block::field_type&
187  access_matrix_element(tags::field_matrix_1_1, Block& b, const RI& ri, const CI& ci, int i, int j)
188  {
189  assert(i == -1);
190  assert(j == -1);
191  return b[0][0];
192  }
193 
194  template<typename RI, typename CI, typename Block>
195  typename Block::field_type&
196  access_matrix_element(tags::field_matrix_n_m, Block& b, const RI& ri, const CI& ci, int i, int j)
197  {
198  assert(i == 0);
199  assert(j == 0);
200  return b[ri[0]][ci[0]];
201  }
202 
203  template<typename RI, typename CI, typename Block>
204  typename Block::field_type&
205  access_matrix_element(tags::field_matrix_1_m, Block& b, const RI& ri, const CI& ci, int i, int j)
206  {
207  assert(i == -1);
208  assert(j == 0);
209  return b[0][ci[0]];
210  }
211 
212  template<typename RI, typename CI, typename Block>
213  typename Block::field_type&
214  access_matrix_element(tags::field_matrix_n_1, Block& b, const RI& ri, const CI& ci, int i, int j)
215  {
216  assert(i == 0);
217  assert(j == -1);
218  return b[ri[0]][0];
219  }
220 
221  template<typename RI, typename CI, typename Block>
222  typename Block::field_type&
223  access_matrix_element(tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
224  {
225  return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
226  }
227 
228 
229  template<typename RI, typename CI, typename Block>
230  const typename Block::field_type&
231  access_matrix_element(tags::field_matrix_1_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
232  {
233  assert(i == -1);
234  assert(j == -1);
235  return b[0][0];
236  }
237 
238  template<typename RI, typename CI, typename Block>
239  const typename Block::field_type&
240  access_matrix_element(tags::field_matrix_n_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
241  {
242  assert(i == 0);
243  assert(j == 0);
244  return b[ri[0]][ci[0]];
245  }
246 
247  template<typename RI, typename CI, typename Block>
248  const typename Block::field_type&
249  access_matrix_element(tags::field_matrix_1_m, const Block& b, const RI& ri, const CI& ci, int i, int j)
250  {
251  assert(i == -1);
252  assert(j == 0);
253  return b[0][ci[0]];
254  }
255 
256  template<typename RI, typename CI, typename Block>
257  const typename Block::field_type&
258  access_matrix_element(tags::field_matrix_n_1, const Block& b, const RI& ri, const CI& ci, int i, int j)
259  {
260  assert(i == 0);
261  assert(j == -1);
262  return b[ri[0]][0];
263  }
264 
265  template<typename RI, typename CI, typename Block>
266  const typename Block::field_type&
267  access_matrix_element(tags::bcrs_matrix, const Block& b, const RI& ri, const CI& ci, int i, int j)
268  {
269  return access_matrix_element(container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
270  }
271 
272 
273 
274  template<typename RI, typename Block>
275  void clear_matrix_row(tags::field_matrix_1_any, Block& b, const RI& ri, int i)
276  {
277  assert(i == -1);
278  b[0] = 0;
279  }
280 
281  template<typename RI, typename Block>
282  void clear_matrix_row(tags::field_matrix_n_any, Block& b, const RI& ri, int i)
283  {
284  assert(i == 0);
285  b[ri[0]] = 0;
286  }
287 
288  template<typename RI, typename Block>
289  void clear_matrix_row(tags::bcrs_matrix, Block& b, const RI& ri, int i)
290  {
291  typedef typename Block::ColIterator col_iterator_type;
292  const col_iterator_type end = b[ri[i]].end();
293  for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
294  clear_matrix_row(container_tag(*cit),*cit,ri,i-1);
295  }
296 
297 
298  template<typename T, typename RI, typename CI, typename Block>
299  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)
300  {
301  assert(i == -1);
302  assert(j == -1);
303  b[0][0] = v;
304  }
305 
306  template<typename T, typename RI, typename CI, typename Block>
307  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)
308  {
309  assert(i == 0);
310  assert(j == 0);
311  b[ri[0]][ci[0]] = v;
312  }
313 
314  template<typename T, typename RI, typename CI, typename Block>
315  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)
316  {
317  assert(i == -1);
318  assert(j == 0);
319  b[0][ci[0]] = v;
320  }
321 
322  template<typename T, typename RI, typename CI, typename Block>
323  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)
324  {
325  assert(i == 0);
326  assert(j == -1);
327  b[ri[0]][0] = v;
328  }
329 
330  template<typename T, typename RI, typename CI, typename Block>
331  void write_matrix_element_if_exists(const T& v, tags::bcrs_matrix, Block& b, const RI& ri, const CI& ci, int i, int j)
332  {
333  if (b.exists(ri[i],ci[j]))
334  write_matrix_element_if_exists(v,container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
335  }
336 
337 
338  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
339  typename enable_if<
342  >::type
343  allocate_matrix(const OrderingV& ordering_v,
344  const OrderingU& ordering_u,
345  const Pattern& p,
346  Container& c)
347  {
348  c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
349  c.setBuildMode(Container::random);
350 
351  for (std::size_t i = 0; i < c.N(); ++i)
352  c.setrowsize(i,p[i].size());
353  c.endrowsizes();
354 
355  for (std::size_t i = 0; i < c.N(); ++i)
356  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
357  c.addindex(i,cit->first);
358  c.endindices();
359 
360  for (std::size_t i = 0; i < c.N(); ++i)
361  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
362  {
363  allocate_matrix(ordering_v.childOrdering(i),
364  ordering_u.childOrdering(cit->first),
365  cit->second,
366  c[i][cit->first]);
367  }
368  }
369 
370  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
371  typename enable_if<
374  >::type
375  allocate_matrix(const OrderingV& ordering_v,
376  const OrderingU& ordering_u,
377  const Pattern& p,
378  Container& c)
379  {
380  for (std::size_t i = 0; i < c.N(); ++i)
381  for (typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
382  {
383  allocate_matrix(ordering_v.childOrdering(i),
384  ordering_u.childOrdering(cit->first),
385  cit->second,
386  c[i][cit->first]);
387  }
388  }
389 
390  template<typename OrderingV, typename OrderingU, typename Pattern, typename Container>
391  typename enable_if<
393  >::type
394  allocate_matrix(const OrderingV& ordering_v,
395  const OrderingU& ordering_u,
396  const Pattern& p,
397  Container& c)
398  {
399  c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),false);
400  c.setBuildMode(Container::random);
401 
402  for (std::size_t i = 0; i < c.N(); ++i)
403  c.setrowsize(i,p[i].size());
404  c.endrowsizes();
405 
406  for (std::size_t i = 0; i < c.N(); ++i)
407  for (typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
408  c.addindex(i,*cit);
409  c.endindices();
410  }
411 
412  } // namespace istl
413 
414 #endif // DOXYGEN
415 
416  } // namespace PDELab
417 } // namespace Dune
418 
419 #endif // DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:246
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: adaptivity.hh:27
const P & p
Definition: constraints.hh:146