Loading...
Searching...
No Matches
cxx_reference_element.cc
Go to the documentation of this file.
//
// reference element connectivity tables:
// - input : #include human-readable C++ code "hexa.h"..."triangle.h" etc
// - output: code efficiently usable internaly
// by the "reference_element" class
//
// author: Pierre.Saramito@imag.fr
//
// date: 24 may 2010
//
#include "rheolef/point.h"
using namespace rheolef;
using namespace std;
// --------------------------------------------------------------------
// global table definitions
// --------------------------------------------------------------------
static const size_t max_size_t = size_t(-1);
static const size_t max_variant = 7;
static const char* table_name [max_variant] = { "p", "e", "t", "q", "T", "P", "H" };
size_t table_dimension [max_variant];
Float table_measure [max_variant];
size_t table_n_vertex [max_variant];
static const size_t max_face = 8; // for hexa : number of faces
static const size_t max_face_vertex = 4; // for hexa : number of vertex on a face
size_t table_n_edge [max_variant];
size_t table_n_face [max_variant];
size_t table_n_face_vertex_max [max_variant];
size_t table_n_face_vertex [max_variant][max_face];
size_t table_fac2edg_idx [max_variant][max_face][max_face_vertex];
int table_fac2edg_ori [max_variant][max_face][max_face_vertex];
// --------------------------------------------------------------------
// generic fcts and specific includes
// --------------------------------------------------------------------
void init_generic_0d (size_t E, size_t d, size_t nv, size_t ne, Float meas) {
table_dimension [E] = d;
table_measure [E] = meas;
table_n_vertex [E] = nv;
table_n_edge [E] = ne;
}
void init_generic_1d (size_t E, size_t d, size_t nv, const point v[], size_t ne, Float meas) {
init_generic_0d (E, d, nv, ne, meas);
}
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) {
init_generic_1d (E, d, nv, v, ne, meas);
}
template <size_t NEdgePerFaceMax>
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)
{
init_generic_2d (E, d, nv, v, nedg, e, meas);
table_n_face [E] = nfac;
table_n_face_vertex_max [E] = 0;
for (size_t ifac = 0; ifac < nfac; ifac++) {
size_t nv_on_fac = 0;
while (nv_on_fac < NEdgePerFaceMax && f[ifac][nv_on_fac] != size_t(-1)) {
nv_on_fac++;
}
table_n_face_vertex [E][ifac] = nv_on_fac;
table_n_face_vertex_max [E] = std::max (nv_on_fac, table_n_face_vertex_max [E]);
for (size_t iv = 0; iv < nv_on_fac; iv++) {
size_t iv2 = (iv+1) % nv_on_fac;
// search (iv,iv2) edge in the edge table e[]
for (size_t iedg = 0; iedg < nedg; iedg++) {
if (f[ifac][iv] == e[iedg][0] && f[ifac][iv2] == e[iedg][1]) {
table_fac2edg_idx [E][ifac][iv] = iedg;
table_fac2edg_ori [E][ifac][iv] = 1;
break;
} else if (f[ifac][iv] == e[iedg][1] && f[ifac][iv2] == e[iedg][0]) {
table_fac2edg_idx [E][ifac][iv] = iedg;
table_fac2edg_ori [E][ifac][iv] = -1;
break;
}
}
}
}
}
void init_p (size_t p) {
#include "point.icc"
init_generic_0d (p, dimension, 1, 0, measure);
}
void init_e (size_t e) {
#include "edge.icc"
init_generic_1d (e, dimension, n_vertex, vertex, 1, measure);
}
void init_t (size_t t) {
#include "triangle.icc"
init_generic_2d (t, dimension, n_vertex, vertex, n_edge, edge, measure);
}
void init_q (size_t q) {
#include "quadrangle.icc"
init_generic_2d (q, dimension, n_vertex, vertex, n_edge, edge, measure);
}
void init_T (size_t T) {
#include "tetrahedron.icc"
init_generic_3d (T, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
}
void init_P (size_t P) {
#include "prism.icc"
init_generic_3d (P, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
}
void init_H (size_t H) {
#include "hexahedron.icc"
init_generic_3d (H, dimension, n_vertex, vertex, n_face, face, n_edge, edge, measure);
}
void licence () {
// --------------------------------------------------------------------
// c++ code generation
// --------------------------------------------------------------------
cout << "// file automaticaly generated by \"cxx_reference_element.cc\"" << endl
<< "//" << endl
<< "///" << endl
<< "/// This file is part of Rheolef." << endl
<< "///" << endl
<< "/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>" << endl
<< "///" << endl
<< "/// Rheolef is free software; you can redistribute it and/or modify" << endl
<< "/// it under the terms of the GNU General Public License as published by" << endl
<< "/// the Free Software Foundation; either version 2 of the License, or" << endl
<< "/// (at your option) any later version." << endl
<< "///" << endl
<< "/// Rheolef is distributed in the hope that it will be useful," << endl
<< "/// but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl
<< "/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl
<< "/// GNU General Public License for more details." << endl
<< "///" << endl
<< "/// You should have received a copy of the GNU General Public License" << endl
<< "/// along with Rheolef; if not, write to the Free Software" << endl
<< "/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" << endl
<< "/// " << endl
<< "/// =========================================================================" << endl
<< "// file automaticaly generated by \"cxx_reference_element.cc\"" << endl
;
}
void cxx_reference_element_header() {
licence ();
// print enum symbol = e, t, q, ..
cout << "static const variant_type" << endl;
for (size_t i = 0; i < max_variant; i++) {
cout << "\t" << table_name[i] << " = " << i << "," << endl;
}
cout << "\tmax_variant = " << max_variant << ";" << endl;
}
void cxx_reference_element_body () {
licence ();
// print char symbol = 'e', 't', ..
cout << "const char" << endl
<< "reference_element::_name [reference_element::max_variant] = {" << endl
;
for (size_t E = 0; E < max_variant; E++) {
cout << "\t'" << table_name[E] << "'";
if (E+1 != max_variant) cout << ",";
cout << endl;
}
cout << "};" << endl;
// element dimension : 1,2 3
cout << "const reference_element::size_type" << endl
<< "reference_element::_dimension [reference_element::max_variant] = {" << endl
;
for (size_t E = 0; E < max_variant; E++) {
cout << "\t" << table_dimension[E];
if (E+1 != max_variant) cout << ",";
cout << endl;
}
cout << "};" << endl;
// first variant by dimension
cout << "const reference_element::size_type" << endl
<< "reference_element::_first_variant_by_dimension [5] = {" << endl
;
for (size_t E = 0, prev_dim = size_t(-1); E < max_variant; E++) {
if (table_dimension[E] == prev_dim) continue;
prev_dim = table_dimension[E];
cout << "\treference_element::" << table_name[E] << "," << endl;
}
cout << "\treference_element::max_variant" << endl
<< "};" << endl;
// n_vertex
cout << "const reference_element::size_type" << endl
<< "reference_element::_n_vertex [reference_element::max_variant] = {" << endl
;
for (size_t E = 0; E < max_variant; E++) {
cout << "\t" << table_n_vertex[E];
if (E+1 != max_variant) cout << ",";
cout << endl;
}
cout << "};" << endl;
// n_subgeo_by_variant
cout << "const reference_element::size_type" << endl
<< "reference_element::_n_subgeo_by_variant [reference_element::max_variant] [reference_element::max_variant] = {" << endl
<< "//";
for (size_t E = 0; E < max_variant; E++) {
cout << " " << table_name[E];
}
cout << endl;
for (size_t E = 0; E < max_variant; E++) {
cout << " {";
for (size_t F = 0; F < max_variant; F++) {
size_t nf;
if (F > E) {
nf = 0;
} else if (F == E) {
nf = 1;
} else if (table_dimension[F] == table_dimension[E]) {
nf = 0;
} else { // table_dimension[F] < table_dimension[E]
switch (table_dimension[F]) {
case 0: nf = table_n_vertex [E]; break;
case 1: nf = table_n_edge [E]; break;
default:
case 2: {
nf = 0;
for (size_t ifac = 0; ifac < table_n_face [E]; ++ifac) {
if (table_n_face_vertex [E][ifac] == table_n_vertex[F]) nf++;
}
break;
}
}
}
cout << " " << nf;
if (F+1 != max_variant) cout << ",";
}
#ifdef TODO
{ 1, 0, 0, 0, 0, 0, 0}, // p
{ 2, 1, 0, 0, 0, 0, 0}, // e
{ 3, 3, 1, 0, 0, 0, 0}, // t
{ 4, 4, 0, 1, 0, 0, 0}, // q
{ 4, 6, 4, 0, 1, 0, 0}, // T
{ 6, 9, 2, 3, 0, 1, 0}, // P
{ 8,12, 0, 6, 0, 0, 1} // H
#endif // TODO
cout << "}" << ((E+1 != max_variant) ? "," : " ") << " // " << table_name[E] << endl;
}
cout << "};" << endl;
// measure
cout << "static const Float" << endl
<< "hat_K_measure [reference_element::max_variant] = {" << endl
<< setprecision(std::numeric_limits<Float>::digits10)
;
for (size_t E = 0; E < max_variant; E++) {
cout << "\t" << table_measure[E];
if (E+1 != max_variant) cout << ",";
cout << endl;
}
cout << "};" << endl;
// 3d faces: edge indexes
for (size_t E = 0; E < max_variant; E++) {
if (table_dimension[E] != 3) continue;
size_t nfac = table_n_face [E];
cout << "const reference_element::size_type" << endl
<< "geo_element_" << table_name[E] << "_fac2edg_idx ["
<<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
;
for (size_t ifac = 0; ifac < nfac; ifac++) {
size_t nedg = table_n_face_vertex [E][ifac];
cout << "\t{";
for (size_t iedg = 0; iedg < nedg; iedg++) {
cout << table_fac2edg_idx [E][ifac][iedg];
if (iedg+1 == nedg) cout << "}"; else cout << ",";
}
if (ifac+1 == nfac) cout << " };"; else cout << ",";
cout << endl;
}
cout << endl;
}
// 3d faces: edge orientation
for (size_t E = 0; E < max_variant; E++) {
if (table_dimension[E] != 3) continue;
size_t nfac = table_n_face [E];
cout << "const int" << endl
<< "geo_element_" << table_name[E] << "_fac2edg_orient ["
<<nfac<<"]["<< table_n_face_vertex_max [E] << "] = {" << endl
;
for (size_t ifac = 0; ifac < nfac; ifac++) {
size_t nedg = table_n_face_vertex [E][ifac];
cout << "\t{";
for (size_t iedg = 0; iedg < nedg; iedg++) {
cout << table_fac2edg_ori [E][ifac][iedg];
if (iedg+1 == nedg) cout << "}"; else cout << ",";
}
if (ifac+1 == nfac) cout << " };"; else cout << ",";
cout << endl;
}
cout << endl;
}
}
// --------------------------------------------------------------------
int main(int argc, char**argv) {
// --------------------------------------------------------------------
size_t E = 0;
init_p (E++);
init_e (E++);
init_t (E++);
init_q (E++);
init_T (E++);
init_P (E++);
init_H (E++);
if (argc > 1 && argv[1] == string("-h")) {
cxx_reference_element_header();
} else {
cxx_reference_element_body();
}
}
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
table_n_face
size_t table_n_face[max_variant]
Definition cxx_reference_element.cc:51
table_dimension
size_t table_dimension[max_variant]
Definition cxx_reference_element.cc:43
table_fac2edg_ori
int table_fac2edg_ori[max_variant][max_face][max_face_vertex]
Definition cxx_reference_element.cc:55
table_n_face_vertex
size_t table_n_face_vertex[max_variant][max_face]
Definition cxx_reference_element.cc:53
table_n_vertex
size_t table_n_vertex[max_variant]
Definition cxx_reference_element.cc:45
init_t
void init_t(size_t t)
Definition cxx_reference_element.cc:114
init_q
void init_q(size_t q)
Definition cxx_reference_element.cc:118
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
table_n_face_vertex_max
size_t table_n_face_vertex_max[max_variant]
Definition cxx_reference_element.cc:52
cxx_reference_element_header
void cxx_reference_element_header()
Definition cxx_reference_element.cc:163
table_fac2edg_idx
size_t table_fac2edg_idx[max_variant][max_face][max_face_vertex]
Definition cxx_reference_element.cc:54
cxx_reference_element_body
void cxx_reference_element_body()
Definition cxx_reference_element.cc:172
init_T
void init_T(size_t T)
Definition cxx_reference_element.cc:122
table_n_edge
size_t table_n_edge[max_variant]
Definition cxx_reference_element.cc:50
init_e
void init_e(size_t e)
Definition cxx_reference_element.cc:110
licence
void licence()
Definition cxx_reference_element.cc:134
init_P
void init_P(size_t P)
Definition cxx_reference_element.cc:126
init_H
void init_H(size_t H)
Definition cxx_reference_element.cc:130
table_measure
Float table_measure[max_variant]
Definition cxx_reference_element.cc:44
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
init_p
void init_p(size_t p)
Definition cxx_reference_element.cc:106
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
edge.icc
edge - reference element
dimension
const size_t dimension
Definition edge.icc:64
vertex
const point vertex[n_vertex]
Definition edge.icc:67
n_vertex
const size_t n_vertex
Definition edge.icc:66
measure
const Float measure
Definition edge.icc:65
main
int main()
Definition field2bb.cc:58
T
Expr1::float_type T
Definition field_expr.h:230
hexahedron.icc
hexahedron - reference element
face
const size_t face[n_face][4]
Definition hexahedron.icc:110
n_face
const size_t n_face
Definition hexahedron.icc:109
n_edge
const size_t n_edge
Definition hexahedron.icc:117
rheolef::edge
Definition reference_element.cc:29
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
point.icc
point - reference element
prism.icc
prism - reference element
quadrangle.icc
quadrangle - reference element
f
Definition cavity_dg.h:29
p
Definition sphere.icc:25
tetrahedron.icc
tetrahedron - reference element
triangle.icc
triangle - reference element