Dune Core Modules (2.6.0)

brezzidouglasmarini1cube3dlocalbasis.hh
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_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
4#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
5
6#include <array>
7#include <bitset>
8#include <numeric>
9#include <vector>
10
11#include <dune/common/fmatrix.hh>
12
13#include "../../common/localbasis.hh"
14
15namespace Dune
16{
27 template<class D, class R>
28 class BDM1Cube3DLocalBasis
29 {
30
31 public:
32 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,
33 R,3,Dune::FieldVector<R,3>,
34 Dune::FieldMatrix<R,3,3> > Traits;
35
37 BDM1Cube3DLocalBasis()
38 {
39 for (size_t i=0; i<6; i++)
40 sign_[i] = 1.0;
41 }
42
48 BDM1Cube3DLocalBasis(std::bitset<6> s)
49 {
50 for (size_t i=0; i<6; i++)
51 sign_[i] = s[i] ? -1.0 : 1.0;
52 }
53
55 unsigned int size() const
56 {
57 return 18;
58 }
59
66 inline void evaluateFunction(const typename Traits::DomainType& in,
67 std::vector<typename Traits::RangeType>& out) const
68 {
69 out.resize(size());
70
71 out[0][0] = sign_[0] * (in[0] - 1.0);
72 out[0][1] = 0.0;
73 out[0][2] = 0.0;
74 out[1][0] = sign_[1] * in[0];
75 out[1][1] = 0.0;
76 out[1][2] = 0.0;
77 out[2][0] = 0.0;
78 out[2][1] = sign_[2] * (in[1] - 1.0);
79 out[2][2] = 0.0;
80 out[3][0] = 0.0;
81 out[3][1] = sign_[3] * in[1];
82 out[3][2] = 0.0;
83 out[4][0] = 0.0;
84 out[4][1] = 0.0;
85 out[4][2] = sign_[4] * (in[2] - 1.0);
86 out[5][0] = 0.0;
87 out[5][1] = 0.0;
88 out[5][2] = sign_[5] * in[2];
89 out[6][0] = 6.0 * in[0] * in[1] - 3 * in[0]-6 * in[1] + 3.0;
90 out[6][1] = -3.0 * in[1] * in[1] + 3 * in[1];
91 out[6][2] = 0.0;
92 out[7][0] = -6.0 * in[0] * in[1] + 3 * in[0];
93 out[7][1] = 3.0 * in[1] * in[1] - 3 * in[1];
94 out[7][2] = 0.0;
95 out[8][0] = 3.0 * in[0] * in[0] - 3 * in[0];
96 out[8][1] = -6.0 * in[0] * in[1] + 3 * in[1]+6 * in[0]-3.0;
97 out[8][2] = 0.0;
98 out[9][0] = -3.0 * in[0] * in[0] + 3 * in[0];
99 out[9][1] = 6.0 * in[0] * in[1] - 3 * in[1];
100 out[9][2] = 0.0;
101 out[10][0] = -3.0 * in[0] * in[0] + 3 * in[0];
102 out[10][1] = 0.0;
103 out[10][2] = 6.0 * in[0] * in[2]-6 * in[0]-3 * in[2] + 3.0;
104 out[11][0] = 3.0 * in[0] * in[0]-3 * in[0];
105 out[11][1] = 0.0;
106 out[11][2] = -6.0 * in[0] * in[2] + 3 * in[2];
107 out[12][0] = -6.0 * in[0] * in[2]+6 * in[2] + 3 * in[0]-3.0;
108 out[12][1] = 0.0;
109 out[12][2] = 3.0 * in[2] * in[2]-3 * in[2];
110 out[13][0] = -3 * in[0]+6 * in[0] * in[2];
111 out[13][1] = 0.0;
112 out[13][2] = -3.0 * in[2] * in[2] + 3 * in[2];
113 out[14][0] = 0.0;
114 out[14][1] = 6.0 * in[1] * in[2]-3 * in[1]-6 * in[2] + 3.0;
115 out[14][2] = -3 * in[2] * in[2] + 3 * in[2];
116 out[15][0] = 0.0;
117 out[15][1] = -6.0 * in[1] * in[2] + 3 * in[1];
118 out[15][2] = 3.0 * in[2] * in[2]-3 * in[2];
119 out[16][0] = 0.0;
120 out[16][1] = 3.0 * in[1] * in[1]-3 * in[1];
121 out[16][2] = -6.0 * in[1] * in[2] + 3 * in[2]+6 * in[1]-3.0;
122 out[17][0] = 0.0;
123 out[17][1] = -3.0 * in[1] * in[1] + 3 * in[1];
124 out[17][2] = 6.0 * in[1] * in[2] - 3.0 * in[2];
125 }
126
133 inline void evaluateJacobian(const typename Traits::DomainType& in,
134 std::vector<typename Traits::JacobianType>& out) const
135 {
136 out.resize(size());
137
138 out[0][0][0] = sign_[0]; out[0][0][1] = 0; out[0][0][2] = 0;
139 out[0][1][0] = 0; out[0][1][1] = 0; out[0][1][2] = 0;
140 out[0][2][0] = 0; out[0][2][1] = 0; out[0][2][2] = 0;
141
142 out[1][0][0] = sign_[1]; out[1][0][1] = 0; out[1][0][2] = 0;
143 out[1][1][0] = 0; out[1][1][1] = 0; out[1][1][2] = 0;
144 out[1][2][0] = 0; out[1][2][1] = 0; out[1][2][2] = 0;
145
146 out[2][0][0] = 0; out[2][0][1] = 0; out[2][0][2] = 0;
147 out[2][1][0] = 0; out[2][1][1] = sign_[2]; out[2][1][2] = 0;
148 out[2][2][0] = 0; out[2][2][1] = 0; out[2][2][2] = 0;
149
150 out[3][0][0] = 0; out[3][0][1] = 0; out[3][0][2] = 0;
151 out[3][1][0] = 0; out[3][1][1] = sign_[3]; out[3][1][2] = 0;
152 out[3][2][0] = 0; out[3][2][1] = 0; out[3][2][2] = 0;
153
154 out[4][0][0] = 0; out[4][0][1] = 0; out[4][0][2] = 0;
155 out[4][1][0] = 0; out[4][1][1] = 0; out[4][1][2] = 0;
156 out[4][2][0] = 0; out[4][2][1] = 0; out[4][2][2] = sign_[4];
157
158 out[5][0][0] = 0; out[5][0][1] = 0; out[5][0][2] = 0;
159 out[5][1][0] = 0; out[5][1][1] = 0; out[5][1][2] = 0;
160 out[5][2][0] = 0; out[5][2][1] = 0; out[5][2][2] = sign_[5];
161
162
163
164 out[6][0][0] = 6*in[1]-3; out[6][0][1] = 6*in[0]-6; out[6][0][2] = 0.0;
165 out[6][1][0] = 0.0; out[6][1][1] =-6*in[1]+3; out[6][1][2] = 0.0;
166 out[6][2][0] = 0.0; out[6][2][1] = 0.0; out[6][2][2] = 0.0;
167
168 out[7][0][0] = -6*in[1]+3; out[7][0][1] =-6*in[0]; out[7][0][2] = 0.0;
169 out[7][1][0] = 0.0; out[7][1][1] = 6*in[1]-3; out[7][1][2] = 0.0;
170 out[7][2][0] = 0.0; out[7][2][1] = 0.0; out[7][2][2] = 0.0;
171
172
173 out[8][0][0] = 6*in[0]-3; out[8][0][1] = 0.0; out[8][0][2] = 0.0;
174 out[8][1][0] =-6*in[1]+6; out[8][1][1] =-6*in[0]+3; out[8][1][2] = 0.0;
175 out[8][2][0] = 0.0; out[8][2][1] = 0.0; out[8][2][2] = 0.0;
176
177 out[9][0][0] = -6*in[0]+3; out[9][0][1] = 0.0; out[9][0][2] = 0.0;
178 out[9][1][0] = 6*in[1]; out[9][1][1] = 6*in[0]-3; out[9][1][2] = 0.0;
179 out[9][2][0] = 0.0; out[9][2][1] = 0.0; out[9][2][2] = 0.0;
180
181 out[10][0][0] = -6*in[0]+3; out[10][0][1] = 0.0; out[10][0][2] = 0.0;
182 out[10][1][0] = 0.0; out[10][1][1] = 0.0; out[10][1][2] = 0.0;
183 out[10][2][0] =6*in[2]-6; out[10][2][1] = 0.0; out[10][2][2] =6*in[0]-3;
184
185 out[11][0][0] =6*in[0]-3; out[11][0][1] = 0.0; out[11][0][2] = 0.0;
186 out[11][1][0] = 0.0; out[11][1][1] = 0.0; out[11][1][2] = 0.0;
187 out[11][2][0] = -6*in[2]; out[11][2][1] = 0.0; out[11][2][2] = -6*in[0]+3;
188
189
190 out[12][0][0] = -6*in[2]+3; out[12][0][1] =0.0; out[12][0][2] = -6*in[0]+6;
191 out[12][1][0] = 0.0; out[12][1][1] = 0.0; out[12][1][2] = 0.0;
192 out[12][2][0] = 0.0; out[12][2][1] = 0.0; out[12][2][2] =6*in[2]-3;
193
194
195 out[13][0][0] =6*in[2]-3; out[13][0][1] = 0.0; out[13][0][2] = 6*in[0];
196 out[13][1][0] = 0.0; out[13][1][1] = 0.0; out[13][1][2] = 0.0;
197 out[13][2][0] = 0.0; out[13][2][1] = 0.0; out[13][2][2] =-6*in[2]+3;
198
199 out[14][0][0] = 0.0; out[14][0][1] = 0.0; out[14][0][2] = 0.0;
200 out[14][1][0] = 0.0; out[14][1][1] = 6*in[2]-3; out[14][1][2] = 6*in[1]-6;
201 out[14][2][0] = 0.0; out[14][2][1] = 0.0; out[14][2][2] =-6*in[2]+3;
202
203 out[15][0][0] = 0.0; out[15][0][1] = 0.0; out[15][0][2] = 0.0;
204 out[15][1][0] = 0.0; out[15][1][1] =-6*in[2]+3; out[15][1][2] =-6*in[1];
205 out[15][2][0] = 0.0; out[15][2][1] = 0.0; out[15][2][2] = 6*in[2]-3;
206
207 out[16][0][0] = 0.0; out[16][0][1] = 0.0; out[16][0][2] = 0.0;
208 out[16][1][0] = 0.0; out[16][1][1] = 6*in[1]-3; out[16][1][2] = 0.0;
209 out[16][2][0] = 0.0; out[16][2][1] =-6*in[2]+6; out[16][2][2] =-6*in[1]+3;
210
211 out[17][0][0] = 0.0; out[17][0][1] = 0.0; out[17][0][2] = 0.0;
212 out[17][1][0] = 0.0; out[17][1][1] =-6*in[1]+3; out[17][1][2] = 0.0;
213 out[17][2][0] = 0.0; out[17][2][1] = 6*in[2]; out[17][2][2] = 6*in[1]-3;
214 }
215
217 void partial (const std::array<unsigned int, 3>& order,
218 const typename Traits::DomainType& in, // position
219 std::vector<typename Traits::RangeType>& out) const // return value
220 {
221 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
222 if (totalOrder == 0) {
223 evaluateFunction(in, out);
224 } else if (totalOrder == 1) {
225 out.resize(size());
226 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
227
228 switch (direction) {
229 case 0:
230 out[0][0] = sign_[0];
231 out[0][1] = 0;
232 out[0][2] = 0;
233
234 out[1][0] = sign_[1];
235 out[1][1] = 0;
236 out[1][2] = 0;
237
238 out[2][0] = 0;
239 out[2][1] = 0;
240 out[2][2] = 0;
241
242 out[3][0] = 0;
243 out[3][1] = 0;
244 out[3][2] = 0;
245
246 out[4][0] = 0;
247 out[4][1] = 0;
248 out[4][2] = 0;
249
250 out[5][0] = 0;
251 out[5][1] = 0;
252 out[5][2] = 0;
253
254 out[6][0] = 6*in[1]-3;
255 out[6][1] = 0.0;
256 out[6][2] = 0.0;
257
258 out[7][0] = -6*in[1]+3;
259 out[7][1] = 0.0;
260 out[7][2] = 0.0;
261
262 out[8][0] = 6*in[0]-3;
263 out[8][1] =-6*in[1]+6;
264 out[8][2] = 0.0;
265
266 out[9][0] = -6*in[0]+3;
267 out[9][1] = 6*in[1];
268 out[9][2] = 0.0;
269
270 out[10][0] = -6*in[0]+3;
271 out[10][1] = 0.0;
272 out[10][2] =6*in[2]-6;
273
274 out[11][0] =6*in[0]-3;
275 out[11][1] = 0.0;
276 out[11][2] = -6*in[2];
277
278 out[12][0] = -6*in[2]+3;
279 out[12][1] = 0.0;
280 out[12][2] = 0.0;
281
282 out[13][0] =6*in[2]-3;
283 out[13][1] = 0.0;
284 out[13][2] = 0.0;
285
286 out[14][0] = 0.0;
287 out[14][1] = 0.0;
288 out[14][2] = 0.0;
289
290 out[15][0] = 0.0;
291 out[15][1] = 0.0;
292 out[15][2] = 0.0;
293
294 out[16][0] = 0.0;
295 out[16][1] = 0.0;
296 out[16][2] = 0.0;
297
298 out[17][0] = 0.0;
299 out[17][1] = 0.0;
300 out[17][2] = 0.0;
301 break;
302 case 1:
303 out[0][0] = 0;
304 out[0][1] = 0;
305 out[0][2] = 0;
306
307 out[1][0] = 0;
308 out[1][1] = 0;
309 out[1][2] = 0;
310
311 out[2][0] = 0;
312 out[2][1] = sign_[2];
313 out[2][2] = 0;
314
315 out[3][0] = 0;
316 out[3][1] = sign_[3];
317 out[3][2] = 0;
318
319 out[4][0] = 0;
320 out[4][1] = 0;
321 out[4][2] = 0;
322
323 out[5][0] = 0;
324 out[5][1] = 0;
325 out[5][2] = 0;
326
327 out[6][0] = 6*in[0]-6;
328 out[6][1] =-6*in[1]+3;
329 out[6][2] = 0.0;
330
331 out[7][0] =-6*in[0];
332 out[7][1] = 6*in[1]-3;
333 out[7][2] = 0.0;
334
335 out[8][0] = 0.0;
336 out[8][1] =-6*in[0]+3;
337 out[8][2] = 0.0;
338
339 out[9][0] = 0.0;
340 out[9][1] = 6*in[0]-3;
341 out[9][2] = 0.0;
342
343 out[10][0] = 0.0;
344 out[10][1] = 0.0;
345 out[10][2] = 0.0;
346
347 out[11][0] = 0.0;
348 out[11][1] = 0.0;
349 out[11][2] = 0.0;
350
351 out[12][0] =0.0;
352 out[12][1] = 0.0;
353 out[12][2] = 0.0;
354
355 out[13][0] = 0.0;
356 out[13][1] = 0.0;
357 out[13][2] = 0.0;
358
359 out[14][0] = 0.0;
360 out[14][1] = 6*in[2]-3;
361 out[14][2] = 0.0;
362
363 out[15][0] = 0.0;
364 out[15][1] =-6*in[2]+3;
365 out[15][2] = 0.0;
366
367 out[16][0] = 0.0;
368 out[16][1] = 6*in[1]-3;
369 out[16][2] =-6*in[2]+6;
370
371 out[17][0] = 0.0;
372 out[17][1] =-6*in[1]+3;
373 out[17][2] = 6*in[2];
374 break;
375 case 2:
376 out[0][0] = 0;
377 out[0][1] = 0;
378 out[0][2] = 0;
379
380 out[1][0] = 0;
381 out[1][1] = 0;
382 out[1][2] = 0;
383
384 out[2][0] = 0;
385 out[2][1] = 0;
386 out[2][2] = 0;
387
388 out[3][0] = 0;
389 out[3][1] = 0;
390 out[3][2] = 0;
391
392 out[4][0] = 0;
393 out[4][1] = 0;
394 out[4][2] = sign_[4];
395
396 out[5][0] = 0;
397 out[5][1] = 0;
398 out[5][2] = sign_[5];
399
400 out[6][0] = 0.0;
401 out[6][1] = 0.0;
402 out[6][2] = 0.0;
403
404 out[7][0] = 0.0;
405 out[7][1] = 0.0;
406 out[7][2] = 0.0;
407
408 out[8][0] = 0.0;
409 out[8][1] = 0.0;
410 out[8][2] = 0.0;
411
412 out[9][0] = 0.0;
413 out[9][1] = 0.0;
414 out[9][2] = 0.0;
415
416 out[10][0] = 0.0;
417 out[10][1] = 0.0;
418 out[10][2] =6*in[0]-3;
419
420 out[11][0] = 0.0;
421 out[11][1] = 0.0;
422 out[11][2] = -6*in[0]+3;
423
424 out[12][0] = -6*in[0]+6;
425 out[12][1] = 0.0;
426 out[12][2] =6*in[2]-3;
427
428 out[13][0] = 6*in[0];
429 out[13][1] = 0.0;
430 out[13][2] =-6*in[2]+3;
431
432 out[14][0] = 0.0;
433 out[14][1] = 6*in[1]-6;
434 out[14][2] =-6*in[2]+3;
435
436 out[15][0] = 0.0;
437 out[15][1] =-6*in[1];
438 out[15][2] = 6*in[2]-3;
439
440 out[16][0] = 0.0;
441 out[16][1] = 0.0;
442 out[16][2] =-6*in[1]+3;
443
444 out[17][0] = 0.0;
445 out[17][1] = 0.0;
446 out[17][2] = 6*in[1]-3;
447 break;
448 default:
449 DUNE_THROW(RangeError, "Component out of range.");
450 }
451 } else {
452 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
453 }
454 }
455
457 unsigned int order() const
458 {
459 return 2;
460 }
461
462 private:
463 std::array<R,6> sign_;
464 };
465} // end namespace Dune
466#endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
Dune::BDM1Cube3DLocalBasis
First order Brezzi-Douglas-Marini shape functions on the reference hexahedron.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:29
Dune::BDM1Cube3DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:133
Dune::BDM1Cube3DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:66
Dune::BDM1Cube3DLocalBasis::partial
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 all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:217
Dune::BDM1Cube3DLocalBasis::BDM1Cube3DLocalBasis
BDM1Cube3DLocalBasis(std::bitset< 6 > s)
Make set number s, where 0 <= s < 64.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:48
Dune::BDM1Cube3DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:457
Dune::BDM1Cube3DLocalBasis::BDM1Cube3DLocalBasis
BDM1Cube3DLocalBasis()
Standard constructor.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:37
Dune::BDM1Cube3DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:55
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:68
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Dune::NotImplemented
Default exception for dummy implementations.
Definition: exceptions.hh:261
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:252
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:331
Dune
Dune namespace.
Definition: alignedallocator.hh:10
Dune::LocalBasisTraits
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)