dune-fem 2.8.0
Loading...
Searching...
No Matches
domainintegral.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_INTEGRAL_HH
2#define DUNE_FEM_INTEGRAL_HH
3
4#include <numeric>
5#include <type_traits>
6
7#include <dune/grid/common/rangegenerators.hh>
8
11
13
14#include <dune/common/bartonnackmanifcheck.hh>
16
20
21namespace Dune
22{
23
24 namespace Fem
25 {
26 // IntegralBase
27 // ----------
28
29 template< class GridPart, class NormImplementation >
31 : public BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >,
32 NormImplementation >
33 {
35 NormImplementation > BaseType ;
37
38 public:
39 typedef GridPart GridPartType;
40
41 protected:
42 using BaseType :: asImp ;
43
44 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
45
46 template <bool uDiscrete, bool vDiscrete>
48 {
49 private:
50 template <class IteratorRange, class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType>
51 static ReturnType forEachLocal ( const ThisType &norm, const IteratorRange& iterators,
52 const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
53 const ReturnType &initialValue,
54 unsigned int order )
55 {
56 static_assert( uDiscrete && vDiscrete, "Distance can only be calculated between GridFunctions" );
57
58 ReturnType sum( 0 );
59 {
62 for( const EntityType &entity : iterators )
63 {
64 auto uGuard = bindGuard( uLocal, entity );
65 auto vGuard = bindGuard( vLocal, entity );
66
67 const unsigned int uOrder = uLocal.order();
68 const unsigned int vOrder = vLocal.order();
69 const unsigned int orderLocal = (order == 0 ? 2*std::max( uOrder, vOrder ) : order);
70 norm.distanceLocal( entity, orderLocal, uLocal, vLocal, sum );
71 }
72 }
73 return sum;
74 }
75
76 public:
77 template <class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet>
78 static ReturnType forEach ( const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
79 const ReturnType &initialValue,
80 const PartitionSet& partitionSet,
81 unsigned int order )
82 {
83 const int nThreads = MPIManager::numThreads();
84 if( nThreads == 1 )
85 return forEachLocal( norm, elements( norm.gridPart_, partitionSet ), u, v, initialValue, order );
86
87 std::vector< ReturnType > sums( nThreads, ReturnType(0) );
88
89 ThreadIterator< GridPartType, PartitionSet::partitionIterator() >
90 iterators( norm.gridPart_ );
91
92 auto doIntegrate = [ &norm, &iterators, &sums, &u, &v, &initialValue, &order ] ()
93 {
94 sums[ MPIManager::thread() ] = forEachLocal( norm, iterators, u, v, initialValue, order );
95 };
96
97 try {
98 // run threaded
99 MPIManager::run( doIntegrate );
100 }
101 catch ( const SingleThreadModeError& e )
102 {
103 // return single threaded variant in this case
104 return forEachLocal( norm, elements( norm.gridPart_, partitionSet ), u, v, initialValue, order );
105 }
106
107 return std::accumulate( sums.begin(), sums.end(), ReturnType(0) );
108 }
109 };
110
111 // this specialization creates a grid function adapter
112 template <bool vDiscrete>
113 struct ForEachCaller<false, vDiscrete>
114 {
115 template <class F,
116 class VDiscreteFunctionType,
117 class ReturnType,
118 class PartitionSet>
119 static
120 ReturnType forEach ( const ThisType& norm,
121 const F& f, const VDiscreteFunctionType& v,
122 const ReturnType& initialValue,
123 const PartitionSet& partitionSet,
124 const unsigned int order )
125 {
126 typedef GridFunctionAdapter< F, GridPartType> GridFunction ;
127 GridFunction u( "Integral::adapter" , f , v.space().gridPart(), v.space().order() );
128
130 forEach( norm, u, v, initialValue, partitionSet, order );
131 }
132 };
133
134 // this specialization simply switches arguments
135 template <bool uDiscrete>
136 struct ForEachCaller<uDiscrete, false>
137 {
138 template <class UDiscreteFunctionType,
139 class F,
140 class ReturnType,
141 class PartitionSet>
142 static
143 ReturnType forEach ( const ThisType& norm,
144 const UDiscreteFunctionType& u,
145 const F& f,
146 const ReturnType& initialValue,
147 const PartitionSet& partitionSet,
148 const unsigned int order )
149 {
151 forEach( norm, f, u, initialValue, partitionSet, order );
152 }
153 };
154
155 template <class IteratorRange, class UDiscreteFunctionType, class ReturnType>
156 ReturnType forEachLocal ( const IteratorRange& iterators,
157 const UDiscreteFunctionType &u,
158 const ReturnType &initialValue,
159 unsigned int order ) const
160 {
161 ReturnType sum( 0 );
162 {
164 for( const EntityType &entity : iterators )
165 {
166 auto uGuard = bindGuard( uLocal, entity );
167
168 const unsigned int uOrder = uLocal.order();
169 const unsigned int orderLocal = (order == 0 ? 2*uOrder : order);
170 normLocal( entity, orderLocal, uLocal, sum );
171 }
172 }
173 return sum;
174 }
175
176 template< class DiscreteFunctionType, class ReturnType, class PartitionSet >
177 ReturnType forEach ( const DiscreteFunctionType &u, const ReturnType &initialValue,
178 const PartitionSet& partitionSet,
179 unsigned int order = 0 ) const
180 {
181 static_assert( (std::is_base_of<Fem::HasLocalFunction, DiscreteFunctionType>::value),
182 "Norm only implemented for quantities implementing a local function!" );
183 const int nThreads = MPIManager::numThreads();
184 if( nThreads == 1 )
185 return forEachLocal( elements( gridPart_, partitionSet ), u, initialValue, order );
186
187 std::vector< ReturnType > sums( nThreads, ReturnType(0) );
188
189 ThreadIterator< GridPartType, PartitionSet::partitionIterator() >
190 iterators( gridPart_ );
191
192 auto doIntegrate = [ this, &iterators, &sums, &u, &initialValue, &order ] ()
193 {
194 sums[ MPIManager::thread() ] = forEachLocal( iterators, u, initialValue, order );
195 };
196
197 try {
198 // run threaded
199 MPIManager::run( doIntegrate );
200 }
201 catch ( const SingleThreadModeError& e )
202 {
203 // return single threaded variant in this case
204 return forEachLocal( elements( gridPart_, partitionSet ), u, initialValue, order );
205 }
206
207 return std::accumulate( sums.begin(), sums.end(), ReturnType(0) );
208 }
209
210 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet >
211 ReturnType forEach ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v,
212 const ReturnType &initialValue, const PartitionSet& partitionSet,
213 unsigned int order = 0 ) const
214 {
215 enum { uDiscrete = std::is_convertible<UDiscreteFunctionType, HasLocalFunction>::value };
216 enum { vDiscrete = std::is_convertible<VDiscreteFunctionType, HasLocalFunction>::value };
217
218 // call forEach depending on which argument is a grid function,
219 // i.e. has a local function
221 forEach( *this, u, v, initialValue, partitionSet, order );
222 }
223
224 public:
226 : gridPart_( gridPart )
227 {}
228
229 protected:
230 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
231 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
232 {
233 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().distanceLocal( entity, order, uLocal, vLocal, sum ) );
234 }
235
236 template< class LocalFunctionType, class ReturnType >
237 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
238 {
239 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().normLocal( entity, order, uLocal, sum ) );
240 }
241
242 const GridPartType &gridPart () const { return gridPart_; }
243
244 const typename GridPartType::CollectiveCommunicationType& comm () const
245 {
246 return gridPart().comm();
247 }
248
249 bool checkCommunicateFlag( bool communicate ) const
250 {
251#ifndef NDEBUG
252 bool commMax = communicate;
253 assert( communicate == comm().max( commMax ) );
254#endif
255 return communicate;
256 }
257
258 private:
259 const GridPartType &gridPart_;
260
261 };
262
263 // Integral
264 // --------
265
266 template< class GridPart >
267 class Integral : public IntegralBase< GridPart, Integral< GridPart > >
268 {
271
272 public:
273 typedef GridPart GridPartType;
274
275 using BaseType :: gridPart ;
276 using BaseType :: comm ;
277
278 protected:
279
280 template< class UFunction, class VFunction >
281 struct FunctionDistance;
282
285
286 const unsigned int order_;
287 const bool communicate_;
288 public:
294 explicit Integral ( const GridPartType &gridPart,
295 const unsigned int order = 0,
296 const bool communicate = true );
297
298
300 template< class DiscreteFunctionType, class PartitionSet >
301 typename DiscreteFunctionType::RangeType
302 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const;
303
305 template< class DiscreteFunctionType >
306 typename DiscreteFunctionType::RangeType
307 norm ( const DiscreteFunctionType &u ) const
308 {
309 return norm( u, Partitions::interior );
310 }
311
313 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
314 typename UDiscreteFunctionType::RangeType
315 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const;
316
318 template< class UDiscreteFunctionType, class VDiscreteFunctionType >
319 typename UDiscreteFunctionType::RangeType
320 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const
321 {
322 return distance( u, v, Partitions::interior );
323 }
324
325 template< class LocalFunctionType, class ReturnType >
326 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const;
327
328 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
329 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const;
330 };
331
332
333
334 // Implementation of Integral
335 // ------------------------
336
337 template< class GridPart >
338 inline Integral< GridPart >::Integral ( const GridPartType &gridPart, const unsigned int order, const bool communicate )
339 : BaseType( gridPart ),
340 order_( order ),
341 communicate_( BaseType::checkCommunicateFlag( communicate ) )
342 {}
343
344
345 template< class GridPart >
346 template< class DiscreteFunctionType, class PartitionSet >
347 inline typename DiscreteFunctionType::RangeType
348 Integral< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const
349 {
350 typedef typename DiscreteFunctionType::RangeType RangeType;
351 typedef RangeType ReturnType ;
352
353 // calculate integral over each element
354 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ );
355
356 // communicate_ indicates global norm
357 if( communicate_ )
358 {
359 sum = comm().sum( sum );
360 }
361
362 return sum;
363 }
364
365
366 template< class GridPart >
367 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet >
368 inline typename UDiscreteFunctionType::RangeType
370 ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const
371 {
372 typedef typename UDiscreteFunctionType::RangeType RangeType;
373 typedef RangeType ReturnType ;
374
375 // calculate integral over each element
376 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ );
377
378 // communicate_ indicates global norm
379 if( communicate_ )
380 {
381 sum = comm().sum( sum );
382 }
383
384 return sum;
385 }
386
387 template< class GridPart >
388 template< class LocalFunctionType, class ReturnType >
389 inline void
390 Integral< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const
391 {
392 Integrator< QuadratureType > integrator( order );
393
394 integrator.integrateAdd( entity, uLocal, sum );
395 }
396
397 template< class GridPart >
398 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType >
399 inline void
400 Integral< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const
401 {
402 Integrator< QuadratureType > integrator( order );
403
405
406 LocalDistanceType dist( uLocal, vLocal );
407
408 integrator.integrateAdd( entity, dist, sum );
409 }
410
411 template< class GridPart >
412 template< class UFunction, class VFunction >
413 struct Integral< GridPart >::FunctionDistance
414 {
415 typedef UFunction UFunctionType;
416 typedef VFunction VFunctionType;
417
418 typedef typename UFunctionType::RangeFieldType RangeFieldType;
419 typedef typename UFunctionType::RangeType RangeType;
420 typedef typename UFunctionType::JacobianRangeType JacobianRangeType;
421
423 : u_( u ), v_( v )
424 {}
425
426 template< class Point >
427 void evaluate ( const Point &x, RangeType &ret ) const
428 {
429 RangeType phi;
430 u_.evaluate( x, ret );
431 v_.evaluate( x, phi );
432 ret -= phi;
433 }
434
435 template< class Point >
436 void jacobian ( const Point &x, JacobianRangeType &ret ) const
437 {
439 u_.jacobian( x, ret );
440 v_.jacobian( x, phi );
441 ret -= phi;
442 }
443
444 private:
445 const UFunctionType &u_;
446 const VFunctionType &v_;
447 };
448
449 } // namespace Fem
450
451} // namespace Dune
452
453#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:604
static double max(const Double &v, const double p)
Definition: double.hh:398
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
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:33
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: domainintegral.hh:44
ReturnType forEachLocal(const IteratorRange &iterators, const UDiscreteFunctionType &u, const ReturnType &initialValue, unsigned int order) const
Definition: domainintegral.hh:156
GridPart GridPartType
Definition: domainintegral.hh:39
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:231
IntegralBase(const GridPartType &gridPart)
Definition: domainintegral.hh:225
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:237
const GridPartType::CollectiveCommunicationType & comm() const
Definition: domainintegral.hh:244
bool checkCommunicateFlag(bool communicate) const
Definition: domainintegral.hh:249
ReturnType forEach(const DiscreteFunctionType &u, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:177
ReturnType forEach(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order=0) const
Definition: domainintegral.hh:211
const GridPartType & gridPart() const
Definition: domainintegral.hh:242
Definition: domainintegral.hh:48
static ReturnType forEach(const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const ReturnType &initialValue, const PartitionSet &partitionSet, unsigned int order)
Definition: domainintegral.hh:78
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:120
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:143
Definition: domainintegral.hh:268
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u) const
|| u ||_L1 on interior partition entities
Definition: domainintegral.hh:307
CachingQuadrature< GridPartType, 0 > QuadratureType
Definition: domainintegral.hh:284
GridPart GridPartType
Definition: domainintegral.hh:273
const bool communicate_
Definition: domainintegral.hh:287
UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType &u, const VDiscreteFunctionType &v) const
|| u - v ||_L2 on interior partition entities
Definition: domainintegral.hh:320
BaseType::EntityType EntityType
Definition: domainintegral.hh:283
DiscreteFunctionType::RangeType norm(const DiscreteFunctionType &u, const PartitionSet &partitionSet) const
|| u ||_L1 on given set of entities (partition set)
Definition: domainintegral.hh:348
Integral(const GridPartType &gridPart, const unsigned int order=0, const bool communicate=true)
constructor
Definition: domainintegral.hh:338
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:370
void normLocal(const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum) const
Definition: domainintegral.hh:390
const unsigned int order_
Definition: domainintegral.hh:286
void distanceLocal(const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum) const
Definition: domainintegral.hh:400
Definition: domainintegral.hh:414
UFunctionType::RangeType RangeType
Definition: domainintegral.hh:419
void evaluate(const Point &x, RangeType &ret) const
Definition: domainintegral.hh:427
UFunction UFunctionType
Definition: domainintegral.hh:415
FunctionDistance(const UFunctionType &u, const VFunctionType &v)
Definition: domainintegral.hh:422
UFunctionType::RangeFieldType RangeFieldType
Definition: domainintegral.hh:418
VFunction VFunctionType
Definition: domainintegral.hh:416
void jacobian(const Point &x, JacobianRangeType &ret) const
Definition: domainintegral.hh:436
UFunctionType::JacobianRangeType JacobianRangeType
Definition: domainintegral.hh:420
BaseType::EntityType EntityType
Definition: lpnorm.hh:143
Exception thrown when a code segment that is supposed to be only accessed in single thread mode is ac...
Definition: mpimanager.hh:43
static void run(F &&f, Args &&... args)
run functor f with given arguments args in threaded mode
Definition: mpimanager.hh:440
static int thread()
return thread number
Definition: mpimanager.hh:424
static int numThreads()
return number of current threads
Definition: mpimanager.hh:421
Thread iterators.
Definition: threaditerator.hh:22
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