dune-localfunctions  2.3.1
monomlocalbasis.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_MONOMLOCALBASIS_HH
4 #define DUNE_MONOMLOCALBASIS_HH
5 
6 #include <cassert>
7 
8 #include <dune/common/fmatrix.hh>
9 
10 #include "../common/localbasis.hh"
11 
12 namespace Dune
13 {
14  namespace MonomImp {
18  template<int d, int k>
19  struct Size {
20  enum { val = Size<d,k-1>::val+Size<d-1,k>::val };
21  };
22  template<int d>
23  struct Size<d, 0> {
24  enum { val = 1 };
25  };
26  template<int k>
27  struct Size<0, k> {
28  enum { val = 1 };
29  };
30  template<>
31  struct Size<0, 0> {
32  enum { val = 1 };
33  };
34 
35  template<class T>
36  T ipow(T base, int exp)
37  {
38  T result(1);
39  while (exp)
40  {
41  if (exp & 1)
42  result *= base;
43  exp >>= 1;
44  base *= base;
45  }
46  return result;
47  }
48 
50  template <typename Traits>
51  class EvalAccess {
52  std::vector<typename Traits::RangeType> &out;
53 #ifndef NDEBUG
54  unsigned int first_unused_index;
55 #endif
56 
57  public:
58  EvalAccess(std::vector<typename Traits::RangeType> &out_)
59  : out(out_)
60 #ifndef NDEBUG
61  , first_unused_index(0)
62 #endif
63  { }
64 #ifndef NDEBUG
66  assert(first_unused_index == out.size());
67  }
68 #endif
69  typename Traits::RangeFieldType &operator[](unsigned int index)
70  {
71  assert(index < out.size());
72 #ifndef NDEBUG
73  if(first_unused_index <= index)
74  first_unused_index = index+1;
75 #endif
76  return out[index][0];
77  }
78  };
79 
81  template <typename Traits>
83  std::vector<typename Traits::JacobianType> &out;
84  unsigned int row;
85 #ifndef NDEBUG
86  unsigned int first_unused_index;
87 #endif
88 
89  public:
90  JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
91  unsigned int row_)
92  : out(out_), row(row_)
93 #ifndef NDEBUG
94  , first_unused_index(0)
95 #endif
96  { }
97 #ifndef NDEBUG
99  assert(first_unused_index == out.size());
100  }
101 #endif
102  typename Traits::RangeFieldType &operator[](unsigned int index)
103  {
104  assert(index < out.size());
105 #ifndef NDEBUG
106  if(first_unused_index <= index)
107  first_unused_index = index+1;
108 #endif
109  return out[index][0][row];
110  }
111  };
112 
125  template <typename Traits, int c>
126  struct Evaluate
127  {
128  enum {
130  d = Traits::dimDomain - c
131  };
138  template <typename Access>
139  static void eval (
140  const typename Traits::DomainType &in,
143  const array<int, Traits::dimDomain> &derivatives,
146  typename Traits::RangeFieldType prod,
148  int bound,
150  int& index,
152  Access &access)
153  {
154  // start with the highest exponent for this dimension, then work down
155  for (int e = bound; e >= 0; --e)
156  {
157  // the rest rest of the available exponents, to be used by the other
158  // dimensions
159  int newbound = bound - e;
160  if(e < derivatives[d])
162  eval(in, derivatives, 0, newbound, index, access);
163  else {
164  int coeff = 1;
165  for(int i = e - derivatives[d] + 1; i <= e; ++i)
166  coeff *= i;
167  // call the evaluator for the next dimension
169  eval( // pass the coordinate and the derivatives unchanged
170  in, derivatives,
171  // also pass the product accumulated so far, but also
172  // include the current dimension
173  prod * ipow(in[d], e-derivatives[d]) * coeff,
174  // pass the number of remaining exponents to the next
175  // dimension
176  newbound,
177  // pass the next index to fill and the output access
178  // wrapper
179  index, access);
180  }
181  }
182  }
183  };
184 
189  template <typename Traits>
190  struct Evaluate<Traits, 1>
191  {
192  enum { d = Traits::dimDomain-1 };
194  template <typename Access>
195  static void eval (const typename Traits::DomainType &in,
196  const array<int, Traits::dimDomain> &derivatives,
197  typename Traits::RangeFieldType prod,
198  int bound, int& index, Access &access)
199  {
200  if(bound < derivatives[d])
201  prod = 0;
202  else {
203  int coeff = 1;
204  for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
205  coeff *= i;
206  prod *= ipow(in[d], bound-derivatives[d]) * coeff;
207  }
208  access[index] = prod;
209  ++index;
210  }
211  };
212 
213  } //namespace MonomImp
214 
229  template<class D, class R, unsigned int d, unsigned int p, unsigned diffOrder = p>
231  {
232  enum { static_size = MonomImp::Size<d,p>::val };
233 
234  public:
236  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,
237  Dune::FieldMatrix<R,1,d>,diffOrder> Traits;
238 
240  unsigned int size () const
241  {
242  return static_size;
243  }
244 
246  inline void evaluateFunction (const typename Traits::DomainType& in,
247  std::vector<typename Traits::RangeType>& out) const
248  {
249  evaluate<0>(array<int, 0>(), in, out);
250  }
251 
253  template<unsigned int k>
254  inline void evaluate (const array<int,k>& directions,
255  const typename Traits::DomainType& in,
256  std::vector<typename Traits::RangeType>& out) const
257  {
258  out.resize(size());
259  int index = 0;
260  array<int, d> derivatives;
261  for(unsigned int i = 0; i < d; ++i) derivatives[i] = 0;
262  for(unsigned int i = 0; i < k; ++i) ++derivatives[directions[i]];
263  MonomImp::EvalAccess<Traits> access(out);
264  for(unsigned int lp = 0; lp <= p; ++lp)
265  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index,
266  access);
267  }
268 
270  inline void
271  evaluateJacobian (const typename Traits::DomainType& in, // position
272  std::vector<typename Traits::JacobianType>& out) const // return value
273  {
274  out.resize(size());
275  array<int, d> derivatives;
276  for(unsigned int i = 0; i < d; ++i)
277  derivatives[i] = 0;
278  for(unsigned int i = 0; i < d; ++i)
279  {
280  derivatives[i] = 1;
281  int index = 0;
282  MonomImp::JacobianAccess<Traits> access(out, i);
283  for(unsigned int lp = 0; lp <= p; ++lp)
284  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
285  derivatives[i] = 0;
286  }
287  }
288 
290  unsigned int order () const
291  {
292  return p;
293  }
294  };
295 
296 }
297 
298 #endif // DUNE_MONOMLOCALBASIS_HH