dune-fem  2.8-git
tensorproduct.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SHAPEFUNCTIONSET_TENSORPRODUCT_HH
2 #define DUNE_FEM_SHAPEFUNCTIONSET_TENSORPRODUCT_HH
3 
4 #include <algorithm>
5 #include <array>
6 #include <cstddef>
7 #include <tuple>
8 #include <utility>
9 
10 #include <dune/common/fmatrix.hh>
12 #include <dune/common/fvector.hh>
13 #include <dune/common/hybridutilities.hh>
14 #include <dune/common/tupleutility.hh>
16 
17 namespace Dune
18 {
19 
20  namespace Fem
21  {
22 
23  // TensorProductShapeFunctionSet
24  // -----------------------------
25 
26  template< class FunctionSpace, class ShapeFunctionSetTuple >
28  {
30 
31  static_assert( (FunctionSpace::dimDomain == std::tuple_size< ShapeFunctionSetTuple >::value),
32  "dimDomain of FunctionSpace must coincide with length of ShapeFunctionSetTuple." );
33  static_assert( (FunctionSpace::dimRange == 1),
34  "FunctionSpace must be scalar (i.e., dimRange = 1)." );
35 
36  struct Assign;
37 
38  template< int i > struct Size;
39  template< int i > struct EvaluateAll;
40  template< int i > struct JacobianAll;
41  template< int i > struct HessianAll;
42 
43  static const int dimension = FunctionSpace::dimDomain;
44 
45  public:
47  typedef ShapeFunctionSetTuple ShapeFunctionSetTupleType;
48 
51 
56 
58 
59  explicit TensorProductShapeFunctionSet ( const ShapeFunctionSetTupleType &shapeFunctionSetTuple );
61 
62  TensorProductShapeFunctionSet ( const ThisType &other );
63  const ThisType &operator= ( const ThisType &other );
64 
65  int order () const;
66 
67  std::size_t size () const;
68 
69  template< class Point, class Functor >
70  void evaluateEach ( const Point &x, Functor functor ) const;
71 
72  template< class Point, class Functor >
73  void jacobianEach ( const Point &x, Functor functor ) const;
74 
75  template< class Point, class Functor >
76  void hessianEach ( const Point &x, Functor functor ) const;
77 
78  private:
79  template< class Functor >
80  void doEvaluateEach ( int d, RangeType value, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const;
81  template< class Functor >
82  void doJacobianEach ( int d, JacobianRangeType jacobian, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const;
83  template< class Functor >
84  void doHessianEach ( int d, HessianRangeType hessian, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const;
85 
86  ShapeFunctionSetTuple shapeFunctionSetTuple_;
87  std::array< std::size_t, dimension > sizes_;
88  RangeFieldType *buffer_;
89  };
90 
91 
92 
93  // TensorProductShapeFunctionSet::Assign
94  // -------------------------------------
95 
96  template< class FunctionSpace, class ShapeFunctionSetTuple >
97  struct TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >::Assign
98  {
99  explicit Assign ( RangeFieldType *buffer ) : buffer_( buffer ) {}
100 
101  void operator() ( const std::size_t i, const RangeFieldType &value )
102  {
103  buffer_[ i ] = value;
104  }
105 
106  template< class T >
107  void operator() ( const std::size_t i, const FieldVector< T, 1 > &value )
108  {
109  (*this)( i, value[ 0 ] );
110  }
111 
112  template< class T >
113  void operator() ( const std::size_t i, const FieldMatrix< T, 1, 1 > &value )
114  {
115  (*this)( i, value[ 0 ][ 0 ] );
116  }
117 
118  private:
119  RangeFieldType *buffer_;
120  };
121 
122 
123 
124  // TensorProductShapeFunctionSet::Size
125  // -----------------------------------
126 
127  template< class FunctionSpace, class ShapeFunctionSetTuple >
128  template< int i >
129  struct TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >::Size
130  {
131  static void apply ( const ShapeFunctionSetTuple &tuple, std::array< std::size_t, FunctionSpace::dimDomain > &size )
132  {
133  size[ i ] = std::get< i >( tuple ).size();
134  }
135  };
136 
137 
138 
139  // TensorProductShapeFunctionSet::EvaluateAll
140  // ------------------------------------------
141 
142  template< class FunctionSpace, class ShapeFunctionSetTuple >
143  template< int i >
144  struct TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >::EvaluateAll
145  {
146  static void apply ( const ShapeFunctionSetTuple &tuple, const DomainType &x, RangeFieldType *&it )
147  {
148  Dune::FieldVector< DomainFieldType, 1 > xi( x[ i ] );
149  std::get< i >( tuple ).evaluateEach( xi, Assign( it ) );
150  it += std::get< i >( tuple ).size();
151  }
152  };
153 
154 
155 
156  // TensorProductShapeFunctionSet::JacobianAll
157  // ------------------------------------------
158 
159  template< class FunctionSpace, class ShapeFunctionSetTuple >
160  template< int i >
161  struct TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >::JacobianAll
162  {
163  static void apply ( const ShapeFunctionSetTuple &tuple, const DomainType &x, RangeFieldType *&it )
164  {
165  Dune::FieldVector< DomainFieldType, 1 > xi( x[ i ] );
166  const std::size_t size = std::get< i >( tuple ).size();
167  std::get< i >( tuple ).evaluateEach( xi, Assign( it ) );
168  std::get< i >( tuple ).jacobianEach( xi, Assign( it+size ) );
169  it += 2*size;
170  }
171  };
172 
173 
174 
175  // TensorProductShapeFunctionSet::HessianAll
176  // -----------------------------------------
177 
178  template< class FunctionSpace, class ShapeFunctionSetTuple >
179  template< int i >
180  struct TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >::HessianAll
181  {
182  static void apply ( const ShapeFunctionSetTuple &tuple, const DomainType &x, RangeFieldType *&it )
183  {
184  Dune::FieldVector< DomainFieldType, 1 > xi( x[ i ] );
185  const std::size_t size = std::get< i >( tuple ).size();
186  std::get< i >( tuple ).evaluateEach( xi, Assign( it ) );
187  std::get< i >( tuple ).jacobianEach( xi, Assign( it+size ) );
188  std::get< i >( tuple ).hessianEach( xi, Assign( it+2*size ) );
189  it += 3*size;
190  }
191  };
192 
193 
194 
195  // Implementation of TensorProductShapeFunctionSet
196  // -----------------------------------------------
197 
198  template< class FunctionSpace, class ShapeFunctionSetTuple >
201  : shapeFunctionSetTuple_( shapeFunctionSetTuple )
202  {
203  Fem::ForLoop< Size, 0, dimension-1 >::apply( shapeFunctionSetTuple_, sizes_ );
204  std::size_t buffer_size = 0;
205  for( int i = 0; i < dimension; ++i )
206  buffer_size += sizes_[ i ];
207  buffer_ = new RangeFieldType[ 3*buffer_size ];
208  }
209 
210 
211  template< class FunctionSpace, class ShapeFunctionSetTuple >
214  {
215  delete[]( buffer_ );
216  }
217 
218 
219  template< class FunctionSpace, class ShapeFunctionSetTuple >
222  : shapeFunctionSetTuple_( other.shapeFunctionSetTuple_ )
223  {
224  std::size_t buffer_size = 0;
225  for( int i = 0; i < dimension; ++i )
226  {
227  sizes_[ i ] = other.sizes_[ i ];
228  buffer_size += sizes_[ i ];
229  }
230  buffer_ = new RangeFieldType[ 3*buffer_size ];
231  }
232 
233 
234  template< class FunctionSpace, class ShapeFunctionSetTuple >
237  ::operator= ( const ThisType &other )
238  {
239  if( this == &other )
240  return *this;
241  delete[]( buffer_ );
242 
243  shapeFunctionSetTuple_ = other.shapeFunctionSetTuple_;
244  std::size_t buffer_size = 0;
245  for( int i = 0; i < dimension; ++i )
246  {
247  sizes_[ i ] = other.sizes_[ i ];
248  buffer_size += sizes_[ i ];
249  }
250  buffer_ = new RangeFieldType[ 3*buffer_size ];
251  return *this;
252  }
253 
254 
255  template< class FunctionSpace, class ShapeFunctionSetTuple >
257  {
258  int value( 0 );
259  Hybrid::forEach( std::make_index_sequence< std::tuple_size< ShapeFunctionSetTupleType >::value >{},
260  [ & ]( auto i ){ value = std::max( value, std::get< i >( shapeFunctionSetTuple_ ).order() ); } );
261  return value;
262  }
263 
264 
265  template< class FunctionSpace, class ShapeFunctionSetTuple >
266  inline std::size_t
268  {
269  std::size_t size( 1 );
270  for( int i = 0; i < dimension; ++i )
271  size *= sizes_[ i ];
272  return size;
273  }
274 
275 
276  template< class FunctionSpace, class ShapeFunctionSetTuple >
277  template< class Point, class Functor >
279  ::evaluateEach ( const Point &x, Functor functor ) const
280  {
281  RangeFieldType *it = buffer_;
282  Fem::ForLoop< EvaluateAll, 0, dimension-1 >::apply( shapeFunctionSetTuple_, coordinate( x ), it );
283 
284  std::size_t index = 0;
285  doEvaluateEach( 0, RangeType( RangeFieldType( 1 ) ), index, buffer_, functor );
286  }
287 
288 
289  template< class FunctionSpace, class ShapeFunctionSetTuple >
290  template< class Point, class Functor >
292  ::jacobianEach ( const Point &x, Functor functor ) const
293  {
294  RangeFieldType *it = buffer_;
295  Fem::ForLoop< JacobianAll, 0, dimension-1 >::apply( shapeFunctionSetTuple_, coordinate( x ), it );
296 
297  std::size_t index = 0;
298  JacobianRangeType jacobian;
299  for( int i = 0; i < dimension; ++i )
300  jacobian[ 0 ][ i ] = RangeFieldType( 1 );
301  doJacobianEach( 0, jacobian, index, buffer_, functor );
302  }
303 
304 
305  template< class FunctionSpace, class ShapeFunctionSetTuple >
306  template< class Point, class Functor >
308  ::hessianEach ( const Point &x, Functor functor ) const
309  {
310  RangeFieldType *it = buffer_;
311  Fem::ForLoop< HessianAll, 0, dimension-1 >::apply( shapeFunctionSetTuple_, coordinate( x ), it );
312 
313  std::size_t index = 0;
314  HessianRangeType hessian;
315  for( int i = 0; i < dimension; ++i )
316  for( int j = 0; j < dimension; ++j )
317  hessian[ 0 ][ i ][ j ] = RangeFieldType( 1 );
318  doHessianEach( 0, hessian, index, buffer_, functor );
319  }
320 
321 
322  template< class FunctionSpace, class ShapeFunctionSetTuple >
323  template< class Functor >
325  ::doEvaluateEach ( int d, RangeType value, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const
326  {
327  if( d < dimension )
328  {
329  for( std::size_t i = 0; i < sizes_[ d ]; ++i )
330  {
331  RangeType v( value );
332  v[ 0 ] *= buffer[ i ];
333  doEvaluateEach( d+1, v, index, buffer+sizes_[ d ], functor );
334  }
335  }
336  else
337  functor( index++, value );
338  }
339 
340 
341  template< class FunctionSpace, class ShapeFunctionSetTuple >
342  template< class Functor >
343  inline void TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >
344  ::doJacobianEach ( int d, JacobianRangeType jacobian, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const
345  {
346  if( d < dimension )
347  {
348  for( std::size_t i = 0; i < sizes_[ d ]; ++i )
349  {
350  JacobianRangeType j( jacobian );
351  j[ 0 ][ d ] *= buffer[ i + sizes_[ d ] ];
352  for( int k = 1; k < dimension; ++k )
353  j[ 0 ][ (d+k)%dimension ] *= buffer[ i ];
354  doJacobianEach( d+1, j, index, buffer+2*sizes_[ d ], functor );
355  }
356  }
357  else
358  functor( index++, jacobian );
359  }
360 
361 
362  template< class FunctionSpace, class ShapeFunctionSetTuple >
363  template< class Functor >
364  inline void TensorProductShapeFunctionSet< FunctionSpace, ShapeFunctionSetTuple >
365  ::doHessianEach ( int d, HessianRangeType hessian, std::size_t &index, const RangeFieldType *buffer, Functor functor ) const
366  {
367  if( d < dimension )
368  {
369  for( std::size_t i = 0; i < sizes_[ d ]; ++i )
370  {
371  HessianRangeType h( hessian );
372  h[ 0 ][ d ][ d ] *= buffer[ i + 2*sizes_[ d ] ];
373  for( int j = 1; j < dimension; ++j )
374  {
375  h[ 0 ][ (d+j)%dimension ][ d ] *= buffer[ i * sizes_[ d ] ];
376  h[ 0 ][ d ][ (d+j)%dimension ] *= buffer[ i * sizes_[ d ] ];
377  for( int k = 1; k < dimension; ++k )
378  h[ 0 ][ (d+j)%dimension ][ (d+k)%dimension ] *= buffer[ i ];
379  }
380  doHessianEach( d+1, h, index, buffer+3*sizes_[ d ], functor );
381  }
382  }
383  else
384  functor( index++, hessian );
385  }
386 
387  } // namespace Fem
388 
389 } // namespace Dune
390 
391 #endif // #ifndef DUNE_FEM_SHAPEFUNCTIONSET_TENSORPRODUCT_HH
Definition: bindguard.hh:11
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static constexpr T max(T a)
Definition: utility.hh:77
static DUNE_PRIVATE void apply(Args &&... args)
Definition: forloop.hh:23
A vector valued function space.
Definition: functionspace.hh:60
ExplicitFieldVector< FieldMatrix< RangeFieldType, dimDomain, dimDomain >, dimRange > HessianRangeType
Intrinsic type used for the hessian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:79
FunctionSpaceTraits::DomainFieldType DomainFieldType
Intrinsic type used for values in the domain field (usually a double)
Definition: functionspaceinterface.hh:60
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
@ dimDomain
dimension of domain vector space
Definition: functionspaceinterface.hh:46
@ dimRange
dimension of range vector space
Definition: functionspaceinterface.hh:48
FunctionSpaceTraits::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
FunctionSpaceTraits::DomainType DomainType
Type of domain vector (using type of domain field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:67
FunctionSpaceTraits::RangeFieldType RangeFieldType
Intrinsic type used for values in the range field (usually a double)
Definition: functionspaceinterface.hh:63
Definition: tensorproduct.hh:28
void evaluateEach(const Point &x, Functor functor) const
Definition: tensorproduct.hh:279
FunctionSpaceType::RangeType RangeType
Definition: tensorproduct.hh:53
void jacobianEach(const Point &x, Functor functor) const
Definition: tensorproduct.hh:292
FunctionSpaceType::HessianRangeType HessianRangeType
Definition: tensorproduct.hh:55
FunctionSpaceType::DomainType DomainType
Definition: tensorproduct.hh:52
void hessianEach(const Point &x, Functor functor) const
Definition: tensorproduct.hh:308
int order() const
Definition: tensorproduct.hh:256
std::size_t size() const
Definition: tensorproduct.hh:267
FunctionSpaceType::JacobianRangeType JacobianRangeType
Definition: tensorproduct.hh:54
const ThisType & operator=(const ThisType &other)
Definition: tensorproduct.hh:237
FunctionSpace FunctionSpaceType
Definition: tensorproduct.hh:46
FunctionSpaceType::RangeFieldType RangeFieldType
Definition: tensorproduct.hh:50
FunctionSpaceType::DomainFieldType DomainFieldType
Definition: tensorproduct.hh:49
~TensorProductShapeFunctionSet()
Definition: tensorproduct.hh:213
ShapeFunctionSetTuple ShapeFunctionSetTupleType
Definition: tensorproduct.hh:47
Assign(RangeFieldType *buffer)
Definition: tensorproduct.hh:99