dune-localfunctions  2.4
monomiallocalbasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
4 #define DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
5 
6 #include <array>
7 #include <cassert>
8 
9 #include <dune/common/fmatrix.hh>
10 
11 #include "../common/localbasis.hh"
12 
13 namespace Dune
14 {
15  namespace MonomImp {
19  template<int d, int k>
20  struct Size {
21  enum { val = Size<d,k-1>::val+Size<d-1,k>::val };
22  };
23  template<int d>
24  struct Size<d, 0> {
25  enum { val = 1 };
26  };
27  template<int k>
28  struct Size<0, k> {
29  enum { val = 1 };
30  };
31  template<>
32  struct Size<0, 0> {
33  enum { val = 1 };
34  };
35 
36  template<class T>
37  T ipow(T base, int exp)
38  {
39  T result(1);
40  while (exp)
41  {
42  if (exp & 1)
43  result *= base;
44  exp >>= 1;
45  base *= base;
46  }
47  return result;
48  }
49 
51  template <typename Traits>
52  class EvalAccess {
53  std::vector<typename Traits::RangeType> &out;
54 #ifndef NDEBUG
55  unsigned int first_unused_index;
56 #endif
57 
58  public:
59  EvalAccess(std::vector<typename Traits::RangeType> &out_)
60  : out(out_)
61 #ifndef NDEBUG
62  , first_unused_index(0)
63 #endif
64  { }
65 #ifndef NDEBUG
67  assert(first_unused_index == out.size());
68  }
69 #endif
70  typename Traits::RangeFieldType &operator[](unsigned int index)
71  {
72  assert(index < out.size());
73 #ifndef NDEBUG
74  if(first_unused_index <= index)
75  first_unused_index = index+1;
76 #endif
77  return out[index][0];
78  }
79  };
80 
82  template <typename Traits>
84  std::vector<typename Traits::JacobianType> &out;
85  unsigned int row;
86 #ifndef NDEBUG
87  unsigned int first_unused_index;
88 #endif
89 
90  public:
91  JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
92  unsigned int row_)
93  : out(out_), row(row_)
94 #ifndef NDEBUG
95  , first_unused_index(0)
96 #endif
97  { }
98 #ifndef NDEBUG
100  assert(first_unused_index == out.size());
101  }
102 #endif
103  typename Traits::RangeFieldType &operator[](unsigned int index)
104  {
105  assert(index < out.size());
106 #ifndef NDEBUG
107  if(first_unused_index <= index)
108  first_unused_index = index+1;
109 #endif
110  return out[index][0][row];
111  }
112  };
113 
126  template <typename Traits, int c>
127  struct Evaluate
128  {
129  enum {
131  d = Traits::dimDomain - c
132  };
139  template <typename Access>
140  static void eval (
141  const typename Traits::DomainType &in,
144  const std::array<int, Traits::dimDomain> &derivatives,
147  typename Traits::RangeFieldType prod,
149  int bound,
151  int& index,
153  Access &access)
154  {
155  // start with the highest exponent for this dimension, then work down
156  for (int e = bound; e >= 0; --e)
157  {
158  // the rest of the available exponents, to be used by the other
159  // dimensions
160  int newbound = bound - e;
161  if(e < derivatives[d])
163  eval(in, derivatives, 0, newbound, index, access);
164  else {
165  int coeff = 1;
166  for(int i = e - derivatives[d] + 1; i <= e; ++i)
167  coeff *= i;
168  // call the evaluator for the next dimension
170  eval( // pass the coordinate and the derivatives unchanged
171  in, derivatives,
172  // also pass the product accumulated so far, but also
173  // include the current dimension
174  prod * ipow(in[d], e-derivatives[d]) * coeff,
175  // pass the number of remaining exponents to the next
176  // dimension
177  newbound,
178  // pass the next index to fill and the output access
179  // wrapper
180  index, access);
181  }
182  }
183  }
184  };
185 
190  template <typename Traits>
191  struct Evaluate<Traits, 1>
192  {
193  enum { d = Traits::dimDomain-1 };
195  template <typename Access>
196  static void eval (const typename Traits::DomainType &in,
197  const std::array<int, Traits::dimDomain> &derivatives,
198  typename Traits::RangeFieldType prod,
199  int bound, int& index, Access &access)
200  {
201  if(bound < derivatives[d])
202  prod = 0;
203  else {
204  int coeff = 1;
205  for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
206  coeff *= i;
207  prod *= ipow(in[d], bound-derivatives[d]) * coeff;
208  }
209  access[index] = prod;
210  ++index;
211  }
212  };
213 
214  } //namespace MonomImp
215 
230  template<class D, class R, unsigned int d, unsigned int p, unsigned diffOrder = p>
232  {
233  public:
235  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,
236  Dune::FieldMatrix<R,1,d>,diffOrder> Traits;
237 
239  unsigned int size () const
240  {
242  }
243 
245  inline void evaluateFunction (const typename Traits::DomainType& in,
246  std::vector<typename Traits::RangeType>& out) const
247  {
248  evaluate<0>(std::array<int, 0>(), in, out);
249  }
250 
252  template<unsigned int k>
253  inline void evaluate (const std::array<int,k>& directions,
254  const typename Traits::DomainType& in,
255  std::vector<typename Traits::RangeType>& out) const
256  {
257  out.resize(size());
258  int index = 0;
259  std::array<int, d> derivatives;
260  for(unsigned int i = 0; i < d; ++i) derivatives[i] = 0;
261  for(unsigned int i = 0; i < k; ++i) ++derivatives[directions[i]];
262  MonomImp::EvalAccess<Traits> access(out);
263  for(unsigned int lp = 0; lp <= p; ++lp)
264  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index,
265  access);
266  }
267 
269  inline void
270  evaluateJacobian (const typename Traits::DomainType& in, // position
271  std::vector<typename Traits::JacobianType>& out) const // return value
272  {
273  out.resize(size());
274  std::array<int, d> derivatives;
275  for(unsigned int i = 0; i < d; ++i)
276  derivatives[i] = 0;
277  for(unsigned int i = 0; i < d; ++i)
278  {
279  derivatives[i] = 1;
280  int index = 0;
281  MonomImp::JacobianAccess<Traits> access(out, i);
282  for(unsigned int lp = 0; lp <= p; ++lp)
283  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
284  derivatives[i] = 0;
285  }
286  }
287 
289  unsigned int order () const
290  {
291  return p;
292  }
293  };
294 
295 }
296 
297 #endif // DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: monomiallocalbasis.hh:245
T ipow(T base, int exp)
Definition: monomiallocalbasis.hh:37
Access output vector of evaluateJacobian()
Definition: monomiallocalbasis.hh:83
unsigned int order() const
Polynomial order of the shape functions.
Definition: monomiallocalbasis.hh:289
static void eval(const typename Traits::DomainType &in, const std::array< int, Traits::dimDomain > &derivatives, typename Traits::RangeFieldType prod, int bound, int &index, Access &access)
Definition: monomiallocalbasis.hh:140
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, diffOrder > Traits
export type traits for function signature
Definition: monomiallocalbasis.hh:236
Definition: monomiallocalbasis.hh:21
~JacobianAccess()
Definition: monomiallocalbasis.hh:99
Definition: monomiallocalbasis.hh:127
Traits::RangeFieldType & operator[](unsigned int index)
Definition: monomiallocalbasis.hh:103
The next dimension to try for factors.
Definition: monomiallocalbasis.hh:131
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: monomiallocalbasis.hh:270
Definition: monomiallocalbasis.hh:20
Constant shape function.
Definition: monomiallocalbasis.hh:231
JacobianAccess(std::vector< typename Traits::JacobianType > &out_, unsigned int row_)
Definition: monomiallocalbasis.hh:91
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
EvalAccess(std::vector< typename Traits::RangeType > &out_)
Definition: monomiallocalbasis.hh:59
static void eval(const typename Traits::DomainType &in, const std::array< int, Traits::dimDomain > &derivatives, typename Traits::RangeFieldType prod, int bound, int &index, Access &access)
Definition: monomiallocalbasis.hh:196
Access output vector of evaluateFunction() and evaluate()
Definition: monomiallocalbasis.hh:52
D DomainType
domain type
Definition: localbasis.hh:49
void evaluate(const std::array< int, k > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
return given derivative of all components
Definition: monomiallocalbasis.hh:253
Traits::RangeFieldType & operator[](unsigned int index)
Definition: monomiallocalbasis.hh:70
~EvalAccess()
Definition: monomiallocalbasis.hh:66
unsigned int size() const
number of shape functions
Definition: monomiallocalbasis.hh:239