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 
9 namespace 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:
44       /** \brief constructor
45        *    \param gridPart     specific gridPart for selection of entities
46        *    \param order        order of integration quadrature (default = 2*space.order())
47        *    \param communicate  if true global (over all ranks) norm is computed (default = true)
48        */
49       explicit L2Norm ( const GridPartType &gridPart,
50                         const unsigned int order = 0,
51                         const bool communicate = true );
52 
53       //! || u ||_L2 on given set of entities (partition set)
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 
58       //! || u ||_L2 on interior partition entities
59       template< class DiscreteFunctionType >
60       typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type
norm(const DiscreteFunctionType & u) const61       norm ( const DiscreteFunctionType &u ) const
62       {
63         return norm( u, Partitions::interior );
64       }
65 
66       //! || u - v ||_L2 on given set of entities (partition set)
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 
73       //! || u - v ||_L2 on interior partition entities
74       template< class UDiscreteFunctionType, class VDiscreteFunctionType >
75       typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type
distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v) const76       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 >
L2Norm(const GridPartType & gridPart,const unsigned int order,const bool communicate)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
norm(const DiscreteFunctionType & u,const PartitionSet & partitionSet) const104     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 >
distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v,const PartitionSet & partitionSet) const127       ::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 >
normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const149     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
distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const162     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 
FunctionSquareDune::Fem::L2Norm::FunctionSquare186       explicit FunctionSquare ( const FunctionType &function )
187       : function_( function )
188       {}
189 
190       template< class Point >
evaluateDune::Fem::L2Norm::FunctionSquare191       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 
FunctionDistanceDune::Fem::L2Norm::FunctionDistance214       FunctionDistance ( const UFunctionType &u, const VFunctionType &v )
215       : u_( u ), v_( v )
216       {}
217 
218       template< class Point >
evaluateDune::Fem::L2Norm::FunctionDistance219       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 >
jacobianDune::Fem::L2Norm::FunctionDistance228       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 >
WeightedL2Norm(const WeightFunctionType & weightFunction,const unsigned int order,const bool communicate)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
normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const309     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
distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const324     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 
WeightedFunctionSquareDune::Fem::WeightedL2Norm::WeightedFunctionSquare350       WeightedFunctionSquare ( const LocalWeightFunctionType &weightFunction,
351                                const FunctionType &function )
352       : weightFunction_( weightFunction ),
353         function_( function )
354       {}
355 
356       template< class Point >
evaluateDune::Fem::WeightedL2Norm::WeightedFunctionSquare357       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
377