CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
umat_single_crystal_creep.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine umat_single_crystal_creep (amat, iel, iint, kode, elconloc, emec, emec0, beta, xokl, voj, xkl, vj, ithermal, t1l, dtime, time, ttime, icmd, ielas, mi, nstate_, xstateini, xstate, stre, stiff, iorien, pgauss, orab, pnewdt)
 

Function/Subroutine Documentation

◆ umat_single_crystal_creep()

subroutine umat_single_crystal_creep ( character*80  amat,
integer  iel,
integer  iint,
integer  kode,
real*8, dimension(21)  elconloc,
real*8, dimension(6)  emec,
real*8, dimension(6)  emec0,
real*8, dimension(6)  beta,
real*8, dimension(3,3)  xokl,
real*8  voj,
real*8, dimension(3,3)  xkl,
real*8  vj,
integer  ithermal,
real*8  t1l,
real*8  dtime,
real*8  time,
real*8  ttime,
integer  icmd,
integer  ielas,
integer, dimension(*)  mi,
integer  nstate_,
real*8, dimension(nstate_,mi(1),*)  xstateini,
real*8, dimension(nstate_,mi(1),*)  xstate,
real*8, dimension(6)  stre,
real*8, dimension(21)  stiff,
integer  iorien,
real*8, dimension(3)  pgauss,
real*8, dimension(7,*)  orab,
real*8  pnewdt 
)
24 !
25 ! calculates stiffness and stresses for a user defined material
26 ! law
27 !
28 ! icmd=3: calcutates stress at mechanical strain
29 ! else: calculates stress at mechanical strain and the stiffness
30 ! matrix
31 !
32 ! INPUT:
33 !
34 ! amat material name
35 ! iel element number
36 ! iint integration point number
37 !
38 ! kode material type (-100-#of constants entered
39 ! under *USER MATERIAL): can be used for materials
40 ! with varying number of constants
41 !
42 ! elconloc(21) user defined constants defined by the keyword
43 ! card *USER MATERIAL (max. 21, actual # =
44 ! -kode-100), interpolated for the
45 ! actual temperature t1l
46 !
47 ! emec(6) Lagrange mechanical strain tensor (component order:
48 ! 11,22,33,12,13,23) at the end of the increment
49 ! (thermal strains are subtracted)
50 ! emec0(6) Lagrange mechanical strain tensor at the start of the
51 ! increment (thermal strains are subtracted)
52 ! beta(6) residual stress tensor (the stress entered under
53 ! the keyword *INITIAL CONDITIONS,TYPE=STRESS)
54 !
55 ! xokl(3,3) deformation gradient at the start of the increment
56 ! voj Jacobian at the start of the increment
57 ! xkl(3,3) deformation gradient at the end of the increment
58 ! vj Jacobian at the end of the increment
59 !
60 ! ithermal 0: no thermal effects are taken into account
61 ! >0: thermal effects are taken into account (triggered
62 ! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
63 ! t1l temperature at the end of the increment
64 ! dtime time length of the increment
65 ! time step time at the end of the current increment
66 ! ttime total time at the start of the current step
67 !
68 ! icmd not equal to 3: calculate stress and stiffness
69 ! 3: calculate only stress
70 ! ielas 0: no elastic iteration: irreversible effects
71 ! are allowed
72 ! 1: elastic iteration, i.e. no irreversible
73 ! deformation allowed
74 !
75 ! mi(1) max. # of integration points per element in the
76 ! model
77 ! nstate_ max. # of state variables in the model
78 !
79 ! xstateini(nstate_,mi(1),# of elements)
80 ! state variables at the start of the increment
81 ! xstate(nstate_,mi(1),# of elements)
82 ! state variables at the end of the increment
83 !
84 ! stre(6) Piola-Kirchhoff stress of the second kind
85 ! at the start of the increment
86 !
87 ! iorien number of the local coordinate axis system
88 ! in the integration point at stake (takes the value
89 ! 0 if no local system applies)
90 ! pgauss(3) global coordinates of the integration point
91 ! orab(7,*) description of all local coordinate systems.
92 ! If a local coordinate system applies the global
93 ! tensors can be obtained by premultiplying the local
94 ! tensors with skl(3,3). skl is determined by calling
95 ! the subroutine transformatrix:
96 ! call transformatrix(orab(1,iorien),pgauss,skl)
97 !
98 !
99 ! OUTPUT:
100 !
101 ! xstate(nstate_,mi(1),# of elements)
102 ! updated state variables at the end of the increment
103 ! stre(6) Piola-Kirchhoff stress of the second kind at the
104 ! end of the increment
105 ! stiff(21): consistent tangent stiffness matrix in the material
106 ! frame of reference at the end of the increment. In
107 ! other words: the derivative of the PK2 stress with
108 ! respect to the Lagrangian strain tensor. The matrix
109 ! is supposed to be symmetric, only the upper half is
110 ! to be given in the same order as for a fully
111 ! anisotropic elastic material (*ELASTIC,TYPE=ANISO).
112 ! Notice that the matrix is an integral part of the
113 ! fourth order material tensor, i.e. the Voigt notation
114 ! is not used.
115 !
116  implicit none
117 !
118  integer convergence
119 !
120  character*80 amat
121 !
122  integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien
123 !
124  integer i,j,k,l,ipiv(18),info,neq,lda,ldb,
125  & nrhs,iplas,icounter
126 !
127  real*8 ep0(6),ep(6),pnewdt,
128  & dg(18),ddg(18),xm(6,18),ck(18),cn(18),
129  & stri(6),htri(18),sg(18),r(6),xmc(6,18),
130  & t(6),gl(18,18),gr(18,18),ee(6),c1111,c1122,c1212,dd,
131  & skl(3,3),xmtran(3,3),ddsdde(6,6),xx(6,18)
132 !
133  real*8 elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
134  & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
135  & elas(21),time,ttime
136 !
137  real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
138 !
139 ! crystallographic slip planes:
140 !
141 ! 1. n=1,1,1 l=1,-1,0
142 ! 2. n=1,1,1 l=1,0,-1
143 ! 3. n=1,1,1 l=0,1,-1
144 ! 4. n=1,-1,1 l=0,1,1
145 ! 5. n=1,-1,1 l=1,0,-1
146 ! 6. n=1,-1,1 l=1,1,0
147 ! 7. n=1,-1,-1 l=0,1,-1
148 ! 8. n=1,-1,-1 l=1,0,1
149 ! 9. n=1,-1,-1 l=1,1,0
150 ! 10. n=1,1,-1 l=0,1,1
151 ! 11. n=1,1,-1 l=1,0,1
152 ! 12. n=1,1,-1 l=1,-1,0
153 ! 13. n=1,0,0 l=0,1,1
154 ! 14. n=1,0,0 l=0,1,-1
155 ! 15. n=0,1,0 l=1,0,1
156 ! 16. n=0,1,0 l=1,0,-1
157 ! 17. n=0,0,1 l=1,1,0
158 ! 18. n=0,0,1 l=1,-1,0
159 !
160  xm=reshape((
161  & /0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
162  & 0.0000000000000e+00, 0.2041241452319e+00,-0.2041241452319e+00,
163  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
164  & 0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
165  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
166  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
167  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
168  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
169  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
170  & -0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
171  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
172  & 0.0000000000000e+00, 0.2041241452319e+00, 0.2041241452319e+00,
173  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
174  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
175  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
176  & -0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
177  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
178  & 0.0000000000000e+00,-0.2041241452319e+00,-0.2041241452319e+00,
179  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
180  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
181  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
182  & 0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
183  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
184  & 0.0000000000000e+00,-0.2041241452319e+00, 0.2041241452319e+00,
185  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
186  & 0.3535533905933e+00, 0.3535533905933e+00, 0.0000000000000e+00,
187  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
188  & 0.3535533905933e+00,-0.3535533905933e+00, 0.0000000000000e+00,
189  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
190  & 0.3535533905933e+00, 0.0000000000000e+00, 0.3535533905933e+00,
191  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
192  & 0.3535533905933e+00, 0.0000000000000e+00,-0.3535533905933e+00,
193  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
194  & 0.0000000000000e+00, 0.3535533905933e+00, 0.3535533905933e+00,
195  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
196  & 0.0000000000000e+00, 0.3535533905933e+00,-0.3535533905933e+00/
197  & ),(/6,18/))
198 !
199  xx=reshape((
200  & /0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
201  & 0.0000000000000e+00, 0.2041241452319e+00,-0.2041241452319e+00,
202  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
203  & 0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
204  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
205  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
206  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
207  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
208  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
209  & -0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
210  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
211  & 0.0000000000000e+00, 0.2041241452319e+00, 0.2041241452319e+00,
212  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
213  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
214  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
215  & -0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
216  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
217  & 0.0000000000000e+00,-0.2041241452319e+00,-0.2041241452319e+00,
218  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
219  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
220  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
221  & 0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
222  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
223  & 0.0000000000000e+00,-0.2041241452319e+00, 0.2041241452319e+00,
224  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
225  & 0.3535533905933e+00, 0.3535533905933e+00, 0.0000000000000e+00,
226  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
227  & 0.3535533905933e+00,-0.3535533905933e+00, 0.0000000000000e+00,
228  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
229  & 0.3535533905933e+00, 0.0000000000000e+00, 0.3535533905933e+00,
230  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
231  & 0.3535533905933e+00, 0.0000000000000e+00,-0.3535533905933e+00,
232  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
233  & 0.0000000000000e+00, 0.3535533905933e+00, 0.3535533905933e+00,
234  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
235  & 0.0000000000000e+00, 0.3535533905933e+00,-0.3535533905933e+00/
236  & ),(/6,18/))
237 !
238 c pnewdt=-1.d0
239 !
240 ! elastic constants
241 !
242  c1111=elconloc(1)
243  c1122=elconloc(2)
244  c1212=elconloc(3)
245 !
246  if(iorien.gt.0) then
247  call transformatrix(orab(1,iorien),pgauss,skl)
248  do k=1,18
249  do i=1,3
250  do j=i,3
251  xmtran(i,j)=skl(i,1)*skl(j,1)*xx(1,k)+
252  & skl(i,2)*skl(j,2)*xx(2,k)+
253  & skl(i,3)*skl(j,3)*xx(3,k)+
254  & (skl(i,1)*skl(j,2)+
255  & skl(i,2)*skl(j,1))*xx(4,k)+
256  & (skl(i,1)*skl(j,3)+
257  & skl(i,3)*skl(j,1))*xx(5,k)+
258  & (skl(i,2)*skl(j,3)+
259  & skl(i,3)*skl(j,2))*xx(6,k)
260  enddo
261  enddo
262  xm(1,k)=xmtran(1,1)
263  xm(2,k)=xmtran(2,2)
264  xm(3,k)=xmtran(3,3)
265  xm(4,k)=xmtran(1,2)
266  xm(5,k)=xmtran(1,3)
267  xm(6,k)=xmtran(2,3)
268  enddo
269 !
270  elas( 1)=
271  & skl(1,1)*skl(1,1)*skl(1,1)*skl(1,1)*c1111+
272  & skl(1,1)*skl(1,1)*skl(1,2)*skl(1,2)*c1122+
273  & skl(1,1)*skl(1,1)*skl(1,3)*skl(1,3)*c1122+
274  & skl(1,1)*skl(1,2)*skl(1,1)*skl(1,2)*c1212+
275  & skl(1,1)*skl(1,2)*skl(1,2)*skl(1,1)*c1212+
276  & skl(1,1)*skl(1,3)*skl(1,1)*skl(1,3)*c1212+
277  & skl(1,1)*skl(1,3)*skl(1,3)*skl(1,1)*c1212+
278  & skl(1,2)*skl(1,1)*skl(1,1)*skl(1,2)*c1212+
279  & skl(1,2)*skl(1,1)*skl(1,2)*skl(1,1)*c1212+
280  & skl(1,2)*skl(1,2)*skl(1,1)*skl(1,1)*c1122+
281  & skl(1,2)*skl(1,2)*skl(1,2)*skl(1,2)*c1111+
282  & skl(1,2)*skl(1,2)*skl(1,3)*skl(1,3)*c1122+
283  & skl(1,2)*skl(1,3)*skl(1,2)*skl(1,3)*c1212+
284  & skl(1,2)*skl(1,3)*skl(1,3)*skl(1,2)*c1212+
285  & skl(1,3)*skl(1,1)*skl(1,1)*skl(1,3)*c1212+
286  & skl(1,3)*skl(1,1)*skl(1,3)*skl(1,1)*c1212+
287  & skl(1,3)*skl(1,2)*skl(1,2)*skl(1,3)*c1212+
288  & skl(1,3)*skl(1,2)*skl(1,3)*skl(1,2)*c1212+
289  & skl(1,3)*skl(1,3)*skl(1,1)*skl(1,1)*c1122+
290  & skl(1,3)*skl(1,3)*skl(1,2)*skl(1,2)*c1122+
291  & skl(1,3)*skl(1,3)*skl(1,3)*skl(1,3)*c1111
292  elas( 2)=
293  & skl(1,1)*skl(1,1)*skl(2,1)*skl(2,1)*c1111+
294  & skl(1,1)*skl(1,1)*skl(2,2)*skl(2,2)*c1122+
295  & skl(1,1)*skl(1,1)*skl(2,3)*skl(2,3)*c1122+
296  & skl(1,1)*skl(1,2)*skl(2,1)*skl(2,2)*c1212+
297  & skl(1,1)*skl(1,2)*skl(2,2)*skl(2,1)*c1212+
298  & skl(1,1)*skl(1,3)*skl(2,1)*skl(2,3)*c1212+
299  & skl(1,1)*skl(1,3)*skl(2,3)*skl(2,1)*c1212+
300  & skl(1,2)*skl(1,1)*skl(2,1)*skl(2,2)*c1212+
301  & skl(1,2)*skl(1,1)*skl(2,2)*skl(2,1)*c1212+
302  & skl(1,2)*skl(1,2)*skl(2,1)*skl(2,1)*c1122+
303  & skl(1,2)*skl(1,2)*skl(2,2)*skl(2,2)*c1111+
304  & skl(1,2)*skl(1,2)*skl(2,3)*skl(2,3)*c1122+
305  & skl(1,2)*skl(1,3)*skl(2,2)*skl(2,3)*c1212+
306  & skl(1,2)*skl(1,3)*skl(2,3)*skl(2,2)*c1212+
307  & skl(1,3)*skl(1,1)*skl(2,1)*skl(2,3)*c1212+
308  & skl(1,3)*skl(1,1)*skl(2,3)*skl(2,1)*c1212+
309  & skl(1,3)*skl(1,2)*skl(2,2)*skl(2,3)*c1212+
310  & skl(1,3)*skl(1,2)*skl(2,3)*skl(2,2)*c1212+
311  & skl(1,3)*skl(1,3)*skl(2,1)*skl(2,1)*c1122+
312  & skl(1,3)*skl(1,3)*skl(2,2)*skl(2,2)*c1122+
313  & skl(1,3)*skl(1,3)*skl(2,3)*skl(2,3)*c1111
314  elas( 3)=
315  & skl(2,1)*skl(2,1)*skl(2,1)*skl(2,1)*c1111+
316  & skl(2,1)*skl(2,1)*skl(2,2)*skl(2,2)*c1122+
317  & skl(2,1)*skl(2,1)*skl(2,3)*skl(2,3)*c1122+
318  & skl(2,1)*skl(2,2)*skl(2,1)*skl(2,2)*c1212+
319  & skl(2,1)*skl(2,2)*skl(2,2)*skl(2,1)*c1212+
320  & skl(2,1)*skl(2,3)*skl(2,1)*skl(2,3)*c1212+
321  & skl(2,1)*skl(2,3)*skl(2,3)*skl(2,1)*c1212+
322  & skl(2,2)*skl(2,1)*skl(2,1)*skl(2,2)*c1212+
323  & skl(2,2)*skl(2,1)*skl(2,2)*skl(2,1)*c1212+
324  & skl(2,2)*skl(2,2)*skl(2,1)*skl(2,1)*c1122+
325  & skl(2,2)*skl(2,2)*skl(2,2)*skl(2,2)*c1111+
326  & skl(2,2)*skl(2,2)*skl(2,3)*skl(2,3)*c1122+
327  & skl(2,2)*skl(2,3)*skl(2,2)*skl(2,3)*c1212+
328  & skl(2,2)*skl(2,3)*skl(2,3)*skl(2,2)*c1212+
329  & skl(2,3)*skl(2,1)*skl(2,1)*skl(2,3)*c1212+
330  & skl(2,3)*skl(2,1)*skl(2,3)*skl(2,1)*c1212+
331  & skl(2,3)*skl(2,2)*skl(2,2)*skl(2,3)*c1212+
332  & skl(2,3)*skl(2,2)*skl(2,3)*skl(2,2)*c1212+
333  & skl(2,3)*skl(2,3)*skl(2,1)*skl(2,1)*c1122+
334  & skl(2,3)*skl(2,3)*skl(2,2)*skl(2,2)*c1122+
335  & skl(2,3)*skl(2,3)*skl(2,3)*skl(2,3)*c1111
336  elas( 4)=
337  & skl(1,1)*skl(1,1)*skl(3,1)*skl(3,1)*c1111+
338  & skl(1,1)*skl(1,1)*skl(3,2)*skl(3,2)*c1122+
339  & skl(1,1)*skl(1,1)*skl(3,3)*skl(3,3)*c1122+
340  & skl(1,1)*skl(1,2)*skl(3,1)*skl(3,2)*c1212+
341  & skl(1,1)*skl(1,2)*skl(3,2)*skl(3,1)*c1212+
342  & skl(1,1)*skl(1,3)*skl(3,1)*skl(3,3)*c1212+
343  & skl(1,1)*skl(1,3)*skl(3,3)*skl(3,1)*c1212+
344  & skl(1,2)*skl(1,1)*skl(3,1)*skl(3,2)*c1212+
345  & skl(1,2)*skl(1,1)*skl(3,2)*skl(3,1)*c1212+
346  & skl(1,2)*skl(1,2)*skl(3,1)*skl(3,1)*c1122+
347  & skl(1,2)*skl(1,2)*skl(3,2)*skl(3,2)*c1111+
348  & skl(1,2)*skl(1,2)*skl(3,3)*skl(3,3)*c1122+
349  & skl(1,2)*skl(1,3)*skl(3,2)*skl(3,3)*c1212+
350  & skl(1,2)*skl(1,3)*skl(3,3)*skl(3,2)*c1212+
351  & skl(1,3)*skl(1,1)*skl(3,1)*skl(3,3)*c1212+
352  & skl(1,3)*skl(1,1)*skl(3,3)*skl(3,1)*c1212+
353  & skl(1,3)*skl(1,2)*skl(3,2)*skl(3,3)*c1212+
354  & skl(1,3)*skl(1,2)*skl(3,3)*skl(3,2)*c1212+
355  & skl(1,3)*skl(1,3)*skl(3,1)*skl(3,1)*c1122+
356  & skl(1,3)*skl(1,3)*skl(3,2)*skl(3,2)*c1122+
357  & skl(1,3)*skl(1,3)*skl(3,3)*skl(3,3)*c1111
358  elas( 5)=
359  & skl(2,1)*skl(2,1)*skl(3,1)*skl(3,1)*c1111+
360  & skl(2,1)*skl(2,1)*skl(3,2)*skl(3,2)*c1122+
361  & skl(2,1)*skl(2,1)*skl(3,3)*skl(3,3)*c1122+
362  & skl(2,1)*skl(2,2)*skl(3,1)*skl(3,2)*c1212+
363  & skl(2,1)*skl(2,2)*skl(3,2)*skl(3,1)*c1212+
364  & skl(2,1)*skl(2,3)*skl(3,1)*skl(3,3)*c1212+
365  & skl(2,1)*skl(2,3)*skl(3,3)*skl(3,1)*c1212+
366  & skl(2,2)*skl(2,1)*skl(3,1)*skl(3,2)*c1212+
367  & skl(2,2)*skl(2,1)*skl(3,2)*skl(3,1)*c1212+
368  & skl(2,2)*skl(2,2)*skl(3,1)*skl(3,1)*c1122+
369  & skl(2,2)*skl(2,2)*skl(3,2)*skl(3,2)*c1111+
370  & skl(2,2)*skl(2,2)*skl(3,3)*skl(3,3)*c1122+
371  & skl(2,2)*skl(2,3)*skl(3,2)*skl(3,3)*c1212+
372  & skl(2,2)*skl(2,3)*skl(3,3)*skl(3,2)*c1212+
373  & skl(2,3)*skl(2,1)*skl(3,1)*skl(3,3)*c1212+
374  & skl(2,3)*skl(2,1)*skl(3,3)*skl(3,1)*c1212+
375  & skl(2,3)*skl(2,2)*skl(3,2)*skl(3,3)*c1212+
376  & skl(2,3)*skl(2,2)*skl(3,3)*skl(3,2)*c1212+
377  & skl(2,3)*skl(2,3)*skl(3,1)*skl(3,1)*c1122+
378  & skl(2,3)*skl(2,3)*skl(3,2)*skl(3,2)*c1122+
379  & skl(2,3)*skl(2,3)*skl(3,3)*skl(3,3)*c1111
380  elas( 6)=
381  & skl(3,1)*skl(3,1)*skl(3,1)*skl(3,1)*c1111+
382  & skl(3,1)*skl(3,1)*skl(3,2)*skl(3,2)*c1122+
383  & skl(3,1)*skl(3,1)*skl(3,3)*skl(3,3)*c1122+
384  & skl(3,1)*skl(3,2)*skl(3,1)*skl(3,2)*c1212+
385  & skl(3,1)*skl(3,2)*skl(3,2)*skl(3,1)*c1212+
386  & skl(3,1)*skl(3,3)*skl(3,1)*skl(3,3)*c1212+
387  & skl(3,1)*skl(3,3)*skl(3,3)*skl(3,1)*c1212+
388  & skl(3,2)*skl(3,1)*skl(3,1)*skl(3,2)*c1212+
389  & skl(3,2)*skl(3,1)*skl(3,2)*skl(3,1)*c1212+
390  & skl(3,2)*skl(3,2)*skl(3,1)*skl(3,1)*c1122+
391  & skl(3,2)*skl(3,2)*skl(3,2)*skl(3,2)*c1111+
392  & skl(3,2)*skl(3,2)*skl(3,3)*skl(3,3)*c1122+
393  & skl(3,2)*skl(3,3)*skl(3,2)*skl(3,3)*c1212+
394  & skl(3,2)*skl(3,3)*skl(3,3)*skl(3,2)*c1212+
395  & skl(3,3)*skl(3,1)*skl(3,1)*skl(3,3)*c1212+
396  & skl(3,3)*skl(3,1)*skl(3,3)*skl(3,1)*c1212+
397  & skl(3,3)*skl(3,2)*skl(3,2)*skl(3,3)*c1212+
398  & skl(3,3)*skl(3,2)*skl(3,3)*skl(3,2)*c1212+
399  & skl(3,3)*skl(3,3)*skl(3,1)*skl(3,1)*c1122+
400  & skl(3,3)*skl(3,3)*skl(3,2)*skl(3,2)*c1122+
401  & skl(3,3)*skl(3,3)*skl(3,3)*skl(3,3)*c1111
402  elas( 7)=
403  & skl(1,1)*skl(1,1)*skl(1,1)*skl(2,1)*c1111+
404  & skl(1,1)*skl(1,1)*skl(1,2)*skl(2,2)*c1122+
405  & skl(1,1)*skl(1,1)*skl(1,3)*skl(2,3)*c1122+
406  & skl(1,1)*skl(1,2)*skl(1,1)*skl(2,2)*c1212+
407  & skl(1,1)*skl(1,2)*skl(1,2)*skl(2,1)*c1212+
408  & skl(1,1)*skl(1,3)*skl(1,1)*skl(2,3)*c1212+
409  & skl(1,1)*skl(1,3)*skl(1,3)*skl(2,1)*c1212+
410  & skl(1,2)*skl(1,1)*skl(1,1)*skl(2,2)*c1212+
411  & skl(1,2)*skl(1,1)*skl(1,2)*skl(2,1)*c1212+
412  & skl(1,2)*skl(1,2)*skl(1,1)*skl(2,1)*c1122+
413  & skl(1,2)*skl(1,2)*skl(1,2)*skl(2,2)*c1111+
414  & skl(1,2)*skl(1,2)*skl(1,3)*skl(2,3)*c1122+
415  & skl(1,2)*skl(1,3)*skl(1,2)*skl(2,3)*c1212+
416  & skl(1,2)*skl(1,3)*skl(1,3)*skl(2,2)*c1212+
417  & skl(1,3)*skl(1,1)*skl(1,1)*skl(2,3)*c1212+
418  & skl(1,3)*skl(1,1)*skl(1,3)*skl(2,1)*c1212+
419  & skl(1,3)*skl(1,2)*skl(1,2)*skl(2,3)*c1212+
420  & skl(1,3)*skl(1,2)*skl(1,3)*skl(2,2)*c1212+
421  & skl(1,3)*skl(1,3)*skl(1,1)*skl(2,1)*c1122+
422  & skl(1,3)*skl(1,3)*skl(1,2)*skl(2,2)*c1122+
423  & skl(1,3)*skl(1,3)*skl(1,3)*skl(2,3)*c1111
424  elas( 8)=
425  & skl(2,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+
426  & skl(2,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+
427  & skl(2,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+
428  & skl(2,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+
429  & skl(2,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+
430  & skl(2,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+
431  & skl(2,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+
432  & skl(2,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+
433  & skl(2,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+
434  & skl(2,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+
435  & skl(2,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+
436  & skl(2,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+
437  & skl(2,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+
438  & skl(2,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+
439  & skl(2,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+
440  & skl(2,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+
441  & skl(2,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+
442  & skl(2,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+
443  & skl(2,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+
444  & skl(2,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+
445  & skl(2,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111
446  elas( 9)=
447  & skl(3,1)*skl(3,1)*skl(1,1)*skl(2,1)*c1111+
448  & skl(3,1)*skl(3,1)*skl(1,2)*skl(2,2)*c1122+
449  & skl(3,1)*skl(3,1)*skl(1,3)*skl(2,3)*c1122+
450  & skl(3,1)*skl(3,2)*skl(1,1)*skl(2,2)*c1212+
451  & skl(3,1)*skl(3,2)*skl(1,2)*skl(2,1)*c1212+
452  & skl(3,1)*skl(3,3)*skl(1,1)*skl(2,3)*c1212+
453  & skl(3,1)*skl(3,3)*skl(1,3)*skl(2,1)*c1212+
454  & skl(3,2)*skl(3,1)*skl(1,1)*skl(2,2)*c1212+
455  & skl(3,2)*skl(3,1)*skl(1,2)*skl(2,1)*c1212+
456  & skl(3,2)*skl(3,2)*skl(1,1)*skl(2,1)*c1122+
457  & skl(3,2)*skl(3,2)*skl(1,2)*skl(2,2)*c1111+
458  & skl(3,2)*skl(3,2)*skl(1,3)*skl(2,3)*c1122+
459  & skl(3,2)*skl(3,3)*skl(1,2)*skl(2,3)*c1212+
460  & skl(3,2)*skl(3,3)*skl(1,3)*skl(2,2)*c1212+
461  & skl(3,3)*skl(3,1)*skl(1,1)*skl(2,3)*c1212+
462  & skl(3,3)*skl(3,1)*skl(1,3)*skl(2,1)*c1212+
463  & skl(3,3)*skl(3,2)*skl(1,2)*skl(2,3)*c1212+
464  & skl(3,3)*skl(3,2)*skl(1,3)*skl(2,2)*c1212+
465  & skl(3,3)*skl(3,3)*skl(1,1)*skl(2,1)*c1122+
466  & skl(3,3)*skl(3,3)*skl(1,2)*skl(2,2)*c1122+
467  & skl(3,3)*skl(3,3)*skl(1,3)*skl(2,3)*c1111
468  elas(10)=
469  & skl(1,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+
470  & skl(1,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+
471  & skl(1,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+
472  & skl(1,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+
473  & skl(1,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+
474  & skl(1,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+
475  & skl(1,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+
476  & skl(1,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+
477  & skl(1,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+
478  & skl(1,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+
479  & skl(1,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+
480  & skl(1,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+
481  & skl(1,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+
482  & skl(1,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+
483  & skl(1,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+
484  & skl(1,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+
485  & skl(1,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+
486  & skl(1,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+
487  & skl(1,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+
488  & skl(1,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+
489  & skl(1,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111
490  elas(11)=
491  & skl(1,1)*skl(1,1)*skl(1,1)*skl(3,1)*c1111+
492  & skl(1,1)*skl(1,1)*skl(1,2)*skl(3,2)*c1122+
493  & skl(1,1)*skl(1,1)*skl(1,3)*skl(3,3)*c1122+
494  & skl(1,1)*skl(1,2)*skl(1,1)*skl(3,2)*c1212+
495  & skl(1,1)*skl(1,2)*skl(1,2)*skl(3,1)*c1212+
496  & skl(1,1)*skl(1,3)*skl(1,1)*skl(3,3)*c1212+
497  & skl(1,1)*skl(1,3)*skl(1,3)*skl(3,1)*c1212+
498  & skl(1,2)*skl(1,1)*skl(1,1)*skl(3,2)*c1212+
499  & skl(1,2)*skl(1,1)*skl(1,2)*skl(3,1)*c1212+
500  & skl(1,2)*skl(1,2)*skl(1,1)*skl(3,1)*c1122+
501  & skl(1,2)*skl(1,2)*skl(1,2)*skl(3,2)*c1111+
502  & skl(1,2)*skl(1,2)*skl(1,3)*skl(3,3)*c1122+
503  & skl(1,2)*skl(1,3)*skl(1,2)*skl(3,3)*c1212+
504  & skl(1,2)*skl(1,3)*skl(1,3)*skl(3,2)*c1212+
505  & skl(1,3)*skl(1,1)*skl(1,1)*skl(3,3)*c1212+
506  & skl(1,3)*skl(1,1)*skl(1,3)*skl(3,1)*c1212+
507  & skl(1,3)*skl(1,2)*skl(1,2)*skl(3,3)*c1212+
508  & skl(1,3)*skl(1,2)*skl(1,3)*skl(3,2)*c1212+
509  & skl(1,3)*skl(1,3)*skl(1,1)*skl(3,1)*c1122+
510  & skl(1,3)*skl(1,3)*skl(1,2)*skl(3,2)*c1122+
511  & skl(1,3)*skl(1,3)*skl(1,3)*skl(3,3)*c1111
512  elas(12)=
513  & skl(2,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+
514  & skl(2,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+
515  & skl(2,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+
516  & skl(2,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+
517  & skl(2,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+
518  & skl(2,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+
519  & skl(2,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+
520  & skl(2,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+
521  & skl(2,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+
522  & skl(2,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+
523  & skl(2,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+
524  & skl(2,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+
525  & skl(2,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+
526  & skl(2,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+
527  & skl(2,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+
528  & skl(2,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+
529  & skl(2,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+
530  & skl(2,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+
531  & skl(2,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+
532  & skl(2,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+
533  & skl(2,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111
534  elas(13)=
535  & skl(3,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+
536  & skl(3,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+
537  & skl(3,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+
538  & skl(3,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+
539  & skl(3,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+
540  & skl(3,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+
541  & skl(3,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+
542  & skl(3,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+
543  & skl(3,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+
544  & skl(3,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+
545  & skl(3,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+
546  & skl(3,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+
547  & skl(3,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+
548  & skl(3,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+
549  & skl(3,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+
550  & skl(3,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+
551  & skl(3,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+
552  & skl(3,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+
553  & skl(3,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+
554  & skl(3,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+
555  & skl(3,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111
556  elas(14)=
557  & skl(1,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+
558  & skl(1,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+
559  & skl(1,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+
560  & skl(1,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+
561  & skl(1,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+
562  & skl(1,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+
563  & skl(1,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+
564  & skl(1,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+
565  & skl(1,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+
566  & skl(1,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+
567  & skl(1,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+
568  & skl(1,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+
569  & skl(1,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+
570  & skl(1,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+
571  & skl(1,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+
572  & skl(1,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+
573  & skl(1,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+
574  & skl(1,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+
575  & skl(1,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+
576  & skl(1,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+
577  & skl(1,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111
578  elas(15)=
579  & skl(1,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+
580  & skl(1,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+
581  & skl(1,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+
582  & skl(1,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+
583  & skl(1,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+
584  & skl(1,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+
585  & skl(1,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+
586  & skl(1,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+
587  & skl(1,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+
588  & skl(1,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+
589  & skl(1,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+
590  & skl(1,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+
591  & skl(1,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+
592  & skl(1,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+
593  & skl(1,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+
594  & skl(1,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+
595  & skl(1,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+
596  & skl(1,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+
597  & skl(1,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+
598  & skl(1,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+
599  & skl(1,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111
600  elas(16)=
601  & skl(1,1)*skl(1,1)*skl(2,1)*skl(3,1)*c1111+
602  & skl(1,1)*skl(1,1)*skl(2,2)*skl(3,2)*c1122+
603  & skl(1,1)*skl(1,1)*skl(2,3)*skl(3,3)*c1122+
604  & skl(1,1)*skl(1,2)*skl(2,1)*skl(3,2)*c1212+
605  & skl(1,1)*skl(1,2)*skl(2,2)*skl(3,1)*c1212+
606  & skl(1,1)*skl(1,3)*skl(2,1)*skl(3,3)*c1212+
607  & skl(1,1)*skl(1,3)*skl(2,3)*skl(3,1)*c1212+
608  & skl(1,2)*skl(1,1)*skl(2,1)*skl(3,2)*c1212+
609  & skl(1,2)*skl(1,1)*skl(2,2)*skl(3,1)*c1212+
610  & skl(1,2)*skl(1,2)*skl(2,1)*skl(3,1)*c1122+
611  & skl(1,2)*skl(1,2)*skl(2,2)*skl(3,2)*c1111+
612  & skl(1,2)*skl(1,2)*skl(2,3)*skl(3,3)*c1122+
613  & skl(1,2)*skl(1,3)*skl(2,2)*skl(3,3)*c1212+
614  & skl(1,2)*skl(1,3)*skl(2,3)*skl(3,2)*c1212+
615  & skl(1,3)*skl(1,1)*skl(2,1)*skl(3,3)*c1212+
616  & skl(1,3)*skl(1,1)*skl(2,3)*skl(3,1)*c1212+
617  & skl(1,3)*skl(1,2)*skl(2,2)*skl(3,3)*c1212+
618  & skl(1,3)*skl(1,2)*skl(2,3)*skl(3,2)*c1212+
619  & skl(1,3)*skl(1,3)*skl(2,1)*skl(3,1)*c1122+
620  & skl(1,3)*skl(1,3)*skl(2,2)*skl(3,2)*c1122+
621  & skl(1,3)*skl(1,3)*skl(2,3)*skl(3,3)*c1111
622  elas(17)=
623  & skl(2,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+
624  & skl(2,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+
625  & skl(2,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+
626  & skl(2,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+
627  & skl(2,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+
628  & skl(2,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+
629  & skl(2,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+
630  & skl(2,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+
631  & skl(2,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+
632  & skl(2,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+
633  & skl(2,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+
634  & skl(2,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+
635  & skl(2,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+
636  & skl(2,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+
637  & skl(2,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+
638  & skl(2,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+
639  & skl(2,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+
640  & skl(2,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+
641  & skl(2,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+
642  & skl(2,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+
643  & skl(2,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111
644  elas(18)=
645  & skl(3,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
646  & skl(3,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
647  & skl(3,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
648  & skl(3,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
649  & skl(3,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
650  & skl(3,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
651  & skl(3,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
652  & skl(3,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
653  & skl(3,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
654  & skl(3,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
655  & skl(3,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
656  & skl(3,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
657  & skl(3,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
658  & skl(3,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
659  & skl(3,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
660  & skl(3,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
661  & skl(3,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
662  & skl(3,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
663  & skl(3,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
664  & skl(3,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
665  & skl(3,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
666  elas(19)=
667  & skl(1,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+
668  & skl(1,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+
669  & skl(1,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+
670  & skl(1,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+
671  & skl(1,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+
672  & skl(1,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+
673  & skl(1,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+
674  & skl(1,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+
675  & skl(1,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+
676  & skl(1,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+
677  & skl(1,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+
678  & skl(1,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+
679  & skl(1,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+
680  & skl(1,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+
681  & skl(1,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+
682  & skl(1,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+
683  & skl(1,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+
684  & skl(1,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+
685  & skl(1,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+
686  & skl(1,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+
687  & skl(1,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111
688  elas(20)=
689  & skl(1,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
690  & skl(1,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
691  & skl(1,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
692  & skl(1,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
693  & skl(1,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
694  & skl(1,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
695  & skl(1,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
696  & skl(1,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
697  & skl(1,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
698  & skl(1,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
699  & skl(1,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
700  & skl(1,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
701  & skl(1,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
702  & skl(1,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
703  & skl(1,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
704  & skl(1,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
705  & skl(1,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
706  & skl(1,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
707  & skl(1,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
708  & skl(1,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
709  & skl(1,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
710  elas(21)=
711  & skl(2,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
712  & skl(2,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
713  & skl(2,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
714  & skl(2,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
715  & skl(2,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
716  & skl(2,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
717  & skl(2,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
718  & skl(2,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
719  & skl(2,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
720  & skl(2,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
721  & skl(2,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
722  & skl(2,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
723  & skl(2,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
724  & skl(2,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
725  & skl(2,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
726  & skl(2,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
727  & skl(2,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
728  & skl(2,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
729  & skl(2,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
730  & skl(2,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
731  & skl(2,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
732  endif
733 !
734  do i=1,6
735  ep0(i)=xstateini(i,iint,iel)
736  enddo
737 !
738 ! elastic strains
739 !
740  do i=1,6
741  ee(i)=emec(i)-ep0(i)
742  enddo
743 !
744 ! (visco)plastic constants: octahedral slip system
745 !
746  do i=1,12
747  ck(i)=elconloc(4)
748  cn(i)=elconloc(5)
749  enddo
750 !
751 ! (visco)plastic constants: cubic slip system
752 !
753  do i=13,18
754  ck(i)=elconloc(6)
755  cn(i)=elconloc(7)
756  enddo
757 !
758 ! global trial stress tensor
759 !
760  if(iorien.gt.0) then
761  stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+
762  & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6))
763  & -beta(1)
764  stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+
765  & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6))
766  & -beta(2)
767  stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+
768  & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6))
769  & -beta(3)
770  stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+
771  & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6))
772  & -beta(4)
773  stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+
774  & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6))
775  & -beta(5)
776  stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+
777  & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6))
778  & -beta(6)
779  else
780  stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1)
781  stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2)
782  stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3)
783  stri(4)=2.d0*c1212*ee(4)-beta(4)
784  stri(5)=2.d0*c1212*ee(5)-beta(5)
785  stri(6)=2.d0*c1212*ee(6)-beta(6)
786  endif
787 !
788 ! stress radius in each slip plane
789 !
790  do i=1,18
791  sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+
792  & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6))
793  enddo
794 !
795 ! evaluation of the yield surface
796 !
797  do i=1,18
798  htri(i)=dabs(sg(i))
799  enddo
800 !
801  if(ielas.eq.1) then
802 !
803 ! elastic stress
804 !
805  do i=1,6
806  stre(i)=stri(i)
807  enddo
808 !
809 ! elastic stiffness
810 !
811  if(icmd.ne.3) then
812  if(iorien.gt.0) then
813  do i=1,21
814  stiff(i)=elas(i)
815  enddo
816  else
817  stiff(1)=c1111
818  stiff(2)=c1122
819  stiff(3)=c1111
820  stiff(4)=c1122
821  stiff(5)=c1122
822  stiff(6)=c1111
823  stiff(7)=0.d0
824  stiff(8)=0.d0
825  stiff(9)=0.d0
826  stiff(10)=c1212
827  stiff(11)=0.d0
828  stiff(12)=0.d0
829  stiff(13)=0.d0
830  stiff(14)=0.d0
831  stiff(15)=c1212
832  stiff(16)=0.d0
833  stiff(17)=0.d0
834  stiff(18)=0.d0
835  stiff(19)=0.d0
836  stiff(20)=0.d0
837  stiff(21)=c1212
838  endif
839 !
840 ! updating the state variables
841 !
842  do i=1,6
843  xstate(i,iint,iel)=ep0(i)
844  enddo
845  do i=7,24
846  xstate(i,iint,iel)=xstateini(i,iint,iel)
847  enddo
848  endif
849 !
850  return
851  endif
852 !
853 ! plastic deformation
854 !
855  nrhs=1
856  lda=18
857  ldb=18
858 !
859 ! initializing the state variables
860 !
861  do i=1,6
862  ep(i)=ep0(i)
863  enddo
864  do i=1,18
865  dg(i)=0.d0
866 c dg(i)=dtime*htri(i)**cn(i)/ck(i)**cn(i)
867  enddo
868 !
869 ! major loop
870 !
871  icounter=0
872  do
873  icounter=icounter+1
874  if(icounter.gt.100) then
875  pnewdt=0.25d0
876  return
877 c write(*,*)
878 c & '*ERROR in umat_single_crystal_creep: no convergence'
879 c call exit(201)
880  endif
881 !
882 ! elastic strains
883 !
884  do i=1,6
885  ee(i)=emec(i)-ep(i)
886 c write(*,*) 'umat_single_crystal_creep ee',icounter,i,ee(i)
887  enddo
888 !
889 ! global trial stress tensor
890 !
891  if(iorien.gt.0) then
892  stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+
893  & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6))
894  & -beta(1)
895  stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+
896  & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6))
897  & -beta(2)
898  stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+
899  & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6))
900  & -beta(3)
901  stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+
902  & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6))
903  & -beta(4)
904  stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+
905  & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6))
906  & -beta(5)
907  stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+
908  & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6))
909  & -beta(6)
910  else
911  stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1)
912  stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2)
913  stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3)
914  stri(4)=2.d0*c1212*ee(4)-beta(4)
915  stri(5)=2.d0*c1212*ee(5)-beta(5)
916  stri(6)=2.d0*c1212*ee(6)-beta(6)
917  endif
918  do i=1,6
919 c write(*,*) 'umat_single_crystal_creep stri',icounter,i,stri(i)
920  enddo
921 !
922 ! stress radius in each slip plane
923 !
924  do i=1,18
925  sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+
926  & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6))
927 c write(*,*) 'umat_single_crystal_creep sg',icounter,i,sg(i)
928  enddo
929 !
930 ! evaluation of the yield surface
931 !
932  do i=1,18
933 c write(*,*) 'umat_single_crystal_creep ck,dtime,cn,dg',icounter,
934 c & i,ck(i),dtime,cn(i),dg(i)
935  htri(i)=dabs(sg(i))-ck(i)*(dg(i)/dtime)**(1.d0/cn(i))
936 c write(*,*) 'umat_single_crystal_creep htri',icounter,i,htri(i)
937  enddo
938 !
939 ! replace sg(i) by sgn(sg(i))
940 !
941  do i=1,18
942  if(sg(i).lt.0.d0) then
943  sg(i)=-1.d0
944  else
945  sg(i)=1.d0
946  endif
947  enddo
948 !
949 ! determining the residual matrix
950 !
951  do i=1,6
952  r(i)=ep0(i)-ep(i)
953  enddo
954  do i=1,18
955  do j=1,6
956  r(j)=r(j)+xm(j,i)*sg(i)*dg(i)
957  enddo
958  enddo
959 !
960 ! check convergence
961 !
962  if(icounter.gt.1) then
963  convergence=1
964  do i=1,18
965 c if(dabs(htri(i)).gt.1.d-10) then
966  if(dabs(htri(i)).gt.1.d-3) then
967 c write(*,*) 'divergence htri dg ',
968 c & icounter,i,htri(i),dg(i)
969  convergence=0
970  exit
971  endif
972  enddo
973  if(convergence.eq.1) exit
974  endif
975 c do i=1,18
976 c if(dabs(htri(i)).gt.1.d-5) then
977 c write(*,*) 'htri ',i,htri(i)
978 c convergence=0
979 c exit
980 c endif
981 c enddo
982 c if(convergence.eq.1) then
983 c dd=0.d0
984 c do i=1,6
985 c dd=dd+r(i)*r(i)
986 c enddo
987 c dd=sqrt(dd)
988 c if(dd.gt.1.d-10) then
989 c convergence=0
990 c write(*,*) 'dd ',dd
991 c else
992 c exit
993 c endif
994 c endif
995 !
996 ! compute xmc=c:xm
997 !
998  do i=1,18
999  if(iorien.gt.0) then
1000  xmc(1,i)=elas(1)*xm(1,i)+elas(2)*xm(2,i)+
1001  & elas(4)*xm(3,i)+2.d0*(elas(7)*xm(4,i)+
1002  & elas(11)*xm(5,i)+elas(16)*xm(6,i))
1003  xmc(2,i)=elas(2)*xm(1,i)+elas(3)*xm(2,i)+
1004  & elas(5)*xm(3,i)+2.d0*(elas(8)*xm(4,i)+
1005  & elas(12)*xm(5,i)+elas(17)*xm(6,i))
1006  xmc(3,i)=elas(4)*xm(1,i)+elas(5)*xm(2,i)+
1007  & elas(6)*xm(3,i)+2.d0*(elas(9)*xm(4,i)+
1008  & elas(13)*xm(5,i)+elas(18)*xm(6,i))
1009  xmc(4,i)=elas(7)*xm(1,i)+elas(8)*xm(2,i)+
1010  & elas(9)*xm(3,i)+2.d0*(elas(10)*xm(4,i)+
1011  & elas(14)*xm(5,i)+elas(19)*xm(6,i))
1012  xmc(5,i)=elas(11)*xm(1,i)+elas(12)*xm(2,i)+
1013  & elas(13)*xm(3,i)+2.d0*(elas(14)*xm(4,i)+
1014  & elas(15)*xm(5,i)+elas(20)*xm(6,i))
1015  xmc(6,i)=elas(16)*xm(1,i)+elas(17)*xm(2,i)+
1016  & elas(18)*xm(3,i)+2.d0*(elas(19)*xm(4,i)+
1017  & elas(20)*xm(5,i)+elas(21)*xm(6,i))
1018  else
1019  xmc(1,i)=c1111*xm(1,i)+c1122*(xm(2,i)+xm(3,i))
1020  xmc(2,i)=c1111*xm(2,i)+c1122*(xm(1,i)+xm(3,i))
1021  xmc(3,i)=c1111*xm(3,i)+c1122*(xm(1,i)+xm(2,i))
1022  xmc(4,i)=2.d0*c1212*xm(4,i)
1023  xmc(5,i)=2.d0*c1212*xm(5,i)
1024  xmc(6,i)=2.d0*c1212*xm(6,i)
1025  endif
1026  enddo
1027 !
1028  neq=18
1029 !
1030 ! filling the LHS
1031 !
1032  do i=1,18
1033  do j=1,18
1034  if(i.ne.j) then
1035  gl(i,j)=(xm(1,i)*xmc(1,j)+
1036  & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0*
1037  & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+
1038  & xm(6,i)*xmc(6,j)))
1039  & *sg(i)*sg(j)
1040  else
1041  gl(i,j)=(xm(1,i)*xmc(1,j)+
1042  & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0*
1043  & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+
1044  & xm(6,i)*xmc(6,j)))
1045  endif
1046  enddo
1047  if(dg(i).gt.1.d-30) then
1048  gl(i,i)=gl(i,i)+
1049  & (dg(i)/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/
1050  & (cn(i)*dtime)
1051  else
1052 !
1053 ! for gamma ein default of 1.d-30 is taken to
1054 ! obtain a finite gradient
1055 !
1056  gl(i,i)=gl(i,i)+
1057  & (1.d-30/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/
1058  & (cn(i)*dtime)
1059  endif
1060  enddo
1061 !
1062 ! filling the RHS
1063 !
1064  do i=1,18
1065 !
1066  gr(i,1)=htri(i)
1067 cccccc & -ck(i)*(dg(i)/dtime)**(1.d0/cn(i))
1068 !
1069  do j=1,6
1070  t(j)=xmc(j,i)*sg(i)
1071  enddo
1072  do j=1,6
1073  gr(i,1)=gr(i,1)-t(j)*r(j)
1074  enddo
1075  gr(i,1)=gr(i,1)
1076  & -t(4)*r(4)-t(5)*r(5)-t(6)*r(6)
1077  enddo
1078 !
1079 ! solve gl*ddg=gr
1080 !
1081  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
1082  if(info.ne.0) then
1083  write(*,*) '*ERROR in umat_single_crystal_creep:'
1084  write(*,*) ' linear equation solver exited'
1085  write(*,*) ' with error: info = ',info
1086  call exit(201)
1087  endif
1088 !
1089  do i=1,18
1090  ddg(i)=gr(i,1)
1091  enddo
1092 !
1093 ! updating the plastic strain
1094 !
1095  do j=1,6
1096  ep(j)=ep(j)+r(j)
1097  do i=1,18
1098  ep(j)=ep(j)+xm(j,i)*sg(i)*ddg(i)
1099  enddo
1100 c write(*,*) 'umat_single_crystal_creep ep',icounter,j,ep(j)
1101  enddo
1102 !
1103 ! updating the consistency parameter
1104 !
1105  do i=1,18
1106  dg(i)=max(dg(i)+ddg(i),0.d0)
1107 c write(*,*) 'umat_single_crystal_creep ddg',icounter,i,ddg(i)
1108  enddo
1109 c write(*,*) 'umat_single_crystal_creep ',icounter,dg(2)
1110 !
1111 ! end of major loop
1112 !
1113  enddo
1114 !
1115 ! inversion of G
1116 !
1117  do i=1,neq
1118  do j=1,neq
1119  gr(i,j)=0.d0
1120  enddo
1121  gr(i,i)=1.d0
1122  enddo
1123  nrhs=neq
1124  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
1125  if(info.ne.0) then
1126  write(*,*) '*ERROR in sc.f: linear equation solver'
1127  write(*,*) ' exited with error: info = ',info
1128  call exit(201)
1129  endif
1130 !
1131 ! storing the stress
1132 !
1133  do i=1,6
1134  stre(i)=stri(i)
1135  enddo
1136 !
1137 ! calculating the tangent stiffness matrix
1138 !
1139  if(icmd.ne.3) then
1140  if(iorien.gt.0) then
1141  ddsdde(1,1)=elas(1)
1142  ddsdde(1,2)=elas(2)
1143  ddsdde(1,3)=elas(4)
1144  ddsdde(1,4)=elas(7)
1145  ddsdde(1,5)=elas(11)
1146  ddsdde(1,6)=elas(16)
1147 c ddsdde(2,1)=elas(2)
1148  ddsdde(2,2)=elas(3)
1149  ddsdde(2,3)=elas(5)
1150  ddsdde(2,4)=elas(8)
1151  ddsdde(2,5)=elas(12)
1152  ddsdde(2,6)=elas(17)
1153 c ddsdde(3,1)=elas(4)
1154 c ddsdde(3,2)=elas(5)
1155  ddsdde(3,3)=elas(6)
1156  ddsdde(3,4)=elas(9)
1157  ddsdde(3,5)=elas(13)
1158  ddsdde(3,6)=elas(18)
1159 c ddsdde(4,1)=elas(7)
1160 c ddsdde(4,2)=elas(8)
1161 c ddsdde(4,3)=elas(9)
1162  ddsdde(4,4)=elas(10)
1163  ddsdde(4,5)=elas(14)
1164  ddsdde(4,6)=elas(19)
1165 c ddsdde(5,1)=elas(11)
1166 c ddsdde(5,2)=elas(12)
1167 c ddsdde(5,3)=elas(13)
1168 c ddsdde(5,4)=elas(14)
1169  ddsdde(5,5)=elas(15)
1170  ddsdde(5,6)=elas(20)
1171 c ddsdde(6,1)=elas(16)
1172 c ddsdde(6,2)=elas(17)
1173 c ddsdde(6,3)=elas(18)
1174 c ddsdde(6,4)=elas(19)
1175 c ddsdde(6,5)=elas(20)
1176  ddsdde(6,6)=elas(21)
1177  else
1178  do i=1,6
1179 c do j=1,6
1180  do j=i,6
1181  ddsdde(i,j)=0.d0
1182  enddo
1183  enddo
1184  do i=1,3
1185  ddsdde(i,i)=c1111
1186  enddo
1187  do i=1,3
1188  do j=i+1,3
1189  ddsdde(i,j)=c1122
1190  enddo
1191 c do j=1,i-1
1192 c ddsdde(i,j)=c1122
1193 c enddo
1194  ddsdde(i+3,i+3)=c1212
1195  enddo
1196  endif
1197  do i=1,18
1198  do j=1,18
1199  do k=1,6
1200 c do l=1,6
1201  do l=k,6
1202  ddsdde(k,l)=ddsdde(k,l)-
1203  & gr(i,j)*xmc(k,i)*sg(i)*xmc(l,j)*sg(j)
1204  enddo
1205  enddo
1206  enddo
1207  enddo
1208 !
1209 ! symmatrizing the stiffness matrix: to be changed: the
1210 ! matrix is symmetric!
1211 !
1212  stiff(1)=ddsdde(1,1)
1213  stiff(2)=ddsdde(1,2)
1214  stiff(3)=ddsdde(2,2)
1215  stiff(4)=ddsdde(1,3)
1216  stiff(5)=ddsdde(2,3)
1217  stiff(6)=ddsdde(3,3)
1218  stiff(7)=ddsdde(1,4)
1219  stiff(8)=ddsdde(2,4)
1220  stiff(9)=ddsdde(3,4)
1221  stiff(10)=ddsdde(4,4)
1222  stiff(11)=ddsdde(1,5)
1223  stiff(12)=ddsdde(2,5)
1224  stiff(13)=ddsdde(3,5)
1225  stiff(14)=ddsdde(4,5)
1226  stiff(15)=ddsdde(5,5)
1227  stiff(16)=ddsdde(1,6)
1228  stiff(17)=ddsdde(2,6)
1229  stiff(18)=ddsdde(3,6)
1230  stiff(19)=ddsdde(4,6)
1231  stiff(20)=ddsdde(5,6)
1232  stiff(21)=ddsdde(6,6)
1233 !
1234  endif
1235 !
1236 ! updating the state variables
1237 !
1238  do i=1,6
1239  xstate(i,iint,iel)=ep(i)
1240  enddo
1241  do i=7,24
1242  xstate(i,iint,iel)=xstateini(i,iint,iel)+dg(i-6)
1243  enddo
1244 !
1245  return
max
#define max(a, b)
Definition: cascade.c:32
dgetrs
subroutine dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:461
dgesv
subroutine dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:58
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)