dune-fem  2.8-git
lpnorm.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_LPNORM_HH
2 #define DUNE_FEM_LPNORM_HH
3 
5 
6 #include "domainintegral.hh"
7 
8 namespace Dune
9 {
10 
11  namespace Fem
12  {
13 
14  // TODO weighte LP norm might be adapted later
15  // LPNorm
16  //
17  // !!!!! It is not cleared which quadrature order have to be applied for p > 2!!!!
18  // !!!!! For p = 1 this norm does not work !!!!
19  // ------
20 
21 
24  {
25  virtual int operator() (const double p)=0;
26  };
27 
29  // can be re-implemented in order to use
30  // a different type of calculation
31  // which can be sepcified in the second template argument of LPNorm
33  {
34  int operator() (const double p)
35  {
36  int ret=0;
37  const double q = p / (p-1);
38  double max = std::max(p,q);
39  assert(max < std::numeric_limits<int>::max()/2. );
40  ret = max +1;
41  return ret;
42  }
43  };
44 
45  template< class GridPart, class OrderCalculator = DefaultOrderCalculator >
46  class LPNorm : public IntegralBase < GridPart, LPNorm< GridPart, OrderCalculator > >
47  {
50 
51  public:
52  typedef GridPart GridPartType;
53 
54  using BaseType::gridPart;
55  using BaseType::comm;
56 
57  protected:
58  template< class Function >
59  struct FunctionMultiplicator;
60 
61  template< class UFunction, class VFunction >
62  struct FunctionDistance;
63 
64  typedef typename BaseType::EntityType EntityType;
67 
68  public:
74  explicit LPNorm ( const GridPartType &gridPart, const double p, const bool communicate = true );
75 
77  template< class DiscreteFunctionType, class PartitionSet >
78  typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
79  norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const;
80 
82  template< class DiscreteFunctionType >
83  typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
84  norm ( const DiscreteFunctionType &u ) const
85  {
86  return norm( u, Partitions::interior );
87  }
88 
90  template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
91  typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
92  distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const;
93 
95  template< class UDiscreteFunctionType, class VDiscreteFunctionType >
96  typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
97  distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const
98  {
99  return distance( u, v, Partitions::interior );
100  }
101 
102  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
103  void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
104 
105  template< class LocalFunctionType, class ReturnType >
106  void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
107 
108  int order ( const int spaceOrder ) const ;
109 
110  protected:
111  double p_ ;
112  OrderCalculator *orderCalc_;
113  const bool communicate_;
114  };
115 
116 
117 
118 
119  // WeightedLPNorm
120  // --------------
121 
122  template< class WeightFunction, class OrderCalculator = DefaultOrderCalculator >
124  : public LPNorm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType,
125  OrderCalculator >
126  {
129 
130  public:
131  typedef WeightFunction WeightFunctionType;
132 
133  typedef typename WeightFunctionType::DiscreteFunctionSpaceType WeightFunctionSpaceType;
134  typedef typename WeightFunctionSpaceType::GridPartType GridPartType;
135 
136  protected:
137  template< class Function >
139 
141  typedef typename WeightFunctionType::RangeType WeightType;
142 
145 
146  using BaseType::gridPart;
147  using BaseType::comm;
148 
149  public:
150  using BaseType::norm;
151  using BaseType::distance;
152 
153  explicit WeightedLPNorm ( const WeightFunctionType &weightFunction, const double p, const bool communicate = true );
154 
155  template< class LocalFunctionType, class ReturnType >
156  void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
157 
158  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
159  void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
160 
161  private:
162  LocalWeightFunctionType wfLocal_;
163  const double p_;
164  };
165 
166 
167  // Implementation of LPNorm
168  // ------------------------
169 
170  template< class GridPart, class OrderCalculator >
171  inline LPNorm< GridPart, OrderCalculator >::LPNorm ( const GridPartType &gridPart, const double p, const bool communicate )
172  : BaseType( gridPart ),
173  p_(p),
174  orderCalc_( new OrderCalculator() ),
175  communicate_( BaseType::checkCommunicateFlag( communicate ) )
176  {
177  }
178 
179  template< class GridPart, class OrderCalculator>
180  inline int LPNorm< GridPart, OrderCalculator>::order(const int spaceOrder) const
181  {
182  return spaceOrder * orderCalc_->operator() (p_);
183  }
184 
185 
186  template< class GridPart, class OrderCalculator>
187  template< class DiscreteFunctionType, class PartitionSet >
188  inline typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
189  LPNorm< GridPart, OrderCalculator >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const
190  {
191  typedef typename DiscreteFunctionType::RangeFieldType RangeFieldType;
192  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
193  typedef FieldVector< RealType, 1 > ReturnType ;
194 
195  // calculate integral over each element
196  ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet );
197 
198  // communicate_ indicates global norm
199  if( communicate_ )
200  {
201  sum[ 0 ] = comm().sum( sum[ 0 ] );
202  }
203 
204  // return result
205  return std::pow ( sum[ 0 ], (1.0 / p_) );
206  }
207 
208  template< class GridPart, class OrderCalculator >
209  template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
210  inline typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
212  ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const
213  {
214  typedef typename UDiscreteFunctionType::RangeFieldType RangeFieldType;
215  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
216  typedef FieldVector< RealType, 1 > ReturnType ;
217 
218  // calculate integral over each element
219  ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet );
220 
221  // communicate_ indicates global norm
222  if( communicate_ )
223  {
224  sum[ 0 ] = comm().sum( sum[ 0 ] );
225  }
226 
227  // return result
228  return std::pow( sum[ 0 ], (1.0/p_) );
229  }
230 
231  template< class GridPart, class OrderCalculator >
232  template< class LocalFunctionType, class ReturnType >
233  inline void
234  LPNorm< GridPart, OrderCalculator >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
235  {
236  IntegratorType integrator( order );
237 
238  FunctionMultiplicator< LocalFunctionType > uLocalp( uLocal, p_ );
239 
240  integrator.integrateAdd( entity, uLocalp, sum );
241  }
242 
243  template< class GridPart, class OrderCalculator >
244  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
245  inline void
246  LPNorm< GridPart, OrderCalculator >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
247  {
249 
250  IntegratorType integrator( order );
251 
252  LocalDistanceType dist( uLocal, vLocal );
254 
255  integrator.integrateAdd( entity, distp, sum );
256  }
257 
258 
259  template< class GridPart, class OrderCalculator >
260  template< class Function >
261  struct LPNorm< GridPart, OrderCalculator >::FunctionMultiplicator
262  {
264 
266  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
267  typedef FieldVector< RealType, 1 > RangeType ;
268 
269  explicit FunctionMultiplicator ( const FunctionType &function, double p )
270  : function_( function ),
271  p_(p)
272  {}
273 
274  template< class Point >
275  void evaluate ( const Point &x, RangeType &ret ) const
276  {
277  typename FunctionType::RangeType phi;
278  function_.evaluate( x, phi );
279  ret = std :: pow ( phi.two_norm(), p_);
280  }
281 
282  private:
283  const FunctionType &function_;
284  double p_;
285  };
286 
287 
288  template< class GridPart, class OrderCalculator >
289  template< class UFunction, class VFunction >
290  struct LPNorm< GridPart, OrderCalculator >::FunctionDistance
291  {
292  typedef UFunction UFunctionType;
293  typedef VFunction VFunctionType;
294 
295  typedef typename UFunctionType::RangeFieldType RangeFieldType;
296  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
297  typedef typename UFunctionType::RangeType RangeType;
298  typedef typename UFunctionType::JacobianRangeType JacobianRangeType;
299 
301  : u_( u ), v_( v )
302  {}
303 
304  template< class Point >
305  void evaluate ( const Point &x, RangeType &ret ) const
306  {
307  RangeType phi;
308  u_.evaluate( x, ret );
309  v_.evaluate( x, phi );
310  ret -= phi;
311  }
312 
313  template< class Point >
314  void jacobian ( const Point &x, JacobianRangeType &ret ) const
315  {
316  JacobianRangeType phi;
317  u_.jacobian( x, ret );
318  v_.jacobian( x, phi );
319  ret -= phi;
320  }
321 
322  private:
323  const UFunctionType &u_;
324  const VFunctionType &v_;
325  };
326 
327 
328  // Implementation of WeightedL2Norm
329  // --------------------------------
330 
331  template< class WeightFunction, class OrderCalculator >
333  ::WeightedLPNorm ( const WeightFunctionType &weightFunction, double p, const bool communicate )
334  : BaseType( weightFunction.space().gridPart(), p, communicate ),
335  wfLocal_( weightFunction ),
336  p_(p)
337  {
338  static_assert( (WeightFunctionSpaceType::dimRange == 1),
339  "Weight function must be scalar." );
340  }
341 
342 
343  template< class WeightFunction, class OrderCalculator >
344  template< class LocalFunctionType, class ReturnType >
345  inline void
346  WeightedLPNorm< WeightFunction, OrderCalculator >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
347  {
348  // !!!! order !!!!
349  IntegratorType integrator( order );
350 
351  wfLocal_.bind( entity );
352  WeightedFunctionMultiplicator< LocalFunctionType > uLocal2( wfLocal_, uLocal, p_ );
353  integrator.integrateAdd( entity, uLocal2, sum );
354  wfLocal_.unbind();
355  }
356 
357 
358  template< class WeightFunction,class OrderCalculator >
359  template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
360  inline void
361  WeightedLPNorm< WeightFunction, OrderCalculator >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
362  {
363  typedef typename BaseType::template FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType;
364 
365  // !!!! order !!!!
366  IntegratorType integrator( order );
367 
368  wfLocal_.bind( entity );
369  LocalDistanceType dist( uLocal, vLocal );
371 
372  integrator.integrateAdd( entity, dist2, sum );
373  wfLocal_.unbind();
374  }
375 
376 
377  template< class WeightFunction, class OrderCalculator>
378  template< class Function >
379  struct WeightedLPNorm< WeightFunction, OrderCalculator>::WeightedFunctionMultiplicator
380  {
382 
384  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
385  typedef FieldVector< RealType, 1 > RangeType;
386 
388  const FunctionType &function,
389  double p )
390  : weightFunction_( weightFunction ),
391  function_( function ),
392  p_(p)
393  {}
394 
395  template< class Point >
396  void evaluate ( const Point &x, RangeType &ret ) const
397  {
398  WeightType weight;
399  weightFunction_.evaluate( x, weight );
400 
401  typename FunctionType::RangeType phi;
402  function_.evaluate( x, phi );
403  ret = weight[ 0 ] * std::pow ( phi.two_norm(), p_);
404  }
405 
406  private:
407  const LocalWeightFunctionType &weightFunction_;
408  const FunctionType &function_;
409  double p_;
410  };
411 
412  } // namespace Fem
413 
414 } // namespace Dune
415 
416 #endif // #ifndef DUNE_FEM_LPNORM_HH
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
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
Abstract class representing a function.
Definition: common/function.hh:50
FunctionSpaceType ::RangeType RangeType
range type
Definition: common/function.hh:68
FunctionSpaceType ::RangeFieldType RangeFieldType
field type of range
Definition: common/function.hh:64
Definition: domainintegral.hh:28
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: domainintegral.hh:39
const GridPartType & gridPart() const
Definition: domainintegral.hh:166
const GridPartType::CollectiveCommunicationType & comm() const
Definition: domainintegral.hh:168
Quadrature Order Interface.
Definition: lpnorm.hh:24
virtual int operator()(const double p)=0
default Quadrature Order Calculator
Definition: lpnorm.hh:33
int operator()(const double p)
Definition: lpnorm.hh:34
Definition: lpnorm.hh:47
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: lpnorm.hh:234
const bool communicate_
Definition: lpnorm.hh:113
LPNorm(const GridPartType &gridPart, const double p, const bool communicate=true)
constructor
Definition: lpnorm.hh:171
BaseType::EntityType EntityType
Definition: lpnorm.hh:62
Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType &u) const
|| u ||_Lp on interior partition entities
Definition: lpnorm.hh:84
GridPart GridPartType
Definition: lpnorm.hh:52
OrderCalculator * orderCalc_
Definition: lpnorm.hh:112
Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet &partitionSet) const
|| u - v ||_Lp on given set of entities (partition set)
Definition: lpnorm.hh:212
int order(const int spaceOrder) const
Definition: lpnorm.hh:180
Integrator< QuadratureType > IntegratorType
Definition: lpnorm.hh:66
double p_
Definition: lpnorm.hh:111
const GridPartType & gridPart() const
Definition: domainintegral.hh:166
Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v) const
|| u - v ||_Lp on interior partition entities
Definition: lpnorm.hh:97
const GridPartType::CollectiveCommunicationType & comm() const
Definition: domainintegral.hh:168
CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: lpnorm.hh:65
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: lpnorm.hh:246
Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType &u, const PartitionSet &partitionSet) const
|| u ||_Lp on given set of entities (partition set)
Definition: lpnorm.hh:189
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: lpnorm.hh:266
FunctionType::RangeFieldType RangeFieldType
Definition: lpnorm.hh:265
void evaluate(const Point &x, RangeType &ret) const
Definition: lpnorm.hh:275
FunctionMultiplicator(const FunctionType &function, double p)
Definition: lpnorm.hh:269
FieldVector< RealType, 1 > RangeType
Definition: lpnorm.hh:267
Function FunctionType
Definition: lpnorm.hh:263
Definition: lpnorm.hh:291
void evaluate(const Point &x, RangeType &ret) const
Definition: lpnorm.hh:305
UFunction UFunctionType
Definition: lpnorm.hh:292
UFunctionType::RangeFieldType RangeFieldType
Definition: lpnorm.hh:295
UFunctionType::RangeType RangeType
Definition: lpnorm.hh:297
VFunction VFunctionType
Definition: lpnorm.hh:293
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: lpnorm.hh:296
void jacobian(const Point &x, JacobianRangeType &ret) const
Definition: lpnorm.hh:314
FunctionDistance(const UFunctionType &u, const VFunctionType &v)
Definition: lpnorm.hh:300
UFunctionType::JacobianRangeType JacobianRangeType
Definition: lpnorm.hh:298
Definition: lpnorm.hh:126
BaseType::EntityType EntityType
Definition: lpnorm.hh:143
WeightFunctionType::DiscreteFunctionSpaceType WeightFunctionSpaceType
Definition: lpnorm.hh:133
WeightFunction WeightFunctionType
Definition: lpnorm.hh:131
ConstLocalFunction< WeightFunctionType > LocalWeightFunctionType
Definition: lpnorm.hh:138
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: lpnorm.hh:346
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: lpnorm.hh:361
WeightedLPNorm(const WeightFunctionType &weightFunction, const double p, const bool communicate=true)
Definition: lpnorm.hh:333
WeightFunctionType::RangeType WeightType
Definition: lpnorm.hh:141
WeightFunctionSpaceType::GridPartType GridPartType
Definition: lpnorm.hh:134
BaseType::IntegratorType IntegratorType
Definition: lpnorm.hh:144
WeightedFunctionMultiplicator(const LocalWeightFunctionType &weightFunction, const FunctionType &function, double p)
Definition: lpnorm.hh:387
FunctionType::RangeFieldType RangeFieldType
Definition: lpnorm.hh:383
void evaluate(const Point &x, RangeType &ret) const
Definition: lpnorm.hh:396
FieldVector< RealType, 1 > RangeType
Definition: lpnorm.hh:385
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: lpnorm.hh:384
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