tensor.h

Go to the documentation of this file.

# ifndef _RHEOLEF_TENSOR_H
# define _RHEOLEF_TENSOR_H

#include "rheolef/point.h"
namespace rheolef {

template<class T>
class tensor_basic {
    public:

        typedef size_t size_type;
        typedef T      element_type;
        typedef T      float_type;


        tensor_basic (const T& init_val = 0);
        tensor_basic (T x[3][3]);
        tensor_basic (const tensor_basic<T>& a);
        static tensor_basic<T> identity (size_type d = 3);


        tensor_basic<T>& operator = (const tensor_basic<T>& a);
        tensor_basic<T>& operator = (const T& val);


        void fill (const T& init_val);
        void reset ();
        void set_row    (const point_basic<T>& r, size_t i, size_t d = 3);
        void set_column (const point_basic<T>& c, size_t j, size_t d = 3);


        T& operator()(size_type i, size_type j);
        T  operator()(size_type i, size_type j) const;
        point_basic<T>  row(size_type i) const;
        point_basic<T>  col(size_type i) const;
        size_t nrow() const; // = 3, for template matrix compatibility
        size_t ncol() const;

// inputs/outputs:

        std::ostream& put (std::ostream& s, size_type d = 3) const;
        std::istream& get (std::istream&);


        bool operator== (const tensor_basic<T>&) const;
        bool operator!= (const tensor_basic<T>& b) const { return ! operator== (b); }
        template <class U>
        friend tensor_basic<U> operator- (const tensor_basic<U>&);
        template <class U>
        friend tensor_basic<U> operator+ (const tensor_basic<U>&, const tensor_basic<U>&);
        template <class U>
        friend tensor_basic<U> operator- (const tensor_basic<U>&, const tensor_basic<U>&);
        template <class U>
        friend tensor_basic<U> operator* (int k, const tensor_basic<U>& a);
        template <class U>
        friend tensor_basic<U> operator* (const U& k, const tensor_basic<U>& a);
        template <class U>
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, int k);
        template <class U>
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, const U& k);
        template <class U>
        friend tensor_basic<U> operator/ (const tensor_basic<U>& a, int k);
        template <class U>
        friend tensor_basic<U> operator/ (const tensor_basic<U>& a, const U& k);
        template <class U>
        friend point_basic<U>  operator* (const tensor_basic<U>&, const point_basic<U>&);
        template <class U>
        friend point_basic<U>  operator* (const point_basic<U>& yt, const tensor_basic<U>& a);
               point_basic<T>  trans_mult (const point_basic<T>& x) const;
        template <class U>
        friend tensor_basic<U> trans      (const tensor_basic<U>& a, size_t d = 3);
        template <class U>
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, const tensor_basic<U>& b);
        template <class U>
        friend void prod (const tensor_basic<U>& a, const tensor_basic<U>& b, tensor_basic<U>& result,
                size_t di=3, size_t dj=3, size_t dk=3);
        template <class U>
        friend tensor_basic<U> inv        (const tensor_basic<U>& a, size_t d = 3);
        template <class U>
        friend tensor_basic<U> diag (const point_basic<U>& d);
        template <class U>
        friend tensor_basic<U> identity (size_t d=3);
        template <class U>
        friend tensor_basic<U> dyadic (const point_basic<U>& u, const point_basic<U>& v, size_t d=3);


        template <class U>
        friend U dotdot (const tensor_basic<U>&, const tensor_basic<U>&);
        template <class U>
        friend U norm2 (const tensor_basic<U>& a) { return dotdot(a,a); }
        template <class U>
        friend U dist2 (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm2(a-b); }
        template <class U>
        friend U norm  (const tensor_basic<U>& a) { return ::sqrt(norm2(a)); }
        template <class U>
        friend U dist  (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm(a-b); }
        T determinant (size_type d = 3) const;
        template <class U>
        friend U determinant (const tensor_basic<U>& A, size_t d = 3);
        template <class U>
        friend bool invert_3x3 (const tensor_basic<U>& A, tensor_basic<U>& result);

        // and   d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
        point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const;
        point_basic<T> eig (size_t dim = 3) const;

        // and   s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
        point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const;

        T _x[3][3];
};
typedef tensor_basic<Float> tensor;

// inputs/outputs:
template<class T>
inline
std::istream& operator>> (std::istream& in, tensor_basic<T>& a)
{
    return a.get (in);
}
template<class T>
inline
std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a)
{
    return a.put (out); 
}
template<class T>
tensor_basic<T> otimes (const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
template<class T>
void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
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);
@endcode

template<class T> struct  float_traits<tensor_basic<T> > { typedef typename float_traits<T>::type type; };
template<class T> struct scalar_traits<tensor_basic<T> > { typedef T type; };

template<class T>
inline
void
tensor_basic<T>::fill (const T& init_val)
{ 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = init_val;
}
template<class T>
inline
void
tensor_basic<T>::reset ()
{
    fill (0);
}
template<class T>
inline
tensor_basic<T>::tensor_basic (const T& init_val)
{ 
    fill (init_val);
}
template<class T>
inline
tensor_basic<T>::tensor_basic (T x[3][3])
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = x[i][j];
}
template<class T>
inline
tensor_basic<T>::tensor_basic (const tensor_basic<T>& a)
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = a._x[i][j];
}
template<class T>
inline
tensor_basic<T>&
tensor_basic<T>::operator= (const tensor_basic<T>& a)
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = a._x[i][j];
    return *this;
}
template<class T>
inline
tensor_basic<T>&
tensor_basic<T>::operator= (const T& val)
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = val;
    return *this;
}
template<class T>
inline
size_t
tensor_basic<T>::nrow() const
{
  return 3;
}
template<class T>
inline
size_t
tensor_basic<T>::ncol() const
{
  return 3;
}
template<class T>
inline
T&
tensor_basic<T>::operator()(size_type i, size_type j)  
{
    return _x[i%3][j%3];
}
template<class T>
inline
T
tensor_basic<T>::operator()(size_type i, size_type j) const
{
    return _x[i%3][j%3];
}
template<class T>
inline
tensor_basic<T>
operator* (int k, const tensor_basic<T>& a)
{
    return T(k)*a;
}
template<class T>
inline
tensor_basic<T>
operator* (const tensor_basic<T>& a, int k)
{
    return T(k)*a;
}
template<class T>
inline
tensor_basic<T>
operator* (const tensor_basic<T>& a, const T& k)
{
    return k*a;
}
template<class T>
inline
tensor_basic<T>
operator/ (const tensor_basic<T>& a, int k)
{
    return (1.0/T(k))*a;
}
template<class T>
inline
tensor_basic<T>
operator/ (const tensor_basic<T>& a, const T& k)
{
    return (1/k)*a;
}
template<class T>
inline
point_basic<T>
tensor_basic<T>::trans_mult (const point_basic<T>& x) const
{
    return x*(*this);
}
template<class T>
inline
void
cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t n)
{
    cumul_otimes (t, a, b, n, n);
}
template<class T>
inline
tensor_basic<T>
otimes (const point_basic<T>& a, const point_basic<T>& b, size_t n)
{
    tensor_basic<T> t;
    cumul_otimes (t, a, b, n, n);
    return t;
}
template<class T>
inline
T
determinant (const tensor_basic<T>& A, size_t d)
{ 
    return A.determinant (d);
}
template<class T>
inline
tensor_basic<T>
diag (const point_basic<T>& d)
{
  tensor_basic<T> a;
  a(0,0) = d[0];
  a(1,1) = d[1];
  a(2,2) = d[2];
  return a;
}
template<class T>
inline
tensor_basic<T>
identity (size_t d)
{
  tensor_basic<T> a;
  for (size_t i = 0; i < d; i++) a(i,i) = 1;
  return a;
}
template<class T>
inline
tensor_basic<T> 
dyadic (const point_basic<T>& u, const point_basic<T>& v, size_t d)
{
  tensor_basic<T> a;
  for (size_t i = 0; i < d; i++) 
    for (size_t j = 0; j < d; j++) 
      a(i,j) = u[i]*v[j];
  return a;
}
template<class T>
inline
void
tensor_basic<T>::set_column (const point_basic<T>& c, size_t j, size_t d)
{
  for (size_t i = 0; i < d; i++)
    operator()(i,j) = c[i];
}
template<class T>
inline
void
tensor_basic<T>::set_row (const point_basic<T>& r, size_t i, size_t d)
{
  for (size_t j = 0; j < d; j++)
    operator()(i,j) = r[j];
}
template<class T>
inline
tensor_basic<T>
tensor_basic<T>::identity (size_type d)
{
  tensor_basic<T> I;
  for (size_t i = 0; i < d; i++)
    I(i,i) = 1;
  return I;
}
}// namespace rheolef
# endif /* _RHEOLEF_TENSOR_H */