Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
dubiner.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_DUBINER_ICC
2#define _RHEOLEF_DUBINER_ICC
23//
24// evaluate the Proriol-Dubiner basis on a point of the reference element
25//
26// author: Pierre.Saramito@imag.fr
27//
28// date: 5 september 2017
29//
30#include "rheolef/quadrature.h"
31
32#include "basis_option.h"
33
34namespace rheolef {
35
36// ----------------------------------------------
37// coefficients in the recurrence formula of the Jacobi polynomials:
38// see War-2010, page 5:4
39// ----------------------------------------------
40namespace details {
41
42namespace legendre {
43template<class T>
44static inline T a (size_t i) {
45 return T(int(i))/sqrt(T(2*int(i)+1)*T(2*int(i)-1));
46}
47} // namespace legendre
48
49namespace jacobi {
50template<class T>
51static inline T a (T alpha, T beta, size_t degree) {
52 return (2*degree + 1 + alpha + beta)*(2*degree + 2 + alpha + beta)
53 /(2*(degree + 1)*(degree + 1 + alpha + beta));
54}
55template<class T>
56static inline T b (T alpha, T beta, size_t degree) {
57 return (sqr(alpha) - sqr(beta))*(2*degree + 1 + alpha + beta)
58 /(2*(degree + 1)*(2*degree + alpha + beta)*(degree + 1 + alpha + beta));
59}
60template<class T>
61static inline T c (T alpha, T beta, size_t degree) {
62 return (degree + alpha)*(degree + beta)*(2*degree + 2 + alpha + beta)
63 /((degree + 1)*(degree + 1 + alpha + beta)*(2*degree + alpha + beta));
64}
65} // namespace jacobi
66} // namespace details
67
68// ----------------------------------------------
69// Legendre/Dubiner basis evaluation:
70// - Dubiner basis [Kir-2010] when K=t,T
71// - Legendre basis when K=e,q,H
72// - product basis when K=P
73// ----------------------------------------------
74// ortho-normalized Legendre polynoms on a segment: see e.g. HesWar-2008, p. 45
75// Note: HesWar-2008 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
76template<class T, class Container>
77static
78void
79eval_dubiner_basis_e (
80 const T& tilde_x, // in [-1,1]
81 size_t degree,
82 const std::vector<size_t>& perm,
83 Container& value,
84 std::vector<point_basic<size_t> >& power_index)
85{
87 typedef typename float_traits<T>::type float_type;
88 value.resize (degree+1);
89 power_index.resize (value.size());
90 size_t loc_idof;
92 power_index[perm[loc_idof]] = point_basic<size_t>(0);
93 value [perm[loc_idof]] = 1/sqrt(float_type(2));
94 if (degree == 0) return;
96 power_index[perm[loc_idof]] = point_basic<size_t>(1);
97 value [perm[loc_idof]] = sqrt(float_type(3)/float_type(2))*tilde_x;
98 for (size_type i = 1; i+1 <= degree; i++) {
99 namespace L = rheolef::details::legendre;
103 power_index[perm[loc_idof]] = point_basic<size_t>(i+1);
104 value [perm[loc_idof]] = (tilde_x*value[perm[loc_idof_i]] - L::a<float_type>(i)*value[perm[loc_idof_im1]])/L::a<float_type>(i+1);
105 }
106}
107// orthogonal Dubiner polynoms on a triangle: see Kir-2010, p. 5:7
108// Notes: polynoms are not normalized => mass matrix is diag but not identity
109// Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
110template<class T, class Container>
111static
112void
113eval_dubiner_basis_t (
114 const point_basic<T>& tilde_x, // in [-1,1]^2
115 size_t degree,
116 const std::vector<size_t>& perm,
117 Container& value,
118 std::vector<point_basic<size_t> >& power_index)
119{
121 typedef typename float_traits<T>::type float_type;
122 value.resize ((degree+1)*(degree+2)/2);
123 power_index.resize (value.size());
124 // Dubiner, by algorithm Kir-2010, page 5:7
125 T xi0 = tilde_x[0],
126 xi1 = tilde_x[1];
127 size_t loc_idof;
128 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
129 power_index[perm[loc_idof]] = point_basic<size_t>(0,0);
130 value [perm[loc_idof]] = 1;
131 if (degree == 0) return;
132 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(1,0));
133 power_index[perm[loc_idof]] = point_basic<size_t>(1,0);
134 value [perm[loc_idof]] = (1 + 2*xi0 + xi1)/float_type(2);
135 for (size_type i = 1; i < degree; i++) {
136 size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i ,0));
137 size_t loc_idof_im1_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i-1,0));
138 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i+1,0));
139 power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0);
140 value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*(1 + 2*xi0 + xi1)/float_type(2)*value[perm[loc_idof_i_0]]
141 - float_type( int(i) )/float_type(int(i)+1)*sqr((1 - xi1)/float_type(2)) *value[perm[loc_idof_im1_0]];
142 }
143 for (size_type i = 0; i < degree; i++) {
144 size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
145 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,1));
146 power_index[perm[loc_idof]] = point_basic<size_t>(i,1);
147 value [perm[loc_idof]] = value[perm[loc_idof_i_0]]*(1 + 2*i + (3+2*i)*xi1)/float_type(2);
148 }
149 for (size_type i = 0; i < degree; i++) {
150 for (size_type j = 1; i+j < degree; j++) {
151 namespace J = rheolef::details::jacobi;
152 size_t loc_idof_i_j = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j ));
153 size_t loc_idof_i_jm1 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j-1));
154 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j+1));
155 power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1);
156 value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*xi1 + J::b<float_type>(2*int(i)+1,0,j))*value[perm[loc_idof_i_j]]
157 - J::c<float_type>(2*int(i)+1,0,j) *value[perm[loc_idof_i_jm1]];
158 }
159 }
160}
161template<class T, class Container>
162static
163void
164eval_dubiner_basis_T (
165 const point_basic<T>& tilde_x, // in [-1,1]^3
166 size_t degree,
167 const std::vector<size_t>& perm,
168 Container& value,
169 std::vector<point_basic<size_t> >& power_index)
170{
172 typedef typename float_traits<T>::type float_type;
173 value.resize (reference_element::n_node(reference_element::T,degree));
174 power_index.resize (value.size());
175 // Dubiner, by algorithm Kir-2010, page 5:7
176 T xi0 = tilde_x[0],
177 xi1 = tilde_x[1],
178 xi2 = tilde_x[2];
179 T F1 = (2 + 2*xi0 + xi1 + xi2)/2,
180 F2 = sqr((xi1 + xi2)/2),
181 F3 = (2 + 3*xi1 + xi2)/2,
182 F4 = (1 + 2*xi1 + xi2)/2,
183 F5 = (1 - xi2)/2;
184 size_t loc_idof;
185 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(0,0,0));
186 power_index[perm[loc_idof]] = point_basic<size_t>(0,0,0);
187 value [perm[loc_idof]] = 1;
188 if (degree == 0) return;
189 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(1,0,0));
190 power_index[perm[loc_idof]] = point_basic<size_t>(1,0,0);
191 value [perm[loc_idof]] = F1;
192 for (size_type i = 1; i+1 <= degree; i++) {
193 size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i ,0,0));
194 size_t loc_idof_im1_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i-1,0,0));
195 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i+1,0,0));
196 power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0,0);
197 value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*F1*value[perm[loc_idof_i_0_0]]
198 - float_type( int(i) )/float_type(int(i)+1)*F2*value[perm[loc_idof_im1_0_0]];
199 }
200 for (size_type i = 0; i+1 <= degree; i++) {
201 size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,0,0));
202 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,1,0));
203 power_index[perm[loc_idof]] = point_basic<size_t>(i,1,0);
204 value [perm[loc_idof]] = (i*(1 + xi1) + F3)*value[perm[loc_idof_i_0_0]];
205 }
206 for (size_type i = 0; i+2 <= degree; i++) {
207 for (size_type j = 1; i+j+1 <= degree; j++) {
208 namespace J = rheolef::details::jacobi;
209 size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j ,0));
210 size_t loc_idof_i_jm1_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j-1,0));
211 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j+1,0));
212 power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1,0);
213 value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*F4 + J::b<float_type>(2*int(i)+1,0,j)*F5)*value[perm[loc_idof_i_j_0]]
214 - J::c<float_type>(2*int(i)+1,0,j)*sqr(F5) *value[perm[loc_idof_i_jm1_0]];
215 }}
216 for (size_type i = 0; i+1 <= degree; i++) {
217 for (size_type j = 0; i+j+1 <= degree; j++) {
218 size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,0));
219 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,1));
220 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,1);
221 value [perm[loc_idof]] = (1+i+j + (2+i+j)*xi2)*value[perm[loc_idof_i_j_0]];
222 }}
223 for (size_type i = 0; i+2 <= degree; i++) {
224 for (size_type j = 0; i+j+2 <= degree; j++) {
225 for (size_type k = 1; i+j+k+1 <= degree; k++) {
226 namespace J = rheolef::details::jacobi;
227 size_t loc_idof_i_j_k = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
228 size_t loc_idof_i_j_km1 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k-1));
229 loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k+1));
230 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k+1);
231 value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+2*int(j)+2,0,k)*xi2 + J::b<float_type>(2*int(i)+2*int(j)+2,0,k))*value[perm[loc_idof_i_j_k]]
232 - J::c<float_type>(2*int(i)+2*int(j)+2,0,k) *value[perm[loc_idof_i_j_km1]];
233 }}}
234}
235// --------------------------------------------------
236// main call: evaluate all the basis at a given point
237// --------------------------------------------------
238template<class T, class Container>
239static
240void
241eval_dubiner_basis (
242 const point_basic<T>& hat_x,
243 reference_element hat_K,
244 size_t degree,
245 const std::vector<size_t>& perm,
246 Container& value,
247 std::vector<point_basic<size_t> >& power_index,
248 bool do_tilde = true)
249{
251 value.resize (reference_element::n_node(hat_K.variant(), degree));
252 power_index.resize (value.size());
253 switch (hat_K.variant()) {
255 power_index[perm[0]] = point_basic<size_t>();
256 value [perm[0]] = 1;
257 break;
258 }
260 // HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
261 T tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
262 eval_dubiner_basis_e (tilde_x, degree, perm, value, power_index);
263 break;
264 }
266 // Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
267 point_basic<T> tilde_x
268 = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1) : hat_x;
269 eval_dubiner_basis_t (tilde_x, degree, perm, value, power_index);
270 break;
271 }
273 // Dubiner, by algorithm Kir-2010, page 5:7
274 // Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
275 point_basic<T> tilde_x
276 = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1, 2*hat_x[2]-1) : hat_x;
277 eval_dubiner_basis_T (tilde_x, degree, perm, value, power_index);
278 break;
279 }
281 Container value0, value1;
282 std::vector<size_t> id_e (degree+1);
283 for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
284 std::vector<point_basic<size_t> > power_index_e;
285 eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
286 eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
287 for (size_type i = 0; i <= degree; i++) {
288 for (size_type j = 0; j <= degree; j++) {
289 size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
290 size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
291 size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
292 power_index[perm[loc_idof]] = point_basic<size_t>(i,j);
293 value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e];
294 }}
295 break;
296 }
298 Container value0, value1, value2;
299 std::vector<size_t> id_e (degree+1);
300 for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
301 std::vector<point_basic<size_t> > power_index_e;
302 eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
303 eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
304 eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
305 for (size_type i = 0; i <= degree; i++) {
306 for (size_type j = 0; j <= degree; j++) {
307 for (size_type k = 0; k <= degree; k++) {
308 size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
309 size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
310 size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
311 size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
312 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
313 value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e]*value2[loc_kdof_e];
314 }}}
315 break;
316 }
318 Container value01, value2;
319 std::vector<size_t> id_e (degree+1), id_t ((degree+1)*(degree+2)/2);
320 for (size_type iloc = 0; iloc < id_e.size(); iloc ++) { id_e [iloc ] = iloc ; }
321 for (size_type iloc = 0; iloc < id_t.size(); iloc ++) { id_t [iloc ] = iloc ; }
322 point_basic<T> tilde_x
323 = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1)
324 : point_basic<T>( hat_x[0], hat_x[1]);
325 std::vector<point_basic<size_t> > power_index_e, power_index_t;
326 eval_dubiner_basis_t (tilde_x, degree, id_t, value01, power_index_t);
327 eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
328 for (size_type i = 0; i <= degree; i++) {
329 for (size_type j = 0; i+j <= degree; j++) {
330 for (size_type k = 0; k <= degree; k++) {
331 size_t loc_ijdof_t = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
332 size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
333 size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
334 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
335 value [perm[loc_idof]] = value01[loc_ijdof_t]*value2[loc_kdof_e];
336 }}}
337 break;
338 }
339 default:
340 error_macro ("unsupported element: "<<hat_K.name());
341 }
342}
343// omit the power_index last arg
344template<class T, class Container>
345static
346void
347eval_dubiner_basis (
348 const point_basic<T>& hat_x,
349 reference_element hat_K,
350 size_t degree,
351 const std::vector<size_t>& perm,
352 Container& value,
353 bool do_tilde = true)
354{
355 std::vector<point_basic<size_t> > power_index;
356 eval_dubiner_basis (hat_x, hat_K, degree, perm, value, power_index, do_tilde);
357}
358
359
360}// namespace rheolef
361#endif // _RHEOLEF_DUBINER_ICC
basis_option - finite element method options
field::size_type size_type
Definition branch.cc:430
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
std::vector< int >::size_type size_type
static size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.