Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
23// shared code, used in:
24// - geo_element_v4.cc
25// - msh2geo.cc
26// in order to avoid code duplication
27//
28namespace 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
35template <typename GeoElement, typename OrientationType, typename ShiftType>
36static
37void
38geo_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
57template <typename GeoElement, typename OrientationType, typename ShiftType>
58static
59void
60geo_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// --------------
86static
87size_t
88geo_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// --------------
100static
101void
102geo_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}
125static
126size_t
127geo_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}
152static
153size_t
154geo_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// --------------
174static
175void
176geo_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}
186static
187size_t
188geo_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}
212static
213size_t
214geo_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
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.