Rheolef  7.1
an efficient C++ finite element environment
tensor4-dexp.cc
Go to the documentation of this file.
1 // dexp(tensor) = d_exp(a)/d a : 2D is explicit while 3D is not yet available
22 //
23 #include "rheolef/tensor4.h"
24 // =========================================================================
25 // 2D case: explicit
26 // =========================================================================
27 //
28 // see tensor_exp_tst.mac
29 // see also maxima: matrixexp.usage allows explicit expm() computation in 2D
30 // http://www.ma.utexas.edu/pipermail/maxima/2006/003031.html
31 // /usr/share/maxima/5.27.0/share/linearalgebra/matrixexp.usage
32 // refs. KanGueFor-2009, p. 50 (is buggy: missing 2 factor in A00 & A11...)
33 namespace rheolef {
34 
35 template <class T>
36 tensor4_basic<T>
37 dexp (const tensor_basic<T>& chi, size_t dim) {
38  check_macro (dim==2, "dexp(d="<<dim<<"): only d=2 case available yet, sorry");
39  static Float eps = 1e3*std::numeric_limits<Float>::epsilon();
40  tensor4 E;
41  Float a = chi(0,0), b = chi(1,0), c = chi(1,1);
42  T a2=sqr(a),
43  b2=sqr(b),
44  c2=sqr(c);
45  Float d2 = sqr(a-c)+4*b2;
46  Float d = sqrt(d2);
47  if (abs(d) < eps) { // chi = a*I, degenerate case, see tensor4_dexp_tst.mac
48  E(0,0, 0,0) =
49  E(1,1, 1,1) =
50  E(0,1, 0,1) =
51  E(0,1, 1,0) =
52  E(1,0, 0,1) =
53  E(1,0, 1,0) = exp(a);
54  return E;
55  }
56 #ifdef TO_CLEAN
57  // --------------------------------------
58  // refs. KanGueFor-2009, p. 50 : bugs...
59  // --------------------------------------
60  Float d3 = d*d2;
61  Float eac = exp((a+c)/2);
62  Float a3 = sinh(d/2), a4 = cosh(d/2);
63  E(0,0, 0,0) = (eac/d3)
64  *((pow(a-c,3) + 4*sqr(b)*(a-c+1))*a3
65  + d*(sqr(a-c) + 2*sqr(b))*a4);
66  E(0,0, 0,1) = E(0,0, 1,0)
67  = ((2*b*eac)/d3)
68  *((sqr(d) + 2*(c-a))*a3 + d*(a-c)*a4);
69  E(0,0, 1,1) = ((2*sqr(b)*eac)/d3)
70  *(-2*a3 + a4*d);
71  E(0,1, 0,0) = E(1,0, 0,0)
72  = ((b*eac)/d3)
73  *((sqr(d) + 2*(c-a))*a3 + d*(a-c)*a4);
74  E(0,1, 0,1) = E(1,0, 0,1) = E(0,1, 1,0) = E(1,0, 1,0)
75  = ((2*eac)/d3)
76  *(sqr(a-c)*a3 + 2*sqr(b)*d*a4);
77  E(0,1, 1,1) = E(1,0, 1,1)
78  = ((b*eac)/d3)
79  *((sqr(d) + 2*(a-c))*a3 + d*(c-a)*a4);
80  E(1,1, 0,0) = ((2*sqr(b)*eac)/d3)
81  *(-2*a3 + d*a4);
82  E(1,1, 0,1) = E(1,1, 1,0)
83  = ((2*b*eac)/d3)
84  *((sqr(d) + 2*(a-c))*a3 + d*(c-a)*a4);
85  E(1,1, 1,1) = ((2*eac)/d3)
86  *((pow(a-c,3) + 4*sqr(b)*(a-c-1))*a3
87  - d*(sqr(a-c) + 2*sqr(b))*a4);
88 #endif // TO_CLEAN
89  // ----------------------------------------
90  // geneated by maxima: see tensor4-exp.mac
91  // ----------------------------------------
92  // TODO: symmetriser !!!!!!!!!!!!!!!!!!!!
93  // E00ij
94  // -----
95  {
96  T k5=1/d2,
97  k7=exp(-d/2.0+c/2.0+a/2.0),
98  k8=2*a-2*c,
99  k9=1/d,
100  k10=exp(d),
101  k11=c*d*k10-a*d*k10-c2*k10+2*a*c*k10-4*b2*k10-a2*k10-c*d+a*d-c2+2*a*c-4*b2-a2;
102  E(0,0,0,0)
103  = -k5*(1.0/2.0-k8*k9/4.0)*k7*k11/2.0+k8*k7*k11/sqr(d2)/2.0-k5*k7*(-d*k10-k8*c2*k9*k10/2.0+a*k8*c*k9*k10+k8*c*k9*k10/2.0-2*b2*k8*k9*k10-a2*k8*k9*k10/2.0-a*k8*k9*k10/2.0+k8*c*k10/2.0+2*c*k10-a*k8*k10/2.0-2*a*k10+d-k8*c*k9/2.0+a*k8*k9/2.0+2*c-2*a)/2.0;
104  }
105  {
106  T k6=exp(-d/2.0+c/2.0+a/2.0),
107  k7=1/d,
108  k8=exp(d),
109  k9=c*d*k8-a*d*k8-c2*k8+2*a*c*k8-4*b2*k8-a2*k8-c*d+a*d-c2+2*a*c-4*b2-a2;
110  E(0,0, 0,1) = E(0,0, 1,0)
111  = b*k6*k9/pow(d,3)+4*b*k6*k9/sqr(d2)-k6*(-4*b*c2*k7*k8+8*a*b*c*k7*k8+4*b*c*k7*k8-16*pow(b,3)*k7*k8-4*a2*b*k7*k8-4*a*b*k7*k8+4*b*c*k8-4*a*b*k8-8*b*k8-4*b*c*k7+4*a*b*k7-8*b)/d2/2.0;
112  }
113  {
114  T k5=1/d2,
115  k7=exp(-d/2.0+c/2.0+a/2.0),
116  k8=2*c-2*a,
117  k9=1/d,
118  k10=exp(d),
119  k11=c*d*k10-a*d*k10-c2*k10+2*a*c*k10-4*b2*k10-a2*k10-c*d+a*d-c2+2*a*c-4*b2-a2;
120  E(0,0,1,1)
121  = -k5*(1.0/2.0-k8*k9/4.0)*k7*k11/2.0+k8*k7*k11/sqr(d2)/2.0-k5*k7*(d*k10-c2*k8*k9*k10/2.0+a*c*k8*k9*k10+c*k8*k9*k10/2.0-2*b2*k8*k9*k10-a2*k8*k9*k10/2.0-a*k8*k9*k10/2.0+c*k8*k10/2.0-a*k8*k10/2.0-2*c*k10+2*a*k10-d-c*k8*k9/2.0+a*k8*k9/2.0-2*c+2*a)/2.0;
122  }
123  // ----------------------------------------
124  // E01ij
125  // ----------------------------------------
126  {
127  T k1=2*a-2*c,
128  k3=a/2.0,
129  k4=c/2.0,
130  k6=exp(-d/2.0+k4+k3),
131  k7=exp(d)-1,
132  k8=1/d;
133  E(0,1, 0,0) = E(1,0, 0,0)
134  = b*k8*(1.0/2.0-k1*k8/4.0)*k6*k7-b*k1*k6*k7/pow(d,3)/2.0+b*k1*exp(d/2.0+k4+k3)/d2/2.0;
135  }
136  {
137  T k3=1/d2,
138  k4=a/2.0,
139  k5=c/2.0,
140  k7=exp(-d/2.0+k5+k4),
141  k8=exp(d)-1;
142  E(0,1, 0,1) = E(1,0, 0,1) = E(0,1, 1,0) = E(1,0, 1,0)
143  = k7*k8/d-2*b2*k3*k7*k8-4*b2*k7*k8/pow(d,3)+4*b2*k3*exp(d/2.0+k5+k4);
144  }
145  {
146  T k1=2*c-2*a,
147  k3=a/2.0,
148  k4=c/2.0,
149  d=sqrt(d2),
150  k6=exp(-d/2.0+k4+k3),
151  k7=exp(d)-1,
152  k8=1/d;
153  E(0,1, 1,1) = E(1,0, 1,1)
154  = b*k8*(1.0/2.0-k1*k8/4.0)*k6*k7-b*k1*k6*k7/pow(d,3)/2.0+b*k1*exp(d/2.0+k4+k3)/d2/2.0;
155  }
156  // ----------------------------------------
157  // E11ij
158  // ----------------------------------------
159  {
160  T k3=4*b2,
161  k4=-2*a*c,
162  k7=1/d2,
163  k9=exp(-d/2.0+c/2.0+a/2.0),
164  k10=2*a,
165  k11=-2*c,
166  k12=k11+k10,
167  k13=1/d,
168  k14=exp(d),
169  k15=c*d*k14-a*d*k14+c2*k14-2*a*c*k14+4*b2*k14+a2*k14-c*d+a*d+c2+k4+k3+a2;
170  E(1,1,0,0) = k7*(1.0/2.0-k12*k13/4.0)*k9*k15/2.0-k12*k9*k15/sqr(d2)/2.0+k7*k9*(-d*k14+k12*c2*k13*k14/2.0-a*k12*c*k13*k14+k12*c*k13*k14/2.0+2*b2*k12*k13*k14+a2*k12*k13*k14/2.0-a*k12*k13*k14/2.0+k12*c*k14/2.0-2*c*k14-a*k12*k14/2.0+2*a*k14+d-k12*c*k13/2.0+a*k12*k13/2.0+k11+k10)/2.0;
171  }
172  {
173  T k3=4*b2,
174  k4=-2*a*c,
175  k8=exp(-d/2.0+c/2.0+a/2.0),
176  k9=1/d,
177  k10=exp(d),
178  k11=c*d*k10-a*d*k10+c2*k10-2*a*c*k10+4*b2*k10+a2*k10-c*d+a*d+c2+k4+k3+a2;
179  E(1,1,0,1) = E(1,1, 1,0)
180  = -b*k8*k11/pow(d,3)-4*b*k8*k11/sqr(d2)+k8*(4*b*c2*k9*k10-8*a*b*c*k9*k10+4*b*c*k9*k10+16*pow(b,3)*k9*k10+4*a2*b*k9*k10-4*a*b*k9*k10+4*b*c*k10-4*a*b*k10+8*b*k10-4*b*c*k9+4*a*b*k9+8*b)/d2/2.0;
181  }
182  {
183  T k3=4*b2,
184  k4=-2*a*c,
185  k7=1/d2,
186  k9=exp(-d/2.0+c/2.0+a/2.0),
187  k10=-2*a,
188  k11=2*c,
189  k12=k11+k10,
190  k13=1/d,
191  k14=exp(d),
192  k15=c*d*k14-a*d*k14+c2*k14-2*a*c*k14+4*b2*k14+a2*k14-c*d+a*d+c2+k4+k3+a2;
193  E(1,1,1,1) = k7*(1.0/2.0-k12*k13/4.0)*k9*k15/2.0-k12*k9*k15/sqr(d2)/2.0+k7*k9*(d*k14+c2*k12*k13*k14/2.0-a*c*k12*k13*k14+c*k12*k13*k14/2.0+2*b2*k12*k13*k14+a2*k12*k13*k14/2.0-a*k12*k13*k14/2.0+c*k12*k14/2.0-a*k12*k14/2.0+2*c*k14-2*a*k14-d-c*k12*k13/2.0+a*k12*k13/2.0+k11+k10)/2.0;
194  }
195  return E;
196 }
197 // ----------------------------------------------------------------------------
198 // instanciation in library
199 // ----------------------------------------------------------------------------
200 #define _RHEOLEF_instanciation(T) \
201 template tensor4_basic<T> dexp (const tensor_basic<T>& a, size_t d); \
202 
204 
205 }// namespace rheolef
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
chi
Definition: transport_tensor_exact.icc:65
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
tensor4
see the tensor4 page for the full documentation
rheolef::tensor_basic
Definition: tensor.h:90
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::dexp
tensor4_basic< T > dexp(const tensor_basic< T > &chi, size_t dim)
Definition: tensor4-dexp.cc:37
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
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
T
Expr1::float_type T
Definition: field_expr.h:261