1#ifndef DUNE_FEM_QUADPROVIDER_HH
2#define DUNE_FEM_QUADPROVIDER_HH
10#include <dune/fem/quadrature/quadratureimp.hh>
11#include <dune/fem/quadrature/idprovider.hh>
12#include <dune/fem/misc/mpimanager.hh>
30 static const int maxFirst = 25 ;
32 static const int maxSecond = 10 ;
42 static const int highest_order = maxFirst * maxSecond ;
48 assert(
first < maxFirst );
49 assert(
second < maxSecond );
56 assert(
first < maxFirst );
60 operator int ()
const {
return id_; }
63 int first ()
const {
return id_ % highest_order ; }
65 int second ()
const {
return id_ / highest_order ; }
78 template<
unsigned int dummy >
82 template<
class QuadImp,
class QuadratureKey,
unsigned int >
83 class QuadratureStorage;
86 template<
class QuadImp,
class QuadratureKey >
87 class QuadratureStorage< QuadImp, QuadratureKey, dummy >
90 typedef QuadImp QuadType;
95 typedef std :: map< QuadratureKey, std::unique_ptr< QuadType > > StorageType;
99 static QuadImp* create(
const GeometryType &geometry,
const QuadratureKey& key )
101 return instance().createImpl( geometry, key );
105 static Storage& instance()
110 QuadImp* createImpl(
const GeometryType &geometry,
const QuadratureKey& key )
112 std::lock_guard< std::mutex > guard( mutex_ );
114 auto& quadPtr = storage_[ key ];
121 return quadPtr.operator->();
125 typedef std :: map< QuadratureKey, std::unique_ptr< QuadType, Dune::null_deleter<QuadType> > > PointerStorageType;
126 PointerStorageType quadPtrs_;
129 QuadratureStorage () {}
131 QuadType &getQuadrature(
const GeometryType &geometry,
const QuadratureKey& key )
133 auto& quadPtr = quadPtrs_[ key ];
137 quadPtr.reset( Storage::create( geometry, key ) );
146 template<
class QuadImp>
147 class QuadratureStorage< QuadImp, int, dummy >
150 typedef QuadImp QuadType;
158 std::vector< std::unique_ptr< QuadType > > storage_;
159 Storage() : storage_( QuadType :: maxOrder() + 1 )
162 static QuadType* create(
const GeometryType &geometry,
int order )
164 return instance().createImpl( geometry, order );
168 static Storage& instance()
173 QuadType* createImpl(
const GeometryType &geometry,
int order )
175 std::lock_guard< std::mutex > guard( mutex_ );
177 auto& quadPtr = storage_[ order ];
183 return quadPtr.operator ->();
187 typedef std::vector< QuadType* > QuadPointerVecType;
189 QuadPointerVecType quadPtrs_;
193 : quadPtrs_( QuadPointerVecType( QuadType :: maxOrder() + 1,
nullptr ))
198 QuadImp &getQuadrature(
const GeometryType &geometry,
int order )
200 if(order >=
int(quadPtrs_.size()) )
203 static thread_local bool showMessage = true ;
206 std::cerr <<
"WARNING: QuadratureStorage::getQuadrature: A quadrature of order " << order
207 <<
" is not implemented!" << std::endl
208 <<
"Choosing maximum order: " << quadPtrs_.size()-1 << std::endl << std::endl;
212 order = quadPtrs_.size() - 1;
215 auto& quadPtr = quadPtrs_[ order ];
219 quadPtr = Storage::create( geometry, order );
228 template<
class QuadImp >
230 :
public QuadratureStorage< QuadImp, int, dummy >
240 template<
class QuadImp,
class QuadratureKey >
242 const QuadratureKey& key )
246 static thread_local QuadratureStorage< QuadImp, QuadratureKey, dummy > storage;
247 return storage.getQuadrature( geometry, key );
256 template<
class QuadImp,
class QuadratureKey >
258 const QuadratureKey& key,
259 const int defaultOrder )
262 assert( geometry.
isNone() );
263 DUNE_THROW(
NotImplemented,
"provideQuad for polyhedral cells (defaultOrder = 0) not implemented for arbitrary QuadratureKey!");
264 QuadImp* ptr =
nullptr;
274 template<
class QuadImp >
277 const int defaultOrder )
279 assert( geometry.
isNone() );
282 static thread_local QuadratureStorage< QuadImp, int, dummy > storage;
283 return storage.getQuadrature( geometry, defaultOrder );
302 template<
typename FieldImp,
int dim,
template<
class,
int >
class IntegrationTraits >
306 typedef FieldImp FieldType;
308 static const int dimension = dim;
313 typedef IntegrationTraits< FieldType, dimension > QuadratureTraits;
316 typedef QuadratureTraits FactoryTraits;
318 template <
class Po
intQuadrature,
class QuadratureKey>
319 class PointQuadratureStorage :
public PointQuadrature
323 PointQuadratureStorage(
const GeometryType &geometry,
const QuadratureKey& quadKey )
334 template <
class FactoryTraits>
337 const typename FactoryTraits::QuadratureKeyType& quadKey )
340 const typename FactoryTraits::QuadratureKeyType& quadKey )
344 typename FactoryTraits::IntegrationPointListType> :: value );
347 if constexpr ( dimension == 0 )
349 typedef typename FactoryTraits :: PointQuadratureType PointQuadratureType;
350 typedef typename FactoryTraits :: QuadratureKeyType QuadratureKeyType;
352 typedef PointQuadratureStorage< PointQuadratureType, QuadratureKeyType > PointQuadratureStorageType;
360 if constexpr ( dimension == 1 )
362 typedef typename FactoryTraits::LineQuadratureType LineQuadratureType;
364 return QuadCreator< 0 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey );
368 if constexpr ( dimension >= 2 )
370 typedef typename FactoryTraits::SimplexQuadratureType SimplexQuadratureType;
371 typedef typename FactoryTraits::CubeQuadratureType CubeQuadratureType;
375 template provideQuad< SimplexQuadratureType > ( geometry, quadKey );
381 template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
387 return QuadCreator< 1 > :: template provideQuad< CubeQuadratureType > ( geometry, 0 );
390 if constexpr ( dimension == 3 )
392 typedef typename FactoryTraits::PrismQuadratureType PrismQuadratureType;
393 typedef typename FactoryTraits::PyramidQuadratureType PyramidQuadratureType;
397 ( geometry, quadKey );
400 ( geometry, quadKey );
406 if constexpr ( dimension == 1 )
408 typedef typename FactoryTraits::LineQuadratureType LineQuadratureType;
409 return QuadCreator< 0 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey, 0 );
413 typedef typename FactoryTraits::SimplexQuadratureType SimplexQuadratureType;
414 return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey, 0 );
421 template <
class FactoryTraits>
425 const typename FactoryTraits::QuadratureKeyType& quadKey )
430 const typename FactoryTraits::QuadratureKeyType& quadKey )
433 if constexpr ( dimension == 0 )
439 if constexpr ( dimension == 1 )
441 typedef typename FactoryTraits::LineQuadratureType LineQuadratureType;
446 QuadCreator< 0 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey ) :
447 QuadCreator< 1 > :: template provideQuad< LineQuadratureType > ( geometry, quadKey ) ;
450 if constexpr ( dimension == 2 )
452 typedef typename FactoryTraits::SimplexQuadratureType SimplexQuadratureType;
453 typedef typename FactoryTraits::CubeQuadratureType CubeQuadratureType;
461 return QuadCreator< 0 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
462 else if( elementGeometry.
isCube() )
463 return QuadCreator< 1 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
464 else if( elementGeometry.
isPrism() )
465 return QuadCreator< 2 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
467 return QuadCreator< 3 > :: template provideQuad< SimplexQuadratureType > ( geometry, quadKey ) ;
476 return QuadCreator< 4 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
477 else if( elementGeometry.
isCube() )
478 return QuadCreator< 5 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
479 else if( elementGeometry.
isPrism() )
480 return QuadCreator< 6 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
482 return QuadCreator< 7 > :: template provideQuad< CubeQuadratureType > ( geometry, quadKey ) ;
488 if constexpr ( dimension >= 3 )
490 typedef typename FactoryTraits::SimplexQuadratureType SimplexQuadratureType;
491 DUNE_THROW(
RangeError,
"QuadProvider::getQuadrature not implemented for 3d face quadratures!" );
494 ( geometry, quadKey, 0 );
A simple quadrature key class for use FemPy.
Definition: quadprovider.hh:26
FemQuadratureKey(const FemQuadratureKey &key)=default
copy constructor
int first() const
return first component
Definition: quadprovider.hh:63
int second() const
return second component
Definition: quadprovider.hh:65
FemQuadratureKey()
empty constructor
Definition: quadprovider.hh:35
FemQuadratureKey(const int first)
constructor taking only order (fallback for standard Fem quadratures)
Definition: quadprovider.hh:53
FemQuadratureKey(const int first, const int second)
constructor taking to ids, like std::pair
Definition: quadprovider.hh:45
static IdProvider & instance()
Access to the singleton object.
Definition: idprovider.hh:21
the actual quadrature storage
Definition: quadprovider.hh:80
static const QuadImp & provideQuad(const GeometryType &geometry, const QuadratureKey &key, const int defaultOrder)
provide quadrature
Definition: quadprovider.hh:257
static const QuadImp & provideQuad(const GeometryType &geometry, const QuadratureKey &key)
provide quadrature
Definition: quadprovider.hh:241
static const QuadImp & provideQuad(const GeometryType &geometry, const int, const int defaultOrder)
provide quadrature
Definition: quadprovider.hh:275
provide a single instance pool of quadratures
Definition: quadprovider.hh:304
static const QuadratureImplementationType & getQuadrature(const GeometryType &geometry, const GeometryType &elementGeometry, const typename FactoryTraits::QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:428
QuadratureTraits::IntegrationPointListType QuadratureImplementationType
type of integration point list implementation
Definition: quadprovider.hh:330
static const QuadratureImplementationType & getQuadrature(const GeometryType &geometry, const typename FactoryTraits::QuadratureKeyType &quadKey)
Access to the quadrature implementations.
Definition: quadprovider.hh:339
static DUNE_EXPORT Object & instance(Args &&... args)
return singleton instance of given Object type.
Definition: singleton.hh:123
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isPyramid() const
Return true if entity is a pyramid.
Definition: type.hh:304
constexpr bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:309
constexpr bool isCube() const
Return true if entity is a cube of any dimension.
Definition: type.hh:324
constexpr bool isNone() const
Return true if entity is a singular of any dimension.
Definition: type.hh:355
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:319
Default exception for dummy implementations.
Definition: exceptions.hh:263
Default exception class for range errors.
Definition: exceptions.hh:254
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
This file implements several utilities related to std::shared_ptr.