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);
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;
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");
432 q(0,0) = 1;
d[0] = operator()(0,0);
return d;
434 case 2 :
return eig2x2 (*
this, q);
435 default:
return eig3x3 (*
this, q);
448 template<
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) \
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); \
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
tensor_basic< T > operator-() const
void fill(const T &init_val)
const tensor_basic< T > & operator+() const
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
point_basic< Float > point
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)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
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)
_RHEOLEF_instanciation(Float) _RHEOLEF_instanciation_evaluate(Float
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
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)