dune-fem  2.8-git
shapefunctionset/legendre.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH
2 #define DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH
3 
4 #include <cassert>
5 #include <cstddef>
6 
7 #include <algorithm>
8 #include <array>
9 #include <iterator>
10 #include <vector>
11 
14 
15 namespace Dune
16 {
17 
18  namespace Fem
19  {
20 
21  // LegendreShapeFunction
22  // ---------------------
23 
30  template< class FunctionSpace >
32  {
33  static_assert( FunctionSpace::dimRange == 1, "FunctionSpace must be scalar (i.e., dimRange = 1)." );
34 
35  public:
38 
43 
52 
53  private:
54  static const int dimDomain = FunctionSpaceType::dimDomain;
55 
56  public:
61  LegendreShapeFunction () = default;
62 
63  template< class MultiIndex >
64  explicit LegendreShapeFunction ( const MultiIndex &multiIndex )
65  {
66  using std::begin;
67  using std::end;
68 
69  auto first = begin( multiIndex );
70  auto last = end( multiIndex );
71  assert( std::distance( first, last ) == dimDomain );
72  std::copy( first, last, multiIndex_.begin() );
73  }
74 
78  int order () const noexcept
79  {
80  return *std::max_element( multiIndex_.begin(), multiIndex_.end() );
81  }
82 
84  const std::array< int, FunctionSpaceType::dimDomain > &orders () const noexcept
85  {
86  return multiIndex_;
87  }
88 
90  void evaluate ( const DomainType &x, RangeType &value ) const noexcept
91  {
92  value[ 0 ] = RangeFieldType( 1 );
93  for( int i = 0; i < dimDomain; ++i )
94  value[ 0 ] *= LegendrePolynomials::evaluate( multiIndex_[ i ], x[ i ] );
95  }
96 
98  void jacobian ( const DomainType &x, JacobianRangeType &jacobian ) const noexcept
99  {
101  for( int k = 0; k < dimDomain; ++k )
102  {
103  const RangeFieldType phi = LegendrePolynomials::evaluate( multiIndex_[ k ], x[ k ] );
104  const RangeFieldType dphi = LegendrePolynomials::jacobian( multiIndex_[ k ], x[ k ]);
105  for( int i = 0; i < dimDomain; ++i )
106  jacobian[ 0 ][ i ] *= ( k == i ) ? dphi : phi;
107  }
108  }
109 
111  void hessian ( const DomainType &x, HessianRangeType &hessian ) const noexcept
112  {
113  hessian = HessianRangeType( typename HessianRangeType::value_type( 1 ) );
114  for( int k = 0; k < dimDomain; ++k )
115  {
116  const RangeFieldType phi = LegendrePolynomials::evaluate( multiIndex_[ k ], x[ k ] );
117  const RangeFieldType dphi = LegendrePolynomials::jacobian( multiIndex_[ k ], x[ k ] );
118  for( int i = 0; i < dimDomain; ++i )
119  {
120  hessian[ 0 ][ i ][ i ] *= ( k == i ) ? LegendrePolynomials::hessian( multiIndex_[ i ], x[ i ]) : phi;
121  for( int j = i+1; j < dimDomain; ++j )
122  {
123  RangeFieldType tmp = ( k == i || k == j ) ? dphi : phi;
124  hessian[ 0 ][ i ][ j ] *= tmp;
125  hessian[ 0 ][ j ][ i ] *= tmp;
126  }
127  }
128  }
129  }
130 
131  private:
132  std::array< int, dimDomain > multiIndex_;
133  };
134 
135 
136 
137 #ifndef DOXYGEN
138 
139  namespace __LegendreShapeFunctionSet
140  {
141 
142  // DefaultFactory
143  // --------------
144 
145  template< class FunctionSpace >
146  class DefaultFactory
147  {
148  typedef LegendreShapeFunction< FunctionSpace > ShapeFunctionType;
149 
150  static const int dimDomain = FunctionSpace::dimDomain;
151  typedef std::array< int, dimDomain > MultiIndex;
152 
153  public:
154  explicit DefaultFactory ( int order )
155  : order_( order )
156  {}
157 
158  int order () const noexcept { return order_; }
159 
160  std::size_t size () const noexcept
161  {
162  std::size_t size = 1;
163  for( int i = 0; i < dimDomain; ++i )
164  size *= order() + 1;
165  return size;
166  }
167 
168  template< class InputIterator >
169  void operator() ( InputIterator first ) const noexcept
170  {
171  auto function = [&first]( const MultiIndex &multiIndex )
172  {
173  *first = ShapeFunctionType( multiIndex );
174  ++first;
175  };
176 
177  MultiIndex multiIndex;
178  fill( multiIndex, function, &multiIndex[ 0 ], dimDomain, order() );
179  }
180 
181  private:
182  template< class Function >
183  static void fill ( MultiIndex &multiIndex, Function function,
184  int *begin, std::size_t n, int order )
185  {
186  if( n > 0u )
187  {
188  for( *begin = 0; *begin <= order; *begin += 1 )
189  fill( multiIndex, function, begin+1, n-1, order );
190  }
191  else
192  function( multiIndex );
193  }
194 
195  int order_;
196  };
197 
198  } // namespace __LegendreShapeFunctionSet
199 
200 #endif // #ifndef DOXYGEN
201 
202 
203 
204  // LegendreShapeFunctionSet
205  // ------------------------
206 
216  template< class FunctionSpace, bool hierarchicalOrdering = false >
218  {
219  protected:
221 
222  public:
225 
234 
235  protected:
236  struct Compare
237  {
238  bool operator() ( const ShapeFunctionType &lhs, const ShapeFunctionType &rhs ) noexcept
239  {
240  if( lhs.order() != rhs.order() )
241  return lhs.order() < rhs.order();
242  else
243  {
244  const auto &a = lhs.orders();
245  const auto &b = rhs.orders();
246  return std::lexicographical_compare( a.begin(), a.end(), b.begin(), b.end() );
247  }
248  }
249  };
250 
251  public:
259 
265  : LegendreShapeFunctionSet( __LegendreShapeFunctionSet::DefaultFactory< FunctionSpaceType >( order ) )
266  {}
267 
287  template< class Factory >
288  LegendreShapeFunctionSet ( const Factory &factory )
289  : shapeFunctions_( factory.size() ),
290  order_( factory.order() )
291  {
292  factory( shapeFunctions_.begin() );
293 
294  // if hierarchical ordering was selected sort shape functions
295  // according to polynomial order
296  if( hierarchicalOrdering )
297  {
298  std::sort( shapeFunctions_.begin(), shapeFunctions_.end(), Compare() );
299  }
300  }
301 
305  int order () const noexcept { return order_; }
306 
308  std::size_t size () const noexcept { return shapeFunctions_.size(); }
309 
311  template< class Point, class Functor >
312  void evaluateEach ( const Point &x, Functor functor ) const noexcept
313  {
314  RangeType value;
315  std::size_t i = 0u;
316  for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
317  {
318  shapeFunction.evaluate( coordinate( x ), value );
319  functor( i++, value );
320  }
321  }
322 
324  template< class Point, class Functor >
325  void jacobianEach ( const Point &x, Functor functor ) const noexcept
326  {
327  JacobianRangeType jacobian;
328  std::size_t i = 0u;
329  for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
330  {
331  shapeFunction.jacobian( coordinate( x ), jacobian );
332  functor( i++, jacobian );
333  }
334  }
335 
337  template< class Point, class Functor >
338  void hessianEach ( const Point &x, Functor functor ) const noexcept
339  {
340  HessianRangeType hessian;
341  std::size_t i = 0u;
342  for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
343  {
344  shapeFunction.hessian( coordinate( x ), hessian );
345  functor( i++, hessian );
346  }
347  }
348 
349  protected:
350  std::vector< ShapeFunctionType > shapeFunctions_;
351  int order_;
352  };
353 
354  } // namespace Fem
355 
356 } // namespace Dune
357 
358 #endif // #ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH
Definition: bindguard.hh:11
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
Definition: explicitfieldvector.hh:75
A vector valued function space.
Definition: functionspace.hh:60
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
implementation of a single scalar-valued Legendre shape function
Definition: shapefunctionset/legendre.hh:32
FunctionSpaceType::HessianRangeType HessianRangeType
hessian type
Definition: shapefunctionset/legendre.hh:51
void evaluate(const DomainType &x, RangeType &value) const noexcept
evaluate the function
Definition: shapefunctionset/legendre.hh:90
void jacobian(const DomainType &x, JacobianRangeType &jacobian) const noexcept
evaluate the Jacobian of the function
Definition: shapefunctionset/legendre.hh:98
void hessian(const DomainType &x, HessianRangeType &hessian) const noexcept
evaluate the hessian of the function
Definition: shapefunctionset/legendre.hh:111
FunctionSpaceType::RangeFieldType RangeFieldType
field type of range
Definition: shapefunctionset/legendre.hh:42
const std::array< int, FunctionSpaceType::dimDomain > & orders() const noexcept
return monomial orders of this function
Definition: shapefunctionset/legendre.hh:84
LegendreShapeFunction(const MultiIndex &multiIndex)
Definition: shapefunctionset/legendre.hh:64
FunctionSpace FunctionSpaceType
type of function space this function belongs to
Definition: shapefunctionset/legendre.hh:33
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian type
Definition: shapefunctionset/legendre.hh:49
int order() const noexcept
return polynomial order of this function
Definition: shapefunctionset/legendre.hh:78
FunctionSpaceType::DomainType DomainType
domain type
Definition: shapefunctionset/legendre.hh:45
FunctionSpaceType::RangeType RangeType
range type
Definition: shapefunctionset/legendre.hh:47
FunctionSpaceType::DomainFieldType DomainFieldType
field type of domain
Definition: shapefunctionset/legendre.hh:40
a Dune::Fem::ShapeFunctionSet of Legendre ansatz polynomials
Definition: shapefunctionset/legendre.hh:218
int order_
Definition: shapefunctionset/legendre.hh:351
FunctionSpaceType::JacobianRangeType JacobianRangeType
jacobian range type
Definition: shapefunctionset/legendre.hh:231
FunctionSpaceType::DomainType DomainType
domain type
Definition: shapefunctionset/legendre.hh:227
void hessianEach(const Point &x, Functor functor) const noexcept
evalute hessian of each shape function
Definition: shapefunctionset/legendre.hh:338
void evaluateEach(const Point &x, Functor functor) const noexcept
evalute each shape function
Definition: shapefunctionset/legendre.hh:312
FunctionSpaceType::RangeType RangeType
range type
Definition: shapefunctionset/legendre.hh:229
LegendreShapeFunctionSet()=default
default constructor resulting in uninitialized shape function set
int order() const noexcept
return order of shape functions
Definition: shapefunctionset/legendre.hh:305
FunctionSpaceType::HessianRangeType HessianRangeType
hessian range type
Definition: shapefunctionset/legendre.hh:233
void jacobianEach(const Point &x, Functor functor) const noexcept
evalute jacobian of each shape function
Definition: shapefunctionset/legendre.hh:325
std::size_t size() const noexcept
return number of shape functions
Definition: shapefunctionset/legendre.hh:308
FunctionSpace FunctionSpaceType
function space type
Definition: shapefunctionset/legendre.hh:224
LegendreShapeFunctionSet(const Factory &factory)
initialize from user-defined factory object
Definition: shapefunctionset/legendre.hh:288
LegendreShapeFunction< FunctionSpace > ShapeFunctionType
Definition: shapefunctionset/legendre.hh:220
LegendreShapeFunctionSet(int order)
initialize with polynomial order
Definition: shapefunctionset/legendre.hh:264
std::vector< ShapeFunctionType > shapeFunctions_
Definition: shapefunctionset/legendre.hh:350
Definition: shapefunctionset/legendre.hh:237
bool operator()(const ShapeFunctionType &lhs, const ShapeFunctionType &rhs) noexcept
Definition: shapefunctionset/legendre.hh:238
static double hessian(const int num, const double x)
Definition: legendrepolynomials.hh:48
static double jacobian(const int num, const double x)
Definition: legendrepolynomials.hh:34
static double evaluate(const int num, const double x)
Definition: legendrepolynomials.hh:24