Rheolef  7.2
an efficient C++ finite element environment
msh2geo_node_renum.icc
Go to the documentation of this file.
1 //
22 // for high order gmsh files:
23 // reorder the node list from gmsh 2 rheolef conventions:
24 // - d=2,3 : for face & volume internal nodes,
25 // change from recursive gmsh numbering style
26 // to rheolef lattice style numbering
27 // - d=3 : change edge & face local index and orientation
28 //
29 // author: Pierre.Saramito@imag.fr
30 //
31 // date: 18 september 2011
32 //
33 namespace rheolef {
34 
35 template<class Iterator>
36 static void check_permutation (Iterator p, size_t n, std::string name)
37 {
38  bool ok = true;
39  for (size_t i = 0; i < n; i++) {
40  if (p[i] >= n) {
41  ok = false;
42  }
43  }
44  vector<size_t> ip (n, numeric_limits<size_t>::max());
45  for (size_t i = 0; i < n; i++) {
46  if (ip[p[i]] < n) {
47  ok = false;
48  }
49  ip[p[i]] = i;
50  }
51  for (size_t i = 0; i < n; i++) {
52  if (ip[p[i]] != i) {
53  ok = false;
54  }
55  }
56  check_macro (ok, "permutation("<<n<<"): "<<name<<" has problem");
57 }
58 static void check_permutation (const vector<size_t>& p, std::string name)
59 {
60  check_permutation (p.begin(), p.size(), name);
61 }
62 // -------------------------------------------------------------------
63 // triangle : from recursive to lattice numbering style
64 // -------------------------------------------------------------------
65 
66 // simple lattice2d class function: (i,j) -> point(i,j)
69  point_basic<size_t> operator() (size_t i, size_t j) const {
70  return point_basic<size_t>(i,j);
71  }
72 };
73 // re-order internal nodes from gmsh/recursive to rheolef/left-right&bottom-top
74 // used also by 3d renum => template functions
75 template <class Lattice2d, class Function>
76 static
77 void
78 t_recursive_run (
79  size_t order,
80  size_t first_level,
81  size_t last_level,
82  const Lattice2d& ilat,
83  Function ilat2loc_inod,
84  const vector<size_t>& p,
85  vector<size_t>& msh_inod2loc_inod,
86  size_t& msh_inod)
87 {
88  for (size_t level = first_level; level < last_level; level++) {
89  // -----------------------------
90  // 3 vertex-nodes at this level:
91  // -----------------------------
92  // lower-left vertex:
93  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, level))];
94  if (level == order - 2*level) continue;
95  // lower-right vertex:
96  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(order-2*level, level))];
97  // upper-left vertex:
98  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, order-2*level))];
99  // -----------------------------
100  // 3 edge-nodes at this level:
101  // -----------------------------
102  size_t first_s = level + 1;
103  size_t last_s = order - 2*level;
104  // lower edge:
105  for (size_t s = first_s; s < last_s; s++) {
106  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(s, level))];
107  }
108  // upper-right edge:
109  for (size_t s = first_s; s < last_s; s++) {
110  size_t s1 = last_s - 1 - (s - first_s);
111  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(s1, s))];
112  }
113  // left edge:
114  for (size_t s = first_s; s < last_s; s++) {
115  size_t s1 = last_s - 1 - (s - first_s);
116  msh_inod2loc_inod [msh_inod++] = p[ilat2loc_inod (order, ilat(level, s1))];
117  }
118  }
119 }
120 static vector<size_t> t_msh_inod2loc_inod;
121 static void t_msh2geo_init_node_renum (size_t order)
122 {
123  static bool has_init = false;
124  if (has_init) return;
125  has_init = true;
126 
127  size_t loc_nnod = (order+1)*(order+2)/2;
128  t_msh_inod2loc_inod.resize (loc_nnod);
129  vector<size_t> id (loc_nnod);
130  for (size_t loc_inod = 0; loc_inod < loc_nnod; loc_inod++) id [loc_inod] = loc_inod;
131  size_t n_recursion = 1 + order/3; // integer division
132  size_t msh_inod = 0;
133  t_recursive_run (order, 0, n_recursion, lattice_simple(), reference_element_t_ilat2loc_inod, id, t_msh_inod2loc_inod, msh_inod);
134  check_permutation (t_msh_inod2loc_inod, "t_msh_inod2loc_inod");
135 }
136 // -------------------------------------------------------------------
137 // quadrangle : from recursive to lattice numbering style
138 // -------------------------------------------------------------------
139 
140 // re-order internal nodes from gmsh/recursive to rheolef/left-right&bottom-top
141 template <class Lattice2d, class Function>
142 static
143 void
144 q_recursive_run (
145  size_t order,
146  size_t first_level,
147  size_t last_level,
148  const Lattice2d& ilat,
149  Function ilat2loc_inod,
150  vector<size_t>& msh_inod2loc_inod,
151  size_t& msh_inod)
152 {
153  for (size_t level = first_level; level < last_level; level++) {
154  // -----------------------------
155  // 3 vertex-nodes at this level:
156  // -----------------------------
157  // lower-left vertex:
158  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, level));
159  if (level == order - level) continue;
160  // lower-right vertex:
161  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, level));
162  // upper-right vertex:
163  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, order-level));
164  // upper-left vertex:
165  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, order-level));
166  // -----------------------------
167  // 4 edge-nodes at this level:
168  // -----------------------------
169  size_t first_s = level + 1;
170  size_t last_s = order - level;
171  // lower edge:
172  for (size_t s = first_s; s < last_s; s++) {
173  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(s, level));
174  }
175  // right edge:
176  for (size_t s = first_s; s < last_s; s++) {
177  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(order-level, s));
178  }
179  // top edge:
180  for (size_t s = first_s; s < last_s; s++) {
181  size_t s1 = last_s - 1 - (s - first_s);
182  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(s1, order-level));
183  }
184  // left edge:
185  for (size_t s = first_s; s < last_s; s++) {
186  size_t s1 = last_s - 1 - (s - first_s);
187  msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (order, ilat(level, s1));
188  }
189  }
190 }
191 static vector<size_t> q_msh_inod2loc_inod;
192 static void q_msh2geo_init_node_renum (size_t order)
193 {
194  static bool has_init = false;
195  if (has_init) return;
196  has_init = true;
197 
198  size_t loc_nnod = (order+1)*(order+1);
199  q_msh_inod2loc_inod.resize (loc_nnod);
200  size_t n_recursion = 1 + order/2; // integer division
201  size_t msh_inod = 0;
202  q_recursive_run (order, 0, n_recursion, lattice_simple(), reference_element_q_ilat2loc_inod, q_msh_inod2loc_inod, msh_inod);
203  check_permutation (q_msh_inod2loc_inod, "q_msh_inod2loc_inod");
204 }
205 // -------------------------------------------------------------------
206 // tetra : in two passes
207 // pass 1: re-index and re-orient edges & face
208 // pass 2: from recursive to lattice numbering style on face and in volume
209 // -------------------------------------------------------------------
210 
211 
212 // some 3d lattice class functions: (i,j) -> (I,J,K) in 3d tetra lattice
213 // tetra: face(02x01)
216  point_basic<size_t> operator() (size_t i, size_t j) const {
217  return point_basic<size_t>(j,i,0);
218  }
219 };
220 // tetra: face(03x02)
223  point_basic<size_t> operator() (size_t i, size_t j) const {
224  return point_basic<size_t>(0,j,i);
225  }
226 };
227 // tetra: face(01x03)
230  point_basic<size_t> operator() (size_t i, size_t j) const {
231  return point_basic<size_t>(i,0,j);
232  }
233 };
234 // tetra: face(12x13)
237  point_basic<size_t> operator() (size_t i, size_t j) const {
238  point_basic<size_t> x(_n-i-j,i,j);
239  return x;
240  }
241  size_t _n;
242 };
243 /* example of edge & face re-indexation
244  for the order=4 tetra (n_node=35)
245 
246 range gmsh rheolef action
247 
248 0-3 4 vert 4 vert ok
249 4-6 e01 e01 -
250 7-9 e12 e12 -
251 10-12 e20 e20 -
252 
253 13-15 e30 e03 reverse orientation sign (rev)
254 
255 16-18 e32 e13
256 19-21 e31 e23 swap 32 :=: 31 and rev orient of both
257 
258 22-24 f021 f021 ok
259 
260 25-27 f013 f032
261 28-30 f032 f013 swap 013 :=: 032
262 
263 31-33 f312 f123 change origin & axis : 31x32 -> 12x13
264 
265 34 vol vol -
266 
267 */
268 static vector<size_t> T_msh_inod2loc_inod;
269 
270 static void T_one_level_run (size_t order, vector<size_t>::iterator msh_inod2loc_inod)
271 {
272  size_t n_edge_node = (order-1);
273  size_t n_face_node = (order-1)*(order-2)/2;
274 #ifdef TO_CLEAN
275  size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
276 #endif // TO_CLEAN
277  size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
278  if (order < 2) { // trivial renumbering
279  for (size_t loc_imsh = 0; loc_imsh < n_bdry_node; loc_imsh++) {
280  msh_inod2loc_inod [loc_imsh] = loc_imsh;
281  }
282  return;
283  }
284  vector<size_t> id (n_bdry_node);
285  for (size_t loc_inod = 0; loc_inod < n_bdry_node; loc_inod++) id [loc_inod] = loc_inod;
286  //
287  // 1) pass 1: change edges and faces idx: from msh 2 pmsh (permuted msh)
288  //
289  vector<size_t> msh_inod2pmsh_inod (n_bdry_node);
290  for (size_t msh_inod = 0; msh_inod < n_bdry_node; msh_inod++) msh_inod2pmsh_inod [msh_inod] = msh_inod;
291  // 1.1) swap edges-nodes between gmsh edge[4]=nodes(3,2) & edge[5]=nodes(3,1) and reverse orientation
292  // => becomes in rheolef 4=nodes(1,3), 5=nodes(2,3)
293  size_t a = 4, b = 5;
294  for (size_t k = 0; order >= 2 && k < order-1; k++) {
295  size_t msh_inod1 = 4+a*(order-1) + k;
296  size_t msh_inod2 = 4+b*(order-1) + (order-2-k);
297  msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
298  msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
299  }
300  // 1.2) reverse orientation for gmsh edge 3 = nodes(3,0) => edge 3 = nodes(0,3)
301  a = 3; b = 3;
302  for (size_t k = 0; order >= 2 && k < order-1; k++) {
303  size_t msh_inod1 = 4+a*(order-1) + k;
304  size_t msh_inod2 = 4+b*(order-1) + (order-2-k);
305  msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
306  }
307  // 1.3) swap nodes-faces between face[1]=nodes(0,3,2) and face[2]=nodes(0,1,3)
308  a = 1; b = 2;
309  for (size_t k = 0; order >= 3 && k < (order-1)*(order-2)/2; k++) {
310  size_t msh_inod1 = 4 + 6*(order-1) + a*(order-1)*(order-2)/2 + k;
311  size_t msh_inod2 = 4 + 6*(order-1) + b*(order-1)*(order-2)/2 + k;
312  msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
313  msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
314  }
315  check_permutation (msh_inod2pmsh_inod, "msh_inod2pmsh_inod");
316  //
317  // 2) pass 2: recursively renumber faces (order >= 5) : pmsh_inod2loc_inod
318  //
319  vector<size_t> pmsh_inod2loc_inod (n_bdry_node);
320  for (size_t pmsh_iloc = 0; pmsh_iloc < n_bdry_node; pmsh_iloc++) pmsh_inod2loc_inod [pmsh_iloc] = pmsh_iloc;
321  size_t n_recursion = 1 + order/3; // integer division
322  size_t pmsh_inod = 4 + 6*n_edge_node;
323 
324  t_recursive_run (order, 1, n_recursion, lattice_T_face_02x01(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
325 
326  t_recursive_run (order, 1, n_recursion, lattice_T_face_03x02(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
327 
328  t_recursive_run (order, 1, n_recursion, lattice_T_face_01x03(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
329 
330  t_recursive_run (order, 1, n_recursion, lattice_T_face_12x13(order), reference_element_T_ilat2loc_inod, id, pmsh_inod2loc_inod, pmsh_inod);
331  check_permutation (pmsh_inod2loc_inod, "pmsh_inod2loc_inod");
332  //
333  // 3) rotate: change origin & axis : 31x32 -> 12x13
334  //
335  vector<size_t> loc_inod2loc_inod2 (n_bdry_node);
336  for (size_t loc_inod = 0; loc_inod < n_bdry_node; loc_inod++) loc_inod2loc_inod2 [loc_inod] = loc_inod;
337  a = 3; // third face
338  for (size_t j = 0, n = order-2; j < n; j++) {
339  for (size_t i = 0; i + j < n; i++) {
340 
341  // lattice rotation and origin change: (i,j) -> (i1,j1)
342  size_t i1 = j;
343  size_t j1 = (order-3) - i - j;
344 
345  // internal face lattice idx (i,j) & (i1,j1) to face lattice idx (ig,jg) & (ir,jr):
346  size_t ig = i + 1;
347  size_t jg = j + 1;
348  size_t ir = i1 + 1;
349  size_t jr = j1 + 1;
350 
351  // ir, then jr:
352  size_t jr1 = order - jr;
353  size_t ijr = (n_face_node - jr1*(jr1-1)/2) + (ir - 1);
354  size_t pmsh_inod = 4 + 6*n_edge_node + a*n_face_node + ijr;
355 
356  // ig, then jg:
357  size_t jg1 = order - jg;
358  size_t ijg = (n_face_node - jg1*(jg1-1)/2) + (ig - 1);
359  size_t msh_inod = 4 + 6*n_edge_node + a*n_face_node + ijg;
360 
361  loc_inod2loc_inod2 [msh_inod] = pmsh_inod;
362  }
363  }
364  check_permutation (loc_inod2loc_inod2, "loc_inod2loc_inod2");
365  //
366  // 4) compose msh_inod2pmsh_inod & pmsh_inod2loc_inod
367  //
368  for (size_t msh_inod = 0; msh_inod < n_bdry_node; msh_inod++) {
369  msh_inod2loc_inod [msh_inod] = loc_inod2loc_inod2[pmsh_inod2loc_inod[msh_inod2pmsh_inod[msh_inod]]];
370  }
371  check_permutation (msh_inod2loc_inod, n_bdry_node, "msh_inod2loc_inod");
372 }
373 static
374 void T_renum_as_lattice (size_t order, vector<size_t>::iterator loc_msh2loc_inod, size_t first_inod)
375 {
376  typedef point_basic<size_t> ilat;
377 
378  size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
379 #ifdef TO_CLEAN
380  size_t n_edge_node = (order-1);
381  size_t n_face_node = (order-1)*(order-2)/2;
382  size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
383 #endif // TO_CLEAN
384 
385  // 1) build inod2ilat permutation
386  vector<size_t> loc_inod2loc_ilat (loc_nnod);
387  for (size_t loc_ilat = 0, k = 0; k < order+1; k++) {
388  for (size_t j = 0; j+k < order+1; j++) {
389  for (size_t i = 0; i+j+k < order+1; i++) {
390  size_t loc_inod = reference_element_T_ilat2loc_inod (order, ilat(i, j, k));
391  loc_inod2loc_ilat [loc_inod] = loc_ilat++;
392  }
393  }
394  }
395  check_permutation (loc_inod2loc_ilat, "loc_inod2loc_ilat");
396 
397  // 2) compose permutation
398  vector<size_t> loc_msh2loc_ilat (loc_nnod);
399  for (size_t loc_imsh = 0; loc_imsh < loc_nnod; loc_imsh++) {
400  loc_msh2loc_ilat [loc_imsh] = loc_inod2loc_ilat [loc_msh2loc_inod [loc_imsh] - first_inod] + first_inod;
401  }
402  // 3) replace it:
403  copy (loc_msh2loc_ilat.begin(), loc_msh2loc_ilat.end(), loc_msh2loc_inod);
404 }
405 static
406 void T_msh2geo_init_node_renum (size_t order)
407 {
408  static bool has_init = false;
409  if (has_init) return;
410  has_init = true;
411 
412  size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
413  T_msh_inod2loc_inod.resize (loc_nnod);
414  for (size_t loc_imsh = 0; loc_imsh < loc_nnod; loc_imsh++) T_msh_inod2loc_inod [loc_imsh] = loc_imsh;
415 
416  size_t n_level = (order + 4)/4; // integer div ; nb of imbricated tetrahedras
417  size_t first_loc_inod = 0; //
418  for (size_t level = 0; level < n_level; level++) {
419  size_t level_order = order - 4*level;
420  T_one_level_run (level_order, T_msh_inod2loc_inod.begin() + first_loc_inod);
421 
422  // count nodes at this level
423  size_t n_edge_node = (level_order-1);
424  size_t n_face_node = (level_order-1)*(level_order-2)/2;
425  size_t n_level_node = (level_order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
426 
427  // shift numbering
428  size_t last_loc_inod = first_loc_inod + n_level_node;
429  for (size_t loc_inod = first_loc_inod; loc_inod < last_loc_inod; loc_inod++) {
430  T_msh_inod2loc_inod [loc_inod] += first_loc_inod;
431  }
432  first_loc_inod += n_level_node;
433  }
434  // renum all internal nodes as lattice:
435  if (order >= 4) {
436  size_t n_edge_node = (order-1);
437  size_t n_face_node = (order-1)*(order-2)/2;
438  size_t n_bdry_node = (order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
439  T_renum_as_lattice (order-4, T_msh_inod2loc_inod.begin() + n_bdry_node, n_bdry_node);
440  }
441  check_permutation (T_msh_inod2loc_inod, "T_msh_inod2loc_inod");
442 }
443 // ---------------------------------------------------------------------------------
444 // the main renumbering function
445 // ---------------------------------------------------------------------------------
446 void
447 msh2geo_node_renum (vector<size_t>& element, size_t variant, size_t order)
448 {
449  if (order < 2) return;
450  vector<size_t> new_element (element.size());
451  copy (element.begin(), element.end(), new_element.begin());
452  if (variant == 't') {
453  t_msh2geo_init_node_renum (order);
454  check_macro (t_msh_inod2loc_inod.size() == element.size(), "invalid element size");
455  for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
456  new_element [t_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
457  }
458  } else if (variant == 'q') {
459  q_msh2geo_init_node_renum (order);
460  check_macro (q_msh_inod2loc_inod.size() == element.size(), "invalid element size");
461  for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
462  new_element [q_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
463  }
464  } else if (variant == 'T') {
465  T_msh2geo_init_node_renum (order);
466  check_macro (T_msh_inod2loc_inod.size() == element.size(), "invalid element size");
467  for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
468  new_element [T_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
469  }
470  } else if (variant == 'H') {
471  check_macro (order <= 2, "unsupported hexaedron order > 2 element");
472  if (order == 2) {
473  static size_t H_p2_msh_inod2loc_inod [27] = {
474  0, 1, 2, 3, 4, 5, 6, 7, // vertices unchanged
475  8, 11, 12, 9, 13, 10, 14, 15, 16, 19, 17, 18, // edges permuted
476  20, 22, 21, 24, 25, 23, // faces permuted
477  26 }; // barycenter unchanged
478  check_macro (27 == element.size(), "invalid element size");
479  check_permutation (H_p2_msh_inod2loc_inod, 27, "H_p2_msh_inod2loc_inod");
480  for (size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
481  new_element [H_p2_msh_inod2loc_inod[msh_inod]] = element[msh_inod];
482  }
483  }
484  }
485  copy (new_element.begin(), new_element.end(), element.begin());
486 }
487 
488 } // namespace rheolef
static vector< size_t > q_msh_inod2loc_inod
static vector< size_t > t_msh_inod2loc_inod
static vector< size_t > T_msh_inod2loc_inod
rheolef::std Function
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
Definition: sphere.icc:25
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const