dune-fem  2.8-git
diagonalpreconditioner.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_DIAGONALPRECONDITIONER_HH
2 #define DUNE_FEM_DIAGONALPRECONDITIONER_HH
3 
4 #include <type_traits>
5 
8 
9 #if HAVE_DUNE_ISTL
10 #include <dune/istl/bvector.hh>
11 #endif
12 
13 namespace Dune
14 {
15 
16  namespace Fem
17  {
18 
19  // DiagonalPreconditionerBase (default non-assembled)
20  // --------------------------------------------------
21  template< class DFImp, class OperatorImp, bool assembled >
23  : public Operator< DFImp, DFImp >
24  {
25  public:
26  typedef DFImp DiscreteFunctionType;
27  typedef OperatorImp OperatorType;
28 
29  typedef typename DiscreteFunctionType :: DofIteratorType DofIteratorType;
30  typedef typename DiscreteFunctionType :: ConstDofIteratorType ConstDofIteratorType;
31 
32  public:
34 
35  virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
36  {
37  DUNE_THROW(NotImplemented,"preconditioning not possible for non-assembled operators");
38  }
39 
40 #if HAVE_DUNE_ISTL
42  template < class YBlock, class XBlock >
43  void applyToISTLBlockVector( const BlockVector< YBlock >& d,
44  BlockVector< XBlock >& v ) const
45  {
46  DUNE_THROW(NotImplemented,"preconditioning not possible for non-assembled operators");
47  }
48 #endif
49 
50  protected:
51  void apply( const DiscreteFunctionType& u, DiscreteFunctionType& res ) const
52  {
53  DUNE_THROW(NotImplemented,"preconditioning not possible for non-assembled operators");
54  }
55  };
56 
57 
58  // DiagonalPreconditionerBase (assembled version)
59  // ----------------------------------------------
60  template< class DFImp, class OperatorImp >
61  class DiagonalPreconditionerBase< DFImp, OperatorImp, true >
62  : public Operator< DFImp, DFImp >
63  {
64  public:
65  typedef DFImp DiscreteFunctionType;
66  typedef OperatorImp OperatorType;
67 
68  typedef typename DiscreteFunctionType :: DofType DofType;
69  typedef typename Dune::FieldTraits< DofType >::real_type RealType;
70  typedef typename DiscreteFunctionType :: DofIteratorType DofIteratorType;
71  typedef typename DiscreteFunctionType :: ConstDofIteratorType ConstDofIteratorType;
72 
73  protected:
75 
76  public:
77  DiagonalPreconditionerBase( const OperatorType& assembledOperator )
78  : diagonalInv_( "diag-preconditioning", assembledOperator.rangeSpace() )
79  {
80  // estract diagonal elements form matrix object
81  assembledOperator.extractDiagonal( diagonalInv_ );
82 
83  // make consistent at border dofs
84  diagonalInv_.communicate();
85 
86  // In general: store 1/diag
87  //
88  // note: We set near-zero entries to 1 to avoid NaNs. Such entries occur
89  // if DoFs are excluded from matrix setup
90  const RealType eps = 16.*std::numeric_limits< RealType >::epsilon();
91  const DofIteratorType dend = diagonalInv_.dend();
92  for( DofIteratorType dit = diagonalInv_.dbegin(); dit != dend; ++dit )
93  *dit = (std::abs( *dit ) < eps ? DofType( 1. ) : DofType( 1. ) / *dit);
94  }
95 
96  virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
97  {
98  apply(u, res);
99  }
100 
101 #if HAVE_DUNE_ISTL
103  template < class YBlock, class XBlock >
104  void applyToISTLBlockVector( const BlockVector< YBlock >& d,
105  BlockVector< XBlock >& v ) const
106  {
107  DiscreteFunctionType vTmp("diag-precon::X", diagonalInv_.space(), v );
108  DiscreteFunctionType dTmp("diag-precon::Y", diagonalInv_.space(), d );
109 
110  // apply 1/diagonal
111  apply( dTmp, vTmp );
112  }
113 #endif
114 
115 
116  protected:
117  void apply( const DiscreteFunctionType& u, DiscreteFunctionType& res ) const
118  {
119  ConstDofIteratorType uIt = u.dbegin();
120  ConstDofIteratorType diagInv = diagonalInv_.dbegin();
121 
122  const DofIteratorType resEnd = res.dend();
123 
124  // apply 1/diagonal
125  for(DofIteratorType resIt = res.dbegin();
126  resIt != resEnd; ++ resIt, ++diagInv, ++ uIt )
127  {
128  assert( diagInv != diagonalInv_.dend() );
129  assert( uIt != u.dend() );
130  (*resIt) = (*uIt) * (*diagInv);
131  }
132  }
133 
134  };
135 
136 
137  // DiagonalPreconditioner
138  // ----------------------
149  template< class DFImp, class Operator>
151  : public DiagonalPreconditionerBase< DFImp, Operator, std::is_base_of< AssembledOperator< DFImp, DFImp >, Operator > :: value >
152  {
154  BaseType;
155  public:
158  : BaseType( op )
159  {}
160  };
161 
162  } // namespace Fem
163 
164 } // namespace Dune
165 
166 #endif // #ifndef DUNE_FEM_DIAGONALPRECONDITIONER_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
abstract operator
Definition: operator.hh:34
Definition: diagonalpreconditioner.hh:24
DFImp DiscreteFunctionType
Definition: diagonalpreconditioner.hh:26
void apply(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
Definition: diagonalpreconditioner.hh:51
DiscreteFunctionType ::ConstDofIteratorType ConstDofIteratorType
Definition: diagonalpreconditioner.hh:30
DiagonalPreconditionerBase(const OperatorType &op)
Definition: diagonalpreconditioner.hh:33
virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
application operator
Definition: diagonalpreconditioner.hh:35
OperatorImp OperatorType
Definition: diagonalpreconditioner.hh:27
DiscreteFunctionType ::DofIteratorType DofIteratorType
Definition: diagonalpreconditioner.hh:29
DiscreteFunctionType ::DofIteratorType DofIteratorType
Definition: diagonalpreconditioner.hh:70
Dune::FieldTraits< DofType >::real_type RealType
Definition: diagonalpreconditioner.hh:69
DiagonalPreconditionerBase(const OperatorType &assembledOperator)
Definition: diagonalpreconditioner.hh:77
DFImp DiscreteFunctionType
Definition: diagonalpreconditioner.hh:65
OperatorImp OperatorType
Definition: diagonalpreconditioner.hh:66
DiscreteFunctionType ::DofType DofType
Definition: diagonalpreconditioner.hh:68
DiscreteFunctionType ::ConstDofIteratorType ConstDofIteratorType
Definition: diagonalpreconditioner.hh:71
virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
application operator
Definition: diagonalpreconditioner.hh:96
DiscreteFunctionType diagonalInv_
Definition: diagonalpreconditioner.hh:74
void apply(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
Definition: diagonalpreconditioner.hh:117
Precondtioner, multiplies with inverse of the diagonal works with.
Definition: diagonalpreconditioner.hh:152
Operator OperatorType
Definition: diagonalpreconditioner.hh:156
DiagonalPreconditioner(const OperatorType &op)
Definition: diagonalpreconditioner.hh:157