DUNE-FEM (unstable)

cachingquadrature.hh
1#ifndef DUNE_FEM_CACHINGQUADRATURE_HH
2#define DUNE_FEM_CACHINGQUADRATURE_HH
3
4//- Dune includes
5#include <dune/common/math.hh>
6
7//- Local includes
8#include <dune/fem/quadrature/elementquadrature.hh>
9#include <dune/fem/quadrature/caching/twistutility.hh>
10#include <dune/fem/quadrature/caching/pointmapper.hh>
11#include <dune/fem/quadrature/caching/cacheprovider.hh>
12
13namespace Dune
14{
15
16 namespace Fem
17 {
18
23 {
24 protected:
25 // do not create instances of this class
27 {
28 }
29
30 public:
32 static constexpr bool twisted () { return false; }
33
35 inline int twistId () const { return 0; }
36
47 inline size_t cachingPoint( const size_t quadraturePoint ) const
48 {
50 "CachingInterface :: cachingPoint must be overloaded!" );
51 }
52
58 inline size_t interpolationPoint( const size_t quadraturePoint ) const
59 {
61 "CachingInterface :: interpolationPoint must be overloaded!" );
62 }
63
70 inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
71 {
73 "CachingInterface :: isInterpolationQuadrature must be overloaded!" );
74 }
75 };
76
77
78
100 template< class GridPartImp, int codim, class IntegrationTraits, bool isQuadrature >
102
103
105 template< class GridPartImp, class IntegrationTraits, bool isQuadrature >
106 class CachingQuadratureImpl< GridPartImp, 0, IntegrationTraits, isQuadrature >
107 : public ElementPointListBase< GridPartImp, 0, IntegrationTraits >,
108 public CachingInterface
109 {
112
113 public:
114 typedef typename Base::GridPartType GridPartType;
115 typedef typename Base::EntityType EntityType;
116 static const int codimension = Base::codimension;
117
119 typedef typename Base::CoordinateType CoordinateType;
120
125
127 static const int pointSetId = SelectQuadraturePointSetId<
128 typename IntegrationTraits::IntegrationPointListType::Traits > :: value;
129
130 protected:
131 using Base::quadImp;
132
133 public:
136 using Base::localPoint;
137 using Base::nop;
138
141 template <class QuadratureKeyType>
142 CachingQuadratureImpl( const EntityType& entity, const QuadratureKeyType& quadKey )
143 : Base( entity.type(), quadKey )
144 {
145 CacheProvider< GridPartType, codimension >::registerQuadrature( quadImp() );
146 }
147
150 template <class QuadratureKeyType>
151 CachingQuadratureImpl( const GeometryType &geometry, const QuadratureKeyType& quadKey )
152 : Base( geometry, quadKey )
153 {
154 CacheProvider< GridPartType, codimension >::registerQuadrature( quadImp() );
155 }
156
157 const QuadraturePointWrapperType operator[] ( const size_t i ) const
158 {
159 return QuadraturePointWrapperType( *this, i );
160 }
161
162 IteratorType begin () const noexcept { return IteratorType( *this, 0 ); }
163 IteratorType end () const noexcept { return IteratorType( *this, nop() ); }
164
166 const CoordinateType &point ( const size_t i ) const
167 {
168 return localPoint( i );
169 }
170
172 auto weight ( std::size_t i ) const
173 {
174 return quadImp().weight( i );
175 }
176
178 inline size_t cachingPoint( const size_t quadraturePoint ) const
179 {
180 return quadraturePoint;
181 }
182
184 inline size_t interpolationPoint( const size_t quadraturePoint ) const
185 {
186 return quadraturePoint;
187 }
188
190 inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
191 {
192 // if pointSetId is not negative then we have an interpolation
193 // quadrature if the number of point are equal to number of shape functions
194 return (pointSetId >= 0) ? (nop() == numShapeFunctions) : false;
195 }
196 };
197
198
199
201 template< typename GridPartImp, class IntegrationTraits, bool isQuadrature >
202 class CachingQuadratureImpl< GridPartImp, 1, IntegrationTraits, isQuadrature >
203 : public ElementPointListBase< GridPartImp, 1, IntegrationTraits >,
204 public CachingInterface
205 {
208
209 public:
211 typedef GridPartImp GridPartType;
212
213 typedef typename Base::RealType RealType;
214 static const int dimension = Base::dimension;
215 static const int codimension = Base::codimension;
216
218 typedef typename Base::CoordinateType CoordinateType;
219
222 typedef typename IntersectionIteratorType::Intersection IntersectionType;
223
227
229 typedef ElementQuadratureImpl< GridPartType, codimension, IntegrationTraits, isQuadrature >
231
232
233 // for compatibility
235 typedef IntersectionIteratorType IntersectionIterator;
236
237 private:
238 static const int quadPointSetId =
240
241 public:
242 // Note: we also exclude GaussLegendre(0) here, because on faces it is not
243 // an interpolation rule
244 static const int pointSetId = (quadPointSetId > 0) ? quadPointSetId :
246
247 protected:
248 typedef typename CachingTraits< RealType, dimension >::MapperPairType MapperPairType;
249 typedef typename CachingTraits< RealType, dimension >::PointVectorType PointVectorType;
250
251 typedef CacheProvider< GridPartType, codimension > CacheProviderType;
252 typedef PointProvider< RealType, dimension, codimension> PointProviderType;
253
254 using Base::quadImp;
255
256 public:
257 using Base::localFaceIndex;
258 using Base::elementGeometry;
259 using Base::nop;
260
273 template <class QuadratureKeyType>
275 const IntersectionType &intersection,
276 const QuadratureKeyType& quadKey, const typename Base :: Side side )
277 : Base( getPointList( intersection, quadKey, side ) ),
278 side_(side),
279 twist_( getTwist( gridPart, intersection, side ) ),
280 mapper_( CacheProviderType::getMapper( quadImp(), elementGeometry(), localFaceIndex(), twist_) ),
281 points_( PointProviderType::getPoints( quadImp().ipList().id(), elementGeometry() ) ),
282 intersection_(intersection)
283 {
284 }
285
286 const QuadraturePointWrapperType operator[] ( const size_t i ) const
287 {
288 return QuadraturePointWrapperType( *this, i );
289 }
290
291 IteratorType begin () const noexcept { return IteratorType( *this, 0 ); }
292 IteratorType end () const noexcept { return IteratorType( *this, nop() ); }
293
294 typename Base :: Side side() const { return side_; }
295 bool isInside() const { return side_ == Base::INSIDE; }
296
299 const CoordinateType &point ( const size_t i ) const
300 {
301 return points_[ cachingPoint( i ) ];
302 }
303
305 auto weight ( std::size_t i ) const
306 {
307 return quadImp().weight( i );
308 }
309
310 const IntersectionType &intersection() const
311 {
312 return intersection_;
313 }
314
316 static constexpr bool twisted() { return true; }
317
319 inline int twistId () const { return twist_ + 4; }
320
322 inline size_t cachingPoint( const size_t quadraturePoint ) const
323 {
324 assert( quadraturePoint < (size_t)nop() );
325 return mapper_.first[ quadraturePoint ];
326 }
327
329 inline size_t interpolationPoint( const size_t quadraturePoint ) const
330 {
331 assert( quadraturePoint < mapper_.second.size() );
332 return mapper_.second[ quadraturePoint ];
333 }
334
336 inline bool isInterpolationQuadrature( const size_t numShapeFunctions ) const
337 {
338 // if pointSetId is not negative then we have an interpolation
339 // quadrature if the number of point are equal to number of shape functions
340 return (pointSetId < 0) ? false :
341 quadImp().ipList().isFaceInterpolationQuadrature( numShapeFunctions );
342 }
343
344 // return local caching point
345 // for debugging issues only
346 size_t localCachingPoint ( const size_t i ) const
347 {
348 const auto& mapper = mapper_.first;
349
350 assert( i < (size_t)nop() );
351
352 assert( mapper[ i ] >= 0 );
353 int faceIndex = localFaceIndex();
354 unsigned int point = mapper[ i ] - faceIndex * mapper.size();
355 assert( point < nop() );
356
357 return point;
358 }
359
360 protected:
361 template <class QuadratureKeyType>
362 Base getPointList ( const IntersectionType &intersection,
363 const QuadratureKeyType& key,
364 const typename Base :: Side side )
365 {
366 switch( side )
367 {
368 case Base :: INSIDE:
369 return Base( TwistUtilityType::elementGeometry( intersection, true ),
370 intersection.indexInInside(), key );
371
372 case Base :: OUTSIDE:
373 return Base( TwistUtilityType::elementGeometry( intersection, false ),
374 intersection.indexInOutside(), key );
375
376 default:
377 DUNE_THROW( InvalidStateException, "ElementIntegrationPointList: side must either be INSIDE or OUTSIDE." );
378 }
379 }
380
381 int getTwist ( const GridPartType &gridPart,
382 const IntersectionType &intersection,
383 const typename Base :: Side side )
384 {
385 switch( side )
386 {
387 case Base :: INSIDE:
388 return TwistUtilityType::twistInSelf( gridPart.grid(), intersection );
389
390 case Base :: OUTSIDE:
391 return TwistUtilityType::twistInNeighbor( gridPart.grid(), intersection );
392
393 default:
394 DUNE_THROW( InvalidStateException, "ElementIntegrationPointList: side must either be INSIDE or OUTSIDE." );
395 }
396 }
397
398 private:
399 const typename Base :: Side side_;
400 const int twist_;
401 const MapperPairType &mapper_;
402 const PointVectorType &points_;
403 const IntersectionType &intersection_;
404 };
405
406
408 template< typename GridPart, int codim, class IntegrationTraits >
409 using CachingPointList = CachingQuadratureImpl< GridPart, codim, IntegrationTraits, false>;
410
412 template< typename GridPart, int codim, template <class, int> class QuadratureTraits = DefaultQuadratureTraits >
413 using CachingQuadrature = CachingQuadratureImpl< GridPart, codim, ElementQuadratureTraits< GridPart, codim, QuadratureTraits >, true>;
414
415 template<class GridPart, class Entity>
416 static inline auto cachingQuadrature(const GridPart& gridPart, const Entity& entity, unsigned quadOrder)
417 {
419 return Quadrature(entity, quadOrder);
420 }
421
422 } //namespace Fem
423
424} //namespace Dune
425
426#endif // #ifndef DUNE_FEM_CACHINGQUADRATURE_HH
Traits::IntersectionIteratorType IntersectionIteratorType
type of intersection iterator
Definition: adaptiveleafgridpart.hh:92
interface a cachable quadrature has to implement
Definition: cachingquadrature.hh:23
size_t interpolationPoint(const size_t quadraturePoint) const
map quadrature points to interpolation points
Definition: cachingquadrature.hh:58
bool isInterpolationQuadrature(const size_t numShapeFunctions) const
check if quadrature is interpolation quadrature
Definition: cachingquadrature.hh:70
static constexpr bool twisted()
returns true if cachingPoint is not the identity mapping
Definition: cachingquadrature.hh:32
size_t cachingPoint(const size_t quadraturePoint) const
map quadrature points to caching points
Definition: cachingquadrature.hh:47
int twistId() const
returns the twistId, i.e. [0,...,7]
Definition: cachingquadrature.hh:35
size_t cachingPoint(const size_t quadraturePoint) const
map quadrature points to caching points
Definition: cachingquadrature.hh:178
auto weight(std::size_t i) const
obtain weight of i-th integration point (if quadrature, else 1.0)
Definition: cachingquadrature.hh:172
Base::CoordinateType CoordinateType
The type of the coordinates in the codim-0 reference element.
Definition: cachingquadrature.hh:119
size_t interpolationPoint(const size_t quadraturePoint) const
map quadrature points to interpolation points
Definition: cachingquadrature.hh:184
QuadraturePointWrapper< This > QuadraturePointWrapperType
the type of the quadrature point
Definition: cachingquadrature.hh:122
const CoordinateType & point(const size_t i) const
obtain coordinates of i-th integration point
Definition: cachingquadrature.hh:166
CachingQuadratureImpl(const EntityType &entity, const QuadratureKeyType &quadKey)
Definition: cachingquadrature.hh:142
bool isInterpolationQuadrature(const size_t numShapeFunctions) const
check if quadrature is interpolation quadrature
Definition: cachingquadrature.hh:190
CachingQuadratureImpl(const GeometryType &geometry, const QuadratureKeyType &quadKey)
Definition: cachingquadrature.hh:151
QuadraturePointIterator< This > IteratorType
type of iterator
Definition: cachingquadrature.hh:124
static constexpr bool twisted()
returns true if cachingPoint is not the identity mapping
Definition: cachingquadrature.hh:316
size_t interpolationPoint(const size_t quadraturePoint) const
map quadrature points to interpolation points
Definition: cachingquadrature.hh:329
GridPartType::IntersectionIteratorType IntersectionIteratorType
Type of the intersection iterator.
Definition: cachingquadrature.hh:221
const CoordinateType & point(const size_t i) const
obtain coordinates of i-th integration point
Definition: cachingquadrature.hh:299
bool isInterpolationQuadrature(const size_t numShapeFunctions) const
check if quadrature is interpolation quadrature
Definition: cachingquadrature.hh:336
size_t cachingPoint(const size_t quadraturePoint) const
map quadrature points to caching points
Definition: cachingquadrature.hh:322
int twistId() const
returns the twistId, i.e. [0,...,7]
Definition: cachingquadrature.hh:319
GridPartImp GridPartType
type of the grid partition
Definition: cachingquadrature.hh:211
ElementQuadratureImpl< GridPartType, codimension, IntegrationTraits, isQuadrature > NonConformingQuadratureType
type of quadrature used for non-conforming intersections
Definition: cachingquadrature.hh:230
CachingQuadratureImpl(const GridPartType &gridPart, const IntersectionType &intersection, const QuadratureKeyType &quadKey, const typename Base ::Side side)
constructor
Definition: cachingquadrature.hh:274
auto weight(std::size_t i) const
obtain weight of i-th integration point (if quadrature, else 1.0)
Definition: cachingquadrature.hh:305
QuadraturePointIterator< This > IteratorType
type of iterator
Definition: cachingquadrature.hh:226
Base::CoordinateType CoordinateType
Type of coordinates in codim-0 reference element.
Definition: cachingquadrature.hh:218
quadrature class supporting base function caching
Definition: cachingquadrature.hh:101
ElementPointListBase.
Definition: elementpointlistbase.hh:193
GridPartType::ctype RealType
coordinate type
Definition: elementpointlistbase.hh:207
iterator over quadrature points
Definition: quadrature.hh:103
wrapper for a (Quadrature,int) pair
Definition: quadrature.hh:42
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Default exception for dummy implementations.
Definition: exceptions.hh:263
actual interface class for quadratures
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
Utility to get twist from IntersectionIterator, if provided by grid (i.e. AlbertaGrid,...
Definition: twistutility.hh:107
selects Obj::pointSetId if available, otherwise defaultValue (default is -1)
Definition: utility.hh:174
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)