dune-localfunctions  2.2.1
pk1dlocalbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_Pk1DLOCALBASIS_HH
2 #define DUNE_Pk1DLOCALBASIS_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};
29 
31  enum {O = k};
32 
33  typedef LocalBasisTraits<D,
34  1,
35  Dune::FieldVector<D,1>,
36  R,
37  1,
38  Dune::FieldVector<R,1>,
39  Dune::FieldMatrix<R,1,1>
40  > Traits;
41 
44  {
45  for (unsigned int i=0; i<=k; i++)
46  pos[i] = (1.0*i)/std::max(k,(unsigned int)1);
47  }
48 
50  unsigned int size () const
51  {
52  return N;
53  }
54 
56  inline void evaluateFunction (const typename Traits::DomainType& x,
57  std::vector<typename Traits::RangeType>& out) const
58  {
59  out.resize(N);
60 
61  for (unsigned int i=0; i<N; i++)
62  {
63  out[i] = 1.0;
64  for (unsigned int alpha=0; alpha<i; alpha++)
65  out[i] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
66  for (unsigned int gamma=i+1; gamma<=k; gamma++)
67  out[i] *= (x[0]-pos[gamma])/(pos[i]-pos[gamma]);
68  }
69  }
70 
72  inline void
73  evaluateJacobian (const typename Traits::DomainType& x, // position
74  std::vector<typename Traits::JacobianType>& out) const // return value
75  {
76  out.resize(N);
77 
78  for (unsigned int i=0; i<=k; i++) {
79 
80  // x_0 derivative
81  out[i][0][0] = 0.0;
82  R factor=1.0;
83  for (unsigned int a=0; a<i; a++)
84  {
85  R product=factor;
86  for (unsigned int alpha=0; alpha<i; alpha++)
87  product *= (alpha==a) ? 1.0/(pos[i]-pos[alpha])
88  : (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
89  for (unsigned int gamma=i+1; gamma<=k; gamma++)
90  product *= (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
91  out[i][0][0] += product;
92  }
93  for (unsigned int c=i+1; c<=k; c++)
94  {
95  R product=factor;
96  for (unsigned int alpha=0; alpha<i; alpha++)
97  product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
98  for (unsigned int gamma=i+1; gamma<=k; gamma++)
99  product *= (gamma==c) ? -1.0/(pos[gamma]-pos[i])
100  : (pos[gamma]-x[0])/(pos[gamma]-pos[i]);
101  out[i][0][0] += product;
102  }
103  }
104 
105  }
106 
108  unsigned int order () const
109  {
110  return k;
111  }
112 
113  private:
114  R pos[k+1]; // positions on the interval
115  };
116 
117 }
118 #endif