1#ifndef _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
2#define _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
32#ifndef _RHEOLEF_REFERENCE_ELEMENT_H
58reference_element_variant (
const char* name_table,
char name)
60 for (
size_t i = 0; i < reference_element__max_variant; i++) {
61 if (name_table [i] == name)
return i;
63 return reference_element__max_variant;
67reference_element_dimension_by_name(
char element_name) {
68 switch (element_name) {
76 default:
error_macro(
"unexpected element name `"<<element_name<<
"' (ascii="<<
int(element_name)<<
")");
82reference_element_dimension_by_variant (
size_t variant)
84 return _reference_element_dimension_by_variant[
variant];
88reference_element_variant (
char element_name)
90 switch (element_name) {
91 case 'p':
return reference_element__p;
92 case 'e':
return reference_element__e;
93 case 't':
return reference_element__t;
94 case 'q':
return reference_element__q;
95 case 'T':
return reference_element__T;
96 case 'P':
return reference_element__P;
97 case 'H':
return reference_element__H;
98 default:
error_macro(
"unexpected element name `"<<element_name<<
"'\"=char("<<
int(element_name)<<
")");
104reference_element_n_vertex (
size_t variant)
107 case reference_element__p:
return 1;
108 case reference_element__e:
return 2;
109 case reference_element__t:
return 3;
110 case reference_element__q:
return 4;
111 case reference_element__T:
return 4;
112 case reference_element__P:
return 6;
113 case reference_element__H:
return 8;
114 default:
error_macro (
"unexpected element variant `"<<variant<<
"'");
return 0;
119reference_element_n_node (
size_t variant,
size_t order)
122 case reference_element__p:
return 1;
123 case reference_element__e:
return order+1;
124 case reference_element__t:
return (order+1)*(
order+2)/2;
125 case reference_element__q:
return (order+1)*(
order+1);
126 case reference_element__T:
return (order+1)*(
order+2)*(order+3)/6;
127 case reference_element__P:
return (order+1)*(
order+1)*(order+2)/2;
128 case reference_element__H:
return (order+1)*(
order+1)*(order+1);
129 default:
error_macro (
"unexpected element variant `"<<variant<<
"'");
return 0;
134reference_element_init_local_nnode_by_variant (
136 std::array<size_t,reference_element__max_variant>& sz)
138 check_macro (order > 0,
"invalid order="<<order<<
": expect order > 0");
139 sz [reference_element__p] = 1;
140 sz [reference_element__e] =
order-1;
141 sz [reference_element__t] = (
order-1)*(order-2)/2;
142 sz [reference_element__q] = (
order-1)*(order-1);
143 sz [reference_element__T] = (
order-1)*(order-2)*(
order-3)/6;
144 sz [reference_element__P] = (
order-1)*(order-1)*(
order-2)/2;
145 sz [reference_element__H] = (
order-1)*(order-1)*(
order-1);
152 reference_element__p,
153 reference_element__e,
154 reference_element__t,
155 reference_element__T,
156 reference_element__max_variant
160reference_element_first_variant_by_dimension (
size_t dim)
162 return _first_variant_by_dimension[dim];
166reference_element_last_variant_by_dimension (
size_t dim)
175reference_element_p_first_inod_by_variant (
177 size_t subgeo_variant)
179 if (subgeo_variant == reference_element__p)
return 0;
184reference_element_e_first_inod_by_variant (
186 size_t subgeo_variant)
188 if (subgeo_variant == reference_element__p)
return 0;
189 if (subgeo_variant == reference_element__e)
return 2;
194reference_element_t_first_inod_by_variant (
196 size_t subgeo_variant)
198 if (subgeo_variant == reference_element__p)
return 0;
199 if (subgeo_variant == reference_element__e)
return 3;
200 if (subgeo_variant == reference_element__t)
return 3 + 3*(
order-1);
201 return (order+1)*(
order+2)/2;
205reference_element_q_first_inod_by_variant (
207 size_t subgeo_variant)
209 if (subgeo_variant == reference_element__p)
return 0;
210 if (subgeo_variant == reference_element__e)
return 4;
211 if (subgeo_variant == reference_element__t)
return 4 + 4*(
order-1);
212 if (subgeo_variant == reference_element__q)
return 4 + 4*(
order-1);
213 return (order+1)*(
order+1);
217reference_element_T_first_inod_by_variant (
219 size_t subgeo_variant)
221 if (subgeo_variant == reference_element__p)
return 0;
222 if (subgeo_variant == reference_element__e)
return 4;
223 if (subgeo_variant == reference_element__t)
return 4 + 6*(
order-1);
224 if (subgeo_variant == reference_element__q)
return 4 + 6*(
order-1) + 4*((order-1)*(
order-2)/2);
225 if (subgeo_variant == reference_element__T)
return 4 + 6*(
order-1) + 4*((order-1)*(
order-2)/2);
226 return (order+1)*(
order+2)*(order+3)/6;
230reference_element_P_first_inod_by_variant (
232 size_t subgeo_variant)
234 if (subgeo_variant == reference_element__p)
return 0;
235 if (subgeo_variant == reference_element__e)
return 6;
236 if (subgeo_variant == reference_element__t)
return 6 + 9*(
order-1);
237 if (subgeo_variant == reference_element__q)
return 6 + 9*(
order-1) + 2*((order-1)*(
order-2)/2);
238 if (subgeo_variant == reference_element__T)
return 6 + 9*(
order-1) + 2*((order-1)*(
order-2)/2) + 3*(
order-1)*(order-1);
239 if (subgeo_variant == reference_element__P)
return 6 + 9*(
order-1) + 2*((order-1)*(
order-2)/2) + 3*(
order-1)*(order-1);
240 return ((order+1)*(
order+1)*(order+2))/2;
244reference_element_H_first_inod_by_variant (
246 size_t subgeo_variant)
248 if (subgeo_variant == reference_element__p)
return 0;
249 if (subgeo_variant == reference_element__e)
return 8;
250 if (subgeo_variant == reference_element__t)
return 8 + 12*(
order-1);
251 if (subgeo_variant == reference_element__q)
return 8 + 12*(
order-1);
252 if (subgeo_variant == reference_element__T)
return 8 + 12*(
order-1) + 6*(order-1)*(
order-1);
253 if (subgeo_variant == reference_element__P)
return 8 + 12*(
order-1) + 6*(order-1)*(
order-1);
254 if (subgeo_variant == reference_element__H)
return 8 + 12*(
order-1) + 6*(order-1)*(
order-1);
255 return (order+1)*(
order+1)*(order+1);
259reference_element_first_inod_by_variant (
262 size_t subgeo_variant)
264#define _RHEOLEF_geo_element_auto_case(NAME) \
265 case reference_element__##NAME: \
266 return reference_element_##NAME##_first_inod_by_variant (order, subgeo_variant);
276 default:
error_macro (
"unexpected element variant `"<<variant<<
"'");
return 0;
278#undef _RHEOLEF_geo_element_auto_case
282reference_element_last_inod_by_variant (
285 size_t subgeo_variant)
287 return reference_element_first_inod_by_variant (variant, order, subgeo_variant+1);
294reference_element_p_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
300reference_element_e_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
303 if (i == 0)
return 0;
304 if (i == order)
return 1;
309reference_element_t_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
314 if (i == 0)
return 0;
315 if (i == order)
return 1;
316 return 3 + 0*(
order-1) + (i - 1);
318 if (i + j == order) {
319 if (j == order)
return 2;
320 return 3 + 1*(
order-1) + (j - 1);
322 size_t j1 =
order - j;
324 return 3 + 2*(
order-1) + (j1 - 1);
327 size_t n_face_node = (
order-1)*(order-2)/2;
328 size_t ij = (n_face_node - (j1-1)*j1/2) + (i - 1);
329 return 3 + 3*(
order-1) + ij;
333reference_element_q_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
338 if (i == 0)
return 0;
339 if (i == order)
return 1;
340 return 4 + 0*(
order-1) + (i - 1);
343 if (j == order)
return 2;
344 return 4 + 1*(
order-1) + (j - 1);
346 size_t i1 =
order - i;
348 if (i == 0)
return 3;
349 return 4 + 2*(
order-1) + (i1 - 1);
351 size_t j1 =
order - j;
353 return 4 + 3*(
order-1) + (j1 - 1);
356 size_t ij = (
order-1)*(j-1) + (i-1);
357 size_t n_bdry_node = 4 + 4*(
order-1);
358 return n_bdry_node + ij;
362reference_element_T_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
367 size_t n_tot_edge_node = 6*(
order-1);
368 size_t n_face_node = (
order-1)*(order-2)/2;
372 if (i == 0)
return 0;
373 if (i == order)
return 1;
374 return 4 + 0*(
order-1) + (i - 1);
376 if (i + j == order) {
377 if (j == order)
return 2;
378 return 4 + 1*(
order-1) + (j - 1);
380 size_t j1 =
order - j;
382 return 4 + 2*(
order-1) + (j1 - 1);
385 size_t i1 =
order - i;
386 size_t ji = (n_face_node - i1*(i1-1)/2) + (j - 1);
387 return 4 + n_tot_edge_node + 0*n_face_node + ji;
392 if (k == order)
return 3;
393 return 4 + 3*(
order-1) + (k - 1);
395 if (j + k == order) {
396 return 4 + 5*(
order-1) + (k - 1);
399 size_t j1 =
order - j;
400 size_t kj = (n_face_node - j1*(j1-1)/2) + (k - 1);
401 return 4 + n_tot_edge_node + 1*n_face_node + kj;
405 if (i + k == order) {
406 return 4 + 4*(
order-1) + (k - 1);
409 size_t k1 =
order - k;
410 size_t ik = (n_face_node - k1*(k1-1)/2) + (i - 1);
411 return 4 + n_tot_edge_node + 2*n_face_node + ik;
413 if (i + j + k == order) {
415 size_t k1 =
order - k;
416 size_t jk = (n_face_node - k1*(k1-1)/2) + (j - 1);
417 return 4 + n_tot_edge_node + 3*n_face_node + jk;
420 size_t n_tot_face_node = 4*n_face_node;
421 size_t n_tot_node = (
order+1)*(order+2)*(
order+3)/6;
422 size_t n_volume_node = (
order-1)*(order-2)*(
order-3)/6;
423 size_t k1 =
order - k;
424 size_t j1 =
order - j - k;
425 size_t n_level_k_node = (k1-1)*(k1-2)/2;
426 size_t ijk = (n_volume_node - k1*(k1-1)*(k1-2)/6)
427 + (n_level_k_node - j1*(j1-1)/2)
429 return 4 + n_tot_edge_node + n_tot_face_node + ijk;
433reference_element_P_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
438 size_t i1 =
order - i;
439 size_t j1 =
order - j;
440 size_t k1 =
order - k;
441 size_t n_node_edge =
order-1;
442 size_t n_tot_edge_node = 9*(
order-1);
443 size_t n_node_face_tri = (
order-1)*(order-2)/2;
444 size_t n_node_face_qua = (
order-1)*(order-1);
448 if (i == 0)
return 0;
449 if (i == order)
return 1;
450 return 6 + 0*(
order-1) + (i - 1);
452 if (i + j == order) {
453 if (j == order)
return 2;
454 return 6 + 1*(
order-1) + (j - 1);
457 return 6 + 2*(
order-1) + (j1 - 1);
460 size_t ji = (n_node_face_tri - i1*(i1-1)/2) + (j - 1);
461 return 6 + n_tot_edge_node + 0*n_node_face_tri + ji;
466 if (i == 0)
return 3;
467 if (i == order)
return 4;
468 return 6 + 6*(
order-1) + (i - 1);
470 if (i + j == order) {
471 if (j == order)
return 5;
472 return 6 + 7*(
order-1) + (j - 1);
475 return 6 + 8*(
order-1) + (j1 - 1);
478 size_t ij = (n_node_face_tri - j1*(j1-1)/2) + (i - 1);
479 return 6 + n_tot_edge_node + 1*n_node_face_tri + ij;
485 return 6 + 3*(
order-1) + (k - 1);
489 return 6 + 4*(
order-1) + (k - 1);
492 size_t ik = n_node_edge*(k-1) + (i-1);
493 return 6 + n_tot_edge_node + 2*n_node_face_tri + 0*n_node_face_qua + ik;
499 return 6 + 5*(
order-1) + (k - 1);
502 size_t kj = n_node_edge*(j-1) + (k-1);
503 return 6 + n_tot_edge_node + 2*n_node_face_tri + 2*n_node_face_qua + kj;
505 if (i + j == order) {
507 size_t jk = n_node_edge*(k-1) + (j-1);
508 return 6 + n_tot_edge_node + 2*n_node_face_tri + 1*n_node_face_qua + jk;
511 size_t n_tot_face_node = 2*n_node_face_tri + 3*n_node_face_qua;
512 size_t ij = (n_node_face_tri - (j1-1)*j1/2) + (i - 1);
513 size_t ijk = n_node_face_tri*(k-1) + ij;
514 return 6 + n_tot_edge_node + n_tot_face_node + ijk;
518reference_element_H_ilat2loc_inod (
size_t order,
const point_basic<size_t>& ilat)
523 size_t i1 =
order - i;
524 size_t j1 =
order - j;
525 size_t k1 =
order - k;
526 size_t n_node_edge =
order-1;
527 size_t n_node_face = (
order-1)*(order-1);
530 if (i == 0)
return 0;
531 if (i == order)
return 1;
532 return 8 + 0*n_node_edge + (i - 1);
535 if (i == 0)
return 3;
536 if (i == order)
return 2;
537 return 8 + 2*n_node_edge + (i1 - 1);
540 return 8 + 3*n_node_edge + (j1 - 1);
543 return 8 + 1*n_node_edge + (j - 1);
546 size_t ji = n_node_edge*(i-1) + (j-1);
547 return 8 + 12*n_node_edge + 0*n_node_face + ji;
551 if (i == 0)
return 4;
552 if (i == order)
return 5;
553 return 8 + 8*n_node_edge + (i - 1);
556 if (i == 0)
return 7;
557 if (i == order)
return 6;
558 return 8 + 10*n_node_edge + (i1 - 1);
561 return 8 + 11*n_node_edge + (j1 - 1);
564 return 8 + 9*n_node_edge + (j - 1);
567 size_t ij = n_node_edge*(j-1) + (i-1);
568 return 8 + 12*n_node_edge + 3*n_node_face + ij;
572 return 8 + 4*n_node_edge + (k - 1);
575 return 8 + 5*n_node_edge + (k - 1);
578 size_t ik = n_node_edge*(k-1) + (i-1);
579 return 8 + 12*n_node_edge + 2*n_node_face + ik;
583 return 8 + 7*n_node_edge + (k - 1);
586 return 8 + 6*n_node_edge + (k - 1);
589 size_t i1k = n_node_edge*(k-1) + (i1-1);
590 return 8 + 12*n_node_edge + 5*n_node_face + i1k;
593 size_t kj = n_node_edge*(j-1) + (k-1);
594 return 8 + 12*n_node_edge + 1*n_node_face + kj;
597 size_t jk = n_node_edge*(k-1) + (j-1);
598 return 8 + 12*n_node_edge + 4*n_node_face + jk;
601 size_t n_node_bdry = 8 + 12*n_node_edge + 6*n_node_face;
602 size_t ijk = (
order-1)*((order-1)*(k-1) + (j-1)) + (i-1);
603 return n_node_bdry + ijk;
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type max_variant
static const variant_type p
static const variant_type T
static const variant_type P
static const variant_type t
constexpr size_t reference_element__H
constexpr size_t reference_element__e
constexpr size_t reference_element__P
static const size_t _first_variant_by_dimension[5]
static size_t _reference_element_n_face_by_variant[reference_element__max_variant]
constexpr size_t reference_element__T
constexpr size_t reference_element__p
constexpr size_t reference_element__t
constexpr size_t reference_element__q
constexpr size_t reference_element__max_variant
static size_t _reference_element_n_edge_by_variant[reference_element__max_variant]
static size_t _reference_element_dimension_by_variant[reference_element__max_variant]
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
#define _RHEOLEF_geo_element_auto_case(VARIANT)