Rheolef  7.1
an efficient C++ finite element environment
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 // -----------------------------------------------------------------------------
34 struct 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 };
117 struct 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 // -----------------------------------------------------------------------------
125 struct 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 };
147 struct 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)))