dune-localfunctions  2.2.1
pk2dlocalbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_PK2DLOCALBASIS_HH
2 #define DUNE_PK2DLOCALBASIS_HH
3 
4 #include <dune/common/fmatrix.hh>
5 
7 
8 namespace Dune
9 {
22  template<class D, class R, unsigned int k>
24  {
25  public:
26 
28  enum {N = (k+1)*(k+2)/2};
29 
33  enum {O = k};
34 
35  typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
36  Dune::FieldMatrix<R,1,2> > Traits;
37 
40  {
41  for (unsigned int i=0; i<=k; i++)
42  pos[i] = (1.0*i)/std::max(k,(unsigned int)1);
43  }
44 
46  unsigned int size () const
47  {
48  return N;
49  }
50 
52  inline void evaluateFunction (const typename Traits::DomainType& x,
53  std::vector<typename Traits::RangeType>& out) const
54  {
55  out.resize(N);
56  // specialization for k==0, not clear whether that is needed
57  if (k==0) {
58  out[0] = 1;
59  return;
60  }
61 
62  int n=0;
63  for (unsigned int j=0; j<=k; j++)
64  for (unsigned int i=0; i<=k-j; i++)
65  {
66  out[n] = 1.0;
67  for (unsigned int alpha=0; alpha<i; alpha++)
68  out[n] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
69  for (unsigned int beta=0; beta<j; beta++)
70  out[n] *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
71  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
72  out[n] *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
73  n++;
74  }
75  }
76 
78  inline void
79  evaluateJacobian (const typename Traits::DomainType& x, // position
80  std::vector<typename Traits::JacobianType>& out) const // return value
81  {
82  out.resize(N);
83 
84  // specialization for k==0, not clear whether that is needed
85  if (k==0) {
86  out[0][0][0] = 0; out[0][0][1] = 0;
87  return;
88  }
89 
90  int n=0;
91  for (unsigned int j=0; j<=k; j++)
92  for (unsigned int i=0; i<=k-j; i++)
93  {
94  // x_0 derivative
95  out[n][0][0] = 0.0;
96  R factor=1.0;
97  for (unsigned int beta=0; beta<j; beta++)
98  factor *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
99  for (unsigned int a=0; a<i; a++)
100  {
101  R product=factor;
102  for (unsigned int alpha=0; alpha<i; alpha++)
103  if (alpha==a)
104  product *= 1.0/(pos[i]-pos[alpha]);
105  else
106  product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
107  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
108  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
109  out[n][0][0] += product;
110  }
111  for (unsigned int c=i+j+1; c<=k; c++)
112  {
113  R product=factor;
114  for (unsigned int alpha=0; alpha<i; alpha++)
115  product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
116  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
117  if (gamma==c)
118  product *= -1.0/(pos[gamma]-pos[i]-pos[j]);
119  else
120  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
121  out[n][0][0] += product;
122  }
123 
124  // x_1 derivative
125  out[n][0][1] = 0.0;
126  factor = 1.0;
127  for (unsigned int alpha=0; alpha<i; alpha++)
128  factor *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
129  for (unsigned int b=0; b<j; b++)
130  {
131  R product=factor;
132  for (unsigned int beta=0; beta<j; beta++)
133  if (beta==b)
134  product *= 1.0/(pos[j]-pos[beta]);
135  else
136  product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
137  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
138  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
139  out[n][0][1] += product;
140  }
141  for (unsigned int c=i+j+1; c<=k; c++)
142  {
143  R product=factor;
144  for (unsigned int beta=0; beta<j; beta++)
145  product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
146  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
147  if (gamma==c)
148  product *= -1.0/(pos[gamma]-pos[i]-pos[j]);
149  else
150  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
151  out[n][0][1] += product;
152  }
153 
154  n++;
155  }
156 
157  }
158 
160  unsigned int order () const
161  {
162  return k;
163  }
164 
165 private:
166  R pos[k+1]; // positions on the interval
167 };
168 
169 }
170 #endif