23#include "rheolef/tensor.h"
24#include "rheolef/compiler_eigen.h"
26#undef _RHEOLEF_EIGEN_HAVE_EXPM
27#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
28#include <unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h>
45exp2d (
const tensor_basic<T>&
chi) {
46 static T eps = 1e3*std::numeric_limits<T>::epsilon();
53 A(0,0) =
A(1,1) =
exp(a);
57 T a2=sqr(a),
c2=sqr(c);
59 T k =
exp((a+c-
d)/2)/(2*
d);
65 A(1,0) =
A(0,1) = k*x3;
68template<
typename T,
int N>
71 using namespace Eigen;
73 SelfAdjointEigenSolver<Matrix<T,N,N> > eig (a);
75 for (
size_t i = 0; i <
N; ++i) {
76 e(i) =
exp(eig.eigenvalues()(i));
78 return eig.eigenvectors()*e.asDiagonal()*eig.eigenvectors().transpose();
80#ifdef _RHEOLEF_EIGEN_HAVE_EXPM
81template<
typename T,
int N>
83expm_pade (
const Eigen::Matrix<T,N,N>& a) {
94 check_macro (a.is_symmetric(
d),
"exp: tensor should be symmetric");
96 using namespace Eigen;
98 if (
d == 2)
return exp2d (a);
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);
109 for (
size_t i = 0; i <
d; ++i)
110 for (
size_t j = 0; j <
d; ++j) b(i,j) = b1(i,j);
116#define _RHEOLEF_instanciation(T) \
117template tensor_basic<T> exp (const tensor_basic<T>& a, size_t d); \
#define _RHEOLEF_instanciation(T, M, A)
see the Float page for the full documentation
see the tensor page for the full documentation
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)
Eigen::Matrix< T, N, N > expm_eig(const Eigen::Matrix< T, N, N > &a)