Loading...
Searching...
No Matches
tensor.h
Go to the documentation of this file.
# 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
//
// =========================================================================
// AUTHOR: Pierre.Saramito@imag.fr
// DATE: 9 october 2003
namespace rheolef {
} // namespace rheolef
#include "rheolef/point.h"
namespace rheolef {
// [verbatim_tensor_basic]
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);
tensor_basic (const std::initializer_list<std::initializer_list<T> >& il);
// 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);
const 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); }
const tensor_basic<T>& operator+ () const { return *this; }
tensor_basic<T> operator- () const;
tensor_basic<T> operator+ (const tensor_basic<T>& b) const;
tensor_basic<T> operator- (const tensor_basic<T>& b) const;
tensor_basic<T> operator* (const tensor_basic<T>& b) const;
tensor_basic<T> operator* (const T& k) const;
tensor_basic<T> operator/ (const T& k) const;
point_basic<T> operator* (const point_basic<T>&) const;
point_basic<T> trans_mult (const point_basic<T>& x) const;
tensor_basic<T>& operator+= (const tensor_basic<T>&);
tensor_basic<T>& operator-= (const tensor_basic<T>&);
tensor_basic<T>& operator*= (const T& k);
tensor_basic<T>& operator/= (const T& k);
T determinant (size_type d = 3) const;
bool is_symmetric (size_type d = 3) const;
// 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;
// [verbatim_tensor_basic]
// data:
T _x[3][3];
// [verbatim_tensor_basic_cont]
};
// [verbatim_tensor_basic_cont]
// [verbatim_tensor]
typedef tensor_basic<Float> tensor;
// [verbatim_tensor]
// [verbatim_tensor_basic_cont2]
template <class U>
point_basic<U> operator* (const point_basic<U>& yt, const tensor_basic<U>& a);
template <class U>
tensor_basic<U> trans (const tensor_basic<U>& a, size_t d = 3);
template <class U>
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>
U tr (const tensor_basic<U>& a, size_t d=3);
template <class U>
U ddot (const tensor_basic<U>&, const tensor_basic<U>&);
// a = u otimes v <==> aij = ui*vj
template <class U>
tensor_basic<U> otimes (const point_basic<U>& u, const point_basic<U>& v, size_t d=3);
template <class U>
tensor_basic<U> inv (const tensor_basic<U>& a, size_t d = 3);
template <class U>
tensor_basic<U> diag (const point_basic<U>& d);
template <class U>
point_basic<U> diag (const tensor_basic<U>& a);
template <class U>
U determinant (const tensor_basic<U>& A, size_t d = 3);
template <class U>
bool invert_3x3 (const tensor_basic<U>& A, tensor_basic<U>& result);
// matrix exponential:
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);
// [verbatim_tensor_basic_cont2]
// -----------------------------------------------------------------------
// 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];
}
template<class T>
tensor_basic<T>::tensor_basic (const std::initializer_list<std::initializer_list<T> >& il) : _x() {
typedef typename std::initializer_list<std::initializer_list<T> >::const_iterator const_iterator;
typedef typename std::initializer_list<T>::const_iterator const_iterator_row;
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;
}
}
}
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
const T&
tensor_basic<T>::operator()(size_type i, size_type j) const
{
return _x[i%3][j%3];
}
template <class T, class U>
inline
typename
std::enable_if<
details::is_rheolef_arithmetic<U>::value
,tensor_basic<T>
>::type
operator* (const U& k, const tensor_basic<T>& a)
{
return a*k;
}
template<class T>
inline
tensor_basic<T>
tensor_basic<T>::operator/ (const T& k) const
{
return operator* (1./k);
}
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)
{
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) {
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
// -------------------------------------------------------------
// serialization
// -------------------------------------------------------------
#ifdef _RHEOLEF_HAVE_MPI
namespace boost {
namespace serialization {
template<class Archive, class T>
void serialize (Archive& ar, class rheolef::tensor_basic<T>& x, const unsigned int version) {
ar & x(0,0);
ar & x(0,1);
ar & x(0,2);
ar & x(1,0);
ar & x(1,1);
ar & x(1,2);
ar & x(2,0);
ar & x(2,1);
ar & x(2,2);
}
} // namespace serialization
} // namespace boost
// Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
// When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
// To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
namespace boost {
namespace mpi {
// TODO: when tensor_basic<T> is not a simple type, such as T=bigfloat or T=gmp, etc
template <>
struct is_mpi_datatype<rheolef::tensor_basic<double> > : mpl::true_ { };
} // namespace mpi
} // namespace boost
#endif // _RHEOLEF_HAVE_MPI
#endif /* _RHEOLEF_TENSOR_H */
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::operator()
T & operator()(size_type i, size_type j)
Definition tensor.h:335
rheolef::tensor_basic::operator=
tensor_basic< T > & operator=(const tensor_basic< T > &a)
Definition tensor.h:303
rheolef::tensor_basic::get
std::istream & get(std::istream &)
Definition tensor.cc:51
rheolef::tensor_basic::trans_mult
point_basic< T > trans_mult(const point_basic< T > &x) const
Definition tensor.h:367
rheolef::tensor_basic::operator+=
tensor_basic< T > & operator+=(const tensor_basic< T > &)
Definition tensor.cc:174
rheolef::tensor_basic::eye
static tensor_basic< T > eye(size_type d=3)
Definition tensor.h:443
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::tensor_basic
tensor_basic(const std::initializer_list< std::initializer_list< T > > &il)
Definition tensor.h:285
rheolef::tensor_basic::operator()
const T & operator()(size_type i, size_type j) const
Definition tensor.h:342
rheolef::tensor_basic::_x
T _x[3][3]
Definition tensor.h:175
rheolef::tensor_basic::operator*
tensor_basic< T > operator*(const tensor_basic< T > &b) const
Definition tensor.cc:270
rheolef::tensor_basic::tensor_basic
tensor_basic(const tensor_basic< T > &a)
Definition tensor.h:279
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::tensor_basic
tensor_basic(T x[3][3])
Definition tensor.h:272
rheolef::tensor_basic::operator/
tensor_basic< T > operator/(const T &k) const
Definition tensor.h:360
rheolef::tensor_basic::ncol
size_t ncol() const
Definition tensor.h:328
rheolef::tensor_basic::reset
void reset()
Definition tensor.h:260
rheolef::tensor_basic::fill
void fill(const T &init_val)
Definition tensor.h:252
rheolef::tensor_basic::operator!=
bool operator!=(const tensor_basic< T > &b) const
Definition tensor.h:135
rheolef::tensor_basic::nrow
size_t nrow() const
Definition tensor.h:321
rheolef::tensor_basic::set_row
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
Definition tensor.h:435
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::element_type
T element_type
Definition tensor.h:94
rheolef::tensor_basic::float_type
T float_type
Definition tensor.h:95
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::tensor_basic
tensor_basic(const T &init_val=0)
Definition tensor.h:266
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
rheolef::tensor_basic::set_column
void set_column(const point_basic< T > &c, size_t j, size_t d=3)
Definition tensor.h:427
rheolef::tensor
tensor_basic< Float > tensor
Definition tensor.h:181
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)")
boost::mpi
Definition mpi_pair_datatype.h:52
boost::serialization::serialize
void serialize(Archive &ar, class rheolef::point_basic< T > &x, const unsigned int version)
Definition point.h:579
boost
Definition mpi_pair_datatype.h:51
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::otimes
tensor_basic< U > otimes(const point_basic< U > &u, const point_basic< U > &v, size_t d=3)
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::tr
U tr(const tensor_basic< U > &a, size_t d=3)
rheolef::exp
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition tensor-exp.cc:92
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
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::dist2
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
rheolef::operator>>
std::istream & operator>>(std::istream &is, const catchmark &m)
Definition catchmark.h:88
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56
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::dist
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
rheolef::invert_3x3
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
d
size_t d
Definition neumann-laplace-lambda.cc:28
A
Definition diffusion_isotropic.h:25
rheolef::const_iterator
Definition field_expr_recursive.h:594
rheolef::float_traits< tensor_basic< T > >::type
float_traits< T >::type type
Definition tensor.h:246
rheolef::float_traits
helper for std::complex<T>: get basic T type
Definition Float.h:93
rheolef::float_traits::type
T type
Definition Float.h:94
rheolef::scalar_traits< tensor_basic< T > >::type
T type
Definition tensor.h:247
rheolef::scalar_traits
helper for point_basic<T> & tensor_basic<T>: get basic T type
Definition point.h:323
u
Definition leveque.h:25