Rheolef  7.2
an efficient C++ finite element environment
tensor.cc
Go to the documentation of this file.
1 //
22 // 2-tensor
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 #include "rheolef/tensor.h"
27 #include "rheolef/compiler_eigen.h"
28 
29 using namespace rheolef;
30 using namespace std;
31 namespace rheolef {
32 
33 
34 // output
35 template<class T>
36 ostream&
38 {
39  switch (d) {
40  case 0 : return out;
41  case 1 : return out << _x[0][0];
42  case 2 : return out << "[" << _x[0][0] << ", " << _x[0][1] << ";\n"
43  << " " << _x[1][0] << ", " << _x[1][1] << "]";
44  default: return out << "[" << _x[0][0] << ", " << _x[0][1] << ", " << _x[0][2] << ";\n"
45  << " " << _x[1][0] << ", " << _x[1][1] << ", " << _x[1][2] << ";\n"
46  << " " << _x[2][0] << ", " << _x[2][1] << ", " << _x[2][2] << "]";
47  }
48 }
49 template<class T>
50 istream&
51 tensor_basic<T>::get (istream& in)
52 {
53  for (size_type i = 0; i < 3; i++)
54  for (size_type j = 0; j < 3; j++)
55  _x[i][j] = 0.;
56  char c;
57  if (!in.good()) return in;
58  in >> ws >> c;
59  // no tensor specifier:
60  // 1d case: "3.14"
61  if (c != '[') {
62  in.unget();
63  return in >> _x[0][0];
64  }
65  // has a tensor specifier, as:
66  // [ 23, 17; 25, 2]
67  in >> _x[0][0] >> ws >> c;
68  if (c == ']') return in; // 1d case
69  if (c != ',') error_macro ("invalid tensor input: expect `,'");
70  in >> _x[0][1];
71  in >> ws >> c;
72  size_type d = 3;
73  if (c == ',') d = 3;
74  else if (c == ';') d = 2;
75  else error_macro ("invalid tensor input: expect `,' or `;'");
76  if (d == 3) {
77  in >> _x[0][2] >> ws >> c;
78  if (c != ';') error_macro ("invalid tensor input: expect `;'");
79  }
80  in >> _x[1][0] >> ws >> c;
81  if (c != ',') error_macro ("invalid tensor input: expect `,'");
82  in >> _x[1][1] >> ws >> c;
83  if (d == 2) {
84  if (c != ']') error_macro ("invalid tensor input: expect `]'");
85  return in;
86  }
87  // d == 3
88  if (c != ',') error_macro ("invalid tensor input: expect `,'");
89  in >> _x[1][2] >> ws >> c;
90  if (c != ';') error_macro ("invalid tensor input: expect `;'");
91  in >> _x[2][0] >> ws >> c;
92  if (c != ',') error_macro ("invalid tensor input: expect `,'");
93  in >> _x[2][1] >> ws >> c;
94  if (c != ',') error_macro ("invalid tensor input: expect `,'");
95  in >> _x[2][2] >> ws >> c;
96  if (c != ']') error_macro ("invalid tensor input: expect `]'");
97  return in;
98 }
99 template<class T>
100 bool
102 {
103  for (size_type i = 0; i < 3; i++)
104  for (size_type j = 0; j < 3; j++)
105  if (_x[i][j] != b._x[i][j]) return false;
106  return true;
107 }
108 template<class T>
111 {
113  typedef typename tensor_basic<T>::size_type size_type;
114  for (size_type i = 0; i < 3; i++)
115  for (size_type j = 0; j < 3; j++)
116  b._x[i][j] = - _x[i][j];
117  return b;
118 }
119 template<class T>
122 {
124  typedef typename tensor_basic<T>::size_type size_type;
125  for (size_type i = 0; i < 3; i++)
126  for (size_type j = 0; j < 3; j++)
127  c._x[i][j] = _x[i][j] + b._x[i][j];
128  return c;
129 }
130 template<class T>
133 {
135  typedef typename tensor_basic<T>::size_type size_type;
136  for (size_type i = 0; i < 3; i++)
137  for (size_type j = 0; j < 3; j++)
138  c._x[i][j] = _x[i][j] - b._x[i][j];
139  return c;
140 }
141 template<class T>
144 {
146  for (size_type i = 0; i < 3; i++)
147  for (size_type j = 0; j < 3; j++)
148  b._x[i][j] = _x[i][j] * k;
149  return b;
150 }
151 template<class T>
154 {
155  point_basic<T> y;
156  typedef typename tensor_basic<T>::size_type size_type;
157  for (size_type i = 0; i < 3; i++)
158  for (size_type j = 0; j < 3; j++)
159  y [i] += _x[i][j] * x[j];
160  return y;
161 }
162 template<class T>
164 {
165  point_basic<T> y;
166  typedef typename tensor_basic<T>::size_type size_type;
167  for (size_type i = 0; i < 3; i++)
168  for (size_type j = 0; j < 3; j++)
169  y [i] += a._x[j][i] * x[j];
170  return y;
171 }
172 template<class T>
175 {
176  for (size_type i = 0; i < 3; i++)
177  for (size_type j = 0; j < 3; j++)
178  _x[i][j] += b._x[i][j];
179  return *this;
180 }
181 template<class T>
184 {
185  for (size_type i = 0; i < 3; i++)
186  for (size_type j = 0; j < 3; j++)
187  _x[i][j] -= b._x[i][j];
188  return *this;
189 }
190 template<class T>
193 {
194  for (size_type i = 0; i < 3; i++)
195  for (size_type j = 0; j < 3; j++)
196  _x[i][j] *= k;
197  return *this;
198 }
199 template<class T>
202 {
203  for (size_type i = 0; i < 3; i++)
204  for (size_type j = 0; j < 3; j++)
205  _x[i][j] /= k;
206  return *this;
207 }
208 template<class T>
210 {
212  typedef typename tensor_basic<T>::size_type size_type;
213  for (size_type i = 0; i < d; i++)
214  for (size_type j = 0; j < d; j++)
215  b._x[i][j] = a._x[j][i];
216  return b;
217 }
218 template<class T>
220 {
222  typedef typename tensor_basic<T>::size_type size_type;
223  switch (d) {
224  case 0:
225  break;
226  case 1:
227  b(0,0) = 1/a(0,0);
228  break;
229  case 2: {
230  T det = a.determinant(2);
231  b(0,0) = a(1,1)/det;
232  b(1,1) = a(0,0)/det;
233  b(0,1) = - a(0,1)/det;
234  b(1,0) = - a(1,0)/det;
235  break;
236  }
237  case 3:
238  default: {
239  T det = a.determinant(3);
240  b(0,0) = (a(1,1)*a(2,2) - a(1,2)*a(2,1))/det;
241  b(0,1) = (a(0,2)*a(2,1) - a(0,1)*a(2,2))/det;
242  b(0,2) = (a(0,1)*a(1,2) - a(0,2)*a(1,1))/det;
243  b(1,0) = (a(1,2)*a(2,0) - a(1,0)*a(2,2))/det;
244  b(1,1) = (a(0,0)*a(2,2) - a(0,2)*a(2,0))/det;
245  b(1,2) = (a(0,2)*a(1,0) - a(0,0)*a(1,2))/det;
246  b(2,0) = (a(1,0)*a(2,1) - a(1,1)*a(2,0))/det;
247  b(2,1) = (a(0,1)*a(2,0) - a(0,0)*a(2,1))/det;
248  b(2,2) = (a(0,0)*a(1,1) - a(0,1)*a(1,0))/det;
249  break;
250  }
251  }
252  return b;
253 }
254 template<class T>
255 void
257  size_t di, size_t dj, size_t dk)
258 {
259  typedef typename tensor_basic<T>::size_type size_type;
260  for (size_type i = 0; i < di; i++)
261  for (size_type j = 0; j < dj; j++) {
262  T sum = 0;
263  for (size_type k = 0; k < dk; k++)
264  sum += a(i,k) * b(k,j);
265  result(i,j) = sum;
266  }
267 }
268 template<class T>
271 {
273  prod (*this,b,c,3,3,3);
274  return c;
275 }
277 template<class T>
279 {
280  T r = 0;
281  typedef typename tensor_basic<T>::size_type size_type;
282  for (size_type i = 0; i < 3; i++)
283  for (size_type j = 0; j < 3; j++)
284  r += a._x[i][j] * b._x[i][j];
285  return r;
286 }
287 template<class T>
289 {
290  switch (d) {
291  case 0 : return 1;
292  case 1 : return _x[0][0];
293  case 2 : return _x[0][0]*_x[1][1] - _x[0][1]*_x[1][0];
294  case 3 : return _x[0][0]*(_x[1][1]*_x[2][2] - _x[1][2]*_x[2][1])
295  - _x[0][1]*(_x[1][0]*_x[2][2] - _x[1][2]*_x[2][0])
296  + _x[0][2]*(_x[1][0]*_x[2][1] - _x[1][1]*_x[2][0]);
297  default: {
298  error_macro ("determinant: unexpected dimension " << d);
299  }
300  }
301 }
302 // t += a otimes b
303 template<class T>
304 void
305 cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb)
306 {
307  for (size_t i = 0; i < na; i++)
308  for (size_t j = 0; j < nb; j++)
309  t(i,j) += a[i] * b[j];
310 }
311 template<class T>
314 {
315  point_basic<T> r;
316  r[0] = _x[i][0];
317  r[1] = _x[i][1];
318  r[2] = _x[i][2];
319  return r;
320 }
321 template<class T>
324 {
326  c[0] = _x[0][j];
327  c[1] = _x[1][j];
328  c[2] = _x[2][j];
329  return c;
330 }
331 template<class T>
332 bool
334 {
335  T det = determinant(A,3);
336  if (1+det == 1) return false;
337  T invdet = 1.0/det;
338  result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
339  result(0,1) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
340  result(0,2) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
341  result(1,0) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
342  result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
343  result(1,2) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
344  result(2,0) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
345  result(2,1) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
346  result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
347  return true;
348 }
349 template<class T>
350 bool
352  if (d <= 1) return true;
353  static T tol = 1e3*std::numeric_limits<T>::epsilon();
354  if (d <= 2) return fabs (_x[0][1] - _x[1][0]) < tol;
355  return (fabs (_x[0][1] - _x[1][0]) < tol)
356  && (fabs (_x[0][2] - _x[2][0]) < tol)
357  && (fabs (_x[1][2] - _x[2][1]) < tol);
358 }
359 // ------------------------------------------
360 // eigenvalues: the 3D case
361 // ------------------------------------------
362 template<class T>
363 static point_basic<T> eig3x3 (const tensor_basic<T>& a, tensor_basic<T>& q)
364 {
365  using namespace Eigen;
366  Matrix<T,3,3> a1;
367  for (size_t i = 0; i < 3; ++i)
368  for (size_t j = 0; j < 3; ++j) a1(i,j) = a(i,j);
369  SelfAdjointEigenSolver<Matrix<T,3,3> > es (a1);
370  point lambda;
371  for (size_t i = 0; i < 3; ++i) {
372  lambda[i] = es.eigenvalues()(i);
373  for (size_t j = 0; j < 3; ++j) {
374  q(i,j) = es.eigenvectors()(i,j);
375  }
376  }
377  return lambda;
378 }
379 // ------------------------------------------
380 // eigenvalues: the 2D case
381 // ------------------------------------------
382 template<class T>
383 static point_basic<T> eig2x2 (const tensor_basic<T>& a, tensor_basic<T>& q)
384 {
385  point_basic<T> d (0,0,0);
386  q.fill (0);
387  if (a(0,1) == 0) {
388  // a is already diagonal
389  if (a(0,0) >= a(1,1)) {
390  d[0] = a(0,0);
391  d[1] = a(1,1);
392  q(0,0) = q(1,1) = 1.;
393  q(0,1) = q(1,0) = 0.;
394  } else { // swap column 0 & 1:
395  d[1] = a(0,0);
396  d[0] = a(1,1);
397  q(0,0) = q(1,1) = 0.;
398  q(0,1) = q(1,0) = 1.;
399  }
400  return d;
401  }
402  // here a(0,1) != 0
403  T discr = sqr(a(1,1) - a(0,0)) + 4*sqr(a(0,1));
404  T trace = a(0,0) + a(1,1);
405  d[0] = (trace + sqrt(discr))/2;
406  d[1] = (trace - sqrt(discr))/2;
407  T lost_precision = 1e-6;
408  if (fabs(d[0]) + lost_precision*fabs(d[1]) == fabs(d[0])) { // d[1] == 0 up to machine precision
409  d[1] = 0.;
410  }
411  T c0 = (d[0]-a(0,0))/a(0,1);
412  T n0 = sqrt(1+sqr(c0));
413  q(0,0) = 1/n0;
414  q(1,0) = c0/n0;
415  T c1 = (d[1]-a(0,0))/a(0,1);
416  T n1 = sqrt(1+sqr(c1));
417  q(0,1) = 1/n1;
418  q(1,1) = c1/n1;
419  return d;
420 }
421 // ------------------------------------------
422 // eigenvalues: general case
423 // ------------------------------------------
424 template<class T>
427 {
428  check_macro (is_symmetric(d), "eig: tensor should be symmetric");
429  switch (d) {
430  case 1 : {
432  q(0,0) = 1; d[0] = operator()(0,0); return d;
433  }
434  case 2 : return eig2x2 (*this, q);
435  default: return eig3x3 (*this, q);
436  }
437 }
438 template<class T>
440 tensor_basic<T>::eig (size_t d) const
441 {
442  tensor_basic<T> q;
443  return eig (q, d);
444 }
445 // ------------------------------------------
446 // svd: the 3D case
447 // ------------------------------------------
448 template<class T, int N>
449 static
451 eigen_svd (const tensor_basic<T>& a, tensor_basic<T>& u, tensor_basic<T>& v)
452 {
453  using namespace Eigen;
454  Matrix<T,N,N> a1;
455  for (size_t i = 0; i < N; ++i)
456  for (size_t j = 0; j < N; ++j) a1(i,j) = a(i,j);
457  JacobiSVD<Matrix<T,N,N> > svd (a1, ComputeFullU | ComputeFullV);
458  point s;
459  for (size_t i = 0; i < N; ++i) {
460  s[i] = svd.singularValues()(i);
461  for (size_t j = 0; j < N; ++j) {
462  u(i,j) = svd.matrixU()(i,j);
463  v(i,j) = svd.matrixV()(i,j);
464  }
465  }
466  return s;
467 }
468 template<class T>
471 {
472  if (d == 1) {
473  u(0,0) = v(0,0) = 1;
474  return point(operator()(0,0));
475  }
476  if (d == 2) return eigen_svd<T,2> (*this, u, v);
477  else return eigen_svd<T,3> (*this, u, v);
478 }
479 // ----------------------------------------------------------------------------
480 // instanciation in library
481 // ----------------------------------------------------------------------------
482 #define _RHEOLEF_instanciation(T) \
483 template class tensor_basic<T>; \
484 template point_basic<T> operator* (const point_basic<T>& x, const tensor_basic<T>& a); \
485 template tensor_basic<T> trans (const tensor_basic<T>&, size_t); \
486 template tensor_basic<T> inv (const tensor_basic<T>&, size_t); \
487 template void prod (const tensor_basic<T>& a, const tensor_basic<T>& b, tensor_basic<T>& result, \
488  size_t di, size_t dj, size_t dk); \
489 template T ddot (const tensor_basic<T>& a, const tensor_basic<T> & b); \
490 template void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb); \
491 template bool invert_3x3 (const tensor_basic<T>& A, tensor_basic<T>& result); \
492 
494 
495 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
see the point page for the full documentation
std::ostream & put(std::ostream &s, size_type d=3) const
Definition: tensor.cc:37
std::istream & get(std::istream &)
Definition: tensor.cc:51
tensor_basic< T > & operator+=(const tensor_basic< T > &)
Definition: tensor.cc:174
point_basic< T > col(size_type i) const
Definition: tensor.cc:323
T determinant(size_type d=3) const
Definition: tensor.cc:288
tensor_basic< T > & operator*=(const T &k)
Definition: tensor.cc:192
tensor_basic< T > operator*(const tensor_basic< T > &b) const
Definition: tensor.cc:270
bool is_symmetric(size_type d=3) const
Definition: tensor.cc:351
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition: tensor.cc:426
tensor_basic< T > operator-() const
Definition: tensor.cc:110
void fill(const T &init_val)
Definition: tensor.h:252
const tensor_basic< T > & operator+() const
Definition: tensor.h:136
tensor_basic< T > & operator-=(const tensor_basic< T > &)
Definition: tensor.cc:183
bool operator==(const tensor_basic< T > &) const
Definition: tensor.cc:101
tensor_basic< T > & operator/=(const T &k)
Definition: tensor.cc:201
point_basic< T > row(size_type i) const
Definition: tensor.cc:313
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition: tensor.cc:470
point_basic< Float > point
Definition: point.h:163
rheolef::space_base_rep< T, M > t
size_t size_type
Definition: basis_get.cc:76
point u(const point &x)
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition: tensor.cc:219
t operator()(const t &a, const t &b)
Definition: space.cc:386
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
U determinant(const tensor_basic< U > &A, size_t d=3)
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:305
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition: tensor.cc:333
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
void prod(const tensor_basic< T > &a, const tensor_basic< T > &b, tensor_basic< T > &result, size_t di, size_t dj, size_t dk)
Definition: tensor.cc:256
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
size_t N
Definition: leveque.h:25
Float epsilon