Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
tensor-exp.cc
Go to the documentation of this file.
1
21// exp(tensor) : 2D is explicit while 3D uses pade' approx
22//
23#include "rheolef/tensor.h"
24#include "rheolef/compiler_eigen.h"
25
26#undef _RHEOLEF_EIGEN_HAVE_EXPM // not yet supported officially, but faster than via eig(a)
27#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
28#include <unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h>
29#endif // _RHEOLEF_EIGEN_HAVE_EXPM
30
31namespace rheolef {
32// =========================================================================
33// 2D case: explicit
34// =========================================================================
35//
36// see tensor_exp_tst.mac
37// see also maxima: matrixexp.usage allows explicit expm() computation in 2D
38// http://www.ma.utexas.edu/pipermail/maxima/2006/003031.html
39// /usr/share/maxima/5.27.0/share/linearalgebra/matrixexp.usage
40// refs. KanGueFor-2009, p. 50 (is buggy: missing 2 factor in A00 & A11...)
41
42template <class T>
43static
44tensor_basic<T>
45exp2d (const tensor_basic<T>& chi) {
46 static T eps = 1e3*std::numeric_limits<T>::epsilon();
47 T a = chi(0,0), b = chi(1,0), c = chi(1,1);
48 T b2=sqr(b);
49 T d2 = sqr(a-c)+4*b2;
50 T d = sqrt(d2);
51 tensor_basic<T> A;
52 if (fabs(d) < eps) { // chi = a*I
53 A(0,0) = A(1,1) = exp(a);
54 A(0,1) = A(1,0) = 0.;
55 return A;
56 }
57 T a2=sqr(a), c2=sqr(c);
58 T ed = exp(d);
59 T k = exp((a+c-d)/2)/(2*d);
60 T x1 = (ed+1)*d;
61 T x2 = (ed-1)*(a-c);
62 T x3 = 2*(ed-1)*b;
63 A(0,0) = k*(x1 + x2);
64 A(1,1) = k*(x1 - x2);
65 A(1,0) = A(0,1) = k*x3;
66 return A;
67}
68template<typename T, int N>
69Eigen::Matrix<T,N,N>
70expm_eig (const Eigen::Matrix<T,N,N>& a) {
71 using namespace Eigen;
72 // eig of sym matrix, then exp
73 SelfAdjointEigenSolver<Matrix<T,N,N> > eig (a);
74 Matrix<T,N,1> e;
75 for (size_t i = 0; i < N; ++i) {
76 e(i) = exp(eig.eigenvalues()(i));
77 }
78 return eig.eigenvectors()*e.asDiagonal()*eig.eigenvectors().transpose();
79}
80#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
81template<typename T, int N>
82Eigen::Matrix<T,N,N>
83expm_pade (const Eigen::Matrix<T,N,N>& a) {
84 return a.exp(); // in "eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h"
85}
86#endif // _RHEOLEF_EIGEN_HAVE_EXPM
87// =========================================================================
88// tensor interface
89// =========================================================================
90template <class T>
91tensor_basic<T>
92exp (const tensor_basic<T>& a, size_t d) {
93 // check that a is symmetric ; otherwise, the result would be complex
94 check_macro (a.is_symmetric(d), "exp: tensor should be symmetric");
95 using namespace std;
96 using namespace Eigen;
97 if (d == 1) return tensor(exp(a(0,0)));
98 if (d == 2) return exp2d (a);
99 // 3D case:
100 Matrix<T,3,3> a1 (d,d);
101 for (size_t i = 0; i < d; ++i)
102 for (size_t j = 0; j < d; ++j) a1(i,j) = a(i,j);
103#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
104 Matrix<T,3,3> b1 = expm_pade(a1);
105#else // _RHEOLEF_EIGEN_HAVE_EXPM
106 Matrix<T,3,3> b1 = expm_eig (a1);
107#endif // _RHEOLEF_EIGEN_HAVE_EXPM
108 tensor b;
109 for (size_t i = 0; i < d; ++i)
110 for (size_t j = 0; j < d; ++j) b(i,j) = b1(i,j);
111 return b;
112}
113// ----------------------------------------------------------------------------
114// instanciation in library
115// ----------------------------------------------------------------------------
116#define _RHEOLEF_instanciation(T) \
117template tensor_basic<T> exp (const tensor_basic<T>& a, size_t d); \
118
120
121}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
see the tensor page for the full documentation
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 > exp(const tensor_basic< T > &a, size_t d)
Definition tensor-exp.cc:92
Eigen::Matrix< T, N, N > expm_eig(const Eigen::Matrix< T, N, N > &a)
Definition tensor-exp.cc:70
STL namespace.
size_t N