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

Go to the source code of this file.

Functions/Subroutines

subroutine umat_single_crystal (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)
 

Function/Subroutine Documentation

◆ umat_single_crystal()

subroutine umat_single_crystal ( 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 
)
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 active(18),convergence,creep
119 !
120  character*80 amat
121 !
122  integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien
123 !
124  integer index(18),i,j,k,l,ipiv(18),info,ichange,neq,lda,ldb,
125  & nrhs,iplas,icounter
126 !
127  real*8 ep0(6),al10(18),al20(18),dg0(18),ep(6),al1(18),
128  & al2(18),dg(18),ddg(18),xm(6,18),h(18,18),ck(18),cn(18),
129  & c(18),d(18),phi(18),delta(18),r0(18),q(18),b(18),cphi(18),
130  & q1(18),q2(18),stri(6),htri(18),sg(18),r(42),xmc(6,18),aux(18),
131  & t(42),gl(18,18),gr(18,18),ee(6),c1111,c1122,c1212,dd,
132  & skl(3,3),xmtran(3,3),ddsdde(6,6),xx(6,18)
133 !
134  real*8 elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
135  & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
136  & elas(21),time,ttime
137 !
138  real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
139 !
140  save ep0,al10,al20,dg0,xx,h
141 !
142 !
143 ! crystallographic slip planes:
144 !
145 ! 1. n=1,1,1 l=1,-1,0
146 ! 2. n=1,1,1 l=1,0,-1
147 ! 3. n=1,1,1 l=0,1,-1
148 ! 4. n=1,-1,1 l=0,1,1
149 ! 5. n=1,-1,1 l=1,0,-1
150 ! 6. n=1,-1,1 l=1,1,0
151 ! 7. n=1,-1,-1 l=0,1,-1
152 ! 8. n=1,-1,-1 l=1,0,1
153 ! 9. n=1,-1,-1 l=1,1,0
154 ! 10. n=1,1,-1 l=0,1,1
155 ! 11. n=1,1,-1 l=1,0,1
156 ! 12. n=1,1,-1 l=1,-1,0
157 ! 13. n=1,0,0 l=0,1,1
158 ! 14. n=1,0,0 l=0,1,-1
159 ! 15. n=0,1,0 l=1,0,1
160 ! 16. n=0,1,0 l=1,0,-1
161 ! 17. n=0,0,1 l=1,1,0
162 ! 18. n=0,0,1 l=1,-1,0
163 !
164  xm=reshape((
165  & /0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
166  & 0.0000000000000e+00, 0.2041241452319e+00,-0.2041241452319e+00,
167  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
168  & 0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
169  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
170  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
171  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
172  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
173  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
174  & -0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
175  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
176  & 0.0000000000000e+00, 0.2041241452319e+00, 0.2041241452319e+00,
177  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
178  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
179  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
180  & -0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
181  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
182  & 0.0000000000000e+00,-0.2041241452319e+00,-0.2041241452319e+00,
183  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
184  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
185  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
186  & 0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
187  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
188  & 0.0000000000000e+00,-0.2041241452319e+00, 0.2041241452319e+00,
189  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
190  & 0.3535533905933e+00, 0.3535533905933e+00, 0.0000000000000e+00,
191  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
192  & 0.3535533905933e+00,-0.3535533905933e+00, 0.0000000000000e+00,
193  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
194  & 0.3535533905933e+00, 0.0000000000000e+00, 0.3535533905933e+00,
195  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
196  & 0.3535533905933e+00, 0.0000000000000e+00,-0.3535533905933e+00,
197  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
198  & 0.0000000000000e+00, 0.3535533905933e+00, 0.3535533905933e+00,
199  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
200  & 0.0000000000000e+00, 0.3535533905933e+00,-0.3535533905933e+00/
201  & ),(/6,18/))
202 !
203  xx=reshape((
204  & /0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
205  & 0.0000000000000e+00, 0.2041241452319e+00,-0.2041241452319e+00,
206  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
207  & 0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
208  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
209  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
210  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
211  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
212  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
213  & -0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
214  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
215  & 0.0000000000000e+00, 0.2041241452319e+00, 0.2041241452319e+00,
216  & 0.0000000000000e+00,-0.4082482904639e+00, 0.4082482904639e+00,
217  & 0.2041241452319e+00,-0.2041241452319e+00, 0.0000000000000e+00,
218  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
219  & -0.2041241452319e+00, 0.0000000000000e+00,-0.2041241452319e+00,
220  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
221  & 0.0000000000000e+00,-0.2041241452319e+00,-0.2041241452319e+00,
222  & 0.0000000000000e+00, 0.4082482904639e+00,-0.4082482904639e+00,
223  & 0.2041241452319e+00, 0.2041241452319e+00, 0.0000000000000e+00,
224  & 0.4082482904639e+00, 0.0000000000000e+00,-0.4082482904639e+00,
225  & 0.2041241452319e+00, 0.0000000000000e+00, 0.2041241452319e+00,
226  & 0.4082482904639e+00,-0.4082482904639e+00, 0.0000000000000e+00,
227  & 0.0000000000000e+00,-0.2041241452319e+00, 0.2041241452319e+00,
228  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
229  & 0.3535533905933e+00, 0.3535533905933e+00, 0.0000000000000e+00,
230  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
231  & 0.3535533905933e+00,-0.3535533905933e+00, 0.0000000000000e+00,
232  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
233  & 0.3535533905933e+00, 0.0000000000000e+00, 0.3535533905933e+00,
234  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
235  & 0.3535533905933e+00, 0.0000000000000e+00,-0.3535533905933e+00,
236  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
237  & 0.0000000000000e+00, 0.3535533905933e+00, 0.3535533905933e+00,
238  & 0.0000000000000e+00, 0.0000000000000e+00, 0.0000000000000e+00,
239  & 0.0000000000000e+00, 0.3535533905933e+00,-0.3535533905933e+00/
240  & ),(/6,18/))
241 !
242  h=reshape((
243  & /0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
244  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00,
245  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00, 0.1e+01,-0.1e+00,
246  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
247  & -0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
248  & 0.0e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,
249  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00,
250  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,
251  & -0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
252  & -0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
253  & 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,
254  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
255  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,
256  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00,
257  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00,
258  & 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
259  & -0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
260  & -0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
261  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
262  & 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00,
263  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,
264  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,
265  & -0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
266  & 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
267  & -0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00, 0.0e+00,
268  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,
269  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
270  & -0.1e+00, 0.1e+01,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
271  & 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
272  & -0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,
273  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
274  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
275  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.1e+01,-0.1e+00,-0.1e+00,
276  & -0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
277  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
278  & 0.0e+00,-0.1e+00, 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
279  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
280  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,
281  & 0.1e+01,-0.1e+00,-0.1e+00,-0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00,
282  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
283  & 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00,
284  & -0.1e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
285  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,
286  & -0.1e+00,-0.1e+00,-0.1e+00, 0.1e+01,-0.1e+00, 0.0e+00, 0.0e+00,
287  & 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00,
288  & 0.0e+00, 0.0e+00, 0.0e+00,-0.1e+00,-0.1e+00,-0.1e+00,-0.1e+00,
289  & -0.1e+00, 0.1e+01/),(/18,18/))
290 !
291 ! elastic constants
292 !
293  c1111=elconloc(1)
294  c1122=elconloc(2)
295  c1212=elconloc(3)
296 !
297  if(iorien.gt.0) then
298  call transformatrix(orab(1,iorien),pgauss,skl)
299  do k=1,18
300  do i=1,3
301  do j=i,3
302  xmtran(i,j)=skl(i,1)*skl(j,1)*xx(1,k)+
303  & skl(i,2)*skl(j,2)*xx(2,k)+
304  & skl(i,3)*skl(j,3)*xx(3,k)+
305  & (skl(i,1)*skl(j,2)+
306  & skl(i,2)*skl(j,1))*xx(4,k)+
307  & (skl(i,1)*skl(j,3)+
308  & skl(i,3)*skl(j,1))*xx(5,k)+
309  & (skl(i,2)*skl(j,3)+
310  & skl(i,3)*skl(j,2))*xx(6,k)
311  enddo
312  enddo
313  xm(1,k)=xmtran(1,1)
314  xm(2,k)=xmtran(2,2)
315  xm(3,k)=xmtran(3,3)
316  xm(4,k)=xmtran(1,2)
317  xm(5,k)=xmtran(1,3)
318  xm(6,k)=xmtran(2,3)
319  enddo
320 !
321  elas( 1)=
322  & skl(1,1)*skl(1,1)*skl(1,1)*skl(1,1)*c1111+
323  & skl(1,1)*skl(1,1)*skl(1,2)*skl(1,2)*c1122+
324  & skl(1,1)*skl(1,1)*skl(1,3)*skl(1,3)*c1122+
325  & skl(1,1)*skl(1,2)*skl(1,1)*skl(1,2)*c1212+
326  & skl(1,1)*skl(1,2)*skl(1,2)*skl(1,1)*c1212+
327  & skl(1,1)*skl(1,3)*skl(1,1)*skl(1,3)*c1212+
328  & skl(1,1)*skl(1,3)*skl(1,3)*skl(1,1)*c1212+
329  & skl(1,2)*skl(1,1)*skl(1,1)*skl(1,2)*c1212+
330  & skl(1,2)*skl(1,1)*skl(1,2)*skl(1,1)*c1212+
331  & skl(1,2)*skl(1,2)*skl(1,1)*skl(1,1)*c1122+
332  & skl(1,2)*skl(1,2)*skl(1,2)*skl(1,2)*c1111+
333  & skl(1,2)*skl(1,2)*skl(1,3)*skl(1,3)*c1122+
334  & skl(1,2)*skl(1,3)*skl(1,2)*skl(1,3)*c1212+
335  & skl(1,2)*skl(1,3)*skl(1,3)*skl(1,2)*c1212+
336  & skl(1,3)*skl(1,1)*skl(1,1)*skl(1,3)*c1212+
337  & skl(1,3)*skl(1,1)*skl(1,3)*skl(1,1)*c1212+
338  & skl(1,3)*skl(1,2)*skl(1,2)*skl(1,3)*c1212+
339  & skl(1,3)*skl(1,2)*skl(1,3)*skl(1,2)*c1212+
340  & skl(1,3)*skl(1,3)*skl(1,1)*skl(1,1)*c1122+
341  & skl(1,3)*skl(1,3)*skl(1,2)*skl(1,2)*c1122+
342  & skl(1,3)*skl(1,3)*skl(1,3)*skl(1,3)*c1111
343  elas( 2)=
344  & skl(1,1)*skl(1,1)*skl(2,1)*skl(2,1)*c1111+
345  & skl(1,1)*skl(1,1)*skl(2,2)*skl(2,2)*c1122+
346  & skl(1,1)*skl(1,1)*skl(2,3)*skl(2,3)*c1122+
347  & skl(1,1)*skl(1,2)*skl(2,1)*skl(2,2)*c1212+
348  & skl(1,1)*skl(1,2)*skl(2,2)*skl(2,1)*c1212+
349  & skl(1,1)*skl(1,3)*skl(2,1)*skl(2,3)*c1212+
350  & skl(1,1)*skl(1,3)*skl(2,3)*skl(2,1)*c1212+
351  & skl(1,2)*skl(1,1)*skl(2,1)*skl(2,2)*c1212+
352  & skl(1,2)*skl(1,1)*skl(2,2)*skl(2,1)*c1212+
353  & skl(1,2)*skl(1,2)*skl(2,1)*skl(2,1)*c1122+
354  & skl(1,2)*skl(1,2)*skl(2,2)*skl(2,2)*c1111+
355  & skl(1,2)*skl(1,2)*skl(2,3)*skl(2,3)*c1122+
356  & skl(1,2)*skl(1,3)*skl(2,2)*skl(2,3)*c1212+
357  & skl(1,2)*skl(1,3)*skl(2,3)*skl(2,2)*c1212+
358  & skl(1,3)*skl(1,1)*skl(2,1)*skl(2,3)*c1212+
359  & skl(1,3)*skl(1,1)*skl(2,3)*skl(2,1)*c1212+
360  & skl(1,3)*skl(1,2)*skl(2,2)*skl(2,3)*c1212+
361  & skl(1,3)*skl(1,2)*skl(2,3)*skl(2,2)*c1212+
362  & skl(1,3)*skl(1,3)*skl(2,1)*skl(2,1)*c1122+
363  & skl(1,3)*skl(1,3)*skl(2,2)*skl(2,2)*c1122+
364  & skl(1,3)*skl(1,3)*skl(2,3)*skl(2,3)*c1111
365  elas( 3)=
366  & skl(2,1)*skl(2,1)*skl(2,1)*skl(2,1)*c1111+
367  & skl(2,1)*skl(2,1)*skl(2,2)*skl(2,2)*c1122+
368  & skl(2,1)*skl(2,1)*skl(2,3)*skl(2,3)*c1122+
369  & skl(2,1)*skl(2,2)*skl(2,1)*skl(2,2)*c1212+
370  & skl(2,1)*skl(2,2)*skl(2,2)*skl(2,1)*c1212+
371  & skl(2,1)*skl(2,3)*skl(2,1)*skl(2,3)*c1212+
372  & skl(2,1)*skl(2,3)*skl(2,3)*skl(2,1)*c1212+
373  & skl(2,2)*skl(2,1)*skl(2,1)*skl(2,2)*c1212+
374  & skl(2,2)*skl(2,1)*skl(2,2)*skl(2,1)*c1212+
375  & skl(2,2)*skl(2,2)*skl(2,1)*skl(2,1)*c1122+
376  & skl(2,2)*skl(2,2)*skl(2,2)*skl(2,2)*c1111+
377  & skl(2,2)*skl(2,2)*skl(2,3)*skl(2,3)*c1122+
378  & skl(2,2)*skl(2,3)*skl(2,2)*skl(2,3)*c1212+
379  & skl(2,2)*skl(2,3)*skl(2,3)*skl(2,2)*c1212+
380  & skl(2,3)*skl(2,1)*skl(2,1)*skl(2,3)*c1212+
381  & skl(2,3)*skl(2,1)*skl(2,3)*skl(2,1)*c1212+
382  & skl(2,3)*skl(2,2)*skl(2,2)*skl(2,3)*c1212+
383  & skl(2,3)*skl(2,2)*skl(2,3)*skl(2,2)*c1212+
384  & skl(2,3)*skl(2,3)*skl(2,1)*skl(2,1)*c1122+
385  & skl(2,3)*skl(2,3)*skl(2,2)*skl(2,2)*c1122+
386  & skl(2,3)*skl(2,3)*skl(2,3)*skl(2,3)*c1111
387  elas( 4)=
388  & skl(1,1)*skl(1,1)*skl(3,1)*skl(3,1)*c1111+
389  & skl(1,1)*skl(1,1)*skl(3,2)*skl(3,2)*c1122+
390  & skl(1,1)*skl(1,1)*skl(3,3)*skl(3,3)*c1122+
391  & skl(1,1)*skl(1,2)*skl(3,1)*skl(3,2)*c1212+
392  & skl(1,1)*skl(1,2)*skl(3,2)*skl(3,1)*c1212+
393  & skl(1,1)*skl(1,3)*skl(3,1)*skl(3,3)*c1212+
394  & skl(1,1)*skl(1,3)*skl(3,3)*skl(3,1)*c1212+
395  & skl(1,2)*skl(1,1)*skl(3,1)*skl(3,2)*c1212+
396  & skl(1,2)*skl(1,1)*skl(3,2)*skl(3,1)*c1212+
397  & skl(1,2)*skl(1,2)*skl(3,1)*skl(3,1)*c1122+
398  & skl(1,2)*skl(1,2)*skl(3,2)*skl(3,2)*c1111+
399  & skl(1,2)*skl(1,2)*skl(3,3)*skl(3,3)*c1122+
400  & skl(1,2)*skl(1,3)*skl(3,2)*skl(3,3)*c1212+
401  & skl(1,2)*skl(1,3)*skl(3,3)*skl(3,2)*c1212+
402  & skl(1,3)*skl(1,1)*skl(3,1)*skl(3,3)*c1212+
403  & skl(1,3)*skl(1,1)*skl(3,3)*skl(3,1)*c1212+
404  & skl(1,3)*skl(1,2)*skl(3,2)*skl(3,3)*c1212+
405  & skl(1,3)*skl(1,2)*skl(3,3)*skl(3,2)*c1212+
406  & skl(1,3)*skl(1,3)*skl(3,1)*skl(3,1)*c1122+
407  & skl(1,3)*skl(1,3)*skl(3,2)*skl(3,2)*c1122+
408  & skl(1,3)*skl(1,3)*skl(3,3)*skl(3,3)*c1111
409  elas( 5)=
410  & skl(2,1)*skl(2,1)*skl(3,1)*skl(3,1)*c1111+
411  & skl(2,1)*skl(2,1)*skl(3,2)*skl(3,2)*c1122+
412  & skl(2,1)*skl(2,1)*skl(3,3)*skl(3,3)*c1122+
413  & skl(2,1)*skl(2,2)*skl(3,1)*skl(3,2)*c1212+
414  & skl(2,1)*skl(2,2)*skl(3,2)*skl(3,1)*c1212+
415  & skl(2,1)*skl(2,3)*skl(3,1)*skl(3,3)*c1212+
416  & skl(2,1)*skl(2,3)*skl(3,3)*skl(3,1)*c1212+
417  & skl(2,2)*skl(2,1)*skl(3,1)*skl(3,2)*c1212+
418  & skl(2,2)*skl(2,1)*skl(3,2)*skl(3,1)*c1212+
419  & skl(2,2)*skl(2,2)*skl(3,1)*skl(3,1)*c1122+
420  & skl(2,2)*skl(2,2)*skl(3,2)*skl(3,2)*c1111+
421  & skl(2,2)*skl(2,2)*skl(3,3)*skl(3,3)*c1122+
422  & skl(2,2)*skl(2,3)*skl(3,2)*skl(3,3)*c1212+
423  & skl(2,2)*skl(2,3)*skl(3,3)*skl(3,2)*c1212+
424  & skl(2,3)*skl(2,1)*skl(3,1)*skl(3,3)*c1212+
425  & skl(2,3)*skl(2,1)*skl(3,3)*skl(3,1)*c1212+
426  & skl(2,3)*skl(2,2)*skl(3,2)*skl(3,3)*c1212+
427  & skl(2,3)*skl(2,2)*skl(3,3)*skl(3,2)*c1212+
428  & skl(2,3)*skl(2,3)*skl(3,1)*skl(3,1)*c1122+
429  & skl(2,3)*skl(2,3)*skl(3,2)*skl(3,2)*c1122+
430  & skl(2,3)*skl(2,3)*skl(3,3)*skl(3,3)*c1111
431  elas( 6)=
432  & skl(3,1)*skl(3,1)*skl(3,1)*skl(3,1)*c1111+
433  & skl(3,1)*skl(3,1)*skl(3,2)*skl(3,2)*c1122+
434  & skl(3,1)*skl(3,1)*skl(3,3)*skl(3,3)*c1122+
435  & skl(3,1)*skl(3,2)*skl(3,1)*skl(3,2)*c1212+
436  & skl(3,1)*skl(3,2)*skl(3,2)*skl(3,1)*c1212+
437  & skl(3,1)*skl(3,3)*skl(3,1)*skl(3,3)*c1212+
438  & skl(3,1)*skl(3,3)*skl(3,3)*skl(3,1)*c1212+
439  & skl(3,2)*skl(3,1)*skl(3,1)*skl(3,2)*c1212+
440  & skl(3,2)*skl(3,1)*skl(3,2)*skl(3,1)*c1212+
441  & skl(3,2)*skl(3,2)*skl(3,1)*skl(3,1)*c1122+
442  & skl(3,2)*skl(3,2)*skl(3,2)*skl(3,2)*c1111+
443  & skl(3,2)*skl(3,2)*skl(3,3)*skl(3,3)*c1122+
444  & skl(3,2)*skl(3,3)*skl(3,2)*skl(3,3)*c1212+
445  & skl(3,2)*skl(3,3)*skl(3,3)*skl(3,2)*c1212+
446  & skl(3,3)*skl(3,1)*skl(3,1)*skl(3,3)*c1212+
447  & skl(3,3)*skl(3,1)*skl(3,3)*skl(3,1)*c1212+
448  & skl(3,3)*skl(3,2)*skl(3,2)*skl(3,3)*c1212+
449  & skl(3,3)*skl(3,2)*skl(3,3)*skl(3,2)*c1212+
450  & skl(3,3)*skl(3,3)*skl(3,1)*skl(3,1)*c1122+
451  & skl(3,3)*skl(3,3)*skl(3,2)*skl(3,2)*c1122+
452  & skl(3,3)*skl(3,3)*skl(3,3)*skl(3,3)*c1111
453  elas( 7)=
454  & skl(1,1)*skl(1,1)*skl(1,1)*skl(2,1)*c1111+
455  & skl(1,1)*skl(1,1)*skl(1,2)*skl(2,2)*c1122+
456  & skl(1,1)*skl(1,1)*skl(1,3)*skl(2,3)*c1122+
457  & skl(1,1)*skl(1,2)*skl(1,1)*skl(2,2)*c1212+
458  & skl(1,1)*skl(1,2)*skl(1,2)*skl(2,1)*c1212+
459  & skl(1,1)*skl(1,3)*skl(1,1)*skl(2,3)*c1212+
460  & skl(1,1)*skl(1,3)*skl(1,3)*skl(2,1)*c1212+
461  & skl(1,2)*skl(1,1)*skl(1,1)*skl(2,2)*c1212+
462  & skl(1,2)*skl(1,1)*skl(1,2)*skl(2,1)*c1212+
463  & skl(1,2)*skl(1,2)*skl(1,1)*skl(2,1)*c1122+
464  & skl(1,2)*skl(1,2)*skl(1,2)*skl(2,2)*c1111+
465  & skl(1,2)*skl(1,2)*skl(1,3)*skl(2,3)*c1122+
466  & skl(1,2)*skl(1,3)*skl(1,2)*skl(2,3)*c1212+
467  & skl(1,2)*skl(1,3)*skl(1,3)*skl(2,2)*c1212+
468  & skl(1,3)*skl(1,1)*skl(1,1)*skl(2,3)*c1212+
469  & skl(1,3)*skl(1,1)*skl(1,3)*skl(2,1)*c1212+
470  & skl(1,3)*skl(1,2)*skl(1,2)*skl(2,3)*c1212+
471  & skl(1,3)*skl(1,2)*skl(1,3)*skl(2,2)*c1212+
472  & skl(1,3)*skl(1,3)*skl(1,1)*skl(2,1)*c1122+
473  & skl(1,3)*skl(1,3)*skl(1,2)*skl(2,2)*c1122+
474  & skl(1,3)*skl(1,3)*skl(1,3)*skl(2,3)*c1111
475  elas( 8)=
476  & skl(2,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+
477  & skl(2,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+
478  & skl(2,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+
479  & skl(2,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+
480  & skl(2,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+
481  & skl(2,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+
482  & skl(2,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+
483  & skl(2,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+
484  & skl(2,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+
485  & skl(2,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+
486  & skl(2,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+
487  & skl(2,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+
488  & skl(2,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+
489  & skl(2,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+
490  & skl(2,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+
491  & skl(2,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+
492  & skl(2,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+
493  & skl(2,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+
494  & skl(2,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+
495  & skl(2,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+
496  & skl(2,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111
497  elas( 9)=
498  & skl(3,1)*skl(3,1)*skl(1,1)*skl(2,1)*c1111+
499  & skl(3,1)*skl(3,1)*skl(1,2)*skl(2,2)*c1122+
500  & skl(3,1)*skl(3,1)*skl(1,3)*skl(2,3)*c1122+
501  & skl(3,1)*skl(3,2)*skl(1,1)*skl(2,2)*c1212+
502  & skl(3,1)*skl(3,2)*skl(1,2)*skl(2,1)*c1212+
503  & skl(3,1)*skl(3,3)*skl(1,1)*skl(2,3)*c1212+
504  & skl(3,1)*skl(3,3)*skl(1,3)*skl(2,1)*c1212+
505  & skl(3,2)*skl(3,1)*skl(1,1)*skl(2,2)*c1212+
506  & skl(3,2)*skl(3,1)*skl(1,2)*skl(2,1)*c1212+
507  & skl(3,2)*skl(3,2)*skl(1,1)*skl(2,1)*c1122+
508  & skl(3,2)*skl(3,2)*skl(1,2)*skl(2,2)*c1111+
509  & skl(3,2)*skl(3,2)*skl(1,3)*skl(2,3)*c1122+
510  & skl(3,2)*skl(3,3)*skl(1,2)*skl(2,3)*c1212+
511  & skl(3,2)*skl(3,3)*skl(1,3)*skl(2,2)*c1212+
512  & skl(3,3)*skl(3,1)*skl(1,1)*skl(2,3)*c1212+
513  & skl(3,3)*skl(3,1)*skl(1,3)*skl(2,1)*c1212+
514  & skl(3,3)*skl(3,2)*skl(1,2)*skl(2,3)*c1212+
515  & skl(3,3)*skl(3,2)*skl(1,3)*skl(2,2)*c1212+
516  & skl(3,3)*skl(3,3)*skl(1,1)*skl(2,1)*c1122+
517  & skl(3,3)*skl(3,3)*skl(1,2)*skl(2,2)*c1122+
518  & skl(3,3)*skl(3,3)*skl(1,3)*skl(2,3)*c1111
519  elas(10)=
520  & skl(1,1)*skl(2,1)*skl(1,1)*skl(2,1)*c1111+
521  & skl(1,1)*skl(2,1)*skl(1,2)*skl(2,2)*c1122+
522  & skl(1,1)*skl(2,1)*skl(1,3)*skl(2,3)*c1122+
523  & skl(1,1)*skl(2,2)*skl(1,1)*skl(2,2)*c1212+
524  & skl(1,1)*skl(2,2)*skl(1,2)*skl(2,1)*c1212+
525  & skl(1,1)*skl(2,3)*skl(1,1)*skl(2,3)*c1212+
526  & skl(1,1)*skl(2,3)*skl(1,3)*skl(2,1)*c1212+
527  & skl(1,2)*skl(2,1)*skl(1,1)*skl(2,2)*c1212+
528  & skl(1,2)*skl(2,1)*skl(1,2)*skl(2,1)*c1212+
529  & skl(1,2)*skl(2,2)*skl(1,1)*skl(2,1)*c1122+
530  & skl(1,2)*skl(2,2)*skl(1,2)*skl(2,2)*c1111+
531  & skl(1,2)*skl(2,2)*skl(1,3)*skl(2,3)*c1122+
532  & skl(1,2)*skl(2,3)*skl(1,2)*skl(2,3)*c1212+
533  & skl(1,2)*skl(2,3)*skl(1,3)*skl(2,2)*c1212+
534  & skl(1,3)*skl(2,1)*skl(1,1)*skl(2,3)*c1212+
535  & skl(1,3)*skl(2,1)*skl(1,3)*skl(2,1)*c1212+
536  & skl(1,3)*skl(2,2)*skl(1,2)*skl(2,3)*c1212+
537  & skl(1,3)*skl(2,2)*skl(1,3)*skl(2,2)*c1212+
538  & skl(1,3)*skl(2,3)*skl(1,1)*skl(2,1)*c1122+
539  & skl(1,3)*skl(2,3)*skl(1,2)*skl(2,2)*c1122+
540  & skl(1,3)*skl(2,3)*skl(1,3)*skl(2,3)*c1111
541  elas(11)=
542  & skl(1,1)*skl(1,1)*skl(1,1)*skl(3,1)*c1111+
543  & skl(1,1)*skl(1,1)*skl(1,2)*skl(3,2)*c1122+
544  & skl(1,1)*skl(1,1)*skl(1,3)*skl(3,3)*c1122+
545  & skl(1,1)*skl(1,2)*skl(1,1)*skl(3,2)*c1212+
546  & skl(1,1)*skl(1,2)*skl(1,2)*skl(3,1)*c1212+
547  & skl(1,1)*skl(1,3)*skl(1,1)*skl(3,3)*c1212+
548  & skl(1,1)*skl(1,3)*skl(1,3)*skl(3,1)*c1212+
549  & skl(1,2)*skl(1,1)*skl(1,1)*skl(3,2)*c1212+
550  & skl(1,2)*skl(1,1)*skl(1,2)*skl(3,1)*c1212+
551  & skl(1,2)*skl(1,2)*skl(1,1)*skl(3,1)*c1122+
552  & skl(1,2)*skl(1,2)*skl(1,2)*skl(3,2)*c1111+
553  & skl(1,2)*skl(1,2)*skl(1,3)*skl(3,3)*c1122+
554  & skl(1,2)*skl(1,3)*skl(1,2)*skl(3,3)*c1212+
555  & skl(1,2)*skl(1,3)*skl(1,3)*skl(3,2)*c1212+
556  & skl(1,3)*skl(1,1)*skl(1,1)*skl(3,3)*c1212+
557  & skl(1,3)*skl(1,1)*skl(1,3)*skl(3,1)*c1212+
558  & skl(1,3)*skl(1,2)*skl(1,2)*skl(3,3)*c1212+
559  & skl(1,3)*skl(1,2)*skl(1,3)*skl(3,2)*c1212+
560  & skl(1,3)*skl(1,3)*skl(1,1)*skl(3,1)*c1122+
561  & skl(1,3)*skl(1,3)*skl(1,2)*skl(3,2)*c1122+
562  & skl(1,3)*skl(1,3)*skl(1,3)*skl(3,3)*c1111
563  elas(12)=
564  & skl(2,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+
565  & skl(2,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+
566  & skl(2,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+
567  & skl(2,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+
568  & skl(2,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+
569  & skl(2,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+
570  & skl(2,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+
571  & skl(2,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+
572  & skl(2,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+
573  & skl(2,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+
574  & skl(2,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+
575  & skl(2,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+
576  & skl(2,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+
577  & skl(2,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+
578  & skl(2,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+
579  & skl(2,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+
580  & skl(2,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+
581  & skl(2,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+
582  & skl(2,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+
583  & skl(2,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+
584  & skl(2,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111
585  elas(13)=
586  & skl(3,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+
587  & skl(3,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+
588  & skl(3,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+
589  & skl(3,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+
590  & skl(3,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+
591  & skl(3,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+
592  & skl(3,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+
593  & skl(3,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+
594  & skl(3,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+
595  & skl(3,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+
596  & skl(3,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+
597  & skl(3,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+
598  & skl(3,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+
599  & skl(3,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+
600  & skl(3,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+
601  & skl(3,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+
602  & skl(3,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+
603  & skl(3,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+
604  & skl(3,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+
605  & skl(3,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+
606  & skl(3,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111
607  elas(14)=
608  & skl(1,1)*skl(2,1)*skl(1,1)*skl(3,1)*c1111+
609  & skl(1,1)*skl(2,1)*skl(1,2)*skl(3,2)*c1122+
610  & skl(1,1)*skl(2,1)*skl(1,3)*skl(3,3)*c1122+
611  & skl(1,1)*skl(2,2)*skl(1,1)*skl(3,2)*c1212+
612  & skl(1,1)*skl(2,2)*skl(1,2)*skl(3,1)*c1212+
613  & skl(1,1)*skl(2,3)*skl(1,1)*skl(3,3)*c1212+
614  & skl(1,1)*skl(2,3)*skl(1,3)*skl(3,1)*c1212+
615  & skl(1,2)*skl(2,1)*skl(1,1)*skl(3,2)*c1212+
616  & skl(1,2)*skl(2,1)*skl(1,2)*skl(3,1)*c1212+
617  & skl(1,2)*skl(2,2)*skl(1,1)*skl(3,1)*c1122+
618  & skl(1,2)*skl(2,2)*skl(1,2)*skl(3,2)*c1111+
619  & skl(1,2)*skl(2,2)*skl(1,3)*skl(3,3)*c1122+
620  & skl(1,2)*skl(2,3)*skl(1,2)*skl(3,3)*c1212+
621  & skl(1,2)*skl(2,3)*skl(1,3)*skl(3,2)*c1212+
622  & skl(1,3)*skl(2,1)*skl(1,1)*skl(3,3)*c1212+
623  & skl(1,3)*skl(2,1)*skl(1,3)*skl(3,1)*c1212+
624  & skl(1,3)*skl(2,2)*skl(1,2)*skl(3,3)*c1212+
625  & skl(1,3)*skl(2,2)*skl(1,3)*skl(3,2)*c1212+
626  & skl(1,3)*skl(2,3)*skl(1,1)*skl(3,1)*c1122+
627  & skl(1,3)*skl(2,3)*skl(1,2)*skl(3,2)*c1122+
628  & skl(1,3)*skl(2,3)*skl(1,3)*skl(3,3)*c1111
629  elas(15)=
630  & skl(1,1)*skl(3,1)*skl(1,1)*skl(3,1)*c1111+
631  & skl(1,1)*skl(3,1)*skl(1,2)*skl(3,2)*c1122+
632  & skl(1,1)*skl(3,1)*skl(1,3)*skl(3,3)*c1122+
633  & skl(1,1)*skl(3,2)*skl(1,1)*skl(3,2)*c1212+
634  & skl(1,1)*skl(3,2)*skl(1,2)*skl(3,1)*c1212+
635  & skl(1,1)*skl(3,3)*skl(1,1)*skl(3,3)*c1212+
636  & skl(1,1)*skl(3,3)*skl(1,3)*skl(3,1)*c1212+
637  & skl(1,2)*skl(3,1)*skl(1,1)*skl(3,2)*c1212+
638  & skl(1,2)*skl(3,1)*skl(1,2)*skl(3,1)*c1212+
639  & skl(1,2)*skl(3,2)*skl(1,1)*skl(3,1)*c1122+
640  & skl(1,2)*skl(3,2)*skl(1,2)*skl(3,2)*c1111+
641  & skl(1,2)*skl(3,2)*skl(1,3)*skl(3,3)*c1122+
642  & skl(1,2)*skl(3,3)*skl(1,2)*skl(3,3)*c1212+
643  & skl(1,2)*skl(3,3)*skl(1,3)*skl(3,2)*c1212+
644  & skl(1,3)*skl(3,1)*skl(1,1)*skl(3,3)*c1212+
645  & skl(1,3)*skl(3,1)*skl(1,3)*skl(3,1)*c1212+
646  & skl(1,3)*skl(3,2)*skl(1,2)*skl(3,3)*c1212+
647  & skl(1,3)*skl(3,2)*skl(1,3)*skl(3,2)*c1212+
648  & skl(1,3)*skl(3,3)*skl(1,1)*skl(3,1)*c1122+
649  & skl(1,3)*skl(3,3)*skl(1,2)*skl(3,2)*c1122+
650  & skl(1,3)*skl(3,3)*skl(1,3)*skl(3,3)*c1111
651  elas(16)=
652  & skl(1,1)*skl(1,1)*skl(2,1)*skl(3,1)*c1111+
653  & skl(1,1)*skl(1,1)*skl(2,2)*skl(3,2)*c1122+
654  & skl(1,1)*skl(1,1)*skl(2,3)*skl(3,3)*c1122+
655  & skl(1,1)*skl(1,2)*skl(2,1)*skl(3,2)*c1212+
656  & skl(1,1)*skl(1,2)*skl(2,2)*skl(3,1)*c1212+
657  & skl(1,1)*skl(1,3)*skl(2,1)*skl(3,3)*c1212+
658  & skl(1,1)*skl(1,3)*skl(2,3)*skl(3,1)*c1212+
659  & skl(1,2)*skl(1,1)*skl(2,1)*skl(3,2)*c1212+
660  & skl(1,2)*skl(1,1)*skl(2,2)*skl(3,1)*c1212+
661  & skl(1,2)*skl(1,2)*skl(2,1)*skl(3,1)*c1122+
662  & skl(1,2)*skl(1,2)*skl(2,2)*skl(3,2)*c1111+
663  & skl(1,2)*skl(1,2)*skl(2,3)*skl(3,3)*c1122+
664  & skl(1,2)*skl(1,3)*skl(2,2)*skl(3,3)*c1212+
665  & skl(1,2)*skl(1,3)*skl(2,3)*skl(3,2)*c1212+
666  & skl(1,3)*skl(1,1)*skl(2,1)*skl(3,3)*c1212+
667  & skl(1,3)*skl(1,1)*skl(2,3)*skl(3,1)*c1212+
668  & skl(1,3)*skl(1,2)*skl(2,2)*skl(3,3)*c1212+
669  & skl(1,3)*skl(1,2)*skl(2,3)*skl(3,2)*c1212+
670  & skl(1,3)*skl(1,3)*skl(2,1)*skl(3,1)*c1122+
671  & skl(1,3)*skl(1,3)*skl(2,2)*skl(3,2)*c1122+
672  & skl(1,3)*skl(1,3)*skl(2,3)*skl(3,3)*c1111
673  elas(17)=
674  & skl(2,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+
675  & skl(2,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+
676  & skl(2,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+
677  & skl(2,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+
678  & skl(2,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+
679  & skl(2,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+
680  & skl(2,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+
681  & skl(2,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+
682  & skl(2,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+
683  & skl(2,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+
684  & skl(2,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+
685  & skl(2,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+
686  & skl(2,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+
687  & skl(2,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+
688  & skl(2,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+
689  & skl(2,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+
690  & skl(2,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+
691  & skl(2,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+
692  & skl(2,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+
693  & skl(2,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+
694  & skl(2,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111
695  elas(18)=
696  & skl(3,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
697  & skl(3,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
698  & skl(3,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
699  & skl(3,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
700  & skl(3,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
701  & skl(3,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
702  & skl(3,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
703  & skl(3,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
704  & skl(3,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
705  & skl(3,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
706  & skl(3,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
707  & skl(3,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
708  & skl(3,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
709  & skl(3,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
710  & skl(3,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
711  & skl(3,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
712  & skl(3,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
713  & skl(3,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
714  & skl(3,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
715  & skl(3,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
716  & skl(3,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
717  elas(19)=
718  & skl(1,1)*skl(2,1)*skl(2,1)*skl(3,1)*c1111+
719  & skl(1,1)*skl(2,1)*skl(2,2)*skl(3,2)*c1122+
720  & skl(1,1)*skl(2,1)*skl(2,3)*skl(3,3)*c1122+
721  & skl(1,1)*skl(2,2)*skl(2,1)*skl(3,2)*c1212+
722  & skl(1,1)*skl(2,2)*skl(2,2)*skl(3,1)*c1212+
723  & skl(1,1)*skl(2,3)*skl(2,1)*skl(3,3)*c1212+
724  & skl(1,1)*skl(2,3)*skl(2,3)*skl(3,1)*c1212+
725  & skl(1,2)*skl(2,1)*skl(2,1)*skl(3,2)*c1212+
726  & skl(1,2)*skl(2,1)*skl(2,2)*skl(3,1)*c1212+
727  & skl(1,2)*skl(2,2)*skl(2,1)*skl(3,1)*c1122+
728  & skl(1,2)*skl(2,2)*skl(2,2)*skl(3,2)*c1111+
729  & skl(1,2)*skl(2,2)*skl(2,3)*skl(3,3)*c1122+
730  & skl(1,2)*skl(2,3)*skl(2,2)*skl(3,3)*c1212+
731  & skl(1,2)*skl(2,3)*skl(2,3)*skl(3,2)*c1212+
732  & skl(1,3)*skl(2,1)*skl(2,1)*skl(3,3)*c1212+
733  & skl(1,3)*skl(2,1)*skl(2,3)*skl(3,1)*c1212+
734  & skl(1,3)*skl(2,2)*skl(2,2)*skl(3,3)*c1212+
735  & skl(1,3)*skl(2,2)*skl(2,3)*skl(3,2)*c1212+
736  & skl(1,3)*skl(2,3)*skl(2,1)*skl(3,1)*c1122+
737  & skl(1,3)*skl(2,3)*skl(2,2)*skl(3,2)*c1122+
738  & skl(1,3)*skl(2,3)*skl(2,3)*skl(3,3)*c1111
739  elas(20)=
740  & skl(1,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
741  & skl(1,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
742  & skl(1,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
743  & skl(1,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
744  & skl(1,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
745  & skl(1,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
746  & skl(1,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
747  & skl(1,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
748  & skl(1,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
749  & skl(1,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
750  & skl(1,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
751  & skl(1,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
752  & skl(1,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
753  & skl(1,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
754  & skl(1,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
755  & skl(1,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
756  & skl(1,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
757  & skl(1,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
758  & skl(1,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
759  & skl(1,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
760  & skl(1,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
761  elas(21)=
762  & skl(2,1)*skl(3,1)*skl(2,1)*skl(3,1)*c1111+
763  & skl(2,1)*skl(3,1)*skl(2,2)*skl(3,2)*c1122+
764  & skl(2,1)*skl(3,1)*skl(2,3)*skl(3,3)*c1122+
765  & skl(2,1)*skl(3,2)*skl(2,1)*skl(3,2)*c1212+
766  & skl(2,1)*skl(3,2)*skl(2,2)*skl(3,1)*c1212+
767  & skl(2,1)*skl(3,3)*skl(2,1)*skl(3,3)*c1212+
768  & skl(2,1)*skl(3,3)*skl(2,3)*skl(3,1)*c1212+
769  & skl(2,2)*skl(3,1)*skl(2,1)*skl(3,2)*c1212+
770  & skl(2,2)*skl(3,1)*skl(2,2)*skl(3,1)*c1212+
771  & skl(2,2)*skl(3,2)*skl(2,1)*skl(3,1)*c1122+
772  & skl(2,2)*skl(3,2)*skl(2,2)*skl(3,2)*c1111+
773  & skl(2,2)*skl(3,2)*skl(2,3)*skl(3,3)*c1122+
774  & skl(2,2)*skl(3,3)*skl(2,2)*skl(3,3)*c1212+
775  & skl(2,2)*skl(3,3)*skl(2,3)*skl(3,2)*c1212+
776  & skl(2,3)*skl(3,1)*skl(2,1)*skl(3,3)*c1212+
777  & skl(2,3)*skl(3,1)*skl(2,3)*skl(3,1)*c1212+
778  & skl(2,3)*skl(3,2)*skl(2,2)*skl(3,3)*c1212+
779  & skl(2,3)*skl(3,2)*skl(2,3)*skl(3,2)*c1212+
780  & skl(2,3)*skl(3,3)*skl(2,1)*skl(3,1)*c1122+
781  & skl(2,3)*skl(3,3)*skl(2,2)*skl(3,2)*c1122+
782  & skl(2,3)*skl(3,3)*skl(2,3)*skl(3,3)*c1111
783  endif
784 !
785  do i=1,6
786  ep0(i)=xstateini(i,iint,iel)
787  enddo
788  do i=1,18
789  q1(i)=xstateini(6+i,iint,iel)
790  q2(i)=xstateini(24+i,iint,iel)
791  dg0(i)=xstateini(42+i,iint,iel)
792  enddo
793 !
794 ! elastic strains
795 !
796  do i=1,6
797  ee(i)=emec(i)-ep0(i)
798  enddo
799 !
800 ! (visco)plastic constants: octahedral slip system
801 !
802  do i=1,12
803  ck(i)=elconloc(4)
804  cn(i)=elconloc(5)
805  c(i)=elconloc(6)
806  d(i)=elconloc(7)
807  phi(i)=elconloc(8)
808  delta(i)=elconloc(9)
809  r0(i)=elconloc(10)
810  q(i)=elconloc(11)
811  b(i)=elconloc(12)
812  enddo
813 !
814 ! (visco)plastic constants: cubic slip system
815 !
816  do i=13,18
817  ck(i)=elconloc(13)
818  cn(i)=elconloc(14)
819  c(i)=elconloc(15)
820  d(i)=elconloc(16)
821  phi(i)=elconloc(17)
822  delta(i)=elconloc(18)
823  r0(i)=elconloc(19)
824  q(i)=elconloc(20)
825  b(i)=elconloc(21)
826  enddo
827 !
828 ! stress state variables q1 and q2
829 !
830  do i=1,18
831  al10(i)=-q1(i)/(b(i)*q(i))
832  al20(i)=-q2(i)/c(i)
833  enddo
834 !
835 ! global trial stress tensor
836 !
837  if(iorien.gt.0) then
838  stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+
839  & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6))
840  & -beta(1)
841  stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+
842  & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6))
843  & -beta(2)
844  stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+
845  & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6))
846  & -beta(3)
847  stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+
848  & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6))
849  & -beta(4)
850  stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+
851  & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6))
852  & -beta(5)
853  stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+
854  & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6))
855  & -beta(6)
856  else
857  stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1)
858  stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2)
859  stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3)
860  stri(4)=2.d0*c1212*ee(4)-beta(4)
861  stri(5)=2.d0*c1212*ee(5)-beta(5)
862  stri(6)=2.d0*c1212*ee(6)-beta(6)
863  endif
864 !
865 ! stress radius in each slip plane
866 !
867  do i=1,18
868  sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+
869  & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6))+q2(i)
870  enddo
871 !
872 ! evaluation of the yield surface
873 !
874  do i=1,18
875  htri(i)=dabs(sg(i))-r0(i)
876  do j=1,18
877  htri(i)=htri(i)+h(i,j)*q1(j)
878  enddo
879  enddo
880 !
881 ! check whether plasticity occurs
882 !
883  iplas=0
884  do i=1,18
885  if(htri(i).gt.0.d0) then
886  iplas=1
887  exit
888  endif
889  enddo
890 c 8 continue
891 !
892  if((iplas.eq.0).or.(ielas.eq.1)) then
893 !
894 ! elastic stress
895 !
896  do i=1,6
897  stre(i)=stri(i)
898  enddo
899 !
900 ! elastic stiffness
901 !
902  if(icmd.ne.3) then
903  if(iorien.gt.0) then
904  do i=1,21
905  stiff(i)=elas(i)
906  enddo
907  else
908  stiff(1)=c1111
909  stiff(2)=c1122
910  stiff(3)=c1111
911  stiff(4)=c1122
912  stiff(5)=c1122
913  stiff(6)=c1111
914  stiff(7)=0.d0
915  stiff(8)=0.d0
916  stiff(9)=0.d0
917  stiff(10)=c1212
918  stiff(11)=0.d0
919  stiff(12)=0.d0
920  stiff(13)=0.d0
921  stiff(14)=0.d0
922  stiff(15)=c1212
923  stiff(16)=0.d0
924  stiff(17)=0.d0
925  stiff(18)=0.d0
926  stiff(19)=0.d0
927  stiff(20)=0.d0
928  stiff(21)=c1212
929  endif
930  endif
931 !
932  return
933  endif
934 !
935 ! plastic deformation
936 !
937  creep=1
938  nrhs=1
939  lda=18
940  ldb=18
941 !
942 ! determining the active slip planes
943 !
944  do i=1,18
945  if(htri(i).gt.0.d0) then
946  active(i)=1
947  else
948  active(i)=0
949  endif
950  enddo
951 !
952 ! initializing the state variables
953 !
954  do i=1,6
955  ep(i)=ep0(i)
956  enddo
957  do i=1,18
958  al1(i)=al10(i)
959  al2(i)=al20(i)
960  dg(i)=0.d0
961  enddo
962 !
963 ! major loop
964 !
965  icounter=0
966  majorloop: do
967  icounter=icounter+1
968  if(icounter.gt.100) then
969  write(*,*) '*ERROR in umat_single_crystal: no convergence'
970  call exit(201)
971  endif
972 !
973 ! elastic strains
974 !
975  do i=1,6
976  ee(i)=emec(i)-ep(i)
977  enddo
978 !
979 ! stress state variables q1 and q2
980 !
981  do i=1,18
982  q1(i)=-b(i)*q(i)*al1(i)
983  q2(i)=-c(i)*al2(i)
984  enddo
985 !
986 ! global trial stress tensor
987 !
988  if(iorien.gt.0) then
989  stri(1)=elas(1)*ee(1)+elas(2)*ee(2)+elas(4)*ee(3)+
990  & 2.d0*(elas(7)*ee(4)+elas(11)*ee(5)+elas(16)*ee(6))
991  & -beta(1)
992  stri(2)=elas(2)*ee(1)+elas(3)*ee(2)+elas(5)*ee(3)+
993  & 2.d0*(elas(8)*ee(4)+elas(12)*ee(5)+elas(17)*ee(6))
994  & -beta(2)
995  stri(3)=elas(4)*ee(1)+elas(5)*ee(2)+elas(6)*ee(3)+
996  & 2.d0*(elas(9)*ee(4)+elas(13)*ee(5)+elas(18)*ee(6))
997  & -beta(3)
998  stri(4)=elas(7)*ee(1)+elas(8)*ee(2)+elas(9)*ee(3)+
999  & 2.d0*(elas(10)*ee(4)+elas(14)*ee(5)+elas(19)*ee(6))
1000  & -beta(4)
1001  stri(5)=elas(11)*ee(1)+elas(12)*ee(2)+elas(13)*ee(3)+
1002  & 2.d0*(elas(14)*ee(4)+elas(15)*ee(5)+elas(20)*ee(6))
1003  & -beta(5)
1004  stri(6)=elas(16)*ee(1)+elas(17)*ee(2)+elas(18)*ee(3)+
1005  & 2.d0*(elas(19)*ee(4)+elas(20)*ee(5)+elas(21)*ee(6))
1006  & -beta(6)
1007  else
1008  stri(1)=c1111*ee(1)+c1122*(ee(2)+ee(3))-beta(1)
1009  stri(2)=c1111*ee(2)+c1122*(ee(1)+ee(3))-beta(2)
1010  stri(3)=c1111*ee(3)+c1122*(ee(1)+ee(2))-beta(3)
1011  stri(4)=2.d0*c1212*ee(4)-beta(4)
1012  stri(5)=2.d0*c1212*ee(5)-beta(5)
1013  stri(6)=2.d0*c1212*ee(6)-beta(6)
1014  endif
1015 !
1016 ! stress radius in each slip plane
1017 !
1018  do i=1,18
1019  sg(i)=xm(1,i)*stri(1)+xm(2,i)*stri(2)+xm(3,i)*stri(3)+
1020  & 2.d0*(xm(4,i)*stri(4)+xm(5,i)*stri(5)+xm(6,i)*stri(6))
1021  & +q2(i)
1022  enddo
1023 !
1024 ! evaluation of the yield surface
1025 !
1026  do i=1,18
1027  htri(i)=dabs(sg(i))-r0(i)-ck(i)*(dg(i)/dtime)**(1.d0/cn(i))
1028  do j=1,18
1029  htri(i)=htri(i)+h(i,j)*q1(j)
1030  enddo
1031  enddo
1032 !
1033 ! replace sg(i) by sgn(sg(i))
1034 !
1035  do i=1,18
1036  if(sg(i).lt.0.d0) then
1037  sg(i)=-1.d0
1038  else
1039  sg(i)=1.d0
1040  endif
1041  enddo
1042 !
1043 ! determining the effect of the accumulated plasticity
1044 !
1045  do i=1,18
1046  cphi(i)=phi(i)+(1.d0-phi(i))*exp(-delta(i)*(dg0(i)+dg(i)))
1047  enddo
1048 !
1049 ! minor loop
1050 !
1051  minorloop: do
1052 !
1053 ! determining the residual matrix
1054 !
1055  do i=1,6
1056  r(i)=ep0(i)-ep(i)
1057  enddo
1058  do i=1,18
1059  r(5+2*i)=al10(i)-al1(i)
1060  r(6+2*i)=al20(i)-al2(i)
1061  enddo
1062  do i=1,18
1063  if(active(i).eq.1) then
1064  do j=1,6
1065  r(j)=r(j)+xm(j,i)*sg(i)*dg(i)
1066  enddo
1067  r(5+2*i)=r(5+2*i)+(1.d0-b(i)*al1(i))*dg(i)
1068  r(6+2*i)=r(6+2*i)+(cphi(i)*sg(i)-d(i)*al2(i))*dg(i)
1069  endif
1070  enddo
1071 !
1072 ! check convergence
1073 !
1074  convergence=1
1075  do i=1,18
1076  if(active(i).eq.0) cycle
1077  if(htri(i).gt.1.d-5) then
1078  convergence=0
1079  exit
1080  endif
1081  enddo
1082 c 9 continue
1083  if(convergence.eq.1) then
1084  dd=0.d0
1085  do i=1,6
1086  dd=dd+r(i)*r(i)
1087  enddo
1088  do i=1,18
1089  if(active(i).eq.0) cycle
1090  dd=dd+r(5+2*i)*r(5+2*i)+r(6+2*i)*r(6+2*i)
1091  enddo
1092  dd=sqrt(dd)
1093  if(dd.gt.1.d-10) then
1094  convergence=0
1095  else
1096  exit majorloop
1097 c go to 12
1098  endif
1099  endif
1100 !
1101 ! compute xmc=c:xm
1102 !
1103  do i=1,18
1104  if(iorien.gt.0) then
1105  xmc(1,i)=elas(1)*xm(1,i)+elas(2)*xm(2,i)+
1106  & elas(4)*xm(3,i)+2.d0*(elas(7)*xm(4,i)+
1107  & elas(11)*xm(5,i)+elas(16)*xm(6,i))
1108  xmc(2,i)=elas(2)*xm(1,i)+elas(3)*xm(2,i)+
1109  & elas(5)*xm(3,i)+2.d0*(elas(8)*xm(4,i)+
1110  & elas(12)*xm(5,i)+elas(17)*xm(6,i))
1111  xmc(3,i)=elas(4)*xm(1,i)+elas(5)*xm(2,i)+
1112  & elas(6)*xm(3,i)+2.d0*(elas(9)*xm(4,i)+
1113  & elas(13)*xm(5,i)+elas(18)*xm(6,i))
1114  xmc(4,i)=elas(7)*xm(1,i)+elas(8)*xm(2,i)+
1115  & elas(9)*xm(3,i)+2.d0*(elas(10)*xm(4,i)+
1116  & elas(14)*xm(5,i)+elas(19)*xm(6,i))
1117  xmc(5,i)=elas(11)*xm(1,i)+elas(12)*xm(2,i)+
1118  & elas(13)*xm(3,i)+2.d0*(elas(14)*xm(4,i)+
1119  & elas(15)*xm(5,i)+elas(20)*xm(6,i))
1120  xmc(6,i)=elas(16)*xm(1,i)+elas(17)*xm(2,i)+
1121  & elas(18)*xm(3,i)+2.d0*(elas(19)*xm(4,i)+
1122  & elas(20)*xm(5,i)+elas(21)*xm(6,i))
1123  else
1124  xmc(1,i)=c1111*xm(1,i)+c1122*(xm(2,i)+xm(3,i))
1125  xmc(2,i)=c1111*xm(2,i)+c1122*(xm(1,i)+xm(3,i))
1126  xmc(3,i)=c1111*xm(3,i)+c1122*(xm(1,i)+xm(2,i))
1127  xmc(4,i)=2.d0*c1212*xm(4,i)
1128  xmc(5,i)=2.d0*c1212*xm(5,i)
1129  xmc(6,i)=2.d0*c1212*xm(6,i)
1130  endif
1131  enddo
1132 !
1133 ! indexing the active slip planes
1134 !
1135  do i=1,18
1136  if(active(i).eq.1) then
1137  index(i)=1.d0
1138  else
1139  index(i)=0.d0
1140  endif
1141  enddo
1142  neq=0
1143  do i=1,18
1144  if(index(i).eq.1) then
1145  neq=neq+1
1146  index(i)=neq
1147  endif
1148  enddo
1149 !
1150 ! filling the LHS
1151 !
1152  do i=1,18
1153  if(active(i).eq.0) cycle
1154  aux(i)=(q(i)+q1(i))/(1.d0/b(i)+dg(i))
1155  enddo
1156 !
1157  do i=1,18
1158  if(active(i).eq.0) cycle
1159  do j=1,18
1160  if(active(j).eq.0) cycle
1161  if(i.ne.j) then
1162  gl(index(i),index(j))=(xm(1,i)*xmc(1,j)+
1163  & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0*
1164  & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+
1165  & xm(6,i)*xmc(6,j)))
1166  & *sg(i)*sg(j)+h(i,j)*aux(j)
1167  else
1168  gl(index(i),index(j))=(xm(1,i)*xmc(1,j)+
1169  & xm(2,i)*xmc(2,j)+xm(3,i)*xmc(3,j)+2.d0*
1170  & (xm(4,i)*xmc(4,j)+xm(5,i)*xmc(5,j)+
1171  & xm(6,i)*xmc(6,j)))
1172  & +h(i,j)*aux(j)+(cphi(j)*c(j)+d(j)*q2(j)*sg(j))
1173  & /(1.d0+dg(j)*d(j))
1174  endif
1175  enddo
1176 c 3 continue
1177  if(creep.eq.1)then
1178  if(dg(i).gt.0.d0) then
1179  gl(index(i),index(i))=gl(index(i),index(i))+
1180  & (dg(i)/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/
1181  & (cn(i)*dtime)
1182  else
1183 !
1184 ! for gamma ein default of 1.d-10 is taken to
1185 ! obtain a finite gradient
1186 !
1187  gl(index(i),index(i))=gl(index(i),index(i))+
1188  & (1.d-10/dtime)**(1.d0/cn(i)-1.d0)*ck(i)/
1189  & (cn(i)*dtime)
1190  endif
1191  endif
1192  enddo
1193 c 2 continue
1194 !
1195 ! filling the RHS
1196 !
1197  do i=1,18
1198  if(active(i).eq.0) cycle
1199  do j=1,6
1200  t(j)=xmc(j,i)*sg(i)
1201  enddo
1202  do j=1,18
1203  t(5+2*j)=h(i,j)*q(j)/(1.d0/b(j)+dg(j))
1204  t(6+2*j)=0.d0
1205  enddo
1206  t(6+2*i)=c(i)*sg(i)/(1.d0+dg(i)*d(i))
1207 !
1208  if(creep.eq.1) then
1209  gr(index(i),1)=htri(i)
1210  else
1211  gr(index(i),1)=htri(i)
1212  & +ck(i)*(dg(i)/dtime)**(1.d0/cn(i))
1213  endif
1214 !
1215  do j=1,42
1216  gr(index(i),1)=gr(index(i),1)-t(j)*r(j)
1217  enddo
1218  gr(index(i),1)=gr(index(i),1)
1219  & -t(4)*r(4)-t(5)*r(5)-t(6)*r(6)
1220  enddo
1221 c 4 continue
1222 !
1223 ! solve gl*ddg=gr
1224 !
1225  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
1226  if(info.ne.0) then
1227  write(*,*) '*ERROR in sc.f: linear equation solver'
1228  write(*,*) ' exited with error: info = ',info
1229  call exit(201)
1230  endif
1231 !
1232  do i=1,18
1233  if(active(i).eq.1) then
1234  ddg(i)=gr(index(i),1)
1235  else
1236  ddg(i)=0.d0
1237  endif
1238  enddo
1239 !
1240 ! check whether active slip planes have changed
1241 !
1242  ichange=0
1243  do i=1,18
1244  if(active(i).eq.0) cycle
1245  if(dg(i)+ddg(i).lt.0.d0) then
1246  active(i)=0
1247  dg(i)=0.d0
1248  al1(i)=al10(i)
1249  al2(i)=al20(i)
1250  ichange=1
1251  endif
1252  enddo
1253 c 5 continue
1254  if(ichange.eq.0) then
1255  exit minorloop
1256 c go to 13
1257  endif
1258 !
1259 ! end of minor loop
1260 !
1261  enddo minorloop
1262 c 13 continue
1263 !
1264 ! updating the residual matrix
1265 !
1266  do i=1,18
1267  if(active(i).eq.1) then
1268  do j=1,6
1269  r(j)=r(j)+xm(j,i)*sg(i)*ddg(i)
1270  enddo
1271  r(5+2*i)=r(5+2*i)+(1.d0-b(i)*al1(i))*ddg(i)
1272  r(6+2*i)=r(6+2*i)+(cphi(i)*sg(i)-d(i)*al2(i))*ddg(i)
1273  endif
1274  enddo
1275 !
1276 ! update the state variables
1277 !
1278  do i=1,6
1279  ep(i)=ep(i)+r(i)
1280  enddo
1281  do i=1,18
1282  if(active(i).eq.1) then
1283  al1(i)=al1(i)+r(5+2*i)/(1.d0+b(i)*dg(i))
1284  al2(i)=al2(i)+r(6+2*i)/(1.d0+d(i)*dg(i))
1285  endif
1286  enddo
1287  do i=1,18
1288  if(active(i).eq.1) then
1289  dg(i)=dg(i)+ddg(i)
1290  endif
1291  enddo
1292 !
1293 ! end of major loop
1294 !
1295  enddo majorloop
1296 c 12 continue
1297 !
1298 ! inversion of G
1299 !
1300  do i=1,neq
1301  do j=1,neq
1302  gr(i,j)=0.d0
1303  enddo
1304  gr(i,i)=1.d0
1305  enddo
1306  nrhs=neq
1307  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
1308  if(info.ne.0) then
1309  write(*,*) '*ERROR in sc.f: linear equation solver'
1310  write(*,*) ' exited with error: info = ',info
1311  call exit(201)
1312  endif
1313 !
1314 ! storing the stress
1315 !
1316  do i=1,6
1317  stre(i)=stri(i)
1318  enddo
1319 !
1320 ! calculating the tangent stiffness matrix
1321 !
1322  if(icmd.ne.3) then
1323  if(iorien.gt.0) then
1324  ddsdde(1,1)=elas(1)
1325  ddsdde(1,2)=elas(2)
1326  ddsdde(1,3)=elas(4)
1327  ddsdde(1,4)=elas(7)
1328  ddsdde(1,5)=elas(11)
1329  ddsdde(1,6)=elas(16)
1330  ddsdde(2,1)=elas(2)
1331  ddsdde(2,2)=elas(3)
1332  ddsdde(2,3)=elas(5)
1333  ddsdde(2,4)=elas(8)
1334  ddsdde(2,5)=elas(12)
1335  ddsdde(2,6)=elas(17)
1336  ddsdde(3,1)=elas(4)
1337  ddsdde(3,2)=elas(5)
1338  ddsdde(3,3)=elas(6)
1339  ddsdde(3,4)=elas(9)
1340  ddsdde(3,5)=elas(13)
1341  ddsdde(3,6)=elas(18)
1342  ddsdde(4,1)=elas(7)
1343  ddsdde(4,2)=elas(8)
1344  ddsdde(4,3)=elas(9)
1345  ddsdde(4,4)=elas(10)
1346  ddsdde(4,5)=elas(14)
1347  ddsdde(4,6)=elas(19)
1348  ddsdde(5,1)=elas(11)
1349  ddsdde(5,2)=elas(12)
1350  ddsdde(5,3)=elas(13)
1351  ddsdde(5,4)=elas(14)
1352  ddsdde(5,5)=elas(15)
1353  ddsdde(5,6)=elas(20)
1354  ddsdde(6,1)=elas(16)
1355  ddsdde(6,2)=elas(17)
1356  ddsdde(6,3)=elas(18)
1357  ddsdde(6,4)=elas(19)
1358  ddsdde(6,5)=elas(20)
1359  ddsdde(6,6)=elas(21)
1360  else
1361  do i=1,6
1362  do j=1,6
1363  ddsdde(i,j)=0.d0
1364  enddo
1365  enddo
1366  do i=1,3
1367  ddsdde(i,i)=c1111
1368  enddo
1369  do i=1,3
1370  do j=i+1,3
1371  ddsdde(i,j)=c1122
1372  enddo
1373  do j=1,i-1
1374  ddsdde(i,j)=c1122
1375  enddo
1376  ddsdde(i+3,i+3)=c1212
1377  enddo
1378  endif
1379  do i=1,18
1380  if(active(i).eq.0) cycle
1381  do j=1,18
1382  if(active(j).eq.0) cycle
1383  do k=1,6
1384  do l=1,6
1385  ddsdde(k,l)=ddsdde(k,l)-
1386  & gr(index(i),index(j))*xmc(k,i)*sg(i)*xmc(l,j)*sg(j)
1387  enddo
1388  enddo
1389  enddo
1390  enddo
1391 c 7 continue
1392 c 6 continue
1393 !
1394 ! symmatrizing the stiffness matrix
1395 !
1396  stiff(1)=ddsdde(1,1)
1397  stiff(2)=(ddsdde(1,2)+ddsdde(2,1))/2.d0
1398  stiff(3)=ddsdde(2,2)
1399  stiff(4)=(ddsdde(1,3)+ddsdde(3,1))/2.d0
1400  stiff(5)=(ddsdde(2,3)+ddsdde(3,2))/2.d0
1401  stiff(6)=ddsdde(3,3)
1402  stiff(7)=(ddsdde(1,4)+ddsdde(4,1))/2.d0
1403  stiff(8)=(ddsdde(2,4)+ddsdde(4,2))/2.d0
1404  stiff(9)=(ddsdde(3,4)+ddsdde(4,3))/2.d0
1405  stiff(10)=ddsdde(4,4)
1406  stiff(11)=(ddsdde(1,5)+ddsdde(5,1))/2.d0
1407  stiff(12)=(ddsdde(2,5)+ddsdde(5,2))/2.d0
1408  stiff(13)=(ddsdde(3,5)+ddsdde(5,3))/2.d0
1409  stiff(14)=(ddsdde(4,5)+ddsdde(5,4))/2.d0
1410  stiff(15)=ddsdde(5,5)
1411  stiff(16)=(ddsdde(1,6)+ddsdde(6,1))/2.d0
1412  stiff(17)=(ddsdde(2,6)+ddsdde(6,2))/2.d0
1413  stiff(18)=(ddsdde(3,6)+ddsdde(6,3))/2.d0
1414  stiff(19)=(ddsdde(4,6)+ddsdde(6,4))/2.d0
1415  stiff(20)=(ddsdde(5,6)+ddsdde(6,5))/2.d0
1416  stiff(21)=ddsdde(6,6)
1417 !
1418  endif
1419 !
1420 ! updating the state variables
1421 !
1422  do i=1,6
1423  xstate(i,iint,iel)=ep(i)
1424  enddo
1425  do i=1,18
1426  xstate(6+i,iint,iel)=q1(i)
1427  xstate(24+i,iint,iel)=q2(i)
1428  xstate(42+i,iint,iel)=dg0(i)+dg(i)
1429  enddo
1430 !
1431  return
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
active
subroutine active(n, l, u, nbd, x, iwhere, iprint, prjctd, cnstnd, boxed)
Definition: lbfgsb.f:1006
creep
subroutine creep(decra, deswa, statev, serd, ec, esw, p, qtild, temp, dtemp, predef, dpred, time, dtime, cmname, leximp, lend, coords, nstatv, noel, npt, layer, kspt, kstep, kinc)
Definition: creep.f:22
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)