DUNE-FEM (unstable)

caching.hh
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
10
11// dune-fem includes
12#include <dune/fem/common/utility.hh>
13#include <dune/fem/misc/functor.hh>
14#include <dune/fem/misc/threads/threadsafevalue.hh>
15#include <dune/fem/quadrature/caching/registry.hh>
16#include <dune/fem/quadrature/cachingquadrature.hh>
17#include <dune/fem/quadrature/quadrature.hh>
18
19#include <dune/fem/storage/dynamicarray.hh>
20
21namespace Dune
22{
23
24 namespace Fem
25 {
26
27 // CachingShapeFunctionSet
28 // -----------------------
29
30 template< class ShapeFunctionSet >
31 class CachingShapeFunctionSet
32 : private QuadratureStorageRegistry::StorageInterface
33 {
34 typedef CachingShapeFunctionSet< ShapeFunctionSet > ThisType;
35
36 public:
37 typedef ShapeFunctionSet ShapeFunctionSetType;
38
40
41 typedef typename ShapeFunctionSet::DomainType DomainType;
42 typedef typename ShapeFunctionSet::RangeType RangeType;
43 typedef typename ShapeFunctionSet::JacobianRangeType JacobianRangeType;
44 typedef typename ShapeFunctionSet::HessianRangeType HessianRangeType;
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 {
61 QuadratureStorageRegistry::registerStorage( *this );
62 }
63
64 ~CachingShapeFunctionSet ();
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,
134 const RangeCacheVectorType&,
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,
162 const JacobianCacheVectorType&,
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,
204 const JacobianRangeVectorType& )
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 >
254 inline CachingShapeFunctionSet< ShapeFunctionSet >::~CachingShapeFunctionSet ()
255 {
256 QuadratureStorageRegistry::unregisterStorage( *this );
257 }
258
259
260 template< class ShapeFunctionSet >
261 template< class Quadrature, class Functor >
262 inline void CachingShapeFunctionSet< ShapeFunctionSet >
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 {
327 if( ! MPIManager::singleThreadMode () )
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
FunctionSpaceTraits::DomainFieldType DomainFieldType
Intrinsic type used for values in the domain field (usually a double)
Definition: functionspaceinterface.hh:60
@ dimDomain
dimension of domain vector space
Definition: functionspaceinterface.hh:46
A vector valued function space.
Definition: functionspace.hh:60
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian range type
Definition: shapefunctionset.hh:43
FunctionSpace FunctionSpaceType
function space type
Definition: shapefunctionset.hh:36
FunctionSpaceType::HessianRangeType HessianRangeType
hessian range type
Definition: shapefunctionset.hh:45
FunctionSpaceType::RangeType RangeType
range type
Definition: shapefunctionset.hh:41
FunctionSpaceType::DomainType DomainType
domain type
Definition: shapefunctionset.hh:39
actual interface class for quadratures
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::bool_constant<(sizeof...(II)==0)> empty(std::integer_sequence< T, II... >)
Checks whether the sequence is empty.
Definition: integersequence.hh:80
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)