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 
7 namespace 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
93  applyInverseMass( w );
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 >
220  explicit MOLAutomaticDifferenceGalerkinOperator ( const GridPartType &gridPart, Args &&... 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
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
BaseType::RangeFunctionType RangeFunctionType
type of discrete function in the operator's range
Definition: differentiableoperator.hh:40
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
RangeFunctionType ::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: molgalerkin.hh:29
bool communicate_
Definition: molgalerkin.hh:102
LocalMassMatrix< DiscreteFunctionSpaceType, InteriorQuadratureType > LocalMassMatrixType
Definition: molgalerkin.hh:34
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
GalerkinOperatorImplType::template QuadratureSelector< DiscreteFunctionSpaceType >::InteriorQuadratureType InteriorQuadratureType
Definition: molgalerkin.hh:32
RangeFunctionType::GridPartType GridPartType
Definition: molgalerkin.hh:24
DomainFunction DomainFunctionType
Definition: molgalerkin.hh:20
Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType
Definition: molgalerkin.hh:28
MOLGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:37
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: molgalerkin.hh:47
const GridPartType & gridPart() const
Definition: molgalerkin.hh:60
Integrands DirichletModelType
Definition: molgalerkin.hh:63
void applyInverseMass(RangeFunctionType &w) const
Definition: molgalerkin.hh:67
Integrands ModelType
Definition: molgalerkin.hh:62
ModelType & model() const
Definition: molgalerkin.hh:64
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
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: molgalerkin.hh:153
void applyInverseMass(JacobianOperatorType &jOp, const bool hasSkeleton) const
Definition: molgalerkin.hh:165
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: molgalerkin.hh:122
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:125
JacobianOperator JacobianOperatorType
Definition: molgalerkin.hh:118
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: molgalerkin.hh:157
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
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