Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
tensor4-dexp.cc
Go to the documentation of this file.
1
21// 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...)
33namespace rheolef {
34
35template <class T>
36tensor4_basic<T>
37dexp (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) \
201template tensor4_basic<T> dexp (const tensor_basic<T>& a, size_t d); \
202
204
205}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
see the tensor4 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
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
tensor4_basic< T > dexp(const tensor_basic< T > &chi, size_t dim)