Rheolef  7.1
an efficient C++ finite element environment
tensor-exp.cc
Go to the documentation of this file.
1 // 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 
31 namespace 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 
42 template <class T>
43 static
44 tensor_basic<T>
45 exp2d (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 }
68 template<typename T, int N>
69 Eigen::Matrix<T,N,N>
70 expm_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
81 template<typename T, int N>
82 Eigen::Matrix<T,N,N>
83 expm_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 // =========================================================================
90 template <class T>
91 tensor_basic<T>
92 exp (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) \
117 template tensor_basic<T> exp (const tensor_basic<T>& a, size_t d); \
118 
120 
121 }// namespace rheolef
rheolef::expm_eig
Eigen::Matrix< T, N, N > expm_eig(const Eigen::Matrix< T, N, N > &a)
Definition: tensor-exp.cc:70
Eigen
Definition: field_eigen.h:35
tensor
see the tensor page for the full documentation
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_contraction.c2
c2
Definition: mkgeo_contraction.sh:199
N
size_t N
Definition: neumann-laplace-check.cc:24
chi
Definition: transport_tensor_exact.icc:65
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::tensor_basic
Definition: tensor.h:90
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::tensor
tensor_basic< Float > tensor
Definition: tensor.h:181
a
Definition: diffusion_isotropic.h:25
rheolef::exp
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:92
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
Float
see the Float page for the full documentation
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::std
Definition: vec_expr_v2.h:402
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
T
Expr1::float_type T
Definition: field_expr.h:261