dune-fem  2.8-git
domainintegral.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_INTEGRAL_HH
2 #define DUNE_FEM_INTEGRAL_HH
3 
4 #include <type_traits>
5 
6 #include <dune/grid/common/rangegenerators.hh>
7 
10 
12 
13 #include <dune/common/bartonnackmanifcheck.hh>
15 
16 namespace Dune
17 {
18 
19  namespace Fem
20  {
21  // IntegralBase
22  // ----------
23 
24  template< class GridPart, class NormImplementation >
26  : public BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >,
27  NormImplementation >
28  {
30  NormImplementation > BaseType ;
32 
33  public:
34  typedef GridPart GridPartType;
35 
36  protected:
37  using BaseType :: asImp ;
38 
39  typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
40 
41  template <bool uDiscrete, bool vDiscrete>
43  {
44  template <class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet>
45  static ReturnType forEach ( const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
46  const ReturnType &initialValue,
47  const PartitionSet& partitionSet,
48  unsigned int order )
49  {
50  static_assert( uDiscrete && vDiscrete, "Distance can only be calculated between GridFunctions" );
51 
52  ReturnType sum( 0 );
55  for( const EntityType &entity : elements( norm.gridPart_, partitionSet ) )
56  {
57  uLocal.bind( entity );
58  vLocal.bind( entity );
59  const unsigned int uOrder = uLocal.order();
60  const unsigned int vOrder = vLocal.order();
61  const unsigned int orderLocal = (order == 0 ? 2*std::max( uOrder, vOrder ) : order);
62  norm.distanceLocal( entity, orderLocal, uLocal, vLocal, sum );
63  uLocal.unbind();
64  vLocal.unbind();
65  }
66  return sum;
67  }
68  };
69 
70  // this specialization creates a grid function adapter
71  template <bool vDiscrete>
72  struct ForEachCaller<false, vDiscrete>
73  {
74  template <class F,
75  class VDiscreteFunctionType,
76  class ReturnType,
77  class PartitionSet>
78  static
79  ReturnType forEach ( const ThisType& norm,
80  const F& f, const VDiscreteFunctionType& v,
81  const ReturnType& initialValue,
82  const PartitionSet& partitionSet,
83  const unsigned int order )
84  {
85  typedef GridFunctionAdapter< F, GridPartType> GridFunction ;
86  GridFunction u( "Integral::adapter" , f , v.space().gridPart(), v.space().order() );
87 
89  forEach( norm, u, v, initialValue, partitionSet, order );
90  }
91  };
92 
93  // this specialization simply switches arguments
94  template <bool uDiscrete>
95  struct ForEachCaller<uDiscrete, false>
96  {
97  template <class UDiscreteFunctionType,
98  class F,
99  class ReturnType,
100  class PartitionSet>
101  static
102  ReturnType forEach ( const ThisType& norm,
103  const UDiscreteFunctionType& u,
104  const F& f,
105  const ReturnType& initialValue,
106  const PartitionSet& partitionSet,
107  const unsigned int order )
108  {
110  forEach( norm, f, u, initialValue, partitionSet, order );
111  }
112  };
113 
114  template< class DiscreteFunctionType, class ReturnType, class PartitionSet >
115  ReturnType forEach ( const DiscreteFunctionType &u, const ReturnType &initialValue,
116  const PartitionSet& partitionSet,
117  unsigned int order = 0 ) const
118  {
119  static_assert( (std::is_base_of<Fem::HasLocalFunction, DiscreteFunctionType>::value),
120  "Norm only implemented for quantities implementing a local function!" );
121 
122  ReturnType sum( 0 );
124  for( const EntityType &entity : elements( gridPart_, partitionSet ) )
125  {
126  uLocal.bind( entity );
127  const unsigned int orderLocal = (order == 0 ? 2*uLocal.order() : order);
128  normLocal( entity, orderLocal, uLocal, sum );
129  uLocal.unbind();
130  }
131  return sum;
132  }
133 
134  template< class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet >
135  ReturnType forEach ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
136  const ReturnType &initialValue, const PartitionSet& partitionSet,
137  unsigned int order = 0 ) const
138  {
139  enum { uDiscrete = std::is_convertible<UDiscreteFunctionType, HasLocalFunction>::value };
140  enum { vDiscrete = std::is_convertible<VDiscreteFunctionType, HasLocalFunction>::value };
141 
142  // call forEach depending on which argument is a grid function,
143  // i.e. has a local function
145  forEach( *this, u, v, initialValue, partitionSet, order );
146  }
147 
148  public:
149  explicit IntegralBase ( const GridPartType &gridPart )
150  : gridPart_( gridPart )
151  {}
152 
153  protected:
154  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
155  void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
156  {
157  CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().distanceLocal( entity, order, uLocal, vLocal, sum ) );
158  }
159 
160  template< class LocalFunctionType, class ReturnType >
161  void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
162  {
163  CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().normLocal( entity, order, uLocal, sum ) );
164  }
165 
166  const GridPartType &gridPart () const { return gridPart_; }
167 
168  const typename GridPartType::CollectiveCommunicationType& comm () const
169  {
170  return gridPart().comm();
171  }
172 
173  bool checkCommunicateFlag( bool communicate ) const
174  {
175 #ifndef NDEBUG
176  bool commMax = communicate;
177  assert( communicate == comm().max( commMax ) );
178 #endif
179  return communicate;
180  }
181 
182  private:
183  const GridPartType &gridPart_;
184  };
185 
186  // Integral
187  // ------
188 
189  template< class GridPart >
190  class Integral : public IntegralBase< GridPart, Integral< GridPart > >
191  {
194 
195  public:
196  typedef GridPart GridPartType;
197 
198  using BaseType :: gridPart ;
199  using BaseType :: comm ;
200 
201  protected:
202 
203  template< class UFunction, class VFunction >
204  struct FunctionDistance;
205 
208 
209  const unsigned int order_;
210  const bool communicate_;
211  public:
217  explicit Integral ( const GridPartType &gridPart,
218  const unsigned int order = 0,
219  const bool communicate = true );
220 
221 
223  template< class DiscreteFunctionType, class PartitionSet >
224  typename DiscreteFunctionType::RangeType
225  norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const;
226 
228  template< class DiscreteFunctionType >
229  typename DiscreteFunctionType::RangeType
230  norm ( const DiscreteFunctionType &u ) const
231  {
232  return norm( u, Partitions::interior );
233  }
234 
236  template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
237  typename UDiscreteFunctionType::RangeType
238  distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const;
239 
241  template< class UDiscreteFunctionType, class VDiscreteFunctionType >
242  typename UDiscreteFunctionType::RangeType
243  distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const
244  {
245  return distance( u, v, Partitions::interior );
246  }
247 
248  template< class LocalFunctionType, class ReturnType >
249  void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
250 
251  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
252  void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
253  };
254 
255 
256 
257  // Implementation of Integral
258  // ------------------------
259 
260  template< class GridPart >
261  inline Integral< GridPart >::Integral ( const GridPartType &gridPart, const unsigned int order, const bool communicate )
262  : BaseType( gridPart ),
263  order_( order ),
264  communicate_( BaseType::checkCommunicateFlag( communicate ) )
265  {}
266 
267 
268  template< class GridPart >
269  template< class DiscreteFunctionType, class PartitionSet >
270  inline typename DiscreteFunctionType::RangeType
271  Integral< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const
272  {
273  typedef typename DiscreteFunctionType::RangeType RangeType;
274  typedef RangeType ReturnType ;
275 
276  // calculate integral over each element
277  ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ );
278 
279  // communicate_ indicates global norm
280  if( communicate_ )
281  {
282  sum = comm().sum( sum );
283  }
284 
285  return sum;
286  }
287 
288 
289  template< class GridPart >
290  template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
291  inline typename UDiscreteFunctionType::RangeType
293  ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const
294  {
295  typedef typename UDiscreteFunctionType::RangeType RangeType;
296  typedef RangeType ReturnType ;
297 
298  // calculate integral over each element
299  ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ );
300 
301  // communicate_ indicates global norm
302  if( communicate_ )
303  {
304  sum = comm().sum( sum );
305  }
306 
307  return sum;
308  }
309 
310  template< class GridPart >
311  template< class LocalFunctionType, class ReturnType >
312  inline void
313  Integral< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
314  {
315  Integrator< QuadratureType > integrator( order );
316 
317  integrator.integrateAdd( entity, uLocal, sum );
318  }
319 
320  template< class GridPart >
321  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
322  inline void
323  Integral< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
324  {
325  Integrator< QuadratureType > integrator( order );
326 
328 
329  LocalDistanceType dist( uLocal, vLocal );
330 
331  integrator.integrateAdd( entity, dist, sum );
332  }
333 
334  template< class GridPart >
335  template< class UFunction, class VFunction >
336  struct Integral< GridPart >::FunctionDistance
337  {
338  typedef UFunction UFunctionType;
339  typedef VFunction VFunctionType;
340 
341  typedef typename UFunctionType::RangeFieldType RangeFieldType;
342  typedef typename UFunctionType::RangeType RangeType;
343  typedef typename UFunctionType::JacobianRangeType JacobianRangeType;
344 
346  : u_( u ), v_( v )
347  {}
348 
349  template< class Point >
350  void evaluate ( const Point &x, RangeType &ret ) const
351  {
352  RangeType phi;
353  u_.evaluate( x, ret );
354  v_.evaluate( x, phi );
355  ret -= phi;
356  }
357 
358  template< class Point >
359  void jacobian ( const Point &x, JacobianRangeType &ret ) const
360  {
361  JacobianRangeType phi;
362  u_.jacobian( x, ret );
363  v_.jacobian( x, phi );
364  ret -= phi;
365  }
366 
367  private:
368  const UFunctionType &u_;
369  const VFunctionType &v_;
370  };
371 
372  } // namespace Fem
373 
374 } // namespace Dune
375 
376 #endif // #ifndef DUNE_FEM_INTEGRAL_HH
Definition: bindguard.hh:11
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:517
static double max(const Double &v, const double p)
Definition: double.hh:398
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr std::decay_t< T > sum(T a)
Definition: utility.hh:33
BasicGridFunctionAdapter provides local functions for a Function.
Definition: gridfunctionadapter.hh:79
Definition: bartonnackmaninterface.hh:17
const NormImplementation & asImp() const
Definition: bartonnackmaninterface.hh:37
Definition: domainintegral.hh:28
GridPart GridPartType
Definition: domainintegral.hh:34
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:155
IntegralBase(const GridPartType &gridPart)
Definition: domainintegral.hh:149
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:161
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: domainintegral.hh:39
bool checkCommunicateFlag(bool communicate) const
Definition: domainintegral.hh:173
ReturnType forEach(const DiscreteFunctionType &u, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:115
const GridPartType & gridPart() const
Definition: domainintegral.hh:166
const GridPartType::CollectiveCommunicationType & comm() const
Definition: domainintegral.hh:168
ReturnType forEach(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:135
Definition: domainintegral.hh:43
static ReturnType forEach(const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order)
Definition: domainintegral.hh:45
static ReturnType forEach(const ThisType &norm, const F &f, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, const unsigned int order)
Definition: domainintegral.hh:79
static ReturnType forEach(const ThisType &norm, const UDiscreteFunctionType &u, const F &f, const ReturnType &initialValue, const PartitionSet &partitionSet, const unsigned int order)
Definition: domainintegral.hh:102
Definition: domainintegral.hh:191
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u) const
|| u ||_L1 on interior partition entities
Definition: domainintegral.hh:230
CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: domainintegral.hh:207
GridPart GridPartType
Definition: domainintegral.hh:196
const bool communicate_
Definition: domainintegral.hh:210
UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v) const
|| u - v ||_L2 on interior partition entities
Definition: domainintegral.hh:243
BaseType::EntityType EntityType
Definition: domainintegral.hh:204
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u, const PartitionSet &partitionSet) const
|| u ||_L1 on given set of entities (partition set)
Definition: domainintegral.hh:271
Integral(const GridPartType &gridPart, const unsigned int order=0, const bool communicate=true)
constructor
Definition: domainintegral.hh:261
UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet &partitionSet) const
|| u - v ||_L2 on given set of entities (partition set)
Definition: domainintegral.hh:293
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:313
const unsigned int order_
Definition: domainintegral.hh:209
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:323
Definition: domainintegral.hh:337
UFunctionType::RangeType RangeType
Definition: domainintegral.hh:342
void evaluate(const Point &x, RangeType &ret) const
Definition: domainintegral.hh:350
UFunction UFunctionType
Definition: domainintegral.hh:338
FunctionDistance(const UFunctionType &u, const VFunctionType &v)
Definition: domainintegral.hh:345
UFunctionType::RangeFieldType RangeFieldType
Definition: domainintegral.hh:341
VFunction VFunctionType
Definition: domainintegral.hh:339
void jacobian(const Point &x, JacobianRangeType &ret) const
Definition: domainintegral.hh:359
UFunctionType::JacobianRangeType JacobianRangeType
Definition: domainintegral.hh:343
BaseType::EntityType EntityType
Definition: lpnorm.hh:143
integrator for arbitrary functions providing evaluate
Definition: integrator.hh:28
void integrateAdd(const EntityType &entity, const Function &function, typename Function ::RangeType &ret) const
add the integral over an entity to a variable
Definition: integrator.hh:67