dune-localfunctions  2.4.1
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 <dune/common/fmatrix.hh>
7 
9 
10 namespace Dune
11 {
24  template<class D, class R, unsigned int k>
26  {
27  public:
28 
30  enum {N = (k+1)*(k+2)/2};
31 
35  enum {O = k};
36 
37  typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
38  Dune::FieldMatrix<R,1,2>, 2 > Traits;
39 
42  {
43  for (unsigned int i=0; i<=k; i++)
44  pos[i] = (1.0*i)/std::max(k,(unsigned int)1);
45  }
46 
48  unsigned int size () const
49  {
50  return N;
51  }
52 
54  inline void evaluateFunction (const typename Traits::DomainType& x,
55  std::vector<typename Traits::RangeType>& out) const
56  {
57  out.resize(N);
58  // specialization for k==0, not clear whether that is needed
59  if (k==0) {
60  out[0] = 1;
61  return;
62  }
63 
64  int n=0;
65  for (unsigned int j=0; j<=k; j++)
66  for (unsigned int i=0; i<=k-j; i++)
67  {
68  out[n] = 1.0;
69  for (unsigned int alpha=0; alpha<i; alpha++)
70  out[n] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
71  for (unsigned int beta=0; beta<j; beta++)
72  out[n] *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
73  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
74  out[n] *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
75  n++;
76  }
77  }
78 
80  inline void
81  evaluateJacobian (const typename Traits::DomainType& x, // position
82  std::vector<typename Traits::JacobianType>& out) const // return value
83  {
84  out.resize(N);
85 
86  // specialization for k==0, not clear whether that is needed
87  if (k==0) {
88  out[0][0][0] = 0; out[0][0][1] = 0;
89  return;
90  }
91 
92  int n=0;
93  for (unsigned int j=0; j<=k; j++)
94  for (unsigned int i=0; i<=k-j; i++)
95  {
96  // x_0 derivative
97  out[n][0][0] = 0.0;
98  R factor=1.0;
99  for (unsigned int beta=0; beta<j; beta++)
100  factor *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
101  for (unsigned int a=0; a<i; a++)
102  {
103  R product=factor;
104  for (unsigned int alpha=0; alpha<i; alpha++)
105  if (alpha==a)
106  product *= D(1)/(pos[i]-pos[alpha]);
107  else
108  product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
109  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
110  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
111  out[n][0][0] += product;
112  }
113  for (unsigned int c=i+j+1; c<=k; c++)
114  {
115  R product=factor;
116  for (unsigned int alpha=0; alpha<i; alpha++)
117  product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
118  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
119  if (gamma==c)
120  product *= -D(1)/(pos[gamma]-pos[i]-pos[j]);
121  else
122  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
123  out[n][0][0] += product;
124  }
125 
126  // x_1 derivative
127  out[n][0][1] = 0.0;
128  factor = 1.0;
129  for (unsigned int alpha=0; alpha<i; alpha++)
130  factor *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
131  for (unsigned int b=0; b<j; b++)
132  {
133  R product=factor;
134  for (unsigned int beta=0; beta<j; beta++)
135  if (beta==b)
136  product *= D(1)/(pos[j]-pos[beta]);
137  else
138  product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
139  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
140  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
141  out[n][0][1] += product;
142  }
143  for (unsigned int c=i+j+1; c<=k; c++)
144  {
145  R product=factor;
146  for (unsigned int beta=0; beta<j; beta++)
147  product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
148  for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
149  if (gamma==c)
150  product *= -D(1)/(pos[gamma]-pos[i]-pos[j]);
151  else
152  product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
153  out[n][0][1] += product;
154  }
155 
156  n++;
157  }
158 
159  }
160 
162  template<unsigned int dOrder> //order of derivative
163  inline void evaluate(const std::array<int,dOrder>& directions, //direction of derivative
164  const typename Traits::DomainType& in, //position
165  std::vector<typename Traits::RangeType>& out) const //return value
166  {
167  out.resize(N);
168 
169  if (dOrder > Traits::diffOrder)
170  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
171 
172  if (dOrder==0)
173  evaluateFunction(in, out);
174  else if (dOrder==1)
175  {
176  int n=0;
177  for (unsigned int j=0; j<=k; j++)
178  for (unsigned int i=0; i<=k-j; i++, n++)
179  {
180  out[n] = 0.0;
181  for (unsigned int no1=0; no1 < k; no1++)
182  {
183  R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
184  for (unsigned int no2=0; no2 < k; no2++)
185  if (no1 != no2)
186  factor *= lagrangianFactor(no2, i, j, in);
187 
188  out[n] += factor;
189  }
190  }
191  }
192  else if (dOrder==2)
193  {
194  // specialization for k<2, not clear whether that is needed
195  if (k<2) {
196  std::fill(out.begin(), out.end(), 0.0);
197  return;
198  }
199 
200  //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
201  int n=0;
202  for (unsigned int j=0; j<=k; j++)
203  for (unsigned int i=0; i<=k-j; i++, n++)
204  {
205  R res = 0.0;
206 
207  for (unsigned int no1=0; no1 < k; no1++)
208  {
209  R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
210  for (unsigned int no2=0; no2 < k; no2++)
211  {
212  if (no1 == no2)
213  continue;
214  R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
215  for (unsigned int no3=0; no3 < k; no3++)
216  {
217  if (no3 == no1 || no3 == no2)
218  continue;
219  factor2 *= lagrangianFactor(no3, i, j, in);
220  }
221  res += factor2;
222  }
223 
224  }
225  out[n] = res;
226  }
227  }
228  }
229 
231  unsigned int order () const
232  {
233  return k;
234  }
235 
236  private:
238  typename Traits::RangeType lagrangianFactor(const int no, const int i, const int j, const typename Traits::DomainType& x) const
239  {
240  if ( no < i)
241  return (x[0]-pos[no])/(pos[i]-pos[no]);
242  if (no < i+j)
243  return (x[1]-pos[no-i])/(pos[j]-pos[no-i]);
244  return (pos[no+1]-x[0]-x[1])/(pos[no+1]-pos[i]-pos[j]);
245  }
246 
250  typename Traits::RangeType lagrangianFactorDerivative(const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x) const
251  {
252  if ( no < i)
253  return (direction == 0) ? 1.0/(pos[i]-pos[no]) : 0;
254 
255  if (no < i+j)
256  return (direction == 0) ? 0: 1.0/(pos[j]-pos[no-i]);
257 
258  return -1.0/(pos[no+1]-pos[i]-pos[j]);
259  }
260 
261  D pos[k+1]; // positions on the interval
262  };
263 
264 }
265 #endif
Lagrange shape functions of arbitrary order on the reference triangle.
Definition: pk2dlocalbasis.hh:25
Definition: pk2dlocalbasis.hh:30
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pk2dlocalbasis.hh:81
unsigned int size() const
number of shape functions
Definition: pk2dlocalbasis.hh:48
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
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:163
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
unsigned int order() const
Polynomial order of the shape functions.
Definition: pk2dlocalbasis.hh:231
Definition: pk2dlocalbasis.hh:35
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition: pk2dlocalbasis.hh:38
number of partial derivatives supported
Definition: localbasis.hh:74
R RangeType
range type
Definition: localbasis.hh:61
D DomainType
domain type
Definition: localbasis.hh:49
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pk2dlocalbasis.hh:54
Pk2DLocalBasis()
Standard constructor.
Definition: pk2dlocalbasis.hh:41