tensor (2rheolef)

The tensor class defines a 3*3 tensor, as the value of a tensorial valued field. Basic algebra with scalars, vectors of R^3 (i.e. the point class) and tensor objects are supported.

template<class T>
class tensor_basic {
    public:

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

// allocators:

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

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
        tensor_basic (const std::initializer_list<std::initializer_list<T> >& il);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

// affectation:

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

// modifiers:

        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);

// accessors:

        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&);

// algebra:

        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);
        // tr(a) = a00 + a11 + a22
        template <class U>
        friend U tr (const tensor_basic<U>& a, size_t d=3);
        // a = u otimes v <==> aij = ui*vj
        template <class U>
        friend tensor_basic<U> otimes (const point_basic<U>& u, const point_basic<U>& v, size_t d=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 point_basic<U> diag (const tensor_basic<U>& a);

// metric and geometric transformations:

        template <class U>
        friend U ddot (const tensor_basic<U>&, const tensor_basic<U>&);
        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);

// spectral:
        // eigenvalues & eigenvectors:
        // a = q*d*q^T
        // a may be symmetric
        // where q=(q1,q2,q3) are eigenvectors in rows (othonormal matrix)
        // and   d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
        // return d
        point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const;
        point_basic<T> eig (size_t dim = 3) const;

        // singular value decomposition:
        // a = u*s*v^T
        // a can be unsymmetric
        // where u=(u1,u2,u3) are left pseudo-eigenvectors in rows (othonormal matrix)
        //       v=(v1,v2,v3) are right pseudo-eigenvectors in rows (othonormal matrix)
        // and   s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
        // return s
        point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const;

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

// nonlinear algebra:
template<class T>
tensor_basic<T> exp (const tensor_basic<T>& a, size_t d = 3);

// 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);
}
// 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 = 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);