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