3 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH 4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXHELPERS_HH 11 #include <unordered_map> 12 #include <unordered_set> 14 #include <dune/common/fmatrix.hh> 15 #include <dune/istl/bcrsmatrix.hh> 27 template<
typename RV,
typename CV,
typename block_type>
28 struct matrix_for_vectors;
30 template<
typename B1,
typename A1,
typename B2,
typename A2,
typename block_type>
31 struct matrix_for_vectors<
Dune::BlockVector<B1,A1>,Dune::BlockVector<B2,A2>,block_type>
33 typedef Dune::BCRSMatrix<block_type> type;
36 template<
typename B1,
int n1,
typename B2,
int n2,
typename block_type>
37 struct matrix_for_vectors<
Dune::FieldVector<B1,n1>,Dune::FieldVector<B2,n2>,block_type>
39 typedef Dune::FieldMatrix<block_type,n1,n2> type;
42 template<
typename E,
typename RV,
typename CV, std::
size_t blocklevel>
43 struct recursive_build_matrix_type
45 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;
48 template<
typename E,
typename RV,
typename CV>
49 struct recursive_build_matrix_type<E,RV,CV,1>
51 typedef Dune::FieldMatrix<E,RV::dimension,CV::dimension> type;
55 template<
typename E,
typename RV,
typename CV>
56 struct build_matrix_type
59 static_assert(static_cast<int>(RV::blocklevel) == static_cast<int>(CV::blocklevel),
"Both vectors must have identical blocking depth");
61 typedef typename recursive_build_matrix_type<E,RV,CV,RV::blocklevel>::type type;
65 template<
typename RowOrdering,
typename ColOrdering,
typename SubPattern_ =
void>
67 :
public std::vector<std::unordered_map<std::size_t,SubPattern_> >
72 typedef SubPattern_ SubPattern;
74 template<
typename RI,
typename CI>
75 void add_link(
const RI& ri,
const CI& ci)
77 recursive_add_entry(ri.view(),ci.view());
80 template<
typename RI,
typename CI>
81 void recursive_add_entry(
const RI& ri,
const CI& ci)
83 this->resize(_row_ordering.blockCount());
84 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()))));
85 r.first->second.recursive_add_entry(ri.back_popped(),ci.back_popped());
88 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
89 : _row_ordering(row_ordering)
90 , _col_ordering(col_ordering)
95 const RowOrdering& _row_ordering;
96 const ColOrdering& _col_ordering;
100 template<
typename RowOrdering,
typename ColOrdering>
101 class Pattern<RowOrdering,ColOrdering,void>
102 :
public std::vector<std::unordered_set<std::size_t> >
107 typedef void SubPattern;
109 template<
typename RI,
typename CI>
110 void add_link(
const RI& ri,
const CI& ci)
112 recursive_add_entry(ri,ci);
115 template<
typename RI,
typename CI>
116 void recursive_add_entry(
const RI& ri,
const CI& ci)
118 this->resize(_row_ordering.blockCount());
119 (*this)[ri.back()].insert(ci.back());
122 Pattern(
const RowOrdering& row_ordering,
const ColOrdering& col_ordering)
123 : _row_ordering(row_ordering)
124 , _col_ordering(col_ordering)
129 const RowOrdering& _row_ordering;
130 const ColOrdering& _col_ordering;
134 template<
typename M,
int blocklevel = M::blocklevel>
135 struct requires_pattern
137 static const bool value = requires_pattern<
typename M::block_type,blocklevel-1>
::value;
141 struct requires_pattern<M,0>
143 static const bool value =
false;
146 template<
typename B,
typename A,
int blocklevel>
147 struct requires_pattern<
Dune::BCRSMatrix<B,A>,blocklevel>
149 static const bool value =
true;
152 template<
typename M,
typename RowOrdering,
typename ColOrdering,
bool pattern>
153 struct _build_pattern_type
158 template<
typename M,
typename RowOrdering,
typename ColOrdering>
159 struct _build_pattern_type<M,RowOrdering,ColOrdering,true>
164 template<
typename M,
typename GFSV,
typename GFSU,
typename Tag>
165 struct build_pattern_type
168 typedef OrderingBase<
169 typename GFSV::Ordering::Traits::DOFIndex,
170 typename GFSV::Ordering::Traits::ContainerIndex
173 typedef OrderingBase<
174 typename GFSU::Ordering::Traits::DOFIndex,
175 typename GFSU::Ordering::Traits::ContainerIndex
181 template<
typename M,
typename GFSV,
typename GFSU>
182 struct build_pattern_type<M,GFSV,GFSU,FlatContainerAllocationTag>
184 typedef Pattern<typename GFSV::Ordering, typename GFSU::Ordering> type;
188 template<
typename RI,
typename CI,
typename Block>
189 typename Block::field_type&
190 access_matrix_element(tags::field_matrix_1_1, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
197 template<
typename RI,
typename CI,
typename Block>
198 typename Block::field_type&
199 access_matrix_element(tags::field_matrix_n_m, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
203 return b[ri[0]][ci[0]];
206 template<
typename RI,
typename CI,
typename Block>
207 typename Block::field_type&
208 access_matrix_element(tags::field_matrix_1_m, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
215 template<
typename RI,
typename CI,
typename Block>
216 typename Block::field_type&
217 access_matrix_element(tags::field_matrix_n_1, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
224 template<
typename RI,
typename CI,
typename Block>
225 typename Block::field_type&
226 access_matrix_element(tags::bcrs_matrix, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
228 return access_matrix_element(
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
232 template<
typename RI,
typename CI,
typename Block>
233 const typename Block::field_type&
234 access_matrix_element(tags::field_matrix_1_1,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
241 template<
typename RI,
typename CI,
typename Block>
242 const typename Block::field_type&
243 access_matrix_element(tags::field_matrix_n_m,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
247 return b[ri[0]][ci[0]];
250 template<
typename RI,
typename CI,
typename Block>
251 const typename Block::field_type&
252 access_matrix_element(tags::field_matrix_1_m,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
259 template<
typename RI,
typename CI,
typename Block>
260 const typename Block::field_type&
261 access_matrix_element(tags::field_matrix_n_1,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
268 template<
typename RI,
typename CI,
typename Block>
269 const typename Block::field_type&
270 access_matrix_element(tags::bcrs_matrix,
const Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
272 return access_matrix_element(
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
277 template<
typename RI,
typename Block>
278 void clear_matrix_row(tags::field_matrix_1_any, Block& b,
const RI& ri,
int i)
284 template<
typename RI,
typename Block>
285 void clear_matrix_row(tags::field_matrix_n_any, Block& b,
const RI& ri,
int i)
291 template<
typename RI,
typename Block>
292 void clear_matrix_row(tags::bcrs_matrix, Block& b,
const RI& ri,
int i)
294 typedef typename Block::ColIterator col_iterator_type;
295 const col_iterator_type end = b[ri[i]].end();
296 for(col_iterator_type cit = b[ri[i]].begin(); cit != end; ++cit)
301 template<
typename T,
typename RI,
typename CI,
typename Block>
302 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)
309 template<
typename T,
typename RI,
typename CI,
typename Block>
310 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)
317 template<
typename T,
typename RI,
typename CI,
typename Block>
318 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)
325 template<
typename T,
typename RI,
typename CI,
typename Block>
326 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)
333 template<
typename T,
typename RI,
typename CI,
typename Block>
334 void write_matrix_element_if_exists(
const T& v, tags::bcrs_matrix, Block& b,
const RI& ri,
const CI& ci,
int i,
int j)
336 if (b.exists(ri[i],ci[j]))
337 write_matrix_element_if_exists(v,
container_tag(b[ri[i]][ci[j]]),b[ri[i]][ci[j]],ri,ci,i-1,j-1);
341 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
342 typename std::enable_if<
346 allocate_matrix(
const OrderingV& ordering_v,
347 const OrderingU& ordering_u,
351 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
352 c.setBuildMode(Container::random);
354 for (std::size_t i = 0; i < c.N(); ++i)
355 c.setrowsize(i,p[i].size());
358 for (std::size_t i = 0; i < c.N(); ++i)
359 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
360 c.addindex(i,cit->first);
363 for (std::size_t i = 0; i < c.N(); ++i)
364 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
366 allocate_matrix(ordering_v.childOrdering(i),
367 ordering_u.childOrdering(cit->first),
373 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
374 typename std::enable_if<
378 allocate_matrix(
const OrderingV& ordering_v,
379 const OrderingU& ordering_u,
383 for (std::size_t i = 0; i < c.N(); ++i)
384 for (
typename Pattern::value_type::iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
386 allocate_matrix(ordering_v.childOrdering(i),
387 ordering_u.childOrdering(cit->first),
393 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
394 typename std::enable_if<
397 allocate_matrix(
const OrderingV& ordering_v,
398 const OrderingU& ordering_u,
402 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
403 c.setBuildMode(Container::random);
405 for (std::size_t i = 0; i < c.N(); ++i)
406 c.setrowsize(i,p[i].size());
409 for (std::size_t i = 0; i < c.N(); ++i)
410 for (
typename Pattern::value_type::const_iterator cit = p[i].begin(); cit != p[i].end(); ++cit)
422 #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:249
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const P & p
Definition: constraints.hh:147
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139