26#include "rheolef/tensor.h"
27#include "rheolef/compiler_eigen.h"
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] <<
"]";
57 if (!in.good())
return in;
63 return in >> _x[0][0];
67 in >> _x[0][0] >> ws >> c;
68 if (c ==
']')
return in;
69 if (c !=
',')
error_macro (
"invalid tensor input: expect `,'");
74 else if (c ==
';')
d = 2;
75 else error_macro (
"invalid tensor input: expect `,' or `;'");
77 in >> _x[0][2] >> ws >> c;
78 if (c !=
';')
error_macro (
"invalid tensor input: expect `;'");
80 in >> _x[1][0] >> ws >> c;
81 if (c !=
',')
error_macro (
"invalid tensor input: expect `,'");
82 in >> _x[1][1] >> ws >> c;
84 if (c !=
']')
error_macro (
"invalid tensor input: expect `]'");
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 `]'");
105 if (_x[i][j] != b._x[i][j])
return false;
116 b._x[i][j] = - _x[i][j];
127 c._x[i][j] = _x[i][j] + b._x[i][j];
138 c._x[i][j] = _x[i][j] - b._x[i][j];
148 b._x[i][j] = _x[i][j] * k;
159 y [i] += _x[i][j] * x[j];
169 y [i] += a._x[j][i] * x[j];
178 _x[i][j] += b._x[i][j];
187 _x[i][j] -= b._x[i][j];
215 b._x[i][j] = a._x[j][i];
230 T det = a.determinant(2);
233 b(0,1) = - a(0,1)/det;
234 b(1,0) = - a(1,0)/det;
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;
257 size_t di,
size_t dj,
size_t dk)
264 sum += a(i,k) * b(k,j);
273 prod (*
this,b,c,3,3,3);
284 r += a._x[i][j] * b._x[i][j];
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]);
298 error_macro (
"determinant: unexpected dimension " <<
d);
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];
336 if (1+det == 1)
return false;
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;
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);
365 using namespace Eigen;
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);
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);
389 if (
a(0,0) >=
a(1,1)) {
392 q(0,0) = q(1,1) = 1.;
393 q(0,1) = q(1,0) = 0.;
397 q(0,0) = q(1,1) = 0.;
398 q(0,1) = q(1,0) = 1.;
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])) {
411 T c0 = (
d[0]-
a(0,0))/
a(0,1);
412 T n0 = sqrt(1+sqr(c0));
415 T c1 = (
d[1]-
a(0,0))/
a(0,1);
416 T n1 = sqrt(1+sqr(c1));
428 check_macro (is_symmetric(
d),
"eig: tensor should be symmetric");
434 case 2 :
return eig2x2 (*
this, q);
435 default:
return eig3x3 (*
this, q);
448template<
class T,
int N>
453 using namespace Eigen;
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);
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);
474 return point(
operator()(0,0));
476 if (
d == 2)
return eigen_svd<T,2> (*
this,
u, v);
477 else return eigen_svd<T,3> (*
this,
u, v);
482#define _RHEOLEF_instanciation(T) \
483template class tensor_basic<T>; \
484template point_basic<T> operator* (const point_basic<T>& x, const tensor_basic<T>& a); \
485template tensor_basic<T> trans (const tensor_basic<T>&, size_t); \
486template tensor_basic<T> inv (const tensor_basic<T>&, size_t); \
487template 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); \
489template T ddot (const tensor_basic<T>& a, const tensor_basic<T> & b); \
490template void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb); \
491template bool invert_3x3 (const tensor_basic<T>& A, tensor_basic<T>& result); \
#define _RHEOLEF_instanciation(T, M, A)
field::size_type size_type
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
std::istream & get(std::istream &)
tensor_basic< T > & operator+=(const tensor_basic< T > &)
point_basic< T > col(size_type i) const
T determinant(size_type d=3) const
tensor_basic< T > & operator*=(const T &k)
tensor_basic< T > operator*(const tensor_basic< T > &b) const
bool is_symmetric(size_type d=3) const
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
const tensor_basic< T > & operator+() const
tensor_basic< T > operator-() const
void fill(const T &init_val)
tensor_basic< T > & operator-=(const tensor_basic< T > &)
bool operator==(const tensor_basic< T > &) const
tensor_basic< T > & operator/=(const T &k)
point_basic< T > row(size_type i) const
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
#define error_macro(message)
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)
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
U determinant(const tensor_basic< U > &A, size_t d=3)
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)
void cumul_otimes(tensor_basic< T > &t, const point_basic< T > &a, const point_basic< T > &b, size_t na, size_t nb)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
t operator()(const t &a, const t &b)
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)