DUNE-FEM (unstable)

lpnorm.hh
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
8namespace 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
23 struct OrderCalculatorInterface
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
32 struct DefaultOrderCalculator : public OrderCalculatorInterface
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 {
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:
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 >
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 >
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
211 LPNorm< GridPart, OrderCalculator >
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 {
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
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
300 FunctionDistance ( const UFunctionType &u, const VFunctionType &v )
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 >
332 inline WeightedLPNorm< WeightFunction, 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 );
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
387 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 >
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
Dune::Fem::CachingQuadratureImpl
quadrature class supporting base function caching
Definition: cachingquadrature.hh:101
Dune::Fem::DiscreteFunctionDefault< TupleDiscreteFunction< DiscreteFunctions ... > >::RangeFieldType
DiscreteFunctionSpaceType::RangeFieldType RangeFieldType
type of range field (usually a float type)
Definition: discretefunction.hh:623
Dune::Fem::Integrator
integrator for arbitrary functions providing evaluate
Definition: integrator.hh:28
Dune::Fem::TupleDiscreteFunction
forward declaration
Definition: discretefunction.hh:51
Dune::Hybrid::forEach
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
Dune::Hybrid::max
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::Fem::DefaultOrderCalculator
default Quadrature Order Calculator
Definition: lpnorm.hh:33
Dune::Fem::OrderCalculatorInterface
Quadrature Order Interface.
Definition: lpnorm.hh:24
Dune::PartitionSet
A set of PartitionType values.
Definition: partitionset.hh:137
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)