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