dune-localfunctions  2.2.1
monomlocalbasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_MONOMLOCALBASIS_HH
3 #define DUNE_MONOMLOCALBASIS_HH
4 
5 #include <cassert>
6 
7 #include <dune/common/fmatrix.hh>
8 
9 #include"../common/localbasis.hh"
10 
11 namespace Dune
12 {
13  namespace MonomImp {
17  template<int d, int k>
18  struct Size {
19  enum { val = Size<d,k-1>::val+Size<d-1,k>::val };
20  };
21  template<int d>
22  struct Size<d, 0> {
23  enum { val = 1 };
24  };
25  template<int k>
26  struct Size<0, k> {
27  enum { val = 1 };
28  };
29  template<>
30  struct Size<0, 0> {
31  enum { val = 1 };
32  };
33 
34  template<class T>
35  T ipow(T base, int exp)
36  {
37  T result(1);
38  while (exp)
39  {
40  if (exp & 1)
41  result *= base;
42  exp >>= 1;
43  base *= base;
44  }
45  return result;
46  }
47 
49  template <typename Traits>
50  class EvalAccess {
51  std::vector<typename Traits::RangeType> &out;
52 #ifndef NDEBUG
53  unsigned int first_unused_index;
54 #endif
55 
56  public:
57  EvalAccess(std::vector<typename Traits::RangeType> &out_)
58  : out(out_)
59 #ifndef NDEBUG
60  , first_unused_index(0)
61 #endif
62  { }
63 #ifndef NDEBUG
65  assert(first_unused_index == out.size());
66  }
67 #endif
68  typename Traits::RangeFieldType &operator[](unsigned int index)
69  {
70  assert(index < out.size());
71 #ifndef NDEBUG
72  if(first_unused_index <= index)
73  first_unused_index = index+1;
74 #endif
75  return out[index][0];
76  }
77  };
78 
80  template <typename Traits>
82  std::vector<typename Traits::JacobianType> &out;
83  unsigned int row;
84 #ifndef NDEBUG
85  unsigned int first_unused_index;
86 #endif
87 
88  public:
89  JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
90  unsigned int row_)
91  : out(out_), row(row_)
92 #ifndef NDEBUG
93  , first_unused_index(0)
94 #endif
95  { }
96 #ifndef NDEBUG
98  assert(first_unused_index == out.size());
99  }
100 #endif
101  typename Traits::RangeFieldType &operator[](unsigned int index)
102  {
103  assert(index < out.size());
104 #ifndef NDEBUG
105  if(first_unused_index <= index)
106  first_unused_index = index+1;
107 #endif
108  return out[index][0][row];
109  }
110  };
111 
124  template <typename Traits, int c>
125  struct Evaluate
126  {
127  enum {
129  d = Traits::dimDomain - c
130  };
137  template <typename Access>
138  static void eval (
139  const typename Traits::DomainType &in,
142  const array<int, Traits::dimDomain> &derivatives,
145  typename Traits::RangeFieldType prod,
147  int bound,
149  int& index,
151  Access &access)
152  {
153  // start with the highest exponent for this dimension, then work down
154  for (int e = bound; e >= 0; --e)
155  {
156  // the rest rest of the available exponents, to be used by the other
157  // dimensions
158  int newbound = bound - e;
159  if(e < derivatives[d])
161  eval(in, derivatives, 0, newbound, index, access);
162  else {
163  int coeff = 1;
164  for(int i = e - derivatives[d] + 1; i <= e; ++i)
165  coeff *= i;
166  // call the evaluator for the next dimension
168  eval(// pass the coordinate and the derivatives unchanged
169  in, derivatives,
170  // also pass the product accumulated so far, but also
171  // include the current dimension
172  prod * ipow(in[d], e-derivatives[d]) * coeff,
173  // pass the number of remaining exponents to the next
174  // dimension
175  newbound,
176  // pass the next index to fill and the output access
177  // wrapper
178  index, access);
179  }
180  }
181  }
182  };
183 
188  template <typename Traits>
189  struct Evaluate<Traits, 1>
190  {
191  enum { d = Traits::dimDomain-1 };
193  template <typename Access>
194  static void eval (const typename Traits::DomainType &in,
195  const array<int, Traits::dimDomain> &derivatives,
196  typename Traits::RangeFieldType prod,
197  int bound, int& index, Access &access)
198  {
199  if(bound < derivatives[d])
200  prod = 0;
201  else {
202  int coeff = 1;
203  for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
204  coeff *= i;
205  prod *= ipow(in[d], bound-derivatives[d]) * coeff;
206  }
207  access[index] = prod;
208  ++index;
209  }
210  };
211 
212  } //namespace MonomImp
213 
228  template<class D, class R, unsigned int d, unsigned int p, unsigned diffOrder = p>
230  {
231  enum { static_size = MonomImp::Size<d,p>::val };
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  {
241  return static_size;
242  }
243 
245  inline void evaluateFunction (const typename Traits::DomainType& in,
246  std::vector<typename Traits::RangeType>& out) const
247  {
248  evaluate<0>(array<int, 0>(), in, out);
249  }
250 
252  template<unsigned int k>
253  inline void evaluate (const 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  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  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_MONOMLOCALBASIS_HH