dune-fem 2.8-git
domainintegral.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_INTEGRAL_HH
2#define DUNE_FEM_INTEGRAL_HH
3
4#include <type_traits>
5
6#include <dune/grid/common/rangegenerators.hh>
7
10
12
13#include <dune/common/bartonnackmanifcheck.hh>
15
16namespace Dune
17{
18
19 namespace Fem
20 {
21 // IntegralBase
22 // ----------
23
24 template< class GridPart, class NormImplementation >
26 : public BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >,
27 NormImplementation >
28 {
30 NormImplementation > BaseType ;
32
33 public:
34 typedef GridPart GridPartType;
35
36 protected:
37 using BaseType :: asImp ;
38
39 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
40
41 template <bool uDiscrete, bool vDiscrete>
43 {
44 template <class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet>
45 static ReturnType forEach ( const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
46 const ReturnType &initialValue,
47 const PartitionSet& partitionSet,
48 unsigned int order )
49 {
50 static_assert( uDiscrete && vDiscrete, "Distance can only be calculated between GridFunctions" );
51
52 ReturnType sum( 0 );
55 for( const EntityType &entity : elements( norm.gridPart_, partitionSet ) )
56 {
57 uLocal.bind( entity );
58 vLocal.bind( entity );
59 const unsigned int uOrder = uLocal.order();
60 const unsigned int vOrder = vLocal.order();
61 const unsigned int orderLocal = (order == 0 ? 2*std::max( uOrder, vOrder ) : order);
62 norm.distanceLocal( entity, orderLocal, uLocal, vLocal, sum );
63 uLocal.unbind();
64 vLocal.unbind();
65 }
66 return sum;
67 }
68 };
69
70 // this specialization creates a grid function adapter
71 template <bool vDiscrete>
72 struct ForEachCaller<false, vDiscrete>
73 {
74 template <class F,
75 class VDiscreteFunctionType,
76 class ReturnType,
77 class PartitionSet>
78 static
79 ReturnType forEach ( const ThisType& norm,
80 const F& f, const VDiscreteFunctionType& v,
81 const ReturnType& initialValue,
82 const PartitionSet& partitionSet,
83 const unsigned int order )
84 {
85 typedef GridFunctionAdapter< F, GridPartType> GridFunction ;
86 GridFunction u( "Integral::adapter" , f , v.space().gridPart(), v.space().order() );
87
89 forEach( norm, u, v, initialValue, partitionSet, order );
90 }
91 };
92
93 // this specialization simply switches arguments
94 template <bool uDiscrete>
95 struct ForEachCaller<uDiscrete, false>
96 {
97 template <class UDiscreteFunctionType,
98 class F,
99 class ReturnType,
100 class PartitionSet>
101 static
102 ReturnType forEach ( const ThisType& norm,
103 const UDiscreteFunctionType& u,
104 const F& f,
105 const ReturnType& initialValue,
106 const PartitionSet& partitionSet,
107 const unsigned int order )
108 {
110 forEach( norm, f, u, initialValue, partitionSet, order );
111 }
112 };
113
114 template< class DiscreteFunctionType, class ReturnType, class PartitionSet >
115 ReturnType forEach ( const DiscreteFunctionType &u, const ReturnType &initialValue,
116 const PartitionSet& partitionSet,
117 unsigned int order = 0 ) const
118 {
119 static_assert( (std::is_base_of<Fem::HasLocalFunction, DiscreteFunctionType>::value),
120 "Norm only implemented for quantities implementing a local function!" );
121
122 ReturnType sum( 0 );
124 for( const EntityType &entity : elements( gridPart_, partitionSet ) )
125 {
126 uLocal.bind( entity );
127 const unsigned int orderLocal = (order == 0 ? 2*uLocal.order() : order);
128 normLocal( entity, orderLocal, uLocal, sum );
129 uLocal.unbind();
130 }
131 return sum;
132 }
133
134 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet >
135 ReturnType forEach ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
136 const ReturnType &initialValue, const PartitionSet& partitionSet,
137 unsigned int order = 0 ) const
138 {
139 enum { uDiscrete = std::is_convertible<UDiscreteFunctionType, HasLocalFunction>::value };
140 enum { vDiscrete = std::is_convertible<VDiscreteFunctionType, HasLocalFunction>::value };
141
142 // call forEach depending on which argument is a grid function,
143 // i.e. has a local function
145 forEach( *this, u, v, initialValue, partitionSet, order );
146 }
147
148 public:
150 : gridPart_( gridPart )
151 {}
152
153 protected:
154 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
155 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
156 {
157 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().distanceLocal( entity, order, uLocal, vLocal, sum ) );
158 }
159
160 template< class LocalFunctionType, class ReturnType >
161 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
162 {
163 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().normLocal( entity, order, uLocal, sum ) );
164 }
165
166 const GridPartType &gridPart () const { return gridPart_; }
167
168 const typename GridPartType::CollectiveCommunicationType& comm () const
169 {
170 return gridPart().comm();
171 }
172
173 bool checkCommunicateFlag( bool communicate ) const
174 {
175#ifndef NDEBUG
176 bool commMax = communicate;
177 assert( communicate == comm().max( commMax ) );
178#endif
179 return communicate;
180 }
181
182 private:
183 const GridPartType &gridPart_;
184 };
185
186 // Integral
187 // ------
188
189 template< class GridPart >
190 class Integral : public IntegralBase< GridPart, Integral< GridPart > >
191 {
194
195 public:
196 typedef GridPart GridPartType;
197
199 using BaseType :: comm ;
200
201 protected:
202
203 template< class UFunction, class VFunction >
204 struct FunctionDistance;
205
208
209 const unsigned int order_;
210 const bool communicate_;
211 public:
217 explicit Integral ( const GridPartType &gridPart,
218 const unsigned int order = 0,
219 const bool communicate = true );
220
221
223 template< class DiscreteFunctionType, class PartitionSet >
224 typename DiscreteFunctionType::RangeType
225 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const;
226
228 template< class DiscreteFunctionType >
229 typename DiscreteFunctionType::RangeType
230 norm ( const DiscreteFunctionType &u ) const
231 {
232 return norm( u, Partitions::interior );
233 }
234
236 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
237 typename UDiscreteFunctionType::RangeType
238 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const;
239
241 template< class UDiscreteFunctionType, class VDiscreteFunctionType >
242 typename UDiscreteFunctionType::RangeType
243 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const
244 {
245 return distance( u, v, Partitions::interior );
246 }
247
248 template< class LocalFunctionType, class ReturnType >
249 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
250
251 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
252 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
253 };
254
255
256
257 // Implementation of Integral
258 // ------------------------
259
260 template< class GridPart >
261 inline Integral< GridPart >::Integral ( const GridPartType &gridPart, const unsigned int order, const bool communicate )
262 : BaseType( gridPart ),
263 order_( order ),
264 communicate_( BaseType::checkCommunicateFlag( communicate ) )
265 {}
266
267
268 template< class GridPart >
269 template< class DiscreteFunctionType, class PartitionSet >
270 inline typename DiscreteFunctionType::RangeType
271 Integral< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const
272 {
273 typedef typename DiscreteFunctionType::RangeType RangeType;
274 typedef RangeType ReturnType ;
275
276 // calculate integral over each element
277 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ );
278
279 // communicate_ indicates global norm
280 if( communicate_ )
281 {
282 sum = comm().sum( sum );
283 }
284
285 return sum;
286 }
287
288
289 template< class GridPart >
290 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
291 inline typename UDiscreteFunctionType::RangeType
293 ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const
294 {
295 typedef typename UDiscreteFunctionType::RangeType RangeType;
296 typedef RangeType ReturnType ;
297
298 // calculate integral over each element
299 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ );
300
301 // communicate_ indicates global norm
302 if( communicate_ )
303 {
304 sum = comm().sum( sum );
305 }
306
307 return sum;
308 }
309
310 template< class GridPart >
311 template< class LocalFunctionType, class ReturnType >
312 inline void
313 Integral< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
314 {
315 Integrator< QuadratureType > integrator( order );
316
317 integrator.integrateAdd( entity, uLocal, sum );
318 }
319
320 template< class GridPart >
321 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
322 inline void
323 Integral< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
324 {
325 Integrator< QuadratureType > integrator( order );
326
328
329 LocalDistanceType dist( uLocal, vLocal );
330
331 integrator.integrateAdd( entity, dist, sum );
332 }
333
334 template< class GridPart >
335 template< class UFunction, class VFunction >
336 struct Integral< GridPart >::FunctionDistance
337 {
338 typedef UFunction UFunctionType;
339 typedef VFunction VFunctionType;
340
341 typedef typename UFunctionType::RangeFieldType RangeFieldType;
342 typedef typename UFunctionType::RangeType RangeType;
343 typedef typename UFunctionType::JacobianRangeType JacobianRangeType;
344
346 : u_( u ), v_( v )
347 {}
348
349 template< class Point >
350 void evaluate ( const Point &x, RangeType &ret ) const
351 {
352 RangeType phi;
353 u_.evaluate( x, ret );
354 v_.evaluate( x, phi );
355 ret -= phi;
356 }
357
358 template< class Point >
359 void jacobian ( const Point &x, JacobianRangeType &ret ) const
360 {
362 u_.jacobian( x, ret );
363 v_.jacobian( x, phi );
364 ret -= phi;
365 }
366
367 private:
368 const UFunctionType &u_;
369 const VFunctionType &v_;
370 };
371
372 } // namespace Fem
373
374} // namespace Dune
375
376#endif // #ifndef DUNE_FEM_INTEGRAL_HH
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
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 std::decay_t< T > sum(T a)
Definition: utility.hh:33
BasicGridFunctionAdapter provides local functions for a Function.
Definition: gridfunctionadapter.hh:79
Definition: bartonnackmaninterface.hh:17
const NormImplementation & asImp() const
Definition: bartonnackmaninterface.hh:37
Definition: domainintegral.hh:28
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: domainintegral.hh:39
GridPart GridPartType
Definition: domainintegral.hh:34
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:155
IntegralBase(const GridPartType &gridPart)
Definition: domainintegral.hh:149
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:161
const GridPartType::CollectiveCommunicationType & comm() const
Definition: domainintegral.hh:168
bool checkCommunicateFlag(bool communicate) const
Definition: domainintegral.hh:173
ReturnType forEach(const DiscreteFunctionType &u, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:115
ReturnType forEach(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:135
const GridPartType & gridPart() const
Definition: domainintegral.hh:166
Definition: domainintegral.hh:43
static ReturnType forEach(const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order)
Definition: domainintegral.hh:45
static ReturnType forEach(const ThisType &norm, const F &f, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, const unsigned int order)
Definition: domainintegral.hh:79
static ReturnType forEach(const ThisType &norm, const UDiscreteFunctionType &u, const F &f, const ReturnType &initialValue, const PartitionSet &partitionSet, const unsigned int order)
Definition: domainintegral.hh:102
Definition: domainintegral.hh:191
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u) const
|| u ||_L1 on interior partition entities
Definition: domainintegral.hh:230
CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: domainintegral.hh:207
GridPart GridPartType
Definition: domainintegral.hh:196
const bool communicate_
Definition: domainintegral.hh:210
UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v) const
|| u - v ||_L2 on interior partition entities
Definition: domainintegral.hh:243
BaseType::EntityType EntityType
Definition: domainintegral.hh:204
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u, const PartitionSet &partitionSet) const
|| u ||_L1 on given set of entities (partition set)
Definition: domainintegral.hh:271
Integral(const GridPartType &gridPart, const unsigned int order=0, const bool communicate=true)
constructor
Definition: domainintegral.hh:261
UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet &partitionSet) const
|| u - v ||_L2 on given set of entities (partition set)
Definition: domainintegral.hh:293
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:313
const unsigned int order_
Definition: domainintegral.hh:209
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:323
Definition: domainintegral.hh:337
UFunctionType::RangeType RangeType
Definition: domainintegral.hh:342
void evaluate(const Point &x, RangeType &ret) const
Definition: domainintegral.hh:350
UFunction UFunctionType
Definition: domainintegral.hh:338
FunctionDistance(const UFunctionType &u, const VFunctionType &v)
Definition: domainintegral.hh:345
UFunctionType::RangeFieldType RangeFieldType
Definition: domainintegral.hh:341
VFunction VFunctionType
Definition: domainintegral.hh:339
void jacobian(const Point &x, JacobianRangeType &ret) const
Definition: domainintegral.hh:359
UFunctionType::JacobianRangeType JacobianRangeType
Definition: domainintegral.hh:343
BaseType::EntityType EntityType
Definition: lpnorm.hh:143
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