dune-localfunctions  2.5.0
pk2dlocalbasis.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_PK2DLOCALBASIS_HH
4 #define DUNE_PK2DLOCALBASIS_HH
5 
6 #include <numeric>
7 
8 #include <dune/common/fmatrix.hh>
9 #include <dune/common/deprecated.hh>
10 
12 
13 namespace Dune
14 {
27  template<class D, class R, unsigned int k>
29  {
30  public:
31 
33  enum {N = (k+1)*(k+2)/2};
34 
38  enum {O = k};
39 
40  typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
41  Dune::FieldMatrix<R,1,2>, 2 > Traits;
42 
45  {
46  for (unsigned int i=0; i<=k; i++)
47  pos_[i] = (1.0*i)/std::max(k,(unsigned int)1);
48  }
49 
51  unsigned int size () const
52  {
53  return N;
54  }
55 
57  inline void evaluateFunction (const typename Traits::DomainType& x,
58  std::vector<typename Traits::RangeType>& out) const
59  {
60  out.resize(N);
61  // specialization for k==0, not clear whether that is needed
62  if (k==0) {
63  out[0] = 1;
64  return;
65  }
66 
67  int n=0;
68  for (unsigned int j=0; j<=k; j++)
69  for (unsigned int i=0; i<=k-j; i++)
70  {
71  out[n] = 1.0;
72  for (unsigned int alpha=0; alpha<i; alpha++)
73  out[n] *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
74  for (unsigned int beta=0; beta<j; beta++)
75  out[n] *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
76  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
77  out[n] *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
78  n++;
79  }
80  }
81 
83  inline void
84  evaluateJacobian (const typename Traits::DomainType& x, // position
85  std::vector<typename Traits::JacobianType>& out) const // return value
86  {
87  out.resize(N);
88 
89  // specialization for k==0, not clear whether that is needed
90  if (k==0) {
91  out[0][0][0] = 0; out[0][0][1] = 0;
92  return;
93  }
94 
95  int n=0;
96  for (unsigned int j=0; j<=k; j++)
97  for (unsigned int i=0; i<=k-j; i++)
98  {
99  // x_0 derivative
100  out[n][0][0] = 0.0;
101  R factor=1.0;
102  for (unsigned int beta=0; beta<j; beta++)
103  factor *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
104  for (unsigned int a=0; a<i; a++)
105  {
106  R product=factor;
107  for (unsigned int alpha=0; alpha<i; alpha++)
108  if (alpha==a)
109  product *= D(1)/(pos_[i]-pos_[alpha]);
110  else
111  product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
112  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
113  product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
114  out[n][0][0] += product;
115  }
116  for (unsigned int c=i+j+1; c<=k; c++)
117  {
118  R product=factor;
119  for (unsigned int alpha=0; alpha<i; alpha++)
120  product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
121  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
122  if (gamma==c)
123  product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
124  else
125  product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
126  out[n][0][0] += product;
127  }
128 
129  // x_1 derivative
130  out[n][0][1] = 0.0;
131  factor = 1.0;
132  for (unsigned int alpha=0; alpha<i; alpha++)
133  factor *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
134  for (unsigned int b=0; b<j; b++)
135  {
136  R product=factor;
137  for (unsigned int beta=0; beta<j; beta++)
138  if (beta==b)
139  product *= D(1)/(pos_[j]-pos_[beta]);
140  else
141  product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
142  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
143  product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
144  out[n][0][1] += product;
145  }
146  for (unsigned int c=i+j+1; c<=k; c++)
147  {
148  R product=factor;
149  for (unsigned int beta=0; beta<j; beta++)
150  product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
151  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
152  if (gamma==c)
153  product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
154  else
155  product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
156  out[n][0][1] += product;
157  }
158 
159  n++;
160  }
161 
162  }
163 
169  void partial(const std::array<unsigned int,2>& order,
170  const typename Traits::DomainType& in,
171  std::vector<typename Traits::RangeType>& out) const
172  {
173  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
174 
175  switch (totalOrder)
176  {
177  case 0:
178  evaluateFunction(in,out);
179  break;
180  case 1:
181  {
182  std::array<int,1> directions;
183  directions[0] = std::find(order.begin(), order.end(), 1)-order.begin();
184  DUNE_NO_DEPRECATED_BEGIN
185  evaluate<1>(directions, in, out);
186  DUNE_NO_DEPRECATED_END
187  break;
188  }
189  case 2:
190  {
191  std::array<int,2> directions;
192  unsigned int counter = 0;
193  auto nonconstOrder = order; // need a copy that I can modify
194  for (int i=0; i<2; i++)
195  {
196  while (nonconstOrder[i])
197  {
198  directions[counter++] = i;
199  nonconstOrder[i]--;
200  }
201  }
202 
203  DUNE_NO_DEPRECATED_BEGIN
204  evaluate<2>(directions, in, out);
205  DUNE_NO_DEPRECATED_END
206  break;
207  }
208  default:
209  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
210  }
211  }
212 
214  template<unsigned int dOrder> //order of derivative
215  inline void DUNE_DEPRECATED_MSG("Use method 'partial' instead!")
216  evaluate(const std::array<int,dOrder>& directions, //direction of derivative
217  const typename Traits::DomainType& in, //position
218  std::vector<typename Traits::RangeType>& out) const //return value
219  {
220  out.resize(N);
221 
222  if (dOrder > Traits::diffOrder)
223  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
224 
225  if (dOrder==0)
226  evaluateFunction(in, out);
227  else if (dOrder==1)
228  {
229  int n=0;
230  for (unsigned int j=0; j<=k; j++)
231  for (unsigned int i=0; i<=k-j; i++, n++)
232  {
233  out[n] = 0.0;
234  for (unsigned int no1=0; no1 < k; no1++)
235  {
236  R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
237  for (unsigned int no2=0; no2 < k; no2++)
238  if (no1 != no2)
239  factor *= lagrangianFactor(no2, i, j, in);
240 
241  out[n] += factor;
242  }
243  }
244  }
245  else if (dOrder==2)
246  {
247  // specialization for k<2, not clear whether that is needed
248  if (k<2) {
249  std::fill(out.begin(), out.end(), 0.0);
250  return;
251  }
252 
253  //f = prod_{i} f_i -> dxa dxb f = sum_{i} {dxa f_i sum_{k \neq i} dxb f_k prod_{l \neq k,i} f_l
254  int n=0;
255  for (unsigned int j=0; j<=k; j++)
256  for (unsigned int i=0; i<=k-j; i++, n++)
257  {
258  R res = 0.0;
259 
260  for (unsigned int no1=0; no1 < k; no1++)
261  {
262  R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
263  for (unsigned int no2=0; no2 < k; no2++)
264  {
265  if (no1 == no2)
266  continue;
267  R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
268  for (unsigned int no3=0; no3 < k; no3++)
269  {
270  if (no3 == no1 || no3 == no2)
271  continue;
272  factor2 *= lagrangianFactor(no3, i, j, in);
273  }
274  res += factor2;
275  }
276 
277  }
278  out[n] = res;
279  }
280  }
281  }
282 
284  unsigned int order () const
285  {
286  return k;
287  }
288 
289  private:
291  typename Traits::RangeType lagrangianFactor(const int no, const int i, const int j, const typename Traits::DomainType& x) const
292  {
293  if ( no < i)
294  return (x[0]-pos_[no])/(pos_[i]-pos_[no]);
295  if (no < i+j)
296  return (x[1]-pos_[no-i])/(pos_[j]-pos_[no-i]);
297  return (pos_[no+1]-x[0]-x[1])/(pos_[no+1]-pos_[i]-pos_[j]);
298  }
299 
303  typename Traits::RangeType lagrangianFactorDerivative(const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x) const
304  {
305  if ( no < i)
306  return (direction == 0) ? 1.0/(pos_[i]-pos_[no]) : 0;
307 
308  if (no < i+j)
309  return (direction == 0) ? 0: 1.0/(pos_[j]-pos_[no-i]);
310 
311  return -1.0/(pos_[no+1]-pos_[i]-pos_[j]);
312  }
313 
314  D pos_[k+1]; // positions on the interval
315  };
316 
317 }
318 #endif
R RangeType
range type
Definition: localbasis.hh:61
unsigned int size() const
number of shape functions
Definition: pk2dlocalbasis.hh:51
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition: pk2dlocalbasis.hh:41
Definition: pk2dlocalbasis.hh:33
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
number of partial derivatives supported
Definition: localbasis.hh:74
D DomainType
domain type
Definition: localbasis.hh:49
Definition: pk2dlocalbasis.hh:38
Lagrange shape functions of arbitrary order on the reference triangle.
Definition: pk2dlocalbasis.hh:28
STL namespace.
Pk2DLocalBasis()
Standard constructor.
Definition: pk2dlocalbasis.hh:44
unsigned int order() const
Polynomial order of the shape functions.
Definition: pk2dlocalbasis.hh:284
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pk2dlocalbasis.hh:84
void partial(const std::array< unsigned int, 2 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition: pk2dlocalbasis.hh:169
void evaluate(const std::array< int, dOrder > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate higher derivatives of all shape functions.
Definition: pk2dlocalbasis.hh:216
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pk2dlocalbasis.hh:57