dune-localfunctions 3.0-git
pk3dlocalbasis.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_PK3DLOCALBASIS_HH
4#define DUNE_PK3DLOCALBASIS_HH
5
6#include <numeric>
7
8#include <dune/common/fmatrix.hh>
9
11
12namespace Dune
13{
26 template<class D, class R, unsigned int k>
28 {
29 public:
30 enum {N = (k+1)*(k+2)*(k+3)/6};
31 enum {O = k};
32
33 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
34 Dune::FieldMatrix<R,1,3> > Traits;
35
38
40 unsigned int size () const
41 {
42 return N;
43 }
44
46 inline void evaluateFunction (const typename Traits::DomainType& x,
47 std::vector<typename Traits::RangeType>& out) const
48 {
49 out.resize(N);
50 typename Traits::DomainType kx = x;
51 kx *= k;
52 unsigned int n = 0;
53 unsigned int i[4];
54 R factor[4];
55 for (i[2] = 0; i[2] <= k; ++i[2])
56 {
57 factor[2] = 1.0;
58 for (unsigned int j = 0; j < i[2]; ++j)
59 factor[2] *= (kx[2]-j) / (i[2]-j);
60 for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
61 {
62 factor[1] = 1.0;
63 for (unsigned int j = 0; j < i[1]; ++j)
64 factor[1] *= (kx[1]-j) / (i[1]-j);
65 for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
66 {
67 factor[0] = 1.0;
68 for (unsigned int j = 0; j < i[0]; ++j)
69 factor[0] *= (kx[0]-j) / (i[0]-j);
70 i[3] = k - i[0] - i[1] - i[2];
71 D kx3 = k - kx[0] - kx[1] - kx[2];
72 factor[3] = 1.0;
73 for (unsigned int j = 0; j < i[3]; ++j)
74 factor[3] *= (kx3-j) / (i[3]-j);
75 out[n++] = factor[0] * factor[1] * factor[2] * factor[3];
76 }
77 }
78 }
79 }
80
82 inline void
83 evaluateJacobian (const typename Traits::DomainType& x, // position
84 std::vector<typename Traits::JacobianType>& out) const // return value
85 {
86 out.resize(N);
87 typename Traits::DomainType kx = x;
88 kx *= k;
89 unsigned int n = 0;
90 unsigned int i[4];
91 R factor[4];
92 for (i[2] = 0; i[2] <= k; ++i[2])
93 {
94 factor[2] = 1.0;
95 for (unsigned int j = 0; j < i[2]; ++j)
96 factor[2] *= (kx[2]-j) / (i[2]-j);
97 for (i[1] = 0; i[1] <= k - i[2]; ++i[1])
98 {
99 factor[1] = 1.0;
100 for (unsigned int j = 0; j < i[1]; ++j)
101 factor[1] *= (kx[1]-j) / (i[1]-j);
102 for (i[0] = 0; i[0] <= k - i[1] - i[2]; ++i[0])
103 {
104 factor[0] = 1.0;
105 for (unsigned int j = 0; j < i[0]; ++j)
106 factor[0] *= (kx[0]-j) / (i[0]-j);
107 i[3] = k - i[0] - i[1] - i[2];
108 D kx3 = k - kx[0] - kx[1] - kx[2];
109 R sum3 = 0.0;
110 factor[3] = 1.0;
111 for (unsigned int j = 0; j < i[3]; ++j)
112 factor[3] /= i[3] - j;
113 R prod_all = factor[0] * factor[1] * factor[2] * factor[3];
114 for (unsigned int j = 0; j < i[3]; ++j)
115 {
116 R prod = prod_all;
117 for (unsigned int l = 0; l < i[3]; ++l)
118 if (j == l)
119 prod *= -R(k);
120 else
121 prod *= kx3 - l;
122 sum3 += prod;
123 }
124 for (unsigned int j = 0; j < i[3]; ++j)
125 factor[3] *= kx3 - j;
126 for (unsigned int m = 0; m < 3; ++m)
127 {
128 out[n][0][m] = sum3;
129 for (unsigned int j = 0; j < i[m]; ++j)
130 {
131 R prod = factor[3];
132 for (unsigned int p = 0; p < 3; ++p)
133 {
134 if (m == p)
135 for (unsigned int l = 0; l < i[p]; ++l)
136 if (j == l)
137 prod *= R(k) / (i[p]-l);
138 else
139 prod *= (kx[p]-l) / (i[p]-l);
140 else
141 prod *= factor[p];
142 }
143 out[n][0][m] += prod;
144 }
145 }
146 n++;
147 }
148 }
149 }
150 }
151
157 void partial(const std::array<unsigned int,3>& order,
158 const typename Traits::DomainType& in,
159 std::vector<typename Traits::RangeType>& out) const
160 {
161 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
162 if (totalOrder == 0) {
163 evaluateFunction(in, out);
164 } else {
165 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
166 }
167 }
168
170 unsigned int order () const
171 {
172 return k;
173 }
174 };
175
176
177 //Specialization for k=0
178 template<class D, class R>
179 class Pk3DLocalBasis<D,R,0>
180 {
181 public:
182 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
183 Dune::FieldMatrix<R,1,3> > Traits;
184
186 enum {N = 1};
187
189 enum {O = 0};
190
191 unsigned int size () const
192 {
193 return 1;
194 }
195
196 inline void evaluateFunction (const typename Traits::DomainType& in,
197 std::vector<typename Traits::RangeType>& out) const
198 {
199 out.resize(1);
200 out[0] = 1;
201 }
202
203 // evaluate derivative of a single component
204 inline void
205 evaluateJacobian (const typename Traits::DomainType& in, // position
206 std::vector<typename Traits::JacobianType>& out) const // return value
207 {
208 out.resize(1);
209 out[0][0][0] = 0;
210 out[0][0][1] = 0;
211 out[0][0][2] = 0;
212 }
213
219 void partial(const std::array<unsigned int,3>& order,
220 const typename Traits::DomainType& in,
221 std::vector<typename Traits::RangeType>& out) const
222 {
223 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
224 if (totalOrder == 0) {
225 evaluateFunction(in, out);
226 } else {
227 out.resize(N);
228 out[0] = 0;
229 }
230 }
231
232 // local interpolation of a function
233 template<typename E, typename F, typename C>
234 void interpolate (const E& e, const F& f, std::vector<C>& out) const
235 {
236 typename Traits::DomainType x;
237 typename Traits::RangeType y;
238 x[0] = 1.0/4.0;
239 x[1] = 1.0/4.0;
240 x[2] = 1.0/4.0;
241 f.eval_local(e,x,y);
242 out[0] = y;
243 }
244
245 unsigned int order () const
246 {
247 return 0;
248 }
249 };
250}
251#endif
Definition brezzidouglasmarini1cube2d.hh:14
Type traits for LocalBasisVirtualInterface.
Definition localbasis.hh:38
R RangeType
range type
Definition localbasis.hh:61
D DomainType
domain type
Definition localbasis.hh:49
Lagrange shape functions of arbitrary order on the reference tetrahedron.
Definition pk3dlocalbasis.hh:28
@ O
Definition pk3dlocalbasis.hh:31
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
Definition pk3dlocalbasis.hh:34
unsigned int order() const
Polynomial order of the shape functions.
Definition pk3dlocalbasis.hh:170
unsigned int size() const
number of shape functions
Definition pk3dlocalbasis.hh:40
@ N
Definition pk3dlocalbasis.hh:30
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition pk3dlocalbasis.hh:157
Pk3DLocalBasis()
Standard constructor.
Definition pk3dlocalbasis.hh:37
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition pk3dlocalbasis.hh:46
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition pk3dlocalbasis.hh:83
unsigned int order() const
Definition pk3dlocalbasis.hh:245
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
Definition pk3dlocalbasis.hh:183
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition pk3dlocalbasis.hh:219
void interpolate(const E &e, const F &f, std::vector< C > &out) const
Definition pk3dlocalbasis.hh:234
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Definition pk3dlocalbasis.hh:196
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Definition pk3dlocalbasis.hh:205
unsigned int size() const
Definition pk3dlocalbasis.hh:191