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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ umat_aniso_plas()

subroutine umat_aniso_plas ( 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 
)
23 !
24 ! calculates stiffness and stresses for a user defined material
25 ! law
26 !
27 ! icmd=3: calcutates 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
60 ! >0: thermal effects are taken into account (triggered
61 ! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
62 ! t1l temperature at the end of the increment
63 ! dtime time length of the increment
64 ! time step time at the end of the current increment
65 ! ttime total time at the start of the current step
66 !
67 ! icmd not equal to 3: calculate stress and stiffness
68 ! 3: calculate only stress
69 ! ielas 0: no elastic iteration: irreversible effects
70 ! are allowed
71 ! 1: elastic iteration, i.e. no irreversible
72 ! deformation allowed
73 !
74 ! mi(1) max. # of integration points per element in the
75 ! model
76 ! nstate_ max. # of state variables in the model
77 !
78 ! xstateini(nstate_,mi(1),# of elements)
79 ! state variables at the start of the increment
80 ! xstate(nstate_,mi(1),# of elements)
81 ! state variables at the end of the increment
82 !
83 ! stre(6) Piola-Kirchhoff stress of the second kind
84 ! at the start of the increment
85 !
86 ! iorien number of the local coordinate axis system
87 ! in the integration point at stake (takes the value
88 ! 0 if no local system applies)
89 ! pgauss(3) global coordinates of the integration point
90 ! orab(7,*) description of all local coordinate systems.
91 ! If a local coordinate system applies the global
92 ! tensors can be obtained by premultiplying the local
93 ! tensors with skl(3,3). skl is determined by calling
94 ! the subroutine transformatrix:
95 ! call transformatrix(orab(1,iorien),pgauss,skl)
96 !
97 !
98 ! OUTPUT:
99 !
100 ! xstate(nstate_,mi(1),# of elements)
101 ! updated state variables at the end of the increment
102 ! stre(6) Piola-Kirchhoff stress of the second kind at the
103 ! end of the increment
104 ! stiff(21): consistent tangent stiffness matrix in the material
105 ! frame of reference at the end of the increment. In
106 ! other words: the derivative of the PK2 stress with
107 ! respect to the Lagrangian strain tensor. The matrix
108 ! is supposed to be symmetric, only the upper half is
109 ! to be given in the same order as for a fully
110 ! anisotropic elastic material (*ELASTIC,TYPE=ANISO).
111 ! Notice that the matrix is an integral part of the
112 ! fourth order material tensor, i.e. the Voigt notation
113 ! is not used.
114 !
115  implicit none
116 !
117  integer visco
118 !
119  character*80 amat
120 !
121  integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien,
122  & i,j,ipiv(6),info,neq,lda,ldb,j1,j2,j3,j4,j5,j6,j7,j8,
123  & nrhs,iplas,kel(4,21),nmethod,ielastic,iloop
124 !
125  real*8 ep0(6),al10,al20(6),eeq,ep(6),al1,b,pn(6),qsn(6),
126  & al2(6),dg,ddg,ca,cn,c(21),r0,x(21),cm1(21),h1,h2,
127  & q1,q2(6),stri(6),htri,sg(6),r(13),au1(21),au2(21),
128  & ee(6),dd,gl(6,6),gr(6,6),c0,c1,c2,c3,c4,c5,c6,pnewdt,
129  & skl(3,3),gcreep,gm1,ya(3,3,3,3),d1,d2,dsg,detc,strinv,
130  & elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
131  & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
132  & time,ttime,xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
133 !
134  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,
135  & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,
136  & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,
137  & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/))
138 !
139  iloop=0
140 !
141  c0=dsqrt(2.d0/3.d0)
142  c1=2.d0/3.d0
143  c2=-1.d0/3.d0
144 !
145 ! elastic constants
146 !
147  if(iorien.gt.0) then
148 !
149  call transformatrix(orab(1,iorien),pgauss,skl)
150 !
151  call orthotropic(elconloc,ya)
152 !
153  do j=1,21
154  j1=kel(1,j)
155  j2=kel(2,j)
156  j3=kel(3,j)
157  j4=kel(4,j)
158  c(j)=0.d0
159  do j5=1,3
160  do j6=1,3
161  do j7=1,3
162  do j8=1,3
163  c(j)=c(j)+ya(j5,j6,j7,j8)*
164  & skl(j1,j5)*skl(j2,j6)*skl(j3,j7)*skl(j4,j8)
165  enddo
166  enddo
167  enddo
168  enddo
169  enddo
170 !
171  else
172  do i=1,9
173  c(i)=elconloc(i)
174  enddo
175  endif
176 !
177 ! state variables
178 !
179 ! equivalent plastic strain
180 !
181  eeq=xstateini(1,iint,iel)
182 !
183 ! plastic strain
184 !
185  do i=1,6
186  ep0(i)=xstateini(1+i,iint,iel)
187  enddo
188 !
189 ! isotropic hardening variable
190 !
191  al10=xstateini(8,iint,iel)
192 !
193 ! kinematic hardening variable
194 !
195  do i=1,6
196  al20(i)=xstateini(8+i,iint,iel)
197  enddo
198 !
199 ! elastic strains
200 !
201  do i=1,6
202  ee(i)=emec(i)-ep0(i)
203  enddo
204 !
205 ! (visco)plastic constants
206 !
207  r0=elconloc(10)
208  d1=elconloc(11)
209  d2=elconloc(12)
210  ca=c0/(elconloc(13)*(ttime+time-dtime)**elconloc(15)*dtime)
211  cn=elconloc(14)
212 !
213  if((ca.lt.0.d0).or.((nmethod.eq.1).and.(ithermal.ne.3))) then
214  visco=0
215  else
216  visco=1
217  endif
218 !
219 ! if the user did not activate a viscous calculation, a pure
220 ! creep calculation (without plasticity) is reduced to an
221 ! elastic calculation
222 !
223  ielastic=0
224  if(visco.eq.0) then
225  if((dabs(r0).lt.1.d-10).and.
226  & (dabs(d1).lt.1.d-10).and.
227  & (dabs(d2).lt.1.d-10)) ielastic=1
228  endif
229 !
230  h1=d1
231  h2=2.d0*d2/3.d0
232 !
233 ! stress state variables q1 and q2
234 !
235  q1=-d1*al10
236  do i=1,6
237  q2(i)=-d2*al20(i)
238  enddo
239 !
240 ! global trial stress tensor
241 !
242  if(iorien.gt.0) then
243  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+
244  & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6))
245  & -beta(1)
246  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+
247  & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6))
248  & -beta(2)
249  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+
250  & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6))
251  & -beta(3)
252  stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+
253  & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6))
254  & -beta(4)
255  stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+
256  & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6))
257  & -beta(5)
258  stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+
259  & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6))
260  & -beta(6)
261  else
262  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1)
263  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1)
264  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1)
265  stri(4)=2.d0*c(7)*ee(4)-beta(4)
266  stri(5)=2.d0*c(8)*ee(5)-beta(5)
267  stri(6)=2.d0*c(9)*ee(6)-beta(6)
268  endif
269 !
270 ! stress radius (only deviatoric part of stress enters)
271 !
272  strinv=(stri(1)+stri(2)+stri(3))/3.d0
273  do i=1,3
274  sg(i)=stri(i)-strinv+q2(i)
275  enddo
276  do i=4,6
277  sg(i)=stri(i)+q2(i)
278  enddo
279  dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+
280  & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6)))
281 !
282 ! evaluation of the yield surface
283 !
284  htri=dsg+c0*(q1-r0)
285 !
286 ! check whether plasticity occurs
287 !
288  if(htri.gt.0.d0) then
289  iplas=1
290  else
291  iplas=0
292  endif
293 !
294  if((iplas.eq.0).or.(ielas.eq.1).or.(ielastic.eq.1)) then
295 !
296 ! elastic stress
297 !
298  do i=1,6
299  stre(i)=stri(i)
300  enddo
301 !
302 ! updating the state variables
303 !
304  xstate(1,iint,iel)=eeq
305  do i=1,6
306  xstate(1+i,iint,iel)=ep0(i)
307  enddo
308  xstate(8,iint,iel)=al10
309  do i=1,6
310  xstate(8+i,iint,iel)=al20(i)
311  enddo
312 !
313 ! elastic stiffness
314 !
315  if(icmd.ne.3) then
316  if(iorien.gt.0) then
317  do i=1,21
318  stiff(i)=c(i)
319  enddo
320  else
321  stiff(1)=c(1)
322  stiff(2)=c(2)
323  stiff(3)=c(3)
324  stiff(4)=c(4)
325  stiff(5)=c(5)
326  stiff(6)=c(6)
327  stiff(7)=0.d0
328  stiff(8)=0.d0
329  stiff(9)=0.d0
330  stiff(10)=c(7)
331  stiff(11)=0.d0
332  stiff(12)=0.d0
333  stiff(13)=0.d0
334  stiff(14)=0.d0
335  stiff(15)=c(8)
336  stiff(16)=0.d0
337  stiff(17)=0.d0
338  stiff(18)=0.d0
339  stiff(19)=0.d0
340  stiff(20)=0.d0
341  stiff(21)=c(9)
342  endif
343  endif
344 !
345  return
346  endif
347 !
348 ! plastic deformation
349 !
350  neq=6
351  nrhs=1
352  lda=6
353  ldb=6
354 !
355 ! initializing the state variables
356 !
357  do i=1,6
358  ep(i)=ep0(i)
359  enddo
360  al1=al10
361  do i=1,6
362  al2(i)=al20(i)
363  enddo
364  dg=0.d0
365  ddg=0.d0
366 !
367 ! determining the inverse of c
368 !
369  if(iorien.gt.0) then
370 !
371 ! solve gl:C=gr
372 !
373  gl(1,1)=c(1)
374  gl(1,2)=c(2)
375  gl(2,2)=c(3)
376  gl(1,3)=c(4)
377  gl(2,3)=c(5)
378  gl(3,3)=c(6)
379  gl(1,4)=c(7)
380  gl(2,4)=c(8)
381  gl(3,4)=c(9)
382  gl(4,4)=c(10)
383  gl(1,5)=c(11)
384  gl(2,5)=c(12)
385  gl(3,5)=c(13)
386  gl(4,5)=c(14)
387  gl(5,5)=c(15)
388  gl(1,6)=c(16)
389  gl(2,6)=c(17)
390  gl(3,6)=c(18)
391  gl(4,6)=c(19)
392  gl(5,6)=c(20)
393  gl(6,6)=c(21)
394  do i=1,6
395  do j=1,i-1
396  gl(i,j)=gl(j,i)
397  enddo
398  enddo
399  do i=1,6
400  do j=4,6
401  gl(i,j)=2.d0*gl(i,j)
402  enddo
403  enddo
404  do i=1,6
405  do j=1,6
406  gr(i,j)=0.d0
407  enddo
408  if(i.le.3) then
409  gr(i,i)=1.d0
410  else
411  gr(i,i)=0.5d0
412  endif
413  enddo
414  nrhs=6
415  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
416  if(info.ne.0) then
417  write(*,*) '*ERROR in sc.f: linear equation solver'
418  write(*,*) ' exited with error: info = ',info
419  call exit(201)
420  endif
421  nrhs=1
422  cm1(1)=gr(1,1)
423  cm1(2)=gr(1,2)
424  cm1(3)=gr(2,2)
425  cm1(4)=gr(1,3)
426  cm1(5)=gr(2,3)
427  cm1(6)=gr(3,3)
428  cm1(7)=gr(1,4)
429  cm1(8)=gr(2,4)
430  cm1(9)=gr(3,4)
431  cm1(10)=gr(4,4)
432  cm1(11)=gr(1,5)
433  cm1(12)=gr(2,5)
434  cm1(13)=gr(3,5)
435  cm1(14)=gr(4,5)
436  cm1(15)=gr(5,5)
437  cm1(16)=gr(1,6)
438  cm1(17)=gr(2,6)
439  cm1(18)=gr(3,6)
440  cm1(19)=gr(4,6)
441  cm1(20)=gr(5,6)
442  cm1(21)=gr(6,6)
443  else
444  detc=c(1)*(c(3)*c(6)-c(5)*c(5))-
445  & c(2)*(c(2)*c(6)-c(4)*c(5))+
446  & c(4)*(c(2)*c(5)-c(4)*c(3))
447  cm1(1)=(c(3)*c(6)-c(5)*c(5))/detc
448  cm1(2)=(c(5)*c(4)-c(2)*c(6))/detc
449  cm1(3)=(c(1)*c(6)-c(4)*c(4))/detc
450  cm1(4)=(c(2)*c(5)-c(3)*c(4))/detc
451  cm1(5)=(c(2)*c(4)-c(1)*c(5))/detc
452  cm1(6)=(c(1)*c(3)-c(2)*c(2))/detc
453  cm1(7)=1.d0/(4.d0*c(7))
454  cm1(8)=1.d0/(4.d0*c(8))
455  cm1(9)=1.d0/(4.d0*c(9))
456  endif
457 !
458 ! loop
459 !
460  do
461 !
462  iloop=iloop+1
463  if(iloop.gt.100) then
464 c NOTE: write statements cause problems for
465 c parallellized execution
466 c write(*,*) '*WARNING in umat_aniso_plas: material loop'
467 c write(*,*) ' did not converge in integration'
468 c write(*,*) ' point',iint,'in element',iel,';'
469 c write(*,*) ' the increment size is reduced'
470 c write(*,*)
471  pnewdt=0.25d0
472  return
473  endif
474 !
475 ! elastic strains
476 !
477  do i=1,6
478  ee(i)=emec(i)-ep(i)
479  enddo
480 !
481 ! stress state variables q1 and q2
482 !
483  q1=-d1*al1
484  do i=1,6
485  q2(i)=-d2*al2(i)
486  enddo
487 !
488 ! global trial stress tensor
489 !
490  if(iorien.gt.0) then
491  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)+
492  & 2.d0*(c(7)*ee(4)+c(11)*ee(5)+c(16)*ee(6))
493  & -beta(1)
494  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)+
495  & 2.d0*(c(8)*ee(4)+c(12)*ee(5)+c(17)*ee(6))
496  & -beta(2)
497  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)+
498  & 2.d0*(c(9)*ee(4)+c(13)*ee(5)+c(18)*ee(6))
499  & -beta(3)
500  stri(4)=c(7)*ee(1)+c(8)*ee(2)+c(9)*ee(3)+
501  & 2.d0*(c(10)*ee(4)+c(14)*ee(5)+c(19)*ee(6))
502  & -beta(4)
503  stri(5)=c(11)*ee(1)+c(12)*ee(2)+c(13)*ee(3)+
504  & 2.d0*(c(14)*ee(4)+c(15)*ee(5)+c(20)*ee(6))
505  & -beta(5)
506  stri(6)=c(16)*ee(1)+c(17)*ee(2)+c(18)*ee(3)+
507  & 2.d0*(c(19)*ee(4)+c(20)*ee(5)+c(21)*ee(6))
508  & -beta(6)
509  else
510  stri(1)=c(1)*ee(1)+c(2)*ee(2)+c(4)*ee(3)-beta(1)
511  stri(2)=c(2)*ee(1)+c(3)*ee(2)+c(5)*ee(3)-beta(1)
512  stri(3)=c(4)*ee(1)+c(5)*ee(2)+c(6)*ee(3)-beta(1)
513  stri(4)=2.d0*c(7)*ee(4)-beta(4)
514  stri(5)=2.d0*c(8)*ee(5)-beta(5)
515  stri(6)=2.d0*c(9)*ee(6)-beta(6)
516  endif
517 !
518 ! stress radius (only deviatoric part of stress enters)
519 !
520  strinv=(stri(1)+stri(2)+stri(3))/3.d0
521  do i=1,3
522  sg(i)=stri(i)-strinv+q2(i)
523  enddo
524  do i=4,6
525  sg(i)=stri(i)+q2(i)
526  enddo
527  dsg=dsqrt(sg(1)*sg(1)+sg(2)*sg(2)+sg(3)*sg(3)+
528  & 2.d0*(sg(4)*sg(4)+sg(5)*sg(5)+sg(6)*sg(6)))
529 !
530 ! evaluation of the yield surface
531 !
532  if(visco.eq.1) then
533  htri=dsg+c0*(q1-r0-(ca*dg)**(1.d0/cn))
534  else
535  htri=dsg+c0*(q1-r0)
536  endif
537 !
538  do i=1,6
539  sg(i)=sg(i)/dsg
540  enddo
541 !
542 ! determining the residual matrix
543 !
544  do i=1,6
545  r(i)=ep0(i)-ep(i)+dg*sg(i)
546  enddo
547  r(7)=al10-al1+dg*c0
548  do i=1,6
549  r(7+i)=al20(i)-al2(i)+dg*sg(i)
550  enddo
551 !
552 ! check convergence
553 !
554  if((htri.le.1.d-5).or.(dabs(ddg).lt.1.d-3*dabs(dg))) then
555  dd=0.d0
556  do i=1,13
557  dd=dd+r(i)*r(i)
558  enddo
559  dd=sqrt(dd)
560  if(dd.le.1.d-10) then
561  exit
562  endif
563  endif
564 !
565 ! determining b.x
566 !
567  b=dg/dsg
568 !
569  x(1)=b*(c1-sg(1)*sg(1))
570  x(2)=b*(c2-sg(1)*sg(2))
571  x(3)=b*(c1-sg(2)*sg(2))
572  x(4)=b*(c2-sg(1)*sg(3))
573  x(5)=b*(c2-sg(2)*sg(3))
574  x(6)=b*(c1-sg(3)*sg(3))
575  x(7)=-b*sg(1)*sg(4)
576  x(8)=-b*sg(2)*sg(4)
577  x(9)=-b*sg(3)*sg(4)
578  x(10)=b*(.5d0-sg(4)*sg(4))
579  x(11)=-b*sg(1)*sg(5)
580  x(12)=-b*sg(2)*sg(5)
581  x(13)=-b*sg(3)*sg(5)
582  x(14)=-b*sg(4)*sg(5)
583  x(15)=b*(.5d0-sg(5)*sg(5))
584  x(16)=-b*sg(1)*sg(6)
585  x(17)=-b*sg(2)*sg(6)
586  x(18)=-b*sg(3)*sg(6)
587  x(19)=-b*sg(4)*sg(6)
588  x(20)=-b*sg(5)*sg(6)
589  x(21)=b*(.5d0-sg(6)*sg(6))
590 !
591  do i=1,21
592  au1(i)=h2*x(i)
593  enddo
594  au1(1)=au1(1)+1.d0
595  au1(3)=au1(3)+1.d0
596  au1(6)=au1(6)+1.d0
597  au1(10)=au1(10)+.5d0
598  au1(15)=au1(15)+.5d0
599  au1(21)=au1(21)+.5d0
600 !
601 ! filling the LHS
602 !
603  if(iorien.gt.0) then
604  gl(1,1)=au1(1)*cm1(1)+au1(2)*cm1(2)+au1(4)*cm1(4)+
605  & 2.d0*(au1(7)*cm1(7)+au1(11)*cm1(11)+au1(16)*cm1(16))+
606  & x(1)
607  gl(1,2)=au1(1)*cm1(2)+au1(2)*cm1(3)+au1(4)*cm1(5)+
608  & 2.d0*(au1(7)*cm1(8)+au1(11)*cm1(12)+au1(16)*cm1(17))+
609  & x(2)
610  gl(2,2)=au1(2)*cm1(2)+au1(3)*cm1(3)+au1(5)*cm1(5)+
611  & 2.d0*(au1(8)*cm1(8)+au1(12)*cm1(12)+au1(17)*cm1(17))+
612  & x(3)
613  gl(1,3)=au1(1)*cm1(4)+au1(2)*cm1(5)+au1(4)*cm1(6)+
614  & 2.d0*(au1(7)*cm1(9)+au1(11)*cm1(13)+au1(16)*cm1(18))+
615  & x(4)
616  gl(2,3)=au1(2)*cm1(4)+au1(3)*cm1(5)+au1(5)*cm1(6)+
617  & 2.d0*(au1(8)*cm1(9)+au1(12)*cm1(13)+au1(17)*cm1(18))+
618  & x(5)
619  gl(3,3)=au1(4)*cm1(4)+au1(5)*cm1(5)+au1(6)*cm1(6)+
620  & 2.d0*(au1(9)*cm1(9)+au1(13)*cm1(13)+au1(18)*cm1(18))+
621  & x(6)
622  gl(1,4)=au1(1)*cm1(7)+au1(2)*cm1(8)+au1(4)*cm1(9)+
623  & 2.d0*(au1(7)*cm1(10)+au1(11)*cm1(14)+au1(16)*cm1(19))+
624  & x(7)
625  gl(2,4)=au1(2)*cm1(7)+au1(3)*cm1(8)+au1(5)*cm1(9)+
626  & 2.d0*(au1(8)*cm1(10)+au1(12)*cm1(14)+au1(17)*cm1(19))+
627  & x(8)
628  gl(3,4)=au1(4)*cm1(7)+au1(5)*cm1(8)+au1(6)*cm1(9)+
629  & 2.d0*(au1(9)*cm1(10)+au1(13)*cm1(14)+au1(18)*cm1(19))+
630  & x(9)
631  gl(4,4)=au1(7)*cm1(7)+au1(8)*cm1(8)+au1(9)*cm1(9)+
632  & 2.d0*(au1(10)*cm1(10)+au1(14)*cm1(14)+au1(19)*cm1(19))+
633  & x(10)
634  gl(1,5)=au1(1)*cm1(11)+au1(2)*cm1(12)+au1(4)*cm1(13)+
635  & 2.d0*(au1(7)*cm1(14)+au1(11)*cm1(15)+au1(16)*cm1(20))+
636  & x(11)
637  gl(2,5)=au1(2)*cm1(11)+au1(3)*cm1(12)+au1(5)*cm1(13)+
638  & 2.d0*(au1(8)*cm1(14)+au1(12)*cm1(15)+au1(17)*cm1(20))+
639  & x(12)
640  gl(3,5)=au1(4)*cm1(11)+au1(5)*cm1(12)+au1(6)*cm1(13)+
641  & 2.d0*(au1(9)*cm1(14)+au1(13)*cm1(15)+au1(18)*cm1(20))+
642  & x(13)
643  gl(4,5)=au1(7)*cm1(11)+au1(8)*cm1(12)+au1(9)*cm1(13)+
644  & 2.d0*(au1(10)*cm1(14)+au1(14)*cm1(15)+au1(19)*cm1(20))+
645  & x(14)
646  gl(5,5)=au1(11)*cm1(11)+au1(12)*cm1(12)+au1(13)*cm1(13)+
647  & 2.d0*(au1(14)*cm1(14)+au1(15)*cm1(15)+au1(20)*cm1(20))+
648  & x(15)
649  gl(1,6)=au1(1)*cm1(16)+au1(2)*cm1(17)+au1(4)*cm1(18)+
650  & 2.d0*(au1(7)*cm1(19)+au1(11)*cm1(20)+au1(16)*cm1(21))+
651  & x(16)
652  gl(2,6)=au1(2)*cm1(16)+au1(3)*cm1(17)+au1(5)*cm1(18)+
653  & 2.d0*(au1(8)*cm1(19)+au1(12)*cm1(20)+au1(17)*cm1(21))+
654  & x(17)
655  gl(3,6)=au1(4)*cm1(16)+au1(5)*cm1(17)+au1(6)*cm1(18)+
656  & 2.d0*(au1(9)*cm1(19)+au1(13)*cm1(20)+au1(18)*cm1(21))+
657  & x(18)
658  gl(4,6)=au1(7)*cm1(16)+au1(8)*cm1(17)+au1(9)*cm1(18)+
659  & 2.d0*(au1(10)*cm1(19)+au1(14)*cm1(20)+au1(19)*cm1(21))+
660  & x(19)
661  gl(5,6)=au1(11)*cm1(16)+au1(12)*cm1(17)+au1(13)*cm1(18)+
662  & 2.d0*(au1(14)*cm1(19)+au1(15)*cm1(20)+au1(20)*cm1(21))+
663  & x(20)
664  gl(6,6)=au1(16)*cm1(16)+au1(17)*cm1(17)+au1(18)*cm1(18)+
665  & 2.d0*(au1(19)*cm1(19)+au1(20)*cm1(20)+au1(21)*cm1(21))+
666  & x(21)
667  do i=1,6
668  do j=1,i-1
669  gl(i,j)=gl(j,i)
670  enddo
671  enddo
672  do i=1,6
673  do j=4,6
674  gl(i,j)=2.d0*gl(i,j)
675  enddo
676  enddo
677  else
678  gl(1,1)=au1(1)*cm1(1)+au1(2)*cm1(2)+au1(4)*cm1(4)+x(1)
679  gl(1,2)=au1(1)*cm1(2)+au1(2)*cm1(3)+au1(4)*cm1(5)+x(2)
680  gl(2,2)=au1(2)*cm1(2)+au1(3)*cm1(3)+au1(5)*cm1(5)+x(3)
681  gl(1,3)=au1(1)*cm1(4)+au1(2)*cm1(5)+au1(4)*cm1(6)+x(4)
682  gl(2,3)=au1(2)*cm1(4)+au1(3)*cm1(5)+au1(5)*cm1(6)+x(5)
683  gl(3,3)=au1(4)*cm1(4)+au1(5)*cm1(5)+au1(6)*cm1(6)+x(6)
684  gl(1,4)=2.d0*au1(7)*cm1(7)+x(7)
685  gl(2,4)=2.d0*au1(8)*cm1(7)+x(8)
686  gl(3,4)=2.d0*au1(9)*cm1(7)+x(9)
687  gl(4,4)=2.d0*au1(10)*cm1(7)+x(10)
688  gl(1,5)=2.d0*au1(11)*cm1(8)+x(11)
689  gl(2,5)=2.d0*au1(12)*cm1(8)+x(12)
690  gl(3,5)=2.d0*au1(13)*cm1(8)+x(13)
691  gl(4,5)=2.d0*au1(14)*cm1(8)+x(14)
692  gl(5,5)=2.d0*au1(15)*cm1(8)+x(15)
693  gl(1,6)=2.d0*au1(16)*cm1(9)+x(16)
694  gl(2,6)=2.d0*au1(17)*cm1(9)+x(17)
695  gl(3,6)=2.d0*au1(18)*cm1(9)+x(18)
696  gl(4,6)=2.d0*au1(19)*cm1(9)+x(19)
697  gl(5,6)=2.d0*au1(20)*cm1(9)+x(20)
698  gl(6,6)=2.d0*au1(21)*cm1(9)+x(21)
699  do i=1,6
700  do j=1,i-1
701  gl(i,j)=gl(j,i)
702  enddo
703  enddo
704  do i=1,6
705  do j=4,6
706  gl(i,j)=2.d0*gl(i,j)
707  enddo
708  enddo
709  endif
710 !
711 ! filling the RHS
712 !
713  do i=1,6
714  gr(i,1)=sg(i)
715  enddo
716 !
717 ! solve gl:(P:n)=gr
718 !
719  call dgesv(neq,nrhs,gl,lda,ipiv,gr,ldb,info)
720  if(info.ne.0) then
721  write(*,*) '*ERROR in sc.f: linear equation solver'
722  write(*,*) ' exited with error: info = ',info
723  call exit(201)
724  endif
725 !
726  do i=1,6
727  pn(i)=gr(i,1)
728  enddo
729 !
730  c3=-h2/(1.d0+b*h2)
731  qsn(1)=c3*(x(1)*pn(1)+x(2)*pn(2)+x(4)*pn(3)+
732  & 2.d0*(x(7)*pn(4)+x(11)*pn(5)+x(16)*pn(6)))+sg(1)*h2
733  qsn(2)=c3*(x(2)*pn(1)+x(3)*pn(2)+x(5)*pn(3)+
734  & 2.d0*(x(8)*pn(4)+x(12)*pn(5)+x(17)*pn(6)))+sg(2)*h2
735  qsn(3)=c3*(x(4)*pn(1)+x(5)*pn(2)+x(6)*pn(3)+
736  & 2.d0*(x(9)*pn(4)+x(13)*pn(5)+x(18)*pn(6)))+sg(3)*h2
737  qsn(4)=c3*(x(7)*pn(1)+x(8)*pn(2)+x(9)*pn(3)+
738  & 2.d0*(x(10)*pn(4)+x(14)*pn(5)+x(19)*pn(6)))+sg(4)*h2
739  qsn(5)=c3*(x(11)*pn(1)+x(12)*pn(2)+x(13)*pn(3)+
740  & 2.d0*(x(14)*pn(4)+x(15)*pn(5)+x(20)*pn(6)))+sg(5)*h2
741  qsn(6)=c3*(x(16)*pn(1)+x(17)*pn(2)+x(18)*pn(3)+
742  & 2.d0*(x(19)*pn(4)+x(20)*pn(5)+x(21)*pn(6)))+sg(6)*h2
743 !
744 ! calculating the creep contribution
745 !
746  if(visco.eq.1) then
747  if(dg.gt.0.d0) then
748  gcreep=c0*ca/cn*(dg*ca)**(1.d0/cn-1.d0)
749  else
750 !
751 ! for gamma ein default of 1.d-10 is taken to
752 ! obtain a finite gradient
753 !
754  gcreep=c0*ca/cn*(1.d-10*ca)**(1.d0/cn-1.d0)
755  endif
756  endif
757 !
758 ! calculating the correction to the consistency parameter
759 !
760  gm1=pn(1)*sg(1)+pn(2)*sg(2)+pn(3)*sg(3)+
761  & 2.d0*(pn(4)*sg(4)+pn(5)*sg(5)+pn(6)*sg(6))+
762  & c1*h1+
763  & qsn(1)*sg(1)+qsn(2)*sg(2)+qsn(3)*sg(3)+
764  & 2.d0*(qsn(4)*sg(4)+qsn(5)*sg(5)+qsn(6)*sg(6))
765  if(visco.eq.1) then
766  gm1=1.d0/(gm1+gcreep)
767  else
768  gm1=1.d0/gm1
769  endif
770  ddg=gm1*(htri-(pn(1)*r(1)+pn(2)*r(2)+pn(3)*r(3)+
771  & 2.d0*(pn(4)*r(4)+pn(5)*r(5)+pn(6)*r(6))+
772  & c0*h1*r(7)+
773  & qsn(1)*r(8)+qsn(2)*r(9)+qsn(3)*r(10)+
774  & 2.d0*(qsn(4)*r(11)+qsn(5)*r(12)+qsn(6)*r(13))))
775 !
776 ! updating the residual matrix
777 !
778  do i=1,6
779  r(i)=r(i)+ddg*sg(i)
780  enddo
781  r(7)=r(7)+ddg*c0
782  do i=1,6
783  r(7+i)=r(7+i)+ddg*sg(i)
784  enddo
785 !
786 ! update the plastic strain
787 !
788  gr(1,1)=au1(1)*r(1)+au1(2)*r(2)+au1(4)*r(3)+
789  & 2.d0*(au1(7)*r(4)+au1(11)*r(5)+au1(16)*r(6))
790  & -h2*(x(1)*r(8)+x(2)*r(9)+x(4)*r(10)+
791  & 2.d0*(x(7)*r(11)+x(11)*r(12)+x(16)*r(13)))
792  gr(2,1)=au1(2)*r(1)+au1(3)*r(2)+au1(5)*r(3)+
793  & 2.d0*(au1(8)*r(4)+au1(12)*r(5)+au1(17)*r(6))
794  & -h2*(x(2)*r(8)+x(3)*r(9)+x(5)*r(10)+
795  & 2.d0*(x(8)*r(11)+x(12)*r(12)+x(17)*r(13)))
796  gr(3,1)=au1(4)*r(1)+au1(5)*r(2)+au1(6)*r(3)+
797  & 2.d0*(au1(9)*r(4)+au1(13)*r(5)+au1(18)*r(6))
798  & -h2*(x(4)*r(8)+x(5)*r(9)+x(6)*r(10)+
799  & 2.d0*(x(9)*r(11)+x(13)*r(12)+x(18)*r(13)))
800  gr(4,1)=au1(7)*r(1)+au1(8)*r(2)+au1(9)*r(3)+
801  & 2.d0*(au1(10)*r(4)+au1(14)*r(5)+au1(19)*r(6))
802  & -h2*(x(7)*r(8)+x(8)*r(9)+x(9)*r(10)+
803  & 2.d0*(x(10)*r(11)+x(14)*r(12)+x(19)*r(13)))
804  gr(5,1)=au1(11)*r(1)+au1(12)*r(2)+au1(13)*r(3)+
805  & 2.d0*(au1(14)*r(4)+au1(15)*r(5)+au1(20)*r(6))
806  & -h2*(x(11)*r(8)+x(12)*r(9)+x(13)*r(10)+
807  & 2.d0*(x(14)*r(11)+x(15)*r(12)+x(20)*r(13)))
808  gr(6,1)=au1(16)*r(1)+au1(17)*r(2)+au1(18)*r(3)+
809  & 2.d0*(au1(19)*r(4)+au1(20)*r(5)+au1(21)*r(6))
810  & -h2*(x(16)*r(8)+x(17)*r(9)+x(18)*r(10)+
811  & 2.d0*(x(19)*r(11)+x(20)*r(12)+x(21)*r(13)))
812 !
813  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
814  if(info.ne.0) then
815  write(*,*) '*ERROR in sc.f: linear equation solver'
816  write(*,*) ' exited with error: info = ',info
817  call exit(201)
818  endif
819 !
820  if(iorien.gt.0) then
821  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1)+
822  & 2.d0*(cm1(7)*gr(4,1)+cm1(11)*gr(5,1)+cm1(16)*gr(6,1))
823  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1)+
824  & 2.d0*(cm1(8)*gr(4,1)+cm1(12)*gr(5,1)+cm1(17)*gr(6,1))
825  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1)+
826  & 2.d0*(cm1(9)*gr(4,1)+cm1(13)*gr(5,1)+cm1(18)*gr(6,1))
827  ep(4)=ep(4)+cm1(7)*gr(1,1)+cm1(8)*gr(2,1)+cm1(9)*gr(3,1)+
828  & 2.d0*(cm1(10)*gr(4,1)+cm1(14)*gr(5,1)+cm1(19)*gr(6,1))
829  ep(5)=ep(5)+cm1(11)*gr(1,1)+cm1(12)*gr(2,1)+cm1(13)*gr(3,1)+
830  & 2.d0*(cm1(14)*gr(4,1)+cm1(15)*gr(5,1)+cm1(20)*gr(6,1))
831  ep(6)=ep(6)+cm1(16)*gr(1,1)+cm1(17)*gr(2,1)+cm1(18)*gr(3,1)+
832  & 2.d0*(cm1(19)*gr(4,1)+cm1(20)*gr(5,1)+cm1(21)*gr(6,1))
833  else
834  ep(1)=ep(1)+cm1(1)*gr(1,1)+cm1(2)*gr(2,1)+cm1(4)*gr(3,1)
835  ep(2)=ep(2)+cm1(2)*gr(1,1)+cm1(3)*gr(2,1)+cm1(5)*gr(3,1)
836  ep(3)=ep(3)+cm1(4)*gr(1,1)+cm1(5)*gr(2,1)+cm1(6)*gr(3,1)
837  ep(4)=ep(4)+2.d0*cm1(7)*gr(4,1)
838  ep(5)=ep(5)+2.d0*cm1(8)*gr(5,1)
839  ep(6)=ep(6)+2.d0*cm1(9)*gr(6,1)
840  endif
841 !
842 ! update the isotropic hardening variable
843 !
844  al1=al1+r(7)
845 !
846 ! update the kinematic hardening variables
847 !
848  c4=1.d0/(1.d0+b*h2)
849  c6=c4*b*h2
850  c5=c6/3.d0
851  au2(1)=c4+c5+c6*sg(1)*sg(1)
852  au2(2)=c5+c6*sg(1)*sg(2)
853  au2(3)=c4+c5+c6*sg(2)*sg(2)
854  au2(4)=c5+c6*sg(1)*sg(3)
855  au2(5)=c5+c6*sg(2)*sg(3)
856  au2(6)=c4+c5+c6*sg(3)*sg(3)
857  au2(7)=c6*sg(1)*sg(4)
858  au2(8)=c6*sg(2)*sg(4)
859  au2(9)=c6*sg(3)*sg(4)
860  au2(10)=c4/2.d0+c6*sg(4)*sg(4)
861  au2(11)=c6*sg(1)*sg(5)
862  au2(12)=c6*sg(2)*sg(5)
863  au2(13)=c6*sg(3)*sg(5)
864  au2(14)=c6*sg(4)*sg(5)
865  au2(15)=c4/2.d0+c6*sg(5)*sg(5)
866  au2(16)=c6*sg(1)*sg(6)
867  au2(17)=c6*sg(2)*sg(6)
868  au2(18)=c6*sg(3)*sg(6)
869  au2(19)=c6*sg(4)*sg(6)
870  au2(20)=c6*sg(5)*sg(6)
871  au2(21)=c4/2.d0+c6*sg(6)*sg(6)
872 !
873  al2(1)=al2(1)+au2(1)*r(8)+au2(2)*r(9)+au2(4)*r(10)+
874  & 2.d0*(au2(7)*r(11)+au2(11)*r(12)+au2(16)*r(13))
875  & -c4*(x(1)*gr(1,1)+x(2)*gr(2,1)+x(4)*gr(3,1)+
876  & 2.d0*(x(7)*gr(4,1)+x(11)*gr(5,1)+x(16)*gr(6,1)))
877  al2(2)=al2(2)+au2(2)*r(8)+au2(3)*r(9)+au2(5)*r(10)+
878  & 2.d0*(au2(8)*r(11)+au2(12)*r(12)+au2(17)*r(13))
879  & -c4*(x(2)*gr(1,1)+x(3)*gr(2,1)+x(5)*gr(3,1)+
880  & 2.d0*(x(8)*gr(4,1)+x(12)*gr(5,1)+x(17)*gr(6,1)))
881  al2(3)=al2(3)+au2(4)*r(8)+au2(5)*r(9)+au2(6)*r(10)+
882  & 2.d0*(au2(9)*r(11)+au2(13)*r(12)+au2(18)*r(13))
883  & -c4*(x(4)*gr(1,1)+x(5)*gr(2,1)+x(6)*gr(3,1)+
884  & 2.d0*(x(9)*gr(4,1)+x(13)*gr(5,1)+x(18)*gr(6,1)))
885  al2(4)=al2(4)+au2(7)*r(8)+au2(8)*r(9)+au2(9)*r(10)+
886  & 2.d0*(au2(10)*r(11)+au2(14)*r(12)+au2(19)*r(13))
887  & -c4*(x(7)*gr(1,1)+x(8)*gr(2,1)+x(9)*gr(3,1)+
888  & 2.d0*(x(10)*gr(4,1)+x(14)*gr(5,1)+x(19)*gr(6,1)))
889  al2(5)=al2(5)+au2(11)*r(8)+au2(12)*r(9)+au2(13)*r(10)+
890  & 2.d0*(au2(14)*r(11)+au2(15)*r(12)+au2(20)*r(13))
891  & -c4*(x(11)*gr(1,1)+x(12)*gr(2,1)+x(13)*gr(3,1)+
892  & 2.d0*(x(14)*gr(4,1)+x(15)*gr(5,1)+x(20)*gr(6,1)))
893  al2(6)=al2(6)+au2(16)*r(8)+au2(17)*r(9)+au2(18)*r(10)+
894  & 2.d0*(au2(19)*r(11)+au2(20)*r(12)+au2(21)*r(13))
895  & -c4*(x(16)*gr(1,1)+x(17)*gr(2,1)+x(18)*gr(3,1)+
896  & 2.d0*(x(19)*gr(4,1)+x(20)*gr(5,1)+x(21)*gr(6,1)))
897 !
898 ! update the consistency parameter
899 !
900  dg=dg+ddg
901 !
902 ! end of major loop
903 !
904  enddo
905 !
906 ! storing the stress
907 !
908  do i=1,6
909  stre(i)=stri(i)
910  enddo
911 !
912 ! calculating the tangent stiffness matrix
913 !
914  if(icmd.ne.3) then
915 !
916 ! determining p
917 !
918  gr(1,1)=au1(1)
919  gr(1,2)=au1(2)
920  gr(2,2)=au1(3)
921  gr(1,3)=au1(4)
922  gr(2,3)=au1(5)
923  gr(3,3)=au1(6)
924  gr(1,4)=au1(7)
925  gr(2,4)=au1(8)
926  gr(3,4)=au1(9)
927  gr(4,4)=au1(10)
928  gr(1,5)=au1(11)
929  gr(2,5)=au1(12)
930  gr(3,5)=au1(13)
931  gr(4,5)=au1(14)
932  gr(5,5)=au1(15)
933  gr(1,6)=au1(16)
934  gr(2,6)=au1(17)
935  gr(3,6)=au1(18)
936  gr(4,6)=au1(19)
937  gr(5,6)=au1(20)
938  gr(6,6)=au1(21)
939  do i=1,6
940  do j=1,i-1
941  gr(i,j)=gr(j,i)
942  enddo
943  enddo
944  nrhs=6
945 !
946  call dgetrs('No transpose',neq,nrhs,gl,lda,ipiv,gr,ldb,info)
947  if(info.ne.0) then
948  write(*,*) '*ERROR in sc.f: linear equation solver'
949  write(*,*) ' exited with error: info = ',info
950  call exit(201)
951  endif
952 !
953  stiff(1)=gr(1,1)-gm1*pn(1)*pn(1)
954  stiff(2)=gr(1,2)-gm1*pn(1)*pn(2)
955  stiff(3)=gr(2,2)-gm1*pn(2)*pn(2)
956  stiff(4)=gr(1,3)-gm1*pn(1)*pn(3)
957  stiff(5)=gr(2,3)-gm1*pn(2)*pn(3)
958  stiff(6)=gr(3,3)-gm1*pn(3)*pn(3)
959  stiff(7)=gr(1,4)-gm1*pn(1)*pn(4)
960  stiff(8)=gr(2,4)-gm1*pn(2)*pn(4)
961  stiff(9)=gr(3,4)-gm1*pn(3)*pn(4)
962  stiff(10)=gr(4,4)-gm1*pn(4)*pn(4)
963  stiff(11)=gr(1,5)-gm1*pn(1)*pn(5)
964  stiff(12)=gr(2,5)-gm1*pn(2)*pn(5)
965  stiff(13)=gr(3,5)-gm1*pn(3)*pn(5)
966  stiff(14)=gr(4,5)-gm1*pn(4)*pn(5)
967  stiff(15)=gr(5,5)-gm1*pn(5)*pn(5)
968  stiff(16)=gr(1,6)-gm1*pn(1)*pn(6)
969  stiff(17)=gr(2,6)-gm1*pn(2)*pn(6)
970  stiff(18)=gr(3,6)-gm1*pn(3)*pn(6)
971  stiff(19)=gr(4,6)-gm1*pn(4)*pn(6)
972  stiff(20)=gr(5,6)-gm1*pn(5)*pn(6)
973  stiff(21)=gr(6,6)-gm1*pn(6)*pn(6)
974 !
975  endif
976 !
977 ! updating the state variables
978 !
979  xstate(1,iint,iel)=eeq+c0*dg
980  do i=1,6
981  xstate(1+i,iint,iel)=ep(i)
982  enddo
983  xstate(8,iint,iel)=al1
984  do i=1,6
985  xstate(8+i,iint,iel)=al2(i)
986  enddo
987 !
988  return
c1
static double * c1
Definition: mafillvcompmain.c:30
dgetrs
subroutine dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:461
au1
static double * au1
Definition: mafillkmain.c:30
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
orthotropic
subroutine orthotropic(orthol, anisox)
Definition: orthotropic.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)