1 #ifndef DUNE_FEM_LPNORM_HH
2 #define DUNE_FEM_LPNORM_HH
3 
4 #include <dune/fem/quadrature/integrator.hh>
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 
22     //! Quadrature Order Interface
23     struct OrderCalculatorInterface
24     {
25       virtual int operator() (const double p)=0;
26     };
27 
28     //! default Quadrature Order Calculator
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
32     struct DefaultOrderCalculator : public OrderCalculatorInterface
33     {
operator ()Dune::Fem::DefaultOrderCalculator34       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     {
48       typedef LPNorm< GridPart, OrderCalculator > ThisType;
49       typedef IntegralBase< GridPart, ThisType > BaseType ;
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;
65       typedef CachingQuadrature< GridPartType, 0 > QuadratureType;
66       typedef Integrator< QuadratureType > IntegratorType;
67 
68     public:
69       /** \brief constructor
70        *    \param gridPart     specific gridPart for selection of entities
71        *    \param p            p in Lp-norm
72        *    \param communicate  if true global (over all ranks) norm is computed (default = true)
73        */
74       explicit LPNorm ( const GridPartType &gridPart, const double p, const bool communicate = true  );
75 
76       //! || u ||_Lp on given set of entities (partition set)
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 
81       //! || u ||_Lp on interior partition entities
82       template< class DiscreteFunctionType >
83       typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
norm(const DiscreteFunctionType & u) const84       norm ( const DiscreteFunctionType &u ) const
85       {
86         return norm( u, Partitions::interior );
87       }
88 
89       //! || u - v ||_Lp on given set of entities (partition set)
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 
94       //! || u - v ||_Lp on interior partition entities
95       template< class UDiscreteFunctionType, class VDiscreteFunctionType >
96       typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v) const97       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 >
123     class WeightedLPNorm
124     : public LPNorm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType,
125                      OrderCalculator >
126     {
127       typedef WeightedLPNorm< WeightFunction, OrderCalculator > ThisType;
128       typedef LPNorm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType, OrderCalculator> BaseType;
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 >
138       struct WeightedFunctionMultiplicator;
139 
140       typedef ConstLocalFunction< WeightFunctionType > LocalWeightFunctionType;
141       typedef typename WeightFunctionType::RangeType WeightType;
142 
143       typedef typename BaseType::EntityType EntityType;
144       typedef typename BaseType::IntegratorType IntegratorType;
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 >
LPNorm(const GridPartType & gridPart,const double p,const bool communicate)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>
order(const int spaceOrder) const180     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
norm(const DiscreteFunctionType & u,const PartitionSet & partitionSet) const189     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
211     LPNorm< GridPart, OrderCalculator >
distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v,const PartitionSet & partitionSet) const212       ::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
normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const234     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
distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const246     LPNorm< GridPart, OrderCalculator >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
247     {
248       typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType;
249 
250       IntegratorType integrator( order );
251 
252       LocalDistanceType dist( uLocal, vLocal );
253       FunctionMultiplicator< LocalDistanceType > distp( dist, p_ );
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     {
263       typedef Function FunctionType;
264 
265       typedef typename FunctionType::RangeFieldType RangeFieldType;
266       typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
267       typedef FieldVector< RealType, 1 > RangeType ;
268 
FunctionMultiplicatorDune::Fem::LPNorm::FunctionMultiplicator269       explicit FunctionMultiplicator ( const FunctionType &function, double p )
270       : function_( function ),
271         p_(p)
272       {}
273 
274       template< class Point >
evaluateDune::Fem::LPNorm::FunctionMultiplicator275       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 
FunctionDistanceDune::Fem::LPNorm::FunctionDistance300       FunctionDistance ( const UFunctionType &u, const VFunctionType &v )
301       : u_( u ), v_( v )
302       {}
303 
304       template< class Point >
evaluateDune::Fem::LPNorm::FunctionDistance305       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 >
jacobianDune::Fem::LPNorm::FunctionDistance314       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 >
332     inline WeightedLPNorm< WeightFunction, OrderCalculator >
WeightedLPNorm(const WeightFunctionType & weightFunction,double p,const bool communicate)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
normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const346     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
distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const361     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 );
370       WeightedFunctionMultiplicator< LocalDistanceType > dist2( wfLocal_, dist );
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     {
381       typedef Function FunctionType;
382 
383       typedef typename FunctionType::RangeFieldType RangeFieldType;
384       typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
385       typedef FieldVector< RealType, 1 > RangeType;
386 
WeightedFunctionMultiplicatorDune::Fem::WeightedLPNorm::WeightedFunctionMultiplicator387       WeightedFunctionMultiplicator ( const LocalWeightFunctionType &weightFunction,
388                                       const FunctionType &function,
389                                       double p )
390       : weightFunction_( weightFunction ),
391         function_( function ),
392         p_(p)
393       {}
394 
395       template< class Point >
evaluateDune::Fem::WeightedLPNorm::WeightedFunctionMultiplicator396       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
417