dune-fem 2.8.0
Loading...
Searching...
No Matches
caching.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_CACHING_HH
2#define DUNE_FEM_SPACE_SHAPEFUNCTIONSET_CACHING_HH
3
4// C++ includes
5#include <cstddef>
6#include <vector>
7#include <type_traits>
8
9#include <dune/geometry/quadraturerules.hh>
10
11// dune-fem includes
18
20
21namespace Dune
22{
23
24 namespace Fem
25 {
26
27 // CachingShapeFunctionSet
28 // -----------------------
29
30 template< class ShapeFunctionSet >
33 {
35
36 public:
38
40
45
46 typedef std::vector< RangeType > RangeVectorType ;
47 typedef std::vector< JacobianRangeType > JacobianRangeVectorType ;
48
49 typedef std::vector< RangeVectorType > RangeCacheVectorType;
50 typedef std::vector< JacobianRangeVectorType > JacobianCacheVectorType;
51
52 public:
53 // point set id if available (otherwise -1)
54 static const int pointSetId = detail::SelectPointSetId< ShapeFunctionSetType >::value;
55
56 explicit CachingShapeFunctionSet ( const GeometryType &type,
57 const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() )
58 : type_( type )
59 , shapeFunctionSet_( shapeFunctionSet )
60 {
62 }
63
65
66 int order () const
67 {
68 return shapeFunctionSet_.order();
69 }
70
71 std::size_t size () const
72 {
73 return shapeFunctionSet_.size();
74 }
75
76 template< class Point, class Functor >
77 void evaluateEach ( const Point &x, Functor functor ) const
78 {
79 return shapeFunctionSet_.evaluateEach( x, functor );
80 }
81
82 template< class Quadrature, class Functor >
83 void evaluateEach ( const QuadraturePointWrapper< Quadrature > &x, Functor functor ) const
84 {
85 const bool cacheable = std::is_convertible< Quadrature, CachingInterface >::value;
86 evaluateEach( x.quadrature(), x.index(), functor, std::integral_constant< bool, cacheable >() );
87 }
88
89 template< class Point, class Functor >
90 void jacobianEach ( const Point &x, Functor functor ) const
91 {
92 return shapeFunctionSet_.jacobianEach( x, functor );
93 }
94
95 template< class Quadrature, class Functor >
96 void jacobianEach ( const QuadraturePointWrapper< Quadrature > &x, Functor functor ) const
97 {
98 const bool cacheable = std::is_convertible< Quadrature, CachingInterface >::value;
99 jacobianEach( x.quadrature(), x.index(), functor, std::integral_constant< bool, cacheable >() );
100 }
101
102 template< class Point, class Functor >
103 void hessianEach ( const Point &x, Functor functor ) const
104 {
105 return shapeFunctionSet_.hessianEach( x, functor );
106 }
107
108 GeometryType type () const { return type_; }
109
110 template < class QuadratureType >
111 const RangeVectorType& rangeCache( const QuadratureType& quadrature ) const
112 {
113 return ReturnCache< QuadratureType, std::is_base_of< CachingInterface, QuadratureType >::value > ::
114 ranges( *this, quadrature, rangeCaches_, localRangeCache_ );
115 }
116
117 template < class QuadratureType >
118 const JacobianRangeVectorType& jacobianCache( const QuadratureType& quadrature ) const
119 {
120 return ReturnCache< QuadratureType, std::is_base_of< CachingInterface, QuadratureType >::value > ::
121 jacobians( *this, quadrature, jacobianCaches_, localJacobianCache_ );
122 }
123
124 const ThisType& scalarShapeFunctionSet() const { return *this; }
125 const ThisType& impl() const { return *this; }
126
127 private:
128 template< class Quad, bool cacheable /* false */ >
129 struct ReturnCache
130 {
131 static const RangeVectorType&
132 ranges( const ThisType& shapeFunctionSet,
133 const Quad& quad,
135 RangeVectorType& storage )
136 {
137 // this feature was disabled in default_codegen.hh
138 // this method should therefore not be called.
139 assert( false );
140 std::abort();
141
142 // evaluate all basis functions and multiply with dof value
143 const unsigned int nop = quad.nop();
144 const unsigned int size = shapeFunctionSet.size();
145
146 // make sure cache has the appropriate size
147 storage.resize( size * nop );
148 RangeType* data = storage.data();
149
150 for( unsigned int qp = 0 ; qp < nop; ++ qp )
151 {
152 const int cacheQp = quad.cachingPoint( qp );
153 AssignFunctor< RangeType* > funztor( data + (cacheQp * size) );
154 shapeFunctionSet.evaluateEach( quad[ qp ], funztor );
155 }
156 return storage;
157 }
158
159 static const JacobianRangeVectorType&
160 jacobians( const ThisType& shapeFunctionSet,
161 const Quad& quad,
163 JacobianRangeVectorType& storage )
164 {
165 // this feature was disabled in default_codegen.hh
166 // this method should therefore not be called.
167 assert( false );
168 std::abort();
169
170 // evaluate all basis functions and multiply with dof value
171 const unsigned int nop = quad.nop();
172 const unsigned int size = shapeFunctionSet.size();
173
174 // make sure cache has the appropriate size
175 storage.resize( size * nop );
176 JacobianRangeType* data = storage.data();
177
178 for( unsigned int qp = 0 ; qp < nop; ++ qp )
179 {
180 const int cacheQp = quad.cachingPoint( qp );
181 AssignFunctor< JacobianRangeType* > funztor( data + ( cacheQp * size ) );
182 shapeFunctionSet.jacobianEach( quad[ qp ], funztor );
183 }
184 return storage;
185 }
186 };
187
188 template< class Quad >
189 struct ReturnCache< Quad, true >
190 {
191 static const RangeVectorType&
192 ranges( const ThisType& shapeFunctionSet,
193 const Quad& quad,
194 const RangeCacheVectorType& cache,
195 const RangeVectorType& )
196 {
197 return cache[ quad.id() ];
198 }
199
200 static const JacobianRangeVectorType&
201 jacobians( const ThisType& shapeFunctionSet,
202 const Quad& quad,
203 const JacobianCacheVectorType& cache,
205 {
206 return cache[ quad.id() ];
207 }
208 };
209
210
211 template< class Quadrature, class Functor >
212 void evaluateEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
213 std::integral_constant< bool, false > ) const
214 {
215 evaluateEach( quadrature.point( pt ), functor );
216 }
217
218 template< class Quadrature, class Functor >
219 void evaluateEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
220 std::integral_constant< bool, true > ) const;
221
222 template< class Quadrature, class Functor >
223 void jacobianEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
224 std::integral_constant< bool, false > ) const
225 {
226 jacobianEach( quadrature.point( pt ), functor );
227 }
228
229 template< class Quadrature, class Functor >
230 void jacobianEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
231 std::integral_constant< bool, true > ) const;
232
233
234 void cacheQuadrature( std::size_t id, std::size_t codim, std::size_t size );
235
236 template< class PointVector >
237 void cachePoints ( std::size_t id, const PointVector &points );
238
239 GeometryType type_;
240 ShapeFunctionSet shapeFunctionSet_;
241 RangeCacheVectorType rangeCaches_;
242 JacobianCacheVectorType jacobianCaches_;
243
244 mutable RangeVectorType localRangeCache_ ;
245 mutable JacobianRangeVectorType localJacobianCache_;
246 };
247
248
249
250 // Implementation of CachingShapeFunctionSet
251 // -----------------------------------------
252
253 template< class ShapeFunctionSet >
255 {
257 }
258
259
260 template< class ShapeFunctionSet >
261 template< class Quadrature, class Functor >
263 ::evaluateEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
264 std::integral_constant< bool, true > ) const
265 {
266 assert( (quadrature.id() < rangeCaches_.size()) && !rangeCaches_[ quadrature.id() ].empty() );
267 const RangeType *cache = rangeCaches_[ quadrature.id() ].data();
268
269 const unsigned int numShapeFunctions = size();
270 const unsigned int cpt = quadrature.cachingPoint( pt );
271
272 // for Lagrange-type basis evaluated on interpolation points
273 // this is the Kronecker delta, there we only need
274 // to evaluate the shapefunction with number 'pt'
275 static const int quadPointSetId = SelectQuadraturePointSetId< Quadrature >::value;
276 //std::cout << "QP:" << quadPointSetId << " " << pointSetId << std::endl;
277
278 if constexpr ( quadPointSetId == pointSetId )
279 {
280 // std::cout << "QP matches: " << quadPointSetId << " " << pointSetId << " " << quadrature.nop() << " " << numShapeFunctions << std::endl;
281 if( quadrature.isInterpolationQuadrature(numShapeFunctions) )
282 {
283 // negative values mean invalid point sets
284 // we should not get here in this case
285 assert( quadPointSetId >= 0 );
286 assert( pointSetId >= 0 );
287
288 //if( Quadrature::codimension == 1 )
289 // std::cout << "using interpolation point " << std::endl;
290 // obtain interpolation point (different for face quadratures)
291 const unsigned int point = quadrature.interpolationPoint( pt );
292#ifndef NDEBUG
293 for( unsigned int i = 0; i < numShapeFunctions; ++i )
294 assert( (cache[ cpt*numShapeFunctions + i ] - RangeType(i==point)).two_norm() < 1e-8 ) ;
295#endif
296 // point should be 1
297 functor( point, RangeType(1) );
298 return;
299 }
300 }
301
302 for( unsigned int i = 0; i < numShapeFunctions; ++i )
303 functor( i, cache[ cpt*numShapeFunctions + i ] );
304 }
305
306
307 template< class ShapeFunctionSet >
308 template< class Quadrature, class Functor >
309 inline void CachingShapeFunctionSet< ShapeFunctionSet >
310 ::jacobianEach ( const Quadrature &quadrature, std::size_t pt, Functor functor,
311 std::integral_constant< bool, true > ) const
312 {
313 assert( (quadrature.id() < jacobianCaches_.size()) && !jacobianCaches_[ quadrature.id() ].empty() );
314 const JacobianRangeType *cache = jacobianCaches_[ quadrature.id() ].data();
315
316 const unsigned int numShapeFunctions = size();
317 const unsigned int cpt = quadrature.cachingPoint( pt );
318 for( unsigned int i = 0; i < numShapeFunctions; ++i )
319 functor( i, cache[ cpt*numShapeFunctions + i ] );
320 }
321
322
323 template< class ShapeFunctionSet >
324 inline void CachingShapeFunctionSet< ShapeFunctionSet >
325 ::cacheQuadrature( std::size_t id, std::size_t codim, std::size_t size )
326 {
328 {
329 DUNE_THROW(SingleThreadModeError,"CachingShapeFunctionSet::cacheQuadrature: only call in single thread mode!");
330 }
331
332 if( id >= rangeCaches_.size() )
333 {
334 rangeCaches_.resize( id+1, RangeVectorType() );
335 jacobianCaches_.resize( id+1, JacobianRangeVectorType() );
336 }
337
338 assert( rangeCaches_[ id ].empty() == jacobianCaches_[ id ].empty() );
339
340 if( rangeCaches_[ id ].empty() )
341 {
342 typedef typename FunctionSpaceType::DomainFieldType ctype;
343 const int dim = FunctionSpaceType::dimDomain;
344 switch( codim )
345 {
346 case 0:
347 cachePoints( id, PointProvider< ctype, dim, 0 >::getPoints( id, type_ ) );
348 break;
349
350 case 1:
351 cachePoints( id, PointProvider< ctype, dim, 1 >::getPoints( id, type_ ) );
352 break;
353
354 default:
355 DUNE_THROW( NotImplemented, "Caching for codim > 1 not implemented." );
356 }
357 }
358 }
359
360
361 template< class ShapeFunctionSet >
362 template< class PointVector >
363 inline void CachingShapeFunctionSet< ShapeFunctionSet >
364 ::cachePoints ( std::size_t id, const PointVector &points )
365 {
366 const unsigned int numShapeFunctions = size();
367 const unsigned int numPoints = points.size();
368
369 RangeVectorType& ranges = rangeCaches_[ id ];
370 ranges.resize( numShapeFunctions * numPoints );
371
372 JacobianRangeVectorType& jacobians = jacobianCaches_[ id ];
373 jacobians.resize( numShapeFunctions * numPoints );
374
375 if( ranges.empty() || jacobians.empty() )
376 DUNE_THROW( OutOfMemoryError, "Unable to allocate shape function set caches." );
377
378 for( unsigned int pt = 0; pt < numPoints; ++pt )
379 {
380 evaluateEach( points[ pt ], AssignFunctor< RangeType * >( ranges.data() + pt*numShapeFunctions ) );
381 jacobianEach( points[ pt ], AssignFunctor< JacobianRangeType * >( jacobians.data() + pt*numShapeFunctions ) );
382 }
383 }
384
385 } // namespace Fem
386
387} // namespace Dune
388
389#endif // #ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_CACHING_HH
Definition: bindguard.hh:11
detail::SelectPointSetId< Quadrature, -Dune::QuadratureType::size > SelectQuadraturePointSetId
Definition: quadrature.hh:541
Definition: explicitfieldvector.hh:75
Definition: misc/functor.hh:31
static bool singleThreadMode()
returns true if program is operating on one thread currently
Definition: mpimanager.hh:436
static void unregisterStorage(StorageInterface &storage)
Definition: registry.hh:85
static void registerStorage(StorageInterface &storage)
Definition: registry.hh:70
wrapper for a (Quadrature,int) pair
Definition: quadrature.hh:43
unsigned int index() const
Definition: quadrature.hh:70
const QuadratureType & quadrature() const
Definition: quadrature.hh:65
size_t id() const
obtain the identifier of the integration point list
Definition: quadrature.hh:327
actual interface class for quadratures
Definition: quadrature.hh:405
A vector valued function space.
Definition: functionspace.hh:60
Definition: caching.hh:33
void evaluateEach(const QuadraturePointWrapper< Quadrature > &x, Functor functor) const
Definition: caching.hh:83
ShapeFunctionSet::HessianRangeType HessianRangeType
Definition: caching.hh:44
const ThisType & impl() const
Definition: caching.hh:125
void jacobianEach(const Point &x, Functor functor) const
Definition: caching.hh:90
const RangeVectorType & rangeCache(const QuadratureType &quadrature) const
Definition: caching.hh:111
std::vector< JacobianRangeType > JacobianRangeVectorType
Definition: caching.hh:47
void jacobianEach(const QuadraturePointWrapper< Quadrature > &x, Functor functor) const
Definition: caching.hh:96
ShapeFunctionSet ShapeFunctionSetType
Definition: caching.hh:37
std::vector< RangeType > RangeVectorType
Definition: caching.hh:46
std::vector< JacobianRangeVectorType > JacobianCacheVectorType
Definition: caching.hh:50
void hessianEach(const Point &x, Functor functor) const
Definition: caching.hh:103
ShapeFunctionSet::RangeType RangeType
Definition: caching.hh:42
CachingShapeFunctionSet(const GeometryType &type, const ShapeFunctionSet &shapeFunctionSet=ShapeFunctionSet())
Definition: caching.hh:56
ShapeFunctionSet::DomainType DomainType
Definition: caching.hh:41
GeometryType type() const
Definition: caching.hh:108
ShapeFunctionSet::JacobianRangeType JacobianRangeType
Definition: caching.hh:43
std::vector< RangeVectorType > RangeCacheVectorType
Definition: caching.hh:49
ShapeFunctionSet::FunctionSpaceType FunctionSpaceType
Definition: caching.hh:39
~CachingShapeFunctionSet()
Definition: caching.hh:254
void evaluateEach(const Point &x, Functor functor) const
Definition: caching.hh:77
std::size_t size() const
Definition: caching.hh:71
int order() const
Definition: caching.hh:66
const ThisType & scalarShapeFunctionSet() const
Definition: caching.hh:124
const JacobianRangeVectorType & jacobianCache(const QuadratureType &quadrature) const
Definition: caching.hh:118
static const int pointSetId
Definition: caching.hh:54
Interface class for shape function sets.
Definition: shapefunctionset/shapefunctionset.hh:33
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian range type
Definition: shapefunctionset/shapefunctionset.hh:43
FunctionSpaceType::RangeType RangeType
range type
Definition: shapefunctionset/shapefunctionset.hh:41
void hessianEach(const Point &x, Functor functor) const
evalute hessian of each shape function
void evaluateEach(const Point &x, Functor functor) const
evalute each shape function
std::size_t size() const
return number of shape functions
int order() const
return order of shape functions
FunctionSpaceType::DomainType DomainType
domain type
Definition: shapefunctionset/shapefunctionset.hh:39
void jacobianEach(const Point &x, Functor functor) const
evalute jacobian of each shape function