dune-localfunctions  2.5.0
qklocalbasis.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 
4 #ifndef DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
5 #define DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
6 
7 #include <numeric>
8 
9 #include <dune/common/deprecated.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 #include <dune/common/power.hh>
13 
14 #include <dune/geometry/type.hh>
15 
18 
19 
20 namespace Dune
21 {
34  template<class D, class R, int k, int d>
36  {
37  enum { n = StaticPower<k+1,d>::power };
38 
39  // ith Lagrange polynomial of degree k in one dimension
40  static R p (int i, D x)
41  {
42  R result(1.0);
43  for (int j=0; j<=k; j++)
44  if (j!=i) result *= (k*x-j)/(i-j);
45  return result;
46  }
47 
48  // derivative of ith Lagrange polynomial of degree k in one dimension
49  static R dp (int i, D x)
50  {
51  R result(0.0);
52 
53  for (int j=0; j<=k; j++)
54  if (j!=i)
55  {
56  R prod( (k*1.0)/(i-j) );
57  for (int l=0; l<=k; l++)
58  if (l!=i && l!=j)
59  prod *= (k*x-l)/(i-l);
60  result += prod;
61  }
62  return result;
63  }
64 
65  // Return i as a d-digit number in the (k+1)-nary system
66  static Dune::FieldVector<int,d> multiindex (int i)
67  {
68  Dune::FieldVector<int,d> alpha;
69  for (int j=0; j<d; j++)
70  {
71  alpha[j] = i % (k+1);
72  i = i/(k+1);
73  }
74  return alpha;
75  }
76 
77  public:
78  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 1> Traits;
79 
81  unsigned int size () const
82  {
83  return StaticPower<k+1,d>::power;
84  }
85 
87  inline void evaluateFunction (const typename Traits::DomainType& in,
88  std::vector<typename Traits::RangeType>& out) const
89  {
90  out.resize(size());
91  for (size_t i=0; i<size(); i++)
92  {
93  // convert index i to multiindex
94  Dune::FieldVector<int,d> alpha(multiindex(i));
95 
96  // initialize product
97  out[i] = 1.0;
98 
99  // dimension by dimension
100  for (int j=0; j<d; j++)
101  out[i] *= p(alpha[j],in[j]);
102  }
103  }
104 
109  inline void
111  std::vector<typename Traits::JacobianType>& out) const
112  {
113  out.resize(size());
114 
115  // Loop over all shape functions
116  for (size_t i=0; i<size(); i++)
117  {
118  // convert index i to multiindex
119  Dune::FieldVector<int,d> alpha(multiindex(i));
120 
121  // Loop over all coordinate directions
122  for (int j=0; j<d; j++)
123  {
124  // Initialize: the overall expression is a product
125  // if j-th bit of i is set to -1, else 1
126  out[i][0][j] = dp(alpha[j],in[j]);
127 
128  // rest of the product
129  for (int l=0; l<d; l++)
130  if (l!=j)
131  out[i][0][j] *= p(alpha[l],in[l]);
132  }
133  }
134  }
135 
141  inline void partial(const std::array<unsigned int,d>& order,
142  const typename Traits::DomainType& in,
143  std::vector<typename Traits::RangeType>& out) const
144  {
145  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
146 
147  switch (totalOrder)
148  {
149  case 0:
150  evaluateFunction(in,out);
151  break;
152  case 1:
153  {
154  out.resize(size());
155 
156  // Loop over all shape functions
157  for (size_t i=0; i<size(); i++)
158  {
159  // convert index i to multiindex
160  Dune::FieldVector<int,d> alpha(multiindex(i));
161 
162  // Initialize: the overall expression is a product
163  out[i][0] = 1.0;
164 
165  // rest of the product
166  for (std::size_t l=0; l<d; l++)
167  out[i][0] *= (order[l]) ? dp(alpha[l],in[l]) : p(alpha[l],in[l]);
168  }
169  break;
170  }
171  default:
172  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
173  }
174  }
175 
181  template<int diffOrder>
182  inline void DUNE_DEPRECATED_MSG("Use method 'partial' instead!")
184  const std::array<int,1>& direction,
185  const typename Traits::DomainType& in,
186  std::vector<typename Traits::RangeType>& out) const
187  {
188  static_assert(diffOrder == 1, "We only can compute first derivatives");
189  out.resize(size());
190 
191  // Loop over all shape functions
192  for (size_t i=0; i<size(); i++)
193  {
194  // convert index i to multiindex
195  Dune::FieldVector<int,d> alpha(multiindex(i));
196 
197  // Loop over all coordinate directions
198  std::size_t j = direction[0];
199 
200  // Initialize: the overall expression is a product
201  // if j-th bit of i is set to -1, else 1
202  out[i][0] = dp(alpha[j],in[j]);
203 
204  // rest of the product
205  for (std::size_t l=0; l<d; l++)
206  if (l!=j)
207  out[i][0] *= p(alpha[l],in[l]);
208  }
209  }
210 
212  unsigned int order () const
213  {
214  return k;
215  }
216  };
217 }
218 
219 #endif
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qklocalbasis.hh:87
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qklocalbasis.hh:110
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
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: qklocalbasis.hh:141
void evaluate(const std::array< int, 1 > &direction, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate derivative in a given direction.
Definition: qklocalbasis.hh:183
D DomainType
domain type
Definition: localbasis.hh:49
unsigned int size() const
number of shape functions
Definition: qklocalbasis.hh:81
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, 1 > Traits
Definition: qklocalbasis.hh:78
STL namespace.
unsigned int order() const
Polynomial order of the shape functions.
Definition: qklocalbasis.hh:212
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
Lagrange shape functions of order k on the reference cube.
Definition: qklocalbasis.hh:35