dune-localfunctions  2.5.1
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 #include <numeric>
9 
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/deprecated.hh>
12 
13 #include "../common/localbasis.hh"
14 
15 namespace Dune
16 {
17  namespace MonomImp {
21  template<int d, int k>
22  struct Size {
23  enum { val = Size<d,k-1>::val+Size<d-1,k>::val };
24  };
25  template<int d>
26  struct Size<d, 0> {
27  enum { val = 1 };
28  };
29  template<int k>
30  struct Size<0, k> {
31  enum { val = 1 };
32  };
33  template<>
34  struct Size<0, 0> {
35  enum { val = 1 };
36  };
37 
38  template<class T>
39  T ipow(T base, int exp)
40  {
41  T result(1);
42  while (exp)
43  {
44  if (exp & 1)
45  result *= base;
46  exp >>= 1;
47  base *= base;
48  }
49  return result;
50  }
51 
53  template <typename Traits>
54  class EvalAccess {
55  std::vector<typename Traits::RangeType> &out;
56 #ifndef NDEBUG
57  unsigned int first_unused_index;
58 #endif
59 
60  public:
61  EvalAccess(std::vector<typename Traits::RangeType> &out_)
62  : out(out_)
63 #ifndef NDEBUG
64  , first_unused_index(0)
65 #endif
66  { }
67 #ifndef NDEBUG
69  assert(first_unused_index == out.size());
70  }
71 #endif
72  typename Traits::RangeFieldType &operator[](unsigned int index)
73  {
74  assert(index < out.size());
75 #ifndef NDEBUG
76  if(first_unused_index <= index)
77  first_unused_index = index+1;
78 #endif
79  return out[index][0];
80  }
81  };
82 
84  template <typename Traits>
86  std::vector<typename Traits::JacobianType> &out;
87  unsigned int row;
88 #ifndef NDEBUG
89  unsigned int first_unused_index;
90 #endif
91 
92  public:
93  JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
94  unsigned int row_)
95  : out(out_), row(row_)
96 #ifndef NDEBUG
97  , first_unused_index(0)
98 #endif
99  { }
100 #ifndef NDEBUG
102  assert(first_unused_index == out.size());
103  }
104 #endif
105  typename Traits::RangeFieldType &operator[](unsigned int index)
106  {
107  assert(index < out.size());
108 #ifndef NDEBUG
109  if(first_unused_index <= index)
110  first_unused_index = index+1;
111 #endif
112  return out[index][0][row];
113  }
114  };
115 
128  template <typename Traits, int c>
129  struct Evaluate
130  {
131  enum {
133  d = Traits::dimDomain - c
134  };
141  template <typename Access>
142  static void eval (
143  const typename Traits::DomainType &in,
146  const std::array<int, Traits::dimDomain> &derivatives,
149  typename Traits::RangeFieldType prod,
151  int bound,
153  int& index,
155  Access &access)
156  {
157  // start with the highest exponent for this dimension, then work down
158  for (int e = bound; e >= 0; --e)
159  {
160  // the rest of the available exponents, to be used by the other
161  // dimensions
162  int newbound = bound - e;
163  if(e < derivatives[d])
165  eval(in, derivatives, 0, newbound, index, access);
166  else {
167  int coeff = 1;
168  for(int i = e - derivatives[d] + 1; i <= e; ++i)
169  coeff *= i;
170  // call the evaluator for the next dimension
172  eval( // pass the coordinate and the derivatives unchanged
173  in, derivatives,
174  // also pass the product accumulated so far, but also
175  // include the current dimension
176  prod * ipow(in[d], e-derivatives[d]) * coeff,
177  // pass the number of remaining exponents to the next
178  // dimension
179  newbound,
180  // pass the next index to fill and the output access
181  // wrapper
182  index, access);
183  }
184  }
185  }
186  };
187 
192  template <typename Traits>
193  struct Evaluate<Traits, 1>
194  {
195  enum { d = Traits::dimDomain-1 };
197  template <typename Access>
198  static void eval (const typename Traits::DomainType &in,
199  const std::array<int, Traits::dimDomain> &derivatives,
200  typename Traits::RangeFieldType prod,
201  int bound, int& index, Access &access)
202  {
203  if(bound < derivatives[d])
204  prod = 0;
205  else {
206  int coeff = 1;
207  for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
208  coeff *= i;
209  prod *= ipow(in[d], bound-derivatives[d]) * coeff;
210  }
211  access[index] = prod;
212  ++index;
213  }
214  };
215 
216  } //namespace MonomImp
217 
232  template<class D, class R, unsigned int d, unsigned int p, unsigned diffOrder = p>
234  {
235  public:
237  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,
238  Dune::FieldMatrix<R,1,d>,diffOrder> Traits;
239 
241  unsigned int size () const
242  {
244  }
245 
247  inline void evaluateFunction (const typename Traits::DomainType& in,
248  std::vector<typename Traits::RangeType>& out) const
249  {
250  DUNE_NO_DEPRECATED_BEGIN
251  evaluate<0>(std::array<int, 0>(), in, out);
252  DUNE_NO_DEPRECATED_END
253  }
254 
260  inline void partial(const std::array<unsigned int,d>& order,
261  const typename Traits::DomainType& in,
262  std::vector<typename Traits::RangeType>& out) const
263  {
264  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
265 
266  switch (totalOrder)
267  {
268  case 0:
269  evaluateFunction(in,out);
270  break;
271  case 1:
272  {
273  std::array<int,1> directions;
274  directions[0] = std::find(order.begin(), order.end(), 1)-order.begin();
275  DUNE_NO_DEPRECATED_BEGIN
276  evaluate<1>(directions, in, out);
277  DUNE_NO_DEPRECATED_END
278  break;
279  }
280  case 2:
281  {
282  std::array<int,2> directions;
283  unsigned int counter = 0;
284  auto nonconstOrder = order; // need a copy that I can modify
285  for (unsigned int i=0; i<d; i++)
286  {
287  while (nonconstOrder[i])
288  {
289  directions[counter++] = i;
290  nonconstOrder[i]--;
291  }
292  }
293 
294  DUNE_NO_DEPRECATED_BEGIN
295  evaluate<2>(directions, in, out);
296  DUNE_NO_DEPRECATED_END
297  break;
298  }
299  default:
300  // \todo The 'evaluate' method implements higher derivatives
301  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
302  }
303  }
304 
306  template<unsigned int k>
307  inline void DUNE_DEPRECATED_MSG("Use method 'partial' instead!")
308  evaluate (const std::array<int,k>& directions,
309  const typename Traits::DomainType& in,
310  std::vector<typename Traits::RangeType>& out) const
311  {
312  out.resize(size());
313  int index = 0;
314  std::array<int, d> derivatives;
315  for(unsigned int i = 0; i < d; ++i) derivatives[i] = 0;
316  for(unsigned int i = 0; i < k; ++i) ++derivatives[directions[i]];
317  MonomImp::EvalAccess<Traits> access(out);
318  for(unsigned int lp = 0; lp <= p; ++lp)
319  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index,
320  access);
321  }
322 
324  inline void
325  evaluateJacobian (const typename Traits::DomainType& in, // position
326  std::vector<typename Traits::JacobianType>& out) const // return value
327  {
328  out.resize(size());
329  std::array<int, d> derivatives;
330  for(unsigned int i = 0; i < d; ++i)
331  derivatives[i] = 0;
332  for(unsigned int i = 0; i < d; ++i)
333  {
334  derivatives[i] = 1;
335  int index = 0;
336  MonomImp::JacobianAccess<Traits> access(out, i);
337  for(unsigned int lp = 0; lp <= p; ++lp)
338  MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
339  derivatives[i] = 0;
340  }
341  }
342 
344  unsigned int order () const
345  {
346  return p;
347  }
348  };
349 
350 }
351 
352 #endif // DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
Access output vector of evaluateFunction() and evaluate()
Definition: monomiallocalbasis.hh:54
EvalAccess(std::vector< typename Traits::RangeType > &out_)
Definition: monomiallocalbasis.hh:61
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:198
T ipow(T base, int exp)
Definition: monomiallocalbasis.hh:39
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:238
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
unsigned int size() const
number of shape functions
Definition: monomiallocalbasis.hh:241
Traits::RangeFieldType & operator[](unsigned int index)
Definition: monomiallocalbasis.hh:105
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: monomiallocalbasis.hh:247
void partial(const std::array< unsigned int, d > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition: monomiallocalbasis.hh:260
~EvalAccess()
Definition: monomiallocalbasis.hh:68
D DomainType
domain type
Definition: localbasis.hh:49
~JacobianAccess()
Definition: monomiallocalbasis.hh:101
Definition: monomiallocalbasis.hh:129
STL namespace.
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:142
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: monomiallocalbasis.hh:325
Definition: monomiallocalbasis.hh:22
Traits::RangeFieldType & operator[](unsigned int index)
Definition: monomiallocalbasis.hh:72
Constant shape function.
Definition: monomiallocalbasis.hh:233
JacobianAccess(std::vector< typename Traits::JacobianType > &out_, unsigned int row_)
Definition: monomiallocalbasis.hh:93
Access output vector of evaluateJacobian()
Definition: monomiallocalbasis.hh:85
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
unsigned int order() const
Polynomial order of the shape functions.
Definition: monomiallocalbasis.hh:344
Definition: monomiallocalbasis.hh:23