Rheolef  7.2
an efficient C++ finite element environment
tensor.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_TENSOR_H
2 # define _RHEOLEF_TENSOR_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // AUTHOR: Pierre.Saramito@imag.fr
24 // DATE: 9 october 2003
25 
26 namespace rheolef {
83 } // namespace rheolef
84 
85 #include "rheolef/point.h"
86 namespace rheolef {
87 
88 // [verbatim_tensor_basic]
89 template<class T>
90 class tensor_basic {
91 public:
92 
93  typedef size_t size_type;
94  typedef T element_type;
95  typedef T float_type;
96 
97 // allocators:
98 
99  tensor_basic (const T& init_val = 0);
100  tensor_basic (T x[3][3]);
103 
104  tensor_basic (const std::initializer_list<std::initializer_list<T> >& il);
105 
106 // affectation:
107 
110 
111 // modifiers:
112 
113  void fill (const T& init_val);
114  void reset ();
115  void set_row (const point_basic<T>& r, size_t i, size_t d = 3);
116  void set_column (const point_basic<T>& c, size_t j, size_t d = 3);
117 
118 // accessors:
119 
121  const T& operator()(size_type i, size_type j) const;
122  point_basic<T> row(size_type i) const;
123  point_basic<T> col(size_type i) const;
124  size_t nrow() const; // = 3, for template matrix compatibility
125  size_t ncol() const;
126 
127 // inputs/outputs:
128 
129  std::ostream& put (std::ostream& s, size_type d = 3) const;
130  std::istream& get (std::istream&);
131 
132 // algebra:
133 
134  bool operator== (const tensor_basic<T>&) const;
135  bool operator!= (const tensor_basic<T>& b) const { return ! operator== (b); }
136  const tensor_basic<T>& operator+ () const { return *this; }
137  tensor_basic<T> operator- () const;
141  tensor_basic<T> operator* (const T& k) const;
142  tensor_basic<T> operator/ (const T& k) const;
147  tensor_basic<T>& operator*= (const T& k);
148  tensor_basic<T>& operator/= (const T& k);
149 
150  T determinant (size_type d = 3) const;
151  bool is_symmetric (size_type d = 3) const;
152 
153 // eigenvalues & eigenvectors:
154 
155  // a = q*d*q^T
156  // a may be symmetric
157  // where q=(q1,q2,q3) are eigenvectors in rows (othonormal matrix)
158  // and d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
159  // return d
160  point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const;
161  point_basic<T> eig (size_t dim = 3) const;
162 
163 // singular value decomposition:
164 
165  // a = u*s*v^T
166  // a can be unsymmetric
167  // where u=(u1,u2,u3) are left pseudo-eigenvectors in rows (othonormal matrix)
168  // v=(v1,v2,v3) are right pseudo-eigenvectors in rows (othonormal matrix)
169  // and s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
170  // return s
171  point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const;
172 // [verbatim_tensor_basic]
173 
174 // data:
175  T _x[3][3];
176 // [verbatim_tensor_basic_cont]
177 };
178 // [verbatim_tensor_basic_cont]
179 
180 // [verbatim_tensor]
182 // [verbatim_tensor]
183 
184 // [verbatim_tensor_basic_cont2]
185 template <class U>
187 
188 template <class U>
189 tensor_basic<U> trans (const tensor_basic<U>& a, size_t d = 3);
190 
191 template <class U>
192 void prod (const tensor_basic<U>& a, const tensor_basic<U>& b, tensor_basic<U>& result,
193  size_t di=3, size_t dj=3, size_t dk=3);
194 
195 // tr(a) = a00 + a11 + a22
196 template <class U>
197 U tr (const tensor_basic<U>& a, size_t d=3);
198 
199 template <class U>
201 
202 // a = u otimes v <==> aij = ui*vj
203 template <class U>
204 tensor_basic<U> otimes (const point_basic<U>& u, const point_basic<U>& v, size_t d=3);
205 
206 template <class U>
207 tensor_basic<U> inv (const tensor_basic<U>& a, size_t d = 3);
208 
209 template <class U>
211 template <class U>
213 
214 template <class U>
215 U determinant (const tensor_basic<U>& A, size_t d = 3);
216 
217 template <class U>
219 
220 // matrix exponential:
221 template<class T>
222 tensor_basic<T> exp (const tensor_basic<T>& a, size_t d = 3);
223 
224 // inputs/outputs:
225 template<class T>
226 inline
227 std::istream& operator>> (std::istream& in, tensor_basic<T>& a) {
228  return a.get (in); }
229 
230 template<class T>
231 inline
232 std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a) {
233  return a.put (out); }
234 
235 // t += a otimes b
236 template<class T>
237 void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
238 
239 template<class T>
240 void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb);
241 // [verbatim_tensor_basic_cont2]
242 
243 // -----------------------------------------------------------------------
244 // inlined
245 // -----------------------------------------------------------------------
246 template<class T> struct float_traits<tensor_basic<T> > { typedef typename float_traits<T>::type type; };
247 template<class T> struct scalar_traits<tensor_basic<T> > { typedef T type; };
248 
249 template<class T>
250 inline
251 void
252 tensor_basic<T>::fill (const T& init_val)
253 {
254  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
255  _x[i][j] = init_val;
256 }
257 template<class T>
258 inline
259 void
261 {
262  fill (0);
263 }
264 template<class T>
265 inline
267 {
268  fill (init_val);
269 }
270 template<class T>
271 inline
273 {
274  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
275  _x[i][j] = x[i][j];
276 }
277 template<class T>
278 inline
280 {
281  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
282  _x[i][j] = a._x[i][j];
283 }
284 template<class T>
285 tensor_basic<T>::tensor_basic (const std::initializer_list<std::initializer_list<T> >& il) : _x() {
286  typedef typename std::initializer_list<std::initializer_list<T> >::const_iterator const_iterator;
287  typedef typename std::initializer_list<T>::const_iterator const_iterator_row;
288  fill (T());
289  check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
290  size_type i = 0;
291  for (const_iterator iter = il.begin(); iter != il.end(); ++iter, ++i) {
292  const std::initializer_list<T>& row = *iter;
293  check_macro (row.size() <= 3, "unexpected initializer list size=" << row.size() << " > 3");
294  size_type j = 0;
295  for (const_iterator_row jter = row.begin(); jter != row.end(); ++jter, ++j) {
296  _x[i][j] = *jter;
297  }
298  }
299 }
300 template<class T>
301 inline
304 {
305  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
306  _x[i][j] = a._x[i][j];
307  return *this;
308 }
309 template<class T>
310 inline
313 {
314  for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
315  _x[i][j] = val;
316  return *this;
317 }
318 template<class T>
319 inline
320 size_t
322 {
323  return 3;
324 }
325 template<class T>
326 inline
327 size_t
329 {
330  return 3;
331 }
332 template<class T>
333 inline
334 T&
336 {
337  return _x[i%3][j%3];
338 }
339 template<class T>
340 inline
341 const T&
343 {
344  return _x[i%3][j%3];
345 }
346 template <class T, class U>
347 inline
348 typename
349 std::enable_if<
350  details::is_rheolef_arithmetic<U>::value
352 >::type
353 operator* (const U& k, const tensor_basic<T>& a)
354 {
355  return a*k;
356 }
357 template<class T>
358 inline
359 tensor_basic<T>
361 {
362  return operator* (1./k);
363 }
364 template<class T>
365 inline
368 {
369  return x*(*this);
370 }
371 template<class T>
372 inline
373 void
375 {
376  cumul_otimes (t, a, b, n, n);
377 }
378 template<class T>
379 inline
380 tensor_basic<T>
381 otimes (const point_basic<T>& u, const point_basic<T>& v, size_t d)
382 {
384  cumul_otimes (a, u, v, d, d);
385  return a;
386 }
387 template<class T>
388 inline
389 T
390 determinant (const tensor_basic<T>& A, size_t d)
391 {
392  return A.determinant (d);
393 }
394 template<class T>
395 inline
396 tensor_basic<T>
398 {
400  a(0,0) = d[0];
401  a(1,1) = d[1];
402  a(2,2) = d[2];
403  return a;
404 }
405 template<class T>
406 inline
409 {
411  d[0] = a(0,0);
412  d[1] = a(1,1);
413  d[2] = a(2,2);
414  return d;
415 }
416 template <class T>
417 inline
418 T
419 tr (const tensor_basic<T>& a, size_t d) {
420  T sum = 0;
421  for (size_t i = 0; i < d; i++) sum += a(i,i);
422  return sum;
423 }
424 template<class T>
425 inline
426 void
427 tensor_basic<T>::set_column (const point_basic<T>& c, size_t j, size_t d)
428 {
429  for (size_t i = 0; i < d; i++)
430  operator()(i,j) = c[i];
431 }
432 template<class T>
433 inline
434 void
435 tensor_basic<T>::set_row (const point_basic<T>& r, size_t i, size_t d)
436 {
437  for (size_t j = 0; j < d; j++)
438  operator()(i,j) = r[j];
439 }
440 template<class T>
441 inline
444 {
445  tensor_basic<T> I;
446  for (size_t i = 0; i < d; i++)
447  I(i,i) = 1;
448  return I;
449 }
450 template <class T>
451 inline
452 T
454 {
455  return ddot(a,a);
456 }
457 template <class T>
458 inline
459 T
461 {
462  return norm2(a-b);
463 }
464 template <class U>
465 inline
466 U
468 {
469  return sqrt(norm2(a));
470 }
471 template <class U>
472 inline
473 U
475 {
476  return norm(a-b);
477 }
478 
479 }// namespace rheolef
480 // -------------------------------------------------------------
481 // serialization
482 // -------------------------------------------------------------
483 #ifdef _RHEOLEF_HAVE_MPI
484 namespace boost {
485  namespace serialization {
486  template<class Archive, class T>
487  void serialize (Archive& ar, class rheolef::tensor_basic<T>& x, const unsigned int version) {
488  ar & x(0,0);
489  ar & x(0,1);
490  ar & x(0,2);
491  ar & x(1,0);
492  ar & x(1,1);
493  ar & x(1,2);
494  ar & x(2,0);
495  ar & x(2,1);
496  ar & x(2,2);
497  }
498  } // namespace serialization
499 } // namespace boost
500 
501 // Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
502 // When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
503 // To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
504 namespace boost {
505  namespace mpi {
506  // TODO: when tensor_basic<T> is not a simple type, such as T=bigfloat or T=gmp, etc
507  template <>
508  struct is_mpi_datatype<rheolef::tensor_basic<double> > : mpl::true_ { };
509  } // namespace mpi
510 } // namespace boost
511 #endif // _RHEOLEF_HAVE_MPI
512 #endif /* _RHEOLEF_TENSOR_H */
std::ostream & put(std::ostream &s, size_type d=3) const
Definition: tensor.cc:37
T & operator()(size_type i, size_type j)
Definition: tensor.h:335
tensor_basic< T > & operator=(const tensor_basic< T > &a)
Definition: tensor.h:303
std::istream & get(std::istream &)
Definition: tensor.cc:51
point_basic< T > trans_mult(const point_basic< T > &x) const
Definition: tensor.h:367
tensor_basic< T > & operator+=(const tensor_basic< T > &)
Definition: tensor.cc:174
static tensor_basic< T > eye(size_type d=3)
Definition: tensor.h:443
point_basic< T > col(size_type i) const
Definition: tensor.cc:323
T determinant(size_type d=3) const
Definition: tensor.cc:288
tensor_basic< T > & operator*=(const T &k)
Definition: tensor.cc:192
tensor_basic(const std::initializer_list< std::initializer_list< T > > &il)
Definition: tensor.h:285
const T & operator()(size_type i, size_type j) const
Definition: tensor.h:342
tensor_basic< T > operator*(const tensor_basic< T > &b) const
Definition: tensor.cc:270
tensor_basic(const tensor_basic< T > &a)
Definition: tensor.h:279
bool is_symmetric(size_type d=3) const
Definition: tensor.cc:351
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition: tensor.cc:426
tensor_basic< T > operator-() const
Definition: tensor.cc:110
tensor_basic(T x[3][3])
Definition: tensor.h:272
tensor_basic< T > operator/(const T &k) const
Definition: tensor.h:360
size_t ncol() const
Definition: tensor.h:328
void fill(const T &init_val)
Definition: tensor.h:252
const tensor_basic< T > & operator+() const
Definition: tensor.h:136
bool operator!=(const tensor_basic< T > &b) const
Definition: tensor.h:135
size_t nrow() const
Definition: tensor.h:321
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
Definition: tensor.h:435
tensor_basic< T > & operator-=(const tensor_basic< T > &)
Definition: tensor.cc:183
bool operator==(const tensor_basic< T > &) const
Definition: tensor.cc:101
tensor_basic< T > & operator/=(const T &k)
Definition: tensor.cc:201
point_basic< T > row(size_type i) const
Definition: tensor.cc:313
tensor_basic(const T &init_val=0)
Definition: tensor.h:266
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition: tensor.cc:470
void set_column(const point_basic< T > &c, size_t j, size_t d=3)
Definition: tensor.h:427
rheolef::space_base_rep< T, M > t
rheolef::std type
tensor_basic< Float > tensor
Definition: tensor.h:181
point_basic< T >
Definition: piola_fem.h:135
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void serialize(Archive &ar, class rheolef::point_basic< T > &x, const unsigned int version)
Definition: point.h:579
This file is part of Rheolef.
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:298
std::istream & operator>>(std::istream &is, const catchmark &m)
Definition: catchmark.h:88
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition: tensor.cc:219
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
U determinant(const tensor_basic< U > &A, size_t d=3)
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:292
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
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:56
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
U tr(const tensor_basic< U > &a, size_t d=3)
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition: tensor.cc:333
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
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
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
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:92
tensor_basic< U > otimes(const point_basic< U > &u, const point_basic< U > &v, size_t d=3)
float_traits< T >::type type
Definition: tensor.h:246
helper for std::complex<T>: get basic T type
Definition: Float.h:93
helper for point_basic<T> & tensor_basic<T>: get basic T type
Definition: point.h:323
Definition: leveque.h:25