← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/geo_element/cxx_reference_element.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/cxx_reference_element.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
/// 
 
20
/// =========================================================================
 
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
 
 
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_face            [max_variant];
 
51
size_t table_n_face_vertex_max [max_variant];
 
52
size_t table_n_face_vertex     [max_variant][max_face];
 
53
size_t table_fac2edg_idx       [max_variant][max_face][max_face_vertex];
 
54
int    table_fac2edg_ori       [max_variant][max_face][max_face_vertex];
 
55
 
 
56
// --------------------------------------------------------------------
 
57
// generic fcts and specific includes
 
58
// --------------------------------------------------------------------
 
59
void init_generic_0d (size_t E, size_t d, size_t nv, Float meas) {
 
60
  table_dimension [E] = d;
 
61
  table_measure   [E] = meas;
 
62
  table_n_vertex  [E] = nv;
 
63
}
 
64
template <int Dimension>
 
65
void init_generic_1d (size_t E, size_t d, size_t nv, const Float v[][Dimension], Float meas) {
 
66
  init_generic_0d (E, d, nv, meas);
 
67
}
 
68
template <int Dimension>
 
69
void init_generic_2d (size_t E, size_t d, size_t nv, const Float v[][Dimension],
 
70
        size_t ne, const size_t e[][2], Float meas) {
 
71
  init_generic_1d (E, d, nv, v, meas);
 
72
}
 
73
template <size_t NEdgePerFaceMax>
 
74
void init_generic_3d (size_t E, size_t d, size_t nv, const Float v[][3],
 
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;
 
81
  table_n_face_vertex_max [E] = 0;
 
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"
 
108
  init_generic_0d (p, dimension, 1, measure);
 
109
}
 
110
void init_e (size_t e) {
 
111
#include "edge.icc"
 
112
  init_generic_1d (e, dimension, n_vertex, vertex, measure);
 
113
}
 
114
void init_t (size_t t) {
 
115
#include "triangle.icc"
 
116
  init_generic_2d (t, dimension, n_vertex, vertex, n_edge, edge, measure);
 
117
}
 
118
void init_q (size_t q) {
 
119
#include "quadrangle.icc"
 
120
  init_generic_2d (q, dimension, n_vertex, vertex, n_edge, edge, measure);
 
121
}
 
122
void init_T (size_t T) {
 
123
#include "tetrahedron.icc"
 
124
  init_generic_3d (T, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
 
125
}
 
126
void init_P (size_t P) {
 
127
#include "prism.icc"
 
128
  init_generic_3d (P, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
 
129
}
 
130
void init_H (size_t H) {
 
131
#include "hexahedron.icc"
 
132
  init_generic_3d (H, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
 
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
       ;
 
161
}
 
162
void cxx_reference_element_header() {
 
163
  licence ();
 
164
  // print enum symbol = e, t, q, ..
 
165
  cout << "static const variant_type" << endl;
 
166
  for (size_t i = 0; i < max_variant; i++) {
 
167
    cout << "\t" << table_name[i] << " = " << i << "," << endl;
 
168
  }
 
169
  cout << "\tmax_variant = " << max_variant << ";" << endl;
 
170
}
 
171
void cxx_reference_element_body () {
 
172
  licence ();
 
173
 
 
174
  // print char symbol = 'e', 't', ..
 
175
  cout << "const char" << endl
 
176
       << "reference_element::_name [reference_element::max_variant] = {" << endl
 
177
    ;
 
178
  for (size_t E = 0; E < max_variant; E++) {
 
179
    cout << "\t'" << table_name[E] << "'";
 
180
    if (E+1 != max_variant) cout << ",";
 
181
    cout << endl;
 
182
  }
 
183
  cout << "};" << endl;
 
184
 
 
185
  // element dimension : 1,2 3
 
186
  cout << "const reference_element::size_type" << endl
 
187
       << "reference_element::_dimension [reference_element::max_variant] = {" << endl
 
188
    ;
 
189
  for (size_t E = 0; E < max_variant; E++) {
 
190
    cout << "\t" << table_dimension[E];
 
191
    if (E+1 != max_variant) cout << ",";
 
192
    cout << endl;
 
193
  }
 
194
  cout << "};" << endl;
 
195
 
 
196
  // first variant by dimension
 
197
  cout << "const reference_element::size_type" << endl
 
198
       << "reference_element::_first_variant_by_dimension [5] = {" << endl
 
199
    ;
 
200
  for (size_t E = 0, prev_dim = size_t(-1); E < max_variant; E++) {
 
201
    if (table_dimension[E] == prev_dim) continue;
 
202
    prev_dim = table_dimension[E];
 
203
    cout << "\treference_element::" << table_name[E] << "," << endl;
 
204
  }
 
205
  cout << "\treference_element::max_variant" << endl
 
206
       << "};" << endl;
 
207
 
 
208
  // n_vertex
 
209
  cout << "const reference_element::size_type" << endl
 
210
       << "reference_element::_n_vertex [reference_element::max_variant] = {" << endl
 
211
    ;
 
212
  for (size_t E = 0; E < max_variant; E++) {
 
213
    cout << "\t" << table_n_vertex[E];
 
214
    if (E+1 != max_variant) cout << ",";
 
215
    cout << endl;
 
216
  }
 
217
  cout << "};" << endl;
 
218
 
 
219
  // measure
 
220
  cout << "static const Float" << endl
 
221
       << "hat_K_measure [reference_element::max_variant] = {" << endl
 
222
       << setprecision(std::numeric_limits<Float>::digits10)
 
223
    ;
 
224
  for (size_t E = 0; E < max_variant; E++) {
 
225
    cout << "\t" << table_measure[E];
 
226
    if (E+1 != max_variant) cout << ",";
 
227
    cout << endl;
 
228
  }
 
229
  cout << "};" << endl;
 
230
 
 
231
  // 3d faces: edge indexes
 
232
  for (size_t E = 0; E < max_variant; E++) {
 
233
    if (table_dimension[E] != 3) continue;
 
234
    size_t nfac = table_n_face [E];
 
235
    cout << "const reference_element::size_type" << endl
 
236
         << "geo_element_" << table_name[E] << "_fac2edg_idx ["
 
237
             <<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
 
238
    ;
 
239
    for (size_t ifac = 0; ifac < nfac; ifac++) {
 
240
      size_t nedg = table_n_face_vertex [E][ifac];
 
241
      cout << "\t{";
 
242
      for (size_t iedg = 0; iedg < nedg; iedg++) {
 
243
        cout << table_fac2edg_idx [E][ifac][iedg];
 
244
        if (iedg+1 == nedg) cout << "}"; else cout << ",";
 
245
      }
 
246
      if (ifac+1 == nfac) cout << " };"; else cout << ",";
 
247
      cout << endl;
 
248
    }
 
249
    cout << endl;
 
250
  }
 
251
 
 
252
  // 3d faces: edge orientation
 
253
  for (size_t E = 0; E < max_variant; E++) {
 
254
    if (table_dimension[E] != 3) continue;
 
255
    size_t nfac = table_n_face [E];
 
256
    cout << "const int" << endl
 
257
         << "geo_element_" << table_name[E] << "_fac2edg_orient ["
 
258
             <<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
 
259
    ;
 
260
    for (size_t ifac = 0; ifac < nfac; ifac++) {
 
261
      size_t nedg = table_n_face_vertex [E][ifac];
 
262
      cout << "\t{";
 
263
      for (size_t iedg = 0; iedg < nedg; iedg++) {
 
264
        cout << table_fac2edg_ori [E][ifac][iedg];
 
265
        if (iedg+1 == nedg) cout << "}"; else cout << ",";
 
266
      }
 
267
      if (ifac+1 == nfac) cout << " };"; else cout << ",";
 
268
      cout << endl;
 
269
    }
 
270
    cout << endl;
 
271
  }
 
272
}
 
273
// --------------------------------------------------------------------
 
274
int main(int argc, char**argv) {
 
275
// --------------------------------------------------------------------
 
276
  size_t E = 0;
 
277
  init_p (E++);
 
278
  init_e (E++);
 
279
  init_t (E++);
 
280
  init_q (E++);
 
281
  init_T (E++);
 
282
  init_P (E++);
 
283
  init_H (E++);
 
284
  if (argc > 1 && argv[1] == string("-h")) {
 
285
    cxx_reference_element_header();
 
286
  } else {
 
287
    cxx_reference_element_body();
 
288
  }
 
289
}