Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
26namespace rheolef {
83} // namespace rheolef
84
85#include "rheolef/point.h"
86namespace rheolef {
87
88// [verbatim_tensor_basic]
89template<class T>
91public:
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;
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; }
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]
185template <class U>
187
188template <class U>
189tensor_basic<U> trans (const tensor_basic<U>& a, size_t d = 3);
190
191template <class U>
192void 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
196template <class U>
197U tr (const tensor_basic<U>& a, size_t d=3);
198
199template <class U>
201
202// a = u otimes v <==> aij = ui*vj
203template <class U>
205
206template <class U>
207tensor_basic<U> inv (const tensor_basic<U>& a, size_t d = 3);
208
209template <class U>
211template <class U>
213
214template <class U>
215U determinant (const tensor_basic<U>& A, size_t d = 3);
216
217template <class U>
219
220// matrix exponential:
221template<class T>
222tensor_basic<T> exp (const tensor_basic<T>& a, size_t d = 3);
223
224// inputs/outputs:
225template<class T>
226inline
227std::istream& operator>> (std::istream& in, tensor_basic<T>& a) {
228 return a.get (in); }
229
230template<class T>
231inline
232std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a) {
233 return a.put (out); }
234
235// t += a otimes b
236template<class T>
237void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
238
239template<class T>
240void 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// -----------------------------------------------------------------------
246template<class T> struct float_traits<tensor_basic<T> > { typedef typename float_traits<T>::type type; };
247template<class T> struct scalar_traits<tensor_basic<T> > { typedef T type; };
248
249template<class T>
250inline
251void
252tensor_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}
257template<class T>
258inline
259void
261{
262 fill (0);
263}
264template<class T>
265inline
267{
268 fill (init_val);
269}
270template<class T>
271inline
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}
277template<class T>
278inline
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}
284template<class T>
285tensor_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}
300template<class T>
301inline
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}
309template<class T>
310inline
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}
318template<class T>
319inline
320size_t
322{
323 return 3;
324}
325template<class T>
326inline
327size_t
329{
330 return 3;
331}
332template<class T>
333inline
334T&
336{
337 return _x[i%3][j%3];
338}
339template<class T>
340inline
341const T&
343{
344 return _x[i%3][j%3];
345}
346template <class T, class U>
347inline
348typename
349std::enable_if<
350 details::is_rheolef_arithmetic<U>::value
352>::type
353operator* (const U& k, const tensor_basic<T>& a)
354{
355 return a*k;
356}
357template<class T>
358inline
359tensor_basic<T>
361{
362 return operator* (1./k);
363}
364template<class T>
365inline
368{
369 return x*(*this);
370}
371template<class T>
372inline
373void
375{
376 cumul_otimes (t, a, b, n, n);
377}
378template<class T>
379inline
380tensor_basic<T>
381otimes (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}
387template<class T>
388inline
389T
391{
392 return A.determinant (d);
393}
394template<class T>
395inline
396tensor_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}
405template<class T>
406inline
407point_basic<T>
409{
411 d[0] = a(0,0);
412 d[1] = a(1,1);
413 d[2] = a(2,2);
414 return d;
415}
416template <class T>
417inline
418T
419tr (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}
424template<class T>
425inline
426void
427tensor_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}
432template<class T>
433inline
434void
435tensor_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}
440template<class T>
441inline
444{
446 for (size_t i = 0; i < d; i++)
447 I(i,i) = 1;
448 return I;
449}
450template <class T>
451inline
452T
454{
455 return ddot(a,a);
456}
457template <class T>
458inline
459T
461{
462 return norm2(a-b);
463}
464template <class U>
465inline
466U
468{
469 return sqrt(norm2(a));
470}
471template <class U>
472inline
473U
475{
476 return norm(a-b);
477}
478
479}// namespace rheolef
480// -------------------------------------------------------------
481// serialization
482// -------------------------------------------------------------
483#ifdef _RHEOLEF_HAVE_MPI
484namespace 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.:
504namespace 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
const tensor_basic< T > & operator+() const
Definition tensor.h:136
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
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
tensor_basic< Float > tensor
Definition tensor.h:181
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.
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition tensor.cc:219
tensor_basic< U > otimes(const point_basic< U > &u, const point_basic< U > &v, size_t d=3)
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
U tr(const tensor_basic< U > &a, size_t d=3)
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition tensor-exp.cc:92
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
U determinant(const tensor_basic< U > &A, size_t d=3)
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
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
std::istream & operator>>(std::istream &is, const catchmark &m)
Definition catchmark.h:88
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56
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, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
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