rheolef  6.3
reference_element.cc
Go to the documentation of this file.
1 
2 #include "reference_element.h"
3 #include "rheolef/point.h"
4 
5 
6 namespace rheolef {
7 
8 namespace prism {
9 #include "prism.icc"
10 } // namespace prism
11 
12 namespace hexahedron {
13 #include "hexahedron.icc"
14 } // namespace hexahedron
15 
16 #include "reference_element_declare.cc"
17 
18 void
19 reference_element::set_name (char name)
20 {
21  _x = variant(name);
22  check_macro (_x != max_variant, "undefined reference element `" << name << "'");
23 }
24 Float
25 measure (reference_element hat_K) {
26  return hat_K_measure [hat_K.variant()];
27 }
28 reference_element::variant_type
29 reference_element::variant (char name)
30 {
31  for (size_type i = 0; i < max_variant; i++) {
32  if (_name [i] == name) return variant_type(i);
33  }
34  return variant_type(max_variant);
35 }
36 reference_element::variant_type
37 reference_element::variant (size_type n_vertex, size_type dim)
38 {
39  size_type variant = 0;
40  for (; variant < max_variant; variant++) {
41  if (_dimension[variant] == dim && _n_vertex[variant] == n_vertex) {
42  return variant_type(variant);
43  }
44  }
45  error_macro ("undefined "<<dim<<"d reference element with "<<n_vertex << " vertices");
46 }
47 reference_element::size_type
48 reference_element::n_sub_edge (variant_type variant)
49 {
50  switch (variant) {
51  case reference_element::p: return 0;
52  case reference_element::e: return 0;
53  case reference_element::t: return 3;
54  case reference_element::q: return 4;
55  case reference_element::T: return 6;
56  case reference_element::P: return 9;
57  case reference_element::H: return 12;
58  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
59  }
60 }
61 reference_element::size_type
62 reference_element::n_sub_face (variant_type variant)
63 {
64  switch (variant) {
65  case reference_element::p: return 0;
66  case reference_element::e: return 0;
67  case reference_element::t: return 0;
68  case reference_element::q: return 0;
69  case reference_element::T: return 4;
70  case reference_element::P: return 5;
71  case reference_element::H: return 6;
72  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
73  }
74 }
75 reference_element::size_type
76 reference_element::n_subgeo (variant_type variant, size_type subgeo_dim)
77 {
78 #define _RHEOLEF_reference_element_case(VARIANT) \
79  case reference_element::VARIANT: \
80  return reference_element_##VARIANT::n_subgeo (subgeo_dim);
81 
82  switch (variant) {
83  _RHEOLEF_reference_element_case(p)
84  _RHEOLEF_reference_element_case(e)
85  _RHEOLEF_reference_element_case(t)
86  _RHEOLEF_reference_element_case(q)
87  _RHEOLEF_reference_element_case(T)
88  _RHEOLEF_reference_element_case(P)
89  _RHEOLEF_reference_element_case(H)
90  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
91  }
92 #undef _RHEOLEF_reference_element_case
93 }
94 reference_element::size_type
95 reference_element::subgeo_n_node (
96  variant_type variant,
97  size_type order,
98  size_type subgeo_dim,
99  size_type loc_isid)
100 {
101 #define _RHEOLEF_geo_element_auto_case(VARIANT) \
102  case reference_element::VARIANT: \
103  return reference_element_##VARIANT::subgeo_n_node (order, subgeo_dim, loc_isid);
104 
105  switch (variant) {
106  _RHEOLEF_geo_element_auto_case(p)
107  _RHEOLEF_geo_element_auto_case(e)
108  _RHEOLEF_geo_element_auto_case(t)
109  _RHEOLEF_geo_element_auto_case(q)
110  _RHEOLEF_geo_element_auto_case(T)
111  _RHEOLEF_geo_element_auto_case(P)
112  _RHEOLEF_geo_element_auto_case(H)
113  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
114  }
115 #undef _RHEOLEF_geo_element_auto_case
116 }
117 reference_element::size_type
118 reference_element::subgeo_local_node (
119  variant_type variant,
120  size_type order,
121  size_type subgeo_dim,
122  size_type loc_isid,
123  size_type loc_jsidnod)
124 {
125 #define _RHEOLEF_geo_element_auto_case(VARIANT) \
126  case reference_element::VARIANT: \
127  return reference_element_##VARIANT::subgeo_local_node (order, subgeo_dim, loc_isid, loc_jsidnod);
128 
129  switch (variant) {
130  _RHEOLEF_geo_element_auto_case(p)
131  _RHEOLEF_geo_element_auto_case(e)
132  _RHEOLEF_geo_element_auto_case(t)
133  _RHEOLEF_geo_element_auto_case(q)
134  _RHEOLEF_geo_element_auto_case(T)
135  _RHEOLEF_geo_element_auto_case(P)
136  _RHEOLEF_geo_element_auto_case(H)
137  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
138  }
139 #undef _RHEOLEF_geo_element_auto_case
140 }
141 reference_element::size_type
142 reference_element::first_inod_by_variant (
143  variant_type variant,
144  size_type order,
145  variant_type subgeo_variant)
146 {
147 #define _RHEOLEF_geo_element_auto_case(VARIANT) \
148  case reference_element::VARIANT: \
149  return reference_element_##VARIANT::first_inod_by_variant (order, subgeo_variant);
150 
151  switch (variant) {
152  _RHEOLEF_geo_element_auto_case(p)
153  _RHEOLEF_geo_element_auto_case(e)
154  _RHEOLEF_geo_element_auto_case(t)
155  _RHEOLEF_geo_element_auto_case(q)
156  _RHEOLEF_geo_element_auto_case(T)
157  _RHEOLEF_geo_element_auto_case(P)
158  _RHEOLEF_geo_element_auto_case(H)
159  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
160  }
161 #undef _RHEOLEF_geo_element_auto_case
162 }
163 reference_element::size_type
164 reference_element::n_node (variant_type variant, size_type order)
165 {
166  switch (variant) {
167  case reference_element::p: return 1;
168  case reference_element::e: return order+1;
169  case reference_element::t: return (order+1)*(order+2)/2;
170  case reference_element::q: return (order+1)*(order+1);
171  case reference_element::T: return (order+1)*(order+2)*(order+3)/6;
172  case reference_element::P: return (order+1)*(order+1)*(order+2)/2;
173  case reference_element::H: return (order+1)*(order+1)*(order+1);
174  default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
175  }
176 }
177 void
178 reference_element::init_local_nnode_by_variant (
179  size_type order,
180  boost::array<size_type,reference_element::max_variant>& sz)
181 {
182  sz [reference_element::p] = 1;
183  sz [reference_element::e] = order-1;
184  sz [reference_element::t] = (order-1)*(order-2)/2;
185  sz [reference_element::q] = (order-1)*(order-1);
186  sz [reference_element::T] = (order-1)*(order-2)*(order-3)/6;
187  sz [reference_element::P] = (order-1)*(order-1)*(order-2)/2;
188  sz [reference_element::H] = (order-1)*(order-1)*(order-1);
189 }
190 Float
191 reference_element::side_measure (size_type loc_isid) const
192 {
193 #define _RHEOLEF_geo_element_auto_case(VARIANT) \
194  case reference_element::VARIANT: \
195  return reference_element_##VARIANT::side_measure (loc_isid);
196 
197  switch (variant()) {
198  _RHEOLEF_geo_element_auto_case(p)
199  _RHEOLEF_geo_element_auto_case(e)
200  _RHEOLEF_geo_element_auto_case(t)
201  _RHEOLEF_geo_element_auto_case(q)
202  _RHEOLEF_geo_element_auto_case(T)
203  _RHEOLEF_geo_element_auto_case(P)
204  _RHEOLEF_geo_element_auto_case(H)
205  default: error_macro ("unexpected element variant `"<<variant()<<"'"); return 0;
206  }
207 #undef _RHEOLEF_geo_element_auto_case
208 }
209 void
210 reference_element::side_normal (size_type loc_isid, point_basic<Float>& hat_n) const
211 {
212 #define _RHEOLEF_geo_element_auto_case(VARIANT) \
213  case reference_element::VARIANT: \
214  reference_element_##VARIANT::side_normal (loc_isid, hat_n); break;
215 
216  switch (variant()) {
217  _RHEOLEF_geo_element_auto_case(p)
218  _RHEOLEF_geo_element_auto_case(e)
219  _RHEOLEF_geo_element_auto_case(t)
220  _RHEOLEF_geo_element_auto_case(q)
221  _RHEOLEF_geo_element_auto_case(T)
222  _RHEOLEF_geo_element_auto_case(P)
223  _RHEOLEF_geo_element_auto_case(H)
224  default: error_macro ("unexpected element variant `"<<variant()<<"'");
225  }
226 #undef _RHEOLEF_geo_element_auto_case
227 }
228 reference_element_p::size_type
229 reference_element_p::n_subgeo (size_type side_dim)
230 {
231  return (side_dim == 0) ? 1 : 0;
232 }
233 reference_element_p::size_type
234 reference_element_p::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
235 {
236  return (side_dim == 0) ? 1 : 0;
237 }
238 reference_element_p::size_type
239 reference_element_p::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
240 {
241  return 0;
242 }
243 // edge 0d-lattice for high order elements, i <= 0
244 reference_element_p::size_type
245 reference_element_p::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
246 {
247  return 0;
248 }
249 reference_element_p::size_type
250 reference_element_p::first_inod_by_variant (
251  size_type order,
252  size_type subgeo_variant)
253 {
254  if (subgeo_variant == reference_element::p) return 0;
255  return 1;
256 }
257 Float
258 reference_element_p::side_measure (size_type loc_isid)
259 {
260  return 0;
261 }
262 void
263 reference_element_p::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
264 {
265 }
266 reference_element_e::size_type
267 reference_element_e::n_subgeo (size_type side_dim)
268 {
269  switch (side_dim) {
270  case 0: return 2;
271  case 1: return 1;
272  default: return 0;
273  }
274 }
275 reference_element_e::size_type
276 reference_element_e::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
277 {
278  switch (side_dim) {
279  case 0: return 1;
280  case 1: return order+1;
281  default: return 0;
282  }
283 }
284 reference_element_e::size_type
285 reference_element_e::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
286 {
287  switch (side_dim) {
288  case 0: return loc_isid;
289  case 1: return loc_jsidnod;
290  default: return 0;
291  }
292 }
293 // edge 1d-lattice (i) for high order elements, i <= order
294 reference_element_e::size_type
295 reference_element_e::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
296 {
297  size_type i = ilat[0];
298  if (i == 0) return 0;
299  if (i == order) return 1;
300  return 2 + (i - 1);
301 }
302 reference_element_e::size_type
303 reference_element_e::first_inod_by_variant (
304  size_type order,
305  size_type subgeo_variant)
306 {
307  if (subgeo_variant == reference_element::p) return 0;
308  if (subgeo_variant == reference_element::e) return 2;
309  return order+1;
310 }
311 Float
312 reference_element_e::side_measure (size_type loc_isid)
313 {
314  return 1;
315 }
316 void
317 reference_element_e::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
318 {
319  hat_n[0] = (loc_isid == 0) ? -1 : 1;
320 }
321 reference_element_t::size_type
322 reference_element_t::n_subgeo (size_type side_dim)
323 {
324  switch (side_dim) {
325  case 0: return 3;
326  case 1: return 3;
327  case 2: return 1;
328  default: return 0;
329  }
330 }
331 reference_element_t::size_type
332 reference_element_t::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
333 {
334  switch (side_dim) {
335  case 0: return 1;
336  case 1: return order+1;
337  case 2: return ((order+1)*(order+2))/2;
338  default: return 0;
339  }
340 }
341 reference_element_t::size_type
342 reference_element_t::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
343 {
344  switch (side_dim) {
345  case 0: return loc_isid;
346  case 1: if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 3; // edge-node is a vertex
347  else return (order-1)*loc_isid + loc_jsidnod + 1; // edge-node is edge-internal
348  case 2: return loc_jsidnod;
349  default: return 0;
350  }
351 }
352 // triangle lattice (i,j) for high order elements, i+j <= order
353 reference_element_t::size_type
354 reference_element_t::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
355 {
356  size_type i = ilat[0];
357  size_type j = ilat[1];
358  if (j == 0) { // horizontal edge [0,1]
359  if (i == 0) return 0;
360  if (i == order) return 1;
361  return 3 + 0*(order-1) + (i - 1);
362  }
363  if (i + j == order) { // oblique edge ]1,2]
364  if (j == order) return 2;
365  return 3 + 1*(order-1) + (j - 1);
366  }
367  size_type j1 = order - j;
368  if (i == 0) { // vertical edge ]2,0[
369  return 3 + 2*(order-1) + (j1 - 1);
370  }
371  size_type n_face_node = (order-1)*(order-2)/2;
372  size_type ij = (n_face_node - (j1-1)*j1/2) + (i - 1);
373  return 3 + 3*(order-1) + ij;
374 }
375 reference_element_t::size_type
376 reference_element_t::first_inod_by_variant (
377  size_type order,
378  size_type subgeo_variant)
379 {
380  if (subgeo_variant == reference_element::p) return 0;
381  if (subgeo_variant == reference_element::e) return 3;
382  if (subgeo_variant == reference_element::t) return 3 + 3*(order-1);
383  return (order+1)*(order+2)/2;
384 }
385 Float
386 reference_element_t::side_measure (size_type loc_isid)
387 {
388  return (loc_isid != 1) ? 1 : sqrt(Float(2.));
389 }
390 void
391 reference_element_t::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
392 {
393  switch (loc_isid) {
394  case 0: hat_n = point_basic<Float>( 0,-1); break;
395  case 1: hat_n = point_basic<Float>( 1, 1)/sqrt(Float(2)); break;
396  case 2:
397  default: hat_n = point_basic<Float>(-1,0); break;
398  }
399 }
400 reference_element_q::size_type
401 reference_element_q::n_subgeo (size_type side_dim)
402 {
403  switch (side_dim) {
404  case 0: return 4;
405  case 1: return 4;
406  case 2: return 1;
407  default: return 0;
408  }
409 }
410 reference_element_q::size_type
411 reference_element_q::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
412 {
413  switch (side_dim) {
414  case 0: return 1;
415  case 1: return order+1;
416  case 2: return (order+1)*(order+1);
417  default: return 0;
418  }
419 }
420 reference_element_q::size_type
421 reference_element_q::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
422 {
423  switch (side_dim) {
424  case 0: return loc_isid;
425  case 1: if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 4; // edge-node is a vertex
426  else return (order-1)*loc_isid + loc_jsidnod + 2; // edge-node is edge-internal
427  case 2: return loc_jsidnod;
428  default: return 0;
429  }
430 }
431 // quadrangle lattice (i,j) for high order elements, 0 <= i,j <= order
432 reference_element_q::size_type
433 reference_element_q::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
434 {
435  size_type i = ilat[0];
436  size_type j = ilat[1];
437  if (j == 0) { // horizontal edge [0,1]
438  if (i == 0) return 0;
439  if (i == order) return 1;
440  return 4 + 0*(order-1) + (i - 1);
441  }
442  if (i == order) { // vertical edge ]1,2]
443  if (j == order) return 2;
444  return 4 + 1*(order-1) + (j - 1);
445  }
446  size_type i1 = order - i;
447  if (j == order) { // horizontal edge ]2,3]
448  if (i == 0) return 3;
449  return 4 + 2*(order-1) + (i1 - 1);
450  }
451  size_type j1 = order - j;
452  if (i == 0) { // vertical edge ]3,0[
453  return 4 + 3*(order-1) + (j1 - 1);
454  }
455  size_type ij = (order-1)*(j-1) + (i-1);
456  size_type n_bdry_node = 4 + 4*(order-1);
457  return n_bdry_node + ij;
458 }
459 reference_element_q::size_type
460 reference_element_q::first_inod_by_variant (
461  size_type order,
462  size_type subgeo_variant)
463 {
464  if (subgeo_variant == reference_element::p) return 0;
465  if (subgeo_variant == reference_element::e) return 4;
466  if (subgeo_variant == reference_element::t) return 4 + 4*(order-1);
467  if (subgeo_variant == reference_element::q) return 4 + 4*(order-1);
468  return (order+1)*(order+1);
469 }
470 Float
471 reference_element_q::side_measure (size_type loc_isid)
472 {
473  return 2;
474 }
475 void
476 reference_element_q::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
477 {
478  switch (loc_isid) {
479  case 0: hat_n = point_basic<Float>( 0,-1); break;
480  case 1: hat_n = point_basic<Float>( 1, 0); break;
481  case 2: hat_n = point_basic<Float>( 0, 1); break;
482  case 3:
483  default: hat_n = point_basic<Float>(-1, 0); break;
484  }
485 }
486 reference_element_T::size_type
487 reference_element_T::n_subgeo (size_type side_dim)
488 {
489  switch (side_dim) {
490  case 0: return 4;
491  case 1: return 6;
492  case 2: return 4;
493  case 3: return 1;
494  default: return 0;
495  }
496 }
497 reference_element_T::size_type
498 reference_element_T::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
499 {
500  switch (side_dim) {
501  case 0: return 1;
502  case 1: return order+1;
503  case 2: return (order+1)*(order+2)/2;
504  case 3: return ((order+1)*(order+2)*(order+3))/6;
505  default: return 0;
506  }
507 }
508 reference_element_T::size_type
509 reference_element_T::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
510 {
511  switch (side_dim) {
512  case 0: return loc_isid;
513  case 1: if (loc_jsidnod < 2) { // edge-node is a vertex
514  if (loc_isid < 3) return (loc_isid + loc_jsidnod) % 3;
515  else return loc_jsidnod == 0 ? loc_isid - 3 : 3;
516  } else { // edge-node is internal to the edge
517  return loc_isid*(order-1) + loc_jsidnod + 2;
518  }
519  case 2: {
520  if (loc_jsidnod < 3) { // face-node is a vertex
521  if (loc_isid == 3) return loc_jsidnod + 1;
522  if (loc_jsidnod == 0) return 0;
523  if (loc_jsidnod == 2) return loc_isid + 1;
524  return ((loc_isid + 1) % 3) + 1;
525  }
526  size_type last_edge_node_iloc = 3 + 3*(order-1);
527  if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
528 #ifdef TO_CLEAN
529  extern const size_type geo_element_T_fac2edg_idx [4][3];
530  extern const int geo_element_T_fac2edg_orient [4][3];
531 #endif // TO_CLEAN
532  size_type order1 = order - 1; // avoid div by zero compiler error
533  size_type loc_jedg = (loc_jsidnod-3) / order1;
534  size_type loc_kedg = (loc_jsidnod-3) % order1;
535  if (geo_element_T_fac2edg_orient [loc_isid][loc_jedg] < 0) {
536  loc_kedg = order - loc_kedg - 2;
537  }
538  return (order-1)*geo_element_T_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 4;
539  }
540  size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
541  return 4 + 6*(order-1) + loc_isid*(order-1)*(order-2)/2 + ij_loc;
542  }
543  case 3: return loc_jsidnod;
544  default: return 0;
545  }
546 }
547 // tetrahedron lattice (i,j,k) for high order elements, i+j+k <= order
548 reference_element_T::size_type
549 reference_element_T::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
550 {
551  size_type i = ilat[0];
552  size_type j = ilat[1];
553  size_type k = ilat[2];
554  size_type n_tot_edge_node = 6*(order-1);
555  size_type n_face_node = (order-1)*(order-2)/2;
556  if (k == 0) {
557  if (j == 0) { // left edge [0,1]
558  if (i == 0) return 0;
559  if (i == order) return 1;
560  return 4 + 0*(order-1) + (i - 1);
561  }
562  if (i + j == order) { // oblique edge ]1,2[
563  if (j == order) return 2;
564  return 4 + 1*(order-1) + (j - 1);
565  }
566  size_type j1 = order - j;
567  if (i == 0) { // right edge ]0,2[
568  return 4 + 2*(order-1) + (j1 - 1);
569  }
570  size_type i1 = order - i;
571  size_type ji = (n_face_node - i1*(i1-1)/2) + (j - 1);
572  return 4 + n_tot_edge_node + 0*n_face_node + ji;
573  }
574  if (i == 0) {
575  if (j == 0) { // vertical edge ]0,3]
576  if (k == order) return 3;
577  return 4 + 3*(order-1) + (k - 1);
578  }
579  if (j + k == order) { // oblique edge ]2,3[
580  return 4 + 5*(order-1) + (k - 1);
581  }
582  size_type j1 = order - j;
583  size_type kj = (n_face_node - j1*(j1-1)/2) + (k - 1);
584  return 4 + n_tot_edge_node + 1*n_face_node + kj;
585  }
586  if (j == 0) {
587  if (i + k == order) { // oblique edge ]1,3[
588  return 4 + 4*(order-1) + (k - 1);
589  }
590  size_type k1 = order - k;
591  size_type ik = (n_face_node - k1*(k1-1)/2) + (i - 1);
592  return 4 + n_tot_edge_node + 2*n_face_node + ik;
593  }
594  if (i + j + k == order) {
595  size_type k1 = order - k;
596  size_type jk = (n_face_node - k1*(k1-1)/2) + (j - 1);
597  return 4 + n_tot_edge_node + 3*n_face_node + jk;
598  }
599  size_type n_tot_face_node = 4*n_face_node;
600  size_type n_tot_node = (order+1)*(order+2)*(order+3)/6;
601  size_type n_volume_node = (order-1)*(order-2)*(order-3)/6;
602  size_type k1 = order - k;
603  size_type j1 = order - j - k;
604  size_type n_level_k_node = (k1-1)*(k1-2)/2;
605  size_type ijk = (n_volume_node - k1*(k1-1)*(k1-2)/6)
606  + (n_level_k_node - j1*(j1-1)/2)
607  + (i - 1);
608  return 4 + n_tot_edge_node + n_tot_face_node + ijk;
609 }
610 reference_element_T::size_type
611 reference_element_T::first_inod_by_variant (
612  size_type order,
613  size_type subgeo_variant)
614 {
615  if (subgeo_variant == reference_element::p) return 0;
616  if (subgeo_variant == reference_element::e) return 4;
617  if (subgeo_variant == reference_element::t) return 4 + 6*(order-1);
618  if (subgeo_variant == reference_element::q) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
619  if (subgeo_variant == reference_element::T) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
620  return (order+1)*(order+2)*(order+3)/6;
621 }
622 reference_element_T::size_type
623 reference_element_T::face2edge (size_type loc_iface, size_type loc_iface_jedg)
624 {
625  return geo_element_T_fac2edg_idx [loc_iface][loc_iface_jedg];
626 }
627 int
628 reference_element_T::face2edge_orient (size_type loc_iface, size_type loc_iface_jedg)
629 {
630  return geo_element_T_fac2edg_orient [loc_iface][loc_iface_jedg];
631 }
632 Float
633 reference_element_T::side_measure (size_type loc_isid)
634 {
635  return (loc_isid != 3) ? 0.5 : sqrt(Float(3.))/2.;
636 }
637 void
638 reference_element_T::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
639 {
640  switch (loc_isid) {
641  case 0: hat_n = point_basic<Float>( 0, 0,-1); break;
642  case 1: hat_n = point_basic<Float>(-1, 0, 0); break;
643  case 2: hat_n = point_basic<Float>( 0,-1, 0); break;
644  case 3:
645  default: hat_n = point_basic<Float>( 1, 1, 1)/sqrt(Float(3)); break;
646  }
647 }
648 reference_element_P::size_type
649 reference_element_P::n_subgeo (size_type side_dim)
650 {
651  switch (side_dim) {
652  case 0: return 6;
653  case 1: return 9;
654  case 2: return 5;
655  case 3: return 1;
656  default: return 0;
657  }
658 }
659 reference_element_P::size_type
660 reference_element_P::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
661 {
662  switch (side_dim) {
663  case 0: return 1;
664  case 1: return order+1;
665  case 2: return (loc_isid < 2) ? (order+1)*(order+2)/2 : (order+1)*(order+1);
666  case 3: return ((order+1)*(order+1)*(order+2))/2;
667  default: return 0;
668  }
669 }
670 reference_element_P::size_type
671 reference_element_P::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
672 {
673  switch (side_dim) {
674  case 0: return loc_isid;
675  case 1: if (loc_jsidnod < 2) { // edge-node is a vertex
676  return prism::edge [loc_isid] [loc_jsidnod];
677  } else { // edge-node is internal to the edge
678  return 6 + loc_isid*(order-1) + (loc_jsidnod - 2);
679  }
680  case 2: {
681  size_type n_vert_on_side = (loc_isid < 2) ? 3 : 4;
682  if (loc_jsidnod < n_vert_on_side) { // face-node is a vertex
683  return prism::face [loc_isid] [loc_jsidnod];
684  }
685  size_type last_edge_node_iloc = n_vert_on_side + n_vert_on_side*(order-1);
686  if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
687 #ifdef TO_CLEAN
688  extern const size_type geo_element_P_fac2edg_idx [5][4];
689  extern const int geo_element_P_fac2edg_orient [5][4];
690 #endif // TO_CLEAN
691  size_type order1 = order - 1; // avoid div by zero compiler error
692  size_type loc_jedg = (loc_jsidnod - n_vert_on_side) / order1;
693  size_type loc_kedg = (loc_jsidnod - n_vert_on_side) % order1;
694  if (geo_element_P_fac2edg_orient [loc_isid][loc_jedg] < 0) {
695  loc_kedg = order - loc_kedg - 2;
696  }
697  return 6 + (order-1)*geo_element_P_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg;
698  }
699  size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
700  size_type shift_prev_faces;
701  if (loc_isid < 2) {
702  shift_prev_faces = loc_isid*(order-1)*(order-2)/2;
703  } else {
704  shift_prev_faces = 2*(order-1)*(order-2)/2 + (loc_isid-2)*(order-1)*(order-1);
705  }
706  return 6 + 9*(order-1) + shift_prev_faces + ij_loc;
707  }
708  case 3: return loc_jsidnod;
709  default: return 0;
710  }
711 }
712 reference_element_P::size_type
713 reference_element_P::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
714 {
715  size_type i = ilat[0];
716  size_type j = ilat[1];
717  size_type k = ilat[2];
718  size_type i1 = order - i;
719  size_type j1 = order - j;
720  size_type k1 = order - k;
721  size_type n_node_edge = order-1;
722  size_type n_tot_edge_node = 9*(order-1);
723  size_type n_node_face_tri = (order-1)*(order-2)/2;
724  size_type n_node_face_qua = (order-1)*(order-1);
725  if (k == 0) {
726  if (j == 0) { // left edge [0,1]
727  if (i == 0) return 0;
728  if (i == order) return 1;
729  return 6 + 0*(order-1) + (i - 1);
730  }
731  if (i + j == order) { // oblique edge ]1,2[
732  if (j == order) return 2;
733  return 6 + 1*(order-1) + (j - 1);
734  }
735  if (i == 0) { // right edge ]0,2[
736  return 6 + 2*(order-1) + (j1 - 1);
737  }
738  size_type ji = (n_node_face_tri - i1*(i1-1)/2) + (j - 1);
739  return 6 + n_tot_edge_node + 0*n_node_face_tri + ji;
740  }
741  if (k == order) {
742  if (j == 0) { // left edge [3,4]
743  if (i == 0) return 3;
744  if (i == order) return 4;
745  return 6 + 6*(order-1) + (i - 1);
746  }
747  if (i + j == order) { // oblique edge ]4,5[
748  if (j == order) return 5;
749  return 6 + 7*(order-1) + (j - 1);
750  }
751  if (i == 0) { // right edge ]3,5[
752  return 6 + 8*(order-1) + (j1 - 1);
753  }
754  size_type ij = (n_node_face_tri - j1*(j1-1)/2) + (i - 1);
755  return 6 + n_tot_edge_node + 1*n_node_face_tri + ij;
756  }
757  if (j == 0) {
758  if (i == 0) {
759  return 6 + 3*(order-1) + (k - 1);
760  }
761  if (i == order) {
762  return 6 + 4*(order-1) + (k - 1);
763  }
764  size_type ik = n_node_edge*(k-1) + (i-1);
765  return 6 + n_tot_edge_node + 2*n_node_face_tri + 0*n_node_face_qua + ik;
766  }
767  if (i == 0) {
768  if (j == order) {
769  return 6 + 5*(order-1) + (k - 1);
770  }
771  size_type kj = n_node_edge*(j-1) + (k-1);
772  return 6 + n_tot_edge_node + 2*n_node_face_tri + 2*n_node_face_qua + kj;
773  }
774  if (i + j == order) {
775  size_type jk = n_node_edge*(k-1) + (j-1);
776  return 6 + n_tot_edge_node + 2*n_node_face_tri + 1*n_node_face_qua + jk;
777  }
778  size_type n_tot_face_node = 2*n_node_face_tri + 3*n_node_face_qua;
779  size_type ij = (n_node_face_tri - (j1-1)*j1/2) + (i - 1);
780  size_type ijk = n_node_face_tri*(k-1) + ij;
781  return 6 + n_tot_edge_node + n_tot_face_node + ijk;
782 }
783 reference_element_P::size_type
784 reference_element_P::first_inod_by_variant (
785  size_type order,
786  size_type subgeo_variant)
787 {
788  if (subgeo_variant == reference_element::p) return 0;
789  if (subgeo_variant == reference_element::e) return 6;
790  if (subgeo_variant == reference_element::t) return 6 + 9*(order-1);
791  if (subgeo_variant == reference_element::q) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2);
792  if (subgeo_variant == reference_element::T) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
793  if (subgeo_variant == reference_element::P) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
794  return ((order+1)*(order+1)*(order+2))/2;
795 }
796 Float
797 reference_element_P::side_measure (size_type loc_isid)
798 {
799  fatal_macro ("side_measure(3d): not yet");
800  return 0;
801 }
802 void
803 reference_element_P::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
804 {
805  fatal_macro ("side_normal(3d): not yet");
806 }
807 
808 reference_element_H::size_type
809 reference_element_H::n_subgeo (size_type side_dim)
810 {
811  switch (side_dim) {
812  case 0: return 8;
813  case 1: return 12;
814  case 2: return 6;
815  case 3: return 1;
816  default: return 0;
817  }
818 }
819 reference_element_H::size_type
820 reference_element_H::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
821 {
822  switch (side_dim) {
823  case 0: return 1;
824  case 1: return order+1;
825  case 2: return (order+1)*(order+1);
826  case 3: return (order+1)*(order+1)*(order+1);
827  default: return 0;
828  }
829 }
830 reference_element_H::size_type
831 reference_element_H::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
832 {
833  switch (side_dim) {
834  case 0: return loc_isid;
835  case 1: if (loc_jsidnod < 2) { // edge-node is a vertex
836  if (loc_isid < 4) return (loc_isid + loc_jsidnod) % 4;
837  else if (loc_isid < 8) return loc_isid - 4 + loc_jsidnod*4;
838  else return loc_jsidnod == 0 ? loc_isid - 4 : (loc_isid - 7)%4 + 4;
839  } else { // edge-node is internal to the edge
840  return loc_isid*(order-1) + loc_jsidnod + 6; // TODO : BUG? +8 a la place de +6 ??
841  }
842  case 2: {
843  if (loc_jsidnod < 4) { // face-node is a vertex
844  return hexahedron::face [loc_isid][loc_jsidnod];
845  }
846  size_type last_edge_node_iloc = 4 + 4*(order-1);
847  if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
848 #ifdef TO_CLEAN
849  extern const size_t geo_element_H_fac2edg_idx [6][4];
850  extern const int geo_element_H_fac2edg_orient [6][4];
851 #endif // TO_CLEAN
852  size_type order1 = order - 1; // avoid div by zero compiler error
853  size_type loc_jedg = (loc_jsidnod-4) / order1;
854  size_type loc_kedg = (loc_jsidnod-4) % order1;
855  if (geo_element_H_fac2edg_orient [loc_isid][loc_jedg] < 0) {
856  loc_kedg = order - loc_kedg - 2;
857  }
858  return (order-1)*geo_element_H_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 8;
859  }
860  return 8 + 12*(order-1) + (order-1)*(order-1)*loc_isid + (loc_jsidnod - last_edge_node_iloc);
861  }
862  case 3: return loc_jsidnod;
863  default: return 0;
864  }
865 }
866 // edge 0d-lattice for high order elements, i <= 0
867 typedef reference_element_H::size_type size_type;
868 static
869 reference_element_H::size_type
870 H_ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
871 {
872  size_type i = ilat[0];
873  size_type j = ilat[1];
874  size_type k = ilat[2];
875  size_type i1 = order - i;
876  size_type j1 = order - j;
877  size_type k1 = order - k;
878  size_type n_node_edge = order-1;
879  size_type n_node_face = (order-1)*(order-1);
880  if (k == 0) { // bottom face [0,3]x[0,1]
881  if (j == 0) { // bottom-left edge [0,1]
882  if (i == 0) return 0;
883  if (i == order) return 1;
884  return 8 + 0*n_node_edge + (i - 1);
885  }
886  if (j == order) { // bottom-right edge [3,2]
887  if (i == 0) return 3;
888  if (i == order) return 2;
889  return 8 + 2*n_node_edge + (i1 - 1);
890  }
891  if (i == 0) { // bottom-back edge ]3,0[
892  return 8 + 3*n_node_edge + (j1 - 1);
893  }
894  if (i == order) { // bottom-front edge ]1,2[
895  return 8 + 1*n_node_edge + (j - 1);
896  }
897  // internal face ]0,3[x]0,1[: j then i (TODO: not checked since Pk(H) only with k <= 2 yet
898  size_type ji = n_node_edge*(i-1) + (j-1);
899  return 8 + 12*n_node_edge + 0*n_node_face + ji;
900  }
901  if (k == order) { // top face [4,5]x[4,7]
902  if (j == 0) { // top-left edge [4,5]
903  if (i == 0) return 4;
904  if (i == order) return 5;
905  return 8 + 8*n_node_edge + (i - 1);
906  }
907  if (j == order) { // top-right edge [7,6]
908  if (i == 0) return 7;
909  if (i == order) return 6;
910  return 8 + 10*n_node_edge + (i1 - 1);
911  }
912  if (i == 0) { // top-back edge ]7,4[
913  return 8 + 11*n_node_edge + (j1 - 1);
914  }
915  if (i == order) { // top-front edge ]5,6[
916  return 8 + 9*n_node_edge + (j - 1);
917  }
918  size_type ij = n_node_edge*(j-1) + (i-1);
919  return 8 + 12*n_node_edge + 3*n_node_face + ij;
920  }
921  if (j == 0) { // left face ]0,1[x[0,4]
922  if (i == 0) { // left-back edge [0,4]
923  return 8 + 4*n_node_edge + (k - 1);
924  }
925  if (i == order) { // left-front edge [1,5]
926  return 8 + 5*n_node_edge + (k - 1);
927  }
928  size_type ik = n_node_edge*(k-1) + (i-1);
929  return 8 + 12*n_node_edge + 2*n_node_face + ik;
930  }
931  if (j == order) { // right face ]2,3[x[2,6]
932  if (i == 0) { // right-back edge [3,7]
933  return 8 + 7*n_node_edge + (k - 1);
934  }
935  if (i == order) { // right-front edge [2,6]
936  return 8 + 6*n_node_edge + (k - 1);
937  }
938  size_type i1k = n_node_edge*(k-1) + (i1-1);
939  return 8 + 12*n_node_edge + 5*n_node_face + i1k;
940  }
941  if (i == 0) { // internal back face ]0,4[x[0,3]: k then j
942  size_type kj = n_node_edge*(j-1) + (k-1);
943  return 8 + 12*n_node_edge + 1*n_node_face + kj;
944  }
945  if (i == order) { // internal front face ]1,2[x[1,5]: j then k
946  size_type jk = n_node_edge*(k-1) + (j-1);
947  return 8 + 12*n_node_edge + 4*n_node_face + jk;
948  }
949  size_type n_node_bdry = 8 + 12*n_node_edge + 6*n_node_face;
950  size_type ijk = (order-1)*((order-1)*(k-1) + (j-1)) + (i-1);
951  return n_node_bdry + ijk;
952 }
953 reference_element_H::size_type
954 reference_element_H::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
955 {
956  size_type loc_inod = H_ilat2loc_inod (order, ilat);
957  return H_ilat2loc_inod (order, ilat);
958 }
959 reference_element_H::size_type
960 reference_element_H::first_inod_by_variant (
961  size_type order,
962  size_type subgeo_variant)
963 {
964  if (subgeo_variant == reference_element::p) return 0;
965  if (subgeo_variant == reference_element::e) return 8;
966  if (subgeo_variant == reference_element::t) return 8 + 12*(order-1);
967  if (subgeo_variant == reference_element::q) return 8 + 12*(order-1);
968  if (subgeo_variant == reference_element::T) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
969  if (subgeo_variant == reference_element::P) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
970  if (subgeo_variant == reference_element::H) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
971  return (order+1)*(order+1)*(order+1);
972 }
973 Float
974 reference_element_H::side_measure (size_type loc_isid)
975 {
976  fatal_macro ("side_measure(3d): not yet");
977  return 0;
978 }
979 void
980 reference_element_H::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
981 {
982  fatal_macro ("side_normal(3d): not yet");
983 }
984 
985 
986 } // namespace rheolef
987