DUNE-FEM (unstable)

l2norm.hh
1#ifndef DUNE_FEM_L2NORM_HH
2#define DUNE_FEM_L2NORM_HH
3
4#include <dune/fem/quadrature/cachingquadrature.hh>
5#include <dune/fem/quadrature/integrator.hh>
6
7#include <dune/fem/misc/domainintegral.hh>
8
9namespace Dune
10{
11
12 namespace Fem
13 {
14
15 // L2Norm
16 // ------
17
18 template< class GridPart >
19 class L2Norm : public IntegralBase< GridPart, L2Norm< GridPart > >
20 {
21 typedef IntegralBase< GridPart, L2Norm< GridPart > > BaseType ;
22 typedef L2Norm< GridPart > ThisType;
23
24 public:
25 typedef GridPart GridPartType;
26
27 using BaseType :: gridPart ;
28 using BaseType :: comm ;
29
30 template< class Function >
31 struct FunctionSquare;
32
33 template< class UFunction, class VFunction >
34 struct FunctionDistance;
35
36 protected:
37 typedef typename BaseType::EntityType EntityType;
38 typedef CachingQuadrature< GridPartType, 0 > QuadratureType;
39 typedef Integrator< QuadratureType > IntegratorType;
40
41 const unsigned int order_;
42 const bool communicate_;
43 public:
49 explicit L2Norm ( const GridPartType &gridPart,
50 const unsigned int order = 0,
51 const bool communicate = true );
52
54 template< class DiscreteFunctionType, class PartitionSet >
55 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
56 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const;
57
59 template< class DiscreteFunctionType >
60 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
61 norm ( const DiscreteFunctionType &u ) const
62 {
63 return norm( u, Partitions::interior );
64 }
65
67 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
68 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
69 distance ( const UDiscreteFunctionType &u,
70 const VDiscreteFunctionType &v,
71 const PartitionSet& partitionSet ) const;
72
74 template< class UDiscreteFunctionType, class VDiscreteFunctionType >
75 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
76 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const
77 {
78 return distance( u, v, Partitions::interior );
79 }
80
81 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
82 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
83
84 template< class LocalFunctionType, class ReturnType >
85 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
86 };
87
88
89 // Implementation of L2Norm
90 // ------------------------
91
92 template< class GridPart >
93 inline L2Norm< GridPart >::L2Norm ( const GridPartType &gridPart, const unsigned int order, const bool communicate )
94 : BaseType( gridPart ),
95 order_( order ),
96 communicate_( BaseType::checkCommunicateFlag( communicate ) )
97 {
98 }
99
100
101 template< class GridPart >
102 template< class DiscreteFunctionType, class PartitionSet >
103 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
104 L2Norm< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const
105 {
106 typedef typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type RealType;
107 typedef FieldVector< RealType, 1 > ReturnType ;
108
109 // calculate integral over each element
110 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ );
111
112 // communicate_ indicates global norm
113 if( communicate_ )
114 {
115 sum[ 0 ] = comm().sum( sum[ 0 ] );
116 }
117
118 // return result, e.g. sqrt of calculated sum
119 return std::sqrt( sum[ 0 ] );
120 }
121
122
123 template< class GridPart >
124 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
125 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
126 L2Norm< GridPart >
127 ::distance ( const UDiscreteFunctionType &u,
128 const VDiscreteFunctionType &v,
129 const PartitionSet& partitionSet ) const
130 {
131 typedef typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type RealType;
132 typedef FieldVector< RealType, 1 > ReturnType ;
133
134 // calculate integral over each element
135 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ );
136
137 // communicate_ indicates global norm
138 if( communicate_ )
139 {
140 sum[ 0 ] = comm().sum( sum[ 0 ] );
141 }
142
143 // return result, e.g. sqrt of calculated sum
144 return std::sqrt( sum[ 0 ] );
145 }
146
147 template< class GridPart >
148 template< class LocalFunctionType, class ReturnType >
149 inline void L2Norm< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
150 {
151 // evaluate norm locally
152 IntegratorType integrator( order );
153
154 FunctionSquare< LocalFunctionType > uLocal2( uLocal );
155
156 integrator.integrateAdd( entity, uLocal2, sum );
157 }
158
159 template< class GridPart >
160 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
161 inline void
162 L2Norm< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
163 {
164 // evaluate norm locally
165 IntegratorType integrator( order );
166
167 typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType;
168
169 LocalDistanceType dist( uLocal, vLocal );
170 FunctionSquare< LocalDistanceType > dist2( dist );
171
172 integrator.integrateAdd( entity, dist2, sum );
173 }
174
175
176 template< class GridPart >
177 template< class Function >
178 struct L2Norm< GridPart >::FunctionSquare
179 {
180 typedef Function FunctionType;
181
182 typedef typename FunctionType::RangeFieldType RangeFieldType;
183 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
184 typedef FieldVector< RealType, 1 > RangeType;
185
186 explicit FunctionSquare ( const FunctionType &function )
187 : function_( function )
188 {}
189
190 template< class Point >
191 void evaluate ( const Point &x, RangeType &ret ) const
192 {
193 typename FunctionType::RangeType phi;
194 function_.evaluate( x, phi );
195 ret = phi.two_norm2();
196 }
197
198 private:
199 const FunctionType &function_;
200 };
201
202
203 template< class GridPart >
204 template< class UFunction, class VFunction >
205 struct L2Norm< GridPart >::FunctionDistance
206 {
207 typedef UFunction UFunctionType;
208 typedef VFunction VFunctionType;
209
210 typedef typename UFunctionType::RangeFieldType RangeFieldType;
211 typedef typename UFunctionType::RangeType RangeType;
212 typedef typename UFunctionType::JacobianRangeType JacobianRangeType;
213
214 FunctionDistance ( const UFunctionType &u, const VFunctionType &v )
215 : u_( u ), v_( v )
216 {}
217
218 template< class Point >
219 void evaluate ( const Point &x, RangeType &ret ) const
220 {
221 RangeType phi;
222 u_.evaluate( x, ret );
223 v_.evaluate( x, phi );
224 ret -= phi;
225 }
226
227 template< class Point >
228 void jacobian ( const Point &x, JacobianRangeType &ret ) const
229 {
230 JacobianRangeType phi;
231 u_.jacobian( x, ret );
232 v_.jacobian( x, phi );
233 ret -= phi;
234 }
235
236 private:
237 const UFunctionType &u_;
238 const VFunctionType &v_;
239 };
240
241
242
243 // WeightedL2Norm
244 // --------------
245
246 template< class WeightFunction >
247 class WeightedL2Norm
248 : public L2Norm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType >
249 {
250 typedef WeightedL2Norm< WeightFunction > ThisType;
251 typedef L2Norm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType > BaseType;
252
253 public:
254 typedef WeightFunction WeightFunctionType;
255
256 typedef typename WeightFunctionType::DiscreteFunctionSpaceType WeightFunctionSpaceType;
257 typedef typename WeightFunctionSpaceType::GridPartType GridPartType;
258
259 protected:
260 template< class Function >
261 struct WeightedFunctionSquare;
262
263 typedef ConstLocalFunction< WeightFunctionType > LocalWeightFunctionType;
264 typedef typename WeightFunctionType::RangeType WeightType;
265
266 typedef typename BaseType::EntityType EntityType;
267 typedef typename BaseType::IntegratorType IntegratorType;
268
269 using BaseType::gridPart;
270 using BaseType::comm;
271
272 public:
273
274 using BaseType::norm;
275 using BaseType::distance;
276
277 explicit WeightedL2Norm ( const WeightFunctionType &weightFunction, const unsigned int order = 0, const bool communicate = true );
278
279 template< class LocalFunctionType, class ReturnType >
280 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
281
282 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
283 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
284
285 private:
286 LocalWeightFunctionType wfLocal_;
287 };
288
289
290
291
292 // Implementation of WeightedL2Norm
293 // --------------------------------
294
295 template< class WeightFunction >
296 inline WeightedL2Norm< WeightFunction >
297 ::WeightedL2Norm ( const WeightFunctionType &weightFunction, const unsigned int order, const bool communicate )
298 : BaseType( weightFunction.space().gridPart(), order, communicate ),
299 wfLocal_( weightFunction )
300 {
301 static_assert( (WeightFunctionSpaceType::dimRange == 1),
302 "Wight function must be scalar." );
303 }
304
305
306 template< class WeightFunction >
307 template< class LocalFunctionType, class ReturnType >
308 inline void
309 WeightedL2Norm< WeightFunction >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
310 {
311 // !!!! order !!!!
312 IntegratorType integrator( order );
313
314 wfLocal_.bind( entity );
315 WeightedFunctionSquare< LocalFunctionType > uLocal2( wfLocal_, uLocal );
316 integrator.integrateAdd( entity, uLocal2, sum );
317 wfLocal_.unbind();
318 }
319
320
321 template< class WeightFunction >
322 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
323 inline void
324 WeightedL2Norm< WeightFunction >::distanceLocal ( const EntityType& entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType& sum ) const
325 {
326 typedef typename BaseType::template FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType;
327
328 // !!!! order !!!!
329 IntegratorType integrator( order );
330
331 wfLocal_.bind( entity );
332 LocalDistanceType dist( uLocal, vLocal );
333 WeightedFunctionSquare< LocalDistanceType > dist2( wfLocal_, dist );
334
335 integrator.integrateAdd( entity, dist2, sum );
336 wfLocal_.unbind();
337 }
338
339
340 template< class WeightFunction >
341 template< class Function >
342 struct WeightedL2Norm< WeightFunction >::WeightedFunctionSquare
343 {
344 typedef Function FunctionType;
345
346 typedef typename FunctionType::RangeFieldType RangeFieldType;
347 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
348 typedef FieldVector< RealType, 1 > RangeType;
349
350 WeightedFunctionSquare ( const LocalWeightFunctionType &weightFunction,
351 const FunctionType &function )
352 : weightFunction_( weightFunction ),
353 function_( function )
354 {}
355
356 template< class Point >
357 void evaluate ( const Point &x, RangeType &ret ) const
358 {
359 WeightType weight;
360 weightFunction_.evaluate( x, weight );
361
362 typename FunctionType::RangeType phi;
363 function_.evaluate( x, phi );
364 ret = weight[ 0 ] * (phi * phi);
365 }
366
367 private:
368 const LocalWeightFunctionType &weightFunction_;
369 const FunctionType &function_;
370 };
371
372 } // end namespace Fem
373
374} // end namespace Dune
375
376#endif // #ifndef DUNE_FEM_L2NORM_HH
Dune::Fem::TupleDiscreteFunction
forward declaration
Definition: discretefunction.hh:51
Dune::PartitionSet::partitionSet
PartitionSet<(1<< p)> partitionSet()
Creates a PartitionSet for the given PartitionType.
Definition: partitionset.hh:221
Dune::Hybrid::forEach
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)