dune-fem  2.8-git
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
13 #include <dune/fem/misc/functor.hh>
18 
20 
21 namespace 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,
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,
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 >
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 >
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( id >= rangeCaches_.size() )
328  {
329  rangeCaches_.resize( id+1, RangeVectorType() );
330  jacobianCaches_.resize( id+1, JacobianRangeVectorType() );
331  }
332 
333  assert( rangeCaches_[ id ].empty() == jacobianCaches_[ id ].empty() );
334 
335  if( rangeCaches_[ id ].empty() )
336  {
337  typedef typename FunctionSpaceType::DomainFieldType ctype;
338  const int dim = FunctionSpaceType::dimDomain;
339  switch( codim )
340  {
341  case 0:
342  cachePoints( id, PointProvider< ctype, dim, 0 >::getPoints( id, type_ ) );
343  break;
344 
345  case 1:
346  cachePoints( id, PointProvider< ctype, dim, 1 >::getPoints( id, type_ ) );
347  break;
348 
349  default:
350  DUNE_THROW( NotImplemented, "Caching for codim > 1 not implemented." );
351  }
352  }
353  }
354 
355 
356  template< class ShapeFunctionSet >
357  template< class PointVector >
358  inline void CachingShapeFunctionSet< ShapeFunctionSet >
359  ::cachePoints ( std::size_t id, const PointVector &points )
360  {
361  const unsigned int numShapeFunctions = size();
362  const unsigned int numPoints = points.size();
363 
364  RangeVectorType& ranges = rangeCaches_[ id ];
365  ranges.resize( numShapeFunctions * numPoints );
366 
367  JacobianRangeVectorType& jacobians = jacobianCaches_[ id ];
368  jacobians.resize( numShapeFunctions * numPoints );
369 
370  if( ranges.empty() || jacobians.empty() )
371  DUNE_THROW( OutOfMemoryError, "Unable to allocate shape function set caches." );
372 
373  for( unsigned int pt = 0; pt < numPoints; ++pt )
374  {
375  evaluateEach( points[ pt ], AssignFunctor< RangeType * >( ranges.data() + pt*numShapeFunctions ) );
376  jacobianEach( points[ pt ], AssignFunctor< JacobianRangeType * >( jacobians.data() + pt*numShapeFunctions ) );
377  }
378  }
379 
380  } // namespace Fem
381 
382 } // namespace Dune
383 
384 #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 void unregisterStorage(StorageInterface &storage)
Definition: registry.hh:82
static void registerStorage(StorageInterface &storage)
Definition: registry.hh:68
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
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
const JacobianRangeVectorType & jacobianCache(const QuadratureType &quadrature) const
Definition: caching.hh:118
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
const RangeVectorType & rangeCache(const QuadratureType &quadrature) const
Definition: caching.hh:111
std::size_t size() const
Definition: caching.hh:71
int order() const
Definition: caching.hh:66
static const int pointSetId
Definition: caching.hh:54
const ThisType & scalarShapeFunctionSet() const
Definition: caching.hh:124
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