DUNE-FEM (unstable)

domainintegral.hh
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
9#include <dune/fem/quadrature/cachingquadrature.hh>
10#include <dune/fem/quadrature/integrator.hh>
11
12#include <dune/fem/function/common/gridfunctionadapter.hh>
13
14#include <dune/common/bartonnackmanifcheck.hh>
15#include <dune/fem/misc/bartonnackmaninterface.hh>
16
17#include <dune/fem/misc/mpimanager.hh>
18#include <dune/fem/misc/threads/threaditerator.hh>
19#include <dune/fem/common/bindguard.hh>
20
21namespace Dune
22{
23
24 namespace Fem
25 {
26 // IntegralBase
27 // ----------
28
29 template< class GridPart, class NormImplementation >
30 class IntegralBase
31 : public BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >,
32 NormImplementation >
33 {
34 typedef BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >,
35 NormImplementation > BaseType ;
36 typedef IntegralBase< GridPart, NormImplementation > ThisType;
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>
47 struct ForEachCaller
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 {
60 ConstLocalFunction< UDiscreteFunctionType > uLocal( u );
61 ConstLocalFunction< VDiscreteFunctionType > vLocal( v );
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
129 return ForEachCaller< true, vDiscrete > ::
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 {
150 return ForEachCaller< false, uDiscrete > ::
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 {
163 ConstLocalFunction< UDiscreteFunctionType > uLocal( u );
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
220 return ForEachCaller< uDiscrete, vDiscrete > ::
221 forEach( *this, u, v, initialValue, partitionSet, order );
222 }
223
224 public:
225 explicit IntegralBase ( const GridPartType &gridPart )
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::CommunicationType& 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 {
269 typedef IntegralBase< GridPart, Integral< GridPart > > BaseType ;
270 typedef Integral< GridPart > ThisType;
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
283 typedef typename BaseType::EntityType EntityType;
284 typedef CachingQuadrature< GridPartType, 0 > QuadratureType;
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
369 Integral< GridPart >
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
404 typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType;
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
422 FunctionDistance ( const UFunctionType &u, const VFunctionType &v )
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 {
438 JacobianRangeType phi;
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
bartonnackmanifcheck.hh
Provides check for implementation of interface methods when using static polymorphism,...
CHECK_AND_CALL_INTERFACE_IMPLEMENTATION
#define CHECK_AND_CALL_INTERFACE_IMPLEMENTATION(__interface_method_to_call__)
Definition: bartonnackmanifcheck.hh:61
Dune::Fem::DiscreteFunctionDefault< TupleDiscreteFunction< DiscreteFunctions ... > >::RangeType
DiscreteFunctionSpaceType::RangeType RangeType
type of range vector
Definition: discretefunction.hh:614
Dune::Fem::GridPartInterface< AdaptiveLeafGridPartTraits< Grid, idxpitype, onlyCodimensionZero > >::CommunicationType
Traits::CommunicationType CommunicationType
Collective communication.
Definition: gridpart.hh:97
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::Hybrid::max
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::PartitionSet::partitionIterator
static constexpr PartitionIteratorType partitionIterator()
Returns the PartitionIteratorType that can be used to iterate over the partitions in the set.
Definition: partitionset.hh:182
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)