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

Go to the source code of this file.

Functions/Subroutines

subroutine umat_abaqusnl (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, istep, kinc, pnewdt, nmethod, iperturb)
 

Function/Subroutine Documentation

◆ umat_abaqusnl()

subroutine umat_abaqusnl ( 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  istep,
integer  kinc,
real*8  pnewdt,
integer  nmethod,
integer, dimension(*)  iperturb 
)
23 !
24 ! calculates stiffness and stresses for a nonlinear material
25 ! defined by an ABAQUS umat routine
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 !
112 ! This routine allows for the use of an ABAQUS umat user subroutine
113 ! in CalculiX.
114 !
115 ! Note that the following fields are not supported
116 ! so far: sse,spd,scd,rpl,ddsddt,drplde,drpldt,predef,
117 ! dpred,pnewdt,celent,layer,kspt
118 !
119 ! Furthermore, the following fields have a different meaning in
120 ! ABAQUS and CalculiX:
121 !
122 ! temp: in CalculiX: temperature at the end of the increment
123 ! in ABAQUS: temperature at the start of the increment
124 ! dtemp: in CalculiX: zero
125 ! in ABAQUS: temperature increment
126 !
127  implicit none
128 !
129  character*80 amat
130 !
131  integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),i,
132  & iorien,nmethod,iperturb(*),istep,
133  & ndi,nshr,ntens,nprops,layer,kspt,jstep(4),kinc,kal(2,6),
134  & kel(4,21),j1,j2,j3,j4,j5,j6,j7,j8,jj,n,ier,j,matz
135 !
136  real*8 elconloc(21),stiff(21),emec(6),emec0(6),beta(6),stre(6),
137  & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
138  & time,ttime,skl(3,3),xa(3,3),ya(3,3,3,3),xstate(nstate_,mi(1),*),
139  & xstateini(nstate_,mi(1),*),w(3),fv1(3),fv2(3),d(6),
140  & v1,v2,v3,c(6),r(3,3),r0(3,3),eln0(6),eln(6),e(3,3),tkl(3,3),
141  & u(6),c2(6),dd,um1(3,3),z(3,3),u0(3,3)
142 !
143  real*8 ddsdde(6,6),sse,spd,scd,rpl,ddsddt(6),drplde(6),
144  & drpldt,stran(6),dstran(6),abqtime(2),predef,temp,dtemp,
145  & dpred,drot(3,3),celent,pnewdt
146 !
147  kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/))
148 !
149  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,
150  & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,
151  & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,
152  & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/))
153 !
154  d=(/1.d0,1.d0,1.d0,0.d0,0.d0,0.d0/)
155 !
156 ! filling field jstep
157 !
158  jstep(1)=istep
159  jstep(2)=nmethod
160  jstep(3)=iperturb(2)
161  if(iperturb(1).eq.1) then
162  jstep(4)=1
163  else
164  jstep(4)=0
165  endif
166 !
167 ! calculating the logarithmic mechanical strain at the
168 ! start of the increment
169 !
170  e(1,1)=emec0(1)
171  e(2,2)=emec0(2)
172  e(3,3)=emec0(3)
173  e(1,2)=emec0(4)
174  e(1,3)=emec0(5)
175  e(2,3)=emec0(6)
176  e(2,1)=emec0(4)
177  e(3,1)=emec0(5)
178  e(3,2)=emec0(6)
179 !
180 ! calculating the eigenvalues and eigenvectors
181 !
182  n=3
183  matz=1
184 !
185  call rs(n,n,e,w,matz,z,fv1,fv2,ier)
186 !
187  if(ier.ne.0) then
188  write(*,*) '
189  & *ERROR calculating the eigenvalues/vectors in umat_abaqusnl'
190  call exit(201)
191  endif
192 !
193 ! calculating the principal stretches at the start of the increment
194 !
195  do i=1,3
196  w(i)=dsqrt(2.d0*w(i)+1.d0)
197  enddo
198 !
199 ! calculating the invariants at the start of the increment
200 !
201  v1=w(1)+w(2)+w(3)
202  v2=w(1)*w(2)+w(2)*w(3)+w(3)*w(1)
203  v3=w(1)*w(2)*w(3)
204 !
205 ! calculating the right Cauchy-Green tensor at the start of the
206 ! increment
207 !
208  do i=1,3
209  c(i)=2.d0*emec0(i)+1.d0
210  enddo
211  do i=4,6
212  c(i)=2.d0*emec0(i)
213  enddo
214 !
215 ! calculating the square of the right Cauchy-Green tensor at the
216 ! start of the increment
217 !
218  c2(1)=c(1)*c(1)+c(4)*c(4)+c(5)*c(5)
219  c2(2)=c(4)*c(4)+c(2)*c(2)+c(6)*c(6)
220  c2(3)=c(5)*c(5)+c(6)*c(6)+c(3)*c(3)
221  c2(4)=c(1)*c(4)+c(4)*c(2)+c(5)*c(6)
222  c2(5)=c(1)*c(5)+c(4)*c(6)+c(5)*c(3)
223  c2(6)=c(4)*c(5)+c(2)*c(6)+c(6)*c(3)
224 !
225 ! calculating the right stretch tensor at the start of the increment
226 ! (cf. Simo and Hughes, Computational Inelasticity)
227 !
228  dd=v1*v2-v3
229  do i=1,6
230  u(i)=(-c2(i)+(v1*v1-v2)*c(i)+v1*v3*d(i))/dd
231  enddo
232 !
233  u0(1,1)=u(1)
234  u0(2,2)=u(2)
235  u0(3,3)=u(3)
236  u0(1,2)=u(4)
237  u0(1,3)=u(5)
238  u0(2,3)=u(6)
239  u0(2,1)=u(4)
240  u0(3,1)=u(5)
241  u0(3,2)=u(6)
242 !
243 ! calculating the inverse of the right stretch tensor at the start
244 ! of the increment
245 !
246  um1(1,1)=(c(1)-v1*u(1)+v2)/v3
247  um1(2,2)=(c(2)-v1*u(2)+v2)/v3
248  um1(3,3)=(c(3)-v1*u(3)+v2)/v3
249  um1(1,2)=(c(4)-v1*u(4))/v3
250  um1(1,3)=(c(5)-v1*u(5))/v3
251  um1(2,3)=(c(6)-v1*u(6))/v3
252  um1(2,1)=um1(1,2)
253  um1(3,1)=um1(1,3)
254  um1(3,2)=um1(2,3)
255 !
256 ! calculation of the local rotation tensor at the start of the
257 ! increment
258 !
259  do i=1,3
260  do j=1,3
261  r0(i,j)=xokl(i,1)*um1(1,j)+xokl(i,2)*um1(2,j)+
262  & xokl(i,3)*um1(3,j)
263  enddo
264  enddo
265 !
266 ! calculating the logarithmic strain at the start of the increment
267 !
268  do i=1,3
269  w(i)=log(w(i))
270  enddo
271 !
272 ! logarithmic strain in global coordinates at the start of the
273 ! increment
274 !
275  eln0(1)=z(1,1)*z(1,1)*w(1)+z(1,2)*z(1,2)*w(2)+
276  & z(1,3)*z(1,3)*w(3)
277  eln0(2)=z(2,1)*z(2,1)*w(1)+z(2,2)*z(2,2)*w(2)+
278  & z(2,3)*z(2,3)*w(3)
279  eln0(3)=z(3,1)*z(3,1)*w(1)+z(3,2)*z(3,2)*w(2)+
280  & z(3,3)*z(3,3)*w(3)
281  eln0(4)=z(1,1)*z(2,1)*w(1)+z(1,2)*z(2,2)*w(2)+
282  & z(1,3)*z(2,3)*w(3)
283  eln0(5)=z(1,1)*z(3,1)*w(1)+z(1,2)*z(3,2)*w(2)+
284  & z(1,3)*z(3,3)*w(3)
285  eln0(6)=z(2,1)*z(3,1)*w(1)+z(2,2)*z(3,2)*w(2)+
286  & z(2,3)*z(3,3)*w(3)
287 !
288 ! calculating the logarithmic mechanical strain at the
289 ! end of the increment
290 !
291  e(1,1)=emec(1)
292  e(2,2)=emec(2)
293  e(3,3)=emec(3)
294  e(1,2)=emec(4)
295  e(1,3)=emec(5)
296  e(2,3)=emec(6)
297  e(2,1)=emec(4)
298  e(3,1)=emec(5)
299  e(3,2)=emec(6)
300 !
301 ! calculating the eigenvalues and eigenvectors
302 !
303  call rs(n,n,e,w,matz,z,fv1,fv2,ier)
304 !
305  if(ier.ne.0) then
306  write(*,*) '
307  & *ERROR calculating the eigenvalues/vectors in umat_abaqusnl'
308  call exit(201)
309  endif
310 !
311 ! calculating the principal stretches at the end of the increment
312 !
313  do i=1,3
314  w(i)=dsqrt(2.d0*w(i)+1.d0)
315  enddo
316 !
317 ! calculating the invariants at the end of the increment
318 !
319  v1=w(1)+w(2)+w(3)
320  v2=w(1)*w(2)+w(2)*w(3)+w(3)*w(1)
321  v3=w(1)*w(2)*w(3)
322 !
323 ! calculating the right Cauchy-Green tensor at the end of the
324 ! increment
325 !
326  do i=1,3
327  c(i)=2.d0*emec(i)+1.d0
328  enddo
329  do i=4,6
330  c(i)=2.d0*emec(i)
331  enddo
332 !
333 ! calculating the square of the right Cauchy-Green tensor at the
334 ! end of the increment
335 !
336  c2(1)=c(1)*c(1)+c(4)*c(4)+c(5)*c(5)
337  c2(2)=c(4)*c(4)+c(2)*c(2)+c(6)*c(6)
338  c2(3)=c(5)*c(5)+c(6)*c(6)+c(3)*c(3)
339  c2(4)=c(1)*c(4)+c(4)*c(2)+c(5)*c(6)
340  c2(5)=c(1)*c(5)+c(4)*c(6)+c(5)*c(3)
341  c2(6)=c(4)*c(5)+c(2)*c(6)+c(6)*c(3)
342 !
343 ! calculating the right stretch tensor at the end of the increment
344 ! (cf. Simo and Hughes, Computational Inelasticity)
345 !
346  dd=v1*v2-v3
347  do i=1,6
348  u(i)=(-c2(i)+(v1*v1-v2)*c(i)+v1*v3*d(i))/dd
349  enddo
350 !
351 ! calculating the inverse of the right stretch tensor at the end
352 ! of the increment
353 !
354  um1(1,1)=(c(1)-v1*u(1)+v2)/v3
355  um1(2,2)=(c(2)-v1*u(2)+v2)/v3
356  um1(3,3)=(c(3)-v1*u(3)+v2)/v3
357  um1(1,2)=(c(4)-v1*u(4))/v3
358  um1(1,3)=(c(5)-v1*u(5))/v3
359  um1(2,3)=(c(6)-v1*u(6))/v3
360  um1(2,1)=um1(1,2)
361  um1(3,1)=um1(1,3)
362  um1(3,2)=um1(2,3)
363 !
364 ! calculation of the local rotation tensor at the end of the
365 ! increment
366 !
367  do i=1,3
368  do j=1,3
369  r(i,j)=xkl(i,1)*um1(1,j)+xkl(i,2)*um1(2,j)+
370  & xkl(i,3)*um1(3,j)
371  enddo
372  enddo
373 !
374 ! calculating the logarithmic strain at the end of the increment
375 ! Elog=Z.ln(w).Z^T
376 !
377  do i=1,3
378  w(i)=log(w(i))
379  enddo
380 !
381 ! logarithmic strain in global coordinates at the end of the
382 ! increment
383 !
384  eln(1)=z(1,1)*z(1,1)*w(1)+z(1,2)*z(1,2)*w(2)+
385  & z(1,3)*z(1,3)*w(3)
386  eln(2)=z(2,1)*z(2,1)*w(1)+z(2,2)*z(2,2)*w(2)+
387  & z(2,3)*z(2,3)*w(3)
388  eln(3)=z(3,1)*z(3,1)*w(1)+z(3,2)*z(3,2)*w(2)+
389  & z(3,3)*z(3,3)*w(3)
390  eln(4)=z(1,1)*z(2,1)*w(1)+z(1,2)*z(2,2)*w(2)+
391  & z(1,3)*z(2,3)*w(3)
392  eln(5)=z(1,1)*z(3,1)*w(1)+z(1,2)*z(3,2)*w(2)+
393  & z(1,3)*z(3,3)*w(3)
394  eln(6)=z(2,1)*z(3,1)*w(1)+z(2,2)*z(3,2)*w(2)+
395  & z(2,3)*z(3,3)*w(3)
396 c write(*,*) 'iel', iel
397 c write(*,*) 'emec',(emec(i),i=1,6)
398 c write(*,*) 'eln',(eln(i),i=1,6)
399 c write(*,*) 'r0',((r0(i,j),j=1,3),i=1,3)
400 c write(*,*) 'r',((r(i,j),j=1,3),i=1,3)
401 !
402 ! calculating the incremental rotation tensor
403 ! drot=r.r0^T
404 !
405  do i=1,3
406  do j=1,3
407  drot(i,j)=r(i,1)*r0(j,1)+r(i,2)*r0(j,2)+r(i,3)*r0(j,3)
408  enddo
409  enddo
410 !
411  ntens=6
412 !
413  do i=1,nstate_
414  xstate(i,iint,iel)=xstateini(i,iint,iel)
415  enddo
416 !
417  abqtime(1)=time-dtime
418  abqtime(2)=ttime+time-dtime
419 !
420  temp=t1l
421  dtemp=0.d0
422 !
423  ndi=3
424  nshr=3
425  ntens=ndi+nshr
426 !
427  nprops=-kode-100
428 !
429 ! taking local material orientations into account
430 !
431  if(iorien.ne.0) then
432  call transformatrix(orab(1,iorien),pgauss,skl)
433 !
434 ! rotating the strain at the start of the increment
435 ! into the local system: Elog'=T^T.Elog.T
436 !
437  xa(1,1)=eln0(1)
438  xa(1,2)=eln0(4)
439  xa(1,3)=eln0(5)
440  xa(2,1)=eln0(4)
441  xa(2,2)=eln0(2)
442  xa(2,3)=eln0(6)
443  xa(3,1)=eln0(5)
444  xa(3,2)=eln0(6)
445  xa(3,3)=eln0(3)
446 !
447  do jj=1,6
448  stran(jj)=0.d0
449  j1=kal(1,jj)
450  j2=kal(2,jj)
451  do j3=1,3
452  do j4=1,3
453  stran(jj)=stran(jj)+
454  & xa(j3,j4)*skl(j3,j1)*skl(j4,j2)
455  enddo
456  enddo
457  enddo
458 !
459 ! rotating the strain at the end of the increment
460 ! into the local system
461 !
462  xa(1,1)=eln(1)
463  xa(1,2)=eln(4)
464  xa(1,3)=eln(5)
465  xa(2,1)=eln(4)
466  xa(2,2)=eln(2)
467  xa(2,3)=eln(6)
468  xa(3,1)=eln(5)
469  xa(3,2)=eln(6)
470  xa(3,3)=eln(3)
471 !
472  do jj=1,6
473  dstran(jj)=-stran(jj)
474  j1=kal(1,jj)
475  j2=kal(2,jj)
476  do j3=1,3
477  do j4=1,3
478  dstran(jj)=dstran(jj)+
479  & xa(j3,j4)*skl(j3,j1)*skl(j4,j2)
480  enddo
481  enddo
482  enddo
483  else
484  do jj=1,6
485  stran(jj)=eln0(jj)
486  dstran(jj)=eln(jj)-eln0(jj)
487  enddo
488  endif
489 !
490 ! rotating the stress into the local system
491 ! s'=J^(-1).U.S.U^T (no orientation card) or
492 ! s'=J^(-1).U.T^T.S.T.U^T (orientation card)
493 !
494  if(iorien.ne.0) then
495  do i=1,3
496  do j=1,3
497  tkl(i,j)=u0(i,1)*skl(j,1)+u0(i,2)*skl(j,2)+
498  & u0(i,3)*skl(j,3)
499  enddo
500  enddo
501  else
502  do i=1,3
503  do j=1,3
504  tkl(i,j)=u0(i,j)
505  enddo
506  enddo
507  endif
508 !
509  xa(1,1)=stre(1)
510  xa(1,2)=stre(4)
511  xa(1,3)=stre(5)
512  xa(2,1)=stre(4)
513  xa(2,2)=stre(2)
514  xa(2,3)=stre(6)
515  xa(3,1)=stre(5)
516  xa(3,2)=stre(6)
517  xa(3,3)=stre(3)
518 !
519  do jj=1,6
520  stre(jj)=0.d0
521  j1=kal(1,jj)
522  j2=kal(2,jj)
523  do j3=1,3
524  do j4=1,3
525  stre(jj)=stre(jj)+
526  & xa(j3,j4)*tkl(j1,j3)*tkl(j2,j4)
527  enddo
528  enddo
529  stre(jj)=stre(jj)/voj
530  enddo
531 !
532 ! changing physical strain into engineering strain
533 ! ABAQUS uses the engineering strain!
534 !
535  do i=4,6
536  stran(i)=2.d0*stran(i)
537  dstran(i)=2.d0*dstran(i)
538  enddo
539 !
540  if(amat(1:1).eq.'@') then
541 !
542  call call_external_umat(stre,xstate(1,iint,iel),ddsdde,
543  & sse,spd,scd,rpl,ddsddt,drplde,drpldt,stran,dstran,
544  & abqtime,dtime,temp,dtemp ,predef,dpred,amat,ndi,nshr,
545  & ntens,nstate_,elconloc,nprops,pgauss,drot,pnewdt,
546  & celent,xokl,xkl,iel,iint,layer,kspt,jstep,kinc)
547 !
548  else
549 !
550  call umat(stre,xstate(1,iint,iel),ddsdde,sse,spd,scd,rpl,ddsddt
551  & ,drplde,drpldt,stran,dstran,abqtime,dtime,temp,dtemp
552  & ,predef,dpred,amat,ndi,nshr,ntens,nstate_,elconloc,nprops
553  & ,pgauss,drot,pnewdt,celent,xokl,xkl,iel,iint,layer,kspt
554  & ,jstep,kinc)
555 !
556  endif
557 !
558 ! rotating the stress into the global system
559 ! S=J.U^(-1).s'.U^(-T) (no orientation card) or
560 ! S=J.T.U^(-1).s'.U^(-T).T^T (orientation card)
561 !
562  if(iorien.ne.0) then
563  do i=1,3
564  do j=1,3
565  tkl(i,j)=skl(i,1)*um1(1,j)+skl(i,2)*um1(2,j)+
566  & skl(i,3)*um1(3,j)
567  enddo
568  enddo
569  else
570  do i=1,3
571  do j=1,3
572  tkl(i,j)=um1(i,j)
573  enddo
574  enddo
575  endif
576 !
577  xa(1,1)=stre(1)
578  xa(1,2)=stre(4)
579  xa(1,3)=stre(5)
580  xa(2,1)=stre(4)
581  xa(2,2)=stre(2)
582  xa(2,3)=stre(6)
583  xa(3,1)=stre(5)
584  xa(3,2)=stre(6)
585  xa(3,3)=stre(3)
586 !
587  do jj=1,6
588  stre(jj)=0.d0
589  j1=kal(1,jj)
590  j2=kal(2,jj)
591  do j3=1,3
592  do j4=1,3
593  stre(jj)=stre(jj)+
594  & xa(j3,j4)*tkl(j1,j3)*tkl(j2,j4)
595  enddo
596  enddo
597  stre(jj)=stre(jj)*vj
598  enddo
599 !
600 ! calculate the stiffness matrix (the matrix is symmetrized)
601 !
602  if(icmd.ne.3) then
603  stiff(1)=ddsdde(1,1)
604  stiff(2)=(ddsdde(1,2)+ddsdde(2,1))/2.d0
605  stiff(3)=ddsdde(2,2)
606  stiff(4)=(ddsdde(1,3)+ddsdde(3,1))/2.d0
607  stiff(5)=(ddsdde(2,3)+ddsdde(3,2))/2.d0
608  stiff(6)=ddsdde(3,3)
609  stiff(7)=(ddsdde(1,4)+ddsdde(4,1))/2.d0
610  stiff(8)=(ddsdde(2,4)+ddsdde(4,2))/2.d0
611  stiff(9)=(ddsdde(3,4)+ddsdde(4,3))/2.d0
612  stiff(10)=ddsdde(4,4)
613  stiff(11)=(ddsdde(1,5)+ddsdde(5,1))/2.d0
614  stiff(12)=(ddsdde(2,5)+ddsdde(5,2))/2.d0
615  stiff(13)=(ddsdde(3,5)+ddsdde(5,3))/2.d0
616  stiff(14)=(ddsdde(4,5)+ddsdde(5,4))/2.d0
617  stiff(15)=ddsdde(5,5)
618  stiff(16)=(ddsdde(1,6)+ddsdde(6,1))/2.d0
619  stiff(17)=(ddsdde(2,6)+ddsdde(6,2))/2.d0
620  stiff(18)=(ddsdde(3,6)+ddsdde(6,3))/2.d0
621  stiff(19)=(ddsdde(4,6)+ddsdde(6,4))/2.d0
622  stiff(20)=(ddsdde(5,6)+ddsdde(6,5))/2.d0
623  stiff(21)=ddsdde(6,6)
624 c stiff(1)=ddsdde(1,1)
625 c stiff(2)=ddsdde(1,2)
626 c stiff(3)=ddsdde(2,2)
627 c stiff(4)=ddsdde(1,3)
628 c stiff(5)=ddsdde(2,3)
629 c stiff(6)=ddsdde(3,3)
630 c stiff(7)=ddsdde(1,4)
631 c stiff(8)=ddsdde(2,4)
632 c stiff(9)=ddsdde(3,4)
633 c stiff(10)=ddsdde(4,4)
634 c stiff(11)=ddsdde(1,5)
635 c stiff(12)=ddsdde(2,5)
636 c stiff(13)=ddsdde(3,5)
637 c stiff(14)=ddsdde(4,5)
638 c stiff(15)=ddsdde(5,5)
639 c stiff(16)=ddsdde(1,6)
640 c stiff(17)=ddsdde(2,6)
641 c stiff(18)=ddsdde(3,6)
642 c stiff(19)=ddsdde(4,6)
643 c stiff(20)=ddsdde(5,6)
644 c stiff(21)=ddsdde(6,6)
645 !
646 ! rotating the stiffness coefficients into the global system
647 !
648  call anisotropic(stiff,ya)
649 !
650  do jj=1,21
651  j1=kel(1,jj)
652  j2=kel(2,jj)
653  j3=kel(3,jj)
654  j4=kel(4,jj)
655  stiff(jj)=0.d0
656  do j5=1,3
657  do j6=1,3
658  do j7=1,3
659  do j8=1,3
660  stiff(jj)=stiff(jj)+ya(j5,j6,j7,j8)*
661  & tkl(j1,j5)*tkl(j2,j6)*tkl(j3,j7)*tkl(j4,j8)
662  enddo
663  enddo
664  enddo
665  enddo
666  enddo
667  endif
668 !
669  return
rs
subroutine rs(nm, n, a, w, matz, z, fv1, fv2, ierr)
Definition: rs.f:27
umat
subroutine umat(stress, statev, ddsdde, sse, spd, scd, rpl, ddsddt, drplde, drpldt, stran, dstran, time, dtime, temp, dtemp, predef, dpred, cmname, ndi, nshr, ntens, nstatv, props, nprops, coords, drot, pnewdt, celent, dfgrd0, dfgrd1, noel, npt, layer, kspt, kstep, kinc)
Definition: umat.f:24
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
v1
static double * v1
Definition: mafillsmmain_se.c:40
anisotropic
subroutine anisotropic(anisol, anisox)
Definition: anisotropic.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)