dune-fem 2.8.0
Loading...
Searching...
No Matches
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
13namespace 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:
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
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
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
DiagonalPreconditionerBase(const OperatorType &op)
Definition: diagonalpreconditioner.hh:33
virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
application operator
Definition: diagonalpreconditioner.hh:35
DiscreteFunctionType::ConstDofIteratorType ConstDofIteratorType
Definition: diagonalpreconditioner.hh:30
DiscreteFunctionType::DofIteratorType DofIteratorType
Definition: diagonalpreconditioner.hh:29
OperatorImp OperatorType
Definition: diagonalpreconditioner.hh:27
DiscreteFunctionType::ConstDofIteratorType ConstDofIteratorType
Definition: diagonalpreconditioner.hh:71
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::DofIteratorType DofIteratorType
Definition: diagonalpreconditioner.hh:70
virtual void operator()(const DiscreteFunctionType &u, DiscreteFunctionType &res) const
application operator
Definition: diagonalpreconditioner.hh:96
Dune::FieldTraits< DofType >::real_type RealType
Definition: diagonalpreconditioner.hh:69
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