dune-fem  2.8-git
quadprovider.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_QUADPROVIDER_HH
2 #define DUNE_FEM_QUADPROVIDER_HH
3 
4 #include <iostream>
5 #include <memory>
6 #include <map>
7 #include <vector>
8 
12 
13 namespace Dune
14 {
15 
16  namespace Fem
17  {
18 
25  {
26  int id_;
27 
28  // something like maxOrder
29  static const int maxFirst = 25 ;
30  // something like different quadratures of the same order
31  static const int maxSecond = 10 ;
32  public:
34  FemQuadratureKey() : id_( -1 ) {}
35 
37  FemQuadratureKey( const FemQuadratureKey& key ) = default;
38 
39  // this may need to be adjusted in case more than 100 quadratures
40  // need to be stored
41  static const int highest_order = maxFirst * maxSecond ;
42 
44  FemQuadratureKey( const int first, const int second )
45  : id_( second * maxFirst + first )
46  {
47  assert( first < maxFirst );
48  assert( second < maxSecond );
49  }
50 
52  FemQuadratureKey( const int first )
53  : id_( first )
54  {
55  assert( first < maxFirst );
56  }
57 
59  operator int () const { return id_; }
60 
62  int first () const { return id_ % highest_order ; }
64  int second () const { return id_ / highest_order ; }
65  };
66 
67 
77  template< unsigned int dummy >
79  {
80  private:
82  template< class QuadImp, class QuadratureKey >
83  class QuadratureStorage
84  {
85  public:
86  typedef QuadImp QuadType;
87 
88  protected:
89  typedef std :: map< QuadratureKey, std::unique_ptr< QuadType > > StorageType;
90  StorageType storage_;
91 
92  public:
93  QuadratureStorage () {}
94 
95  QuadImp &getQuadrature( const GeometryType &geometry, const QuadratureKey& key )
96  {
97  QuadType* quadPtr = nullptr;
98  auto it = storage_.find( key );
99  if( it == storage_.end() )
100  {
101  // make sure we work in single thread mode when quadrature is created
102  assert( Fem :: ThreadManager:: singleThreadMode() );
103  quadPtr = new QuadImp( geometry, key, IdProvider :: instance().newId() );
104  storage_[ key ].reset( quadPtr );
105  }
106  else
107  {
108  quadPtr = it->second.operator ->();
109  }
110 
111  assert( quadPtr != nullptr );
112  return *quadPtr;
113  }
114  }; // end class QuadratureStorage
115 
117  template< class QuadImp>
118  class QuadratureStorage< QuadImp, int > // QuadratureKey == int
119  {
120  public:
121  typedef QuadImp QuadType;
122 
123  protected:
124  std::vector< std::unique_ptr< QuadType > > storage_;
125 
126  public:
127  QuadratureStorage ()
128  : storage_( QuadType :: maxOrder() + 1 )
129  {
130  }
131 
132  QuadImp &getQuadrature( const GeometryType &geometry, unsigned int order )
133  {
134  if(order >= storage_.size() )
135  {
136 #ifndef NDEBUG
137  static bool showMessage = true ;
138  if( showMessage )
139  {
140  std::cerr << "WARNING: QuadratureStorage::getQuadrature: A quadrature of order " << order
141  << " is not implemented!" << std::endl
142  << "Choosing maximum order: " << storage_.size()-1 << std::endl << std::endl;
143  showMessage = false;
144  }
145 #endif
146  order = storage_.size() - 1;
147  }
148 
149  auto& quadPtr = storage_[ order ];
150  if( ! quadPtr )
151  {
152  // make sure we work in single thread mode when quadrature is created
153  assert( Fem :: ThreadManager:: singleThreadMode() );
154  quadPtr.reset( new QuadImp( geometry, int(order), IdProvider :: instance().newId() ) );
155  }
156 
157  assert( quadPtr );
158  return *quadPtr;
159  }
160  }; // end class QuadratureStorage
161 
163  template< class QuadImp >
164  class QuadratureStorage< QuadImp, FemQuadratureKey >
165  : public QuadratureStorage< QuadImp, int >
166  {
167  };
168 
169  public:
175  template< class QuadImp, class QuadratureKey >
176  static const QuadImp &provideQuad( const GeometryType& geometry,
177  const QuadratureKey& key )
178  {
179  static QuadratureStorage< QuadImp, QuadratureKey > storage;
180  return storage.getQuadrature( geometry, key );
181  }
182 
189  template< class QuadImp, class QuadratureKey >
190  static const QuadImp &provideQuad( const GeometryType& geometry,
191  const QuadratureKey& key,
192  const int defaultOrder )
193  {
194  // this function should only be called for geometry types equal to none
195  assert( geometry.isNone() );
196  DUNE_THROW(NotImplemented,"provideQuad for polyhedral cells (defaultOrder = 0) not implemented for arbitrary QuadratureKey!");
197  QuadImp* ptr = nullptr;
198  return *ptr;
199  }
200 
207  template< class QuadImp >
208  static const QuadImp &provideQuad( const GeometryType& geometry,
209  const int ,
210  const int defaultOrder )
211  {
212  assert( geometry.isNone() );
213  static QuadratureStorage< QuadImp, int > storage;
214  return storage.getQuadrature( geometry, defaultOrder );
215  }
216  };
217 
233  template< typename FieldImp, int dim, template< class, int > class QuadratureTraits >
235  {
236  public:
237  typedef FieldImp FieldType;
238 
239  enum { dimension = dim };
240 
241  private:
243 
244  typedef QuadratureTraits< FieldType, dimension > QuadratureTraitsType;
245 
246  public:
248  typedef typename QuadratureTraitsType :: CubeQuadratureType CubeQuadratureType;
249 
251  typedef typename QuadratureTraitsType :: IntegrationPointListType
253 
255  typedef typename QuadratureTraitsType :: QuadratureKeyType QuadratureKeyType;
256 
258  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
259  const QuadratureKeyType& quadKey )
260  {
261  assert( geometry.isCube() );
262  return QuadCreator< 0 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey );
263  }
265  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
266  const GeometryType &elementGeometry,
267  const QuadratureKeyType& quadKey )
268  {
269  return getQuadrature( geometry, quadKey );
270  }
271 
272  QuadratureProvider() = delete;
273  QuadratureProvider( const ThisType& ) = delete;
274  QuadratureProvider &operator=( const ThisType& ) = delete;
275  };
276 
277 
278 
280  template< typename FieldImp, template< class, int > class QuadratureTraits >
281  class QuadratureProvider< FieldImp, 0, QuadratureTraits >
282  {
283  public:
284  typedef FieldImp FieldType;
285 
286  enum { dimension = 0 };
287 
288  private:
290 
291  typedef QuadratureTraits< FieldType, dimension > QuadratureTraitsType;
292 
293  public:
295  typedef typename QuadratureTraitsType :: PointQuadratureType PointQuadratureType;
296 
298  typedef typename QuadratureTraitsType :: IntegrationPointListType IntegrationPointListType;
299 
301  typedef typename QuadratureTraitsType :: QuadratureKeyType QuadratureKeyType;
302 
304  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
305  const QuadratureKeyType& quadKey )
306  {
307  assert( geometry.isCube() || geometry.isSimplex() );
308  static PointQuadratureType quad( geometry,
309  quadKey,
310  IdProvider ::instance().newId());
311  return quad;
312  }
313 
315  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
316  const GeometryType &elementGeometry,
317  const QuadratureKeyType& quadKey )
318  {
319  return getQuadrature(geometry, quadKey);
320  }
321 
322  QuadratureProvider() = delete;
323  QuadratureProvider( const ThisType& ) = delete;
324  QuadratureProvider &operator=( const ThisType& ) = delete;
325  };
326 
327 
328 
330  template< class FieldImp, template< class, int > class QuadratureTraits >
331  class QuadratureProvider< FieldImp, 1, QuadratureTraits >
332  {
333  public:
334  typedef FieldImp FieldType;
335 
336  enum { dimension = 1 };
337 
338  private:
340 
341  typedef QuadratureTraits< FieldType, dimension > QuadratureTraitsType;
342 
343  public:
345  typedef typename QuadratureTraitsType :: LineQuadratureType LineQuadratureType;
346 
348  typedef typename QuadratureTraitsType :: IntegrationPointListType IntegrationPointListType;
349 
351  typedef typename QuadratureTraitsType :: QuadratureKeyType QuadratureKeyType;
352 
354  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
355  const QuadratureKeyType& quadKey )
356  {
357  assert( geometry.isCube() || geometry.isSimplex() );
358  return QuadCreator< 0 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey );
359  }
360 
362  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
363  const GeometryType &elementGeometry,
364  const QuadratureKeyType& quadKey )
365  {
366  assert( geometry.isCube() || geometry.isSimplex() );
367  // we need here to distinguish between the basic types
368  // otherwise the this won't work for UGGrid
369  return ( elementGeometry.isSimplex() ) ?
370  QuadCreator< 0 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey ) :
371  QuadCreator< 1 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey ) ;
372  }
373 
374  QuadratureProvider() = delete;
375  QuadratureProvider( const ThisType& ) = delete;
376  QuadratureProvider &operator=( const ThisType& ) = delete;
377  };
378 
379 
380 
382  template< class FieldImp, template< class, int > class QuadratureTraits >
383  class QuadratureProvider< FieldImp, 2, QuadratureTraits >
384  {
385  public:
386  typedef FieldImp FieldType;
387 
388  enum { dimension = 2 };
389 
390  private:
392 
393  typedef QuadratureTraits< FieldType, dimension > QuadratureTraitsType;
394 
395  public:
397  typedef typename QuadratureTraitsType :: SimplexQuadratureType SimplexQuadratureType;
399  typedef typename QuadratureTraitsType :: CubeQuadratureType CubeQuadratureType;
400 
402  typedef typename QuadratureTraitsType :: IntegrationPointListType IntegrationPointListType;
403 
405  typedef typename QuadratureTraitsType :: QuadratureKeyType QuadratureKeyType;
406 
408  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
409  const QuadratureKeyType& quadKey )
410  {
411  assert( geometry.isCube() || geometry.isSimplex() || geometry.isNone() );
412 
413  if( geometry.isSimplex() )
414  {
415  return QuadCreator< 0 > ::
416  template provideQuad< SimplexQuadratureType > ( geometry, quadKey );
417  }
418  else if( geometry.isCube() )
419  {
420  return QuadCreator< 1 > ::
421  template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
422  }
423  else // type == None
424  {
425  // dummy return for polygonal grid cells, i.e. geometry type none
426  return QuadCreator< 1 > :: template provideQuad< CubeQuadratureType > ( geometry, 0 );
427  }
428  }
429 
431  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
432  const GeometryType &elementGeometry,
433  const QuadratureKeyType& quadKey )
434  {
435  assert( geometry.isCube() || geometry.isSimplex() );
436 
437  // if geometry is simplex return simplex quadrature
438  if ( geometry.isSimplex() )
439  {
440  // check element geometry to provide quadratures with different ids
441  if( elementGeometry.isSimplex() )
442  return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
443  else if( elementGeometry.isCube() )
444  return QuadCreator< 1 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
445  else if( elementGeometry.isPrism() )
446  return QuadCreator< 2 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
447  else if( elementGeometry.isPyramid() )
448  return QuadCreator< 3 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
449  else
450  DUNE_THROW( RangeError, "Element type not available for dimension 3" );
451  }
452  else
453  {
454  // return cube quadrature
455  // check element geometry to provide quadratures with different ids
456  if( elementGeometry.isSimplex() )
457  return QuadCreator< 4 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
458  else if( elementGeometry.isCube() )
459  return QuadCreator< 5 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
460  else if( elementGeometry.isPrism() )
461  return QuadCreator< 6 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
462  else if( elementGeometry.isPyramid() )
463  return QuadCreator< 7 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
464  else
465  DUNE_THROW( RangeError, "Element type not available for dimension 3" );
466  }
467 
468  DUNE_THROW( RangeError, "Element type not available for dimension 2" );
469  // dummy return
470  return QuadCreator< 0 > ::
471  template provideQuad< SimplexQuadratureType >( geometry, quadKey, 0 );
472  }
473 
474  QuadratureProvider() = delete;
475  QuadratureProvider( const ThisType& ) = delete;
476  QuadratureProvider &operator=( const ThisType& ) = delete;
477  };
478 
479 
480 
482  template< class FieldImp, template< class, int > class QuadratureTraits >
483  class QuadratureProvider< FieldImp, 3, QuadratureTraits >
484  {
485  public:
486  typedef FieldImp FieldType;
487 
488  enum { dimension = 3 };
489 
490  private:
492 
493  typedef QuadratureTraits< FieldType, dimension > QuadratureTraitsType;
494 
495  public:
497  typedef typename QuadratureTraitsType :: SimplexQuadratureType SimplexQuadratureType;
499  typedef typename QuadratureTraitsType :: CubeQuadratureType CubeQuadratureType;
501  typedef typename QuadratureTraitsType :: PrismQuadratureType PrismQuadratureType;
503  typedef typename QuadratureTraitsType :: PyramidQuadratureType PyramidQuadratureType;
504 
506  typedef typename QuadratureTraitsType :: IntegrationPointListType IntegrationPointListType;
507 
509  typedef typename QuadratureTraitsType :: QuadratureKeyType QuadratureKeyType;
510 
512  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
513  const QuadratureKeyType& quadKey )
514  {
515  assert( geometry.isCube() || geometry.isSimplex() || geometry.isNone()
516  || geometry.isPrism() || geometry.isPyramid() );
517 
518  if( geometry.isSimplex() )
519  return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType >
520  ( geometry, quadKey );
521  if( geometry.isCube() )
522  return QuadCreator< 1 > :: template provideQuad< CubeQuadratureType >
523  ( geometry, quadKey );
524 
525  if( geometry.isPrism() )
526  return QuadCreator< 2 > :: template provideQuad< PrismQuadratureType >
527  ( geometry, quadKey );
528  if( geometry.isPyramid() )
529  return QuadCreator< 3 > :: template provideQuad< PyramidQuadratureType >
530  ( geometry, quadKey );
531 
532  if( geometry.isNone() )
533  {
534  // dummy return for polyhedral grid cells
535  return QuadCreator< 1 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey, 0 );
536  }
537 
538  DUNE_THROW( RangeError, "Element type not available for dimension 3" );
539  // dummy return
540  return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType >
541  ( geometry, quadKey, 0 );
542  }
543 
544  static const IntegrationPointListType &getQuadrature( const GeometryType &geometry,
545  const GeometryType &elementGeometry,
546  const QuadratureKeyType& quadKey )
547  {
548  DUNE_THROW( RangeError, "QuadProvider::getQuadrature not implemented for 3d face quadratures!" );
549  // dummy return
550  return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType >
551  ( geometry, quadKey, 0 );
552  }
553 
554  QuadratureProvider() = delete;
555  QuadratureProvider( const ThisType& ) = delete;
556  QuadratureProvider &operator=( const ThisType& ) = delete;
557  };
558 
559  } // namespace Fem
560 
561 } // namespace Dune
562 
563 #endif // #ifndef DUNE_FEM_QUADPROVIDER_HH
Definition: bindguard.hh:11
Definition: pointmapper.hh:18
static IdProvider & instance()
Access to the singleton object.
Definition: idprovider.hh:21
A simple quadrature key class for use FemPy.
Definition: quadprovider.hh:25
FemQuadratureKey(const FemQuadratureKey &key)=default
copy constructor
int first() const
return first component
Definition: quadprovider.hh:62
int second() const
return second component
Definition: quadprovider.hh:64
FemQuadratureKey()
empty constructor
Definition: quadprovider.hh:34
FemQuadratureKey(const int first)
constructor taking only order (fallback for standard Fem quadratures)
Definition: quadprovider.hh:52
static const int highest_order
Definition: quadprovider.hh:41
FemQuadratureKey(const int first, const int second)
constructor taking to ids, like std::pair
Definition: quadprovider.hh:44
the actual quadrature storage
Definition: quadprovider.hh:79
static const QuadImp & provideQuad(const GeometryType &geometry, const QuadratureKey &key)
provide quadrature
Definition: quadprovider.hh:176
static const QuadImp & provideQuad(const GeometryType &geometry, const int, const int defaultOrder)
provide quadrature
Definition: quadprovider.hh:208
static const QuadImp & provideQuad(const GeometryType &geometry, const QuadratureKey &key, const int defaultOrder)
provide quadrature
Definition: quadprovider.hh:190
provide a single instance pool of quadratures
Definition: quadprovider.hh:235
@ dimension
Definition: quadprovider.hh:239
FieldImp FieldType
Definition: quadprovider.hh:237
QuadratureProvider & operator=(const ThisType &)=delete
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:265
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:258
QuadratureTraitsType ::CubeQuadratureType CubeQuadratureType
type for cube quadrature
Definition: quadprovider.hh:248
QuadratureTraitsType ::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:255
QuadratureTraitsType ::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:252
QuadratureProvider(const ThisType &)=delete
QuadratureProvider & operator=(const ThisType &)=delete
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:304
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:315
QuadratureTraitsType ::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:298
QuadratureTraitsType ::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:301
QuadratureTraitsType ::PointQuadratureType PointQuadratureType
type of point quadrature
Definition: quadprovider.hh:295
FieldImp FieldType
Definition: quadprovider.hh:284
FieldImp FieldType
Definition: quadprovider.hh:334
QuadratureTraitsType ::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:348
QuadratureTraitsType ::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:351
QuadratureProvider & operator=(const ThisType &)=delete
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:362
QuadratureTraitsType ::LineQuadratureType LineQuadratureType
type of line quadrature
Definition: quadprovider.hh:345
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:354
QuadratureTraitsType ::CubeQuadratureType CubeQuadratureType
type of cube quadrature
Definition: quadprovider.hh:399
QuadratureTraitsType ::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:402
FieldImp FieldType
Definition: quadprovider.hh:386
QuadratureTraitsType ::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:405
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:408
QuadratureProvider & operator=(const ThisType &)=delete
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:431
QuadratureTraitsType ::SimplexQuadratureType SimplexQuadratureType
type of simplex quadrature
Definition: quadprovider.hh:397
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:512
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Definition: quadprovider.hh:544
QuadratureTraitsType ::CubeQuadratureType CubeQuadratureType
type of cube quadrature
Definition: quadprovider.hh:499
FieldImp FieldType
Definition: quadprovider.hh:486
QuadratureTraitsType ::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:509
QuadratureTraitsType ::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:506
QuadratureProvider & operator=(const ThisType &)=delete
QuadratureTraitsType ::PrismQuadratureType PrismQuadratureType
type of prims quadrature
Definition: quadprovider.hh:501
QuadratureTraitsType ::SimplexQuadratureType SimplexQuadratureType
type of simplex quadrature
Definition: quadprovider.hh:497
QuadratureTraitsType ::PyramidQuadratureType PyramidQuadratureType
type of pyramid quadrature
Definition: quadprovider.hh:503