rheolef  6.5
piola.cc
Go to the documentation of this file.
1 #include "rheolef/piola.h"
2 
3 namespace rheolef {
4 
5 // F_K : hat_K --> K
6 // hat_x --> x
7 template<class T, class M>
8 point_basic<T>
10  const geo_basic<T,M>& omega,
11  const basis_on_pointset<T>& piola_on_pointset,
12  reference_element hat_K,
13  const std::vector<size_t>& dis_inod,
14  size_t q)
15 {
16  typedef typename geo_basic<T,M>::size_type size_type;
17  size_type d = omega.dimension();
18  size_type loc_inod = 0;
19  point_basic<T> xq;
20  piola_on_pointset.evaluate (hat_K, q);
22  iter = piola_on_pointset.begin(),
23  last = piola_on_pointset.end();
24  iter != last;
25  iter++, loc_inod++) {
26  const T& cnod = *iter;
27  const point_basic<T>& xnod = omega.dis_node (dis_inod[loc_inod]);
28  for (size_type i = 0; i < d; i++) {
29  xq[i] += cnod*xnod[i];
30  }
31  }
32  return xq;
33 }
34 template<class T, class M>
35 void
37  const geo_basic<T,M>& omega,
38  const basis_on_pointset<T>& piola_on_pointset,
39  reference_element hat_K,
40  const std::vector<size_t>& dis_inod,
41  const point_basic<T>& hat_x,
42  tensor_basic<T>& DF)
43 {
44  typedef typename geo_basic<T,M>::size_type size_type;
45  size_type d = omega.dimension();
46  size_type map_d = hat_K.dimension();
47  size_type loc_inod = 0;
48  DF.fill (0);
49  piola_on_pointset.evaluate_grad (hat_K, hat_x);
51  first_grad_tr = piola_on_pointset.begin_grad(),
52  last_grad_tr = piola_on_pointset.end_grad();
53  first_grad_tr != last_grad_tr;
54  first_grad_tr++, loc_inod++) {
55  const point_basic<T>& xnod = omega.dis_node (dis_inod[loc_inod]);
56  cumul_otimes (DF, xnod, *first_grad_tr, d, map_d);
57  }
58 }
59 template<class T, class M>
60 void
62  const geo_basic<T,M>& omega,
63  const basis_on_pointset<T>& piola_on_pointset,
64  reference_element hat_K,
65  const std::vector<size_t>& dis_inod,
66  size_t q,
67  tensor_basic<T>& DF)
68 {
69  typedef typename geo_basic<T,M>::size_type size_type;
70  size_type d = omega.dimension();
71  size_type map_d = hat_K.dimension();
72  size_type loc_inod = 0;
73  DF.fill (0);
74  piola_on_pointset.evaluate_grad (hat_K, q);
76  first_grad_tr = piola_on_pointset.begin_grad(),
77  last_grad_tr = piola_on_pointset.end_grad();
78  first_grad_tr != last_grad_tr;
79  first_grad_tr++, loc_inod++) {
80  const point_basic<T>& xnod = omega.dis_node (dis_inod[loc_inod]);
81  cumul_otimes (DF, xnod, *first_grad_tr, d, map_d);
82  }
83 }
84 template <class T>
85 T
86 det_jacobian_piola_transformation (const tensor_basic<T>& DF, size_t d , size_t map_d)
87 {
88  if (d == map_d) {
89  return DF.determinant (map_d);
90  }
91  /* surface jacobian: references:
92  * Spectral/hp element methods for CFD
93  * G. E. M. Karniadakis and S. J. Sherwin
94  * Oxford university press
95  * 1999
96  * page 165
97  */
98  switch (map_d) {
99  case 0: return 1;
100  case 1: return norm(DF.col(0));
101  case 2: return norm(vect(DF.col(0), DF.col(1)));
102  default:
103  error_macro ("det_jacobian_piola_transformation: unsupported element dimension "
104  << map_d << " in " << d << "D mesh.");
105  return 0;
106  }
107 }
108 template<class T, class M>
109 point_basic<T>
111  const geo_basic<T,M>& omega,
112  const geo_element& S,
113  const tensor_basic<T>& DF,
114  size_t d)
115 {
116  switch (d) {
117  case 1: { // point in 1D: DF[1][0] is empty, so scan S[0] node connectivity
118  typedef typename geo_basic<T,M>::size_type size_type;
119  if (S.dimension() + 1 == omega.map_dimension()) {
120  size_type dis_ie = S.master(0);
121  check_macro (dis_ie != std::numeric_limits<size_type>::max(), "normal: requires neighbours initialization");
122  const geo_element& K = omega.dis_get_geo_element (S.dimension()+1, dis_ie);
123  Float sign = (S[0] == K[1]) ? 1 : -1;
124  return point_basic<T>(sign);
125  }
126  if (omega.variant() == geo_abstract_base_rep<T>::geo_domain_indirect) {
127  size_type dis_isid = S.dis_ie();
128  size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
129  size_type isid = dis_isid - first_dis_isid;
130  check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
131  const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
132  const geo_basic<T,M>& bgd_dom = omega.get_background_domain();
133  const geo_element& bgd_S = bgd_omega.get_geo_element(0, isid);
134  size_type bgd_dis_ie = bgd_S.master(0);
136  "normal: bgd_S={"<<bgd_S<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
137  const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
138  Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
139  return point_basic<T>(sign);
140  } else {
141  size_type dis_isid = S.dis_ie();
142  size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
143  size_type isid = dis_isid - first_dis_isid;
144  check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
145  const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
146  const geo_basic<T,M>& bgd_dom = omega.get_background_domain();
147  const geo_element& bgd_S = bgd_dom[isid]; // TODO: pas clair, differe du cas precedent ?
148  size_type bgd_dis_ie = bgd_S.master(0);
150  "normal: bgd_S={"<<bgd_S<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
151  const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
152  Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
153  return point_basic<T>(sign);
154  }
155  }
156  case 2: { // edge in 2D
157  // 2D: S=edge(a,b) then t=b-a, DF=[t] and n = (t1,-t0) => det(n,1)=1 : (n,t) is direct
158  point_basic<T> t = DF.col(0);
159  t /= norm(t);
160  return point_basic<T>(t[1], -t[0]);
161  }
162  case 3: { // 3D: S=triangle(a,b,c) then t0=b-a, t1=c-a, DF=[t0,t1] and n = t0^t1/|t0^t1|.
163  point_basic<T> t0 = DF.col(0);
164  point_basic<T> t1 = DF.col(1);
165  point_basic<T> n = vect(t0,t1);
166  n /= norm(n);
167  return n;
168  }
169  default:
170  error_macro ("normal: unsupported " << d << "D mesh.");
171  return point_basic<T>();
172  }
173 }
174 // if K=triangle(a,b,c) then u=ab=b-a, v=ac=c-a and w = n = u^v/|u^v|.
175 // and inv(DF)^T = [ac^n/|ab^ac|, -ab^n/|ab^ac|, n]
176 template<class T>
177 tensor_basic<T>
179  const tensor_basic<T>& DF,
180  size_t d,
181  size_t map_d)
182 {
183  if (d == map_d) {
184  return inv(DF, map_d);
185  }
186  tensor_basic<T> invDF;
187  switch (map_d) {
188  case 0: { // point in 1D
189  invDF(0,0) = 1;
190  return invDF;
191  }
192  case 1: { // segment in 2D
193  point_basic<T> t = DF.col(0);
194  invDF.set_row (t/norm2(t), 0, d);
195  return invDF;
196  }
197  case 2: {
198  point_basic<T> t0 = DF.col(0);
199  point_basic<T> t1 = DF.col(1);
200  point_basic<T> n = vect(t0,t1);
201  T det2 = norm2(n);
202  point_basic<T> v0 = vect(t1,n)/det2;
203  point_basic<T> v1 = - vect(t0,n)/det2;
204  invDF.set_row (v0, 0, d);
205  invDF.set_row (v1, 1, d);
206  return invDF;
207  }
208  default:
209  error_macro ("pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
210  << map_d << " in " << d << "D mesh.");
211  return invDF;
212  }
213 }
214 // point_basic<T> xq = rheolef::piola_transformation (_omega, _piola_table, K, dis_inod, q);
215 template<class T>
216 T
218 {
219  switch (sys_coord) {
220  case space_constant::axisymmetric_rz: return xq[0];
221  case space_constant::axisymmetric_zr: return xq[1];
222  case space_constant::cartesian: return 1;
223  default: {
224  fatal_macro ("unsupported coordinate system `"
225  << space_constant::coordinate_system_name (sys_coord) << "'");
226  return 0;
227  }
228  }
229 }
230 template<class T, class M>
231 T
233  const geo_basic<T,M>& omega,
234  const basis_on_pointset<T>& piola_table,
235  const geo_element& K,
236  const std::vector<size_t>& dis_inod,
237  size_t q)
238 {
239  if (omega.coordinate_system() == space_constant::cartesian) return 1;
240  point_basic<T> xq = rheolef::piola_transformation (omega, piola_table, K, dis_inod, q);
241  return weight_coordinate_system (omega.coordinate_system(), xq);
242 }
243 template<class T>
244 void
246  const tensor_basic<T>& DF,
247  size_t d,
248  size_t map_d,
249  tensor_basic<T>& P)
250 {
252  switch (map_d) {
253  case 1: {
254  point_basic<Float> t = DF.col(0);
255  check_macro (d == map_d+1, "unexpected dimension map_d="<<map_d<<" and d="<<d);
256  n = point_basic<T>(t[1],-t[0]);
257  break;
258  }
259  case 2: {
260  point_basic<Float> t0 = DF.col(0);
261  point_basic<Float> t1 = DF.col(1);
262  n = vect(t0,t1);
263  break;
264  }
265  default:
266  error_macro ("unexpected dimension "<<map_d);
267  }
268  n = n/norm(n);
269  for (size_t l = 0; l < d; l++) {
270  for (size_t m = 0; m < d; m++) {
271  P(l,m) = - n[l]*n[m];
272  }
273  P(l,l) += 1;
274  }
275 }
276 // F_K^{-1} : K --> hat(K)
277 // x --> hat(x)
278 template<class T>
279 static
280 inline
281 point_basic<T>
283  const point_basic<T>& x,
284  const point_basic<T>& a,
285  const point_basic<T>& b)
286 {
287  return point_basic<T>((x[0]-a[0])/(b[0]-a[0]));
288 }
289 template<class T>
290 static
291 inline
292 point_basic<T>
294  const point_basic<T>& x,
295  const point_basic<T>& a,
296  const point_basic<T>& b,
297  const point_basic<T>& c)
298 {
299  T t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
300  T t11 = -a[0]+x[0];
301  T t15 = -a[1]+x[1];
302  return point_basic<T>((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
303  ( b[1]-a[1])*t9*t11-( b[0]-a[0])*t9*t15);
304 }
305 template<class T>
306 static
307 inline
308 point_basic<T>
310  const point_basic<T>& x,
311  const point_basic<T>& a,
312  const point_basic<T>& b,
313  const point_basic<T>& c,
314  const point_basic<T>& d)
315 {
316  tensor_basic<T> A;
317  point_basic<T> ax;
318  for (size_t i = 0; i < 3; i++) {
319  ax[i] = x[i]-a[i];
320  A(i,0) = b[i]-a[i];
321  A(i,1) = c[i]-a[i];
322  A(i,2) = d[i]-a[i];
323  }
324  tensor_basic<T> inv_A;
325  bool is_singular = ! invert_3x3 (A, inv_A);
326  check_macro(!is_singular, "inv_piola: singular transformation in a tetrahedron");
327  point_basic<T> hat_x = inv_A*ax;
328  return hat_x;
329 }
330 template<class T, class M>
331 point_basic<T>
333  const geo_basic<T,M>& omega,
334  reference_element hat_K,
335  const std::vector<size_t>& dis_inod,
336  const point_basic<T>& x)
337 {
338  check_macro (omega.order() == 1, "inverse piola: mesh order > 1: not yet");
339  switch (hat_K.variant()) {
340  case reference_element::e: return inv_piola_e (x, omega.dis_node(dis_inod [0]),
341  omega.dis_node(dis_inod [1]));
342  case reference_element::t: return inv_piola_t (x, omega.dis_node(dis_inod [0]),
343  omega.dis_node(dis_inod [1]),
344  omega.dis_node(dis_inod [2]));
345  case reference_element::T: return inv_piola_T (x, omega.dis_node(dis_inod [0]),
346  omega.dis_node(dis_inod [1]),
347  omega.dis_node(dis_inod [2]),
348  omega.dis_node(dis_inod [3]));
349  default: {
350  fatal_macro ("inverse piola: non-linear transformation on '"<<hat_K.name()<<"' element: not yet");
351  return point_basic<T>();
352  }
353  }
354 }
355 #define _RHEOLEF_instanciation1(T) \
356 template \
357 T \
358 det_jacobian_piola_transformation ( \
359  const tensor_basic<T>& DF, \
360  size_t d, \
361  size_t map_d); \
362 template \
363 tensor_basic<T> \
364 pseudo_inverse_jacobian_piola_transformation ( \
365  const tensor_basic<T>& DF, \
366  size_t d, \
367  size_t map_d); \
368 template \
369 T \
370 weight_coordinate_system ( \
371  space_constant::coordinate_type sys_coord, \
372  const point_basic<T>& xq); \
373 template \
374 void \
375 map_projector ( \
376  const tensor_basic<T>& DF, \
377  size_t d, \
378  size_t map_d, \
379  tensor_basic<T>& P); \
380 
381 
382 #define _RHEOLEF_instanciation2(T,M) \
383 template \
384 point_basic<T> \
385 piola_transformation ( \
386  const geo_basic<T,M>& omega, \
387  const basis_on_pointset<T>& piola_on_pointset, \
388  reference_element hat_K, \
389  const std::vector<size_t>& dis_inod, \
390  size_t q); \
391 template \
392 point_basic<T> \
393 inverse_piola_transformation ( \
394  const geo_basic<T,M>& omega, \
395  reference_element hat_K, \
396  const std::vector<size_t>& dis_inod, \
397  const point_basic<T>& x); \
398 template \
399 void \
400 jacobian_piola_transformation ( \
401  const geo_basic<T,M>& omega, \
402  const basis_on_pointset<T>& piola_on_pointset, \
403  reference_element hat_K, \
404  const std::vector<size_t>& dis_inod, \
405  const point_basic<T>& hat_x, \
406  tensor_basic<T>& DF); \
407 template \
408 void \
409 jacobian_piola_transformation ( \
410  const geo_basic<T,M>& omega, \
411  const basis_on_pointset<T>& piola_on_pointset, \
412  reference_element hat_K, \
413  const std::vector<size_t>& dis_inod, \
414  size_t q, \
415  tensor_basic<T>& DF); \
416 template \
417 point_basic<T> \
418 normal_from_piola_transformation ( \
419  const geo_basic<T,M>& omega, \
420  const geo_element& S, \
421  const tensor_basic<T>& DF, \
422  size_t d); \
423 template \
424 T \
425 weight_coordinate_system ( \
426  const geo_basic<T,M>& omega, \
427  const basis_on_pointset<T>& piola_table, \
428  const geo_element& K, \
429  const std::vector<size_t>& dis_inod, \
430  size_t q); \
431 
432 
433 
435 _RHEOLEF_instanciation2(Float,sequential)
436 #ifdef _RHEOLEF_HAVE_MPI
438 #endif // _RHEOLEF_HAVE_MPI
439 
440 } // namespace rheolef
441