librheolef-dev 6.5-1+b1 / usr / include / rheolef / tensor.h

# ifndef _RHEOLEF_TENSOR_H
# define _RHEOLEF_TENSOR_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

/*Class:tensor
NAME: @code{tensor} - a N*N tensor, N=1,2,3
@cindex tensor
@clindex tensor
@clindex point
@clindex field
SYNOPSYS:
@noindent
The @code{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 @code{point} class) and @code{tensor} objects
are supported.
AUTHOR: Pierre.Saramito@imag.fr
DATE:   9 october 2003
End:
*/

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

//<tensor:
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);
//>tensor:

// -----------------------------------------------------------------------
// inlined
// -----------------------------------------------------------------------
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];
}
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
template<class T>
tensor_basic<T>::tensor_basic (const std::initializer_list<std::initializer_list<T> >& il) : _x() {
#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef typename std::initializer_list<std::initializer_list<T> >::const_iterator const_iterator;
    typedef typename std::initializer_list<T>::const_iterator                         const_iterator_row;
#else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    typedef const std::initializer_list<T>* const_iterator;
    typedef const T*                        const_iterator_row;
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
    fill (T());
    check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
    size_type i = 0;
    for (const_iterator iter = il.begin(); iter != il.end(); ++iter, ++i) {
      const std::initializer_list<T>& row = *iter;
      check_macro (row.size() <= 3, "unexpected initializer list size=" << row.size() << " > 3");
      size_type j = 0;
      for (const_iterator_row jter = row.begin(); jter != row.end(); ++jter, ++j) {
        _x[i][j] = *jter;
      }
    }
}
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
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>& u, const point_basic<T>& v, size_t d=3)
{
    tensor_basic<T> a;
    cumul_otimes (a, u, v, d, d);
    return a;
}
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
point_basic<T>
diag (const tensor_basic<T>& a)
{
  point_basic<T> d;
  d[0] = a(0,0);
  d[1] = a(1,1);
  d[2] = a(2,2);
  return d;
}
template <class T>
inline
T
tr (const tensor_basic<T>& a, size_t d = 3) {
  T sum = 0;
  for (size_t i = 0; i < d; i++) sum += a(i,i);
  return sum;
} 
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>::eye (size_type d)
{
  tensor_basic<T> I;
  for (size_t i = 0; i < d; i++)
    I(i,i) = 1;
  return I;
}
template <class T>
inline
T
norm2 (const tensor_basic<T>& a)
{
  return ddot(a,a);
}
template <class T>
inline
T
dist2 (const tensor_basic<T>& a, const tensor_basic<T>& b)
{
  return norm2(a-b);
}
template <class U>
inline
U
norm  (const tensor_basic<U>& a)
{
  return ::sqrt(norm2(a));
}
template <class U>
inline
U
dist (const tensor_basic<U>& a, const tensor_basic<U>& b)
{
  return norm(a-b);
}

}// namespace rheolef
# endif /* _RHEOLEF_TENSOR_H */