dune-fem  2.8-git
space/fourier/functionset.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_FOURIER_FUNCTIONSET_HH
2 #define DUNE_FEM_SPACE_FOURIER_FUNCTIONSET_HH
3 
4 #include <cassert>
5 #include <cstddef>
6 #include <limits>
7 #include <array>
8 
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/power.hh>
12 
14 #include <dune/fem/version.hh>
15 
16 
17 namespace Dune
18 {
19 
20  namespace Fem
21  {
22 
23  // FourierFunctionSetSize
24  // ----------------------
25 
26  template< int dimension, int Order >
28  {
29  static const int v = StaticPower< (2*Order+1), dimension >::power;
30  };
31 
32 
33 
34  // FourierFunctionSet
35  // ------------------
36 
37  template< class FunctionSpace, int Order >
39 
40 
41 
42  // Template specialization for dimDomain = dimRange = 1
43  // ----------------------------------------------------
44 
45  template< class DomainFieldType, class RangeFieldType, int Order >
46  class FourierFunctionSet< FunctionSpace< DomainFieldType, RangeFieldType, 1, 1 >, Order >
47  {
49 
50  public:
52 
57 
58  typedef std::size_t SizeType;
59 
60  explicit FourierFunctionSet ( int order ) : order_( order ) {}
61 
62  int order () const { return order_; }
63 
65 
66  template< class Functor >
67  static void evaluateEach ( const DomainType &x, Functor functor )
68  {
69  using std::sin;
70  using std::cos;
71  functor( 0, RangeType( RangeFieldType( 1 ) / RangeFieldType( 2 ) ) );
72  // use recursion:
73  // sin((n+1)*x) = sin(n*x)*cos(x) + cos(n*x)*sin(x)
74  // cos((n+1)*x) = cos(n*x)*cos(x) - sin(n*x)*sin(x)
75  SizeType basisFunction = 1;
76  for( int n = 1; n <= Order; ++n )
77  {
78  functor( basisFunction++, RangeType( cos( n*x[ 0 ] ) ) );
79  functor( basisFunction++, RangeType( sin( n*x[ 0 ] ) ) );
80  }
81  }
82 
83  template< class Functor >
84  static void jacobianEach ( const DomainType &x, Functor functor )
85  {
86  using std::sin;
87  using std::cos;
88  functor( 0, JacobianRangeType( RangeFieldType( 0 ) ) );
89  SizeType basisFunction = 1;
90  for( int n = 1; n <= Order; ++n )
91  {
92  functor( basisFunction++, JacobianRangeType( -n*sin( n*x[ 0 ] ) ) );
93  functor( basisFunction++, JacobianRangeType( n*cos( n*x[ 0 ] ) ) );
94  }
95  }
96 
97  template< class Functor >
98  static void hessianEach ( const DomainType &x, Functor functor )
99  {
100  using std::sin;
101  using std::cos;
102  functor( 0, HessianRangeType( RangeFieldType( 0 ) ) );
103  SizeType basisFunction = 1;
104  for( int n = 1; n <= Order; ++n )
105  {
106  functor( basisFunction++, HessianRangeType( -(n*n)*cos( n*x[ 0 ] ) ) );
107  functor( basisFunction++, HessianRangeType( -(n*n)*sin( n*x[ 0 ] ) ) );
108  }
109  }
110  private:
111  int order_;
112  };
113 
114 
115 
116  // Template specialization for dimDomain > 1, dimRange = 1
117  // -------------------------------------------------------
118 
119  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
120  class FourierFunctionSet< FunctionSpace< DomainFieldType, RangeFieldType, dimDomain, 1 >, Order >
121  {
123 
124  public:
126 
127  typedef typename FunctionSpaceType::DomainType DomainType;
131 
132  typedef std::size_t SizeType;
133 
134  private:
135  // number of Fourier basis function for dimDomain = 1
136  static const int buffer_size = FourierFunctionSetSize< 1, Order >::v;
137 
138  // tags used for building cache
139  struct Evaluate {}; //< evaluate basis functions
140  struct Jacobian {}; //< evaluate basis functions and jacobians
141  struct Hessian {}; //< evaluate basis functions, jacobians, and hessians
142 
143  struct Assign;
144 
145  protected:
146  // multi index used for tensor product basis functions
147  typedef Dune::FieldVector< int, dimDomain > MultiIndexType;
148 
149  // iterator type and methods for accessing multi indices
150  struct MultiIndexIterator;
151  typedef MultiIndexIterator IteratorType;
152  static IteratorType begin () { return IteratorType::begin(); }
153  static IteratorType end () { return IteratorType::end(); }
154 
155  public:
156  explicit FourierFunctionSet ( int order ) : order_( order ) {}
157 
158  int order () const { return order_; }
159 
161 
162  template< class Functor >
163  void evaluateEach ( const DomainType &x, Functor functor ) const
164  {
165  prepare( Evaluate(), x );
166  SizeType index( 0 );
167  const IteratorType end = ThisType::end();
168  for( IteratorType it = ThisType::begin(); it != end; ++it, ++index )
169  {
170  RangeType value;
171  evaluate( *it, value );
172  assert( index == IteratorType::index( *it ) );
173  functor( index, value );
174  }
175  assert( index == size() );
176  }
177 
178  template< class Functor >
179  void jacobianEach ( const DomainType &x, Functor functor ) const
180  {
181  prepare( Jacobian(), x );
182  SizeType index( 0 );
183  const IteratorType end = ThisType::end();
184  for( IteratorType it = ThisType::begin(); it != end; ++it, ++index )
185  {
186  JacobianRangeType jacobian;
187  evaluate( *it, jacobian );
188  assert( index == IteratorType::index( *it ) );
189  functor( index, jacobian );
190  }
191  assert( index == size() );
192  }
193 
194  template< class Functor >
195  void hessianEach ( const DomainType &x, Functor functor ) const
196  {
197  prepare( Hessian(), x );
198  SizeType index( 0 );
199  const IteratorType end = ThisType::end();
200  for( IteratorType it = ThisType::begin(); it != end; ++it, ++index )
201  {
202  HessianRangeType hessian;
203  evaluate( *it, hessian );
204  assert( index == IteratorType::index( *it ) );
205  functor( index, hessian );
206  }
207  assert( index == size() );
208  }
209 
210  protected:
211  // evaluate tensor product basis function
212  void evaluate ( const MultiIndexType &multiIndex, RangeType &value ) const
213  {
214  value = RangeType( RangeFieldType( 1 ) );
215  for( SizeType i = 0; i < dimDomain; ++i )
216  value *= buffer_[ i ][ multiIndex[ i ] ];
217  }
218 
219  // evaluate jacobian of tensor product basis function
220  void evaluate ( const MultiIndexType &multiIndex, JacobianRangeType &jacobian ) const
221  {
222  jacobian = JacobianRangeType( 1 );
223  for( int k = 0; k < dimDomain; ++k )
224  {
225  const RangeFieldType phi = buffer_[ k ][ multiIndex[ k ] ];
226  const RangeFieldType dphi = buffer_[ k ][ buffer_size + multiIndex[ k ] ];
227  for( int i = 0; i < dimDomain; ++i )
228  jacobian[ 0 ][ i ] *= (k == i ? dphi : phi);
229  }
230  }
231 
232  // evaluate hessian of tensor product basis function
233  void evaluate ( const MultiIndexType &multiIndex, HessianRangeType &hessian ) const
234  {
235  for( int i = 0; i < dimDomain; ++i )
236  for( int j = 0; j < dimDomain; ++j )
237  hessian[ 0 ][ i ][ j ] = RangeFieldType( 1 );
238 
239  for( int k = 0; k < dimDomain; ++k )
240  {
241  const RangeFieldType phi = buffer_[ k ][ multiIndex[ k ] ];
242  const RangeFieldType dphi = buffer_[ k ][ buffer_size + multiIndex[ k ] ];
243  for( int i = 0; i < dimDomain; ++i )
244  {
245  hessian[ 0 ][ i ][ i ] *= (k == i ? buffer_[ i ][ 2*buffer_size + multiIndex[ i ] ] : phi);
246  for( int j = i+1; j < dimDomain; ++j )
247  {
248  RangeFieldType tmp = ( k == i || k == j ) ? dphi : phi;
249  hessian[ 0 ][ i ][ j ] *= tmp;
250  hessian[ 0 ][ j ][ i ] *= tmp;
251  }
252  }
253  }
254  }
255 
256  // methods for building cache
257  void prepare ( const Evaluate &, const DomainType &x ) const;
258  void prepare ( const Jacobian &, const DomainType &x ) const;
259  void prepare ( const Hessian &, const DomainType &x ) const;
260 
261  private:
262  int order_;
263  // cache for evaluation of basis functions
264  mutable std::array< std::array< RangeFieldType, 3*buffer_size>, dimDomain > buffer_;
265  };
266 
267 
268 
269  // Implementation of FourierFunctionSet::Assign
270  // --------------------------------------------
271 
272  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
273  struct FourierFunctionSet< FunctionSpace< DomainFieldType, RangeFieldType, dimDomain, 1 >, Order >::Assign
274  {
275  explicit Assign ( RangeFieldType *buffer ) : buffer_( buffer ) {}
276 
277  void operator() ( const std::size_t i, const RangeFieldType &value )
278  {
279  buffer_[ i ] = value;
280  }
281 
282  template< class T >
283  void operator() ( const std::size_t i, const FieldVector< T, 1 > &value )
284  {
285  (*this)( i, value[ 0 ] );
286  }
287 
288  template< class T >
289  void operator() ( const std::size_t i, const FieldMatrix< T, 1, 1 > &value )
290  {
291  (*this)( i, value[ 0 ][ 0 ] );
292  }
293 
294  private:
295  RangeFieldType *buffer_;
296  };
297 
298 
299 
300  // Implementation of FourierFunctionSet::MultiIndexIterator
301  // --------------------------------------------------------
302 
303  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
304  struct FourierFunctionSet< FunctionSpace< DomainFieldType, RangeFieldType, dimDomain, 1 >, Order >::MultiIndexIterator
305  {
306  typedef MultiIndexIterator ThisType;
307 
308  protected:
309  typedef int IndexType;
310 
311  explicit MultiIndexIterator ( IndexType n ) : multiIndex_( n ) {}
312 
314 
316 
317  public:
318  static ThisType begin () { return ThisType( 0 ); }
319  static ThisType end () { return ThisType( invalidIndex() ); }
320 
321  ThisType operator++ ()
322  {
323  // try to increment and leave eventually
324  for( int i = 0; i < dimDomain; ++i )
325  {
326  const int j = dimDomain-i-1;
327  if( ++multiIndex_[ j ] < N )
328  return *this;
329  multiIndex_[ j ] = 0;
330  }
331 
332  // otherwise, reset this iterator to end iterator
333  *this = end();
334  return *this;
335  }
336 
337  bool operator== ( const ThisType &other ) const { return ( multiIndex_ == other.multiIndex_ ); }
338 
339  bool operator!= ( const ThisType &other ) const { return !( *this == other ); }
340 
341  const MultiIndexType &operator* () const { return multiIndex_; }
342 
343  const MultiIndexType *operator-> () const { return &multiIndex_; }
344 
345  SizeType index () const { return index( multiIndex_ ); }
346 
347  static SizeType index ( const MultiIndexType &multiIndex )
348  {
349  SizeType index = 0, factor = 1;
350  for( int i = dimDomain-1; i >= 0; --i )
351  {
352  index += multiIndex[ i ]*factor;
353  factor *= N;
354  }
355  assert( index < size() );
356  return index;
357  }
358 
359  private:
360  MultiIndexType multiIndex_;
361  };
362 
363 
364 
365  // Implementation of FourierFunctionSet
366  // ------------------------------------
367 
368  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
370  ::prepare ( const Evaluate &, const DomainType &x ) const
371  {
373  for( SizeType i = 0; i < dimDomain; ++i )
374  {
375  RangeFieldType *it = &buffer_[ i ][ 0 ];
376  Dune::FieldVector< DomainFieldType, 1 > y( x[ i ] );
377  FunctionSetImp::evaluateEach( y, Assign( it ) );
378  }
379  };
380 
381 
382  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
384  ::prepare ( const Jacobian &, const DomainType &x ) const
385  {
387  for( SizeType i = 0; i < dimDomain; ++i )
388  {
389  RangeFieldType *it = &buffer_[ i ][ 0 ];
390  Dune::FieldVector< DomainFieldType, 1 > y( x[ i ] );
391  FunctionSetImp::evaluateEach( y, Assign( it ) );
392  FunctionSetImp::jacobianEach( y, Assign( it+buffer_size ) );
393  }
394  };
395 
396 
397  template< class DomainFieldType, class RangeFieldType, int dimDomain, int Order >
399  ::prepare ( const Hessian &, const DomainType &x ) const
400  {
402  for( SizeType i = 0; i < dimDomain; ++i )
403  {
404  RangeFieldType *it = &buffer_[ i ][ 0 ];
405  Dune::FieldVector< DomainFieldType, 1 > y( x[ i ] );
406  FunctionSetImp::evaluateEach( y, Assign( it ) );
407  FunctionSetImp::jacobianEach( y, Assign( it+buffer_size) );
408  FunctionSetImp::hessianEach( y, Assign( it+2*buffer_size ) );
409  }
410  };
411 
412  } // namespace Fem
413 
414 } // namespace Dune
415 
416 #endif // #ifndef DUNE_FEM_SPACE_FOURIER_FUNCTIONSET_HH
Definition: bindguard.hh:11
Double operator*(const Double &a, const Double &b)
Definition: double.hh:506
static double cos(const Double &v)
Definition: double.hh:896
static double sin(const Double &v)
Definition: double.hh:891
bool operator==(const Double &a, const Double &b)
Definition: double.hh:600
bool operator!=(const Double &a, const Double &b)
Definition: double.hh:640
static constexpr T max(T a)
Definition: utility.hh:77
Definition: explicitfieldvector.hh:75
A vector valued function space.
Definition: functionspace.hh:60
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
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
Definition: space/fourier/functionset.hh:28
static const int v
Definition: space/fourier/functionset.hh:29
Definition: space/fourier/functionset.hh:38
FunctionSpaceType::RangeType RangeType
Definition: space/fourier/functionset.hh:54
static void jacobianEach(const DomainType &x, Functor functor)
Definition: space/fourier/functionset.hh:84
FunctionSpace< DomainFieldType, RangeFieldType, 1, 1 > FunctionSpaceType
Definition: space/fourier/functionset.hh:51
static void hessianEach(const DomainType &x, Functor functor)
Definition: space/fourier/functionset.hh:98
static void evaluateEach(const DomainType &x, Functor functor)
Definition: space/fourier/functionset.hh:67
FunctionSpaceType::HessianRangeType HessianRangeType
Definition: space/fourier/functionset.hh:56
FunctionSpaceType::JacobianRangeType JacobianRangeType
Definition: space/fourier/functionset.hh:55
FunctionSpaceType::DomainType DomainType
Definition: space/fourier/functionset.hh:53
void evaluate(const MultiIndexType &multiIndex, HessianRangeType &hessian) const
Definition: space/fourier/functionset.hh:233
void evaluateEach(const DomainType &x, Functor functor) const
Definition: space/fourier/functionset.hh:163
void evaluate(const MultiIndexType &multiIndex, RangeType &value) const
Definition: space/fourier/functionset.hh:212
FunctionSpaceType::HessianRangeType HessianRangeType
Definition: space/fourier/functionset.hh:130
void jacobianEach(const DomainType &x, Functor functor) const
Definition: space/fourier/functionset.hh:179
void hessianEach(const DomainType &x, Functor functor) const
Definition: space/fourier/functionset.hh:195
void evaluate(const MultiIndexType &multiIndex, JacobianRangeType &jacobian) const
Definition: space/fourier/functionset.hh:220
FunctionSpaceType::DomainType DomainType
Definition: space/fourier/functionset.hh:127
Dune::FieldVector< int, dimDomain > MultiIndexType
Definition: space/fourier/functionset.hh:143
FunctionSpace< DomainFieldType, RangeFieldType, dimDomain, 1 > FunctionSpaceType
Definition: space/fourier/functionset.hh:125
FunctionSpaceType::JacobianRangeType JacobianRangeType
Definition: space/fourier/functionset.hh:129
FunctionSpaceType::RangeType RangeType
Definition: space/fourier/functionset.hh:128
static SizeType index(const MultiIndexType &multiIndex)
Definition: space/fourier/functionset.hh:347