Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
interpolate_RTk_polynom.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_INTERPOLATE_RTK_POLYNOM_ICC
2#define _RHEOLEF_INTERPOLATE_RTK_POLYNOM_ICC
23//
24// check the RT1 basis interpolatation on a triangle or a quadrangle
25//
26// author: Pierre.Saramito@imag.fr
27//
28// date: 29 april 2019
29//
30#include "rheolef/basis.h"
31// -----------------------------------------------------------------------------
32// 1) polynomial function
33// -----------------------------------------------------------------------------
34struct psi {
35 point operator() (const point& x) const {
36 switch (hat_K_name) {
37 case 't': {
38 switch (i) {
39 // P0^2 + (x,y)*P0
40 case 0: return point(1,0);
41 case 1: return point(0,1);
42 case 2: return point(x[0],x[1]);
43 // P1^2 + (x,y)*P1
44 case 3: return point(x[0],-x[1]);
45 case 4: return point(0,x[0]);
46 case 5: return point(x[1],0);
47 case 6: return point(sqr(x[0]),x[0]*x[1]);
48 default: return point(x[0]*x[1],sqr(x[1]));
49 }
50 }
51 case 'q':
52 default: {
53 switch (i) {
54 // P_{1,0}*P_{0,1}
55 case 0: return point(1,0);
56 case 1: return point(0,1);
57 case 2: return point(x[0],0);
58 case 3: return point(0,x[1]);
59 // P_{2,1}*P_{1,2}
60 case 4: return point(0,x[0]);
61 case 5: return point(x[1],0);
62 case 6: return point(x[0]*x[1],0);
63 case 7: return point(0,x[0]*x[1]);
64 case 8: return point(sqr(x[0]),0);
65 case 9: return point(0,sqr(x[1]));
66 case 10: return point(sqr(x[0])*x[1],0);
67 default: return point(0,x[0]*sqr(x[1]));
68 }
69 }
70 }
71 }
72 Float div (const point& x) const {
73 switch (hat_K_name) {
74 case 't': {
75 switch (i) {
76 // P0^2 + (x,y)*P0
77 case 0: return 0;
78 case 1: return 0;
79 case 2: return 2;
80 // P1^2 + (x,y)*P1
81 case 3: return 0;
82 case 4: return 0;
83 case 5: return 0;
84 case 6: return 3*x[0];
85 default: return 3*x[1];
86 }
87 }
88 case 'q':
89 default: {
90 switch (i) {
91 // P_{2,1}*P_{1,2}
92 case 0: return 0;
93 case 1: return 0;
94 case 2: return 1;
95 case 3: return 1;
96 // P_{1,0}*P_{0,1}
97 case 4: return 0;
98 case 5: return 0;
99 case 6: return x[1];
100 case 7: return x[0];
101 case 8: return 2*x[0];
102 case 9: return 2*x[1];
103 case 10: return 2*x[0]*x[1];
104 default: return 2*x[0]*x[1];
105 }
106 }
107 }
108 }
109 psi (char hat_K_name1, size_t i1) : hat_K_name(hat_K_name1), i(i1) {}
110 static size_t n_index (char hat_K_name, size_t k) {
111 return (hat_K_name == 't') ?
112 ((k == 0) ? 3 : 8) :
113 ((k == 0) ? 4 : 12);
114 }
115 const char hat_K_name; const size_t i;
116};
117struct div_psi {
118 Float operator() (const point& x) const { return _psi.div (x); }
119 div_psi (char hat_K_name, size_t i) : _psi(hat_K_name,i) {}
121};
122// -----------------------------------------------------------------------------
123// 2) non-polynomial function
124// -----------------------------------------------------------------------------
125struct u_exact {
126 point operator() (const point& x) const {
127 switch (d) {
128 case 2: return point(cos(w*(x[0]+2*x[1])),
129 sin(w*(x[0]-2*x[1])));
130 default: return point(cos(w*(x[0]+2*x[1]+x[2])),
131 sin(w*(x[0]-2*x[1]-x[2])),
132 sin(w*(x[0]+2*x[1]-x[2])));
133 }
134 }
135 Float div (const point& x) const {
136 switch (d) {
137 case 2: return -w*( sin(w*( x[0]+2*x[1]))
138 + 2*cos(w*(-x[0]+2*x[1])));
139 default: return -w*( sin(w*( x[0]+2*x[1]+x[2]))
140 + 2*cos(w*(-x[0]+2*x[1]+x[2]))
141 + cos(w*(-x[0]-2*x[1]+x[2])));
142 }
143 }
144 u_exact(size_t d1, Float w1 = acos(Float(-1))) : d(d1), w(w1) {}
145 size_t d; Float w;
146};
147struct div_u {
148 Float operator() (const point& x) const { return _u.div(x); }
149 div_u(size_t d, Float w = acos(Float(-1))) : _u(d,w) {}
151};
152
153#endif // _RHEOLEF_INTERPOLATE_RTK_POLYNOM_ICC
see the Float page for the full documentation
see the point page for the full documentation
div_psi(char hat_K_name, size_t i)
Float operator()(const point &x) const
Float operator()(const point &x) const
div_u(size_t d, Float w=acos(Float(-1)))
point operator()(const point &x) const
Float div(const point &x) const
const size_t i
const char hat_K_name
static size_t n_index(char hat_K_name, size_t k)
psi(char hat_K_name1, size_t i1)
point operator()(const point &x) const
Float div(const point &x) const
u_exact(size_t d1, Float w1=acos(Float(-1)))