Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
cxx_reference_element.cc
Go to the documentation of this file.
1
21//
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
33using namespace rheolef;
34using namespace std;
35
36// --------------------------------------------------------------------
37// global table definitions
38// --------------------------------------------------------------------
39static const size_t max_size_t = size_t(-1);
40static const size_t max_variant = 7;
41static const char* table_name [max_variant] = { "p", "e", "t", "q", "T", "P", "H" };
42
43size_t table_dimension [max_variant];
44Float table_measure [max_variant];
45size_t table_n_vertex [max_variant];
46
47static const size_t max_face = 8; // for hexa : number of faces
48static const size_t max_face_vertex = 4; // for hexa : number of vertex on a face
49
50size_t table_n_edge [max_variant];
51size_t table_n_face [max_variant];
52size_t table_n_face_vertex_max [max_variant];
53size_t table_n_face_vertex [max_variant][max_face];
54size_t table_fac2edg_idx [max_variant][max_face][max_face_vertex];
55int table_fac2edg_ori [max_variant][max_face][max_face_vertex];
56
57// --------------------------------------------------------------------
58// generic fcts and specific includes
59// --------------------------------------------------------------------
60void 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}
66void 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}
69void 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}
73template <size_t NEdgePerFaceMax>
74void 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}
106void init_p (size_t p) {
107#include "point.icc"
109}
110void init_e (size_t e) {
111#include "edge.icc"
113}
114void init_t (size_t t) {
115#include "triangle.icc"
117}
118void init_q (size_t q) {
119#include "quadrangle.icc"
121}
122void init_T (size_t T) {
123#include "tetrahedron.icc"
125}
126void init_P (size_t P) {
127#include "prism.icc"
129}
130void init_H (size_t H) {
131#include "hexahedron.icc"
133}
134void 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// --------------------------------------------------------------------
322int 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}
see the Float page for the full documentation
see the edge page for the full documentation
see the point page for the full documentation
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)
size_t table_n_face[max_variant]
size_t table_dimension[max_variant]
int table_fac2edg_ori[max_variant][max_face][max_face_vertex]
size_t table_n_face_vertex[max_variant][max_face]
size_t table_n_vertex[max_variant]
void init_t(size_t t)
void init_q(size_t q)
void init_generic_0d(size_t E, size_t d, size_t nv, size_t ne, Float meas)
size_t table_n_face_vertex_max[max_variant]
void cxx_reference_element_header()
size_t table_fac2edg_idx[max_variant][max_face][max_face_vertex]
void cxx_reference_element_body()
void init_T(size_t T)
size_t table_n_edge[max_variant]
void init_e(size_t e)
void licence()
void init_P(size_t P)
void init_H(size_t H)
Float table_measure[max_variant]
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)
void init_p(size_t p)
void init_generic_1d(size_t E, size_t d, size_t nv, const point v[], size_t ne, Float meas)
edge - reference element
const size_t dimension
Definition edge.icc:64
const point vertex[n_vertex]
Definition edge.icc:67
const size_t n_vertex
Definition edge.icc:66
const Float measure
Definition edge.icc:65
int main()
Definition field2bb.cc:58
Expr1::float_type T
Definition field_expr.h:230
hexahedron - reference element
const size_t face[n_face][4]
const size_t n_face
const size_t n_edge
This file is part of Rheolef.
STL namespace.
point - reference element
prism - reference element
quadrangle - reference element
Definition cavity_dg.h:29
Definition sphere.icc:25
tetrahedron - reference element
triangle - reference element