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

Go to the source code of this file.

Functions/Subroutines

subroutine umat_aniso_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, nmethod, pnewdt, depvisc)
 

Function/Subroutine Documentation

◆ umat_aniso_creep()

subroutine umat_aniso_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,
integer  nmethod,
real*8  pnewdt,
real*8  depvisc 
)
23 !
24 ! calculates stiffness and stresses for a user defined material
25 ! law
26 !
27 ! icmd=3: calculates stress at mechanical strain
28 ! else: calculates stress at mechanical strain and the stiffness
29 ! matrix
30 !
31 ! INPUT:
32 !
33 ! amat material name
34 ! iel element number
35 ! iint integration point number
36 !
37 ! kode material type (-100-#of constants entered
38 ! under *USER MATERIAL): can be used for materials
39 ! with varying number of constants
40 !
41 ! elconloc(21) user defined constants defined by the keyword
42 ! card *USER MATERIAL (max. 21, actual # =
43 ! -kode-100), interpolated for the
44 ! actual temperature t1l
45 !
46 ! emec(6) Lagrange mechanical strain tensor (component order:
47 ! 11,22,33,12,13,23) at the end of the increment
48 ! (thermal strains are subtracted)
49 ! emec0(6) Lagrange mechanical strain tensor at the start of the
50 ! increment (thermal strains are subtracted)
51 ! beta(6) residual stress tensor (the stress entered under
52 ! the keyword *INITIAL CONDITIONS,TYPE=STRESS)
53 !
54 ! xokl(3,3) deformation gradient at the start of the increment
55 ! voj Jacobian at the start of the increment
56 ! xkl(3,3) deformation gradient at the end of the increment
57 ! vj Jacobian at the end of the increment
58 !
59 ! ithermal 0: no thermal effects are taken into account: for
60 ! creep this does not make sense.
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 exitcriterion
119 !
120  character*80 amat
121 !
122  integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien,
123  & i,j,ipiv(6),info,neq,lda,ldb,j1,j2,j3,j4,j5,j6,j7,j8,
124  & nrhs,iplas,kel(4,21),iloop,leximp,lend,layer,kspt,kstep,
125  & kinc,ii,nmethod
126 !
127  real*8 ep0(6),epqini,ep(6),b,pn(6),dg,ddg,c(21),x(21),cm1(21),
128  & stri(6),htri,sg(6),r(13),ee(6),dd,gl(6,6),gr(6,6),c0,c1,c2,
129  & skl(3,3),gcreep,gm1,ya(3,3,3,3),dsg,detc,strinv,depvisc,
130  & depq,svm,dsvm,dg1,dg2,fu,fu1,fu2,expon,ec(2),pnewdt,
131  & timeabq(2),r1(13),ep1(6),gl1(6,6),sg1(6),ckl(3,3),
132  & elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
133  & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
134  & time,ttime,decra(5),deswa(5),serd,esw(2),p,predef(1),dpred(1),
135  & dtemp,xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
136 !
137  kel=reshape((/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
138  & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,
139  & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,
140  & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/))
141 !
142  leximp=1
143  lend=3
144 !
145  if(ithermal.eq.0) then
146  write(*,*)'*ERROR in umat_aniso_creep: no temperature defined;'
147  write(*,*) ' a creep calculation without temperature'
148  write(*,*) ' does not make sense'
149  write(*,*)
150  call exit(201)
151  endif
152 !
153  iloop=0
154  exitcriterion=0
155 !
156  c0=dsqrt(2.d0/3.d0)
157  c1=2.d0/3.d0
158  c2=-1.d0/3.d0
159 !
160 ! elastic constants
161 !
162  if(iorien.gt.0) then
163 !
164  call transformatrix(orab(1,iorien),pgauss,skl)
165 !
166  call orthotropic(elconloc,ya)
167 !
168  do j=1,21
169  j1=kel(1,j)
170  j2=kel(2,j)
171  j3=kel(3,j)
172  j4=kel(4,j)
173  c(j)=0.d0
174  do j5=1,3
175  do j6=1,3
176  do j7=1,3
177  do j8=1,3
178  c(j)=c(j)+ya(j5,j6,j7,j8)*
179  & skl(j1,j5)*skl(j2,j6)*skl(j3,j7)*skl(j4,j8)
180  enddo
181  enddo
182  enddo
183  enddo
184  enddo
185 !
186  else
187  do i=1,9
188  c(i)=elconloc(i)
189  enddo
190  endif
191 !
192 ! state variables
193 !
194 ! equivalent plastic strain
195 !
196  epqini=xstateini(1,iint,iel)
197 !
198 ! plastic strain
199 !
200  do i=1,6
201  ep0(i)=xstateini(1+i,iint,iel)
202  enddo
203 ! elastic strains
204 !
205  do i=1,6
206  ee(i)=emec(i)-ep0(i)
207  enddo
208 !
209 ! global trial stress tensor
210 !
211  if(iorien.gt.0) then
212  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+
213  & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6))
214  & -beta(1)
215  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+
216  & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6))
217  & -beta(2)
218  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+
219  & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6))
220  & -beta(3)
221  stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+
222  & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6))
223  & -beta(4)
224  stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+
225  & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6))
226  & -beta(5)
227  stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+
228  & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6))
229  & -beta(6)
230  else
231  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1)
232  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1)
233  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1)
234  stri(4)=2.d0*c(7)*ee(4)-beta(4)
235  stri(5)=2.d0*c(8)*ee(5)-beta(5)
236  stri(6)=2.d0*c(9)*ee(6)-beta(6)
237  endif
238 !
239 ! stress radius (only deviatoric part of stress enters)
240 !
241  strinv=(stri(1)+stri(2)+stri(3))/3.d0
242  do i=1,3
243  sg(i)=stri(i)-strinv
244  enddo
245  do i=4,6
246  sg(i)=stri(i)
247  enddo
248  dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+
249  & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6)))
250 !
251 ! evaluation of the yield surface
252 !
253  ec(1)=epqini
254 !
255  htri=dsg
256 !
257 ! check whether plasticity occurs
258 !
259  if(htri.gt.1.d-10) then
260  iplas=1
261  else
262  iplas=0
263  endif
264 !
265  if((iplas.eq.0).or.(ielas.eq.1).or.(dtime.lt.1.d-30).or.
266  & ((nmethod.eq.1).and.(ithermal.ne.3))) then
267 !
268 ! elastic stress
269 !
270  do i=1,6
271  stre(i)=stri(i)
272  enddo
273 !
274 ! updating the state variables
275 !
276  xstate(1,iint,iel)=epqini
277  do i=1,6
278  xstate(1+i,iint,iel)=ep0(i)
279  enddo
280 !
281 ! elastic stiffness
282 !
283  if(icmd.ne.3) then
284  if(iorien.gt.0) then
285  do i=1,21
286  stiff(i)=c(i)
287  enddo
288  else
289  stiff(1)=c(1)
290  stiff(2)=c(2)
291  stiff(3)=c(3)
292  stiff(4)=c(4)
293  stiff(5)=c(5)
294  stiff(6)=c(6)
295  stiff(7)=0.d0
296  stiff(8)=0.d0
297  stiff(9)=0.d0
298  stiff(10)=c(7)
299  stiff(11)=0.d0
300  stiff(12)=0.d0
301  stiff(13)=0.d0
302  stiff(14)=0.d0
303  stiff(15)=c(8)
304  stiff(16)=0.d0
305  stiff(17)=0.d0
306  stiff(18)=0.d0
307  stiff(19)=0.d0
308  stiff(20)=0.d0
309  stiff(21)=c(9)
310  endif
311  endif
312 !
313  return
314  endif
315 !
316 ! plastic deformation
317 !
318  neq=6
319  nrhs=1
320  lda=6
321  ldb=6
322 !
323 ! initializing the state variables
324 !
325  do i=1,6
326  ep(i)=ep0(i)
327  enddo
328  dg=0.d0
329 !
330 ! determining the inverse of c
331 !
332  if(iorien.gt.0) then
333 !
334 ! solve gl:C=gr
335 !
336  gl(1,1)=c(1)
337  gl(1,2)=c(2)
338  gl(2,2)=c(3)
339  gl(1,3)=c(4)
340  gl(2,3)=c(5)
341  gl(3,3)=c(6)
342  gl(1,4)=c(7)
343  gl(2,4)=c(8)
344  gl(3,4)=c(9)
345  gl(4,4)=c(10)
346  gl(1,5)=c(11)
347  gl(2,5)=c(12)
348  gl(3,5)=c(13)
349  gl(4,5)=c(14)
350  gl(5,5)=c(15)
351  gl(1,6)=c(16)
352  gl(2,6)=c(17)
353  gl(3,6)=c(18)
354  gl(4,6)=c(19)
355  gl(5,6)=c(20)
356  gl(6,6)=c(21)
357  do i=1,6
358  do j=1,i-1
359  gl(i,j)=gl(j,i)
360  enddo
361  enddo
362  do i=1,6
363  do j=1,6
364  gr(i,j)=0.d0
365  enddo
366  gr(i,i)=1.d0
367  enddo
368  nrhs=6
369  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
370  if(info.ne.0) then
371  write(*,*) '*ERROR in sc.f: linear equation solver'
372  write(*,*) ' exited with error: info = ',info
373  call exit(201)
374  endif
375  nrhs=1
376  cm1(1)=gr(1,1)
377  cm1(2)=gr(1,2)
378  cm1(3)=gr(2,2)
379  cm1(4)=gr(1,3)
380  cm1(5)=gr(2,3)
381  cm1(6)=gr(3,3)
382  cm1(7)=gr(1,4)/2.d0
383  cm1(8)=gr(2,4)/2.d0
384  cm1(9)=gr(3,4)/2.d0
385  cm1(10)=gr(4,4)/4.d0
386  cm1(11)=gr(1,5)/2.d0
387  cm1(12)=gr(2,5)/2.d0
388  cm1(13)=gr(3,5)/2.d0
389  cm1(14)=gr(4,5)/4.d0
390  cm1(15)=gr(5,5)/4.d0
391  cm1(16)=gr(1,6)/2.d0
392  cm1(17)=gr(2,6)/2.d0
393  cm1(18)=gr(3,6)/2.d0
394  cm1(19)=gr(4,6)/4.d0
395  cm1(20)=gr(5,6)/4.d0
396  cm1(21)=gr(6,6)/4.d0
397  else
398  detc=c(1)*(c(3)*c(6)-c(5)*c(5))-
399  & c(2)*(c(2)*c(6)-c(4)*c(5))+
400  & c(4)*(c(2)*c(5)-c(4)*c(3))
401  cm1(1)=(c(3)*c(6)-c(5)*c(5))/detc
402  cm1(2)=(c(5)*c(4)-c(2)*c(6))/detc
403  cm1(3)=(c(1)*c(6)-c(4)*c(4))/detc
404  cm1(4)=(c(2)*c(5)-c(3)*c(4))/detc
405  cm1(5)=(c(2)*c(4)-c(1)*c(5))/detc
406  cm1(6)=(c(1)*c(3)-c(2)*c(2))/detc
407  cm1(7)=1.d0/(4.d0*c(7))
408  cm1(8)=1.d0/(4.d0*c(8))
409  cm1(9)=1.d0/(4.d0*c(9))
410  endif
411 !
412 ! first attempt: root search with Newton-Raphson
413 !
414  loop: do
415 !
416  iloop=iloop+1
417 !
418 ! elastic strains
419 !
420  do i=1,6
421  ee(i)=emec(i)-ep(i)
422  enddo
423 !
424 ! global trial stress tensor
425 !
426  if(iorien.gt.0) then
427  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+
428  & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6))
429  & -beta(1)
430  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+
431  & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6))
432  & -beta(2)
433  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+
434  & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6))
435  & -beta(3)
436  stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+
437  & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6))
438  & -beta(4)
439  stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+
440  & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6))
441  & -beta(5)
442  stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+
443  & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6))
444  & -beta(6)
445  else
446  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1)
447  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1)
448  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1)
449  stri(4)=2.d0*c(7)*ee(4)-beta(4)
450  stri(5)=2.d0*c(8)*ee(5)-beta(5)
451  stri(6)=2.d0*c(9)*ee(6)-beta(6)
452  endif
453 !
454 ! stress radius (only deviatoric part of stress enters)
455 !
456  strinv=(stri(1)+stri(2)+stri(3))/3.d0
457  do i=1,3
458  sg(i)=stri(i)-strinv
459  enddo
460  do i=4,6
461  sg(i)=stri(i)
462  enddo
463  dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+
464  & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6)))
465 !
466 ! evaluation of the yield surface
467 !
468  ec(1)=epqini
469  decra(1)=c0*dg
470  timeabq(1)=time
471  timeabq(2)=ttime+time
472  call creep(decra,deswa,xstateini(1,iint,iel),serd,ec,
473  & esw,p,svm,t1l,dtemp,predef,dpred,timeabq,dtime,
474  & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt,
475  & kstep,kinc)
476 !
477 ! if the creep routine returns an increased value of decra(1)
478 ! it means that there is a lower cut-off for decra(1);
479 ! if the routine stays in a range lower than this cut-off,
480 ! it will never leave it and the exit conditions are
481 ! assumed to be satisfied.
482 !
483  if(decra(1).gt.c0*dg) then
484  dg=decra(1)/c0
485  if(iloop.gt.1) exitcriterion=1
486  endif
487 !
488  htri=dsg-c0*svm
489 !
490  do i=1,6
491  sg(i)=sg(i)/dsg
492  enddo
493 !
494 ! determining the residual matrix
495 !
496  do i=1,6
497  r(i)=ep0(i)-ep(i)+dg*sg(i)
498  enddo
499 !
500 ! check convergence
501 !
502  if(exitcriterion.eq.1) exit
503  if((dabs(htri).le.1.d-3).and.
504  & ((iloop.gt.1).and.((dabs(ddg).lt.1.d-10).or.
505  & (dabs(ddg).lt.1.d-3*dabs(dg))))) then
506  dd=0.d0
507  do i=1,6
508  dd=dd+r(i)*r(i)
509  enddo
510  dd=sqrt(dd)
511  if(dd.le.1.d-10) then
512  exit
513  endif
514  endif
515 !
516 ! determining b.x
517 !
518  b=dg/dsg
519 !
520  x(1)=b*(c1-sg(1)*sg(1))
521  x(2)=b*(c2-sg(1)*sg(2))
522  x(3)=b*(c1-sg(2)*sg(2))
523  x(4)=b*(c2-sg(1)*sg(3))
524  x(5)=b*(c2-sg(2)*sg(3))
525  x(6)=b*(c1-sg(3)*sg(3))
526  x(7)=-b*sg(1)*sg(4)
527  x(8)=-b*sg(2)*sg(4)
528  x(9)=-b*sg(3)*sg(4)
529  x(10)=b*(.5d0-sg(4)*sg(4))
530  x(11)=-b*sg(1)*sg(5)
531  x(12)=-b*sg(2)*sg(5)
532  x(13)=-b*sg(3)*sg(5)
533  x(14)=-b*sg(4)*sg(5)
534  x(15)=b*(.5d0-sg(5)*sg(5))
535  x(16)=-b*sg(1)*sg(6)
536  x(17)=-b*sg(2)*sg(6)
537  x(18)=-b*sg(3)*sg(6)
538  x(19)=-b*sg(4)*sg(6)
539  x(20)=-b*sg(5)*sg(6)
540  x(21)=b*(.5d0-sg(6)*sg(6))
541 !
542 ! filling the LHS
543 !
544  if(iorien.gt.0) then
545  gl(1,1)=cm1(1)+x(1)
546  gl(1,2)=cm1(2)+x(2)
547  gl(2,2)=cm1(3)+x(3)
548  gl(1,3)=cm1(4)+x(4)
549  gl(2,3)=cm1(5)+x(5)
550  gl(3,3)=cm1(6)+x(6)
551  gl(1,4)=cm1(7)+x(7)
552  gl(2,4)=cm1(8)+x(8)
553  gl(3,4)=cm1(9)+x(9)
554  gl(4,4)=cm1(10)+x(10)
555  gl(1,5)=cm1(11)+x(11)
556  gl(2,5)=cm1(12)+x(12)
557  gl(3,5)=cm1(13)+x(13)
558  gl(4,5)=cm1(14)+x(14)
559  gl(5,5)=cm1(15)+x(15)
560  gl(1,6)=cm1(16)+x(16)
561  gl(2,6)=cm1(17)+x(17)
562  gl(3,6)=cm1(18)+x(18)
563  gl(4,6)=cm1(19)+x(19)
564  gl(5,6)=cm1(20)+x(20)
565  gl(6,6)=cm1(21)+x(21)
566  do i=1,6
567  do j=1,i-1
568  gl(i,j)=gl(j,i)
569  enddo
570  enddo
571  else
572  gl(1,1)=cm1(1)+x(1)
573  gl(1,2)=cm1(2)+x(2)
574  gl(2,2)=cm1(3)+x(3)
575  gl(1,3)=cm1(4)+x(4)
576  gl(2,3)=cm1(5)+x(5)
577  gl(3,3)=cm1(6)+x(6)
578  gl(1,4)=x(7)
579  gl(2,4)=x(8)
580  gl(3,4)=x(9)
581  gl(4,4)=cm1(7)+x(10)
582  gl(1,5)=x(11)
583  gl(2,5)=x(12)
584  gl(3,5)=x(13)
585  gl(4,5)=x(14)
586  gl(5,5)=cm1(8)+x(15)
587  gl(1,6)=x(16)
588  gl(2,6)=x(17)
589  gl(3,6)=x(18)
590  gl(4,6)=x(19)
591  gl(5,6)=x(20)
592  gl(6,6)=cm1(9)+x(21)
593  do i=1,6
594  do j=1,i-1
595  gl(i,j)=gl(j,i)
596  enddo
597  enddo
598  endif
599 !
600 ! filling the RHS
601 !
602  do i=1,6
603  gr(i,1)=sg(i)
604  enddo
605 !
606 ! solve gl:(P:n)=gr
607 !
608  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
609  if(info.ne.0) then
610  write(*,*) '*ERROR in sc.f: linear equation solver'
611  write(*,*) ' exited with error: info = ',info
612  call exit(201)
613  endif
614 !
615  do i=1,6
616  pn(i)=gr(i,1)
617  enddo
618 !
619 ! calculating the creep contribution
620 !
621  gcreep=c1/decra(5)
622 !
623 ! calculating the correction to the consistency parameter
624 !
625  gm1=pn(1)*sg(1)+pn(2)*sg(2)+pn(3)*sg(3)+
626  & (pn(4)*sg(4)+pn(5)*sg(5)+pn(6)*sg(6))
627  gm1=1.d0/(gm1+gcreep)
628  ddg=gm1*(htri-(pn(1)*r(1)+pn(2)*r(2)+pn(3)*r(3)+
629  & (pn(4)*r(4)+pn(5)*r(5)+pn(6)*r(6))))
630 !
631 ! updating the residual matrix
632 !
633  do i=1,6
634  r(i)=r(i)+ddg*sg(i)
635  enddo
636 !
637 ! update the plastic strain
638 !
639  gr(1,1)=r(1)
640  gr(2,1)=r(2)
641  gr(3,1)=r(3)
642  gr(4,1)=r(4)
643  gr(5,1)=r(5)
644  gr(6,1)=r(6)
645 !
646  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
647  if(info.ne.0) then
648  write(*,*) '*ERROR in sc.f: linear equation solver'
649  write(*,*) ' exited with error: info = ',info
650  call exit(201)
651  endif
652 !
653  if(iorien.gt.0) then
654  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1)+
655  & (cm1(7)*gr(4,1)+cm1(11)*gr(5,1)+cm1(16)*gr(6,1))
656  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1)+
657  & (cm1(8)*gr(4,1)+cm1(12)*gr(5,1)+cm1(17)*gr(6,1))
658  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1)+
659  & (cm1(9)*gr(4,1)+cm1(13)*gr(5,1)+cm1(18)*gr(6,1))
660  ep(4)=ep(4)+cm1(7)*gr(1,1)+cm1(8)*gr(2,1)+cm1(9)*gr(3,1)+
661  & (cm1(10)*gr(4,1)+cm1(14)*gr(5,1)+cm1(19)*gr(6,1))
662  ep(5)=ep(5)+cm1(11)*gr(1,1)+cm1(12)*gr(2,1)+cm1(13)*gr(3,1)+
663  & (cm1(14)*gr(4,1)+cm1(15)*gr(5,1)+cm1(20)*gr(6,1))
664  ep(6)=ep(6)+cm1(16)*gr(1,1)+cm1(17)*gr(2,1)+cm1(18)*gr(3,1)+
665  & (cm1(19)*gr(4,1)+cm1(20)*gr(5,1)+cm1(21)*gr(6,1))
666  else
667  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1)
668  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1)
669  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1)
670  ep(4)=ep(4)+cm1(7)*gr(4,1)
671  ep(5)=ep(5)+cm1(8)*gr(5,1)
672  ep(6)=ep(6)+cm1(9)*gr(6,1)
673  endif
674 !
675 ! update the consistency parameter
676 !
677  dg=dg+ddg
678 !
679 ! end of major loop
680 !
681  if((iloop.gt.15).or.(dg.le.0.d0)) then
682  iloop=1
683  dg=0.d0
684  do i=1,6
685  ep(i)=ep0(i)
686  enddo
687 !
688 ! second attempt: root search through interval division
689 !
690  do
691  if(iloop.gt.100) then
692 c NOTE: write statements cause problems for
693 c parallellized execution
694 c write(*,*)
695 c & '*WARNING in umat_aniso_creep: material loop'
696 c write(*,*) ' did not converge in integration'
697 c write(*,*) ' point',iint,'in element',iel,';'
698 c write(*,*) ' the increment size is reduced'
699 c write(*,*)
700  pnewdt=0.25d0
701  return
702  endif
703 !
704 ! elastic strains
705 !
706  do i=1,6
707  ee(i)=emec(i)-ep(i)
708  enddo
709 !
710 ! global trial stress tensor
711 !
712  if(iorien.gt.0) then
713  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+
714  & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6))
715  & -beta(1)
716  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+
717  & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6))
718  & -beta(2)
719  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+
720  & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6))
721  & -beta(3)
722  stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+
723  & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6))
724  & -beta(4)
725  stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+
726  & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6))
727  & -beta(5)
728  stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+
729  & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6))
730  & -beta(6)
731  else
732  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1)
733  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1)
734  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1)
735  stri(4)=2.d0*c(7)*ee(4)-beta(4)
736  stri(5)=2.d0*c(8)*ee(5)-beta(5)
737  stri(6)=2.d0*c(9)*ee(6)-beta(6)
738  endif
739 !
740 ! stress radius (only deviatoric part of stress enters)
741 !
742  strinv=(stri(1)+stri(2)+stri(3))/3.d0
743  do i=1,3
744  sg(i)=stri(i)-strinv
745  enddo
746  do i=4,6
747  sg(i)=stri(i)
748  enddo
749  dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+
750  & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6)))
751 !
752 ! evaluation of the yield surface
753 !
754  ec(1)=epqini
755  decra(1)=c0*dg
756  timeabq(1)=time
757  timeabq(2)=ttime+time
758  call creep(decra,deswa,xstateini(1,iint,iel),serd,ec,
759  & esw,p,svm,t1l,dtemp,predef,dpred,timeabq,dtime,
760  & amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt,
761  & kstep,kinc)
762  if(decra(1).gt.c0*dg) then
763  dg=decra(1)/c0
764  if(abs(iloop).gt.2) exitcriterion=1
765  endif
766 !
767 ! needed in case decra(1) was changed in subroutine creep,
768 ! for instance because it is too small
769 !
770  dg=decra(1)/c0
771 !
772  htri=dsg-c0*svm
773 !
774  do i=1,6
775  sg(i)=sg(i)/dsg
776  enddo
777 !
778 ! determining the residual matrix
779 !
780  do i=1,6
781  r(i)=ep0(i)-ep(i)+dg*sg(i)
782  enddo
783 !
784 ! check convergence
785 !
786  if(exitcriterion.eq.1) exit loop
787  if((dabs(htri).le.1.d-3).and.
788  & ((iloop.gt.2).and.((dabs(ddg).lt.1.d-10).or.
789  & (dabs(ddg).lt.1.d-3*dabs(dg))))) then
790  dd=0.d0
791  do i=1,6
792  dd=dd+r(i)*r(i)
793  enddo
794  dd=sqrt(dd)
795  if(dd.le.1.d-10) then
796  exit loop
797  endif
798  endif
799  if(iloop.gt.100) then
800 c write(*,*)
801 c & '*ERROR: no convergence in umat_aniso_creep'
802 c write(*,*) ' iloop>100'
803 c write(*,*) 'htri,dd ',htri,dd
804  exit loop
805  endif
806 !
807 ! determining b.x
808 !
809  b=dg/dsg
810 !
811  x(1)=b*(c1-sg(1)*sg(1))
812  x(2)=b*(c2-sg(1)*sg(2))
813  x(3)=b*(c1-sg(2)*sg(2))
814  x(4)=b*(c2-sg(1)*sg(3))
815  x(5)=b*(c2-sg(2)*sg(3))
816  x(6)=b*(c1-sg(3)*sg(3))
817  x(7)=-b*sg(1)*sg(4)
818  x(8)=-b*sg(2)*sg(4)
819  x(9)=-b*sg(3)*sg(4)
820  x(10)=b*(.5d0-sg(4)*sg(4))
821  x(11)=-b*sg(1)*sg(5)
822  x(12)=-b*sg(2)*sg(5)
823  x(13)=-b*sg(3)*sg(5)
824  x(14)=-b*sg(4)*sg(5)
825  x(15)=b*(.5d0-sg(5)*sg(5))
826  x(16)=-b*sg(1)*sg(6)
827  x(17)=-b*sg(2)*sg(6)
828  x(18)=-b*sg(3)*sg(6)
829  x(19)=-b*sg(4)*sg(6)
830  x(20)=-b*sg(5)*sg(6)
831  x(21)=b*(.5d0-sg(6)*sg(6))
832 !
833 ! filling the LHS
834 !
835  if(iorien.gt.0) then
836  gl(1,1)=cm1(1)+x(1)
837  gl(1,2)=cm1(2)+x(2)
838  gl(2,2)=cm1(3)+x(3)
839  gl(1,3)=cm1(4)+x(4)
840  gl(2,3)=cm1(5)+x(5)
841  gl(3,3)=cm1(6)+x(6)
842  gl(1,4)=cm1(7)+x(7)
843  gl(2,4)=cm1(8)+x(8)
844  gl(3,4)=cm1(9)+x(9)
845  gl(4,4)=cm1(10)+x(10)
846  gl(1,5)=cm1(11)+x(11)
847  gl(2,5)=cm1(12)+x(12)
848  gl(3,5)=cm1(13)+x(13)
849  gl(4,5)=cm1(14)+x(14)
850  gl(5,5)=cm1(15)+x(15)
851  gl(1,6)=cm1(16)+x(16)
852  gl(2,6)=cm1(17)+x(17)
853  gl(3,6)=cm1(18)+x(18)
854  gl(4,6)=cm1(19)+x(19)
855  gl(5,6)=cm1(20)+x(20)
856  gl(6,6)=cm1(21)+x(21)
857  do i=1,6
858  do j=1,i-1
859  gl(i,j)=gl(j,i)
860  enddo
861  enddo
862  else
863  gl(1,1)=cm1(1)+x(1)
864  gl(1,2)=cm1(2)+x(2)
865  gl(2,2)=cm1(3)+x(3)
866  gl(1,3)=cm1(4)+x(4)
867  gl(2,3)=cm1(5)+x(5)
868  gl(3,3)=cm1(6)+x(6)
869  gl(1,4)=x(7)
870  gl(2,4)=x(8)
871  gl(3,4)=x(9)
872  gl(4,4)=cm1(7)+x(10)
873  gl(1,5)=x(11)
874  gl(2,5)=x(12)
875  gl(3,5)=x(13)
876  gl(4,5)=x(14)
877  gl(5,5)=cm1(8)+x(15)
878  gl(1,6)=x(16)
879  gl(2,6)=x(17)
880  gl(3,6)=x(18)
881  gl(4,6)=x(19)
882  gl(5,6)=x(20)
883  gl(6,6)=cm1(9)+x(21)
884  do i=1,6
885  do j=1,i-1
886  gl(i,j)=gl(j,i)
887  enddo
888  enddo
889  endif
890 !
891 ! filling the RHS
892 !
893  do i=1,6
894  gr(i,1)=sg(i)
895  enddo
896 !
897 ! solve gl:(P:n)=gr
898 !
899  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
900  if(info.ne.0) then
901  write(*,*) '*ERROR in sc.f: linear equation solver'
902  write(*,*) ' exited with error: info = ',info
903  call exit(201)
904  endif
905 !
906  do i=1,6
907  pn(i)=gr(i,1)
908  enddo
909 !
910 ! calculating the creep contribution
911 !
912  gcreep=c1/decra(5)
913 !
914 ! calculating the correction to the consistency parameter
915 !
916  gm1=pn(1)*sg(1)+pn(2)*sg(2)+pn(3)*sg(3)+
917  & (pn(4)*sg(4)+pn(5)*sg(5)+pn(6)*sg(6))
918  gm1=1.d0/(gm1+gcreep)
919  fu=(htri-(pn(1)*r(1)+pn(2)*r(2)+pn(3)*r(3)+
920  & (pn(4)*r(4)+pn(5)*r(5)+pn(6)*r(6))))
921 !
922  if(iloop.eq.1) then
923 c write(*,*) 'iloop,dg,fu ',iloop,dg,fu
924  dg1=0.d0
925  fu1=fu
926  iloop=2
927  dg=1.d-10
928  ddg=dg
929  do i=1,6
930  ep1(i)=ep(i)
931  r1(i)=r(i)
932  sg1(i)=sg(i)
933  do j=1,6
934  gl1(i,j)=gl(i,j)
935  enddo
936  enddo
937  elseif((iloop.eq.2).or.(iloop.lt.0)) then
938  if(fu*fu1.lt.0.d0) then
939 c write(*,*) 'iloop,dg,fu ',iloop,dg,fu
940  if(iloop.eq.2) then
941  iloop=3
942  else
943  iloop=-iloop+1
944  endif
945  fu2=fu
946  dg2=dg
947  dg=(dg1+dg2)/2.d0
948  ddg=(dg2-dg1)/2.d0
949  do i=1,6
950  ep(i)=ep1(i)
951  r(i)=r1(i)
952  sg(i)=sg1(i)
953  do j=1,6
954  gl(i,j)=gl1(i,j)
955  enddo
956  enddo
957  else
958 c write(*,*) 'iloop,dg,fu ',iloop,dg,fu
959 c dg1=dg
960 c fu1=fu
961  if(iloop.eq.2) then
962  if(dabs(fu).gt.dabs(fu1)) exitcriterion=1
963  dg1=dg
964  fu1=fu
965  ddg=dg*9.d0
966  dg=dg*10.d0
967  else
968  dg1=dg
969  fu1=fu
970  dg=dg+ddg
971  iloop=iloop-1
972  endif
973  if(dg.gt.10.1d0) then
974  write(*,*)
975  & '*ERROR: no convergence in umat_aniso_creep'
976  write(*,*) ' dg>10.'
977  call exit(201)
978  endif
979  do i=1,6
980  ep1(i)=ep(i)
981  r1(i)=r(i)
982  sg1(i)=sg(i)
983  do j=1,6
984  gl1(i,j)=gl(i,j)
985  enddo
986  enddo
987  endif
988  else
989 c write(*,*) 'iloop,dg,fu ',iloop,dg,fu
990  if(fu*fu1.ge.0.d0) then
991  dg1=dg
992  fu1=fu
993  dg=(dg1+dg2)/2.d0
994  ddg=(dg2-dg1)/2.d0
995  do i=1,6
996  ep1(i)=ep(i)
997  r1(i)=r(i)
998  sg1(i)=sg(i)
999  do j=1,6
1000  gl1(i,j)=gl(i,j)
1001  enddo
1002  enddo
1003  iloop=-iloop-1
1004  else
1005  dg2=dg
1006  fu2=fu
1007  dg=(dg1+dg2)/2.d0
1008  ddg=(dg2-dg1)/2.d0
1009  do i=1,6
1010  ep(i)=ep1(i)
1011  r(i)=r1(i)
1012  sg(i)=sg1(i)
1013  do j=1,6
1014  gl(i,j)=gl1(i,j)
1015  enddo
1016  enddo
1017  iloop=iloop+1
1018  endif
1019  endif
1020 !
1021 ! updating the residual matrix
1022 !
1023  do i=1,6
1024  r(i)=r(i)+ddg*sg(i)
1025  enddo
1026 !
1027 ! update the plastic strain
1028 !
1029  gr(1,1)=r(1)
1030  gr(2,1)=r(2)
1031  gr(3,1)=r(3)
1032  gr(4,1)=r(4)
1033  gr(5,1)=r(5)
1034  gr(6,1)=r(6)
1035 !
1036  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,
1037  & info)
1038  if(info.ne.0) then
1039  write(*,*) '*ERROR in sc.f: linear equation solver'
1040  write(*,*) ' exited with error: info = ',info
1041  call exit(201)
1042  endif
1043 !
1044  if(iorien.gt.0) then
1045  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+
1046  & cm1(4)*gr(3,1)+
1047  & (cm1(7)*gr(4,1)+cm1(11)*gr(5,1)+
1048  & cm1(16)*gr(6,1))
1049  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+
1050  & cm1(5)*gr(3,1)+
1051  & (cm1(8)*gr(4,1)+cm1(12)*gr(5,1)+
1052  & cm1(17)*gr(6,1))
1053  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)
1054  & +cm1(6)*gr(3,1)+
1055  & (cm1(9)*gr(4,1)+cm1(13)*gr(5,1)+
1056  & cm1(18)*gr(6,1))
1057  ep(4)=ep(4)+cm1(7)*gr(1,1)+cm1(8)*gr(2,1)+
1058  & cm1(9)*gr(3,1)+
1059  & (cm1(10)*gr(4,1)+cm1(14)*gr(5,1)+
1060  & cm1(19)*gr(6,1))
1061  ep(5)=ep(5)+cm1(11)*gr(1,1)+cm1(12)*gr(2,1)+
1062  & cm1(13)*gr(3,1)+
1063  & (cm1(14)*gr(4,1)+cm1(15)*gr(5,1)+
1064  & cm1(20)*gr(6,1))
1065  ep(6)=ep(6)+cm1(16)*gr(1,1)+cm1(17)*gr(2,1)+
1066  & cm1(18)*gr(3,1)+
1067  & (cm1(19)*gr(4,1)+cm1(20)*gr(5,1)+
1068  & cm1(21)*gr(6,1))
1069  else
1070  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+
1071  & cm1(4)*gr(3,1)
1072  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+
1073  & cm1(5)*gr(3,1)
1074  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+
1075  & cm1(6)*gr(3,1)
1076  ep(4)=ep(4)+cm1(7)*gr(4,1)
1077  ep(5)=ep(5)+cm1(8)*gr(5,1)
1078  ep(6)=ep(6)+cm1(9)*gr(6,1)
1079  endif
1080 !
1081 ! end of major loop
1082 !
1083  enddo
1084 !
1085  endif
1086 !
1087  enddo loop
1088 !
1089 ! storing the stress
1090 !
1091  do i=1,6
1092  stre(i)=stri(i)
1093  enddo
1094 !
1095 ! calculating the tangent stiffness matrix
1096 !
1097  if(icmd.ne.3) then
1098 !
1099 ! determining p
1100 !
1101  gr(1,1)=1.d0
1102  gr(1,2)=0.
1103  gr(2,2)=1.d0
1104  gr(1,3)=0.
1105  gr(2,3)=0.
1106  gr(3,3)=1.d0
1107  gr(1,4)=0.
1108  gr(2,4)=0.
1109  gr(3,4)=0.
1110  gr(4,4)=1.d0
1111  gr(1,5)=0.
1112  gr(2,5)=0.
1113  gr(3,5)=0.
1114  gr(4,5)=0.
1115  gr(5,5)=1.d0
1116  gr(1,6)=0.
1117  gr(2,6)=0.
1118  gr(3,6)=0.
1119  gr(4,6)=0.
1120  gr(5,6)=0.
1121  gr(6,6)=1.d0
1122  do i=1,6
1123  do j=1,i-1
1124  gr(i,j)=gr(j,i)
1125  enddo
1126  enddo
1127  nrhs=6
1128 !
1129  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
1130  if(info.ne.0) then
1131  write(*,*) '*ERROR in sc.f: linear equation solver'
1132  write(*,*) ' exited with error: info = ',info
1133  call exit(201)
1134  endif
1135 !
1136  stiff(1)=gr(1,1)-gm1*pn(1)*pn(1)
1137  stiff(2)=gr(1,2)-gm1*pn(1)*pn(2)
1138  stiff(3)=gr(2,2)-gm1*pn(2)*pn(2)
1139  stiff(4)=gr(1,3)-gm1*pn(1)*pn(3)
1140  stiff(5)=gr(2,3)-gm1*pn(2)*pn(3)
1141  stiff(6)=gr(3,3)-gm1*pn(3)*pn(3)
1142  stiff(7)=(gr(1,4)-gm1*pn(1)*pn(4))/2.d0
1143  stiff(8)=(gr(2,4)-gm1*pn(2)*pn(4))/2.d0
1144  stiff(9)=(gr(3,4)-gm1*pn(3)*pn(4))/2.d0
1145  stiff(10)=(gr(4,4)-gm1*pn(4)*pn(4))/4.d0
1146  stiff(11)=(gr(1,5)-gm1*pn(1)*pn(5))/2.d0
1147  stiff(12)=(gr(2,5)-gm1*pn(2)*pn(5))/2.d0
1148  stiff(13)=(gr(3,5)-gm1*pn(3)*pn(5))/2.d0
1149  stiff(14)=(gr(4,5)-gm1*pn(4)*pn(5))/4.d0
1150  stiff(15)=(gr(5,5)-gm1*pn(5)*pn(5))/4.d0
1151  stiff(16)=(gr(1,6)-gm1*pn(1)*pn(6))/2.d0
1152  stiff(17)=(gr(2,6)-gm1*pn(2)*pn(6))/2.d0
1153  stiff(18)=(gr(3,6)-gm1*pn(3)*pn(6))/2.d0
1154  stiff(19)=(gr(4,6)-gm1*pn(4)*pn(6))/4.d0
1155  stiff(20)=(gr(5,6)-gm1*pn(5)*pn(6))/4.d0
1156  stiff(21)=(gr(6,6)-gm1*pn(6)*pn(6))/4.d0
1157  endif
1158 !
1159 ! updating the state variables
1160 !
1161  xstate(1,iint,iel)=epqini+c0*dg
1162  do i=1,6
1163  xstate(1+i,iint,iel)=ep(i)
1164  enddo
1165 !
1166 ! maximum equivalent viscoplastic strain in this increment
1167 !
1168  depvisc=max(depvisc,c0*dg)
1169 !
1170  return
max
#define max(a, b)
Definition: cascade.c:32
c1
static double * c1
Definition: mafillvcompmain.c:30
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
r1
static double * r1
Definition: filtermain.c:48
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
orthotropic
subroutine orthotropic(orthol, anisox)
Definition: orthotropic.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)