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
3 //
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
psi::n_index
static size_t n_index(char hat_K_name, size_t k)
Definition: interpolate_RTk_polynom.icc:110
u_exact::d
size_t d
Definition: interpolate_RTk_polynom.icc:145
div_u::operator()
Float operator()(const point &x) const
Definition: interpolate_RTk_polynom.icc:148
u_exact::w
Float w
Definition: interpolate_RTk_polynom.icc:145
div_u::_u
u_exact _u
Definition: interpolate_RTk_polynom.icc:150
u_exact::operator()
point operator()(const point &x) const
Definition: interpolate_RTk_polynom.icc:126
psi
Definition: interpolate_RTk_polynom.icc:34
div_psi::operator()
Float operator()(const point &x) const
Definition: interpolate_RTk_polynom.icc:118
u_exact::u_exact
u_exact(size_t d1, Float w1=acos(Float(-1)))
Definition: interpolate_RTk_polynom.icc:144
Float
see the Float page for the full documentation
u_exact
Definition: interpolate_RTk_polynom.icc:125
point
see the point page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
psi::psi
psi(char hat_K_name1, size_t i1)
Definition: interpolate_RTk_polynom.icc:109
div_psi::_psi
psi _psi
Definition: interpolate_RTk_polynom.icc:120
psi::i
const size_t i
Definition: interpolate_RTk_polynom.icc:115
div_u
Definition: interpolate_RTk_polynom.icc:147
psi::div
Float div(const point &x) const
Definition: interpolate_RTk_polynom.icc:72
u_exact::div
Float div(const point &x) const
Definition: interpolate_RTk_polynom.icc:135
div_psi
Definition: interpolate_RTk_polynom.icc:117
div_psi::div_psi
div_psi(char hat_K_name, size_t i)
Definition: interpolate_RTk_polynom.icc:119
psi::hat_K_name
const char hat_K_name
Definition: interpolate_RTk_polynom.icc:115
div_u::div_u
div_u(size_t d, Float w=acos(Float(-1)))
Definition: interpolate_RTk_polynom.icc:149
psi::operator()
point operator()(const point &x) const
Definition: interpolate_RTk_polynom.icc:35