msh2geo_node_renum.icc
Go to the documentation of this file.
//
// for high order gmsh files:
// reorder the node list from gmsh 2 rheolef conventions:
// - d=2,3 : for face & volume internal nodes,
// change from recursive gmsh numbering style
// to rheolef lattice style numbering
// - d=3 : change edge & face local index and orientation
//
// author: Pierre.Saramito@imag.fr
//
// date: 18 september 2011
//
namespace rheolef {
template<class Iterator>
static void check_permutation (Iterator p, size_t n, std::string name)
{
bool ok = true;
for (size_t i = 0; i < n; i++) {
if (p[i] >= n) {
ok = false;
}
}
vector<size_t> ip (n, numeric_limits<size_t>::max());
for (size_t i = 0; i < n; i++) {
if (ip[p[i]] < n) {
ok = false;
}
ip[p[i]] = i;
}
for (size_t i = 0; i < n; i++) {
if (ip[p[i]] != i) {
ok = false;
}
}
check_macro (ok, "permutation("<<n<<"): "<<name<<" has problem");
}
static void check_permutation (const vector<size_t>& p, std::string name)
{
check_permutation (p.begin(), p.size(), name);
}
// -------------------------------------------------------------------
// triangle : from recursive to lattice numbering style
// -------------------------------------------------------------------
// simple lattice2d class function: (i,j) -> point(i,j)
struct lattice_simple {
lattice_simple() {}
point_basic<size_t> operator() (size_t i, size_t j) const {
return point_basic<size_t>(i,j);
}
};
// re-order internal nodes from gmsh/recursive to rheolef/left-right&bottom-top
// used also by 3d renum => template functions
template <class Lattice2d, class Function>
static
void
t_recursive_run (
size_t order,
size_t first_level,
size_t last_level,
const Lattice2d& ilat,
Function ilat2loc_inod,
const vector<size_t>& p,
vector<size_t>& msh_inod2loc_inod,
size_t& msh_inod)
{
for (size_t level = first_level; level < last_level; level++) {
// -----------------------------
// 3 vertex-nodes at this level:
// -----------------------------
// lower-left vertex:
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, level))];
if (level == order - 2*level) continue;
// lower-right vertex:
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(order-2*level, level))];
// upper-left vertex:
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, order-2*level))];
// -----------------------------
// 3 edge-nodes at this level:
// -----------------------------
size_t first_s = level + 1;
size_t last_s = order - 2*level;
// lower edge:
for (size_t s = first_s; s < last_s; s++) {
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(s, level))];
}
// upper-right edge:
for (size_t s = first_s; s < last_s; s++) {
size_t s1 = last_s - 1 - (s - first_s);
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(s1, s))];
}
// left edge:
for (size_t s = first_s; s < last_s; s++) {
size_t s1 = last_s - 1 - (s - first_s);
msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, s1))];
}
}
}
static vector<size_t> t_msh_inod2loc_inod;
static void t_msh2geo_init_node_renum (size_t order)
{
static bool has_init = false;
if (has_init) return;
has_init = true;
size_t loc_nnod = (order+1)*(order+2)/2;
t_msh_inod2loc_inod.resize (loc_nnod);
vector<size_t> id (loc_nnod);
for (size_t loc_inod = 0; loc_inod < loc_nnod; loc_inod++) id [loc_inod] = loc_inod;
size_t n_recursion = 1 + order/3; // integer division
size_t msh_inod = 0;
t_recursive_run (order, 0, n_recursion, lattice_simple(), reference_element_t_ilat2loc_inod, id, t_msh_inod2loc_inod, msh_inod);
check_permutation (t_msh_inod2loc_inod, "t_msh_inod2loc_inod");
}
// -------------------------------------------------------------------
// quadrangle : from recursive to lattice numbering style
// -------------------------------------------------------------------
// re-order internal nodes from gmsh/recursive to rheolef/left-right&bottom-top
template <class Lattice2d, class Function>
static
void
q_recursive_run (
size_t order,
size_t first_level,
size_t last_level,
const Lattice2d& ilat,
Function ilat2loc_inod,
vector<size_t>& msh_inod2loc_inod,
size_t& msh_inod)
{
for (size_t level = first_level; level < last_level; level++) {
// -----------------------------
// 3 vertex-nodes at this level:
// -----------------------------
// lower-left vertex:
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, level));
if (level == order - level) continue;
// lower-right vertex:
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, level));
// upper-right vertex:
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, order-level));
// upper-left vertex:
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, order-level));
// -----------------------------
// 4 edge-nodes at this level:
// -----------------------------
size_t first_s = level + 1;
size_t last_s = order - level;
// lower edge:
for (size_t s = first_s; s < last_s; s++) {
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(s, level));
}
// right edge:
for (size_t s = first_s; s < last_s; s++) {
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, s));
}
// top edge:
for (size_t s = first_s; s < last_s; s++) {
size_t s1 = last_s - 1 - (s - first_s);
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(s1, order-level));
}
// left edge:
for (size_t s = first_s; s < last_s; s++) {
size_t s1 = last_s - 1 - (s - first_s);
msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, s1));
}
}
}
static vector<size_t> q_msh_inod2loc_inod;
static void q_msh2geo_init_node_renum (size_t order)
{
static bool has_init = false;
if (has_init) return;
has_init = true;
size_t loc_nnod = (order+1)*(order+1);
q_msh_inod2loc_inod.resize (loc_nnod);
size_t n_recursion = 1 + order/2; // integer division
size_t msh_inod = 0;
q_recursive_run (order, 0, n_recursion, lattice_simple(), reference_element_q_ilat2loc_inod, q_msh_inod2loc_inod, msh_inod);
check_permutation (q_msh_inod2loc_inod, "q_msh_inod2loc_inod");
}
// -------------------------------------------------------------------
// tetra : in two passes
// pass 1: re-index and re-orient edges & face
// pass 2: from recursive to lattice numbering style on face and in volume
// -------------------------------------------------------------------
// some 3d lattice class functions: (i,j) -> (I,J,K) in 3d tetra lattice
// tetra: face(02x01)
struct lattice_T_face_02x01 {
lattice_T_face_02x01 (size_t order) {}
point_basic<size_t> operator() (size_t i, size_t j) const {
return point_basic<size_t>(j,i,0);
}
};
// tetra: face(03x02)
struct lattice_T_face_03x02 {
lattice_T_face_03x02 (size_t order) {}
point_basic<size_t> operator() (size_t i, size_t j) const {
return point_basic<size_t>(0,j,i);
}
};
// tetra: face(01x03)
struct lattice_T_face_01x03 {
lattice_T_face_01x03 (size_t order) {}
point_basic<size_t> operator() (size_t i, size_t j) const {
return point_basic<size_t>(i,0,j);
}
};
// tetra: face(12x13)
struct lattice_T_face_12x13 {
lattice_T_face_12x13 (size_t order) : _n(order) {}
point_basic<size_t> operator() (size_t i, size_t j) const {
point_basic<size_t> x(_n-i-j,i,j);
return x;
}
size_t _n;
};
/* example of edge & face re-indexation
for the order=4 tetra (n_node=35)
range gmsh rheolef action
0-3 4 vert 4 vert ok
4-6 e01 e01 -
7-9 e12 e12 -
10-12 e20 e20 -
13-15 e30 e03 reverse orientation sign (rev)
16-18 e32 e13
19-21 e31 e23 swap 32 :=: 31 and rev orient of both
22-24 f021 f021 ok
25-27 f013 f032
28-30 f032 f013 swap 013 :=: 032
31-33 f312 f123 change origin & axis : 31x32 -> 12x13
34 vol vol -
*/
static vector<size_t> T_msh_inod2loc_inod;
static void T_one_level_run (size_t order, vector<size_t>::iterator msh_inod2loc_inod)
{
size_t n_edge_node = (order-1);
size_t n_face_node = (order-1)*(order-2)/2;
#ifdef TO_CLEAN
size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
#endif // TO_CLEAN
size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
if (order < 2) { // trivial renumbering
for (size_t loc_imsh = 0; loc_imsh < n_bdry_node; loc_imsh++) {
msh_inod2loc_inod [loc_imsh] = loc_imsh;
}
return;
}
vector<size_t> id (n_bdry_node);
for (size_t loc_inod = 0; loc_inod < n_bdry_node; loc_inod++) id [loc_inod] = loc_inod;
//
// 1) pass 1: change edges and faces idx: from msh 2 pmsh (permuted msh)
//
vector<size_t> msh_inod2pmsh_inod (n_bdry_node);
for (size_t msh_inod = 0; msh_inod < n_bdry_node; msh_inod++) msh_inod2pmsh_inod [msh_inod] = msh_inod;
// 1.1) swap edges-nodes between gmsh edge[4]=nodes(3,2) & edge[5]=nodes(3,1) and reverse orientation
// => becomes in rheolef 4=nodes(1,3), 5=nodes(2,3)
size_t a = 4, b = 5;
for (size_t k = 0; order >= 2 && k < order-1; k++) {
size_t msh_inod1 = 4+a*(order-1) + k;
size_t msh_inod2 = 4+b*(order-1) + (order-2-k);
msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
}
// 1.2) reverse orientation for gmsh edge 3 = nodes(3,0) => edge 3 = nodes(0,3)
a = 3; b = 3;
for (size_t k = 0; order >= 2 && k < order-1; k++) {
size_t msh_inod1 = 4+a*(order-1) + k;
size_t msh_inod2 = 4+b*(order-1) + (order-2-k);
msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
}
// 1.3) swap nodes-faces between face[1]=nodes(0,3,2) and face[2]=nodes(0,1,3)
a = 1; b = 2;
for (size_t k = 0; order >= 3 && k < (order-1)*(order-2)/2; k++) {
size_t msh_inod1 = 4 + 6*(order-1) + a*(order-1)*(order-2)/2 + k;
size_t msh_inod2 = 4 + 6*(order-1) + b*(order-1)*(order-2)/2 + k;
msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
}
check_permutation (msh_inod2pmsh_inod, "msh_inod2pmsh_inod");
//
// 2) pass 2: recursively renumber faces (order >= 5) : pmsh_inod2loc_inod
//
vector<size_t> pmsh_inod2loc_inod (n_bdry_node);
for (size_t pmsh_iloc = 0; pmsh_iloc < n_bdry_node; pmsh_iloc++) pmsh_inod2loc_inod [pmsh_iloc] = pmsh_iloc;
size_t n_recursion = 1 + order/3; // integer division
size_t pmsh_inod = 4 + 6*n_edge_node;
t_recursive_run (order, 1, n_recursion, lattice_T_face_02x01(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
t_recursive_run (order, 1, n_recursion, lattice_T_face_03x02(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
t_recursive_run (order, 1, n_recursion, lattice_T_face_01x03(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
t_recursive_run (order, 1, n_recursion, lattice_T_face_12x13(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
check_permutation (pmsh_inod2loc_inod, "pmsh_inod2loc_inod");
//
// 3) rotate: change origin & axis : 31x32 -> 12x13
//
vector<size_t> loc_inod2loc_inod2 (n_bdry_node);
for (size_t loc_inod = 0; loc_inod < n_bdry_node; loc_inod++) loc_inod2loc_inod2 [loc_inod] = loc_inod;
a = 3; // third face
for (size_t j = 0, n = order-2; j < n; j++) {
for (size_t i = 0; i + j < n; i++) {
// lattice rotation and origin change: (i,j) -> (i1,j1)
size_t i1 = j;
size_t j1 = (order-3) - i - j;
// internal face lattice idx (i,j) & (i1,j1) to face lattice idx (ig,jg) & (ir,jr):
size_t ig = i + 1;
size_t jg = j + 1;
size_t ir = i1 + 1;
size_t jr = j1 + 1;
// ir, then jr:
size_t jr1 = order - jr;
size_t ijr = (n_face_node - jr1*(jr1-1)/2) + (ir - 1);
size_t pmsh_inod = 4 + 6*n_edge_node + a*n_face_node + ijr;
// ig, then jg:
size_t jg1 = order - jg;
size_t ijg = (n_face_node - jg1*(jg1-1)/2) + (ig - 1);
size_t msh_inod = 4 + 6*n_edge_node + a*n_face_node + ijg;
loc_inod2loc_inod2 [msh_inod] = pmsh_inod;
}
}
check_permutation (loc_inod2loc_inod2, "loc_inod2loc_inod2");
//
// 4) compose msh_inod2pmsh_inod & pmsh_inod2loc_inod
//
for (size_t msh_inod = 0; msh_inod < n_bdry_node; msh_inod++) {
msh_inod2loc_inod [msh_inod] = loc_inod2loc_inod2[pmsh_inod2loc_inod[msh_inod2pmsh_inod[msh_inod]]];
}
check_permutation (msh_inod2loc_inod, n_bdry_node, "msh_inod2loc_inod");
}
static
void T_renum_as_lattice (size_t order, vector<size_t>::iterator loc_msh2loc_inod, size_t first_inod)
{
typedef point_basic<size_t> ilat;
size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
#ifdef TO_CLEAN
size_t n_edge_node = (order-1);
size_t n_face_node = (order-1)*(order-2)/2;
size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
#endif // TO_CLEAN
// 1) build inod2ilat permutation
vector<size_t> loc_inod2loc_ilat (loc_nnod);
for (size_t loc_ilat = 0, k = 0; k < order+1; k++) {
for (size_t j = 0; j+k < order+1; j++) {
for (size_t i = 0; i+j+k < order+1; i++) {
size_t loc_inod = reference_element_T_ilat2loc_inod (order, ilat(i, j, k));
loc_inod2loc_ilat [loc_inod] = loc_ilat++;
}
}
}
check_permutation (loc_inod2loc_ilat, "loc_inod2loc_ilat");
// 2) compose permutation
vector<size_t> loc_msh2loc_ilat (loc_nnod);
for (size_t loc_imsh = 0; loc_imsh < loc_nnod; loc_imsh++) {
loc_msh2loc_ilat [loc_imsh] = loc_inod2loc_ilat [loc_msh2loc_inod [loc_imsh] - first_inod] + first_inod;
}
// 3) replace it:
copy (loc_msh2loc_ilat.begin(), loc_msh2loc_ilat.end(), loc_msh2loc_inod);
}
static
void T_msh2geo_init_node_renum (size_t order)
{
static bool has_init = false;
if (has_init) return;
has_init = true;
size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
T_msh_inod2loc_inod.resize (loc_nnod);
for (size_t loc_imsh = 0; loc_imsh < loc_nnod; loc_imsh++) T_msh_inod2loc_inod [loc_imsh] = loc_imsh;
size_t n_level = (order + 4)/4; // integer div ; nb of imbricated tetrahedras
size_t first_loc_inod = 0; //
for (size_t level = 0; level < n_level; level++) {
size_t level_order = order - 4*level;
T_one_level_run (level_order, T_msh_inod2loc_inod.begin() + first_loc_inod);
// count nodes at this level
size_t n_edge_node = (level_order-1);
size_t n_face_node = (level_order-1)*(level_order-2)/2;
size_t n_level_node = (level_order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
// shift numbering
size_t last_loc_inod = first_loc_inod + n_level_node;
for (size_t loc_inod = first_loc_inod; loc_inod < last_loc_inod; loc_inod++) {
T_msh_inod2loc_inod [loc_inod] += first_loc_inod;
}
first_loc_inod += n_level_node;
}
// renum all internal nodes as lattice:
if (order >= 4) {
size_t n_edge_node = (order-1);
size_t n_face_node = (order-1)*(order-2)/2;
size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
T_renum_as_lattice (order-4, T_msh_inod2loc_inod.begin() + n_bdry_node, n_bdry_node);
}
check_permutation (T_msh_inod2loc_inod, "T_msh_inod2loc_inod");
}
// ---------------------------------------------------------------------------------
// the main renumbering function
// ---------------------------------------------------------------------------------
void
msh2geo_node_renum (vector<size_t>& element, size_t variant, size_t order)
{
if (order < 2) return;
vector<size_t> new_element (element.size());
copy (element.begin(), element.end(), new_element.begin());
if (variant == 't') {
t_msh2geo_init_node_renum (order);
check_macro (t_msh_inod2loc_inod.size() == element.size(), "invalid element size");
for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
new_element [t_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
}
} else if (variant == 'q') {
q_msh2geo_init_node_renum (order);
check_macro (q_msh_inod2loc_inod.size() == element.size(), "invalid element size");
for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
new_element [q_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
}
} else if (variant == 'T') {
T_msh2geo_init_node_renum (order);
check_macro (T_msh_inod2loc_inod.size() == element.size(), "invalid element size");
for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
new_element [T_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
}
} else if (variant == 'H') {
check_macro (order <= 2, "unsupported hexaedron order > 2 element");
if (order == 2) {
static size_t H_p2_msh_inod2loc_inod [27] = {
0, 1, 2, 3, 4, 5, 6, 7, // vertices unchanged
8, 11, 12, 9, 13, 10, 14, 15, 16, 19, 17, 18, // edges permuted
20, 22, 21, 24, 25, 23, // faces permuted
26 }; // barycenter unchanged
check_macro (27 == element.size(), "invalid element size");
check_permutation (H_p2_msh_inod2loc_inod, 27, "H_p2_msh_inod2loc_inod");
for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
new_element [H_p2_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
}
}
}
copy (new_element.begin(), new_element.end(), element.begin());
}
} // namespace rheolef
rheolef::lattice_T_face_02x01::operator()
point_basic< size_t > operator()(size_t i, size_t j) const
Definition: msh2geo_node_renum.icc:216
rheolef::lattice_T_face_02x01
Definition: msh2geo_node_renum.icc:214
rheolef::point_basic< size_t >
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::t_msh_inod2loc_inod
static vector< size_t > t_msh_inod2loc_inod
Definition: msh2geo_node_renum.icc:120
rheolef::lattice_simple::lattice_simple
lattice_simple()
Definition: msh2geo_node_renum.icc:68
rheolef::Function
rheolef::std Function
rheolef::T_msh_inod2loc_inod
static vector< size_t > T_msh_inod2loc_inod
Definition: msh2geo_node_renum.icc:268
rheolef::lattice_T_face_03x02::lattice_T_face_03x02
lattice_T_face_03x02(size_t order)
Definition: msh2geo_node_renum.icc:222
rheolef::lattice_T_face_12x13
Definition: msh2geo_node_renum.icc:235
rheolef::lattice_T_face_01x03::lattice_T_face_01x03
lattice_T_face_01x03(size_t order)
Definition: msh2geo_node_renum.icc:229
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::lattice_simple
Definition: msh2geo_node_renum.icc:67
rheolef::msh2geo_node_renum
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
Definition: msh2geo_node_renum.icc:447
p
Definition: sphere.icc:25
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
a
Definition: diffusion_isotropic.h:25
rheolef::lattice_T_face_01x03::operator()
point_basic< size_t > operator()(size_t i, size_t j) const
Definition: msh2geo_node_renum.icc:230
rheolef::lattice_T_face_03x02::operator()
point_basic< size_t > operator()(size_t i, size_t j) const
Definition: msh2geo_node_renum.icc:223
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::lattice_T_face_02x01::lattice_T_face_02x01
lattice_T_face_02x01(size_t order)
Definition: msh2geo_node_renum.icc:215
rheolef::lattice_T_face_12x13::_n
size_t _n
Definition: msh2geo_node_renum.icc:241
rheolef::lattice_T_face_12x13::operator()
point_basic< size_t > operator()(size_t i, size_t j) const
Definition: msh2geo_node_renum.icc:237
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
rheolef::lattice_simple::operator()
point_basic< size_t > operator()(size_t i, size_t j) const
Definition: msh2geo_node_renum.icc:69
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::lattice_T_face_03x02
Definition: msh2geo_node_renum.icc:221
rheolef::lattice_T_face_01x03
Definition: msh2geo_node_renum.icc:228
rheolef::lattice_T_face_12x13::lattice_T_face_12x13
lattice_T_face_12x13(size_t order)
Definition: msh2geo_node_renum.icc:236
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::q_msh_inod2loc_inod
static vector< size_t > q_msh_inod2loc_inod
Definition: msh2geo_node_renum.icc:191