Rheolef  7.1
an efficient C++ finite element environment
cxx_reference_element.cc
Go to the documentation of this file.
1 //
22 // reference element connectivity tables:
23 // - input : #include human-readable C++ code "hexa.h"..."triangle.h" etc
24 // - output: code efficiently usable internaly
25 // by the "reference_element" class
26 //
27 // author: Pierre.Saramito@imag.fr
28 //
29 // date: 24 may 2010
30 //
31 #include "rheolef/point.h"
32 
33 using namespace rheolef;
34 using namespace std;
35 
36 // --------------------------------------------------------------------
37 // global table definitions
38 // --------------------------------------------------------------------
39 static const size_t max_size_t = size_t(-1);
40 static const size_t max_variant = 7;
41 static const char* table_name [max_variant] = { "p", "e", "t", "q", "T", "P", "H" };
42 
43 size_t table_dimension [max_variant];
44 Float table_measure [max_variant];
45 size_t table_n_vertex [max_variant];
46 
47 static const size_t max_face = 8; // for hexa : number of faces
48 static const size_t max_face_vertex = 4; // for hexa : number of vertex on a face
49 
50 size_t table_n_edge [max_variant];
51 size_t table_n_face [max_variant];
52 size_t table_n_face_vertex_max [max_variant];
53 size_t table_n_face_vertex [max_variant][max_face];
54 size_t table_fac2edg_idx [max_variant][max_face][max_face_vertex];
55 int table_fac2edg_ori [max_variant][max_face][max_face_vertex];
56 
57 // --------------------------------------------------------------------
58 // generic fcts and specific includes
59 // --------------------------------------------------------------------
60 void init_generic_0d (size_t E, size_t d, size_t nv, size_t ne, Float meas) {
61  table_dimension [E] = d;
62  table_measure [E] = meas;
63  table_n_vertex [E] = nv;
64  table_n_edge [E] = ne;
65 }
66 void init_generic_1d (size_t E, size_t d, size_t nv, const point v[], size_t ne, Float meas) {
67  init_generic_0d (E, d, nv, ne, meas);
68 }
69 void init_generic_2d (size_t E, size_t d, size_t nv, const point v[],
70  size_t ne, const size_t e[][2], Float meas) {
71  init_generic_1d (E, d, nv, v, ne, meas);
72 }
73 template <size_t NEdgePerFaceMax>
74 void init_generic_3d (size_t E, size_t d, size_t nv, const point v[],
75  size_t nfac, const size_t f[][NEdgePerFaceMax],
76  size_t nedg, const size_t e[][2],
77  Float meas)
78 {
79  init_generic_2d (E, d, nv, v, nedg, e, meas);
80  table_n_face [E] = nfac;
82  for (size_t ifac = 0; ifac < nfac; ifac++) {
83  size_t nv_on_fac = 0;
84  while (nv_on_fac < NEdgePerFaceMax && f[ifac][nv_on_fac] != size_t(-1)) {
85  nv_on_fac++;
86  }
87  table_n_face_vertex [E][ifac] = nv_on_fac;
88  table_n_face_vertex_max [E] = std::max (nv_on_fac, table_n_face_vertex_max [E]);
89  for (size_t iv = 0; iv < nv_on_fac; iv++) {
90  size_t iv2 = (iv+1) % nv_on_fac;
91  // search (iv,iv2) edge in the edge table e[]
92  for (size_t iedg = 0; iedg < nedg; iedg++) {
93  if (f[ifac][iv] == e[iedg][0] && f[ifac][iv2] == e[iedg][1]) {
94  table_fac2edg_idx [E][ifac][iv] = iedg;
95  table_fac2edg_ori [E][ifac][iv] = 1;
96  break;
97  } else if (f[ifac][iv] == e[iedg][1] && f[ifac][iv2] == e[iedg][0]) {
98  table_fac2edg_idx [E][ifac][iv] = iedg;
99  table_fac2edg_ori [E][ifac][iv] = -1;
100  break;
101  }
102  }
103  }
104  }
105 }
106 void init_p (size_t p) {
107 #include "point.icc"
109 }
110 void init_e (size_t e) {
111 #include "edge.icc"
113 }
114 void init_t (size_t t) {
115 #include "triangle.icc"
117 }
118 void init_q (size_t q) {
119 #include "quadrangle.icc"
121 }
122 void init_T (size_t T) {
123 #include "tetrahedron.icc"
125 }
126 void init_P (size_t P) {
127 #include "prism.icc"
129 }
130 void init_H (size_t H) {
131 #include "hexahedron.icc"
133 }
134 void licence () {
135 // --------------------------------------------------------------------
136 // c++ code generation
137 // --------------------------------------------------------------------
138  cout << "// file automaticaly generated by \"cxx_reference_element.cc\"" << endl
139  << "//" << endl
140  << "///" << endl
141  << "/// This file is part of Rheolef." << endl
142  << "///" << endl
143  << "/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>" << endl
144  << "///" << endl
145  << "/// Rheolef is free software; you can redistribute it and/or modify" << endl
146  << "/// it under the terms of the GNU General Public License as published by" << endl
147  << "/// the Free Software Foundation; either version 2 of the License, or" << endl
148  << "/// (at your option) any later version." << endl
149  << "///" << endl
150  << "/// Rheolef is distributed in the hope that it will be useful," << endl
151  << "/// but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl
152  << "/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl
153  << "/// GNU General Public License for more details." << endl
154  << "///" << endl
155  << "/// You should have received a copy of the GNU General Public License" << endl
156  << "/// along with Rheolef; if not, write to the Free Software" << endl
157  << "/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" << endl
158  << "/// " << endl
159  << "/// =========================================================================" << endl
160  << "// file automaticaly generated by \"cxx_reference_element.cc\"" << endl
161  ;
162 }
164  licence ();
165  // print enum symbol = e, t, q, ..
166  cout << "static const variant_type" << endl;
167  for (size_t i = 0; i < max_variant; i++) {
168  cout << "\t" << table_name[i] << " = " << i << "," << endl;
169  }
170  cout << "\tmax_variant = " << max_variant << ";" << endl;
171 }
173  licence ();
174 
175  // print char symbol = 'e', 't', ..
176  cout << "const char" << endl
177  << "reference_element::_name [reference_element::max_variant] = {" << endl
178  ;
179  for (size_t E = 0; E < max_variant; E++) {
180  cout << "\t'" << table_name[E] << "'";
181  if (E+1 != max_variant) cout << ",";
182  cout << endl;
183  }
184  cout << "};" << endl;
185 
186  // element dimension : 1,2 3
187  cout << "const reference_element::size_type" << endl
188  << "reference_element::_dimension [reference_element::max_variant] = {" << endl
189  ;
190  for (size_t E = 0; E < max_variant; E++) {
191  cout << "\t" << table_dimension[E];
192  if (E+1 != max_variant) cout << ",";
193  cout << endl;
194  }
195  cout << "};" << endl;
196 
197  // first variant by dimension
198  cout << "const reference_element::size_type" << endl
199  << "reference_element::_first_variant_by_dimension [5] = {" << endl
200  ;
201  for (size_t E = 0, prev_dim = size_t(-1); E < max_variant; E++) {
202  if (table_dimension[E] == prev_dim) continue;
203  prev_dim = table_dimension[E];
204  cout << "\treference_element::" << table_name[E] << "," << endl;
205  }
206  cout << "\treference_element::max_variant" << endl
207  << "};" << endl;
208 
209  // n_vertex
210  cout << "const reference_element::size_type" << endl
211  << "reference_element::_n_vertex [reference_element::max_variant] = {" << endl
212  ;
213  for (size_t E = 0; E < max_variant; E++) {
214  cout << "\t" << table_n_vertex[E];
215  if (E+1 != max_variant) cout << ",";
216  cout << endl;
217  }
218  cout << "};" << endl;
219 
220  // n_subgeo_by_variant
221  cout << "const reference_element::size_type" << endl
222  << "reference_element::_n_subgeo_by_variant [reference_element::max_variant] [reference_element::max_variant] = {" << endl
223  << "//";
224  for (size_t E = 0; E < max_variant; E++) {
225  cout << " " << table_name[E];
226  }
227  cout << endl;
228  for (size_t E = 0; E < max_variant; E++) {
229  cout << " {";
230  for (size_t F = 0; F < max_variant; F++) {
231  size_t nf;
232  if (F > E) {
233  nf = 0;
234  } else if (F == E) {
235  nf = 1;
236  } else if (table_dimension[F] == table_dimension[E]) {
237  nf = 0;
238  } else { // table_dimension[F] < table_dimension[E]
239  switch (table_dimension[F]) {
240  case 0: nf = table_n_vertex [E]; break;
241  case 1: nf = table_n_edge [E]; break;
242  default:
243  case 2: {
244  nf = 0;
245  for (size_t ifac = 0; ifac < table_n_face [E]; ++ifac) {
246  if (table_n_face_vertex [E][ifac] == table_n_vertex[F]) nf++;
247  }
248  break;
249  }
250  }
251  }
252  cout << " " << nf;
253  if (F+1 != max_variant) cout << ",";
254  }
255 #ifdef TODO
256  { 1, 0, 0, 0, 0, 0, 0}, // p
257  { 2, 1, 0, 0, 0, 0, 0}, // e
258  { 3, 3, 1, 0, 0, 0, 0}, // t
259  { 4, 4, 0, 1, 0, 0, 0}, // q
260  { 4, 6, 4, 0, 1, 0, 0}, // T
261  { 6, 9, 2, 3, 0, 1, 0}, // P
262  { 8,12, 0, 6, 0, 0, 1} // H
263 #endif // TODO
264  cout << "}" << ((E+1 != max_variant) ? "," : " ") << " // " << table_name[E] << endl;
265  }
266  cout << "};" << endl;
267  // measure
268  cout << "static const Float" << endl
269  << "hat_K_measure [reference_element::max_variant] = {" << endl
270  << setprecision(std::numeric_limits<Float>::digits10)
271  ;
272  for (size_t E = 0; E < max_variant; E++) {
273  cout << "\t" << table_measure[E];
274  if (E+1 != max_variant) cout << ",";
275  cout << endl;
276  }
277  cout << "};" << endl;
278 
279  // 3d faces: edge indexes
280  for (size_t E = 0; E < max_variant; E++) {
281  if (table_dimension[E] != 3) continue;
282  size_t nfac = table_n_face [E];
283  cout << "const reference_element::size_type" << endl
284  << "geo_element_" << table_name[E] << "_fac2edg_idx ["
285  <<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
286  ;
287  for (size_t ifac = 0; ifac < nfac; ifac++) {
288  size_t nedg = table_n_face_vertex [E][ifac];
289  cout << "\t{";
290  for (size_t iedg = 0; iedg < nedg; iedg++) {
291  cout << table_fac2edg_idx [E][ifac][iedg];
292  if (iedg+1 == nedg) cout << "}"; else cout << ",";
293  }
294  if (ifac+1 == nfac) cout << " };"; else cout << ",";
295  cout << endl;
296  }
297  cout << endl;
298  }
299 
300  // 3d faces: edge orientation
301  for (size_t E = 0; E < max_variant; E++) {
302  if (table_dimension[E] != 3) continue;
303  size_t nfac = table_n_face [E];
304  cout << "const int" << endl
305  << "geo_element_" << table_name[E] << "_fac2edg_orient ["
306  <<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
307  ;
308  for (size_t ifac = 0; ifac < nfac; ifac++) {
309  size_t nedg = table_n_face_vertex [E][ifac];
310  cout << "\t{";
311  for (size_t iedg = 0; iedg < nedg; iedg++) {
312  cout << table_fac2edg_ori [E][ifac][iedg];
313  if (iedg+1 == nedg) cout << "}"; else cout << ",";
314  }
315  if (ifac+1 == nfac) cout << " };"; else cout << ",";
316  cout << endl;
317  }
318  cout << endl;
319  }
320 }
321 // --------------------------------------------------------------------
322 int main(int argc, char**argv) {
323 // --------------------------------------------------------------------
324  size_t E = 0;
325  init_p (E++);
326  init_e (E++);
327  init_t (E++);
328  init_q (E++);
329  init_T (E++);
330  init_P (E++);
331  init_H (E++);
332  if (argc > 1 && argv[1] == string("-h")) {
334  } else {
336  }
337 }
init_q
void init_q(size_t q)
Definition: cxx_reference_element.cc:118
n_face
const size_t n_face
Definition: hexahedron.icc:109
init_t
void init_t(size_t t)
Definition: cxx_reference_element.cc:114
n_edge
const size_t n_edge
Definition: hexahedron.icc:117
init_generic_0d
void init_generic_0d(size_t E, size_t d, size_t nv, size_t ne, Float meas)
Definition: cxx_reference_element.cc:60
triangle.icc
triangle - reference element
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
table_fac2edg_idx
size_t table_fac2edg_idx[max_variant][max_face][max_face_vertex]
Definition: cxx_reference_element.cc:54
init_generic_1d
void init_generic_1d(size_t E, size_t d, size_t nv, const point v[], size_t ne, Float meas)
Definition: cxx_reference_element.cc:66
table_n_vertex
size_t table_n_vertex[max_variant]
Definition: cxx_reference_element.cc:45
cxx_reference_element_header
void cxx_reference_element_header()
Definition: cxx_reference_element.cc:163
p
Definition: sphere.icc:25
table_n_face_vertex_max
size_t table_n_face_vertex_max[max_variant]
Definition: cxx_reference_element.cc:52
quadrangle.icc
quadrangle - reference element
prism.icc
prism - reference element
init_H
void init_H(size_t H)
Definition: cxx_reference_element.cc:130
vertex
const point vertex[n_vertex]
Definition: edge.icc:67
table_n_edge
size_t table_n_edge[max_variant]
Definition: cxx_reference_element.cc:50
table_dimension
size_t table_dimension[max_variant]
Definition: cxx_reference_element.cc:43
rheolef::measure
Float measure(reference_element hat_K)
Definition: reference_element.cc:65
init_T
void init_T(size_t T)
Definition: cxx_reference_element.cc:122
cxx_reference_element_body
void cxx_reference_element_body()
Definition: cxx_reference_element.cc:172
tetrahedron.icc
tetrahedron - reference element
table_fac2edg_ori
int table_fac2edg_ori[max_variant][max_face][max_face_vertex]
Definition: cxx_reference_element.cc:55
point.icc
point - reference element
table_measure
Float table_measure[max_variant]
Definition: cxx_reference_element.cc:44
dimension
const size_t dimension
Definition: edge.icc:64
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
edge.icc
edge - reference element
face
const size_t face[n_face][4]
Definition: hexahedron.icc:110
table_n_face
size_t table_n_face[max_variant]
Definition: cxx_reference_element.cc:51
Float
see the Float page for the full documentation
point
see the point page for the full documentation
init_generic_3d
void init_generic_3d(size_t E, size_t d, size_t nv, const point v[], size_t nfac, const size_t f[][NEdgePerFaceMax], size_t nedg, const size_t e[][2], Float meas)
Definition: cxx_reference_element.cc:74
init_P
void init_P(size_t P)
Definition: cxx_reference_element.cc:126
hexahedron.icc
hexahedron - reference element
init_generic_2d
void init_generic_2d(size_t E, size_t d, size_t nv, const point v[], size_t ne, const size_t e[][2], Float meas)
Definition: cxx_reference_element.cc:69
licence
void licence()
Definition: cxx_reference_element.cc:134
main
int main(int argc, char **argv)
Definition: cxx_reference_element.cc:322
table_n_face_vertex
size_t table_n_face_vertex[max_variant][max_face]
Definition: cxx_reference_element.cc:53
init_p
void init_p(size_t p)
Definition: cxx_reference_element.cc:106
f
Definition: cavity_dg.h:29
n_vertex
const size_t n_vertex
Definition: edge.icc:66
rheolef::std
Definition: vec_expr_v2.h:402
T
Expr1::float_type T
Definition: field_expr.h:261
init_e
void init_e(size_t e)
Definition: cxx_reference_element.cc:110
edge
see the edge page for the full documentation