35 template<
class Iterator>
36 static void check_permutation (Iterator
p,
size_t n, std::string
name)
39 for (
size_t i = 0; i <
n; i++) {
44 vector<size_t> ip (
n, numeric_limits<size_t>::max());
45 for (
size_t i = 0; i <
n; i++) {
51 for (
size_t i = 0; i <
n; i++) {
58 static void check_permutation (
const vector<size_t>&
p, std::string
name)
60 check_permutation (
p.begin(),
p.size(),
name);
75 template <
class Lattice2d,
class Function>
82 const Lattice2d& ilat,
84 const vector<size_t>&
p,
85 vector<size_t>& msh_inod2loc_inod,
88 for (
size_t level = first_level; level < last_level; level++) {
93 msh_inod2loc_inod [msh_inod++] =
p[ilat2loc_inod (
order, ilat(level, level))];
94 if (level ==
order - 2*level)
continue;
96 msh_inod2loc_inod [msh_inod++] =
p[ilat2loc_inod (
order, ilat(
order-2*level, level))];
98 msh_inod2loc_inod [msh_inod++] =
p[ilat2loc_inod (
order, ilat(level,
order-2*level))];
102 size_t first_s = level + 1;
103 size_t last_s =
order - 2*level;
105 for (
size_t s = first_s; s < last_s; s++) {
106 msh_inod2loc_inod [msh_inod++] =
p[ilat2loc_inod (
order, ilat(s, level))];
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))];
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))];
121 static void t_msh2geo_init_node_renum (
size_t order)
123 static bool has_init =
false;
124 if (has_init)
return;
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;
133 t_recursive_run (
order, 0, n_recursion, lattice_simple(), reference_element_t_ilat2loc_inod,
id,
t_msh_inod2loc_inod, msh_inod);
141 template <
class Lattice2d,
class Function>
148 const Lattice2d& ilat,
150 vector<size_t>& msh_inod2loc_inod,
153 for (
size_t level = first_level; level < last_level; level++) {
158 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(level, level));
159 if (level ==
order - level)
continue;
161 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(
order-level, level));
163 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(
order-level,
order-level));
165 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(level,
order-level));
169 size_t first_s = level + 1;
170 size_t last_s =
order - level;
172 for (
size_t s = first_s; s < last_s; s++) {
173 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(s, level));
176 for (
size_t s = first_s; s < last_s; s++) {
177 msh_inod2loc_inod [msh_inod++] = ilat2loc_inod (
order, ilat(
order-level, s));
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));
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));
192 static void q_msh2geo_init_node_renum (
size_t order)
194 static bool has_init =
false;
195 if (has_init)
return;
200 size_t n_recursion = 1 +
order/2;
202 q_recursive_run (
order, 0, n_recursion, lattice_simple(), reference_element_q_ilat2loc_inod,
q_msh_inod2loc_inod, msh_inod);
270 static void T_one_level_run (
size_t order, vector<size_t>::iterator msh_inod2loc_inod)
272 size_t n_edge_node = (
order-1);
277 size_t n_bdry_node = (
order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
279 for (
size_t loc_imsh = 0; loc_imsh < n_bdry_node; loc_imsh++) {
280 msh_inod2loc_inod [loc_imsh] = loc_imsh;
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;
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;
294 for (
size_t k = 0;
order >= 2 && k <
order-1; k++) {
295 size_t msh_inod1 = 4+
a*(
order-1) + k;
297 msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
298 msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
302 for (
size_t k = 0;
order >= 2 && k <
order-1; k++) {
303 size_t msh_inod1 = 4+
a*(
order-1) + k;
305 msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
312 msh_inod2pmsh_inod[msh_inod1] = msh_inod2;
313 msh_inod2pmsh_inod[msh_inod2] = msh_inod1;
315 check_permutation (msh_inod2pmsh_inod,
"msh_inod2pmsh_inod");
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;
322 size_t pmsh_inod = 4 + 6*n_edge_node;
324 t_recursive_run (
order, 1, n_recursion, lattice_T_face_02x01(
order), reference_element_T_ilat2loc_inod,
id, pmsh_inod2loc_inod, pmsh_inod);
326 t_recursive_run (
order, 1, n_recursion, lattice_T_face_03x02(
order), reference_element_T_ilat2loc_inod,
id, pmsh_inod2loc_inod, pmsh_inod);
328 t_recursive_run (
order, 1, n_recursion, lattice_T_face_01x03(
order), reference_element_T_ilat2loc_inod,
id, pmsh_inod2loc_inod, pmsh_inod);
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");
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;
338 for (
size_t j = 0,
n =
order-2; j <
n; j++) {
339 for (
size_t i = 0; i + j <
n; i++) {
343 size_t j1 = (
order-3) - i - j;
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;
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;
361 loc_inod2loc_inod2 [msh_inod] = pmsh_inod;
364 check_permutation (loc_inod2loc_inod2,
"loc_inod2loc_inod2");
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]]];
371 check_permutation (msh_inod2loc_inod, n_bdry_node,
"msh_inod2loc_inod");
374 void T_renum_as_lattice (
size_t order, vector<size_t>::iterator loc_msh2loc_inod,
size_t first_inod)
376 typedef point_basic<size_t> ilat;
380 size_t n_edge_node = (
order-1);
382 size_t n_bdry_node = (
order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
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++;
395 check_permutation (loc_inod2loc_ilat,
"loc_inod2loc_ilat");
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;
403 copy (loc_msh2loc_ilat.begin(), loc_msh2loc_ilat.end(), loc_msh2loc_inod);
406 void T_msh2geo_init_node_renum (
size_t order)
408 static bool has_init =
false;
409 if (has_init)
return;
414 for (
size_t loc_imsh = 0; loc_imsh < loc_nnod; loc_imsh++)
T_msh_inod2loc_inod [loc_imsh] = loc_imsh;
416 size_t n_level = (
order + 4)/4;
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;
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;
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++) {
432 first_loc_inod += n_level_node;
436 size_t n_edge_node = (
order-1);
438 size_t n_bdry_node = (
order == 0) ? 1 : 4 + 6*n_edge_node + 4*n_face_node;
449 if (
order < 2)
return;
450 vector<size_t> new_element (element.size());
451 copy (element.begin(), element.end(), new_element.begin());
453 t_msh2geo_init_node_renum (
order);
455 for (
size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
459 q_msh2geo_init_node_renum (
order);
461 for (
size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
465 T_msh2geo_init_node_renum (
order);
467 for (
size_t msh_inod = 0; msh_inod < element.size(); msh_inod++) {
471 check_macro (order <= 2, "unsupported hexaedron order > 2 element
");
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];
485 copy (new_element.begin(), new_element.end(), element.begin());
488 } // namespace rheolef