1#ifndef _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
2#define _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
35template <
typename GeoElement,
typename OrientationType,
typename ShiftType>
38geo_element_get_orientation_and_shift (
40 size_t dis_iv0,
size_t dis_iv1,
size_t dis_iv2,
41 OrientationType& orient,
44 if (S[0] == dis_iv0) {
45 orient = (S[1] == dis_iv1) ? 1 : -1;
47 }
else if (S[0] == dis_iv1) {
48 orient = (S[1] == dis_iv2) ? 1 : -1;
51 orient = (S[1] == dis_iv0) ? 1 : -1;
57template <
typename GeoElement,
typename OrientationType,
typename ShiftType>
60geo_element_get_orientation_and_shift (
62 size_t dis_iv0,
size_t dis_iv1,
size_t dis_iv2,
size_t dis_iv3,
63 OrientationType& orient,
66 if (S[0] == dis_iv0) {
67 orient = (S[1] == dis_iv1) ? 1 : -1;
69 }
else if (S[0] == dis_iv1) {
70 orient = (S[1] == dis_iv2) ? 1 : -1;
72 }
else if (S[0] == dis_iv2) {
73 orient = (S[1] == dis_iv3) ? 1 : -1;
76 orient = (S[1] == dis_iv0) ? 1 : -1;
88geo_element_fix_edge_indirect (
94 if (K.dimension() < 2 || order <= 2 || K.edge_indirect(loc_iedg).orientation() > 0)
return loc_iedg_j;
95 return order - 2 - loc_iedg_j;
102geo_element_loc_tri_inod2lattice (
105 point_basic<size_t>& ij_lattice)
112 size_t n_face_node = (order-1)*(order-2)/2;
115 double c = -2*int(n_face_node - loc_tri_inod);
116 double delta = b*b - 4*a*c;
117 check_macro (
delta >= 0,
"loc_tri_inod2lattice(loc_tri_inod="<<loc_tri_inod<<
",order="<<order<<
"): unexpected delta="<<
delta<<
" < 0");
118 double j1_plus = (-b + sqrt(
delta))/(2*a);
120 size_t j = floor (order - j1_plus);
121 size_t j1 = order - j;
122 size_t i = loc_tri_inod - (n_face_node - (j1-1)*j1/2) + 1 ;
123 ij_lattice = point_basic<size_t>(i,j);
127geo_element_fix_triangle_indirect (
134 point_basic<size_t> ij_lattice;
135 geo_element_loc_tri_inod2lattice (loc_itri_j, order, ij_lattice);
136 size_t i = ij_lattice[0];
137 size_t j = ij_lattice[1];
142 coord [2] = order - i - j;
143 size_t i_fix = coord [shift%3];
144 size_t j_fix = coord [(shift+1)%3];
145 if (orient < 0) std::swap (i_fix, j_fix);
147 size_t loc_tri_nnod = (order-1)*(order-2)/2;
148 size_t j1_fix = order - j_fix;
149 size_t loc_itri_j_fix = (loc_tri_nnod - (j1_fix-1)*j1_fix/2) + (i_fix-1);
150 return loc_itri_j_fix;
154geo_element_fix_triangle_indirect (
155 const geo_element& K,
161 if (K.dimension() < 3 ||
162 (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
165 return geo_element_fix_triangle_indirect (
166 K.face_indirect(loc_ifac).orientation(),
167 K.face_indirect(loc_ifac).shift(),
176geo_element_loc_qua_inod2lattice (
179 point_basic<size_t>& ij_lattice)
183 ij_lattice[0] = (loc_qua_inod % (order-1)) + 1;
184 ij_lattice[1] = (loc_qua_inod / (order-1)) + 1;
188geo_element_fix_quadrangle_indirect (
195 point_basic<size_t> ij_lattice;
196 geo_element_loc_qua_inod2lattice (loc_iqua_j, order, ij_lattice);
197 size_t i = ij_lattice[0];
198 size_t j = ij_lattice[1];
203 coord [2] = order - i;
204 coord [3] = order - j;
205 size_t i_fix = coord [shift%4];
206 size_t j_fix = coord [(shift+1)%4];
207 if (orient < 0) std::swap (i_fix, j_fix);
209 size_t loc_iqua_j_fix = (order-1)*(j_fix-1) + (i_fix-1);
210 return loc_iqua_j_fix;
214geo_element_fix_quadrangle_indirect (
215 const geo_element& K,
221 if (K.dimension() < 3 || order <= 2 ||
222 (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
225 return geo_element_fix_quadrangle_indirect (
226 K.face_indirect(loc_ifac).orientation(),
227 K.face_indirect(loc_ifac).shift(),
short int orientation_type
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Float delta(Float f, Float g)
This file is part of Rheolef.