Dune Core Modules (2.8.0)

hierarchicalsimplexp2withelementbubble.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_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
4#define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
5
10#include <numeric>
11#include <vector>
12
15
16#include <dune/localfunctions/common/localbasis.hh>
17#include <dune/localfunctions/common/localkey.hh>
18#include <dune/localfunctions/common/localinterpolation.hh>
19
20namespace Dune
21{
22 template<class D, class R, int dim>
23 class HierarchicalSimplexP2WithElementBubbleLocalBasis
24 {
25 public:
26 HierarchicalSimplexP2WithElementBubbleLocalBasis()
27 {
28 DUNE_THROW(Dune::NotImplemented,"HierarchicalSimplexP2LocalBasis not implemented for dim > 3.");
29 }
30 };
31
46 template<class D, class R>
47 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,1>
48 {
49 public:
53
55 unsigned int size () const
56 {
57 return 3;
58 }
59
61 inline void evaluateFunction (const typename Traits::DomainType& in,
62 std::vector<typename Traits::RangeType>& out) const
63 {
64 out.resize(3);
65
66 out[0] = 1-in[0];
67 out[1] = in[0];
68 out[2] = 1-4*(in[0]-0.5)*(in[0]-0.5);
69 }
70
72 inline void
73 evaluateJacobian (const typename Traits::DomainType& in, // position
74 std::vector<typename Traits::JacobianType>& out) const // return value
75 {
76 out.resize(3);
77
78 out[0][0][0] = -1;
79 out[1][0][0] = 1;
80 out[2][0][0] = 4-8*in[0];
81 }
82
84 void partial (const std::array<unsigned int, 1>& order,
85 const typename Traits::DomainType& in, // position
86 std::vector<typename Traits::RangeType>& out) const // return value
87 {
88 auto totalOrder = order[0];
89 if (totalOrder == 0) {
90 evaluateFunction(in, out);
91 } else if (totalOrder == 1) {
92 out.resize(size());
93 out[0] = -1;
94 out[1] = 1;
95 out[2] = 4-8*in[0];
96 } else if (totalOrder == 2) {
97 out.resize(size());
98 out[0] = 0;
99 out[1] = 0;
100 out[2] =-8;
101 } else {
102 out.resize(size());
103 out[0] = out[1] = out[2] = 0;
104 }
105 }
106
109 unsigned int order () const
110 {
111 return 2;
112 }
113
114 };
115
136 template<class D, class R>
137 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,2>
138 {
139 public:
143
145 unsigned int size () const
146 {
147 return 7;
148 }
149
151 inline void evaluateFunction (const typename Traits::DomainType& in,
152 std::vector<typename Traits::RangeType>& out) const
153 {
154 out.resize(7);
155
156 out[0] = 1 - in[0] - in[1];
157 out[1] = 4*in[0]*(1-in[0]-in[1]);
158 out[2] = in[0];
159 out[3] = 4*in[1]*(1-in[0]-in[1]);
160 out[4] = 4*in[0]*in[1];
161 out[5] = in[1];
162 out[6] = 27*in[0]*in[1]*(1-in[0]-in[1]);
163
164 }
165
167 inline void
168 evaluateJacobian (const typename Traits::DomainType& in, // position
169 std::vector<typename Traits::JacobianType>& out) const // return value
170 {
171 out.resize(7);
172
173 out[0][0][0] = -1; out[0][0][1] = -1;
174 out[1][0][0] = 4-8*in[0]-4*in[1]; out[1][0][1] = -4*in[0];
175 out[2][0][0] = 1; out[2][0][1] = 0;
176 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1];
177 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0];
178 out[5][0][0] = 0; out[5][0][1] = 1;
179
180 // Cubic bubble
181 out[6][0][0] = 27 * in[1] * (1 - 2*in[0] - in[1]);
182 out[6][0][1] = 27 * in[0] * (1 - 2*in[1] - in[0]);
183
184 }
185
187 void partial (const std::array<unsigned int, 2>& order,
188 const typename Traits::DomainType& in, // position
189 std::vector<typename Traits::RangeType>& out) const // return value
190 {
191 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
192 if (totalOrder == 0) {
193 evaluateFunction(in, out);
194 } else if (totalOrder == 1) {
195 out.resize(size());
196 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
197
198 switch (direction) {
199 case 0:
200 out[0] = -1;
201 out[1] = 4-8*in[0]-4*in[1];
202 out[2] = 1;
203 out[3] = -4*in[1];
204 out[4] = 4*in[1];
205 out[5] = 0;
206 out[6] = 27 * in[1] * (1 - 2*in[0] - in[1]);
207 break;
208 case 1:
209 out[0] = -1;
210 out[1] = -4*in[0];
211 out[2] = 0;
212 out[3] = 4-4*in[0]-8*in[1];
213 out[4] = 4*in[0];
214 out[5] = 1;
215 out[6] = 27 * in[0] * (1 - 2*in[1] - in[0]);
216 break;
217 default:
218 DUNE_THROW(RangeError, "Component out of range.");
219 }
220 } else {
221 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
222 }
223 }
224
227 unsigned int order () const
228 {
229 return 3;
230 }
231
232 };
233
258 template<class D, class R>
259 class HierarchicalSimplexP2WithElementBubbleLocalBasis<D,R,3>
260 {
261 public:
265
267 unsigned int size () const
268 {
269 return 11;
270 }
271
273 void evaluateFunction (const typename Traits::DomainType& in,
274 std::vector<typename Traits::RangeType>& out) const
275 {
276 out.resize(10);
277
278 out[0] = 1 - in[0] - in[1] - in[2];
279 out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]);
280 out[2] = in[0];
281 out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]);
282 out[4] = 4 * in[0] * in[1];
283 out[5] = in[1];
284 out[6] = 4 * in[2] * (1 - in[0] - in[1] - in[2]);
285 out[7] = 4 * in[0] * in[2];
286 out[8] = 4 * in[1] * in[2];
287 out[9] = in[2];
288
289 // quartic element bubble
290 out[10] = 81*in[0]*in[1]*in[2]*(1-in[0]-in[1]-in[2]);
291 }
292
294 void evaluateJacobian (const typename Traits::DomainType& in, // position
295 std::vector<typename Traits::JacobianType>& out) const // return value
296 {
297 out.resize(10);
298
299 out[0][0][0] = -1; out[0][0][1] = -1; out[0][0][2] = -1;
300 out[1][0][0] = 4-8*in[0]-4*in[1]-4*in[2]; out[1][0][1] = -4*in[0]; out[1][0][2] = -4*in[0];
301 out[2][0][0] = 1; out[2][0][1] = 0; out[2][0][2] = 0;
302 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1]-4*in[2]; out[3][0][2] = -4*in[1];
303 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; out[4][0][2] = 0;
304 out[5][0][0] = 0; out[5][0][1] = 1; out[5][0][2] = 0;
305 out[6][0][0] = -4*in[2]; out[6][0][1] = -4*in[2]; out[6][0][2] = 4-4*in[0]-4*in[1]-8*in[2];
306 out[7][0][0] = 4*in[2]; out[7][0][1] = 0; out[7][0][2] = 4*in[0];
307 out[8][0][0] = 0; out[8][0][1] = 4*in[2]; out[8][0][2] = 4*in[1];
308 out[9][0][0] = 0; out[9][0][1] = 0; out[9][0][2] = 1;
309
310 out[10][0][0] = 81 * in[1] * in[2] * (1 - 2*in[0] - in[1] - in[2]);
311 out[10][0][1] = 81 * in[0] * in[2] * (1 - in[0] - 2*in[1] - in[2]);
312 out[10][0][2] = 81 * in[0] * in[1] * (1 - in[0] - in[1] - 2*in[2]);
313 }
314
316 void partial (const std::array<unsigned int, 3>& order,
317 const typename Traits::DomainType& in, // position
318 std::vector<typename Traits::RangeType>& out) const // return value
319 {
320 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
321 if (totalOrder == 0) {
322 evaluateFunction(in, out);
323 } else if (totalOrder == 1) {
324 out.resize(size());
325 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
326
327 switch (direction) {
328 case 0:
329 out[0] = -1;
330 out[1] = 4-8*in[0]-4*in[1]-4*in[2];
331 out[2] = 1;
332 out[3] = -4*in[1];
333 out[4] = 4*in[1];
334 out[5] = 0;
335 out[6] = -4*in[2];
336 out[7] = 4*in[2];
337 out[8] = 0;
338 out[9] = 0;
339 out[10] = 81 * in[1] * in[2] * (1 - 2*in[0] - in[1] - in[2]);
340 break;
341 case 1:
342 out[0] = -1;
343 out[1] = -4*in[0];
344 out[2] = 0;
345 out[3] = 4-4*in[0]-8*in[1]-4*in[2];
346 out[4] = 4*in[0];
347 out[5] = 1;
348 out[6] = -4*in[2];
349 out[7] = 0;
350 out[8] = 4*in[2];
351 out[9] = 0;
352 out[10] = 81 * in[0] * in[2] * (1 - in[0] - 2*in[1] - in[2]);
353 break;
354 case 2:
355 out[0] = -1;
356 out[1] = -4*in[0];
357 out[2] = 0;
358 out[3] = -4*in[1];
359 out[4] = 0;
360 out[5] = 0;
361 out[6] = 4-4*in[0]-4*in[1]-8*in[2];
362 out[7] = 4*in[0];
363 out[8] = 4*in[1];
364 out[9] = 1;
365 out[10] = 81 * in[0] * in[1] * (1 - in[0] - in[1] - 2*in[2]);
366 break;
367 default:
368 DUNE_THROW(RangeError, "Component out of range.");
369 }
370 } else {
371 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
372 }
373 }
374
377 unsigned int order () const
378 {
379 return 4;
380 }
381
382 };
383
384
410 template <int dim>
412 {
413 // The binomial coefficient: dim+1 over 1
414 static const int numVertices = dim+1;
415
416 // The binomial coefficient: dim+1 over 2
417 static const int numEdges = (dim+1)*dim / 2;
418
419 public:
422 : li(numVertices+numEdges + 1)
423 {
424 if (dim!=2)
425 DUNE_THROW(NotImplemented, "only for 2d");
426
427 li[0] = Dune::LocalKey(0,2,0); // Vertex (0,0)
428 li[1] = Dune::LocalKey(0,1,0); // Edge (0.5, 0)
429 li[2] = Dune::LocalKey(1,2,0); // Vertex (1,0)
430 li[3] = Dune::LocalKey(1,1,0); // Edge (0, 0.5)
431 li[4] = Dune::LocalKey(2,1,0); // Edge (0.5, 0.5)
432 li[5] = Dune::LocalKey(2,2,0); // Vertex (0,1)
433 li[6] = Dune::LocalKey(0,0,0); // Element (1/3, 1/3)
434 }
435
437 size_t size () const
438 {
439 return numVertices+numEdges + 1;
440 }
441
443 const Dune::LocalKey& localKey (size_t i) const
444 {
445 return li[i];
446 }
447
448 private:
449 std::vector<Dune::LocalKey> li;
450 };
451
452 template<class LB>
453 class HierarchicalSimplexP2WithElementBubbleLocalInterpolation
454 {
455 public:
456
458 template<typename F, typename C>
459 void interpolate (const F& ff, std::vector<C>& out) const
460 {
461 typename LB::Traits::DomainType x;
462 typename LB::Traits::RangeType y;
463
464 out.resize(7);
465
466 auto&& f = Impl::makeFunctionWithCallOperator<decltype(x)>(ff);
467
468 // vertices
469 x[0] = 0.0; x[1] = 0.0; out[0] = f(x);
470 x[0] = 1.0; x[1] = 0.0; out[2] = f(x);
471 x[0] = 0.0; x[1] = 1.0; out[5] = f(x);
472
473 // edge bubbles
474 x[0] = 0.5; x[1] = 0.0; y = f(x);
475 out[1] = y - out[0]*(1-x[0]) - out[2]*x[0];
476
477 x[0] = 0.0; x[1] = 0.5; y = f(x);
478 out[3] = y - out[0]*(1-x[1]) - out[5]*x[1];
479
480 x[0] = 0.5; x[1] = 0.5; y = f(x);
481 out[4] = y - out[2]*x[0] - out[5]*x[1];
482
483 // element bubble
484 x[0] = 1.0/3; x[1] = 1.0/3; y = f(x);
485
487 HierarchicalSimplexP2WithElementBubbleLocalBasis<double,double,2> shapeFunctions;
488 std::vector<typename LB::Traits::RangeType> sfValues;
489 shapeFunctions.evaluateFunction(x, sfValues);
490
491 out[6] = y;
492 for (int i=0; i<6; i++)
493 out[6] -= out[i]*sfValues[i];
494
495 }
496
497 };
498
499
500}
501#endif
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:73
LocalBasisTraits< D, 1, Dune::FieldVector< D, 1 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 1 > > Traits
export type traits for function signature
Definition: hierarchicalsimplexp2withelementbubble.hh:52
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:61
unsigned int order() const
Polynomial order of the shape functions (2, in this case)
Definition: hierarchicalsimplexp2withelementbubble.hh:109
unsigned int size() const
number of shape functions
Definition: hierarchicalsimplexp2withelementbubble.hh:55
void partial(const std::array< unsigned int, 1 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:84
unsigned int size() const
number of shape functions
Definition: hierarchicalsimplexp2withelementbubble.hh:145
unsigned int order() const
Polynomial order of the shape functions (3 in this case)
Definition: hierarchicalsimplexp2withelementbubble.hh:227
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:151
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:168
void partial(const std::array< unsigned int, 2 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:187
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 > > Traits
export type traits for function signature
Definition: hierarchicalsimplexp2withelementbubble.hh:142
unsigned int size() const
number of shape functions
Definition: hierarchicalsimplexp2withelementbubble.hh:267
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:273
unsigned int order() const
Polynomial order of the shape functions (4 in this case)
Definition: hierarchicalsimplexp2withelementbubble.hh:377
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: hierarchicalsimplexp2withelementbubble.hh:294
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: hierarchicalsimplexp2withelementbubble.hh:316
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
export type traits for function signature
Definition: hierarchicalsimplexp2withelementbubble.hh:264
The local finite element needed for the Zou-Kornhuber estimator for Signorini problems.
Definition: hierarchicalsimplexp2withelementbubble.hh:412
size_t size() const
number of coefficients
Definition: hierarchicalsimplexp2withelementbubble.hh:437
const Dune::LocalKey & localKey(size_t i) const
get i'th index
Definition: hierarchicalsimplexp2withelementbubble.hh:443
HierarchicalSimplexP2WithElementBubbleLocalCoefficients()
Standard constructor.
Definition: hierarchicalsimplexp2withelementbubble.hh:421
Describe position of one degree of freedom.
Definition: localkey.hh:21
Default exception for dummy implementations.
Definition: exceptions.hh:261
Default exception class for range errors.
Definition: exceptions.hh:252
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Dune namespace.
Definition: alignedallocator.hh:11
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)