librheolef-dev 5.93-2 / 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:
class tensor {
    public:

	typedef size_t size_type;
	typedef Float  element_type;

// allocators:

	tensor (const Float& init_val = 0);
	tensor (Float x[3][3]);
	tensor (const tensor& a);

// affectation:

	tensor& operator = (const tensor& a);
	tensor& operator = (const Float& val);

// modifiers:

	void fill (const Float& init_val);
	void reset ();
	void set_row    (const point& r, size_t i, size_t d = 3);
	void set_column (const point& c, size_t j, size_t d = 3);

// accessors:

	Float& operator()(size_type i, size_type j);
	Float  operator()(size_type i, size_type j) const;
	point  row(size_type i) const;
	point  col(size_type i) const;
	size_t nrow() const; // = 3, for template matrix compatibility
	size_t ncol() const;

// inputs/outputs:

	friend std::istream& operator >> (std::istream& in, tensor& a);
	friend std::ostream& operator << (std::ostream& out, const tensor& a);
	std::ostream& put (std::ostream& s, size_type d = 3) const;

// algebra:

	friend bool  operator == (const tensor&, const tensor&);
	friend tensor operator - (const tensor&);
	friend tensor operator + (const tensor&, const tensor&);
	friend tensor operator - (const tensor&, const tensor&);
	friend tensor operator * (Float, const tensor&);
	friend tensor operator * (const tensor& a, Float k);
	friend tensor operator / (const tensor& a, Float k);
	friend point  operator * (const tensor& a, const point& x);
	friend point  operator * (const point& yt, const tensor& a);
	       point  trans_mult (const point& x) const;
	friend tensor trans      (const tensor& a, size_type d = 3);
	friend tensor operator * (const tensor& a, const tensor& b);
	friend tensor inv        (const tensor& a, size_type d = 3);
	friend tensor diag (const point& d);

// metric and geometric transformations:

	friend Float dotdot (const tensor&, const tensor&);
	friend Float norm2 (const tensor& a) { return dotdot(a,a); }
	friend Float dist2 (const tensor& a, const tensor& b) { return norm2(a-b); }
	friend Float norm  (const tensor& a) { return ::sqrt(norm2(a)); }
	friend Float dist  (const tensor& a, const tensor& b) { return norm(a-b); }
	Float determinant (size_type d = 3) const;
	friend Float determinant (const tensor& A, size_t d = 3);
	friend bool invert_3x3 (const tensor& A, tensor& 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 eig (tensor& q, size_t dim = 3) const;
	point 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 svd (tensor& u, tensor& v, size_t dim = 3) const;

// data:
	Float _x[3][3];
// internal:
	std::istream& get (std::istream&);
};
// t = a otimes b
tensor otimes (const point& a, const point& b, size_t na = 3);
// t += a otimes b
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na = 3);
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na, size_t nb);
//>tensor:

// -----------------------------------------------------------------------
// inlined
// -----------------------------------------------------------------------

inline
void
tensor::fill (const Float& init_val)
{ 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = init_val;
}
inline
void
tensor::reset ()
{
    fill (0);
}
inline
tensor::tensor (const Float& init_val)
{ 
    fill (init_val);
}
inline
tensor::tensor (Float 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];
}
inline
tensor::tensor (const tensor& a)
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = a._x[i][j];
}
inline
tensor&
tensor::operator = (const tensor& 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;
}
inline
tensor&
tensor::operator = (const Float& val)
{
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
      _x[i][j] = val;
    return *this;
}
inline
size_t
tensor::nrow() const
{
  return 3;
}
inline
size_t
tensor::ncol() const
{
  return 3;
}
inline
Float&
tensor::operator()(size_type i, size_type j)  
{
    return _x[i%3][j%3];
}
inline
Float
tensor::operator()(size_type i, size_type j) const
{
    return _x[i%3][j%3];
}
inline
std::istream& operator >> (std::istream& in, tensor& a)
{
    return a.get (in);
}
inline
std::ostream& operator << (std::ostream& out, const tensor& a)
{
    return a.put (out); 
}
inline
tensor
operator * (const tensor& a, Float k)
{
    return k*a;
}
inline
tensor
operator / (const tensor& a, Float k)
{
    return (1/k)*a;
}
inline
point
tensor::trans_mult (const point& x) const
{
    return x*(*this);
}
inline
void
cumul_otimes (tensor& t, const point& a, const point& b, size_t n)
{
    cumul_otimes (t, a, b, n, n);
}
inline
tensor
otimes (const point& a, const point& b, size_t n)
{
    tensor t;
    cumul_otimes (t, a, b, n, n);
    return t;
}
inline
Float
determinant (const tensor& A, size_t d)
{ 
    return A.determinant (d);
}
inline
tensor
diag (const point& d)
{
  tensor a;
  a(0,0) = d[0];
  a(1,1) = d[1];
  a(2,2) = d[2];
  return a;
}
inline
void
tensor::set_column (const point& c, size_t j, size_t d)
{
  for (size_t i = 0; i < d; i++)
    operator()(i,j) = c[i];
}
inline
void
tensor::set_row (const point& r, size_t i, size_t d)
{
  for (size_t j = 0; j < d; j++)
    operator()(i,j) = r[j];
}
}// namespace rheolef
# endif /* _RHEOLEF_TENSOR_H */