1#ifndef DUNE_FEM_OPERATOR_LINEAR_HIERARCHICAL_HH
2#define DUNE_FEM_OPERATOR_LINEAR_HIERARCHICAL_HH
11#include <dune/common/exceptions.hh>
12#include <dune/common/fmatrix.hh>
15#include <dune/istl/bcrsmatrix.hh>
16#include <dune/istl/bvector.hh>
17#include <dune/istl/multitypeblockmatrix.hh>
18#include <dune/istl/multitypeblockvector.hh>
34 template<
class DomainFunction,
class RangeFunction >
35 class HierarchicalLinearOperator;
42 template<
class Dof,
class DomainIndices,
class RangeIndices >
43 struct HierarchicalMatrixChooser;
46 template<
class Dof,
int rows,
int cols >
47 struct HierarchicalMatrixChooser< Dof, Hybrid::IndexRange< int, cols >, Hybrid::IndexRange< int, rows > >
49 typedef BCRSMatrix< FieldMatrix< Dof, rows, cols > > Type;
52 template<
class Dof,
class... SD,
class... SR >
53 struct HierarchicalMatrixChooser< Dof, Hybrid::CompositeIndexRange< SD... >, Hybrid::CompositeIndexRange< SR... > >
57 using Row = MultiTypeBlockVector< typename HierarchicalMatrixChooser< Dof, SD, R >::Type... >;
60 typedef MultiTypeBlockMatrix< Row< SR >... > Type;
71 template<
class DomainFunction,
class RangeFunction >
79 typedef std::common_type_t< typename DomainFunction::DofType, typename RangeFunction::DofType >
DofType;
85 typedef typename RangeFunctionType::DiscreteFunctionSpaceType
RangeSpaceType;
93 typedef typename Impl::HierarchicalMatrixChooser< DofType, typename DomainSpaceType::LocalBlockIndices, typename RangeSpaceType::LocalBlockIndices >::Type
MatrixType;
96 template<
class Functor >
97 static auto blockFunctor ( Functor &&functor )
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 ) );
120 virtual void operator() (
const DomainFunction &u, RangeFunction &w )
const
123 umv( matrix_, u.dofVector().array(), w.dofVector().array() );
134 template<
class LocalMatrix >
137 auto f = blockFunctor( [
this, &localMatrix ] (
auto local,
auto global ) {
138 ThisType::entry( matrix_, global.first, global.second ) += localMatrix[ local.first ][ local.second ];
143 template<
class LocalMatrix,
class Scalar >
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 ];
152 template<
class LocalMatrix >
155 auto f = blockFunctor( [
this, &localMatrix ] (
auto local,
auto global ) {
156 localMatrix[ local.first ][ local.second ] = ThisType::entry( matrix_, global.first, global.second );
161 template<
class LocalMatrix >
164 auto f = blockFunctor( [
this, &localMatrix ] (
auto local,
auto global ) {
165 ThisType::entry( matrix_, global.first, global.second ) = localMatrix[ local.first ][ local.second ];
172 void unitRow(
const I localRow,
const double diag = 1.0 )
174 DUNE_THROW(NotImplemented,
"unitRow not implemented on HierarchicalLinearOperator");
177 template<
class Stencil >
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 )
188 Hybrid::forEach( std::index_sequence_for< C... >(), [ &row, &u, &w ] (
auto &&j ) {
189 ThisType::umv( row[ j ], u[ j ], w );
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 )
197 ThisType::umv(
matrix[ i ], u, w[ i ] );
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 )
207 template<
class... C >
208 static void clear ( MultiTypeBlockVector< C... > &row )
210 Hybrid::forEach( std::index_sequence_for< C... >(), [ &row ] (
auto &&j ) {
215 template<
class... R >
216 static void clear ( MultiTypeBlockMatrix< R... > &
matrix )
223 template<
class K,
int m,
int n,
class A >
224 static void clear ( BCRSMatrix< FieldMatrix< K, m, n >, A > &
matrix )
227 std::fill( row.begin(), row.end(), FieldMatrix< K, m, n >( K( 0 ) ) );
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 )
233 return entry( row[ std::integral_constant< std::size_t, component >() ], i, std::make_pair( j.first, j.second.subIndex() ) );
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 )
239 return entry( row[ std::integral_constant< std::size_t, component >() ], i, std::make_pair( j.first, j.second.subIndex() ) );
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 )
245 return entry(
matrix[ std::integral_constant< std::size_t, component >() ], std::make_pair( i.first, i.second.subIndex() ), j );
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 )
251 return entry(
matrix[ std::integral_constant< std::size_t, component >() ], std::make_pair( i.first, i.second.subIndex() ), j );
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 )
257 return matrix[ i.first ][ j.first ][ i.second ][ j.second ];
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 )
263 return matrix[ i.first ][ j.first ][ i.second ][ j.second ];
266 template<
class... C,
class Stencil >
267 static void reserve ( MultiTypeBlockVector< C... > &row,
const Stencil &stencil )
269 Hybrid::forEach( std::index_sequence_for< C... >(), [ &row, &stencil ] (
auto &&j ) {
274 template<
class... R,
class Stencil >
275 static void reserve ( MultiTypeBlockMatrix< R... > &
matrix,
const Stencil &stencil )
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 )
288 matrix.setSize( stencil.rows(), stencil.cols() );
291 const auto& globalStencil = stencil.globalStencil();
292 const std::size_t nRows = globalStencil.size();
293 for( std::size_t row = 0; row<nRows; ++row )
296 const auto& columns = globalStencil.at( row );
297 matrix.setrowsize( row, columns.size() );
299 catch (
const std::out_of_range& e )
307 for( std::size_t row = 0; row<nRows; ++row )
310 const auto& columns = globalStencil.at( row );
311 matrix.setIndices( row, columns.begin(), columns.end() );
313 catch (
const std::out_of_range& e )
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
RangeFunction 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
const RangeSpaceType & rangeSpace() const
Definition: operator/linear/hierarchical.hh:130
void clear()
Definition: operator/linear/hierarchical.hh:170
void reserve(const Stencil &stencil)
Definition: operator/linear/hierarchical.hh:178
MatrixType & exportMatrix() const
Definition: operator/linear/hierarchical.hh:132
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
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
Impl::HierarchicalMatrixChooser< DofType, typenameDomainSpaceType::LocalBlockIndices, typenameRangeSpaceType::LocalBlockIndices >::Type MatrixType
Definition: operator/linear/hierarchical.hh:93
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
MatrixType & matrix()
Definition: operator/linear/hierarchical.hh:113
const DomainSpaceType & domainSpace() const
Definition: operator/linear/hierarchical.hh:129
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
application operator
Definition: operator/linear/hierarchical.hh:120
BaseType::DomainFunctionType DomainFunctionType
Definition: operator/linear/hierarchical.hh:81