Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
msh2geo_node_renum.icc
Go to the documentation of this file.
1
21//
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//
33namespace rheolef {
34
35template<class Iterator>
36static 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}
58static 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
75template <class Lattice2d, class Function>
76static
77void
78t_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}
120static vector<size_t> t_msh_inod2loc_inod;
121static 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
141template <class Lattice2d, class Function>
142static
143void
144q_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}
191static vector<size_t> q_msh_inod2loc_inod;
192static 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)
215 lattice_T_face_02x01 (size_t order) {}
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)
222 lattice_T_face_03x02 (size_t order) {}
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)
229 lattice_T_face_01x03 (size_t order) {}
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)
236 lattice_T_face_12x13 (size_t order) : _n(order) {}
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
246range gmsh rheolef action
247
2480-3 4 vert 4 vert ok
2494-6 e01 e01 -
2507-9 e12 e12 -
25110-12 e20 e20 -
252
25313-15 e30 e03 reverse orientation sign (rev)
254
25516-18 e32 e13
25619-21 e31 e23 swap 32 :=: 31 and rev orient of both
257
25822-24 f021 f021 ok
259
26025-27 f013 f032
26128-30 f032 f013 swap 013 :=: 032
262
26331-33 f312 f123 change origin & axis : 31x32 -> 12x13
264
26534 vol vol -
266
267*/
268static vector<size_t> T_msh_inod2loc_inod;
269
270static 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}
373static
374void 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}
405static
406void 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// ---------------------------------------------------------------------------------
446void
447msh2geo_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
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