Rheolef  7.1
an efficient C++ finite element environment
geo_element_aux.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
2 #define _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
3 // shared code, used in:
24 // - geo_element_v4.cc
25 // - msh2geo.cc
26 // in order to avoid code duplication
27 //
28 namespace rheolef {
29 
30 // =============================================================================
31 // compute orientation & shift
32 // =============================================================================
33 // for a triangular side (dis_iv0,dis_iv1,dis_iv2)
34 // assume that current geo_element is a triangle that contains the same vertices
35 template <typename GeoElement, typename OrientationType, typename ShiftType>
36 static
37 void
38 geo_element_get_orientation_and_shift (
39  const GeoElement& S,
40  size_t dis_iv0, size_t dis_iv1, size_t dis_iv2,
41  OrientationType& orient,
42  ShiftType& shift)
43 {
44  if (S[0] == dis_iv0) {
45  orient = (S[1] == dis_iv1) ? 1 : -1;
46  shift = 0;
47  } else if (S[0] == dis_iv1) {
48  orient = (S[1] == dis_iv2) ? 1 : -1;
49  shift = 1;
50  } else {
51  orient = (S[1] == dis_iv0) ? 1 : -1;
52  shift = 2;
53  }
54 }
55 // for a quadrangular side (dis_iv0,dis_iv1,dis_iv2,dis_iv3)
56 // assume that current geo_element is a quadrangular that contains the same vertices
57 template <typename GeoElement, typename OrientationType, typename ShiftType>
58 static
59 void
60 geo_element_get_orientation_and_shift (
61  const GeoElement& S,
62  size_t dis_iv0, size_t dis_iv1, size_t dis_iv2, size_t dis_iv3,
63  OrientationType& orient,
64  ShiftType& shift)
65 {
66  if (S[0] == dis_iv0) {
67  orient = (S[1] == dis_iv1) ? 1 : -1;
68  shift = 0;
69  } else if (S[0] == dis_iv1) {
70  orient = (S[1] == dis_iv2) ? 1 : -1;
71  shift = 1;
72  } else if (S[0] == dis_iv2) {
73  orient = (S[1] == dis_iv3) ? 1 : -1;
74  shift = 2;
75  } else {
76  orient = (S[1] == dis_iv0) ? 1 : -1;
77  shift = 3;
78  }
79 }
80 // =============================================================================
81 // fix rotation and orientation on a 2d edge or 3d face
82 // =============================================================================
83 // --------------
84 // 1) edge
85 // --------------
86 static
87 size_t
88 geo_element_fix_edge_indirect (
89  const geo_element& K,
90  size_t loc_iedg,
91  size_t loc_iedg_j,
92  size_t order)
93 {
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;
96 }
97 // --------------
98 // 2) triangles
99 // --------------
100 static
101 void
102 geo_element_loc_tri_inod2lattice (
103  size_t loc_tri_inod,
104  size_t order,
105  point_basic<size_t>& ij_lattice)
106 {
107  // search (i,j) integers such that
108  // loc_tri_inod = (n_face_node - (j1-1)*j1/2) + i-1, where j1=order-1-j
109  // first, search j1 integer such that for i=1:
110  // i - 1 = 0 = loc_tri_inod - (n_face_node - (j1-1)*j1/2);
111  // <=> j1^2 - j1 - 2*(n_face_node - loc_tri_inod) = 0
112  size_t n_face_node = (order-1)*(order-2)/2;
113  double a = 1;
114  double b = -1;
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);
119  // j1_minus = (-b - sqrt(delta))/(2*a); is negative always
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);
124 }
125 static
126 size_t
127 geo_element_fix_triangle_indirect (
130  size_t order,
131  size_t loc_itri_j)
132 {
133  // 1) compute lattice coord (i,j) on this face, from current loc_itri_j
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];
138  // 2) then shift and/or inverse-orient (i,j)
139  size_t coord[3];
140  coord [0] = i;
141  coord [1] = j;
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);
146  // 3) re-compute the fixed loc_itri_j face index
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;
151 }
152 static
153 size_t
154 geo_element_fix_triangle_indirect (
155  const geo_element& K,
156  size_t loc_ifac,
157  size_t loc_itri_j,
158  size_t order)
159 {
160  // shift and/or inverse-orient the lattice(i,j)
161  if (K.dimension() < 3 ||
162  (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
163  return loc_itri_j;
164  }
165  return geo_element_fix_triangle_indirect (
166  K.face_indirect(loc_ifac).orientation(),
167  K.face_indirect(loc_ifac).shift(),
168  order,
169  loc_itri_j);
170 }
171 // --------------
172 // 3) quadrangles
173 // --------------
174 static
175 void
176 geo_element_loc_qua_inod2lattice (
177  size_t loc_qua_inod,
178  size_t order,
179  point_basic<size_t>& ij_lattice)
180 {
181  // search (i,j) integers in [0:order-1[ such that
182  // loc_qua_inod = (order-1)*(j-1) + (i-1);
183  ij_lattice[0] = (loc_qua_inod % (order-1)) + 1;
184  ij_lattice[1] = (loc_qua_inod / (order-1)) + 1;
185 }
186 static
187 size_t
188 geo_element_fix_quadrangle_indirect (
191  size_t order,
192  size_t loc_iqua_j)
193 {
194  // 1) compute lattice coord (i,j) on this face, from current loc_itri_j
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];
199  // 2) then shift and/or inverse-orient (i,j)
200  size_t coord [4];
201  coord [0] = i;
202  coord [1] = j;
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);
208  // 3) re-compute the fixed loc_iqua_j face index
209  size_t loc_iqua_j_fix = (order-1)*(j_fix-1) + (i_fix-1);
210  return loc_iqua_j_fix;
211 }
212 static
213 size_t
214 geo_element_fix_quadrangle_indirect (
215  const geo_element& K,
216  size_t loc_ifac,
217  size_t loc_iqua_j,
218  size_t order)
219 {
220  // shift and/or inverse-orient the lattice(i,j)
221  if (K.dimension() < 3 || order <= 2 ||
222  (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
223  return loc_iqua_j;
224  }
225  return geo_element_fix_quadrangle_indirect (
226  K.face_indirect(loc_ifac).orientation(),
227  K.face_indirect(loc_ifac).shift(),
228  order,
229  loc_iqua_j);
230 }
231 
232 } // namespace rheolef
233 #endif // _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
rheolef::geo_element_indirect::shift_type
short int shift_type
Definition: geo_element_indirect.h:39
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::geo_element_indirect::orientation_type
short int orientation_type
Definition: geo_element_indirect.h:38
a
Definition: diffusion_isotropic.h:25
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
delta
Float delta(Float f, Float g)
Definition: mosolov_error_yield_surface.cc:29
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153