Rheolef  7.2
an efficient C++ finite element environment
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 
49 namespace rheolef {
50 
51 // ====================================================================================
52 // quadrangle (a,b,c,d) with curved edge (a,b)
53 // ====================================================================================
54 template<class T>
56 public:
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:
69 protected:
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  }
102 protected:
103 // data:
106 };
107 // ====================================================================================
108 // triangle (a,b,c) with i-th curved edge
109 // ====================================================================================
110 template<class T>
112 public:
113 // allocator:
114  curved_ball_t (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0,
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  }
149 protected:
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 // ====================================================================================
173 template<class T>
175 public:
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:
284 protected:
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  }
329  point_basic<T> f_abc (const point_basic<T>& x) const { return project_on_boundary (x); }
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 // ====================================================================================
350 template<class T>
352 public:
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:
364 protected:
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 
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 > f_bc(const point_basic< T > &x) const
point_basic< T > f_bd(const point_basic< T > &x) const
point_basic< T > f_ad(const point_basic< T > &x) const
std::array< bool, 4 > is_curved_fac
point_basic< T > x_edge(size_t loc_iedg, const T &hat_x) const
std::array< bool, 4 > is_bdry_fac
point_basic< T > f_ab(const point_basic< T > &x) const
std::array< bool, 6 > is_bdry_edg
void set_boundary_face(size_t loc_ifac_curved)
point_basic< T > project_on_boundary(const point_basic< T > &x) const
point_basic< T > f_cd(const point_basic< T > &x) const
point_basic< T > x_face(size_t loc_ifac, const point_basic< T > &hat_x) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
void set_boundary_edge(size_t loc_iedg_curved)
point_basic< T > f_abc(const point_basic< T > &x) const
std::array< point_basic< T >, 4 > node
point_basic< T > f_acd(const point_basic< T > &x) const
point_basic< T > f_bcd(const point_basic< T > &x) const
point_basic< T > f_ac(const point_basic< T > &x) const
point_basic< T > f_abd(const point_basic< T > &x) const
point_basic< T > f_small_carre(const T &x, const T &y) 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_bc(const T &s) const
point_basic< T > f_ab(const T &s) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
point_basic< T > f_ad(const T &s) const
point_basic< T > f_dc(const T &s) 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 > project_on_boundary(const point_basic< T > &x) const
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)
point_basic< T >
Definition: piola_fem.h:135
#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