UMAT-ABAQUS
A general framework to develop material models (UMAT) in ABAQUS
umat_general.for
Go to the documentation of this file.
1  SUBROUTINE pull2(PK,SIG,FINV,DET,NDI)
2 C> PULL-BACK TIMES DET OF A 2ND ORDER TENSOR
3  IMPLICIT NONE
4  include 'PARAM_UMAT.INC'
5 C
6  INTEGER I1,J1,II1,JJ1,NDI
7  DOUBLE PRECISION PK(ndi,ndi),FINV(ndi,ndi)
8  DOUBLE PRECISION SIG(ndi,ndi)
9  DOUBLE PRECISION AUX,DET
10 C
11  DO i1=1,ndi
12  DO j1=1,ndi
13  aux=zero
14  DO ii1=1,ndi
15  DO jj1=1,ndi
16  aux=aux+det*finv(i1,ii1)*finv(j1,jj1)*sig(ii1,jj1)
17  END DO
18  END DO
19  pk(i1,j1)=aux
20  END DO
21  END DO
22 C
23  RETURN
24  END SUBROUTINE pull2
25  SUBROUTINE pinvariants(A,INV4,NDI,ST,LAMBDA,BARLAMBDA,DET)
26 C> AND 4TH PSEUDO-INVARIANTS OF A TENSOR
27  IMPLICIT NONE
28  include 'PARAM_UMAT.INC'
29 C
30  INTEGER NDI,I,J
31  DOUBLE PRECISION A(ndi,ndi),DET,SCALE1,ST(ndi,ndi),LAMBDA
32  DOUBLE PRECISION BARLAMBDA,INV4
33 C
34  inv4=zero
35  DO i=1,ndi
36  DO j=1, ndi
37  inv4=inv4+a(i,j)*st(i,j)
38  ENDDO
39  ENDDO
40 C STRETCH
41  scale1=det**(-one /three)
42  barlambda=dsqrt(inv4)
43  lambda=barlambda/scale1
44 C
45  RETURN
46  END SUBROUTINE pinvariants
47  SUBROUTINE metvol(CVOL,C,PV,PPV,DET,NDI)
48 C> VOLUMETRIC MATERIAL ELASTICITY TENSOR
49  IMPLICIT NONE
50  include 'PARAM_UMAT.INC'
51 C
52  INTEGER NDI,I1,J1,K1,L1
53  DOUBLE PRECISION C(ndi,ndi),CINV(ndi,ndi),
54  1 cvol(ndi,ndi,ndi,ndi)
55  DOUBLE PRECISION PV,PPV,DET
56 C
57  CALL matinv3d(c,cinv,ndi)
58 C
59  DO i1 = 1, ndi
60  DO j1 = 1, ndi
61  DO k1 = 1, ndi
62  DO l1 = 1, ndi
63  cvol(i1,j1,k1,l1)=
64  1 det*ppv*cinv(i1,j1)*cinv(k1,l1)
65  2 -det*pv*(cinv(i1,k1)*cinv(j1,l1)
66  3 +cinv(i1,l1)*cinv(j1,k1))
67  END DO
68  END DO
69  END DO
70  END DO
71 C
72  RETURN
73  END SUBROUTINE metvol
74  SUBROUTINE contraction42(S,LT,RT,NDI)
75 C> DOUBLE CONTRACTION BETWEEN 4TH ORDER AND 2ND ORDER TENSOR
76 C> INPUT:
77 C> LT - RIGHT 4TH ORDER TENSOR
78 C> RT - LEFT 2ND ODER TENSOR
79 C> OUTPUT:
80 C> S - DOUBLE CONTRACTED TENSOR (2ND ORDER)
81  IMPLICIT NONE
82  include 'PARAM_UMAT.INC'
83 C
84  INTEGER I1,J1,K1,L1,NDI
85 C
86  DOUBLE PRECISION RT(ndi,ndi),LT(ndi,ndi,ndi,ndi)
87  DOUBLE PRECISION S(ndi,ndi)
88  DOUBLE PRECISION AUX
89 C
90  DO i1=1,ndi
91  DO j1=1,ndi
92  aux=zero
93  DO k1=1,ndi
94  DO l1=1,ndi
95  aux=aux+lt(i1,j1,k1,l1)*rt(k1,l1)
96  END DO
97  END DO
98  s(i1,j1)=aux
99  END DO
100  END DO
101  RETURN
102  END SUBROUTINE contraction42
103 C>********************************************************************
104 C> Record of revisions: |
105 C> Date Programmer Description of change |
106 C> ==== ========== ===================== |
107 C> 15/11/2017 Joao Ferreira cont mech general eqs |
108 C> 01/11/2018 Joao Ferreira comments added |
109 C>--------------------------------------------------------------------
110 C> Description:
111 C> UMAT: IMPLEMENTATION OF THE CONSTITUTIVE EQUATIONS BASED UPON
112 C> A STRAIN-ENERGY FUNCTION (SEF).
113 C> THIS CODE, AS IS, EXPECTS A SEF BASED ON THE INVARIANTS OF THE
114 C> CAUCHY-GREEN TENSORS. A VISCOELASTIC COMPONENT IS ALSO
115 C> INCLUDED IF NEEDED.
116 C> YOU CAN CHOOSE TO COMPUTE AT THE MATERIAL FRAME AND THEN
117 C> PUSHFORWARD OR COPUTE AND THE SPATIAL FRAME DIRECTLY.
118 C>--------------------------------------------------------------------
119 C> IF YOU WANT TO ADAPT THE CODE ACCORDING TO YOUR SEF:
120 C> ISOMAT - DERIVATIVES OF THE SEF IN ORDER TO THE INVARIANTS
121 C> ADD OTHER CONTRIBUTIONS: STRESS AND TANGENT MATRIX
122 C>--------------------------------------------------------------------
123 C STATE VARIABLES: CHECK ROUTINES - INITIALIZE, WRITESDV, READSDV.
124 C>--------------------------------------------------------------------
125 C> UEXTERNALDB: READ FILAMENTS ORIENTATION AND PREFERED DIRECTION
126 C>--------------------------------------------------------------------
127 C>---------------------------------------------------------------------
128  SUBROUTINE umat(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
129  1 rpl,ddsddt,drplde,drpldt,
130  2 stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname,
131  3 ndi,nshr,ntens,nstatev,props,nprops,coords,drot,pnewdt,
132  4 celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)
133 C
134 C----------------------------------------------------------------------
135 C--------------------------- DECLARATIONS -----------------------------
136 C----------------------------------------------------------------------
137  IMPLICIT NONE
138  include 'PARAM_UMAT.INC'
139 C ADD COMMON BLOCKS HERE IF NEEDED (and in uexternal)
140 C COMMON /KBLOCK/KBLOCK
141 
142  COMMON /kfib/fibori
143 C
144  CHARACTER*8 CMNAME
145 C
146  INTEGER NDI, NSHR, NTENS, NSTATEV, NPROPS, NOEL, NPT,
147  1 layer, kspt, kstep, kinc
148 C
149  DOUBLE PRECISION STRESS(ntens),STATEV(nstatev),
150  1 ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),
151  2 stran(ntens),dstran(ntens),time(2),predef(1),dpred(1),
152  3 props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3),
153  4 fibori(nelem,4)
154 C
155  DOUBLE PRECISION SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP,
156  1 dtemp,pnewdt,celent
157 C
158  INTEGER NTERM
159 C
160 C FLAGS
161 C INTEGER FLAG1
162 C UTILITY TENSORS
163  DOUBLE PRECISION UNIT2(ndi,ndi),UNIT4(ndi,ndi,ndi,ndi),
164  1 unit4s(ndi,ndi,ndi,ndi),
165  2 proje(ndi,ndi,ndi,ndi),projl(ndi,ndi,ndi,ndi)
166 C KINEMATICS
167  DOUBLE PRECISION DISTGR(ndi,ndi),C(ndi,ndi),B(ndi,ndi),
168  1 cbar(ndi,ndi),bbar(ndi,ndi),distgrinv(ndi,ndi),
169  2 ubar(ndi,ndi),vbar(ndi,ndi),rot(ndi,ndi),
170  3 dfgrd1inv(ndi,ndi)
171  DOUBLE PRECISION DET,CBARI1,CBARI2
172 C VOLUMETRIC CONTRIBUTION
173  DOUBLE PRECISION PKVOL(ndi,ndi),SVOL(ndi,ndi),
174  1 cvol(ndi,ndi,ndi,ndi),cmvol(ndi,ndi,ndi,ndi)
175  DOUBLE PRECISION KBULK,PV,PPV,SSEV
176 C ISOCHORIC CONTRIBUTION
177  DOUBLE PRECISION SISO(ndi,ndi),PKISO(ndi,ndi),PK2(ndi,ndi),
178  1 ciso(ndi,ndi,ndi,ndi),cmiso(ndi,ndi,ndi,ndi),
179  2 sfic(ndi,ndi),cfic(ndi,ndi,ndi,ndi),
180  3 pkfic(ndi,ndi),cmfic(ndi,ndi,ndi,ndi)
181 C ISOCHORIC ISOTROPIC CONTRIBUTION
182  DOUBLE PRECISION C10,C01,SSEISO,DISO(5),PKMATFIC(ndi,ndi),
183  1 smatfic(ndi,ndi),sisomatfic(ndi,ndi),
184  2 cmisomatfic(ndi,ndi,ndi,ndi),
185  3 cisomatfic(ndi,ndi,ndi,ndi)
186 C ISOCHORIC ANISOTROPIC CONTRIBUTION
187  DOUBLE PRECISION K1,K2,KDISP,SSEANISO,
188  1 daniso(4),
189  2 pkmatficaniso(ndi,ndi),
190  3 sanisomatfic(ndi,ndi),
191  4 cmanisomatfic(ndi,ndi,ndi,ndi),
192  6 canisomatfic(ndi,ndi,ndi,ndi),
193  8 lambda,barlambda,
194  9 cbari4
195  DOUBLE PRECISION VORIF(3),VD(3),M0(3,3),MM(3,3),
196  1 vorif2(3),vd2(3),n0(3,3),nn(3,3)
197 C LIST VARS OF OTHER CONTRIBUTIONS HERE
198 C
199 C VISCOUS PROPERTIES (GENERALIZED MAXWEL DASHPOTS)
200  DOUBLE PRECISION VSCPROPS(6)
201  INTEGER VV
202 C JAUMMAN RATE CONTRIBUTION (REQUIRED FOR ABAQUS UMAT)
203  DOUBLE PRECISION CJR(ndi,ndi,ndi,ndi)
204 C CAUCHY STRESS AND ELASTICITY TENSOR
205  DOUBLE PRECISION SIGMA(ndi,ndi),DDSIGDDE(ndi,ndi,ndi,ndi),
206  1 ddpkdde(ndi,ndi,ndi,ndi)
207 C TESTING/DEBUG VARS
208  DOUBLE PRECISION STEST(ndi,ndi), CTEST(ndi,ndi,ndi,ndi)
209 C----------------------------------------------------------------------
210 C-------------------------- INITIALIZATIONS ---------------------------
211 C----------------------------------------------------------------------
212 C IDENTITY AND PROJECTION TENSORS
213  unit2=zero
214  unit4=zero
215  unit4s=zero
216  proje=zero
217  projl=zero
218 C KINEMATICS
219  distgr=zero
220  c=zero
221  b=zero
222  cbar=zero
223  bbar=zero
224  ubar=zero
225  vbar=zero
226  rot=zero
227  det=zero
228  cbari1=zero
229  cbari2=zero
230 C VOLUMETRIC
231  pkvol=zero
232  svol=zero
233  cvol=zero
234  kbulk=zero
235  pv=zero
236  ppv=zero
237  ssev=zero
238 C ISOCHORIC
239  siso=zero
240  pkiso=zero
241  pk2=zero
242  ciso=zero
243  cfic=zero
244  sfic=zero
245  pkfic=zero
246 C ISOTROPIC
247  c10=zero
248  c01=zero
249  sseiso=zero
250  diso=zero
251  pkmatfic=zero
252  smatfic=zero
253  sisomatfic=zero
254  cmisomatfic=zero
255  cisomatfic=zero
256 C INITIALIZE OTHER CONT HERE
257 C
258 C JAUMANN RATE
259  cjr=zero
260 C TOTAL CAUCHY STRESS AND ELASTICITY TENSORS
261  sigma=zero
262  ddsigdde=zero
263 C
264 C----------------------------------------------------------------------
265 C------------------------ IDENTITY TENSORS ----------------------------
266 C----------------------------------------------------------------------
267  CALL onem(unit2,unit4,unit4s,ndi)
268 C----------------------------------------------------------------------
269 C------------------- MATERIAL CONSTANTS AND DATA ----------------------
270 C----------------------------------------------------------------------
271 C VOLUMETRIC
272  kbulk = props(1)
273 C ISOCHORIC ISOTROPIC NEO HOOKE
274  c10 = props(2)
275  c01 = props(3)
276 C ISOCHORIC ANISOTROPIC GHO
277  k1 = props(4)
278  k2 = props(5)
279  kdisp = props(6)
280 C VISCOUS EFFECTS: MAXWELL ELEMENTS (MAX:3)
281 C VV = INT(PROPS(7))
282 C VSCPROPS = PROPS(8:13)
283 C NUMERICAL COMPUTATIONS
284  nterm = 60
285 C
286 C STATE VARIABLES
287 C
288  IF ((time(1).EQ.zero).AND.(kstep.EQ.1)) THEN
289  CALL initialize(statev)
290  ENDIF
291 C READ STATEV
292  CALL sdvread(statev)
293 C
294 C----------------------------------------------------------------------
295 C---------------------------- KINEMATICS ------------------------------
296 C----------------------------------------------------------------------
297 C DISTORTION GRADIENT
298  CALL fslip(dfgrd1,distgr,det,ndi)
299 C INVERSE OF DISTORTION GRADIENT
300  CALL matinv3d(dfgrd1,dfgrd1inv,ndi)
301 C INVERSE OF DISTORTION GRADIENT
302  CALL matinv3d(distgr,distgrinv,ndi)
303 C CAUCHY-GREEN DEFORMATION TENSORS
304  CALL deformation(dfgrd1,c,b,ndi)
305  CALL deformation(distgr,cbar,bbar,ndi)
306 C FIBER UNIT VECTOR AND STRUCTURAL TENSOR
307  CALL fibdir(fibori,m0,mm,nelem,noel,ndi,vorif,vd,distgr,dfgrd1)
308 C INVARIANTS OF DEVIATORIC DEFORMATION TENSORS
309  CALL invariants(cbar,cbari1,cbari2,ndi)
310 C
311  CALL pinvariants(cbar,cbari4,ndi,m0,lambda,barlambda,det)
312 C
313 C STRETCH TENSORS
314  CALL stretch(cbar,bbar,ubar,vbar,ndi)
315 C ROTATION TENSORS
316  CALL rotation(distgr,rot,ubar,ndi)
317 C DEVIATORIC PROJECTION TENSORS
318  CALL projeul(unit2,unit4s,proje,ndi)
319 C
320  CALL projlag(c,unit4,projl,ndi)
321 C----------------------------------------------------------------------
322 C--------------------- CONSTITUTIVE RELATIONS ------------------------
323 C----------------------------------------------------------------------
324 C
325 C---- VOLUMETRIC ------------------------------------------------------
326 C STRAIN-ENERGY AND DERIVATIVES (CHANGE HERE ACCORDING TO YOUR MODEL)
327  CALL vol(ssev,pv,ppv,kbulk,det)
328  CALL isomat(sseiso,diso,c10,c01,cbari1,cbari2)
329  CALL anisomat(sseaniso,daniso,diso,k1,k2,kdisp,cbari4,cbari1)
330 C
331 C---- ISOCHORIC ISOTROPIC ---------------------------------------------
332 C PK2 'FICTICIOUS' STRESS TENSOR
333  CALL pk2isomatfic(pkmatfic,diso,cbar,cbari1,unit2,ndi)
334 C CAUCHY 'FICTICIOUS' STRESS TENSOR
335  CALL sigisomatfic(sisomatfic,pkmatfic,distgr,det,ndi)
336 C 'FICTICIOUS' MATERIAL ELASTICITY TENSOR
337  CALL cmatisomatfic(cmisomatfic,cbar,cbari1,cbari2,
338  1 diso,unit2,unit4,det,ndi)
339 C 'FICTICIOUS' SPATIAL ELASTICITY TENSOR
340  CALL csisomatfic(cisomatfic,cmisomatfic,distgr,det,ndi)
341 C
342 C---- FIBERS (ONE FAMILY) -------------------------------------------
343 C
344  CALL pk2anisomatfic(pkmatficaniso,daniso,cbar,cbari4,m0,ndi)
345  CALL push2(sanisomatfic,pkmatficaniso,distgr,det,ndi)
346 C
347  CALL cmatanisomatfic(cmanisomatfic,m0,daniso,unit2,det,ndi)
348  CALL push4(canisomatfic,cmanisomatfic,distgr,det,ndi)
349 C----------------------------------------------------------------------
350 C SUM OF ALL ELASTIC CONTRIBUTIONS
351 C----------------------------------------------------------------------
352 C STRAIN-ENERGY
353  sse=ssev+sseiso+sseaniso
354 C PK2 'FICTICIOUS' STRESS
355  pkfic=pkmatfic+pkmatficaniso
356 C CAUCHY 'FICTICIOUS' STRESS
357  sfic=sisomatfic+sanisomatfic
358 C MATERIAL 'FICTICIOUS' ELASTICITY TENSOR
359  cmfic=cmisomatfic+cmanisomatfic
360 C SPATIAL 'FICTICIOUS' ELASTICITY TENSOR
361  cfic=cisomatfic+canisomatfic
362 C
363 C----------------------------------------------------------------------
364 C-------------------------- STRESS MEASURES ---------------------------
365 C----------------------------------------------------------------------
366 C
367 C---- VOLUMETRIC ------------------------------------------------------
368 C PK2 STRESS
369  CALL pk2vol(pkvol,pv,c,ndi)
370 C CAUCHY STRESS
371  CALL sigvol(svol,pv,unit2,ndi)
372 C
373 C---- ISOCHORIC -------------------------------------------------------
374 C PK2 STRESS
375  CALL pk2iso(pkiso,pkfic,projl,det,ndi)
376 C CAUCHY STRESS
377  CALL sigiso(siso,sfic,proje,ndi)
378 C
379 C---- VOLUMETRIC + ISOCHORIC ------------------------------------------
380 C PK2 STRESS
381  pk2 = pkvol + pkiso
382 C CAUCHY STRESS
383  sigma = svol + siso
384 C
385 C----------------------------------------------------------------------
386 C-------------------- MATERIAL ELASTICITY TENSOR ----------------------
387 C----------------------------------------------------------------------
388 C
389 C---- VOLUMETRIC ------------------------------------------------------
390 C
391  CALL metvol(cmvol,c,pv,ppv,det,ndi)
392 C
393 C---- ISOCHORIC -------------------------------------------------------
394 C
395  CALL metiso(cmiso,cmfic,projl,pkiso,pkfic,c,unit2,det,ndi)
396 C
397 C----------------------------------------------------------------------
398 C
399  ddpkdde= cmvol + cmiso
400 C
401 C----------------------------------------------------------------------
402 C--------------------- SPATIAL ELASTICITY TENSOR ----------------------
403 C----------------------------------------------------------------------
404 C
405 C---- VOLUMETRIC ------------------------------------------------------
406 C
407  CALL setvol(cvol,pv,ppv,unit2,unit4s,ndi)
408 C
409 C---- ISOCHORIC -------------------------------------------------------
410 C
411  CALL setiso(ciso,cfic,proje,siso,sfic,unit2,ndi)
412 C
413 C-----JAUMMAN RATE ----------------------------------------------------
414 C
415  CALL setjr(cjr,sigma,unit2,ndi)
416 C
417 C----------------------------------------------------------------------
418 C
419 C ELASTICITY TENSOR
420  ddsigdde=cvol+ciso+cjr
421 C
422 C----------------------------------------------------------------------
423 C------------------------- INDEX ALLOCATION ---------------------------
424 C----------------------------------------------------------------------
425 C VOIGT NOTATION - FULLY SIMMETRY IMPOSED
426  CALL indexx(stress,ddsdde,sigma,ddsigdde,ntens,ndi)
427 C
428 C----------------------------------------------------------------------
429 C--------------------------- STATE VARIABLES --------------------------
430 C----------------------------------------------------------------------
431 C DO K1 = 1, NTENS
432 C STATEV(1:27) = VISCOUS TENSORS
433  CALL sdvwrite(det,lambda,statev)
434 C END DO
435 C----------------------------------------------------------------------
436  RETURN
437  END
438 C----------------------------------------------------------------------
439 C--------------------------- END OF UMAT ------------------------------
440 C----------------------------------------------------------------------
441 C
442  SUBROUTINE projeul(A,AA,PE,NDI)
443 C> EULERIAN PROJECTION TENSOR
444 C INPUTS:
445 C IDENTITY TENSORS - A, AA
446 C OUTPUTS:
447 C 4TH ORDER SYMMETRIC EULERIAN PROJECTION TENSOR - PE
448 C
449  IMPLICIT NONE
450  include 'PARAM_UMAT.INC'
451 C
452  INTEGER I,J,K,L,NDI
453 C
454  DOUBLE PRECISION A(ndi,ndi),AA(ndi,ndi,ndi,ndi),
455  1 pe(ndi,ndi,ndi,ndi)
456 C
457  DO i=1,ndi
458  DO j=1,ndi
459  DO k=1,ndi
460  DO l=1,ndi
461  pe(i,j,k,l)=aa(i,j,k,l)-(one/three)*(a(i,j)*a(k,l))
462  END DO
463  END DO
464  END DO
465  END DO
466 C
467  RETURN
468  END SUBROUTINE projeul
469  SUBROUTINE sigvol(SVOL,PV,UNIT2,NDI)
470 C> VOLUMETRIC CAUCHY STRESS
471  IMPLICIT NONE
472  include 'PARAM_UMAT.INC'
473 C
474  INTEGER NDI,I1,J1
475  DOUBLE PRECISION UNIT2(ndi,ndi),SVOL(ndi,ndi)
476  DOUBLE PRECISION PV
477 C
478  DO i1=1,ndi
479  DO j1=1,ndi
480  svol(i1,j1)=pv*unit2(i1,j1)
481  END DO
482  END DO
483 C
484  RETURN
485  END SUBROUTINE sigvol
486  SUBROUTINE isomat(SSEISO,DISO,C10,C01,CBARI1,CBARI2)
487 C> ISOTROPIC MATRIX : ISOCHORIC SEF AND DERIVATIVES
488  IMPLICIT NONE
489  include 'PARAM_UMAT.INC'
490 C
491  DOUBLE PRECISION SSEISO,DISO(5)
492  DOUBLE PRECISION C10,C01,CBARI1,CBARI2
493 C
494  sseiso=c10*(cbari1-three)+c01*(cbari2-three)
495 C
496  !FIRST DERIVATIVE OF SSEISO IN ORDER TO I1
497  diso(1)=c10
498  !FIRST DERIVATIVE OF SSEISO IN ORDER TO I2
499  diso(2)=c01
500  !SECOND DERIVATIVE OF SSEISO IN ORDER TO I1
501  diso(3)=zero
502  !SECOND DERIVATIVE OF SSEISO IN ORDER TO I2
503  diso(4)=zero
504  !SECOND DERIVATIVE OF SSEISO IN ORDER TO I1 AND I2
505  diso(5)=zero
506 C
507  RETURN
508  END SUBROUTINE isomat
509  SUBROUTINE sdvread(STATEV)
510 C> VISCOUS DISSIPATION: READ STATE VARS
511  IMPLICIT NONE
512  include 'PARAM_UMAT.INC'
513 C
514  DOUBLE PRECISION STATEV(nsdv)
515 C read your sdvs here. they should be allocated.
516 C after the viscous terms (only if you use viscosity check hvread)
517 ! POS1=9*VV
518 ! DO I1=1,NCH
519 ! POS2=POS1+I1
520 ! FRAC(I1)=STATEV(POS2)
521 ! ENDDO
522 C
523 
524 C
525  RETURN
526 C
527  END SUBROUTINE sdvread
528  SUBROUTINE anisomat(SSEANISO,DANISO,DISO,K1,K2,KDISP,I4,I1)
529 C> ANISOTROPIC PART : ISOCHORIC SEF AND DERIVATIVES
530  IMPLICIT NONE
531  include 'PARAM_UMAT.INC'
532 C
533  DOUBLE PRECISION SSEISO,DANISO(4),DISO(5)
534  DOUBLE PRECISION K1,K2,KDISP,I4,I1
535  DOUBLE PRECISION DUDI1,D2UD2I1,SSEANISO
536  DOUBLE PRECISION E1,EE2,EE3,DUDI4,D2UD2I4,D2DUDI1DI4,D2DUDI2DI4
537 C
538  dudi1=diso(1)
539  d2ud2i1=diso(3)
540 C
541  e1=i4*(one-three*kdisp)+i1*kdisp-one
542 C
543  sseaniso=(k1/k2)*(dexp(k1*e1*e1)-one)
544 
545  IF(e1.GT.zero) THEN
546 C
547  ee2=dexp(k2*e1*e1)
548  ee3=(one+two*k2*e1*e1)
549 C
550  dudi1=dudi1+k1*kdisp*e1*ee2
551  d2ud2i1=d2ud2i1+k1*kdisp*kdisp*ee3*ee2
552 C
553  dudi4=k1*(one-three*kdisp)*e1*ee2
554 C
555  d2ud2i4=k1*((one-three*kdisp)**two)*ee3*ee2
556 
557  d2dudi1di4=k1*(one-three*kdisp)*kdisp*ee3*ee2
558  d2dudi2di4=zero
559 
560 
561 C
562  ELSE
563  dudi4=zero
564  d2ud2i4=zero
565  d2dudi1di4=zero
566  d2dudi2di4=zero
567 
568  d2ud2i1=zero
569 C
570  END IF
571  !FIRST DERIVATIVE OF SSEANISO IN ORDER TO I1
572  daniso(1)=dudi4
573  !FIRST DERIVATIVE OF SSEANISO IN ORDER TO I2
574  daniso(2)=d2ud2i4
575  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I1
576  daniso(3)=d2dudi1di4
577  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I2
578  daniso(4)=d2dudi2di4
579 C
580  diso(1)=dudi1
581  diso(3)=d2ud2i1
582 C
583  RETURN
584  END SUBROUTINE anisomat
585  SUBROUTINE pk2vol(PKVOL,PV,C,NDI)
586 C> VOLUMETRIC PK2 STRESS
587  IMPLICIT NONE
588  include 'PARAM_UMAT.INC'
589 C
590  INTEGER NDI,I1,J1
591  DOUBLE PRECISION PKVOL(ndi,ndi),C(ndi,ndi),CINV(ndi,ndi)
592  DOUBLE PRECISION PV
593 C
594  CALL matinv3d(c,cinv,ndi)
595 C
596  DO i1=1,ndi
597  DO j1=1,ndi
598  pkvol(i1,j1)=pv*cinv(i1,j1)
599  END DO
600  END DO
601 C
602  RETURN
603  END SUBROUTINE pk2vol
604  SUBROUTINE pull4(MAT,SPATIAL,FINV,DET,NDI)
605 C> PULL-BACK TIMES DET OF 4TH ORDER TENSOR
606  IMPLICIT NONE
607  include 'PARAM_UMAT.INC'
608 C
609  INTEGER I1,J1,K1,L1,II1,JJ1,KK1,LL1,NDI
610  DOUBLE PRECISION MAT(ndi,ndi,ndi,ndi),FINV(ndi,ndi)
611  DOUBLE PRECISION SPATIAL(ndi,ndi,ndi,ndi)
612  DOUBLE PRECISION AUX,DET
613 C
614  DO i1=1,ndi
615  DO j1=1,ndi
616  DO k1=1,ndi
617  DO l1=1,ndi
618  aux=zero
619  DO ii1=1,ndi
620  DO jj1=1,ndi
621  DO kk1=1,ndi
622  DO ll1=1,ndi
623  aux=aux+det*
624  + finv(i1,ii1)*finv(j1,jj1)*
625  + finv(k1,kk1)*finv(l1,ll1)*spatial(ii1,jj1,kk1,ll1)
626  END DO
627  END DO
628  END DO
629  END DO
630  mat(i1,j1,k1,l1)=aux
631  END DO
632  END DO
633  END DO
634  END DO
635 C
636  RETURN
637  END SUBROUTINE pull4
638  SUBROUTINE setjr(CJR,SIGMA,UNIT2,NDI)
639 C> JAUMAN RATE CONTRIBUTION FOR THE SPATIAL ELASTICITY TENSOR
640  IMPLICIT NONE
641  include 'PARAM_UMAT.INC'
642 C
643  INTEGER NDI,I1,J1,K1,L1
644  DOUBLE PRECISION UNIT2(ndi,ndi),
645  1 cjr(ndi,ndi,ndi,ndi),sigma(ndi,ndi)
646 C
647  DO i1 = 1, ndi
648  DO j1 = 1, ndi
649  DO k1 = 1, ndi
650  DO l1 = 1, ndi
651  cjr(i1,j1,k1,l1)=
652  1 (one/two)*(unit2(i1,k1)*sigma(j1,l1)
653  2 +sigma(i1,k1)*unit2(j1,l1)+unit2(i1,l1)*sigma(j1,k1)
654  3 +sigma(i1,l1)*unit2(j1,k1))
655  END DO
656  END DO
657  END DO
658  END DO
659 C
660  RETURN
661  END SUBROUTINE setjr
662  SUBROUTINE contraction44(S,LT,RT,NDI)
663 C> DOUBLE CONTRACTION BETWEEN 4TH ORDER TENSORS
664 C> INPUT:
665 C> LT - RIGHT 4TH ORDER TENSOR
666 C> RT - LEFT 4TH ORDER TENSOR
667 C> OUTPUT:
668 C> S - DOUBLE CONTRACTED TENSOR (4TH ORDER)
669  IMPLICIT NONE
670  include 'PARAM_UMAT.INC'
671 C
672  INTEGER I1,J1,K1,L1,M1,N1,NDI
673 C
674  DOUBLE PRECISION LT(ndi,ndi,ndi,ndi),RT(ndi,ndi,ndi,ndi)
675  DOUBLE PRECISION S(ndi,ndi,ndi,ndi)
676  DOUBLE PRECISION AUX
677 C
678  DO i1=1,ndi
679  DO j1=1,ndi
680  DO k1=1,ndi
681  DO l1=1,ndi
682  aux=zero
683  DO m1=1,ndi
684  DO n1=1,ndi
685  aux=aux+lt(i1,j1,m1,n1)*rt(m1,n1,k1,l1)
686  END DO
687  END DO
688  s(i1,j1,k1,l1)=aux
689  END DO
690  END DO
691  END DO
692  END DO
693 C
694  RETURN
695  END SUBROUTINE contraction44
696  SUBROUTINE initialize(STATEV)
697 C
698  IMPLICIT NONE
699  include 'PARAM_UMAT.INC'
700 C
701 C COMMON /KCOMMON/KBLOCK
702 C
703 C DOUBLE PRECISION TIME(2),KSTEP
704  INTEGER I1,POS,POS1,POS2,POS3
705  DOUBLE PRECISION STATEV(nsdv)
706 C DETERMINANT
707  statev(1)=one
708 C
709  RETURN
710 C
711  END SUBROUTINE initialize
712  SUBROUTINE setvol(CVOL,PV,PPV,UNIT2,UNIT4S,NDI)
713 C> VOLUMETRIC SPATIAL ELASTICITY TENSOR
714  IMPLICIT NONE
715  include 'PARAM_UMAT.INC'
716 C
717  INTEGER NDI,I1,J1,K1,L1
718  DOUBLE PRECISION UNIT2(ndi,ndi),UNIT4S(ndi,ndi,ndi,ndi),
719  1 cvol(ndi,ndi,ndi,ndi)
720  DOUBLE PRECISION PV,PPV
721 C
722  DO i1 = 1, ndi
723  DO j1 = 1, ndi
724  DO k1 = 1, ndi
725  DO l1 = 1, ndi
726  cvol(i1,j1,k1,l1)=
727  1 ppv*unit2(i1,j1)*unit2(k1,l1)
728  2 -two*pv*unit4s(i1,j1,k1,l1)
729  END DO
730  END DO
731  END DO
732  END DO
733 C
734  RETURN
735  END SUBROUTINE setvol
736  SUBROUTINE pk2isomatfic(FIC,DISO,CBAR,CBARI1,UNIT2,NDI)
737 C> ISOTROPIC MATRIX: 2PK 'FICTICIOUS' STRESS TENSOR
738 C INPUT:
739 C DISO - STRAIN-ENERGY DERIVATIVES
740 C CBAR - DEVIATORIC LEFT CAUCHY-GREEN TENSOR
741 C CBARI1,CBARI2 - CBAR INVARIANTS
742 C UNIT2 - 2ND ORDER IDENTITY TENSOR
743 C OUTPUT:
744 C FIC - 2ND PIOLA KIRCHOOF 'FICTICIOUS' STRESS TENSOR
745  IMPLICIT NONE
746  include 'PARAM_UMAT.INC'
747 C
748  INTEGER I1,J1,NDI
749  DOUBLE PRECISION FIC(ndi,ndi),DISO(5),CBAR(ndi,ndi),UNIT2(ndi,ndi)
750  DOUBLE PRECISION DUDI1,DUDI2,CBARI1
751  DOUBLE PRECISION AUX1,AUX2
752 C
753  dudi1=diso(1)
754  dudi2=diso(2)
755 C
756  aux1=two*(dudi1+cbari1*dudi2)
757  aux2=-two*dudi2
758 C
759  DO i1=1,ndi
760  DO j1=1,ndi
761  fic(i1,j1)=aux1*unit2(i1,j1)+aux2*cbar(i1,j1)
762  END DO
763  END DO
764 C
765  RETURN
766  END SUBROUTINE pk2isomatfic
767  SUBROUTINE uexternaldb(LOP,LRESTART,TIME,DTIME,KSTEP,KINC)
768 C> READ MESH DATA
769  include 'ABA_PARAM.INC'
770  include 'PARAM_UMAT.INC'
771 C
772 C UEXTERNAL just called once; work in parallel computing
773 C ADD COMMON BLOCKS HERE IF NEEDED (and in UMAT)
774 C COMMON /KBLOCK/KBLOCK
775  COMMON /kfib/fibori
776 C
777  REAL*8 DTIME
778  dimension time(2)
779  CHARACTER(256) FILENAME
780  CHARACTER(256) JOBDIR
781  INTEGER LENJOBDIR
782 
783  REAL*8 FIBORI(nelem,4)
784 
785 C LOP=0 --> START OF THE ANALYSIS
786  IF(lop.EQ.0.OR.lop.EQ.4) THEN
787 C
788  CALL getoutdir(jobdir,lenjobdir)
789 C DIR1 DEFNIED IN PARAM_UMAT.INC
790  filename=jobdir(:lenjobdir)//'/'//dir1
791 C
792  OPEN(15,file=filename)
793  DO i=1,nelem
794  READ(15,*) (fibori(i,j),j=1,4)
795  END DO
796  CLOSE(15)
797 !C
798  END IF
799 C
800  RETURN
801 C
802  END SUBROUTINE uexternaldb
803  SUBROUTINE getoutdir(OUTDIR, LENOUTDIR)
804 C> GET CURRENT WORKING DIRECTORY
805  include 'aba_param.inc'
806 C
807  CHARACTER*256 OUTDIR
808  INTEGER LENOUTDIR
809 C
810  CALL getcwd(outdir)
811 c OUTDIR=OUTDIR(1:SCAN(OUTDIR,'\',BACK=.TRUE.)-1)
812  lenoutdir=len_trim(outdir)
813 C
814  RETURN
815  END SUBROUTINE getoutdir
816  SUBROUTINE sigisomatfic(SFIC,PKFIC,F,DET,NDI)
817 C> ISOTROPIC MATRIX: ISOCHORIC CAUCHY STRESS
818  IMPLICIT NONE
819  include 'PARAM_UMAT.INC'
820 C
821  INTEGER NDI
822  DOUBLE PRECISION SFIC(ndi,ndi),F(ndi,ndi),
823  1 pkfic(ndi,ndi)
824  DOUBLE PRECISION DET
825 C
826  CALL push2(sfic,pkfic,f,det,ndi)
827 C
828  RETURN
829  END SUBROUTINE sigisomatfic
830  SUBROUTINE pk2anisomatfic(AFIC,DANISO,CBAR,INV4,ST0,NDI)
831 C> ANISOTROPIC MATRIX: 2PK 'FICTICIOUS' STRESS TENSOR
832 C INPUT:
833 C DANISO - ANISOTROPIC STRAIN-ENERGY DERIVATIVES
834 C CBAR - DEVIATORIC LEFT CAUCHY-GREEN TENSOR
835 C INV1,INV4 -CBAR INVARIANTS
836 C UNIT2 - 2ND ORDER IDENTITY TENSOR
837 C OUTPUT:
838 C AFIC - 2ND PIOLA KIRCHOOF 'FICTICIOUS' STRESS TENSOR
839  IMPLICIT NONE
840  include 'PARAM_UMAT.INC'
841 C
842  INTEGER NDI
843  DOUBLE PRECISION AFIC(ndi,ndi),DANISO(3),CBAR(3,3)
844  DOUBLE PRECISION DUDI4,DI4DC(3,3),INV4
845  DOUBLE PRECISION ST0(3,3)
846 C
847 C
848 C-----------------------------------------------------------------------------
849  !FIRST DERIVATIVE OF SSEANISO IN ORDER TO I4
850  dudi4=daniso(1)
851 C
852  di4dc=st0
853 C
854  afic=two*(dudi4*di4dc)
855 C
856  RETURN
857  END SUBROUTINE pk2anisomatfic
858  SUBROUTINE stretch(C,B,U,V,NDI)
859 C> STRETCH TENSORS
860  IMPLICIT NONE
861  include 'PARAM_UMAT.INC'
862 C
863  INTEGER NDI
864  DOUBLE PRECISION C(ndi,ndi),B(ndi,ndi),U(ndi,ndi),V(ndi,ndi)
865  DOUBLE PRECISION EIGVAL(ndi),OMEGA(ndi),EIGVEC(ndi,ndi)
866 C
867  CALL spectral(c,omega,eigvec)
868 C
869  eigval(1) = dsqrt(omega(1))
870  eigval(2) = dsqrt(omega(2))
871  eigval(3) = dsqrt(omega(3))
872 C
873  u(1,1) = eigval(1)
874  u(2,2) = eigval(2)
875  u(3,3) = eigval(3)
876 C
877  u = matmul(matmul(eigvec,u),transpose(eigvec))
878 C
879  CALL spectral(b,omega,eigvec)
880 C
881  eigval(1) = dsqrt(omega(1))
882  eigval(2) = dsqrt(omega(2))
883  eigval(3) = dsqrt(omega(3))
884 C write(*,*) eigvec(1,1),eigvec(2,1),eigvec(3,1)
885 C
886  v(1,1) = eigval(1)
887  v(2,2) = eigval(2)
888  v(3,3) = eigval(3)
889 C
890  v = matmul(matmul(eigvec,v),transpose(eigvec))
891  RETURN
892  END SUBROUTINE stretch
893  SUBROUTINE matinv3d(A,A_INV,NDI)
894 C> INVERSE OF A 3X3 MATRIX
895 C RETURN THE INVERSE OF A(3,3) - A_INV
896  IMPLICIT NONE
897 C
898  INTEGER NDI
899 C
900  DOUBLE PRECISION A(ndi,ndi),A_INV(ndi,ndi),DET_A,DET_A_INV
901 C
902  det_a = a(1,1)*(a(2,2)*a(3,3) - a(3,2)*a(2,3)) -
903  + a(2,1)*(a(1,2)*a(3,3) - a(3,2)*a(1,3)) +
904  + a(3,1)*(a(1,2)*a(2,3) - a(2,2)*a(1,3))
905 
906  IF (det_a .LE. 0.d0) THEN
907  WRITE(*,*) 'WARNING: SUBROUTINE MATINV3D:'
908  WRITE(*,*) 'WARNING: DET OF MAT=',det_a
909  RETURN
910  END IF
911 C
912  det_a_inv = 1.d0/det_a
913 C
914  a_inv(1,1) = det_a_inv*(a(2,2)*a(3,3)-a(3,2)*a(2,3))
915  a_inv(1,2) = det_a_inv*(a(3,2)*a(1,3)-a(1,2)*a(3,3))
916  a_inv(1,3) = det_a_inv*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
917  a_inv(2,1) = det_a_inv*(a(3,1)*a(2,3)-a(2,1)*a(3,3))
918  a_inv(2,2) = det_a_inv*(a(1,1)*a(3,3)-a(3,1)*a(1,3))
919  a_inv(2,3) = det_a_inv*(a(2,1)*a(1,3)-a(1,1)*a(2,3))
920  a_inv(3,1) = det_a_inv*(a(2,1)*a(3,2)-a(3,1)*a(2,2))
921  a_inv(3,2) = det_a_inv*(a(3,1)*a(1,2)-a(1,1)*a(3,2))
922  a_inv(3,3) = det_a_inv*(a(1,1)*a(2,2)-a(2,1)*a(1,2))
923 C
924  RETURN
925  END SUBROUTINE matinv3d
926  SUBROUTINE spectral(A,D,V)
927 C> EIGENVALUES AND EIGENVECTOR OF A 3X3 MATRIX
928 C THIS SUBROUTINE CALCULATES THE EIGENVALUES AND EIGENVECTORS OF
929 C A SYMMETRIC 3X3 MATRIX A.
930 C
931 C THE OUTPUT CONSISTS OF A VECTOR D CONTAINING THE THREE
932 C EIGENVALUES IN ASCENDING ORDER, AND A MATRIX V WHOSE
933 C COLUMNS CONTAIN THE CORRESPONDING EIGENVECTORS.
934 C
935  IMPLICIT NONE
936 C
937  INTEGER NP,NROT
938  parameter(np=3)
939 C
940  DOUBLE PRECISION D(3),V(3,3),A(3,3),E(3,3)
941 C
942  e = a
943 C
944  CALL jacobi(e,3,np,d,v,nrot)
945  CALL eigsrt(d,v,3,np)
946 C
947  RETURN
948  END SUBROUTINE spectral
949 
950 C***********************************************************************
951 
952  SUBROUTINE jacobi(A,N,NP,D,V,NROT)
953 C
954 C COMPUTES ALL EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
955 C MATRIX A, WHICH IS OF SIZE N BY N, STORED IN A PHYSICAL
956 C NP BY NP ARRAY. ON OUTPUT, ELEMENTS OF A ABOVE THE DIAGONAL
957 C ARE DESTROYED, BUT THE DIAGONAL AND SUB-DIAGONAL ARE UNCHANGED
958 C AND GIVE FULL INFORMATION ABOUT THE ORIGINAL SYMMETRIC MATRIX.
959 C VECTOR D RETURNS THE EIGENVALUES OF A IN ITS FIRST N ELEMENTS.
960 C V IS A MATRIX WITH THE SAME LOGICAL AND PHYSICAL DIMENSIONS AS
961 C A WHOSE COLUMNS CONTAIN, UPON OUTPUT, THE NORMALIZED
962 C EIGENVECTORS OF A. NROT RETURNS THE NUMBER OF JACOBI ROTATION
963 C WHICH WERE REQUIRED.
964 C
965 C THIS SUBROUTINE IS TAKEN FROM 'NUMERICAL RECIPES.'
966 C
967  IMPLICIT NONE
968  include 'PARAM_UMAT.INC'
969 C
970  INTEGER IP,IQ,N,NMAX,NP,NROT,I,J
971  parameter(nmax=100)
972 C
973  DOUBLE PRECISION A(np,np),D(np),V(np,np),B(nmax),Z(nmax),
974  + sm,tresh,g,t,h,theta,s,c,tau
975 
976 
977 C INITIALIZE V TO THE IDENTITY MATRIX
978  DO i=1,3
979  v(i,i)=one
980  DO j=1,3
981  IF (i.NE.j)THEN
982  v(i,j)=zero
983  ENDIF
984  END DO
985  END DO
986 C INITIALIZE B AND D TO THE DIAGONAL OF A, AND Z TO ZERO.
987 C THE VECTOR Z WILL ACCUMULATE TERMS OF THE FORM T*A_PQ AS
988 C IN EQUATION (11.1.14)
989 C
990  DO ip = 1,n
991  b(ip) = a(ip,ip)
992  d(ip) = b(ip)
993  z(ip) = 0.d0
994  END DO
995 
996 
997 C BEGIN ITERATION
998 C
999  nrot = 0
1000  DO i=1,50
1001 C
1002 C SUM OFF-DIAGONAL ELEMENTS
1003 C
1004  sm = 0.d0
1005  DO ip=1,n-1
1006  DO iq=ip+1,n
1007  sm = sm + dabs(a(ip,iq))
1008  END DO
1009  END DO
1010 C
1011 C IF SM = 0., THEN RETURN. THIS IS THE NORMAL RETURN,
1012 C WHICH RELIES ON QUADRATIC CONVERGENCE TO MACHINE
1013 C UNDERFLOW.
1014 C
1015  IF (sm.EQ.0.d0) RETURN
1016 C
1017 C IN THE FIRST THREE SWEEPS CARRY OUT THE PQ ROTATION ONLY IF
1018 C |A_PQ| > TRESH, WHERE TRESH IS SOME THRESHOLD VALUE,
1019 C SEE EQUATION (11.1.25). THEREAFTER TRESH = 0.
1020 C
1021  IF (i.LT.4) THEN
1022  tresh = 0.2d0*sm/n**2
1023  ELSE
1024  tresh = 0.d0
1025  END IF
1026 C
1027  DO ip=1,n-1
1028  DO iq=ip+1,n
1029  g = 100.d0*dabs(a(ip,iq))
1030 C
1031 C AFTER FOUR SWEEPS, SKIP THE ROTATION IF THE
1032 C OFF-DIAGONAL ELEMENT IS SMALL.
1033 C
1034  IF ((i.GT.4).AND.(dabs(d(ip))+g.EQ.dabs(d(ip)))
1035  + .AND.(dabs(d(iq))+g.EQ.dabs(d(iq)))) THEN
1036  a(ip,iq) = 0.d0
1037  ELSE IF (dabs(a(ip,iq)).GT.tresh) THEN
1038  h = d(iq) - d(ip)
1039  IF (dabs(h)+g.EQ.dabs(h)) THEN
1040 C
1041 C T = 1./(2.*THETA), EQUATION (11.1.10)
1042 C
1043  t =a(ip,iq)/h
1044  ELSE
1045  theta = 0.5d0*h/a(ip,iq)
1046  t =1.d0/(dabs(theta)+dsqrt(1.d0+theta**2.d0))
1047  IF (theta.LT.0.d0) t = -t
1048  END IF
1049  c = 1.d0/dsqrt(1.d0 + t**2.d0)
1050  s = t*c
1051  tau = s/(1.d0 + c)
1052  h = t*a(ip,iq)
1053  z(ip) = z(ip) - h
1054  z(iq) = z(iq) + h
1055  d(ip) = d(ip) - h
1056  d(iq) = d(iq) + h
1057  a(ip,iq) = 0.d0
1058 C
1059 C CASE OF ROTATIONS 1 <= J < P
1060 C
1061  DO j=1,ip-1
1062  g = a(j,ip)
1063  h = a(j,iq)
1064  a(j,ip) = g - s*(h + g*tau)
1065  a(j,iq) = h + s*(g - h*tau)
1066  END DO
1067 C
1068 C CASE OF ROTATIONS P < J < Q
1069 C
1070  DO j=ip+1,iq-1
1071  g = a(ip,j)
1072  h = a(j,iq)
1073  a(ip,j) = g - s*(h + g*tau)
1074  a(j,iq) = h + s*(g - h*tau)
1075  END DO
1076 C
1077 C CASE OF ROTATIONS Q < J <= N
1078 C
1079  DO j=iq+1,n
1080  g = a(ip,j)
1081  h = a(iq,j)
1082  a(ip,j) = g - s*(h + g*tau)
1083  a(iq,j) = h + s*(g - h*tau)
1084  END DO
1085  DO j = 1,n
1086  g = v(j,ip)
1087  h = v(j,iq)
1088  v(j,ip) = g - s*(h + g*tau)
1089  v(j,iq) = h + s*(g - h*tau)
1090  END DO
1091  nrot = nrot + 1
1092  END IF
1093  END DO
1094  END DO
1095 C
1096 C UPDATE D WITH THE SUM OF T*A_PQ, AND REINITIALIZE Z
1097 C
1098  DO ip=1,n
1099  b(ip) = b(ip) + z(ip)
1100  d(ip) = b(ip)
1101  z(ip) = 0.d0
1102  END DO
1103  END DO
1104 C
1105 C IF THE ALGORITHM HAS REACHED THIS STAGE, THEN THERE
1106 C ARE TOO MANY SWEEPS. PRINT A DIAGNOSTIC AND CUT THE
1107 C TIME INCREMENT.
1108 C
1109  WRITE (*,'(/1X,A/)') '50 ITERATIONS IN JACOBI SHOULD NEVER HAPPEN'
1110 C
1111  RETURN
1112  END SUBROUTINE jacobi
1113 
1114 C**********************************************************************
1115  SUBROUTINE eigsrt(D,V,N,NP)
1116 C
1117 C GIVEN THE EIGENVALUES D AND EIGENVECTORS V AS OUTPUT FROM
1118 C JACOBI, THIS SUBROUTINE SORTS THE EIGENVALUES INTO ASCENDING
1119 C ORDER AND REARRANGES THE COLMNS OF V ACCORDINGLY.
1120 C
1121 C THE SUBROUTINE WAS TAKEN FROM 'NUMERICAL RECIPES.'
1122 C
1123  IMPLICIT NONE
1124 C
1125  INTEGER N,NP,I,J,K
1126 C
1127  DOUBLE PRECISION D(np),V(np,np),P
1128 C
1129  DO i=1,n-1
1130  k = i
1131  p = d(i)
1132  DO j=i+1,n
1133  IF (d(j).GE.p) THEN
1134  k = j
1135  p = d(j)
1136  END IF
1137  END DO
1138  IF (k.NE.i) THEN
1139  d(k) = d(i)
1140  d(i) = p
1141  DO j=1,n
1142  p = v(j,i)
1143  v(j,i) = v(j,k)
1144  v(j,k) = p
1145  END DO
1146  END IF
1147  END DO
1148 C
1149  RETURN
1150  END SUBROUTINE eigsrt
1151  SUBROUTINE tensorprod2(A,B,C,NDI)
1152 C
1153  Implicit None
1154 C
1155  INTEGER I,J,K,L,NDI
1156 C
1157  DOUBLE PRECISION A(ndi,ndi),B(ndi,ndi),C(ndi,ndi,ndi,ndi)
1158 C
1159  DO i=1,ndi
1160  DO j=1,ndi
1161  DO k=1,ndi
1162  DO l=1,ndi
1163  c(i,j,k,l)=a(i,j)*b(k,l)
1164  END DO
1165  END DO
1166  END DO
1167  END DO
1168 C
1169  RETURN
1170 C
1171  end SUBROUTINE tensorprod2
1172  SUBROUTINE contraction24(S,LT,RT,NDI)
1173 C> DOUBLE CONTRACTION BETWEEN 4TH ORDER AND 2ND ORDER TENSOR
1174 C> INPUT:
1175 C> LT - RIGHT 2ND ORDER TENSOR
1176 C> RT - LEFT 4TH ODER TENSOR
1177 C> OUTPUT:
1178 C> S - DOUBLE CONTRACTED TENSOR (2ND ORDER)
1179  IMPLICIT NONE
1180  include 'PARAM_UMAT.INC'
1181 C
1182  INTEGER I1,J1,K1,L1,NDI
1183 C
1184  DOUBLE PRECISION LT(ndi,ndi),RT(ndi,ndi,ndi,ndi)
1185  DOUBLE PRECISION S(ndi,ndi)
1186  DOUBLE PRECISION AUX
1187 C
1188  DO k1=1,ndi
1189  DO l1=1,ndi
1190  aux=zero
1191  DO i1=1,ndi
1192  DO j1=1,ndi
1193  aux=aux+lt(k1,l1)*rt(i1,j1,k1,l1)
1194  END DO
1195  END DO
1196  s(k1,l1)=aux
1197  END DO
1198  END DO
1199  RETURN
1200  END SUBROUTINE contraction24
1201  SUBROUTINE invariants(A,INV1,INV2,NDI)
1202 C> 1ST AND 2ND INVARIANTS OF A TENSOR
1203  IMPLICIT NONE
1204  include 'PARAM_UMAT.INC'
1205 C
1206  INTEGER NDI,I1
1207  DOUBLE PRECISION A(ndi,ndi),AA(ndi,ndi)
1208  DOUBLE PRECISION INV1,INV1AA, INV2
1209 C
1210  inv1=zero
1211  inv1aa=zero
1212  aa=matmul(a,a)
1213  DO i1=1,ndi
1214  inv1=inv1+a(i1,i1)
1215  inv1aa=inv1aa+aa(i1,i1)
1216  END DO
1217  inv2=(one/two)*(inv1*inv1-inv1aa)
1218 C
1219  RETURN
1220  END SUBROUTINE invariants
1221  SUBROUTINE sigiso(SISO,SFIC,PE,NDI)
1222 C> ISOCHORIC CAUCHY STRESS
1223  IMPLICIT NONE
1224  include 'PARAM_UMAT.INC'
1225 C
1226  INTEGER NDI
1227  DOUBLE PRECISION SISO(ndi,ndi),
1228  1 pe(ndi,ndi,ndi,ndi),sfic(ndi,ndi)
1229 C
1230  CALL contraction42(siso,pe,sfic,ndi)
1231 C
1232  RETURN
1233  END SUBROUTINE sigiso
1234  SUBROUTINE vol(SSEV,PV,PPV,K,DET)
1235 C> VOLUMETRIC CONTRIBUTION :STRAIN ENERGY FUNCTION AND DERIVATIVES
1236  IMPLICIT NONE
1237  include 'PARAM_UMAT.INC'
1238 C
1239  DOUBLE PRECISION SSEV,PV,PPV
1240  DOUBLE PRECISION K,G,DET,AUX
1241 C
1242  g=(one/four)*(det*det-one-two*log(det))
1243 C
1244  ssev=k*g
1245 C
1246  pv=k*(one/two)*(det-one/det)
1247  aux=k*(one/two)*(one+one/(det*det))
1248  ppv=pv+det*aux
1249 C
1250  RETURN
1251  END SUBROUTINE vol
1252  SUBROUTINE push4(SPATIAL,MAT,F,DET,NDI)
1253 C> PIOLA TRANSFORMATION
1254 C> INPUT:
1255 C> MAT - MATERIAL ELASTICITY TENSOR
1256 C> F - DEFORMATION GRADIENT
1257 C> DET - DEFORMATION DETERMINANT
1258 C> OUTPUT:
1259 C> SPATIAL - SPATIAL ELASTICITY TENSOR
1260  IMPLICIT NONE
1261  include 'PARAM_UMAT.INC'
1262 C
1263  INTEGER I1,J1,K1,L1,II1,JJ1,KK1,LL1,NDI
1264 C
1265  DOUBLE PRECISION MAT(ndi,ndi,ndi,ndi),F(ndi,ndi)
1266  DOUBLE PRECISION SPATIAL(ndi,ndi,ndi,ndi)
1267  DOUBLE PRECISION AUX,DET
1268 C
1269  DO i1=1,ndi
1270  DO j1=1,ndi
1271  DO k1=1,ndi
1272  DO l1=1,ndi
1273  aux=zero
1274  DO ii1=1,ndi
1275  DO jj1=1,ndi
1276  DO kk1=1,ndi
1277  DO ll1=1,ndi
1278  aux=aux+(det**(-one))*
1279  + f(i1,ii1)*f(j1,jj1)*
1280  + f(k1,kk1)*f(l1,ll1)*mat(ii1,jj1,kk1,ll1)
1281  END DO
1282  END DO
1283  END DO
1284  END DO
1285  spatial(i1,j1,k1,l1)=aux
1286  END DO
1287  END DO
1288  END DO
1289  END DO
1290 C
1291  RETURN
1292  END SUBROUTINE push4
1293  SUBROUTINE pk2iso(PKISO,PKFIC,PL,DET,NDI)
1294 C> ISOCHORIC PK2 STRESS TENSOR
1295  IMPLICIT NONE
1296  include 'PARAM_UMAT.INC'
1297 C
1298  INTEGER NDI,I1,J1
1299  DOUBLE PRECISION PKISO(ndi,ndi),
1300  1 pl(ndi,ndi,ndi,ndi),pkfic(ndi,ndi)
1301  DOUBLE PRECISION DET,SCALE2
1302 C
1303  CALL contraction42(pkiso,pl,pkfic,ndi)
1304 C
1305  scale2=det**(-two/three)
1306  DO i1=1,ndi
1307  DO j1=1,ndi
1308  pkiso(i1,j1)=scale2*pkiso(i1,j1)
1309  END DO
1310  END DO
1311 C
1312  RETURN
1313  END SUBROUTINE pk2iso
1314  SUBROUTINE onem(A,AA,AAS,NDI)
1315 C
1316 C> THIS SUBROUTINE GIVES:
1317 C> 2ND ORDER IDENTITY TENSORS - A
1318 C> 4TH ORDER IDENTITY TENSOR - AA
1319 C> 4TH ORDER SYMMETRIC IDENTITY TENSOR - AAS
1320 C
1321  IMPLICIT NONE
1322  include 'PARAM_UMAT.INC'
1323 C
1324  INTEGER I,J,K,L,NDI
1325 C
1326  DOUBLE PRECISION A(ndi,ndi),AA(ndi,ndi,ndi,ndi),
1327  1 aas(ndi,ndi,ndi,ndi)
1328 C
1329  DO i=1,ndi
1330  DO j=1,ndi
1331  IF (i .EQ. j) THEN
1332  a(i,j) = one
1333  ELSE
1334  a(i,j) = zero
1335  END IF
1336  END DO
1337  END DO
1338 C
1339  DO i=1,ndi
1340  DO j=1,ndi
1341  DO k=1,ndi
1342  DO l=1,ndi
1343  aa(i,j,k,l)=a(i,k)*a(j,l)
1344  aas(i,j,k,l)=(one/two)*(a(i,k)*a(j,l)+a(i,l)*a(j,k))
1345  END DO
1346  END DO
1347  END DO
1348  END DO
1349 C
1350  RETURN
1351  END SUBROUTINE onem
1352  SUBROUTINE metiso(CMISO,CMFIC,PL,PKISO,PKFIC,C,UNIT2,DET,NDI)
1353 C> ISOCHORIC MATERIAL ELASTICITY TENSOR
1354  IMPLICIT NONE
1355  include 'PARAM_UMAT.INC'
1356 C
1357  INTEGER NDI,I1,J1,K1,L1
1358  DOUBLE PRECISION UNIT2(ndi,ndi),PL(ndi,ndi,ndi,ndi),
1359  1 cmiso(ndi,ndi,ndi,ndi),pkiso(ndi,ndi),
1360  2 cmfic(ndi,ndi,ndi,ndi),pkfic(ndi,ndi),
1361  3 cisoaux(ndi,ndi,ndi,ndi),
1362  4 cisoaux1(ndi,ndi,ndi,ndi),c(ndi,ndi),
1363  5 plt(ndi,ndi,ndi,ndi),cinv(ndi,ndi),
1364  6 pll(ndi,ndi,ndi,ndi)
1365  DOUBLE PRECISION TRFIC,XX,YY,ZZ,DET,AUX,AUX1
1366 C
1367  CALL matinv3d(c,cinv,ndi)
1368  cisoaux1=zero
1369  cisoaux=zero
1370  CALL contraction44(cisoaux1,pl,cmfic,ndi)
1371 C
1372 C transpose of lagrangian projection tensor
1373  DO i1=1,ndi
1374  DO j1=1,ndi
1375  DO k1=1,ndi
1376  DO l1=1,ndi
1377  plt(i1,j1,k1,l1)=pl(k1,l1,i1,j1)
1378  END DO
1379  END DO
1380  END DO
1381  END DO
1382 C
1383  CALL contraction44(cisoaux,cisoaux1,plt,ndi)
1384 C
1385  trfic=zero
1386  aux=det**(-two/three)
1387  aux1=aux**two
1388  DO i1=1,ndi
1389  trfic=trfic+aux*pkfic(i1,i1)*c(i1,i1)
1390  END DO
1391 C
1392  DO i1=1,ndi
1393  DO j1=1,ndi
1394  DO k1=1,ndi
1395  DO l1=1,ndi
1396  xx=aux1*cisoaux(i1,j1,k1,l1)
1397  pll(i1,j1,k1,l1)=(one/two)*(cinv(i1,k1)*cinv(j1,l1)+
1398  1 cinv(i1,l1)*cinv(j1,k1))-
1399  2 (one/three)*cinv(i1,j1)*cinv(k1,l1)
1400  yy=trfic*pll(i1,j1,k1,l1)
1401  zz=pkiso(i1,j1)*cinv(k1,l1)+cinv(i1,j1)*pkiso(k1,l1)
1402 C
1403  cmiso(i1,j1,k1,l1)=xx+(two/three)*yy-(two/three)*zz
1404  END DO
1405  END DO
1406  END DO
1407  END DO
1408 C
1409  RETURN
1410  END SUBROUTINE metiso
1411  SUBROUTINE setiso(CISO,CFIC,PE,SISO,SFIC,UNIT2,NDI)
1412 C> ISOCHORIC SPATIAL ELASTICITY TENSOR
1413  IMPLICIT NONE
1414  include 'PARAM_UMAT.INC'
1415 C
1416  INTEGER NDI,I1,J1,K1,L1
1417  DOUBLE PRECISION UNIT2(ndi,ndi),PE(ndi,ndi,ndi,ndi),
1418  1 ciso(ndi,ndi,ndi,ndi),siso(ndi,ndi),
1419  2 cfic(ndi,ndi,ndi,ndi),sfic(ndi,ndi),
1420  3 cisoaux(ndi,ndi,ndi,ndi),
1421  4 cisoaux1(ndi,ndi,ndi,ndi)
1422  DOUBLE PRECISION TRFIC,XX,YY,ZZ
1423 C
1424  cisoaux1=zero
1425  cisoaux=zero
1426 
1427  CALL contraction44(cisoaux1,pe,cfic,ndi)
1428  CALL contraction44(cisoaux,cisoaux1,pe,ndi)
1429 C
1430  trfic=zero
1431  DO i1=1,ndi
1432  trfic=trfic+sfic(i1,i1)
1433  END DO
1434 C
1435  DO i1=1,ndi
1436  DO j1=1,ndi
1437  DO k1=1,ndi
1438  DO l1=1,ndi
1439  xx=cisoaux(i1,j1,k1,l1)
1440  yy=trfic*pe(i1,j1,k1,l1)
1441  zz=siso(i1,j1)*unit2(k1,l1)+unit2(i1,j1)*siso(k1,l1)
1442 C
1443  ciso(i1,j1,k1,l1)=xx+(two/three)*yy-(two/three)*zz
1444  END DO
1445  END DO
1446  END DO
1447  END DO
1448 C
1449  RETURN
1450  END SUBROUTINE setiso
1451  SUBROUTINE push2(SIG,PK,F,DET,NDI)
1452 C> PIOLA TRANSFORMATION
1453 C> INPUT:
1454 C> PK - 2ND PIOLA KIRCHOOF STRESS TENSOR
1455 C> F - DEFORMATION GRADIENT
1456 C> DET - DEFORMATION DETERMINANT
1457 C> OUTPUT:
1458 C> SIG - CAUCHY STRESS TENSOR
1459  IMPLICIT NONE
1460  include 'PARAM_UMAT.INC'
1461 C
1462  INTEGER I1,J1,II1,JJ1,NDI
1463  DOUBLE PRECISION PK(ndi,ndi),F(ndi,ndi)
1464  DOUBLE PRECISION SIG(ndi,ndi)
1465  DOUBLE PRECISION AUX,DET
1466 C
1467  DO i1=1,ndi
1468  DO j1=1,ndi
1469  aux=zero
1470  DO ii1=1,ndi
1471  DO jj1=1,ndi
1472  aux=aux+(det**(-one))*f(i1,ii1)*f(j1,jj1)*pk(ii1,jj1)
1473  END DO
1474  END DO
1475  sig(i1,j1)=aux
1476  END DO
1477  END DO
1478 C
1479  RETURN
1480  END SUBROUTINE push2
1481  SUBROUTINE deformation(F,C,B,NDI)
1482 C> RIGHT AND LEFT CAUCHY-GREEN DEFORMATION TENSORS
1483  IMPLICIT NONE
1484  include 'PARAM_UMAT.INC'
1485 C
1486  INTEGER NDI
1487  DOUBLE PRECISION F(ndi,ndi),C(ndi,ndi),B(ndi,ndi)
1488 C RIGHT CAUCHY-GREEN DEFORMATION TENSOR
1489  c=matmul(transpose(f),f)
1490 C LEFT CAUCHY-GREEN DEFORMATION TENSOR
1491  b=matmul(f,transpose(f))
1492  RETURN
1493  END SUBROUTINE deformation
1494  SUBROUTINE fslip(F,FBAR,DET,NDI)
1495 C> DISTORTION GRADIENT
1496  IMPLICIT NONE
1497  include 'PARAM_UMAT.INC'
1498 C
1499  INTEGER NDI,I1,J1
1500  DOUBLE PRECISION F(ndi,ndi),FBAR(ndi,ndi)
1501  DOUBLE PRECISION DET,SCALE1
1502 C
1503 C JACOBIAN DETERMINANT
1504  det = f(1,1) * f(2,2) * f(3,3)
1505  1 - f(1,2) * f(2,1) * f(3,3)
1506 C
1507  IF (ndi .EQ. 3) THEN
1508  det = det + f(1,2) * f(2,3) * f(3,1)
1509  1 + f(1,3) * f(3,2) * f(2,1)
1510  2 - f(1,3) * f(3,1) * f(2,2)
1511  3 - f(2,3) * f(3,2) * f(1,1)
1512  END IF
1513 C
1514  scale1=det**(-one /three)
1515 C
1516  DO i1=1,ndi
1517  DO j1=1,ndi
1518  fbar(i1,j1)=scale1*f(i1,j1)
1519  END DO
1520  END DO
1521 C
1522  RETURN
1523  END SUBROUTINE fslip
1524  SUBROUTINE cmatanisomatfic(CMANISOMATFIC,M0,DANISO,UNIT2,DET,NDI)
1525 C
1526 C> ANISOTROPIC MATRIX: MATERIAL 'FICTICIOUS' ELASTICITY TENSOR
1527  IMPLICIT NONE
1528  include 'PARAM_UMAT.INC'
1529 C
1530  INTEGER NDI,I,J,K,L
1531  DOUBLE PRECISION CMANISOMATFIC(ndi,ndi,ndi,ndi),UNIT2(ndi,ndi),
1532  1 m0(ndi,ndi),daniso(3),det
1533  DOUBLE PRECISION CINV4(ndi,ndi,ndi,ndi),CINV14(ndi,ndi,ndi,ndi)
1534  DOUBLE PRECISION D2UDI4,D2UDI1DI4
1535  DOUBLE PRECISION IMM(ndi,ndi,ndi,ndi),MMI(ndi,ndi,ndi,ndi),
1536  1 mm0(ndi,ndi,ndi,ndi)
1537 C
1538 C-----------------------------------------------------------------------------
1539  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I4
1540  d2udi4=daniso(2)
1541  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I1 AND I4
1542  d2udi1di4=daniso(3)
1543 C
1544  CALL tensorprod2(m0,m0,mm0,ndi)
1545  CALL tensorprod2(unit2,m0,imm,ndi)
1546  CALL tensorprod2(m0,unit2,mmi,ndi)
1547 C
1548  DO i=1,ndi
1549  DO j=1,ndi
1550  DO k=1,ndi
1551  DO l=1,ndi
1552  cinv4(i,j,k,l)=d2udi4*mm0(i,j,k,l)
1553  cinv14(i,j,k,l)=d2udi1di4*(imm(i,j,k,l)+mmi(i,j,k,l))
1554  cmanisomatfic(i,j,k,l)=four*(cinv4(i,j,k,l)+cinv14(i,j,k,l))
1555  END DO
1556  END DO
1557  END DO
1558  END DO
1559 C
1560  RETURN
1561  END SUBROUTINE cmatanisomatfic
1562  SUBROUTINE sdvwrite(DET,LAMBDA,STATEV)
1563 C> VISCOUS DISSIPATION: WRITE STATE VARS
1564  IMPLICIT NONE
1565  include 'PARAM_UMAT.INC'
1566 C
1567  DOUBLE PRECISION STATEV(nsdv),DET,LAMBDA
1568 C write your sdvs here. they should be allocated
1569 C after the viscous terms (check hvwrite)
1570  statev(1)=det
1571  statev(2)=lambda
1572 
1573  RETURN
1574 C
1575  END SUBROUTINE sdvwrite
1576  SUBROUTINE resetdfgrd(DFGRD,NDI)
1577 C
1578  IMPLICIT NONE
1579  include 'PARAM_UMAT.INC'
1580 
1581  INTEGER NDI
1582  DOUBLE PRECISION DFGRD(ndi,ndi)
1583 
1584  dfgrd(1,1)= one
1585  dfgrd(1,2)= zero
1586  dfgrd(1,3)= zero
1587  dfgrd(2,1)= zero
1588  dfgrd(2,2)= one
1589  dfgrd(2,3)= zero
1590  dfgrd(3,1)= zero
1591  dfgrd(3,2)= zero
1592  dfgrd(3,3)= one
1593 
1594  END SUBROUTINE resetdfgrd
1595  SUBROUTINE fibdir(FIB,ST0,ST,NE,NOEL,NDI,VORIF,VD,DISTGR,DFGRD1)
1596 C
1597  IMPLICIT NONE
1598  include 'PARAM_UMAT.INC'
1599 C
1600  INTEGER NDI, NE, NOEL,INOEL,I,J,I1,J1
1601  DOUBLE PRECISION SUM1, DFGRD1(3,3), DNORM
1602  DOUBLE PRECISION VORIF(3),ST(3,3),VD(3),ST0(3,3),DISTGR(3,3)
1603  DOUBLE PRECISION FIB(ne,4)
1604 C
1605  inoel=0
1606  i=0
1607  DO i=1,ne
1608 C ELEMENT IDENTIFICATION
1609  IF(noel.EQ.int(fib(i,1))) THEN
1610  inoel=i
1611  ENDIF
1612  ENDDO
1613 C
1614 C FIB - FIBER ORIENTATION
1615  dnorm=dsqrt(fib(inoel,2)*fib(inoel,2)+
1616  1 fib(inoel,3)*fib(inoel,3)+
1617  2 fib(inoel,4)*fib(inoel,4))
1618 C
1619 C UNDERFORMED FIBER ORIENTATION TENSOR
1620 C
1621  DO i1=1,ndi
1622  j1=i1+1
1623 C FIBER ORIENTATION NORMALIZED VECTOR - FAMILY 1
1624  vorif(i1)=fib(inoel,j1)/dnorm
1625  END DO
1626 C
1627 
1628  DO i=1,ndi
1629  sum1=zero
1630  DO j=1,ndi
1631  sum1=sum1+dfgrd1(i,j)*vorif(j)
1632  ENDDO
1633 C FIBER DIRECTIONS IN THE DEFORMED CONFIGURATION
1634 C -FAMILY 1
1635  vd(i)=sum1
1636  ENDDO
1637  dnorm=dsqrt(vd(1)*vd(1)+
1638  1 vd(2)*vd(2)+
1639  2 vd(3)*vd(3))
1640 C COSINE OF THE ANGLE BETWEEN FIBERS
1641 C
1642 C
1643 C--------------------------------------------------------------------------
1644  DO i=1,ndi
1645  DO j=1,ndi
1646 C STRUCTURAL TENSOR - FAMILY 1
1647  st0(i,j)=vorif(i)*vorif(j)
1648  END DO
1649  END DO
1650 C
1651 C STRUCTURE TENSOR IN THE DEFORMED CONFIGURATION - FAMILY 1
1652  st=matmul(st0,transpose(distgr))
1653  st=matmul(distgr,st)
1654 C
1655 C
1656  RETURN
1657  END SUBROUTINE fibdir
1658  SUBROUTINE indexx(STRESS,DDSDDE,SIG,TNG,NTENS,NDI)
1659 C> INDEXATION: FULL SIMMETRY IN STRESSES AND ELASTICITY TENSORS
1660  IMPLICIT NONE
1661  include 'PARAM_UMAT.INC'
1662 C
1663  INTEGER II1(6),II2(6),NTENS,NDI,I1,J1
1664  DOUBLE PRECISION STRESS(ntens),DDSDDE(ntens,ntens)
1665  DOUBLE PRECISION SIG(ndi,ndi),TNG(ndi,ndi,ndi,ndi)
1666  DOUBLE PRECISION PP1,PP2
1667 C
1668  ii1(1)=1
1669  ii1(2)=2
1670  ii1(3)=3
1671  ii1(4)=1
1672  ii1(5)=1
1673  ii1(6)=2
1674 C
1675  ii2(1)=1
1676  ii2(2)=2
1677  ii2(3)=3
1678  ii2(4)=2
1679  ii2(5)=3
1680  ii2(6)=3
1681 C
1682  DO i1=1,ntens
1683 C STRESS VECTOR
1684  stress(i1)=sig(ii1(i1),ii2(i1))
1685  DO j1=1,ntens
1686 C DDSDDE - FULLY SIMMETRY IMPOSED
1687  pp1=tng(ii1(i1),ii2(i1),ii1(j1),ii2(j1))
1688  pp2=tng(ii1(i1),ii2(i1),ii2(j1),ii1(j1))
1689  ddsdde(i1,j1)=(one/two)*(pp1+pp2)
1690  END DO
1691  END DO
1692 C
1693  RETURN
1694 C
1695  END SUBROUTINE indexx
1696  SUBROUTINE rotation(F,R,U,NDI)
1697 C> COMPUTES ROTATION TENSOR
1698  IMPLICIT NONE
1699  include 'PARAM_UMAT.INC'
1700 C
1701  INTEGER NDI
1702  DOUBLE PRECISION F(ndi,ndi),R(ndi,ndi),U(ndi,ndi),UINV(ndi,ndi)
1703 C
1704  CALL matinv3d(u,uinv,ndi)
1705 C
1706  r = matmul(f,uinv)
1707  RETURN
1708  END SUBROUTINE rotation
1709  SUBROUTINE cmatisomatfic(CMISOMATFIC,CBAR,CBARI1,CBARI2,
1710  1 diso,unit2,unit4,det,ndi)
1711 C> ISOTROPIC MATRIX: MATERIAL 'FICTICIOUS' ELASTICITY TENSOR
1712  IMPLICIT NONE
1713  include 'PARAM_UMAT.INC'
1714 C
1715  INTEGER NDI,I1,J1,K1,L1
1716  DOUBLE PRECISION CMISOMATFIC(ndi,ndi,ndi,ndi),UNIT2(ndi,ndi),
1717  1 cbar(ndi,ndi),diso(5),
1718  2 unit4(ndi,ndi,ndi,ndi)
1719  DOUBLE PRECISION CBARI1,CBARI2
1720  DOUBLE PRECISION DUDI1,DUDI2,D2UD2I1,D2UD2I2,D2UDI1I2
1721  DOUBLE PRECISION AUX,AUX1,AUX2,AUX3,AUX4,DET
1722  DOUBLE PRECISION UIJ,UKL,CIJ,CKL
1723 C
1724  dudi1=diso(1)
1725  dudi2=diso(2)
1726  d2ud2i1=diso(3)
1727  d2ud2i2=diso(4)
1728  d2udi1i2=diso(5)
1729 C
1730  aux1=four*(d2ud2i1+two*cbari1*d2udi1i2+
1731  1 dudi2+cbari1*cbari1*d2ud2i2)
1732  aux2=-four*(d2udi1i2+cbari1*d2ud2i2)
1733  aux3=four*d2ud2i2
1734  aux4=-four*dudi2
1735 
1736  DO i1=1,ndi
1737  DO j1=1,ndi
1738  DO k1=1,ndi
1739  DO l1=1,ndi
1740  uij=unit2(i1,j1)
1741  ukl=unit2(k1,l1)
1742  cij=cbar(i1,j1)
1743  ckl=cbar(k1,l1)
1744  aux=aux1*uij*ukl+
1745  1 aux2*(uij*ckl+cij*ukl)+aux3*cij*ckl+
1746  3 aux4*unit4(i1,j1,k1,l1)
1747  cmisomatfic(i1,j1,k1,l1)=aux
1748  END DO
1749  END DO
1750  END DO
1751  END DO
1752 C
1753  RETURN
1754  END SUBROUTINE cmatisomatfic
1755  SUBROUTINE projlag(C,AA,PL,NDI)
1756 C> LAGRANGIAN PROJECTION TENSOR
1757 C INPUTS:
1758 C IDENTITY TENSORS - A, AA
1759 C ISOCHORIC LEFT CAUCHY GREEN TENSOR - C
1760 C INVERSE OF C - CINV
1761 C OUTPUTS:
1762 C 4TH ORDER SYMMETRIC LAGRANGIAN PROJECTION TENSOR - PL
1763 C
1764  IMPLICIT NONE
1765  include 'PARAM_UMAT.INC'
1766 C
1767  INTEGER I,J,K,L,NDI
1768 C
1769  DOUBLE PRECISION CINV(ndi,ndi),AA(ndi,ndi,ndi,ndi),
1770  1 pl(ndi,ndi,ndi,ndi),c(ndi,ndi)
1771 C
1772  CALL matinv3d(c,cinv,ndi)
1773 C
1774  DO i=1,ndi
1775  DO j=1,ndi
1776  DO k=1,ndi
1777  DO l=1,ndi
1778  pl(i,j,k,l)=aa(i,j,k,l)-(one/three)*(cinv(i,j)*c(k,l))
1779  END DO
1780  END DO
1781  END DO
1782  END DO
1783 C
1784  RETURN
1785  END SUBROUTINE projlag
1786  SUBROUTINE csisomatfic(CISOMATFIC,CMISOMATFIC,DISTGR,DET,NDI)
1787 C> ISOTROPIC MATRIX: SPATIAL 'FICTICIOUS' ELASTICITY TENSOR
1788  IMPLICIT NONE
1789  include 'PARAM_UMAT.INC'
1790 C
1791  INTEGER NDI
1792  DOUBLE PRECISION CMISOMATFIC(ndi,ndi),DISTGR(ndi,ndi),
1793  1 cisomatfic(ndi,ndi,ndi,ndi)
1794  DOUBLE PRECISION DET
1795 C
1796  CALL push4(cisomatfic,cmisomatfic,distgr,det,ndi)
1797 C
1798  RETURN
1799  END SUBROUTINE csisomatfic
pull4
subroutine pull4(MAT, SPATIAL, FINV, DET, NDI)
Definition: pull4.for:2
uexternaldb
subroutine uexternaldb(LOP, LRESTART, TIME, DTIME, KSTEP, KINC)
Definition: uexternaldb.for:2
sdvread
subroutine sdvread(STATEV)
Definition: sdvread.for:2
anisomat
subroutine anisomat(SSEANISO, DANISO, DISO, K1, K2, KDISP, I4, I1)
Definition: anisomat.for:2
initialize
subroutine initialize(STATEV)
Definition: initialize.for:2
contraction44
subroutine contraction44(S, LT, RT, NDI)
Definition: contraction44.for:2
deformation
subroutine deformation(F, C, B, NDI)
Definition: deformation.for:2
push2
subroutine push2(SIG, PK, F, DET, NDI)
Definition: push2.for:2
sdvwrite
subroutine sdvwrite(DET, LAMBDA, STATEV)
Definition: sdvwrite.for:2
setiso
subroutine setiso(CISO, CFIC, PE, SISO, SFIC, UNIT2, NDI)
Definition: setiso.for:2
cmatanisomatfic
subroutine cmatanisomatfic(CMANISOMATFIC, M0, DANISO, UNIT2, DET, NDI)
Definition: cmatanisomatfic.for:2
jacobi
subroutine jacobi(A, N, NP, D, V, NROT)
Definition: spectral.for:28
stretch
subroutine stretch(C, B, U, V, NDI)
Definition: stretch.for:2
matinv3d
subroutine matinv3d(A, A_INV, NDI)
Definition: minverse3d.for:2
onem
subroutine onem(A, AA, AAS, NDI)
Definition: identity_tensors.for:2
fibdir
subroutine fibdir(FIB, ST0, ST, NE, NOEL, NDI, VORIF, VD, DISTGR, DFGRD1)
Definition: fibdir.for:2
push4
subroutine push4(SPATIAL, MAT, F, DET, NDI)
Definition: push4.for:2
metvol
subroutine metvol(CVOL, C, PV, PPV, DET, NDI)
Definition: metvol.for:2
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, NSTATEV, PROPS, NPROPS, COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, KSTEP, KINC)
Record of revisions: | Date Programmer Description of change | ==== ========== =====================...
Definition: _umat.for:31
csisomatfic
subroutine csisomatfic(CISOMATFIC, CMISOMATFIC, DISTGR, DET, NDI)
Definition: cisomatfic.for:2
sigisomatfic
subroutine sigisomatfic(SFIC, PKFIC, F, DET, NDI)
Definition: sigisomatfic.for:2
invariants
subroutine invariants(A, INV1, INV2, NDI)
Definition: invariants.for:2
sigvol
subroutine sigvol(SVOL, PV, UNIT2, NDI)
Definition: sigvol.for:2
projeul
subroutine projeul(A, AA, PE, NDI)
Definition: proj_eulerian.for:2
getoutdir
subroutine getoutdir(OUTDIR, LENOUTDIR)
Definition: GETOUTDIR.for:2
contraction42
subroutine contraction42(S, LT, RT, NDI)
Definition: contraction42.for:2
indexx
subroutine indexx(STRESS, DDSDDE, SIG, TNG, NTENS, NDI)
Definition: index.for:2
pk2iso
subroutine pk2iso(PKISO, PKFIC, PL, DET, NDI)
Definition: pk2iso.for:2
fslip
subroutine fslip(F, FBAR, DET, NDI)
Definition: fslip.for:2
setjr
subroutine setjr(CJR, SIGMA, UNIT2, NDI)
Definition: setjr.for:2
pull2
subroutine pull2(PK, SIG, FINV, DET, NDI)
Definition: pull2.for:2
pinvariants
subroutine pinvariants(A, INV4, NDI, ST, LAMBDA, BARLAMBDA, DET)
Definition: pinvariants.for:2
cmatisomatfic
subroutine cmatisomatfic(CMISOMATFIC, CBAR, CBARI1, CBARI2, DISO, UNIT2, UNIT4, DET, NDI)
Definition: cmatisomatfic.for:3
isomat
subroutine isomat(SSEISO, DISO, C10, C01, CBARI1, CBARI2)
Definition: isomat.for:2
spectral
subroutine spectral(A, D, V)
Definition: spectral.for:2
tensorprod2
subroutine tensorprod2(A, B, C, NDI)
Definition: tensorprod22.for:2
eigsrt
subroutine eigsrt(D, V, N, NP)
Definition: spectral.for:191
pk2anisomatfic
subroutine pk2anisomatfic(AFIC, DANISO, CBAR, INV4, ST0, NDI)
Definition: pk2anisomatfic.for:2
contraction24
subroutine contraction24(S, LT, RT, NDI)
Definition: contraction24.for:2
projlag
subroutine projlag(C, AA, PL, NDI)
Definition: proj_lagrangian.for:2
rotation
subroutine rotation(F, R, U, NDI)
Definition: rotation.for:2
metiso
subroutine metiso(CMISO, CMFIC, PL, PKISO, PKFIC, C, UNIT2, DET, NDI)
Definition: metiso.for:2
resetdfgrd
subroutine resetdfgrd(DFGRD, NDI)
Definition: resetdfgr.for:2
vol
subroutine vol(SSEV, PV, PPV, K, DET)
Definition: vol.for:2
pk2vol
subroutine pk2vol(PKVOL, PV, C, NDI)
Definition: pk2vol.for:2
pk2isomatfic
subroutine pk2isomatfic(FIC, DISO, CBAR, CBARI1, UNIT2, NDI)
Definition: pk2isomatfic.for:2
sigiso
subroutine sigiso(SISO, SFIC, PE, NDI)
Definition: sigiso.for:2
setvol
subroutine setvol(CVOL, PV, PPV, UNIT2, UNIT4S, NDI)
Definition: setvol.for:2
Generated by   doxygen 1.8.11