Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
tensor.cc
Go to the documentation of this file.
1
21//
22// 2-tensor
23//
24// author: Pierre.Saramito@imag.fr
25//
26#include "rheolef/tensor.h"
27#include "rheolef/compiler_eigen.h"
28
29using namespace rheolef;
30using namespace std;
31namespace rheolef {
32
33
34// output
35template<class T>
36ostream&
37tensor_basic<T>::put (ostream& out, size_type d) const
38{
39 switch (d) {
40 case 0 : return out;
41 case 1 : return out << _x[0][0];
42 case 2 : return out << "[" << _x[0][0] << ", " << _x[0][1] << ";\n"
43 << " " << _x[1][0] << ", " << _x[1][1] << "]";
44 default: return out << "[" << _x[0][0] << ", " << _x[0][1] << ", " << _x[0][2] << ";\n"
45 << " " << _x[1][0] << ", " << _x[1][1] << ", " << _x[1][2] << ";\n"
46 << " " << _x[2][0] << ", " << _x[2][1] << ", " << _x[2][2] << "]";
47 }
48}
49template<class T>
50istream&
52{
53 for (size_type i = 0; i < 3; i++)
54 for (size_type j = 0; j < 3; j++)
55 _x[i][j] = 0.;
56 char c;
57 if (!in.good()) return in;
58 in >> ws >> c;
59 // no tensor specifier:
60 // 1d case: "3.14"
61 if (c != '[') {
62 in.unget();
63 return in >> _x[0][0];
64 }
65 // has a tensor specifier, as:
66 // [ 23, 17; 25, 2]
67 in >> _x[0][0] >> ws >> c;
68 if (c == ']') return in; // 1d case
69 if (c != ',') error_macro ("invalid tensor input: expect `,'");
70 in >> _x[0][1];
71 in >> ws >> c;
72 size_type d = 3;
73 if (c == ',') d = 3;
74 else if (c == ';') d = 2;
75 else error_macro ("invalid tensor input: expect `,' or `;'");
76 if (d == 3) {
77 in >> _x[0][2] >> ws >> c;
78 if (c != ';') error_macro ("invalid tensor input: expect `;'");
79 }
80 in >> _x[1][0] >> ws >> c;
81 if (c != ',') error_macro ("invalid tensor input: expect `,'");
82 in >> _x[1][1] >> ws >> c;
83 if (d == 2) {
84 if (c != ']') error_macro ("invalid tensor input: expect `]'");
85 return in;
86 }
87 // d == 3
88 if (c != ',') error_macro ("invalid tensor input: expect `,'");
89 in >> _x[1][2] >> ws >> c;
90 if (c != ';') error_macro ("invalid tensor input: expect `;'");
91 in >> _x[2][0] >> ws >> c;
92 if (c != ',') error_macro ("invalid tensor input: expect `,'");
93 in >> _x[2][1] >> ws >> c;
94 if (c != ',') error_macro ("invalid tensor input: expect `,'");
95 in >> _x[2][2] >> ws >> c;
96 if (c != ']') error_macro ("invalid tensor input: expect `]'");
97 return in;
98}
99template<class T>
100bool
102{
103 for (size_type i = 0; i < 3; i++)
104 for (size_type j = 0; j < 3; j++)
105 if (_x[i][j] != b._x[i][j]) return false;
106 return true;
107}
108template<class T>
111{
113 typedef typename tensor_basic<T>::size_type size_type;
114 for (size_type i = 0; i < 3; i++)
115 for (size_type j = 0; j < 3; j++)
116 b._x[i][j] = - _x[i][j];
117 return b;
118}
119template<class T>
124 typedef typename tensor_basic<T>::size_type size_type;
125 for (size_type i = 0; i < 3; i++)
126 for (size_type j = 0; j < 3; j++)
127 c._x[i][j] = _x[i][j] + b._x[i][j];
128 return c;
130template<class T>
133{
135 typedef typename tensor_basic<T>::size_type size_type;
136 for (size_type i = 0; i < 3; i++)
137 for (size_type j = 0; j < 3; j++)
138 c._x[i][j] = _x[i][j] - b._x[i][j];
139 return c;
141template<class T>
144{
146 for (size_type i = 0; i < 3; i++)
147 for (size_type j = 0; j < 3; j++)
148 b._x[i][j] = _x[i][j] * k;
149 return b;
151template<class T>
154{
156 typedef typename tensor_basic<T>::size_type size_type;
157 for (size_type i = 0; i < 3; i++)
158 for (size_type j = 0; j < 3; j++)
159 y [i] += _x[i][j] * x[j];
160 return y;
162template<class T>
164{
166 typedef typename tensor_basic<T>::size_type size_type;
167 for (size_type i = 0; i < 3; i++)
168 for (size_type j = 0; j < 3; j++)
169 y [i] += a._x[j][i] * x[j];
170 return y;
172template<class T>
175{
176 for (size_type i = 0; i < 3; i++)
177 for (size_type j = 0; j < 3; j++)
178 _x[i][j] += b._x[i][j];
179 return *this;
180}
181template<class T>
184{
185 for (size_type i = 0; i < 3; i++)
186 for (size_type j = 0; j < 3; j++)
187 _x[i][j] -= b._x[i][j];
188 return *this;
189}
190template<class T>
193{
194 for (size_type i = 0; i < 3; i++)
195 for (size_type j = 0; j < 3; j++)
196 _x[i][j] *= k;
197 return *this;
198}
199template<class T>
202{
203 for (size_type i = 0; i < 3; i++)
204 for (size_type j = 0; j < 3; j++)
205 _x[i][j] /= k;
206 return *this;
207}
208template<class T>
210{
212 typedef typename tensor_basic<T>::size_type size_type;
213 for (size_type i = 0; i < d; i++)
214 for (size_type j = 0; j < d; j++)
215 b._x[i][j] = a._x[j][i];
216 return b;
217}
218template<class T>
220{
222 typedef typename tensor_basic<T>::size_type size_type;
223 switch (d) {
224 case 0:
225 break;
226 case 1:
227 b(0,0) = 1/a(0,0);
228 break;
229 case 2: {
230 T det = a.determinant(2);
231 b(0,0) = a(1,1)/det;
232 b(1,1) = a(0,0)/det;
233 b(0,1) = - a(0,1)/det;
234 b(1,0) = - a(1,0)/det;
235 break;
236 }
237 case 3:
238 default: {
239 T det = a.determinant(3);
240 b(0,0) = (a(1,1)*a(2,2) - a(1,2)*a(2,1))/det;
241 b(0,1) = (a(0,2)*a(2,1) - a(0,1)*a(2,2))/det;
242 b(0,2) = (a(0,1)*a(1,2) - a(0,2)*a(1,1))/det;
243 b(1,0) = (a(1,2)*a(2,0) - a(1,0)*a(2,2))/det;
244 b(1,1) = (a(0,0)*a(2,2) - a(0,2)*a(2,0))/det;
245 b(1,2) = (a(0,2)*a(1,0) - a(0,0)*a(1,2))/det;
246 b(2,0) = (a(1,0)*a(2,1) - a(1,1)*a(2,0))/det;
247 b(2,1) = (a(0,1)*a(2,0) - a(0,0)*a(2,1))/det;
248 b(2,2) = (a(0,0)*a(1,1) - a(0,1)*a(1,0))/det;
249 break;
250 }
251 }
252 return b;
253}
254template<class T>
255void
257 size_t di, size_t dj, size_t dk)
258{
259 typedef typename tensor_basic<T>::size_type size_type;
260 for (size_type i = 0; i < di; i++)
261 for (size_type j = 0; j < dj; j++) {
262 T sum = 0;
263 for (size_type k = 0; k < dk; k++)
264 sum += a(i,k) * b(k,j);
265 result(i,j) = sum;
266 }
267}
268template<class T>
271{
273 prod (*this,b,c,3,3,3);
274 return c;
275}
277template<class T>
279{
280 T r = 0;
281 typedef typename tensor_basic<T>::size_type size_type;
282 for (size_type i = 0; i < 3; i++)
283 for (size_type j = 0; j < 3; j++)
284 r += a._x[i][j] * b._x[i][j];
285 return r;
286}
287template<class T>
289{
290 switch (d) {
291 case 0 : return 1;
292 case 1 : return _x[0][0];
293 case 2 : return _x[0][0]*_x[1][1] - _x[0][1]*_x[1][0];
294 case 3 : return _x[0][0]*(_x[1][1]*_x[2][2] - _x[1][2]*_x[2][1])
295 - _x[0][1]*(_x[1][0]*_x[2][2] - _x[1][2]*_x[2][0])
296 + _x[0][2]*(_x[1][0]*_x[2][1] - _x[1][1]*_x[2][0]);
297 default: {
298 error_macro ("determinant: unexpected dimension " << d);
299 }
300 }
301}
302// t += a otimes b
303template<class T>
304void
305cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb)
306{
307 for (size_t i = 0; i < na; i++)
308 for (size_t j = 0; j < nb; j++)
309 t(i,j) += a[i] * b[j];
310}
311template<class T>
314{
316 r[0] = _x[i][0];
317 r[1] = _x[i][1];
318 r[2] = _x[i][2];
319 return r;
320}
321template<class T>
324{
326 c[0] = _x[0][j];
327 c[1] = _x[1][j];
328 c[2] = _x[2][j];
329 return c;
330}
331template<class T>
332bool
334{
335 T det = determinant(A,3);
336 if (1+det == 1) return false;
337 T invdet = 1.0/det;
338 result(0,0) = (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
339 result(0,1) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
340 result(0,2) = (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
341 result(1,0) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
342 result(1,1) = (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
343 result(1,2) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
344 result(2,0) = (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
345 result(2,1) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
346 result(2,2) = (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
347 return true;
348}
349template<class T>
350bool
352 if (d <= 1) return true;
353 static T tol = 1e3*std::numeric_limits<T>::epsilon();
354 if (d <= 2) return fabs (_x[0][1] - _x[1][0]) < tol;
355 return (fabs (_x[0][1] - _x[1][0]) < tol)
356 && (fabs (_x[0][2] - _x[2][0]) < tol)
357 && (fabs (_x[1][2] - _x[2][1]) < tol);
358}
359// ------------------------------------------
360// eigenvalues: the 3D case
361// ------------------------------------------
362template<class T>
363static point_basic<T> eig3x3 (const tensor_basic<T>& a, tensor_basic<T>& q)
364{
365 using namespace Eigen;
366 Matrix<T,3,3> a1;
367 for (size_t i = 0; i < 3; ++i)
368 for (size_t j = 0; j < 3; ++j) a1(i,j) = a(i,j);
369 SelfAdjointEigenSolver<Matrix<T,3,3> > es (a1);
371 for (size_t i = 0; i < 3; ++i) {
372 lambda[i] = es.eigenvalues()(i);
373 for (size_t j = 0; j < 3; ++j) {
374 q(i,j) = es.eigenvectors()(i,j);
375 }
376 }
377 return lambda;
378}
379// ------------------------------------------
380// eigenvalues: the 2D case
381// ------------------------------------------
382template<class T>
383static point_basic<T> eig2x2 (const tensor_basic<T>& a, tensor_basic<T>& q)
384{
385 point_basic<T> d (0,0,0);
386 q.fill (0);
387 if (a(0,1) == 0) {
388 // a is already diagonal
389 if (a(0,0) >= a(1,1)) {
390 d[0] = a(0,0);
391 d[1] = a(1,1);
392 q(0,0) = q(1,1) = 1.;
393 q(0,1) = q(1,0) = 0.;
394 } else { // swap column 0 & 1:
395 d[1] = a(0,0);
396 d[0] = a(1,1);
397 q(0,0) = q(1,1) = 0.;
398 q(0,1) = q(1,0) = 1.;
399 }
400 return d;
401 }
402 // here a(0,1) != 0
403 T discr = sqr(a(1,1) - a(0,0)) + 4*sqr(a(0,1));
404 T trace = a(0,0) + a(1,1);
405 d[0] = (trace + sqrt(discr))/2;
406 d[1] = (trace - sqrt(discr))/2;
407 T lost_precision = 1e-6;
408 if (fabs(d[0]) + lost_precision*fabs(d[1]) == fabs(d[0])) { // d[1] == 0 up to machine precision
409 d[1] = 0.;
410 }
411 T c0 = (d[0]-a(0,0))/a(0,1);
412 T n0 = sqrt(1+sqr(c0));
413 q(0,0) = 1/n0;
414 q(1,0) = c0/n0;
415 T c1 = (d[1]-a(0,0))/a(0,1);
416 T n1 = sqrt(1+sqr(c1));
417 q(0,1) = 1/n1;
418 q(1,1) = c1/n1;
419 return d;
420}
421// ------------------------------------------
422// eigenvalues: general case
423// ------------------------------------------
424template<class T>
427{
428 check_macro (is_symmetric(d), "eig: tensor should be symmetric");
429 switch (d) {
430 case 1 : {
432 q(0,0) = 1; d[0] = operator()(0,0); return d;
433 }
434 case 2 : return eig2x2 (*this, q);
435 default: return eig3x3 (*this, q);
436 }
437}
438template<class T>
441{
443 return eig (q, d);
444}
445// ------------------------------------------
446// svd: the 3D case
447// ------------------------------------------
448template<class T, int N>
449static
451eigen_svd (const tensor_basic<T>& a, tensor_basic<T>& u, tensor_basic<T>& v)
452{
453 using namespace Eigen;
454 Matrix<T,N,N> a1;
455 for (size_t i = 0; i < N; ++i)
456 for (size_t j = 0; j < N; ++j) a1(i,j) = a(i,j);
457 JacobiSVD<Matrix<T,N,N> > svd (a1, ComputeFullU | ComputeFullV);
458 point s;
459 for (size_t i = 0; i < N; ++i) {
460 s[i] = svd.singularValues()(i);
461 for (size_t j = 0; j < N; ++j) {
462 u(i,j) = svd.matrixU()(i,j);
463 v(i,j) = svd.matrixV()(i,j);
464 }
465 }
466 return s;
467}
468template<class T>
471{
472 if (d == 1) {
473 u(0,0) = v(0,0) = 1;
474 return point(operator()(0,0));
475 }
476 if (d == 2) return eigen_svd<T,2> (*this, u, v);
477 else return eigen_svd<T,3> (*this, u, v);
478}
479// ----------------------------------------------------------------------------
480// instanciation in library
481// ----------------------------------------------------------------------------
482#define _RHEOLEF_instanciation(T) \
483template class tensor_basic<T>; \
484template point_basic<T> operator* (const point_basic<T>& x, const tensor_basic<T>& a); \
485template tensor_basic<T> trans (const tensor_basic<T>&, size_t); \
486template tensor_basic<T> inv (const tensor_basic<T>&, size_t); \
487template void prod (const tensor_basic<T>& a, const tensor_basic<T>& b, tensor_basic<T>& result, \
488 size_t di, size_t dj, size_t dk); \
489template T ddot (const tensor_basic<T>& a, const tensor_basic<T> & b); \
490template void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb); \
491template bool invert_3x3 (const tensor_basic<T>& A, tensor_basic<T>& result); \
492
494
495}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the point page for the full documentation
std::ostream & put(std::ostream &s, size_type d=3) const
Definition tensor.cc:37
std::istream & get(std::istream &)
Definition tensor.cc:51
tensor_basic< T > & operator+=(const tensor_basic< T > &)
Definition tensor.cc:174
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< T > operator*(const tensor_basic< T > &b) const
Definition tensor.cc:270
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
void fill(const T &init_val)
Definition tensor.h:252
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
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition tensor.cc:470
#define error_macro(message)
Definition dis_macros.h:49
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)")
This file is part of Rheolef.
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition tensor.cc:219
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 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
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
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
t operator()(const t &a, const t &b)
Definition space.cc:386
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
STL namespace.
size_t N
Definition leveque.h:25