Loading...
Searching...
No Matches
tensor.cc
Go to the documentation of this file.
//
// 2-tensor
//
// author: Pierre.Saramito@imag.fr
//
#include "rheolef/tensor.h"
#include "rheolef/compiler_eigen.h"
using namespace rheolef;
using namespace std;
namespace rheolef {
// output
template<class T>
ostream&
tensor_basic<T>::put (ostream& out, size_type d) const
{
switch (d) {
case 0 : return out;
case 1 : return out << _x[0][0];
case 2 : return out << "[" << _x[0][0] << ", " << _x[0][1] << ";\n"
<< " " << _x[1][0] << ", " << _x[1][1] << "]";
default: return out << "[" << _x[0][0] << ", " << _x[0][1] << ", " << _x[0][2] << ";\n"
<< " " << _x[1][0] << ", " << _x[1][1] << ", " << _x[1][2] << ";\n"
<< " " << _x[2][0] << ", " << _x[2][1] << ", " << _x[2][2] << "]";
}
}
template<class T>
istream&
tensor_basic<T>::get (istream& in)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
_x[i][j] = 0.;
char c;
if (!in.good()) return in;
in >> ws >> c;
// no tensor specifier:
// 1d case: "3.14"
if (c != '[') {
in.unget();
return in >> _x[0][0];
}
// has a tensor specifier, as:
// [ 23, 17; 25, 2]
in >> _x[0][0] >> ws >> c;
if (c == ']') return in; // 1d case
if (c != ',') error_macro ("invalid tensor input: expect `,'");
in >> _x[0][1];
in >> ws >> c;
size_type d = 3;
if (c == ',') d = 3;
else if (c == ';') d = 2;
else error_macro ("invalid tensor input: expect `,' or `;'");
if (d == 3) {
in >> _x[0][2] >> ws >> c;
if (c != ';') error_macro ("invalid tensor input: expect `;'");
}
in >> _x[1][0] >> ws >> c;
if (c != ',') error_macro ("invalid tensor input: expect `,'");
in >> _x[1][1] >> ws >> c;
if (d == 2) {
if (c != ']') error_macro ("invalid tensor input: expect `]'");
return in;
}
// d == 3
if (c != ',') error_macro ("invalid tensor input: expect `,'");
in >> _x[1][2] >> ws >> c;
if (c != ';') error_macro ("invalid tensor input: expect `;'");
in >> _x[2][0] >> ws >> c;
if (c != ',') error_macro ("invalid tensor input: expect `,'");
in >> _x[2][1] >> ws >> c;
if (c != ',') error_macro ("invalid tensor input: expect `,'");
in >> _x[2][2] >> ws >> c;
if (c != ']') error_macro ("invalid tensor input: expect `]'");
return in;
}
template<class T>
bool
tensor_basic<T>::operator== (const tensor_basic<T>& b) const
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
if (_x[i][j] != b._x[i][j]) return false;
return true;
}
template<class T>
tensor_basic<T>
tensor_basic<T>::operator- () const
{
tensor_basic<T> b;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
b._x[i][j] = - _x[i][j];
return b;
}
template<class T>
tensor_basic<T>
tensor_basic<T>::operator+ (const tensor_basic<T>& b) const
{
tensor_basic<T> c;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
c._x[i][j] = _x[i][j] + b._x[i][j];
return c;
}
template<class T>
tensor_basic<T>
tensor_basic<T>::operator- (const tensor_basic<T>& b) const
{
tensor_basic<T> c;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
c._x[i][j] = _x[i][j] - b._x[i][j];
return c;
}
template<class T>
tensor_basic<T>
tensor_basic<T>::operator* (const T& k) const
{
tensor_basic<T> b;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
b._x[i][j] = _x[i][j] * k;
return b;
}
template<class T>
point_basic<T>
tensor_basic<T>::operator* (const point_basic<T>& x) const
{
point_basic<T> y;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
y [i] += _x[i][j] * x[j];
return y;
}
template<class T>
point_basic<T> operator* (const point_basic<T>& x, const tensor_basic<T>& a)
{
point_basic<T> y;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
y [i] += a._x[j][i] * x[j];
return y;
}
template<class T>
tensor_basic<T>&
tensor_basic<T>::operator+= (const tensor_basic<T>& b)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
_x[i][j] += b._x[i][j];
return *this;
}
template<class T>
tensor_basic<T>&
tensor_basic<T>::operator-= (const tensor_basic<T>& b)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
_x[i][j] -= b._x[i][j];
return *this;
}
template<class T>
tensor_basic<T>&
tensor_basic<T>::operator*= (const T& k)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
_x[i][j] *= k;
return *this;
}
template<class T>
tensor_basic<T>&
tensor_basic<T>::operator/= (const T& k)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
_x[i][j] /= k;
return *this;
}
template<class T>
tensor_basic<T> trans (const tensor_basic<T>& a, size_t d)
{
tensor_basic<T> b;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < d; i++)
for (size_type j = 0; j < d; j++)
b._x[i][j] = a._x[j][i];
return b;
}
template<class T>
tensor_basic<T> inv (const tensor_basic<T>& a, size_t d)
{
tensor_basic<T> b;
typedef typename tensor_basic<T>::size_type size_type;
switch (d) {
case 0:
break;
case 1:
b(0,0) = 1/a(0,0);
break;
case 2: {
T det = a.determinant(2);
b(0,0) = a(1,1)/det;
b(1,1) = a(0,0)/det;
b(0,1) = - a(0,1)/det;
b(1,0) = - a(1,0)/det;
break;
}
case 3:
default: {
T det = a.determinant(3);
b(0,0) = (a(1,1)*a(2,2) - a(1,2)*a(2,1))/det;
b(0,1) = (a(0,2)*a(2,1) - a(0,1)*a(2,2))/det;
b(0,2) = (a(0,1)*a(1,2) - a(0,2)*a(1,1))/det;
b(1,0) = (a(1,2)*a(2,0) - a(1,0)*a(2,2))/det;
b(1,1) = (a(0,0)*a(2,2) - a(0,2)*a(2,0))/det;
b(1,2) = (a(0,2)*a(1,0) - a(0,0)*a(1,2))/det;
b(2,0) = (a(1,0)*a(2,1) - a(1,1)*a(2,0))/det;
b(2,1) = (a(0,1)*a(2,0) - a(0,0)*a(2,1))/det;
b(2,2) = (a(0,0)*a(1,1) - a(0,1)*a(1,0))/det;
break;
}
}
return b;
}
template<class T>
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)
{
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < di; i++)
for (size_type j = 0; j < dj; j++) {
T sum = 0;
for (size_type k = 0; k < dk; k++)
sum += a(i,k) * b(k,j);
result(i,j) = sum;
}
}
template<class T>
tensor_basic<T>
tensor_basic<T>::operator* (const tensor_basic<T>& b) const
{
tensor_basic<T> c;
prod (*this,b,c,3,3,3);
return c;
}
template<class T>
T ddot (const tensor_basic<T>& a, const tensor_basic<T> & b)
{
T r = 0;
typedef typename tensor_basic<T>::size_type size_type;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
r += a._x[i][j] * b._x[i][j];
return r;
}
template<class T>
T tensor_basic<T>::determinant (size_type d) const
{
switch (d) {
case 0 : return 1;
case 1 : return _x[0][0];
case 2 : return _x[0][0]*_x[1][1] - _x[0][1]*_x[1][0];
case 3 : return _x[0][0]*(_x[1][1]*_x[2][2] - _x[1][2]*_x[2][1])
- _x[0][1]*(_x[1][0]*_x[2][2] - _x[1][2]*_x[2][0])
+ _x[0][2]*(_x[1][0]*_x[2][1] - _x[1][1]*_x[2][0]);
default: {
error_macro ("determinant: unexpected dimension " << d);
}
}
}
// t += a otimes b
template<class T>
void
cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb)
{
for (size_t i = 0; i < na; i++)
for (size_t j = 0; j < nb; j++)
t(i,j) += a[i] * b[j];
}
template<class T>
point_basic<T>
tensor_basic<T>::row(size_type i) const
{
point_basic<T> r;
r[0] = _x[i][0];
r[1] = _x[i][1];
r[2] = _x[i][2];
return r;
}
template<class T>
point_basic<T>
tensor_basic<T>::col(size_type j) const
{
point_basic<T> c;
c[0] = _x[0][j];
c[1] = _x[1][j];
c[2] = _x[2][j];
return c;
}
template<class T>
bool
invert_3x3 (const tensor_basic<T>& A, tensor_basic<T>& result)
{
T det = determinant(A,3);
if (1+det == 1) return false;
T invdet = 1.0/det;
result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
result(0,1) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
result(0,2) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
result(1,0) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
result(1,2) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
result(2,0) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
result(2,1) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
return true;
}
template<class T>
bool
tensor_basic<T>::is_symmetric (size_t d) const {
if (d <= 1) return true;
static T tol = 1e3*std::numeric_limits<T>::epsilon();
if (d <= 2) return fabs (_x[0][1] - _x[1][0]) < tol;
return (fabs (_x[0][1] - _x[1][0]) < tol)
&& (fabs (_x[0][2] - _x[2][0]) < tol)
&& (fabs (_x[1][2] - _x[2][1]) < tol);
}
// ------------------------------------------
// eigenvalues: the 3D case
// ------------------------------------------
template<class T>
static point_basic<T> eig3x3 (const tensor_basic<T>& a, tensor_basic<T>& q)
{
using namespace Eigen;
Matrix<T,3,3> a1;
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j) a1(i,j) = a(i,j);
SelfAdjointEigenSolver<Matrix<T,3,3> > es (a1);
point lambda;
for (size_t i = 0; i < 3; ++i) {
lambda[i] = es.eigenvalues()(i);
for (size_t j = 0; j < 3; ++j) {
q(i,j) = es.eigenvectors()(i,j);
}
}
return lambda;
}
// ------------------------------------------
// eigenvalues: the 2D case
// ------------------------------------------
template<class T>
static point_basic<T> eig2x2 (const tensor_basic<T>& a, tensor_basic<T>& q)
{
point_basic<T> d (0,0,0);
q.fill (0);
if (a(0,1) == 0) {
// a is already diagonal
if (a(0,0) >= a(1,1)) {
d[0] = a(0,0);
d[1] = a(1,1);
q(0,0) = q(1,1) = 1.;
q(0,1) = q(1,0) = 0.;
} else { // swap column 0 & 1:
d[1] = a(0,0);
d[0] = a(1,1);
q(0,0) = q(1,1) = 0.;
q(0,1) = q(1,0) = 1.;
}
return d;
}
// here a(0,1) != 0
T discr = sqr(a(1,1) - a(0,0)) + 4*sqr(a(0,1));
T trace = a(0,0) + a(1,1);
d[0] = (trace + sqrt(discr))/2;
d[1] = (trace - sqrt(discr))/2;
T lost_precision = 1e-6;
if (fabs(d[0]) + lost_precision*fabs(d[1]) == fabs(d[0])) { // d[1] == 0 up to machine precision
d[1] = 0.;
}
T c0 = (d[0]-a(0,0))/a(0,1);
T n0 = sqrt(1+sqr(c0));
q(0,0) = 1/n0;
q(1,0) = c0/n0;
T c1 = (d[1]-a(0,0))/a(0,1);
T n1 = sqrt(1+sqr(c1));
q(0,1) = 1/n1;
q(1,1) = c1/n1;
return d;
}
// ------------------------------------------
// eigenvalues: general case
// ------------------------------------------
template<class T>
point_basic<T>
tensor_basic<T>::eig (tensor_basic<T>& q, size_t d) const
{
check_macro (is_symmetric(d), "eig: tensor should be symmetric");
switch (d) {
case 1 : {
point_basic<T> d;
q(0,0) = 1; d[0] = operator()(0,0); return d;
}
case 2 : return eig2x2 (*this, q);
default: return eig3x3 (*this, q);
}
}
template<class T>
point_basic<T>
tensor_basic<T>::eig (size_t d) const
{
tensor_basic<T> q;
return eig (q, d);
}
// ------------------------------------------
// svd: the 3D case
// ------------------------------------------
template<class T, int N>
static
point_basic<T>
eigen_svd (const tensor_basic<T>& a, tensor_basic<T>& u, tensor_basic<T>& v)
{
using namespace Eigen;
Matrix<T,N,N> a1;
for (size_t i = 0; i < N; ++i)
for (size_t j = 0; j < N; ++j) a1(i,j) = a(i,j);
JacobiSVD<Matrix<T,N,N> > svd (a1, ComputeFullU | ComputeFullV);
point s;
for (size_t i = 0; i < N; ++i) {
s[i] = svd.singularValues()(i);
for (size_t j = 0; j < N; ++j) {
u(i,j) = svd.matrixU()(i,j);
v(i,j) = svd.matrixV()(i,j);
}
}
return s;
}
template<class T>
point_basic<T>
tensor_basic<T>::svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t d) const
{
if (d == 1) {
u(0,0) = v(0,0) = 1;
return point(operator()(0,0));
}
if (d == 2) return eigen_svd<T,2> (*this, u, v);
else return eigen_svd<T,3> (*this, u, v);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class tensor_basic<T>; \
template point_basic<T> operator* (const point_basic<T>& x, const tensor_basic<T>& a); \
template tensor_basic<T> trans (const tensor_basic<T>&, size_t); \
template tensor_basic<T> inv (const tensor_basic<T>&, size_t); \
template 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); \
template T ddot (const tensor_basic<T>& a, const tensor_basic<T> & b); \
template void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb); \
template bool invert_3x3 (const tensor_basic<T>& A, tensor_basic<T>& result); \
_RHEOLEF_instanciation(Float)
}// namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
point
see the point page for the full documentation
rheolef::point_basic
Definition point.h:86
rheolef::tensor_basic
Definition tensor.h:90
rheolef::tensor_basic::put
std::ostream & put(std::ostream &s, size_type d=3) const
Definition tensor.cc:37
rheolef::tensor_basic::get
std::istream & get(std::istream &)
Definition tensor.cc:51
rheolef::tensor_basic::operator+=
tensor_basic< T > & operator+=(const tensor_basic< T > &)
Definition tensor.cc:174
rheolef::tensor_basic::col
point_basic< T > col(size_type i) const
Definition tensor.cc:323
rheolef::tensor_basic::determinant
T determinant(size_type d=3) const
Definition tensor.cc:288
rheolef::tensor_basic::size_type
size_t size_type
Definition tensor.h:93
rheolef::tensor_basic::operator*=
tensor_basic< T > & operator*=(const T &k)
Definition tensor.cc:192
rheolef::tensor_basic::operator*
tensor_basic< T > operator*(const tensor_basic< T > &b) const
Definition tensor.cc:270
rheolef::tensor_basic::is_symmetric
bool is_symmetric(size_type d=3) const
Definition tensor.cc:351
rheolef::tensor_basic::eig
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition tensor.cc:426
rheolef::tensor_basic::operator+
const tensor_basic< T > & operator+() const
Definition tensor.h:136
rheolef::tensor_basic::operator-
tensor_basic< T > operator-() const
Definition tensor.cc:110
rheolef::tensor_basic::fill
void fill(const T &init_val)
Definition tensor.h:252
rheolef::tensor_basic::operator-=
tensor_basic< T > & operator-=(const tensor_basic< T > &)
Definition tensor.cc:183
rheolef::tensor_basic::operator==
bool operator==(const tensor_basic< T > &) const
Definition tensor.cc:101
rheolef::tensor_basic::operator/=
tensor_basic< T > & operator/=(const T &k)
Definition tensor.cc:201
rheolef::tensor_basic::row
point_basic< T > row(size_type i) const
Definition tensor.cc:313
rheolef::tensor_basic::svd
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition tensor.cc:470
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Eigen
Definition field_eigen.h:35
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::inv
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition tensor.cc:219
rheolef::ddot
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
rheolef::determinant
U determinant(const tensor_basic< U > &A, size_t d=3)
rheolef::prod
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::cumul_otimes
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
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
rheolef::operator()
t operator()(const t &a, const t &b)
Definition space.cc:386
rheolef::invert_3x3
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
std
STL namespace.
N
size_t N
Definition neumann-laplace-check.cc:24
d
size_t d
Definition neumann-laplace-lambda.cc:28
A
Definition diffusion_isotropic.h:25
lambda
Definition yield_slip_circle.h:34
u
Definition leveque.h:25