1#ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_CACHING_HH
2#define DUNE_FEM_SPACE_SHAPEFUNCTIONSET_CACHING_HH
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>
19#include <dune/fem/storage/dynamicarray.hh>
30 template<
class ShapeFunctionSet >
31 class CachingShapeFunctionSet
32 :
private QuadratureStorageRegistry::StorageInterface
34 typedef CachingShapeFunctionSet< ShapeFunctionSet > ThisType;
37 typedef ShapeFunctionSet ShapeFunctionSetType;
46 typedef std::vector< RangeType > RangeVectorType ;
47 typedef std::vector< JacobianRangeType > JacobianRangeVectorType ;
49 typedef std::vector< RangeVectorType > RangeCacheVectorType;
50 typedef std::vector< JacobianRangeVectorType > JacobianCacheVectorType;
54 static const int pointSetId = detail::SelectPointSetId< ShapeFunctionSetType >::value;
56 explicit CachingShapeFunctionSet (
const GeometryType &type,
57 const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() )
59 , shapeFunctionSet_( shapeFunctionSet )
61 QuadratureStorageRegistry::registerStorage( *
this );
64 ~CachingShapeFunctionSet ();
68 return shapeFunctionSet_.order();
71 std::size_t
size ()
const
73 return shapeFunctionSet_.size();
76 template<
class Po
int,
class Functor >
77 void evaluateEach (
const Point &x, Functor functor )
const
79 return shapeFunctionSet_.evaluateEach( x, functor );
82 template<
class Quadrature,
class Functor >
83 void evaluateEach (
const QuadraturePointWrapper< Quadrature > &x, Functor functor )
const
85 const bool cacheable = std::is_convertible< Quadrature, CachingInterface >::value;
86 evaluateEach( x.quadrature(), x.index(), functor, std::integral_constant< bool, cacheable >() );
89 template<
class Po
int,
class Functor >
90 void jacobianEach (
const Point &x, Functor functor )
const
92 return shapeFunctionSet_.jacobianEach( x, functor );
95 template<
class Quadrature,
class Functor >
96 void jacobianEach (
const QuadraturePointWrapper< Quadrature > &x, Functor functor )
const
98 const bool cacheable = std::is_convertible< Quadrature, CachingInterface >::value;
99 jacobianEach( x.quadrature(), x.index(), functor, std::integral_constant< bool, cacheable >() );
102 template<
class Po
int,
class Functor >
103 void hessianEach (
const Point &x, Functor functor )
const
105 return shapeFunctionSet_.hessianEach( x, functor );
110 template <
class QuadratureType >
111 const RangeVectorType& rangeCache(
const QuadratureType& quadrature )
const
113 return ReturnCache< QuadratureType, std::is_base_of< CachingInterface, QuadratureType >::value > ::
114 ranges( *
this, quadrature, rangeCaches_, localRangeCache_ );
117 template <
class QuadratureType >
118 const JacobianRangeVectorType& jacobianCache(
const QuadratureType& quadrature )
const
120 return ReturnCache< QuadratureType, std::is_base_of< CachingInterface, QuadratureType >::value > ::
121 jacobians( *
this, quadrature, jacobianCaches_, localJacobianCache_ );
124 const ThisType& scalarShapeFunctionSet()
const {
return *
this; }
125 const ThisType& impl()
const {
return *
this; }
128 template<
class Quad,
bool cacheable >
131 static const RangeVectorType&
132 ranges(
const ThisType& shapeFunctionSet,
134 const RangeCacheVectorType&,
135 RangeVectorType& storage )
143 const unsigned int nop = quad.nop();
144 const unsigned int size = shapeFunctionSet.size();
147 storage.resize(
size * nop );
148 RangeType* data = storage.data();
150 for(
unsigned int qp = 0 ; qp < nop; ++ qp )
152 const int cacheQp = quad.cachingPoint( qp );
153 AssignFunctor< RangeType* > funztor( data + (cacheQp *
size) );
154 shapeFunctionSet.evaluateEach( quad[ qp ], funztor );
159 static const JacobianRangeVectorType&
160 jacobians(
const ThisType& shapeFunctionSet,
162 const JacobianCacheVectorType&,
163 JacobianRangeVectorType& storage )
171 const unsigned int nop = quad.nop();
172 const unsigned int size = shapeFunctionSet.size();
175 storage.resize(
size * nop );
176 JacobianRangeType* data = storage.data();
178 for(
unsigned int qp = 0 ; qp < nop; ++ qp )
180 const int cacheQp = quad.cachingPoint( qp );
181 AssignFunctor< JacobianRangeType* > funztor( data + ( cacheQp *
size ) );
182 shapeFunctionSet.jacobianEach( quad[ qp ], funztor );
188 template<
class Quad >
189 struct ReturnCache< Quad, true >
191 static const RangeVectorType&
192 ranges(
const ThisType& shapeFunctionSet,
194 const RangeCacheVectorType& cache,
195 const RangeVectorType& )
197 return cache[ quad.id() ];
200 static const JacobianRangeVectorType&
201 jacobians(
const ThisType& shapeFunctionSet,
203 const JacobianCacheVectorType& cache,
204 const JacobianRangeVectorType& )
206 return cache[ quad.id() ];
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
215 evaluateEach( quadrature.point( pt ), functor );
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;
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
226 jacobianEach( quadrature.point( pt ), functor );
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;
234 void cacheQuadrature( std::size_t
id, std::size_t codim, std::size_t
size );
236 template<
class Po
intVector >
237 void cachePoints ( std::size_t
id,
const PointVector &points );
240 ShapeFunctionSet shapeFunctionSet_;
241 RangeCacheVectorType rangeCaches_;
242 JacobianCacheVectorType jacobianCaches_;
244 mutable RangeVectorType localRangeCache_ ;
245 mutable JacobianRangeVectorType localJacobianCache_;
253 template<
class ShapeFunctionSet >
254 inline CachingShapeFunctionSet< ShapeFunctionSet >::~CachingShapeFunctionSet ()
256 QuadratureStorageRegistry::unregisterStorage( *
this );
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
266 assert( (quadrature.id() < rangeCaches_.size()) && !rangeCaches_[ quadrature.id() ].empty() );
267 const RangeType *cache = rangeCaches_[ quadrature.id() ].data();
269 const unsigned int numShapeFunctions =
size();
270 const unsigned int cpt = quadrature.cachingPoint( pt );
278 if constexpr ( quadPointSetId == pointSetId )
281 if( quadrature.isInterpolationQuadrature(numShapeFunctions) )
285 assert( quadPointSetId >= 0 );
286 assert( pointSetId >= 0 );
291 const unsigned int point = quadrature.interpolationPoint( pt );
293 for(
unsigned int i = 0; i < numShapeFunctions; ++i )
294 assert( (cache[ cpt*numShapeFunctions + i ] - RangeType(i==point)).two_norm() < 1e-8 ) ;
297 functor( point, RangeType(1) );
302 for(
unsigned int i = 0; i < numShapeFunctions; ++i )
303 functor( i, cache[ cpt*numShapeFunctions + i ] );
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
313 assert( (quadrature.id() < jacobianCaches_.size()) && !jacobianCaches_[ quadrature.id() ].empty() );
314 const JacobianRangeType *cache = jacobianCaches_[ quadrature.id() ].data();
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 ] );
323 template<
class ShapeFunctionSet >
324 inline void CachingShapeFunctionSet< ShapeFunctionSet >
325 ::cacheQuadrature( std::size_t
id, std::size_t codim, std::size_t
size )
327 if( ! MPIManager::singleThreadMode () )
329 DUNE_THROW(SingleThreadModeError,
"CachingShapeFunctionSet::cacheQuadrature: only call in single thread mode!");
332 if(
id >= rangeCaches_.size() )
334 rangeCaches_.resize(
id+1, RangeVectorType() );
335 jacobianCaches_.resize(
id+1, JacobianRangeVectorType() );
338 assert( rangeCaches_[
id ].
empty() == jacobianCaches_[
id ].
empty() );
340 if( rangeCaches_[
id ].
empty() )
347 cachePoints(
id, PointProvider< ctype, dim, 0 >::getPoints(
id, type_ ) );
351 cachePoints(
id, PointProvider< ctype, dim, 1 >::getPoints(
id, type_ ) );
355 DUNE_THROW( NotImplemented,
"Caching for codim > 1 not implemented." );
361 template<
class ShapeFunctionSet >
362 template<
class Po
intVector >
363 inline void CachingShapeFunctionSet< ShapeFunctionSet >
364 ::cachePoints ( std::size_t
id,
const PointVector &points )
366 const unsigned int numShapeFunctions =
size();
367 const unsigned int numPoints = points.size();
369 RangeVectorType& ranges = rangeCaches_[ id ];
370 ranges.resize( numShapeFunctions * numPoints );
372 JacobianRangeVectorType& jacobians = jacobianCaches_[ id ];
373 jacobians.resize( numShapeFunctions * numPoints );
375 if( ranges.empty() || jacobians.empty() )
376 DUNE_THROW( OutOfMemoryError,
"Unable to allocate shape function set caches." );
378 for(
unsigned int pt = 0; pt < numPoints; ++pt )
380 evaluateEach( points[ pt ], AssignFunctor< RangeType * >( ranges.data() + pt*numShapeFunctions ) );
381 jacobianEach( points[ pt ], AssignFunctor< JacobianRangeType * >( jacobians.data() + pt*numShapeFunctions ) );
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