dune-fem 2.8-git
molgalerkin.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SCHEMES_MOLGALERKIN_HH
2#define DUNE_FEM_SCHEMES_MOLGALERKIN_HH
3
4// fem includes
6
7namespace Dune
8{
9
10 namespace Fem
11 {
12
13 // GalerkinOperator
14 // ----------------
15
16 template< class Integrands, class DomainFunction, class RangeFunction = DomainFunction >
18 : public virtual Operator< DomainFunction, RangeFunction >
19 {
20 typedef DomainFunction DomainFunctionType;
21 typedef RangeFunction RangeFunctionType;
22
23
24 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value, "DomainFunction and RangeFunction must be defined on the same grid part." );
25
26 typedef typename RangeFunctionType::GridPartType GridPartType;
27
28 typedef Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType;
29 typedef typename RangeFunctionType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
30
31 typedef typename GalerkinOperatorImplType::template QuadratureSelector<
33
35
36 template< class... Args >
37 explicit MOLGalerkinOperator ( const GridPartType &gridPart, Args &&... args )
38 : impl_( gridPart, std::forward< Args >( args )... ),
39 communicate_( true )
40 {
41 // disable communicate in Impl::GalerkinOperator
42 // since applyInverseMass has to be applied first
43 impl_.setCommunicate( false );
44 }
45
46 void setCommunicate( const bool communicate ) { communicate_ = communicate; }
47 void setQuadratureOrders(unsigned int interior, unsigned int surface) { impl_.setQuadratureOrders(interior,surface); }
48
49 virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const final override
50 {
51 evaluate( u, w );
52 }
53
54 template< class GridFunction >
55 void operator() ( const GridFunction &u, RangeFunctionType &w ) const
56 {
57 evaluate( u, w );
58 }
59
60 const GridPartType &gridPart () const { return impl_.gridPart(); }
61
62 typedef Integrands ModelType;
63 typedef Integrands DirichletModelType;
64 ModelType &model() const { return impl_.model(); }
65
66 protected:
68 {
69 // get-set local contribution
71
72 LocalMassMatrixType localMassMatrix( w.space(), impl_.interiorQuadratureOrder( w.space().order() ) );
73
74 // iterate over all elements
75 for( const auto& entity : elements( gridPart(), Partitions::interiorBorder ) )
76 {
77 // fill local contribution
78 auto guard = bindGuard( wLocal, entity );
79
80 // apply inverse mass matrix
81 // TODO: add mass term if needed (from ufl expression)
82 localMassMatrix.applyInverse( entity, wLocal );
83 }
84 }
85
86 template< class GridFunction >
87 void evaluate( const GridFunction &u, RangeFunctionType &w ) const
88 {
89 // Impl::GalerkinOperator::evaluate without communicate
90 impl_.evaluate( u, w );
91
92 // method of lines
94
95 // synchronize data
96 if( communicate_ )
97 w.communicate();
98 }
99
100 // GalerkinOperator implementation (see galerkin.hh)
103 };
104
105
106
107 // DifferentiableGalerkinOperator
108 // ------------------------------
109
110 template< class Integrands, class JacobianOperator >
112 : public MOLGalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
113 public DifferentiableOperator< JacobianOperator >
114 {
116
117 public:
118 typedef JacobianOperator JacobianOperatorType;
119
122 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
123 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
124
126
127 template< class... Args >
129 const RangeDiscreteFunctionSpaceType &rSpace,
130 Args &&... args )
131 : BaseType( rSpace.gridPart(), std::forward< Args >( args )... ),
132 dSpace_(dSpace),
133 rSpace_(rSpace)
134 {}
135
136 virtual void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const final override
137 {
138 // assemble Jacobian, same as GalerkinOperator
139 impl_.assemble( u, jOp );
140 // apply inverse mass
141 applyInverseMass( jOp, impl_.model().hasSkeleton() );
142 }
143
144 template< class GridFunction >
145 void jacobian ( const GridFunction &u, JacobianOperatorType &jOp ) const
146 {
147 // assemble Jacobian, same as GalerkinOperator
148 impl_.assemble( u, jOp );
149 // apply inverse mass
150 applyInverseMass( jOp, impl_.model().hasSkeleton() );
151 }
152
154 {
155 return dSpace_;
156 }
158 {
159 return rSpace_;
160 }
161
162 using BaseType::gridPart;
163
164 protected:
165 void applyInverseMass ( JacobianOperatorType &jOp, const bool hasSkeleton ) const
166 {
168
169 LocalMassMatrixType localMassMatrix( jOp.rangeSpace(), impl_.interiorQuadratureOrder( jOp.rangeSpace().order() ) );
170
172
173 // multiply with inverse mass matrix
174 for( const auto& inside : elements( gridPart(), Partitions::interiorBorder ) )
175 {
176 // scale diagonal
177 {
178 auto guard = bindGuard( jOpLocal, inside, inside );
179 localMassMatrix.leftMultiplyInverse( jOpLocal );
180 }
181
182 if( hasSkeleton )
183 {
184 for( const auto &intersection : intersections( gridPart(), inside ) )
185 {
186 // scale off-diagonal
187 if( intersection.neighbor() )
188 {
189 const auto& outside = intersection.outside();
190 auto guard = bindGuard( jOpLocal, outside, inside );
191 localMassMatrix.leftMultiplyInverse( jOpLocal );
192 }
193 }
194 }
195 }
196 }
197
198 using BaseType::impl_;
201 };
202
203
204
205 // AutomaticDifferenceGalerkinOperator
206 // -----------------------------------
207
208 template< class Integrands, class DomainFunction, class RangeFunction >
210 : public MOLGalerkinOperator< Integrands, DomainFunction, RangeFunction >,
211 public AutomaticDifferenceOperator< DomainFunction, RangeFunction >
212 {
215
216 public:
218
219 template< class... Args >
221 : BaseType( gridPart, std::forward< Args >( args )... ), AutomaticDifferenceOperatorType()
222 {}
223 };
224
225
226
227 // ModelDifferentiableGalerkinOperator
228 // -----------------------------------
229
230 template < class LinearOperator, class ModelIntegrands >
232 : public MOLDifferentiableGalerkinOperator< ModelIntegrands, LinearOperator >
233 {
235
236 typedef typename ModelIntegrands::ModelType ModelType;
237
239 typedef typename LinearOperator::RangeSpaceType DiscreteFunctionSpaceType;
240
242 : BaseType( dfSpace.gridPart(), model )
243 {}
244
245 template< class GridFunction >
246 void apply ( const GridFunction &u, RangeFunctionType &w ) const
247 {
248 (*this)( u, w );
249 }
250
251 template< class GridFunction >
252 void apply ( const GridFunction &u, LinearOperator &jOp ) const
253 {
254 (*this).jacobian( u, jOp );
255 }
256 };
257
258
259 // MethodOfLinesScheme
260 // -------------------
261
262 template< class Integrands, class LinearOperator, class InverseOperator, bool addDirichletBC>
263 using MethodOfLinesScheme = Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC,
265
266 } // namespace Fem
267
268} // namespace Dune
269
270#endif // #ifndef DUNE_FEM_SCHEMES_MOLGALERKIN_HH
STL namespace.
Definition: bindguard.hh:11
Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC, MOLDifferentiableGalerkinOperator > MethodOfLinesScheme
Definition: molgalerkin.hh:264
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
Definition: common/localcontribution.hh:14
void applyInverse(MassCaller &caller, const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:285
void leftMultiplyInverse(LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:339
Local Mass Matrix for arbitrary spaces.
Definition: localmassmatrix.hh:909
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
JacobianOperator JacobianOperatorType
type of linear operator modelling the operator's Jacobian
Definition: differentiableoperator.hh:35
BaseType::DomainFunctionType DomainFunctionType
type of discrete function in the operator's domain
Definition: differentiableoperator.hh:38
abstract operator
Definition: operator.hh:34
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
abstract affine-linear operator
Definition: operator.hh:87
Definition: molgalerkin.hh:19
bool communicate_
Definition: molgalerkin.hh:102
LocalMassMatrix< DiscreteFunctionSpaceType, InteriorQuadratureType > LocalMassMatrixType
Definition: molgalerkin.hh:34
ModelType & model() const
Definition: molgalerkin.hh:64
const GridPartType & gridPart() const
Definition: molgalerkin.hh:60
void setCommunicate(const bool communicate)
Definition: molgalerkin.hh:46
RangeFunction RangeFunctionType
Definition: molgalerkin.hh:21
void evaluate(const GridFunction &u, RangeFunctionType &w) const
Definition: molgalerkin.hh:87
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const final override
Definition: molgalerkin.hh:49
GalerkinOperatorImplType impl_
Definition: molgalerkin.hh:101
RangeFunctionType::GridPartType GridPartType
Definition: molgalerkin.hh:24
DomainFunction DomainFunctionType
Definition: molgalerkin.hh:20
Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType
Definition: molgalerkin.hh:28
RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: molgalerkin.hh:29
MOLGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:37
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: molgalerkin.hh:47
Integrands DirichletModelType
Definition: molgalerkin.hh:63
void applyInverseMass(RangeFunctionType &w) const
Definition: molgalerkin.hh:67
Integrands ModelType
Definition: molgalerkin.hh:62
GalerkinOperatorImplType::template QuadratureSelector< DiscreteFunctionSpaceType >::InteriorQuadratureType InteriorQuadratureType
Definition: molgalerkin.hh:32
Definition: molgalerkin.hh:114
MOLDifferentiableGalerkinOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, Args &&... args)
Definition: molgalerkin.hh:128
const RangeDiscreteFunctionSpaceType & rSpace_
Definition: molgalerkin.hh:200
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const final override
obtain linearization
Definition: molgalerkin.hh:136
const DomainDiscreteFunctionSpaceType & dSpace_
Definition: molgalerkin.hh:199
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: molgalerkin.hh:123
void applyInverseMass(JacobianOperatorType &jOp, const bool hasSkeleton) const
Definition: molgalerkin.hh:165
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: molgalerkin.hh:153
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: molgalerkin.hh:122
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:125
JacobianOperator JacobianOperatorType
Definition: molgalerkin.hh:118
BaseType::RangeFunctionType RangeFunctionType
Definition: molgalerkin.hh:121
BaseType::DomainFunctionType DomainFunctionType
Definition: molgalerkin.hh:120
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: molgalerkin.hh:145
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: molgalerkin.hh:157
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:217
MOLAutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:220
MOLModelDifferentiableGalerkinOperator(ModelType &model, const DiscreteFunctionSpaceType &dfSpace)
Definition: molgalerkin.hh:241
LinearOperator::DomainFunctionType RangeFunctionType
Definition: molgalerkin.hh:238
void apply(const GridFunction &u, RangeFunctionType &w) const
Definition: molgalerkin.hh:246
MOLDifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType
Definition: molgalerkin.hh:234
LinearOperator::RangeSpaceType DiscreteFunctionSpaceType
Definition: molgalerkin.hh:239
void apply(const GridFunction &u, LinearOperator &jOp) const
Definition: molgalerkin.hh:252
ModelIntegrands::ModelType ModelType
Definition: molgalerkin.hh:236