Rheolef  7.1
an efficient C++ finite element environment
reference_element_aux.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
2 #define _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
23 // minimal auto-contained rheolef helpers for msh2geo
24 // used in:
25 // - reference_element.cc
26 // - msh2geo.cc
27 // in order to avoid code duplication
28 //
29 namespace rheolef {
30 
31 // used by msh2geo.cc that do not include reference_element.h
32 #ifndef _RHEOLEF_REFERENCE_ELEMENT_H
33 constexpr size_t reference_element__p = 0;
34 constexpr size_t reference_element__e = 1;
35 constexpr size_t reference_element__t = 2;
36 constexpr size_t reference_element__q = 3;
37 constexpr size_t reference_element__T = 4;
38 constexpr size_t reference_element__P = 5;
39 constexpr size_t reference_element__H = 6;
40 constexpr size_t reference_element__max_variant = 7;
41 #else // _RHEOLEF_REFERENCE_ELEMENT_H
50 #endif // _RHEOLEF_REFERENCE_ELEMENT_H
51 
55 
56 static
57 size_t
58 reference_element_variant (const char* name_table, char name)
59 {
60  for (size_t i = 0; i < reference_element__max_variant; i++) {
61  if (name_table [i] == name) return i;
62  }
64 }
65 static
66 size_t
67 reference_element_dimension_by_name(char element_name) {
68  switch (element_name) {
69  case 'p': return 0;
70  case 'e': return 1;
71  case 't':
72  case 'q': return 2;
73  case 'T':
74  case 'P':
75  case 'H': return 3;
76  default: error_macro("unexpected element name `"<<element_name<<"' (ascii="<<int(element_name)<<")");
77  return 0;
78  }
79 }
80 static
81 size_t
82 reference_element_dimension_by_variant (size_t variant)
83 {
85 }
86 static
87 size_t
88 reference_element_variant (char element_name)
89 {
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)<<")");
99  return 0;
100  }
101 }
102 static
103 size_t
104 reference_element_n_vertex (size_t variant)
105 {
106  switch (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;
115  }
116 }
117 static
118 size_t
119 reference_element_n_node (size_t variant, size_t order)
120 {
121  switch (variant) {
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;
130  }
131 }
132 static
133 void
134 reference_element_init_local_nnode_by_variant (
135  size_t order,
136  std::array<size_t,reference_element__max_variant>& sz)
137 {
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);
146 }
147 // -----------------------------------------------------------------------------
148 // reference_element::first_inod_by_dimension
149 // -----------------------------------------------------------------------------
150 static const size_t
157 };
158 static
159 size_t
160 reference_element_first_variant_by_dimension (size_t dim)
161 {
163 }
164 static
165 size_t
166 reference_element_last_variant_by_dimension (size_t dim)
167 {
169 }
170 // -----------------------------------------------------------------------------
171 // reference_element::first_inod_by_variant
172 // -----------------------------------------------------------------------------
173 static
174 size_t
175 reference_element_p_first_inod_by_variant (
176  size_t order,
177  size_t subgeo_variant)
178 {
179  if (subgeo_variant == reference_element__p) return 0;
180  return 1;
181 }
182 static
183 size_t
184 reference_element_e_first_inod_by_variant (
185  size_t order,
186  size_t subgeo_variant)
187 {
188  if (subgeo_variant == reference_element__p) return 0;
189  if (subgeo_variant == reference_element__e) return 2;
190  return order+1;
191 }
192 static
193 size_t
194 reference_element_t_first_inod_by_variant (
195  size_t order,
196  size_t subgeo_variant)
197 {
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;
202 }
203 static
204 size_t
205 reference_element_q_first_inod_by_variant (
206  size_t order,
207  size_t subgeo_variant)
208 {
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);
214 }
215 static
216 size_t
217 reference_element_T_first_inod_by_variant (
218  size_t order,
219  size_t subgeo_variant)
220 {
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;
227 }
228 static
229 size_t
230 reference_element_P_first_inod_by_variant (
231  size_t order,
232  size_t subgeo_variant)
233 {
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;
241 }
242 static
243 size_t
244 reference_element_H_first_inod_by_variant (
245  size_t order,
246  size_t subgeo_variant)
247 {
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);
256 }
257 static
258 size_t
259 reference_element_first_inod_by_variant (
260  size_t variant,
261  size_t order,
262  size_t subgeo_variant)
263 {
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);
267 
268  switch (variant) {
276  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
277  }
278 #undef _RHEOLEF_geo_element_auto_case
279 }
280 static
281 size_t
282 reference_element_last_inod_by_variant (
283  size_t variant,
284  size_t order,
285  size_t subgeo_variant)
286 {
287  return reference_element_first_inod_by_variant (variant, order, subgeo_variant+1);
288 }
289 // -----------------------------------------------------------------------------
290 // reference_element::ilat2loc_inod
291 // -----------------------------------------------------------------------------
292 static
293 size_t
294 reference_element_p_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
295 {
296  return 0;
297 }
298 static
299 size_t
300 reference_element_e_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
301 {
302  size_t i = ilat[0];
303  if (i == 0) return 0;
304  if (i == order) return 1;
305  return 2 + (i - 1);
306 }
307 static
308 size_t
309 reference_element_t_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
310 {
311  size_t i = ilat[0];
312  size_t j = ilat[1];
313  if (j == 0) { // horizontal edge [0,1]
314  if (i == 0) return 0;
315  if (i == order) return 1;
316  return 3 + 0*(order-1) + (i - 1);
317  }
318  if (i + j == order) { // oblique edge ]1,2]
319  if (j == order) return 2;
320  return 3 + 1*(order-1) + (j - 1);
321  }
322  size_t j1 = order - j;
323  if (i == 0) { // vertical edge ]2,0[
324  return 3 + 2*(order-1) + (j1 - 1);
325  }
326  // internal face ]0,1[x]0,2[: i, then j
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;
330 }
331 static
332 size_t
333 reference_element_q_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
334 {
335  size_t i = ilat[0];
336  size_t j = ilat[1];
337  if (j == 0) { // horizontal edge [0,1]
338  if (i == 0) return 0;
339  if (i == order) return 1;
340  return 4 + 0*(order-1) + (i - 1);
341  }
342  if (i == order) { // vertical edge ]1,2]
343  if (j == order) return 2;
344  return 4 + 1*(order-1) + (j - 1);
345  }
346  size_t i1 = order - i;
347  if (j == order) { // horizontal edge ]2,3]
348  if (i == 0) return 3;
349  return 4 + 2*(order-1) + (i1 - 1);
350  }
351  size_t j1 = order - j;
352  if (i == 0) { // vertical edge ]3,0[
353  return 4 + 3*(order-1) + (j1 - 1);
354  }
355  // internal face ]0,1[x]0,3[: i, then j
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;
359 }
360 static
361 size_t
362 reference_element_T_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
363 {
364  size_t i = ilat[0];
365  size_t j = ilat[1];
366  size_t k = ilat[2];
367  size_t n_tot_edge_node = 6*(order-1);
368  size_t n_face_node = (order-1)*(order-2)/2;
369  if (k == 0) {
370  // horizontal face 021 = [0,2]x[0,1]
371  if (j == 0) { // left edge [0,1]
372  if (i == 0) return 0;
373  if (i == order) return 1;
374  return 4 + 0*(order-1) + (i - 1);
375  }
376  if (i + j == order) { // oblique edge ]1,2[
377  if (j == order) return 2;
378  return 4 + 1*(order-1) + (j - 1);
379  }
380  size_t j1 = order - j;
381  if (i == 0) { // right edge ]0,2[
382  return 4 + 2*(order-1) + (j1 - 1);
383  }
384  // internal face 021 = ]0,2[x]0,1[: j, then i
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;
388  }
389  if (i == 0) {
390  // back face 032 =]0,3]x]0,2[
391  if (j == 0) { // vertical edge ]0,3]
392  if (k == order) return 3;
393  return 4 + 3*(order-1) + (k - 1);
394  }
395  if (j + k == order) { // oblique edge ]2,3[
396  return 4 + 5*(order-1) + (k - 1);
397  }
398  // internal face 032 =]0,3]x]0,2[, k, then j
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;
402  }
403  if (j == 0) {
404  // left face 013 = ]0,1[x]0,3[
405  if (i + k == order) { // oblique edge ]1,3[
406  return 4 + 4*(order-1) + (k - 1);
407  }
408  // internal face 013 = ]0,1[x]0,3[: i, then k
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;
412  }
413  if (i + j + k == order) {
414  // internal to oblique face 123 = ]1,2[x]1,3], j, then k
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;
418  }
419  // internal to volume: i, then j, then k
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)
428  + (i - 1);
429  return 4 + n_tot_edge_node + n_tot_face_node + ijk;
430 }
431 static
432 size_t
433 reference_element_P_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
434 {
435  size_t i = ilat[0];
436  size_t j = ilat[1];
437  size_t k = ilat[2];
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);
445  if (k == 0) {
446  // horizontal face 021 = [0,1]x[0,2]
447  if (j == 0) { // left edge [0,1]
448  if (i == 0) return 0;
449  if (i == order) return 1;
450  return 6 + 0*(order-1) + (i - 1);
451  }
452  if (i + j == order) { // oblique edge ]1,2[
453  if (j == order) return 2;
454  return 6 + 1*(order-1) + (j - 1);
455  }
456  if (i == 0) { // right edge ]0,2[
457  return 6 + 2*(order-1) + (j1 - 1);
458  }
459  // internal face 021 = ]0,2[x]0,1[: j, then i
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;
462  }
463  if (k == order) {
464  // horizontal face 354 = [3,4]x[3,5]
465  if (j == 0) { // left edge [3,4]
466  if (i == 0) return 3;
467  if (i == order) return 4;
468  return 6 + 6*(order-1) + (i - 1);
469  }
470  if (i + j == order) { // oblique edge ]4,5[
471  if (j == order) return 5;
472  return 6 + 7*(order-1) + (j - 1);
473  }
474  if (i == 0) { // right edge ]3,5[
475  return 6 + 8*(order-1) + (j1 - 1);
476  }
477  // internal face 354 = ]3,4[x]3,5[: i, then j
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;
480  }
481  if (j == 0) {
482  // left face 0143 = [0,1]x]0,3[
483  if (i == 0) {
484  // vertical edge ]0,3[
485  return 6 + 3*(order-1) + (k - 1);
486  }
487  if (i == order) {
488  // vertical edge ]1,4[
489  return 6 + 4*(order-1) + (k - 1);
490  }
491  // internal face 01243 = ]0,1[x]0,3[: i then k
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;
494  }
495  if (i == 0) {
496  // back face 0352 = ]0,2]x]0,3[
497  if (j == order) {
498  // vertical edge ]2,5[
499  return 6 + 5*(order-1) + (k - 1);
500  }
501  // internal face 0352 = ]0,2[x]0,3[: k then j
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;
504  }
505  if (i + j == order) {
506  // internal face 1254 = ]1,2[x]1,4[: j then k
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;
509  }
510  // internal volume ]0,1[x]0,2[x][0,3[: i, then j, then k
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;
515 }
516 static
517 size_t
518 reference_element_H_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
519 {
520  size_t i = ilat[0];
521  size_t j = ilat[1];
522  size_t k = ilat[2];
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);
528  if (k == 0) { // bottom face [0,3]x[0,1]
529  if (j == 0) { // bottom-left edge [0,1]
530  if (i == 0) return 0;
531  if (i == order) return 1;
532  return 8 + 0*n_node_edge + (i - 1);
533  }
534  if (j == order) { // bottom-right edge [3,2]
535  if (i == 0) return 3;
536  if (i == order) return 2;
537  return 8 + 2*n_node_edge + (i1 - 1);
538  }
539  if (i == 0) { // bottom-back edge ]3,0[
540  return 8 + 3*n_node_edge + (j1 - 1);
541  }
542  if (i == order) { // bottom-front edge ]1,2[
543  return 8 + 1*n_node_edge + (j - 1);
544  }
545  // internal face ]0,3[x]0,1[: j then i (TODO: not checked since Pk(H) only with k <= 2 yet
546  size_t ji = n_node_edge*(i-1) + (j-1);
547  return 8 + 12*n_node_edge + 0*n_node_face + ji;
548  }
549  if (k == order) { // top face [4,5]x[4,7]
550  if (j == 0) { // top-left edge [4,5]
551  if (i == 0) return 4;
552  if (i == order) return 5;
553  return 8 + 8*n_node_edge + (i - 1);
554  }
555  if (j == order) { // top-right edge [7,6]
556  if (i == 0) return 7;
557  if (i == order) return 6;
558  return 8 + 10*n_node_edge + (i1 - 1);
559  }
560  if (i == 0) { // top-back edge ]7,4[
561  return 8 + 11*n_node_edge + (j1 - 1);
562  }
563  if (i == order) { // top-front edge ]5,6[
564  return 8 + 9*n_node_edge + (j - 1);
565  }
566  // internal face ]4,5[x]4,7[: i then j
567  size_t ij = n_node_edge*(j-1) + (i-1);
568  return 8 + 12*n_node_edge + 3*n_node_face + ij;
569  }
570  if (j == 0) { // left face ]0,1[x[0,4]
571  if (i == 0) { // left-back edge [0,4]
572  return 8 + 4*n_node_edge + (k - 1);
573  }
574  if (i == order) { // left-front edge [1,5]
575  return 8 + 5*n_node_edge + (k - 1);
576  }
577  // internal face ]0,1[x]0,4[: i then k
578  size_t ik = n_node_edge*(k-1) + (i-1);
579  return 8 + 12*n_node_edge + 2*n_node_face + ik;
580  }
581  if (j == order) { // right face ]2,3[x[2,6]
582  if (i == 0) { // right-back edge [3,7]
583  return 8 + 7*n_node_edge + (k - 1);
584  }
585  if (i == order) { // right-front edge [2,6]
586  return 8 + 6*n_node_edge + (k - 1);
587  }
588  // internal face ]2,3[x]2,6[: i1 then k
589  size_t i1k = n_node_edge*(k-1) + (i1-1);
590  return 8 + 12*n_node_edge + 5*n_node_face + i1k;
591  }
592  if (i == 0) { // internal back face ]0,4[x[0,3]: k then j
593  size_t kj = n_node_edge*(j-1) + (k-1);
594  return 8 + 12*n_node_edge + 1*n_node_face + kj;
595  }
596  if (i == order) { // internal front face ]1,2[x[1,5]: j then k
597  size_t jk = n_node_edge*(k-1) + (j-1);
598  return 8 + 12*n_node_edge + 4*n_node_face + jk;
599  }
600  // internal volume: i then j then k
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;
604 }
605 
606 } // namespace rheolef
607 #endif // _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
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)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_geo_element_auto_case(NAME)
This file is part of Rheolef.
Definition: sphere.icc:25