rheolef  7.0
tensor.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_TENSOR_H
2 # define _RHEOLEF_TENSOR_H
3 
16 #include "rheolef/point.h"
17 namespace rheolef {
18 
19 template<class T>
20 class tensor_basic {
21  public:
22 
23  typedef size_t size_type;
24  typedef T element_type;
25  typedef T float_type;
26 
27 
28  tensor_basic (const T& init_val = 0);
29  tensor_basic (T x[3][3]);
31  static tensor_basic<T> eye (size_type d = 3);
32 
33 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
34  tensor_basic (const std::initializer_list<std::initializer_list<T> >& il);
35 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
36 
37 
39  tensor_basic<T>& operator= (const T& val);
40 
41 
42  void fill (const T& init_val);
43  void reset ();
44  void set_row (const point_basic<T>& r, size_t i, size_t d = 3);
45  void set_column (const point_basic<T>& c, size_t j, size_t d = 3);
46 
47 
48  T& operator()(size_type i, size_type j);
49  const T& operator()(size_type i, size_type j) const;
50  point_basic<T> row(size_type i) const;
51  point_basic<T> col(size_type i) const;
52  size_t nrow() const; // = 3, for template matrix compatibility
53  size_t ncol() const;
54 
55 // inputs/outputs:
56 
57  std::ostream& put (std::ostream& s, size_type d = 3) const;
58  std::istream& get (std::istream&);
59 
60 
61  bool operator== (const tensor_basic<T>&) const;
62  bool operator!= (const tensor_basic<T>& b) const { return ! operator== (b); }
63  const tensor_basic<T>& operator+ () const { return *this; }
68  tensor_basic<T> operator* (const T& k) const;
69  tensor_basic<T> operator/ (const T& k) const;
71  point_basic<T> trans_mult (const point_basic<T>& x) const;
72 
73 
74  T determinant (size_type d = 3) const;
75 
76  // and d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
77  point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const;
78  point_basic<T> eig (size_t dim = 3) const;
79 
80  // and s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
81  point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const;
82 
83  T _x[3][3];
84 };
86 
87 
88 template <class U>
90 template <class U>
91 tensor_basic<U> trans (const tensor_basic<U>& a, size_t d = 3);
92 template <class U>
93 void prod (const tensor_basic<U>& a, const tensor_basic<U>& b, tensor_basic<U>& result,
94  size_t di=3, size_t dj=3, size_t dk=3);
95 template <class U>
96 U tr (const tensor_basic<U>& a, size_t d=3);
97 template <class U>
98 U ddot (const tensor_basic<U>&, const tensor_basic<U>&);
99 // a = u otimes v <==> aij = ui*vj
100 template <class U>
101 tensor_basic<U> otimes (const point_basic<U>& u, const point_basic<U>& v, size_t d=3);
102 template <class U>
103 tensor_basic<U> inv (const tensor_basic<U>& a, size_t d = 3);
104 template <class U>
106 template <class U>
108 template <class U>
109 U determinant (const tensor_basic<U>& A, size_t d = 3);
110 template <class U>
111 bool invert_3x3 (const tensor_basic<U>& A, tensor_basic<U>& result);
112 
113 template<class T>
114 tensor_basic<T> exp (const tensor_basic<T>& a, size_t d = 3);
115 
116 // inputs/outputs:
117 template<class T>
118 inline
119 std::istream& operator>> (std::istream& in, tensor_basic<T>& a)
120 {
121  return a.get (in);
122 }
123 template<class T>
124 inline
125 std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a)
126 {
127  return a.put (out);
128 }
129 template<class T>
130 void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
131 template<class T>
132 void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb);
133 @endcode
134 
135 template<class T> struct float_traits<tensor_basic<T> > { typedef typename float_traits<T>::type type; };
136 template<class T> struct scalar_traits<tensor_basic<T> > { typedef T type; };
137 
138 template<class T>
139 inline
140 void
141 tensor_basic<T>::fill (const T& init_val)
142 {
143  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
144  _x[i][j] = init_val;
145 }
146 template<class T>
147 inline
148 void
150 {
151  fill (0);
152 }
153 template<class T>
154 inline
155 tensor_basic<T>::tensor_basic (const T& init_val)
156 {
157  fill (init_val);
158 }
159 template<class T>
160 inline
162 {
163  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
164  _x[i][j] = x[i][j];
165 }
166 template<class T>
167 inline
169 {
170  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
171  _x[i][j] = a._x[i][j];
172 }
173 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
174 template<class T>
175 tensor_basic<T>::tensor_basic (const std::initializer_list<std::initializer_list<T> >& il) : _x() {
176 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
177  typedef typename std::initializer_list<std::initializer_list<T> >::const_iterator const_iterator;
178  typedef typename std::initializer_list<T>::const_iterator const_iterator_row;
179 #else // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
180  typedef const std::initializer_list<T>* const_iterator;
181  typedef const T* const_iterator_row;
182 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_ITERATOR
183  fill (T());
184  check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
185  size_type i = 0;
186  for (const_iterator iter = il.begin(); iter != il.end(); ++iter, ++i) {
187  const std::initializer_list<T>& row = *iter;
188  check_macro (row.size() <= 3, "unexpected initializer list size=" << row.size() << " > 3");
189  size_type j = 0;
190  for (const_iterator_row jter = row.begin(); jter != row.end(); ++jter, ++j) {
191  _x[i][j] = *jter;
192  }
193  }
194 }
195 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
196 template<class T>
197 inline
200 {
201  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
202  _x[i][j] = a._x[i][j];
203  return *this;
204 }
205 template<class T>
206 inline
209 {
210  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
211  _x[i][j] = val;
212  return *this;
213 }
214 template<class T>
215 inline
216 size_t
218 {
219  return 3;
220 }
221 template<class T>
222 inline
223 size_t
225 {
226  return 3;
227 }
228 template<class T>
229 inline
230 T&
232 {
233  return _x[i%3][j%3];
234 }
235 template<class T>
236 inline
237 const T&
239 {
240  return _x[i%3][j%3];
241 }
242 template <class T, class U>
243 inline
244 typename
245 std::enable_if<
248 >::type
249 operator* (const U& k, const tensor_basic<T>& a)
250 {
251  return a*k;
252 }
253 template<class T>
254 inline
257 {
258  return operator* (1./k);
259 }
260 template<class T>
261 inline
264 {
265  return x*(*this);
266 }
267 template<class T>
268 inline
269 void
271 {
272  cumul_otimes (t, a, b, n, n);
273 }
274 template<class T>
275 inline
277 otimes (const point_basic<T>& u, const point_basic<T>& v, size_t d)
278 {
280  cumul_otimes (a, u, v, d, d);
281  return a;
282 }
283 template<class T>
284 inline
285 T
286 determinant (const tensor_basic<T>& A, size_t d)
287 {
288  return A.determinant (d);
289 }
290 template<class T>
291 inline
294 {
296  a(0,0) = d[0];
297  a(1,1) = d[1];
298  a(2,2) = d[2];
299  return a;
300 }
301 template<class T>
302 inline
305 {
307  d[0] = a(0,0);
308  d[1] = a(1,1);
309  d[2] = a(2,2);
310  return d;
311 }
312 template <class T>
313 inline
314 T
315 tr (const tensor_basic<T>& a, size_t d) {
316  T sum = 0;
317  for (size_t i = 0; i < d; i++) sum += a(i,i);
318  return sum;
319 }
320 template<class T>
321 inline
322 void
323 tensor_basic<T>::set_column (const point_basic<T>& c, size_t j, size_t d)
324 {
325  for (size_t i = 0; i < d; i++)
326  operator()(i,j) = c[i];
327 }
328 template<class T>
329 inline
330 void
331 tensor_basic<T>::set_row (const point_basic<T>& r, size_t i, size_t d)
332 {
333  for (size_t j = 0; j < d; j++)
334  operator()(i,j) = r[j];
335 }
336 template<class T>
337 inline
340 {
341  tensor_basic<T> I;
342  for (size_t i = 0; i < d; i++)
343  I(i,i) = 1;
344  return I;
345 }
346 template <class T>
347 inline
348 T
350 {
351  return ddot(a,a);
352 }
353 template <class T>
354 inline
355 T
357 {
358  return norm2(a-b);
359 }
360 template <class U>
361 inline
362 U
364 {
365  return sqrt(norm2(a));
366 }
367 template <class U>
368 inline
369 U
371 {
372  return norm(a-b);
373 }
374 
375 }// namespace rheolef
376 # endif /* _RHEOLEF_TENSOR_H */
T determinant(const tensor_basic< T > &A, size_t d)
Definition: tensor.h:286
tensor_basic< Float > tensor
Definition: tensor.h:85
void fill(const T &init_val)
Definition: tensor.h:141
tensor_basic< U > otimes(const point_basic< U > &u, const point_basic< U > &v, size_t d=3)
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition: tensor.cc:348
tensor_basic< T > & operator=(const tensor_basic< T > &a)
Definition: tensor.h:199
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
Definition: tensor.h:331
U tr(const tensor_basic< U > &a, size_t d=3)
std::ostream & put(std::ostream &s, size_type d=3) const
Definition: tensor.cc:11
tensor_basic(const T &init_val=0)
Definition: tensor.h:155
point - vertex of a mesh
Definition: point.h:22
point_basic< T > row(size_type i) const
Definition: tensor.cc:244
bool operator==(const tensor_basic< T > &) const
Definition: tensor.cc:70
point_basic< T > trans_mult(const point_basic< T > &x) const
Definition: tensor.h:263
irheostream, orheostream - large data streams
Definition: compiler.h:7
static tensor_basic< T > eye(size_type d=3)
Definition: tensor.h:339
T norm2(const tensor_basic< T > &a)
Definition: tensor.h:349
tensor_basic< T > operator/(const T &k) const
Definition: tensor.h:256
size_t ncol() const
Definition: tensor.h:224
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:33
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition: dia.h:118
#define check_macro(ok_condition, message)
Definition: compiler.h:104
size_t d
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
Definition: tensor.cc:211
tensor_basic< T > operator*(const tensor_basic< T > &b) const
Definition: tensor.cc:204
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:147
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition: tensor.cc:328
void set_column(const point_basic< T > &c, size_t j, size_t d=3)
Definition: tensor.h:323
std::istream & get(std::istream &)
Definition: tensor.cc:25
Definition: rotating-hill.h:1
T & operator()(size_type i, size_type j)
Definition: tensor.h:231
point_basic< T > col(size_type i) const
Definition: tensor.cc:254
size_t nrow() const
Definition: tensor.h:217
T tr(const tensor_basic< T > &a, size_t d)
Definition: tensor.h:315
bool operator!=(const tensor_basic< T > &b) const
Definition: tensor.h:62
U dist(const tensor_basic< U > &a, const tensor_basic< U > &b)
Definition: tensor.h:370
const tensor_basic< T > & operator+() const
Definition: tensor.h:63
void cumul_otimes(tensor_basic< T > &t, const point_basic< T > &a, const point_basic< T > &b, size_t na=3)
Definition: tensor.h:270
solver_basic< Float > eye()
Definition: solver.h:162
T dist2(const tensor_basic< T > &a, const tensor_basic< T > &b)
Definition: tensor.h:356
bool invert_3x3(const ublas::matrix< T > &A, ublas::matrix< T > &result)
helper for point_basic<T> & tensor_basic<T>: get basic T type
Definition: point.h:234
return a operator*(xh)
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition: tensor.cc:153
U ddot(const tensor_basic< U > &, const tensor_basic< U > &)
std::istream & operator>>(std::istream &is, const catchmark &m)
Definition: catchmark.h:48
rheolef::std type
tensor_basic< T > otimes(const point_basic< T > &u, const point_basic< T > &v, size_t d)
Definition: tensor.h:277
T determinant(size_type d=3) const
Definition: tensor.cc:221
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:190
U norm(const tensor_basic< U > &a)
Definition: tensor.h:363
point_basic< T > diag(const tensor_basic< T > &a)
Definition: tensor.h:304
csr< T, sequential > trans(const csr< T, sequential > &a)
Definition: csr.h:381
float_traits< T >::type type
Definition: tensor.h:135
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:236
tensor_basic< T > operator-() const
Definition: tensor.cc:79