35template<
class Iterator>
36static 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++) {
56 check_macro (ok,
"permutation("<<n<<
"): "<<name<<
" has problem");
58static void check_permutation (
const vector<size_t>&
p, std::string name)
60 check_permutation (
p.begin(),
p.size(), name);
75template <
class Lattice2d,
class Function>
82 const Lattice2d& ilat,
83 Function ilat2loc_inod,
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))];
121static void t_msh2geo_init_node_renum (
size_t order)
123 static bool has_init =
false;
124 if (has_init)
return;
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;
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");
141template <
class Lattice2d,
class Function>
148 const Lattice2d& ilat,
149 Function ilat2loc_inod,
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));
192static void q_msh2geo_init_node_renum (
size_t order)
194 static bool has_init =
false;
195 if (has_init)
return;
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;
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");
270static void T_one_level_run (
size_t order, vector<size_t>::iterator msh_inod2loc_inod)
272 size_t n_edge_node = (order-1);
273 size_t n_face_node = (order-1)*(order-2)/2;
275 size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
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;
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;
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;
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;
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");
374void 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;
378 size_t loc_nnod = (order+1)*(order+2)*(order+3)/6;
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;
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);
406void T_msh2geo_init_node_renum (
size_t order)
408 static bool has_init =
false;
409 if (has_init)
return;
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;
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;
420 T_one_level_run (level_order, T_msh_inod2loc_inod.begin() + first_loc_inod);
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++) {
430 T_msh_inod2loc_inod [loc_inod] += first_loc_inod;
432 first_loc_inod += n_level_node;
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);
441 check_permutation (T_msh_inod2loc_inod,
"T_msh_inod2loc_inod");
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];
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];
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];
470 }
else if (variant ==
'H') {
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
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)
lattice_T_face_01x03(size_t order)
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const
lattice_T_face_02x01(size_t order)
point_basic< size_t > operator()(size_t i, size_t j) const
lattice_T_face_03x02(size_t order)
lattice_T_face_12x13(size_t order)
point_basic< size_t > operator()(size_t i, size_t j) const
point_basic< size_t > operator()(size_t i, size_t j) const