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