3 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH 4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH 8 #include <unordered_map> 9 #include <unordered_set> 11 #include <dune/common/fmatrix.hh> 12 #include <dune/istl/bcrsmatrix.hh> 24 template<
typename RV,
typename CV,
typename block_type>
25 struct matrix_for_vectors;
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>
30 typedef Dune::BCRSMatrix<block_type> type;
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>
36 typedef Dune::FieldMatrix<block_type,n1,n2> type;
39 template<
typename E,
typename RV,
typename CV, std::
size_t blocklevel>
40 struct recursive_build_matrix_type
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;
45 template<
typename E,
typename RV,
typename CV>
46 struct recursive_build_matrix_type<E,RV,CV,1>
48 typedef Dune::FieldMatrix<E,RV::dimension,CV::dimension> type;
52 template<
typename E,
typename RV,
typename CV>
53 struct build_matrix_type
56 static_assert(static_cast<int>(RV::blocklevel) == static_cast<int>(CV::blocklevel),
"Both vectors must have identical blocking depth");
58 typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
62 template<
typename RowOrdering,
typename ColOrdering,
typename SubPattern_ =
void>
64 :
public std::vector<std::unordered_map<std::size_t,SubPattern_> >
69 typedef SubPattern_ SubPattern;
71 template<
typename RI,
typename CI>
72 void add_link(
const RI& ri,
const CI& ci)
74 recursive_add_entry(ri.view(),ci.view());
77 template<
typename RI,
typename CI>
78 void recursive_add_entry(
const RI& ri,
const CI& ci)
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());
85 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
86 : _row_ordering(row_ordering)
87 , _col_ordering(col_ordering)
92 const RowOrdering& _row_ordering;
93 const ColOrdering& _col_ordering;
97 template<
typename RowOrdering,
typename ColOrdering>
98 class Pattern<RowOrdering,ColOrdering,void>
99 :
public std::vector<std::unordered_set<std::size_t> >
104 typedef void SubPattern;
106 template<
typename RI,
typename CI>
107 void add_link(
const RI& ri,
const CI& ci)
109 recursive_add_entry(ri,ci);
112 template<
typename RI,
typename CI>
113 void recursive_add_entry(
const RI& ri,
const CI& ci)
115 this->resize(_row_ordering.blockCount());
116 (*this)[ri.back()].insert(ci.back());
119 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
120 : _row_ordering(row_ordering)
121 , _col_ordering(col_ordering)
126 const RowOrdering& _row_ordering;
127 const ColOrdering& _col_ordering;
131 template<
typename M,
int blocklevel = M::blocklevel>
132 struct requires_pattern
134 static const bool value = requires_pattern<
typename M::block_type,blocklevel-1>
::value;
138 struct requires_pattern<M,0>
140 static const bool value =
false;
143 template<
typename B,
typename A,
int blocklevel>
144 struct requires_pattern<
Dune::BCRSMatrix<B,A>,blocklevel>
146 static const bool value =
true;
149 template<
typename M,
typename RowOrdering,
typename ColOrdering,
bool pattern>
150 struct _build_pattern_type
155 template<
typename M,
typename RowOrdering,
typename ColOrdering>
156 struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
161 template<
typename M,
typename GFSV,
typename GFSU,
typename Tag>
162 struct build_pattern_type
165 typedef OrderingBase<
166 typename GFSV::Ordering::Traits::DOFIndex,
167 typename GFSV::Ordering::Traits::ContainerIndex
170 typedef OrderingBase<
171 typename GFSU::Ordering::Traits::DOFIndex,
172 typename GFSU::Ordering::Traits::ContainerIndex
178 template<
typename M,
typename GFSV,
typename GFSU>
179 struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
181 typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
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)
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)
200 return b[ri[0]][ci[0]];
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)
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)
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)
225 return access_matrix_element(
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
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)
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)
244 return b[ri[0]][ci[0]];
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)
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)
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)
269 return access_matrix_element(
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
274 template<
typename RI,
typename Block>
275 void clear_matrix_row(tags::field_matrix_1_any, Block& b,
const RI& ri,
int i)
281 template<
typename RI,
typename Block>
282 void clear_matrix_row(tags::field_matrix_n_any, Block& b,
const RI& ri,
int i)
288 template<
typename RI,
typename Block>
289 void clear_matrix_row(tags::bcrs_matrix, Block& b,
const RI& ri,
int i)
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)
298 template<
typename RI,
typename Block>
299 void clear_matrix_row_block(tags::field_matrix_1_1, Block& b,
const RI& ri,
int i)
305 template<
typename RI,
typename Block>
306 void clear_matrix_row_block(tags::field_matrix_1_any, Block& b,
const RI& ri,
int i)
308 DUNE_THROW(Dune::Exception,
"Should never get here!");
311 template<
typename RI,
typename Block>
312 void clear_matrix_row_block(tags::field_matrix_n_any, Block& b,
const RI& ri,
int i)
318 template<
typename RI,
typename Block>
319 void clear_matrix_row_block(tags::bcrs_matrix, Block& b,
const RI& ri,
int i)
321 typedef typename Block::ColIterator col_iterator_type;
322 const col_iterator_type end = b[ri[i]].end();
323 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
329 template<
typename T,
typename RI,
typename CI,
typename Block>
330 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)
337 template<
typename T,
typename RI,
typename CI,
typename Block>
338 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)
345 template<
typename T,
typename RI,
typename CI,
typename Block>
346 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)
353 template<
typename T,
typename RI,
typename CI,
typename Block>
354 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)
361 template<
typename T,
typename RI,
typename CI,
typename Block>
362 void write_matrix_element_if_exists(
const T& v, tags::bcrs_matrix, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
364 if (b.exists(ri[i],ci[j]))
365 write_matrix_element_if_exists(v,
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
371 template<
typename T,
typename RI,
typename CI,
typename Block>
372 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)
379 template<
typename T,
typename RI,
typename CI,
typename Block>
380 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)
384 for (std::size_t i = 0; i < b.rows; ++i)
388 template<
typename T,
typename RI,
typename CI,
typename Block>
389 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)
391 DUNE_THROW(Dune::Exception,
"Should never get here!");
394 template<
typename T,
typename RI,
typename CI,
typename Block>
395 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)
397 DUNE_THROW(Dune::Exception,
"Should never get here!");
400 template<
typename T,
typename RI,
typename CI,
typename Block>
401 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)
403 if (b.exists(ri[i],ci[j]))
404 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);
408 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
409 typename std::enable_if<
413 allocate_matrix(
const OrderingV& ordering_v,
414 const OrderingU& ordering_u,
418 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
419 c.setBuildMode(Container::random);
421 for (std::size_t i = 0; i < c.N(); ++i)
422 c.setrowsize(i,p[i].size());
425 for (std::size_t i = 0; i < c.N(); ++i)
426 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
427 c.addindex(i,cit->first);
430 for (std::size_t i = 0; i < c.N(); ++i)
431 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
433 allocate_matrix(ordering_v.childOrdering(i),
434 ordering_u.childOrdering(cit->first),
440 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
441 typename std::enable_if<
445 allocate_matrix(
const OrderingV& ordering_v,
446 const OrderingU& ordering_u,
450 for (std::size_t i = 0; i < c.N(); ++i)
451 for (
typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
453 allocate_matrix(ordering_v.childOrdering(i),
454 ordering_u.childOrdering(cit->first),
460 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
461 typename std::enable_if<
464 allocate_matrix(
const OrderingV& ordering_v,
465 const OrderingU& ordering_u,
469 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
470 c.setBuildMode(Container::random);
472 for (std::size_t i = 0; i < c.N(); ++i)
473 c.setrowsize(i,p[i].size());
476 for (std::size_t i = 0; i < c.N(); ++i)
477 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
489 #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:139
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
const P & p
Definition: constraints.hh:147