1 #ifndef MASSOPERATOR_HH
2 #define MASSOPERATOR_HH
3
4 #include <dune/common/dynvector.hh>
5
6 #include <dune/grid/common/rangegenerators.hh>
7
8 #include <dune/fem/common/bindguard.hh>
9 #include <dune/fem/function/common/localcontribution.hh>
10 #include <dune/fem/function/localfunction/const.hh>
11 #include <dune/fem/operator/common/stencil.hh>
12 #include <dune/fem/operator/common/operator.hh>
13 #include <dune/fem/operator/common/localcontribution.hh>
14 #include <dune/fem/operator/common/differentiableoperator.hh>
15
16 #include <dune/fem/quadrature/cachingquadrature.hh>
17
18
19 template< class DiscreteFunction, class LinearOperator >
20 class MassOperator
21 : public LinearOperator
22 {
23 typedef MassOperator< DiscreteFunction, LinearOperator > ThisType;
24 typedef LinearOperator BaseType;
25
26 typedef DiscreteFunction DiscreteFunctionType;
27 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
28
29 typedef typename DiscreteFunctionSpaceType::DomainType DomainType;
30 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
31
32 typedef Dune::Fem::CachingQuadrature< typename DiscreteFunctionSpaceType::GridPartType, 0 > QuadratureType;
33 public:
MassOperator(const DiscreteFunctionSpaceType & dfSpace)34 explicit MassOperator ( const DiscreteFunctionSpaceType &dfSpace )
35 : BaseType( "mass operator", dfSpace , dfSpace ),
36 dfSpace_( dfSpace )
37 {
38 assemble(*this);
39 }
40
41 template< class Function >
42 void assembleRHS( const Function &u, DiscreteFunctionType &w ) const;
43 void assemble (LinearOperator &matrix) const;
44 private:
45 const DiscreteFunctionSpaceType &dfSpace_;
46 };
47
48 template< class DiscreteFunction, class LinearOperator >
49 class AffineMassOperator
50 : public Dune::Fem::DifferentiableOperator<LinearOperator>
51 {
52 public:
53 typedef DiscreteFunction DiscreteFunctionType;
54 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
55 typedef LinearOperator JacobianOperatorType;
56 template< class Function >
AffineMassOperator(const DiscreteFunctionSpaceType & dfSpace,const Function & u)57 explicit AffineMassOperator ( const DiscreteFunctionSpaceType &dfSpace, const Function &u )
58 : matrix_(dfSpace), rhs_("rhs", dfSpace)
59 {
60 matrix_.assembleRHS(u, rhs_);
61 }
rhs() const62 const DiscreteFunctionType rhs() const { return rhs_; }
operator ()(const DiscreteFunctionType & u,DiscreteFunctionType & w) const63 virtual void operator() ( const DiscreteFunctionType &u, DiscreteFunctionType &w ) const
64 {
65 matrix_(u,w);
66 w -= rhs_;
67 }
jacobian(const DiscreteFunctionType & u,JacobianOperatorType & jOp) const68 virtual void jacobian ( const DiscreteFunctionType &u, JacobianOperatorType &jOp ) const
69 {
70 matrix_.assemble( jOp );
71 }
72 private:
73 MassOperator<DiscreteFunction,LinearOperator> matrix_;
74 DiscreteFunction rhs_;
75 };
76
77 template< class DiscreteFunction, class LinearOperator >
78 template< class Function >
79 void MassOperator< DiscreteFunction, LinearOperator >
assembleRHS(const Function & u,DiscreteFunctionType & w) const80 ::assembleRHS( const Function &u, DiscreteFunctionType &w ) const
81 {
82 // clear result
83 w.clear();
84
85 Dune::Fem::ConstLocalFunction< Function > uLocal( u );
86 Dune::Fem::AddLocalContribution< DiscreteFunctionType > wLocal( w );
87
88 // run over entities
89 for( const auto &entity : elements( dfSpace_.gridPart(), Dune::Partitions::interiorBorder ) )
90 {
91 const auto geometry = entity.geometry();
92
93 auto guard = bindGuard( std::tie( uLocal, wLocal ), entity );
94
95 // run over quadrature points
96 for( const auto qp : QuadratureType( entity, uLocal.order()+wLocal.order() ) )
97 {
98 // evaluate u
99 RangeType uValue = uLocal.evaluate( qp );
100
101 // apply quadrature weight
102 uValue *= qp.weight() * geometry.integrationElement( qp.position() );
103
104 // add to local function
105 wLocal.axpy( qp, uValue );
106 }
107 }
108 }
109
110
111 template< class DiscreteFunction, class LinearOperator >
assemble(LinearOperator & matrix) const112 void MassOperator< DiscreteFunction, LinearOperator >::assemble (LinearOperator &matrix) const
113 {
114 matrix.reserve( Dune::Fem::DiagonalStencil<DiscreteFunctionSpaceType,DiscreteFunctionSpaceType>( dfSpace_, dfSpace_ ) );
115 matrix.clear();
116
117 Dune::Fem::AddLocalContribution< LinearOperator > localMatrix( matrix );
118 Dune::DynamicVector< typename DiscreteFunctionSpaceType::RangeType > values;
119
120 // run over entities
121 for( const auto &entity : elements( dfSpace_.gridPart(), Dune::Partitions::interiorBorder ) )
122 {
123 const auto geometry = entity.geometry();
124
125 auto guard = bindGuard( localMatrix, entity, entity );
126
127 const auto &basis = localMatrix.domainBasisFunctionSet();
128 const unsigned int numBasisFunctions = basis.size();
129 values.resize( numBasisFunctions );
130
131 // run over quadrature points
132 for( const auto qp : QuadratureType( entity, localMatrix.order() ) )
133 {
134 // evaluate base functions
135 basis.evaluateAll( qp, values );
136
137 // apply quadrature weight
138 const auto weight = qp.weight() * geometry.integrationElement( qp.position() );
139 std::for_each( values.begin(), values.end(), [ weight ] ( auto &a ) { a *= weight; } );
140
141 // update system matrix
142 localMatrix.axpy( qp, values );
143 }
144 }
145 }
146
147 #endif // #ifndef MASSOPERATOR_HH
148