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 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)
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)
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)
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)
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)
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);
338 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
343 allocate_matrix(
const OrderingV& ordering_v,
344 const OrderingU& ordering_u,
348 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
349 c.setBuildMode(Container::random);
351 for (std::size_t i = 0; i < c.N(); ++i)
352 c.setrowsize(i,p[i].size());
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);
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)
363 allocate_matrix(ordering_v.childOrdering(i),
364 ordering_u.childOrdering(cit->first),
370 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
375 allocate_matrix(
const OrderingV& ordering_v,
376 const OrderingU& ordering_u,
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)
383 allocate_matrix(ordering_v.childOrdering(i),
384 ordering_u.childOrdering(cit->first),
390 template<
typename OrderingV,
typename OrderingU,
typename Pattern,
typename Container>
394 allocate_matrix(
const OrderingV& ordering_v,
395 const OrderingU& ordering_u,
399 c.setSize(ordering_v.blockCount(),ordering_u.blockCount(),
false);
400 c.setBuildMode(Container::random);
402 for (std::size_t i = 0; i < c.N(); ++i)
403 c.setrowsize(i,p[i].size());
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)
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