Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_element_curved_ball.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
2#define _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
23//
24// transfinite piola transformaton from
25// * a triangle (a,b,c) or a quadrangle(a,b,c,d)
26// to a curved triangle or quadrangle
27// where edge (a,b) is an arc of a circle.
28// * a tetra (a,b,c,d) or an hexa(a,b,c,d,e,f,g,h)
29// to a curved tetra or hexa
30// where face (a,b,c), or (a,b,c,d) for hexa,
31// is a part of sphere boundary
32//
33// Remark: can be generalized from a circle/ball
34// to any curved boundary, by sending a function
35// that project or replace internal nodes onto the boundary.
36// See SherKar-2005, page 224.
37//
38// author: Pierre.Saramito@imag.fr
39//
40// date: 22 november 2011
41//
42// TODO:
43// * tetra, hexa
44// * ellipse(r1,r2,r3,c)
45// * parametrized by a function that project on the bdry ?
46//
47#include "point.h"
48
49namespace rheolef {
50
51// ====================================================================================
52// quadrangle (a,b,c,d) with curved edge (a,b)
53// ====================================================================================
54template<class T>
56public:
57// allocator:
58 curved_ball_q (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
59 const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
60 : a(a0), b(b0), c(c0), d(d0), center(center0), radius(radius0) {}
61// accessor:
63 // transform large square [-1:1]^2 into small one [0:1]^2
64 T x = 0.5*(hat_x[0]+1);
65 T y = 0.5*(hat_x[1]+1);
66 return f_small_carre (x, y);
67 }
68// internals:
69protected:
70 point_basic<T> f_ab (const T& s) const {
71 point_basic<T> x = (1-s)*a + s*b;
72 return center + radius*(x-center)/norm(x-center);
73 }
74 point_basic<T> f_bc (const T& s) const { return (1-s)*b + s*c; }
75 point_basic<T> f_dc (const T& s) const { return (1-s)*d + s*c; }
76 point_basic<T> f_ad (const T& s) const { return (1-s)*a + s*d; }
77 /*
78 Transformation for a square: from eqn (17) in:
79 @article{GorHal-1972,
80 author = {W. J. Gordon and C. A. Hall},
81 title = {Transfinite element methods:
82 blending-function interpolation over
83 arbitrary curved element domains},
84 journal = {Numer. Math.},
85 volume = 21,
86 pages = {109--129},
87 year = 1973,
88 url = {http://www.springerlink.com/content/p5517206551228j6/fulltext.pdf},
89 keywords = {method!fem!curved}
90 }
91 */
92 point_basic<T> f_small_carre (const T& x, const T& y) const {
93 // small carre: 0 <= x,y <= 1
94 return (1-x)*f_ad(y) + x*f_bc(y)
95 + (1-y)*f_ab(x) + y*f_dc(x)
96 - (1-x)*(1-y)*a
97 - x *(1-y)*b
98 - x * y *c
99 - (1-x)* y *d
100 ;
101 }
102protected:
103// data:
106};
107// ====================================================================================
108// triangle (a,b,c) with i-th curved edge
109// ====================================================================================
110template<class T>
112public:
113// allocator:
115 size_t loc_curved_iedg, const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
116 : node(), center(center0), radius(radius0), is_bdry_edg()
117 {
118 node[0] = a0;
119 node[1] = b0;
120 node[2] = c0;
121 is_bdry_edg.fill (false);
122 is_bdry_edg [loc_curved_iedg] = true;
123 }
124// accessor:
126 /* Transform for a triangle, from eqn (27) in:
127 @article{DeySheFla-1997,
128 title = {Geometry representation issues associated with p-version finite element computations},
129 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
130 journal = {Comput. Meth. Appl. Mech. Engrg.},
131 volume = 150,
132 year = 1997,
133 pages = {39--55}
134 }
135 Remark : SheKar-1999, p. 162-163 : cite GorHal-1973 for a quadrangle
136 but says that the triangle-to-square transform + the quadrangle formulae
137 does not work here (its true, I have checked) because of a singular jacobian.
138 Triangle requires thus a specific formulae as it is done here.
139 Also: similar ref: eqn (3.17), page 169 in SolSegDol-2004 is incorrect.
140 */
141 T la = 1 - hat_x[0] - hat_x[1]; // barycentric coords
142 T lb = hat_x[0];
143 T lc = hat_x[1];
144 return 0.5*( la/(1-lb)*edge(0,lb) + lb/(1-la)*edge(0,1-la)
145 + lb/(1-lc)*edge(1,lc) + lc/(1-lb)*edge(1,1-lb)
146 + lc/(1-la)*edge(2,la) + la/(1-lc)*edge(2,1-lc)
147 - la*node[0] - lb*node[1] - lc*node[2]);
148 }
149protected:
150// internals:
152 return center + radius*(x-center)/norm(x-center);
153 }
154 point_basic<T> edge (size_t loc_iedg, const T& hat_x) const {
155 const point_basic<T>& a = node [loc_iedg];
156 const point_basic<T>& b = node [(loc_iedg+1)%3];
157 point_basic<T> x = (1-hat_x)*a + hat_x*b;
158 if (is_bdry_edg [loc_iedg]) {
159 return project_on_boundary (x);
160 } else {
161 return x;
162 }
163 }
164// data:
165 std::array<point_basic<T>,3> node;
168 std::array<bool,3> is_bdry_edg;
169};
170// ====================================================================================
171// tetrahedron (a,b,c,d) with i-th curved face
172// ====================================================================================
173template<class T>
175public:
176// allocator:
177 curved_ball_T (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
178 const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
179 : node(), center(center0), radius(radius0), is_bdry_edg(), is_bdry_fac(), is_curved_fac()
180 {
181 node[0] = a0;
182 node[1] = b0;
183 node[2] = c0;
184 node[3] = d0;
185 is_bdry_fac.fill(false);
186 is_bdry_edg.fill(false);
187 is_curved_fac.fill(false);
188 }
189// modifiers:
190 void set_boundary_face (size_t loc_ifac_curved)
191 {
192 is_bdry_fac[loc_ifac_curved] = true;
193 // then, all faces are either on boundary or have a curved edge:
194 for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
195 if (!is_bdry_fac [loc_ifac]) {
196 is_curved_fac[loc_ifac] = true;
197 }
198 }
199 // the three edges of the bdry face are also on the boundary
200 for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
201 size_t loc_jedg = reference_element_T::face2edge (loc_ifac_curved, loc_ifac_jedg);
202 is_bdry_edg [loc_jedg] = true;
203 }
204 }
205 void set_boundary_edge (size_t loc_iedg_curved)
206 {
207 is_bdry_edg [loc_iedg_curved] = true;
208 // also, the two faces that contains this edge are then curved (if not yet on the boundary)
209 for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
210 for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
211 size_t loc_jedg = reference_element_T::face2edge (loc_ifac, loc_ifac_jedg);
212 if (loc_jedg != loc_iedg_curved) continue;
213 // this face contains the boundary edge:
214 if (!is_bdry_fac [loc_ifac]) {
215 is_curved_fac[loc_ifac] = true;
216 break;
217 }
218 }
219 }
220 }
221// accessor:
223 /* Transform for a tetra, from eqn (30) in:
224 @article{DeySheFla-1997,
225 title = {Geometry representation issues associated with p-version finite element computations},
226 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
227 journal = {Comput. Meth. Appl. Mech. Engrg.},
228 volume = 150,
229 pages = {39--55},
230 year = 1997
231 }
232 */
233 // 1) compute the baycentric coordinates: lambda_i(xj) = delta_ij
234 std::array<T,4> l;
235 l[0] = 1 - hat_x[0] - hat_x[1] - hat_x[2];
236 l[1] = hat_x[0];
237 l[2] = hat_x[1];
238 l[3] = hat_x[2];
239
240 // 2) compute the contribution of all faces
241 std::array<size_t,3> S; // table of local vertices indexes
242 std::array<bool,4> is_on_S; // when vertex is on S
243 point_basic<T> x_all_faces (0,0,0);
244 for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
245 is_on_S.fill (false);
246 for (size_t loc_ifac_jv = 0; loc_ifac_jv < 3; loc_ifac_jv++) {
247 S[loc_ifac_jv] = reference_element_T::subgeo_local_node (1, 2, loc_ifac, loc_ifac_jv);
248 is_on_S [S[loc_ifac_jv]] = true;
249 }
250 // jv_3 is the only index in 0..3 that is not on the face ifac:
251 size_t jv_3 = size_t(-1);
252 for (size_t loc_jv = 0; loc_jv < 4; loc_jv++) {
253 if (!is_on_S [loc_jv]) jv_3 = loc_jv;
254 }
255 T sum = l[S[0]] + l[S[1]] + l[S[2]];
256 point_basic<T> hat_xi (l[S[1]]/sum, l[S[2]]/sum);
257 x_all_faces += (1 - l[jv_3])*x_face (loc_ifac, hat_xi);
258 }
259 // 2) compute the contribution of all edges
260 point_basic<T> x_all_edges (0,0,0);
261 for (size_t loc_iedg = 0; loc_iedg < 6; loc_iedg++) {
262 // edge = [a,b] and c & d are the two others verices of the tetra
263 size_t ia = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 0);
264 size_t ib = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 1);
265 size_t ic = size_t(-1);
266 size_t id = size_t(-1);
267 for (size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
268 if (loc_iv == ia || loc_iv == ib) continue;
269 if (ic == size_t(-1)) { ic = loc_iv; continue; }
270 id = loc_iv;
271 }
272 assert_macro (ic != size_t(-1) && id != size_t(-1), "unexpected");
273 T hat_xi = l[ib]/(l[ia] + l[ib]);
274 x_all_edges += (1 - l[ic] - l[id])*x_edge (loc_iedg, hat_xi);
275 }
276 // 3) compute the contribution of all vertices
277 point_basic<T> x_all_vertices (0,0,0);
278 for (size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
279 x_all_vertices += l[loc_iv]*node[loc_iv];
280 }
281 return x_all_faces - x_all_edges + x_all_vertices;
282 }
283// internals:
284protected:
285 point_basic<T> x_face (size_t loc_ifac, const point_basic<T>& hat_x) const {
286 size_t ia = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 0);
287 size_t ib = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 1);
288 size_t ic = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 2);
289 const point_basic<T>& a = node[ia];
290 const point_basic<T>& b = node[ib];
291 const point_basic<T>& c = node[ic];
292 if (is_curved_fac [loc_ifac]) {
293 // search the edge who lives on the boundary
294 size_t loc_ifac_jedg_bdry = size_t(-1);
295 for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
296 size_t loc_jedg = reference_element_T::face2edge (loc_ifac, loc_ifac_jedg);
297 if (is_bdry_edg [loc_jedg]) {
298 loc_ifac_jedg_bdry = loc_ifac_jedg;
299 break;
300 }
301 }
302 assert_macro (loc_ifac_jedg_bdry != size_t(-1), "unexpected");
303 curved_ball_t<T> F (a, b, c, loc_ifac_jedg_bdry, center, radius);
304 return F (hat_x);
305 } else {
306 point_basic<T> x = (1 - hat_x[0] - hat_x[1])*a + hat_x[0]*b + hat_x[1]*c;
307 if (is_bdry_fac [loc_ifac]) {
308 return project_on_boundary (x);
309 } else { // strait face
310 return x;
311 }
312 }
313 }
314 point_basic<T> x_edge (size_t loc_iedg, const T& hat_x) const {
315 size_t ia = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 0);
316 size_t ib = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 1);
317 const point_basic<T>& a = node[ia];
318 const point_basic<T>& b = node[ib];
319 point_basic<T> x = (1 - hat_x)*a + hat_x*b;
320 if (is_bdry_edg [loc_iedg]) {
321 return project_on_boundary (x);
322 } else {
323 return x;
324 }
325 }
327 return center + radius*(x-center)/norm(x-center);
328 }
330 point_basic<T> f_abd (const point_basic<T>& x) const { return x; }
331 point_basic<T> f_acd (const point_basic<T>& x) const { return x; }
332 point_basic<T> f_bcd (const point_basic<T>& x) const { return x; }
333 point_basic<T> f_ab (const point_basic<T>& x) const { return project_on_boundary (x); }
334 point_basic<T> f_ac (const point_basic<T>& x) const { return project_on_boundary (x); }
335 point_basic<T> f_bc (const point_basic<T>& x) const { return project_on_boundary (x); }
336 point_basic<T> f_ad (const point_basic<T>& x) const { return x; }
337 point_basic<T> f_bd (const point_basic<T>& x) const { return x; }
338 point_basic<T> f_cd (const point_basic<T>& x) const { return x; }
339// data:
340 std::array<point_basic<T>,4> node;
343 std::array<bool,6> is_bdry_edg;
344 std::array<bool,4> is_bdry_fac;
345 std::array<bool,4> is_curved_fac;
346};
347// ====================================================================================
348// TODO: hexahedron with a curved face
349// ====================================================================================
350template<class T>
352public:
353// allocator:
354 curved_ball_H (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
355 const point_basic<T>& e0, const point_basic<T>& f0, const point_basic<T>& g0, const point_basic<T>& h0,
356 const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
357 : a(a0), b(b0), c(c0), d(d0), e(e0), f(f0), g(g0), h(h0), center(center0), radius(radius0) {}
358// accessor:
360 error_macro ("curved H: not yet");
361 return point_basic<T>();
362 }
363// internals:
364protected:
365 /* Transform for an hexa, from annex A in:
366 @article{KirSza-1997,
367 author = {G. Kir\'alyfalvi and B. A. Szab\'o},
368 title = {Quasi-regional mapping for the $p$-version
369 of the finite element method},
370 journal = {Finite Elements in Analysis and Design},
371 volume = 27,
372 pages = {85--97},
373 year = 1997,
374 keywords = {method!fem!curved!blending}
375 }
376 */
377// data:
380};
381
382} // namespace rheolef
383#endif // _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
384
see the edge page for the full documentation
point_basic< T > operator()(const point_basic< T > &z) const
curved_ball_H(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > &e0, const point_basic< T > &f0, const point_basic< T > &g0, const point_basic< T > &h0, const point_basic< T > &center0=point_basic< T >(0, 0), const T &radius0=1)
curved_ball_T(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > &center0=point_basic< T >(0, 0), const T &radius0=1)
point_basic< T > project_on_boundary(const point_basic< T > &x) const
point_basic< T > f_bcd(const point_basic< T > &x) const
std::array< bool, 4 > is_curved_fac
point_basic< T > f_cd(const point_basic< T > &x) const
std::array< bool, 4 > is_bdry_fac
point_basic< T > f_acd(const point_basic< T > &x) const
point_basic< T > f_ab(const point_basic< T > &x) const
std::array< bool, 6 > is_bdry_edg
point_basic< T > f_ac(const point_basic< T > &x) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
point_basic< T > f_bc(const point_basic< T > &x) const
point_basic< T > f_bd(const point_basic< T > &x) const
void set_boundary_face(size_t loc_ifac_curved)
point_basic< T > f_abd(const point_basic< T > &x) const
point_basic< T > x_edge(size_t loc_iedg, const T &hat_x) const
point_basic< T > f_abc(const point_basic< T > &x) const
point_basic< T > x_face(size_t loc_ifac, const point_basic< T > &hat_x) const
void set_boundary_edge(size_t loc_iedg_curved)
point_basic< T > f_ad(const point_basic< T > &x) const
std::array< point_basic< T >, 4 > node
point_basic< T > f_small_carre(const T &x, const T &y) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
curved_ball_q(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > &center0=point_basic< T >(0, 0), const T &radius0=1)
point_basic< T > f_ab(const T &s) const
point_basic< T > f_bc(const T &s) const
point_basic< T > f_ad(const T &s) const
point_basic< T > f_dc(const T &s) const
point_basic< T > project_on_boundary(const point_basic< T > &x) const
curved_ball_t(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, size_t loc_curved_iedg, const point_basic< T > &center0=point_basic< T >(0, 0), const T &radius0=1)
std::array< bool, 3 > is_bdry_edg
point_basic< T > operator()(const point_basic< T > &hat_x) const
point_basic< T > edge(size_t loc_iedg, const T &hat_x) const
std::array< point_basic< T >, 3 > node
static size_type subgeo_local_node(size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
static size_type face2edge(size_type loc_iface, size_type loc_iface_jedg)
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
point - d-dimensional physical point or vector
Definition cavity_dg.h:29
Definition cavity_dg.h:25