dune-fem  2.8-git
operator/linear/hierarchical.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_OPERATOR_LINEAR_HIERARCHICAL_HH
2 #define DUNE_FEM_OPERATOR_LINEAR_HIERARCHICAL_HH
3 
4 #include <cstddef>
5 
6 #include <algorithm>
7 #include <string>
8 #include <type_traits>
9 #include <utility>
10 
11 #include <dune/common/exceptions.hh>
12 #include <dune/common/fmatrix.hh>
13 
14 #if HAVE_DUNE_ISTL
15 #include <dune/istl/bcrsmatrix.hh>
16 #include <dune/istl/bvector.hh>
17 #include <dune/istl/multitypeblockmatrix.hh>
18 #include <dune/istl/multitypeblockvector.hh>
19 #endif // #if HAVE_DUNE_ISTL
20 
24 
25 namespace Dune
26 {
27 
28  namespace Fem
29  {
30 
31  // Internal Forward Declaration
32  // ----------------------------
33 
34  template< class DomainFunction, class RangeFunction >
35  class HierarchicalLinearOperator;
36 
37 
38 
39  namespace Impl
40  {
41 
42  template< class Dof, class DomainIndices, class RangeIndices >
43  struct HierarchicalMatrixChooser;
44 
45 #if HAVE_DUNE_ISTL
46  template< class Dof, int rows, int cols >
47  struct HierarchicalMatrixChooser< Dof, Hybrid::IndexRange< int, cols >, Hybrid::IndexRange< int, rows > >
48  {
49  typedef BCRSMatrix< FieldMatrix< Dof, rows, cols > > Type;
50  };
51 
52  template< class Dof, class... SD, class... SR >
53  struct HierarchicalMatrixChooser< Dof, Hybrid::CompositeIndexRange< SD... >, Hybrid::CompositeIndexRange< SR... > >
54  {
55  private:
56  template< class R >
57  using Row = MultiTypeBlockVector< typename HierarchicalMatrixChooser< Dof, SD, R >::Type... >;
58 
59  public:
60  typedef MultiTypeBlockMatrix< Row< SR >... > Type;
61  };
62 #endif // #if HAVE_DUNE_ISTL
63 
64  } // namespace Impl
65 
66 
67 
68  // HierarchicalLinearOperator
69  // --------------------------
70 
71  template< class DomainFunction, class RangeFunction >
73  : public AssembledOperator< DomainFunction, RangeFunction >
74  {
77 
78  public:
79  typedef std::common_type_t< typename DomainFunction::DofType, typename RangeFunction::DofType > DofType;
80 
83 
84  typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
85  typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
86 
87  typedef typename DomainSpaceType::EntityType DomainEntityType;
88  typedef typename RangeSpaceType::EntityType RangeEntityType;
89 
90  // typedef DomainEntityType ColumnEntityType;
91  // typedef RangeEntityType RowEntityType;
92 
93  typedef typename Impl::HierarchicalMatrixChooser< DofType, typename DomainSpaceType::LocalBlockIndices, typename RangeSpaceType::LocalBlockIndices >::Type MatrixType;
94 
95  private:
96  template< class Functor >
97  static auto blockFunctor ( Functor &&functor )
98  {
99  return [ functor ] ( std::pair< std::size_t, std::size_t > local, const auto &global ) {
100  local.first *= Hybrid::size( typename RangeSpaceType::LocalBlockIndices() );
101  local.second *= Hybrid::size( typename DomainSpaceType::LocalBlockIndices() );
102  Hybrid::forEach( typename RangeSpaceType::LocalBlockIndices(), [ functor, local, global ] ( auto &&i ) {
103  Hybrid::forEach( typename DomainSpaceType::LocalBlockIndices(), [ functor, local, global, i ] ( auto &&j ) {
104  const auto iGlobal = std::make_pair( static_cast< std::size_t >( global.first ), i );
105  const auto jGlobal = std::make_pair( static_cast< std::size_t >( global.second ), j );
106  functor( std::make_pair( local.first + i, local.second + j ), std::make_pair( iGlobal, jGlobal ) );
107  } );
108  } );
109  };
110  }
111 
112  protected:
113  MatrixType &matrix () { return matrix_; }
114 
115  public:
117  : domainSpace_( domainSpace ), rangeSpace_( rangeSpace )
118  {}
119 
120  virtual void operator() ( const DomainFunction &u, RangeFunction &w ) const
121  {
122  w.clear();
123  umv( matrix_, u.dofVector().array(), w.dofVector().array() );
124  w.communicate();
125  }
126 
127  void communicate () {}
128 
129  const DomainSpaceType &domainSpace () const { return domainSpace_; }
130  const RangeSpaceType &rangeSpace () const { return domainSpace_; }
131 
132  MatrixType &exportMatrix () const { return matrix_; }
133 
134  template< class LocalMatrix >
135  void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix )
136  {
137  auto f = blockFunctor( [ this, &localMatrix ] ( auto local, auto global ) {
138  ThisType::entry( matrix_, global.first, global.second ) += localMatrix[ local.first ][ local.second ];
139  } );
140  rangeSpace().blockMapper().mapEach( rangeEntity, makePairFunctor( domainSpace().blockMapper(), domainEntity, std::move( f ) ) );
141  }
142 
143  template< class LocalMatrix, class Scalar >
144  void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix, const Scalar &scalar )
145  {
146  auto f = blockFunctor( [ this, &localMatrix, &scalar ] ( auto local, auto global ) {
147  ThisType::entry( matrix_, global.first, global.second ) += scalar * localMatrix[ local.first ][ local.second ];
148  } );
149  rangeSpace().blockMapper().mapEach( rangeEntity, makePairFunctor( domainSpace().blockMapper(), domainEntity, std::move( f ) ) );
150  }
151 
152  template< class LocalMatrix >
153  void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMatrix ) const
154  {
155  auto f = blockFunctor( [ this, &localMatrix ] ( auto local, auto global ) {
156  localMatrix[ local.first ][ local.second ] = ThisType::entry( matrix_, global.first, global.second );
157  } );
158  rangeSpace().blockMapper().mapEach( rangeEntity, makePairFunctor( domainSpace().blockMapper(), domainEntity, std::move( f ) ) );
159  }
160 
161  template< class LocalMatrix >
162  void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix )
163  {
164  auto f = blockFunctor( [ this, &localMatrix ] ( auto local, auto global ) {
165  ThisType::entry( matrix_, global.first, global.second ) = localMatrix[ local.first ][ local.second ];
166  } );
167  rangeSpace().blockMapper().mapEach( rangeEntity, makePairFunctor( domainSpace().blockMapper(), domainEntity, std::move( f ) ) );
168  }
169 
170  void clear () { clear( matrix_ ); }
171  template <class I>
172  void unitRow( const I localRow, const double diag = 1.0 )
173  {
174  DUNE_THROW(NotImplemented,"unitRow not implemented on HierarchicalLinearOperator");
175  }
176 
177  template< class Stencil >
178  void reserve ( const Stencil &stencil )
179  {
180  reserve( matrix_, stencil );
181  }
182 
183  private:
184 #if HAVE_DUNE_ISTL
185  template< class... C, class... B, class RangeVector >
186  static std::enable_if_t< sizeof...( C ) == sizeof...( B ) > umv ( const MultiTypeBlockVector< C... > &row, const MultiTypeBlockVector< B... > &u, RangeVector &w )
187  {
188  Hybrid::forEach( std::index_sequence_for< C... >(), [ &row, &u, &w ] ( auto &&j ) {
189  ThisType::umv( row[ j ], u[ j ], w );
190  } );
191  }
192 
193  template< class... R, class DomainVector, class... B >
194  static std::enable_if_t< sizeof...( R ) == sizeof...( B ) > umv ( const MultiTypeBlockMatrix< R... > &matrix, const DomainVector &u, MultiTypeBlockVector< B... > &w )
195  {
196  Hybrid::forEach( std::index_sequence_for< R... >(), [ &matrix, &u, &w ] ( auto &&i ) {
197  ThisType::umv( matrix[ i ], u, w[ i ] );
198  } );
199  }
200 
201  template< class K, int m, int n, class AM, class AU, class AW >
202  static void umv ( const BCRSMatrix< FieldMatrix< K, m, n >, AM > &matrix, const BlockVector< FieldVector< K, n >, AU > &u, BlockVector< FieldVector< K, m >, AW > &w )
203  {
204  matrix.umv( u, w );
205  }
206 
207  template< class... C >
208  static void clear ( MultiTypeBlockVector< C... > &row )
209  {
210  Hybrid::forEach( std::index_sequence_for< C... >(), [ &row ] ( auto &&j ) {
211  ThisType::clear( row[ j ] );
212  } );
213  }
214 
215  template< class... R >
216  static void clear ( MultiTypeBlockMatrix< R... > &matrix )
217  {
218  Hybrid::forEach( std::index_sequence_for< R... >(), [ &matrix ] ( auto &&i ) {
219  ThisType::clear( matrix[ i ] );
220  } );
221  }
222 
223  template< class K, int m, int n, class A >
224  static void clear ( BCRSMatrix< FieldMatrix< K, m, n >, A > &matrix )
225  {
226  for( auto &row : matrix )
227  std::fill( row.begin(), row.end(), FieldMatrix< K, m, n >( K( 0 ) ) );
228  }
229 
230  template< class... C, class I, std::size_t component, class J, J offset, class SJ >
231  static decltype( auto ) entry ( const MultiTypeBlockVector< C... > &row, I i, std::pair< std::size_t, Hybrid::CompositeIndex< component, J, offset, SJ > > j )
232  {
233  return entry( row[ std::integral_constant< std::size_t, component >() ], i, std::make_pair( j.first, j.second.subIndex() ) );
234  }
235 
236  template< class... C, class I, std::size_t component, class J, J offset, class SJ >
237  static decltype( auto ) entry ( MultiTypeBlockVector< C... > &row, I i, std::pair< std::size_t, Hybrid::CompositeIndex< component, J, offset, SJ > > j )
238  {
239  return entry( row[ std::integral_constant< std::size_t, component >() ], i, std::make_pair( j.first, j.second.subIndex() ) );
240  }
241 
242  template< class... R, std::size_t component, class I, I offset, class SI, class J >
243  static decltype( auto ) entry ( const MultiTypeBlockMatrix< R... > &matrix, std::pair< std::size_t, Hybrid::CompositeIndex< component, I, offset, SI > > i, J j )
244  {
245  return entry( matrix[ std::integral_constant< std::size_t, component >() ], std::make_pair( i.first, i.second.subIndex() ), j );
246  }
247 
248  template< class... R, std::size_t component, class I, I offset, class SI, class J >
249  static decltype( auto ) entry ( MultiTypeBlockMatrix< R... > &matrix, std::pair< std::size_t, Hybrid::CompositeIndex< component, I, offset, SI > > i, J j )
250  {
251  return entry( matrix[ std::integral_constant< std::size_t, component >() ], std::make_pair( i.first, i.second.subIndex() ), j );
252  }
253 
254  template< class K, int m, int n, class A >
255  static const K &entry ( const BCRSMatrix< FieldMatrix< K, m, n >, A > &matrix, std::pair< std::size_t, int > i, std::pair< std::size_t, int > j )
256  {
257  return matrix[ i.first ][ j.first ][ i.second ][ j.second ];
258  }
259 
260  template< class K, int m, int n, class A >
261  static K &entry ( BCRSMatrix< FieldMatrix< K, m, n >, A > &matrix, std::pair< std::size_t, int > i, std::pair< std::size_t, int > j )
262  {
263  return matrix[ i.first ][ j.first ][ i.second ][ j.second ];
264  }
265 
266  template< class... C, class Stencil >
267  static void reserve ( MultiTypeBlockVector< C... > &row, const Stencil &stencil )
268  {
269  Hybrid::forEach( std::index_sequence_for< C... >(), [ &row, &stencil ] ( auto &&j ) {
270  ThisType::reserve( row[ j ], stencil );
271  } );
272  }
273 
274  template< class... R, class Stencil >
275  static void reserve ( MultiTypeBlockMatrix< R... > &matrix, const Stencil &stencil )
276  {
277  Hybrid::forEach( std::index_sequence_for< R... >(), [ &matrix, &stencil ] ( auto &&i ) {
278  ThisType::reserve( matrix[ i ], stencil );
279  } );
280  }
281 
282  template< class K, int m, int n, class A, class Stencil >
283  static void reserve ( BCRSMatrix< FieldMatrix< K, m, n >, A > &matrix, const Stencil &stencil )
284  {
285  // reallocate matrix of correct size (destroys previously allocated matrix)
286  if( matrix.buildMode() == matrix.unknown )
287  matrix.setBuildMode( matrix.random );
288  matrix.setSize( stencil.rows(), stencil.cols() );
289 
290  // setup sparsity pattern (using random build mode)
291  for( const auto &row : stencil.globalStencil() )
292  matrix.setrowsize( row.first, row.second.size() );
293  matrix.endrowsizes();
294 
295  for( const auto &row : stencil.globalStencil() )
296  matrix.setIndices( row.first, row.second.begin(), row.second.end() );
297  matrix.endindices();
298  }
299 #endif // #if HAVE_DUNE_ISTL
300 
301  const DomainSpaceType &domainSpace_;
302  const RangeSpaceType &rangeSpace_;
303  mutable MatrixType matrix_;
304  };
305 
306  } // namespace Fem
307 
308 } // namespace Dune
309 
310 #endif // #ifndef DUNE_FEM_OPERATOR_LINEAR_HIERARCHICAL_HH
Definition: bindguard.hh:11
PairFunctor< Mapper, Entity, Functor > makePairFunctor(const Mapper &mapper, const Entity &entity, Functor functor)
Definition: operator/matrix/functor.hh:65
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
DomainFunction RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:38
abstract matrix operator
Definition: operator.hh:124
default implementation for a general operator stencil
Definition: stencil.hh:35
Definition: operator/linear/hierarchical.hh:74
void clear()
Definition: operator/linear/hierarchical.hh:170
MatrixType & matrix()
Definition: operator/linear/hierarchical.hh:113
Impl::HierarchicalMatrixChooser< DofType, typename DomainSpaceType::LocalBlockIndices, typename RangeSpaceType::LocalBlockIndices >::Type MatrixType
Definition: operator/linear/hierarchical.hh:93
MatrixType & exportMatrix() const
Definition: operator/linear/hierarchical.hh:132
void reserve(const Stencil &stencil)
Definition: operator/linear/hierarchical.hh:178
RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType
Definition: operator/linear/hierarchical.hh:85
void addLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix)
Definition: operator/linear/hierarchical.hh:135
const DomainSpaceType & domainSpace() const
Definition: operator/linear/hierarchical.hh:129
const RangeSpaceType & rangeSpace() const
Definition: operator/linear/hierarchical.hh:130
void unitRow(const I localRow, const double diag=1.0)
Definition: operator/linear/hierarchical.hh:172
RangeSpaceType::EntityType RangeEntityType
Definition: operator/linear/hierarchical.hh:88
void setLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix)
Definition: operator/linear/hierarchical.hh:162
DomainSpaceType::EntityType DomainEntityType
Definition: operator/linear/hierarchical.hh:87
void getLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMatrix) const
Definition: operator/linear/hierarchical.hh:153
BaseType::RangeFunctionType RangeFunctionType
Definition: operator/linear/hierarchical.hh:82
DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType
Definition: operator/linear/hierarchical.hh:84
HierarchicalLinearOperator(const std::string &, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace)
Definition: operator/linear/hierarchical.hh:116
void addScaledLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMatrix, const Scalar &scalar)
Definition: operator/linear/hierarchical.hh:144
std::common_type_t< typename DomainFunction::DofType, typename RangeFunction::DofType > DofType
Definition: operator/linear/hierarchical.hh:79
void communicate()
Definition: operator/linear/hierarchical.hh:127
virtual void operator()(const DomainFunction &u, RangeFunction &w) const
Definition: operator/linear/hierarchical.hh:120
BaseType::DomainFunctionType DomainFunctionType
Definition: operator/linear/hierarchical.hh:81