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
13namespace 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
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
273 QuadratureProvider( 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
323 QuadratureProvider( 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
375 QuadratureProvider( 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
475 QuadratureProvider( 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
555 QuadratureProvider( 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, const int defaultOrder)
provide quadrature
Definition: quadprovider.hh:190
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
provide a single instance pool of quadratures
Definition: quadprovider.hh:235
@ dimension
Definition: quadprovider.hh:239
FieldImp FieldType
Definition: quadprovider.hh:237
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
QuadratureProvider & operator=(const ThisType &)=delete
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
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
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:315
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:304
QuadratureTraitsType::PointQuadratureType PointQuadratureType
type of point quadrature
Definition: quadprovider.hh:295
QuadratureProvider & operator=(const ThisType &)=delete
FieldImp FieldType
Definition: quadprovider.hh:284
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:354
FieldImp FieldType
Definition: quadprovider.hh:334
QuadratureTraitsType::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:351
QuadratureProvider & operator=(const ThisType &)=delete
QuadratureTraitsType::LineQuadratureType LineQuadratureType
type of line quadrature
Definition: quadprovider.hh:345
QuadratureTraitsType::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:348
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:362
QuadratureProvider & operator=(const ThisType &)=delete
FieldImp FieldType
Definition: quadprovider.hh:386
QuadratureTraitsType::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:402
QuadratureTraitsType::CubeQuadratureType CubeQuadratureType
type of cube quadrature
Definition: quadprovider.hh:399
QuadratureTraitsType::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:405
QuadratureTraitsType::SimplexQuadratureType SimplexQuadratureType
type of simplex quadrature
Definition: quadprovider.hh:397
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:431
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:408
QuadratureProvider & operator=(const ThisType &)=delete
QuadratureTraitsType::PrismQuadratureType PrismQuadratureType
type of prims quadrature
Definition: quadprovider.hh:501
QuadratureTraitsType::CubeQuadratureType CubeQuadratureType
type of cube quadrature
Definition: quadprovider.hh:499
FieldImp FieldType
Definition: quadprovider.hh:486
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const QuadratureKeyType &quadKey)
Definition: quadprovider.hh:544
QuadratureTraitsType::QuadratureKeyType QuadratureKeyType
key for access of quadratures in the storage
Definition: quadprovider.hh:509
QuadratureTraitsType::SimplexQuadratureType SimplexQuadratureType
type of simplex quadrature
Definition: quadprovider.hh:497
QuadratureTraitsType::PyramidQuadratureType PyramidQuadratureType
type of pyramid quadrature
Definition: quadprovider.hh:503
QuadratureTraitsType::IntegrationPointListType IntegrationPointListType
type of integration point list implementation
Definition: quadprovider.hh:506
static const IntegrationPointListType & getQuadrature(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:512