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 sigisomatfic(SFIC,PKFIC,F,DET,NDI)
804 C> ISOTROPIC MATRIX: ISOCHORIC CAUCHY STRESS
805  IMPLICIT NONE
806  include 'PARAM_UMAT.INC'
807 C
808  INTEGER NDI
809  DOUBLE PRECISION SFIC(ndi,ndi),F(ndi,ndi),
810  1 pkfic(ndi,ndi)
811  DOUBLE PRECISION DET
812 C
813  CALL push2(sfic,pkfic,f,det,ndi)
814 C
815  RETURN
816  END SUBROUTINE sigisomatfic
817  SUBROUTINE pk2anisomatfic(AFIC,DANISO,CBAR,INV4,ST0,NDI)
818 C> ANISOTROPIC MATRIX: 2PK 'FICTICIOUS' STRESS TENSOR
819 C INPUT:
820 C DANISO - ANISOTROPIC STRAIN-ENERGY DERIVATIVES
821 C CBAR - DEVIATORIC LEFT CAUCHY-GREEN TENSOR
822 C INV1,INV4 -CBAR INVARIANTS
823 C UNIT2 - 2ND ORDER IDENTITY TENSOR
824 C OUTPUT:
825 C AFIC - 2ND PIOLA KIRCHOOF 'FICTICIOUS' STRESS TENSOR
826  IMPLICIT NONE
827  include 'PARAM_UMAT.INC'
828 C
829  INTEGER NDI
830  DOUBLE PRECISION AFIC(ndi,ndi),DANISO(3),CBAR(3,3)
831  DOUBLE PRECISION DUDI4,DI4DC(3,3),INV4
832  DOUBLE PRECISION ST0(3,3)
833 C
834 C
835 C-----------------------------------------------------------------------------
836  !FIRST DERIVATIVE OF SSEANISO IN ORDER TO I4
837  dudi4=daniso(1)
838 C
839  di4dc=st0
840 C
841  afic=two*(dudi4*di4dc)
842 C
843  RETURN
844  END SUBROUTINE pk2anisomatfic
845  SUBROUTINE stretch(C,B,U,V,NDI)
846 C> STRETCH TENSORS
847  IMPLICIT NONE
848  include 'PARAM_UMAT.INC'
849 C
850  INTEGER NDI
851  DOUBLE PRECISION C(ndi,ndi),B(ndi,ndi),U(ndi,ndi),V(ndi,ndi)
852  DOUBLE PRECISION EIGVAL(ndi),OMEGA(ndi),EIGVEC(ndi,ndi)
853 C
854  CALL spectral(c,omega,eigvec)
855 C
856  eigval(1) = dsqrt(omega(1))
857  eigval(2) = dsqrt(omega(2))
858  eigval(3) = dsqrt(omega(3))
859 C
860  u(1,1) = eigval(1)
861  u(2,2) = eigval(2)
862  u(3,3) = eigval(3)
863 C
864  u = matmul(matmul(eigvec,u),transpose(eigvec))
865 C
866  CALL spectral(b,omega,eigvec)
867 C
868  eigval(1) = dsqrt(omega(1))
869  eigval(2) = dsqrt(omega(2))
870  eigval(3) = dsqrt(omega(3))
871 C write(*,*) eigvec(1,1),eigvec(2,1),eigvec(3,1)
872 C
873  v(1,1) = eigval(1)
874  v(2,2) = eigval(2)
875  v(3,3) = eigval(3)
876 C
877  v = matmul(matmul(eigvec,v),transpose(eigvec))
878  RETURN
879  END SUBROUTINE stretch
880  SUBROUTINE matinv3d(A,A_INV,NDI)
881 C> INVERSE OF A 3X3 MATRIX
882 C RETURN THE INVERSE OF A(3,3) - A_INV
883  IMPLICIT NONE
884 C
885  INTEGER NDI
886 C
887  DOUBLE PRECISION A(ndi,ndi),A_INV(ndi,ndi),DET_A,DET_A_INV
888 C
889  det_a = a(1,1)*(a(2,2)*a(3,3) - a(3,2)*a(2,3)) -
890  + a(2,1)*(a(1,2)*a(3,3) - a(3,2)*a(1,3)) +
891  + a(3,1)*(a(1,2)*a(2,3) - a(2,2)*a(1,3))
892 
893  IF (det_a .LE. 0.d0) THEN
894  WRITE(*,*) 'WARNING: SUBROUTINE MATINV3D:'
895  WRITE(*,*) 'WARNING: DET OF MAT=',det_a
896  RETURN
897  END IF
898 C
899  det_a_inv = 1.d0/det_a
900 C
901  a_inv(1,1) = det_a_inv*(a(2,2)*a(3,3)-a(3,2)*a(2,3))
902  a_inv(1,2) = det_a_inv*(a(3,2)*a(1,3)-a(1,2)*a(3,3))
903  a_inv(1,3) = det_a_inv*(a(1,2)*a(2,3)-a(2,2)*a(1,3))
904  a_inv(2,1) = det_a_inv*(a(3,1)*a(2,3)-a(2,1)*a(3,3))
905  a_inv(2,2) = det_a_inv*(a(1,1)*a(3,3)-a(3,1)*a(1,3))
906  a_inv(2,3) = det_a_inv*(a(2,1)*a(1,3)-a(1,1)*a(2,3))
907  a_inv(3,1) = det_a_inv*(a(2,1)*a(3,2)-a(3,1)*a(2,2))
908  a_inv(3,2) = det_a_inv*(a(3,1)*a(1,2)-a(1,1)*a(3,2))
909  a_inv(3,3) = det_a_inv*(a(1,1)*a(2,2)-a(2,1)*a(1,2))
910 C
911  RETURN
912  END SUBROUTINE matinv3d
913  SUBROUTINE spectral(A,D,V)
914 C> EIGENVALUES AND EIGENVECTOR OF A 3X3 MATRIX
915 C THIS SUBROUTINE CALCULATES THE EIGENVALUES AND EIGENVECTORS OF
916 C A SYMMETRIC 3X3 MATRIX A.
917 C
918 C THE OUTPUT CONSISTS OF A VECTOR D CONTAINING THE THREE
919 C EIGENVALUES IN ASCENDING ORDER, AND A MATRIX V WHOSE
920 C COLUMNS CONTAIN THE CORRESPONDING EIGENVECTORS.
921 C
922  IMPLICIT NONE
923 C
924  INTEGER NP,NROT
925  parameter(np=3)
926 C
927  DOUBLE PRECISION D(3),V(3,3),A(3,3),E(3,3)
928 C
929  e = a
930 C
931  CALL jacobi(e,3,np,d,v,nrot)
932  CALL eigsrt(d,v,3,np)
933 C
934  RETURN
935  END SUBROUTINE spectral
936 
937 C***********************************************************************
938 
939  SUBROUTINE jacobi(A,N,NP,D,V,NROT)
940 C
941 C COMPUTES ALL EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
942 C MATRIX A, WHICH IS OF SIZE N BY N, STORED IN A PHYSICAL
943 C NP BY NP ARRAY. ON OUTPUT, ELEMENTS OF A ABOVE THE DIAGONAL
944 C ARE DESTROYED, BUT THE DIAGONAL AND SUB-DIAGONAL ARE UNCHANGED
945 C AND GIVE FULL INFORMATION ABOUT THE ORIGINAL SYMMETRIC MATRIX.
946 C VECTOR D RETURNS THE EIGENVALUES OF A IN ITS FIRST N ELEMENTS.
947 C V IS A MATRIX WITH THE SAME LOGICAL AND PHYSICAL DIMENSIONS AS
948 C A WHOSE COLUMNS CONTAIN, UPON OUTPUT, THE NORMALIZED
949 C EIGENVECTORS OF A. NROT RETURNS THE NUMBER OF JACOBI ROTATION
950 C WHICH WERE REQUIRED.
951 C
952 C THIS SUBROUTINE IS TAKEN FROM 'NUMERICAL RECIPES.'
953 C
954  IMPLICIT NONE
955  include 'PARAM_UMAT.INC'
956 C
957  INTEGER IP,IQ,N,NMAX,NP,NROT,I,J
958  parameter(nmax=100)
959 C
960  DOUBLE PRECISION A(np,np),D(np),V(np,np),B(nmax),Z(nmax),
961  + sm,tresh,g,t,h,theta,s,c,tau
962 
963 
964 C INITIALIZE V TO THE IDENTITY MATRIX
965  DO i=1,3
966  v(i,i)=one
967  DO j=1,3
968  IF (i.NE.j)THEN
969  v(i,j)=zero
970  ENDIF
971  END DO
972  END DO
973 C INITIALIZE B AND D TO THE DIAGONAL OF A, AND Z TO ZERO.
974 C THE VECTOR Z WILL ACCUMULATE TERMS OF THE FORM T*A_PQ AS
975 C IN EQUATION (11.1.14)
976 C
977  DO ip = 1,n
978  b(ip) = a(ip,ip)
979  d(ip) = b(ip)
980  z(ip) = 0.d0
981  END DO
982 
983 
984 C BEGIN ITERATION
985 C
986  nrot = 0
987  DO i=1,50
988 C
989 C SUM OFF-DIAGONAL ELEMENTS
990 C
991  sm = 0.d0
992  DO ip=1,n-1
993  DO iq=ip+1,n
994  sm = sm + dabs(a(ip,iq))
995  END DO
996  END DO
997 C
998 C IF SM = 0., THEN RETURN. THIS IS THE NORMAL RETURN,
999 C WHICH RELIES ON QUADRATIC CONVERGENCE TO MACHINE
1000 C UNDERFLOW.
1001 C
1002  IF (sm.EQ.0.d0) RETURN
1003 C
1004 C IN THE FIRST THREE SWEEPS CARRY OUT THE PQ ROTATION ONLY IF
1005 C |A_PQ| > TRESH, WHERE TRESH IS SOME THRESHOLD VALUE,
1006 C SEE EQUATION (11.1.25). THEREAFTER TRESH = 0.
1007 C
1008  IF (i.LT.4) THEN
1009  tresh = 0.2d0*sm/n**2
1010  ELSE
1011  tresh = 0.d0
1012  END IF
1013 C
1014  DO ip=1,n-1
1015  DO iq=ip+1,n
1016  g = 100.d0*dabs(a(ip,iq))
1017 C
1018 C AFTER FOUR SWEEPS, SKIP THE ROTATION IF THE
1019 C OFF-DIAGONAL ELEMENT IS SMALL.
1020 C
1021  IF ((i.GT.4).AND.(dabs(d(ip))+g.EQ.dabs(d(ip)))
1022  + .AND.(dabs(d(iq))+g.EQ.dabs(d(iq)))) THEN
1023  a(ip,iq) = 0.d0
1024  ELSE IF (dabs(a(ip,iq)).GT.tresh) THEN
1025  h = d(iq) - d(ip)
1026  IF (dabs(h)+g.EQ.dabs(h)) THEN
1027 C
1028 C T = 1./(2.*THETA), EQUATION (11.1.10)
1029 C
1030  t =a(ip,iq)/h
1031  ELSE
1032  theta = 0.5d0*h/a(ip,iq)
1033  t =1.d0/(dabs(theta)+dsqrt(1.d0+theta**2.d0))
1034  IF (theta.LT.0.d0) t = -t
1035  END IF
1036  c = 1.d0/dsqrt(1.d0 + t**2.d0)
1037  s = t*c
1038  tau = s/(1.d0 + c)
1039  h = t*a(ip,iq)
1040  z(ip) = z(ip) - h
1041  z(iq) = z(iq) + h
1042  d(ip) = d(ip) - h
1043  d(iq) = d(iq) + h
1044  a(ip,iq) = 0.d0
1045 C
1046 C CASE OF ROTATIONS 1 <= J < P
1047 C
1048  DO j=1,ip-1
1049  g = a(j,ip)
1050  h = a(j,iq)
1051  a(j,ip) = g - s*(h + g*tau)
1052  a(j,iq) = h + s*(g - h*tau)
1053  END DO
1054 C
1055 C CASE OF ROTATIONS P < J < Q
1056 C
1057  DO j=ip+1,iq-1
1058  g = a(ip,j)
1059  h = a(j,iq)
1060  a(ip,j) = g - s*(h + g*tau)
1061  a(j,iq) = h + s*(g - h*tau)
1062  END DO
1063 C
1064 C CASE OF ROTATIONS Q < J <= N
1065 C
1066  DO j=iq+1,n
1067  g = a(ip,j)
1068  h = a(iq,j)
1069  a(ip,j) = g - s*(h + g*tau)
1070  a(iq,j) = h + s*(g - h*tau)
1071  END DO
1072  DO j = 1,n
1073  g = v(j,ip)
1074  h = v(j,iq)
1075  v(j,ip) = g - s*(h + g*tau)
1076  v(j,iq) = h + s*(g - h*tau)
1077  END DO
1078  nrot = nrot + 1
1079  END IF
1080  END DO
1081  END DO
1082 C
1083 C UPDATE D WITH THE SUM OF T*A_PQ, AND REINITIALIZE Z
1084 C
1085  DO ip=1,n
1086  b(ip) = b(ip) + z(ip)
1087  d(ip) = b(ip)
1088  z(ip) = 0.d0
1089  END DO
1090  END DO
1091 C
1092 C IF THE ALGORITHM HAS REACHED THIS STAGE, THEN THERE
1093 C ARE TOO MANY SWEEPS. PRINT A DIAGNOSTIC AND CUT THE
1094 C TIME INCREMENT.
1095 C
1096  WRITE (*,'(/1X,A/)') '50 ITERATIONS IN JACOBI SHOULD NEVER HAPPEN'
1097 C
1098  RETURN
1099  END SUBROUTINE jacobi
1100 
1101 C**********************************************************************
1102  SUBROUTINE eigsrt(D,V,N,NP)
1103 C
1104 C GIVEN THE EIGENVALUES D AND EIGENVECTORS V AS OUTPUT FROM
1105 C JACOBI, THIS SUBROUTINE SORTS THE EIGENVALUES INTO ASCENDING
1106 C ORDER AND REARRANGES THE COLMNS OF V ACCORDINGLY.
1107 C
1108 C THE SUBROUTINE WAS TAKEN FROM 'NUMERICAL RECIPES.'
1109 C
1110  IMPLICIT NONE
1111 C
1112  INTEGER N,NP,I,J,K
1113 C
1114  DOUBLE PRECISION D(np),V(np,np),P
1115 C
1116  DO i=1,n-1
1117  k = i
1118  p = d(i)
1119  DO j=i+1,n
1120  IF (d(j).GE.p) THEN
1121  k = j
1122  p = d(j)
1123  END IF
1124  END DO
1125  IF (k.NE.i) THEN
1126  d(k) = d(i)
1127  d(i) = p
1128  DO j=1,n
1129  p = v(j,i)
1130  v(j,i) = v(j,k)
1131  v(j,k) = p
1132  END DO
1133  END IF
1134  END DO
1135 C
1136  RETURN
1137  END SUBROUTINE eigsrt
1138  SUBROUTINE tensorprod2(A,B,C,NDI)
1139 C
1140  Implicit None
1141 C
1142  INTEGER I,J,K,L,NDI
1143 C
1144  DOUBLE PRECISION A(ndi,ndi),B(ndi,ndi),C(ndi,ndi,ndi,ndi)
1145 C
1146  DO i=1,ndi
1147  DO j=1,ndi
1148  DO k=1,ndi
1149  DO l=1,ndi
1150  c(i,j,k,l)=a(i,j)*b(k,l)
1151  END DO
1152  END DO
1153  END DO
1154  END DO
1155 C
1156  RETURN
1157 C
1158  end SUBROUTINE tensorprod2
1159  SUBROUTINE contraction24(S,LT,RT,NDI)
1160 C> DOUBLE CONTRACTION BETWEEN 4TH ORDER AND 2ND ORDER TENSOR
1161 C> INPUT:
1162 C> LT - RIGHT 2ND ORDER TENSOR
1163 C> RT - LEFT 4TH ODER TENSOR
1164 C> OUTPUT:
1165 C> S - DOUBLE CONTRACTED TENSOR (2ND ORDER)
1166  IMPLICIT NONE
1167  include 'PARAM_UMAT.INC'
1168 C
1169  INTEGER I1,J1,K1,L1,NDI
1170 C
1171  DOUBLE PRECISION LT(ndi,ndi),RT(ndi,ndi,ndi,ndi)
1172  DOUBLE PRECISION S(ndi,ndi)
1173  DOUBLE PRECISION AUX
1174 C
1175  DO k1=1,ndi
1176  DO l1=1,ndi
1177  aux=zero
1178  DO i1=1,ndi
1179  DO j1=1,ndi
1180  aux=aux+lt(k1,l1)*rt(i1,j1,k1,l1)
1181  END DO
1182  END DO
1183  s(k1,l1)=aux
1184  END DO
1185  END DO
1186  RETURN
1187  END SUBROUTINE contraction24
1188  SUBROUTINE invariants(A,INV1,INV2,NDI)
1189 C> 1ST AND 2ND INVARIANTS OF A TENSOR
1190  IMPLICIT NONE
1191  include 'PARAM_UMAT.INC'
1192 C
1193  INTEGER NDI,I1
1194  DOUBLE PRECISION A(ndi,ndi),AA(ndi,ndi)
1195  DOUBLE PRECISION INV1,INV1AA, INV2
1196 C
1197  inv1=zero
1198  inv1aa=zero
1199  aa=matmul(a,a)
1200  DO i1=1,ndi
1201  inv1=inv1+a(i1,i1)
1202  inv1aa=inv1aa+aa(i1,i1)
1203  END DO
1204  inv2=(one/two)*(inv1*inv1-inv1aa)
1205 C
1206  RETURN
1207  END SUBROUTINE invariants
1208  SUBROUTINE sigiso(SISO,SFIC,PE,NDI)
1209 C> ISOCHORIC CAUCHY STRESS
1210  IMPLICIT NONE
1211  include 'PARAM_UMAT.INC'
1212 C
1213  INTEGER NDI
1214  DOUBLE PRECISION SISO(ndi,ndi),
1215  1 pe(ndi,ndi,ndi,ndi),sfic(ndi,ndi)
1216 C
1217  CALL contraction42(siso,pe,sfic,ndi)
1218 C
1219  RETURN
1220  END SUBROUTINE sigiso
1221  SUBROUTINE vol(SSEV,PV,PPV,K,DET)
1222 C> VOLUMETRIC CONTRIBUTION :STRAIN ENERGY FUNCTION AND DERIVATIVES
1223  IMPLICIT NONE
1224  include 'PARAM_UMAT.INC'
1225 C
1226  DOUBLE PRECISION SSEV,PV,PPV
1227  DOUBLE PRECISION K,G,DET,AUX
1228 C
1229  g=(one/four)*(det*det-one-two*log(det))
1230 C
1231  ssev=k*g
1232 C
1233  pv=k*(one/two)*(det-one/det)
1234  aux=k*(one/two)*(one+one/(det*det))
1235  ppv=pv+det*aux
1236 C
1237  RETURN
1238  END SUBROUTINE vol
1239  SUBROUTINE push4(SPATIAL,MAT,F,DET,NDI)
1240 C> PIOLA TRANSFORMATION
1241 C> INPUT:
1242 C> MAT - MATERIAL ELASTICITY TENSOR
1243 C> F - DEFORMATION GRADIENT
1244 C> DET - DEFORMATION DETERMINANT
1245 C> OUTPUT:
1246 C> SPATIAL - SPATIAL ELASTICITY TENSOR
1247  IMPLICIT NONE
1248  include 'PARAM_UMAT.INC'
1249 C
1250  INTEGER I1,J1,K1,L1,II1,JJ1,KK1,LL1,NDI
1251 C
1252  DOUBLE PRECISION MAT(ndi,ndi,ndi,ndi),F(ndi,ndi)
1253  DOUBLE PRECISION SPATIAL(ndi,ndi,ndi,ndi)
1254  DOUBLE PRECISION AUX,DET
1255 C
1256  DO i1=1,ndi
1257  DO j1=1,ndi
1258  DO k1=1,ndi
1259  DO l1=1,ndi
1260  aux=zero
1261  DO ii1=1,ndi
1262  DO jj1=1,ndi
1263  DO kk1=1,ndi
1264  DO ll1=1,ndi
1265  aux=aux+(det**(-one))*
1266  + f(i1,ii1)*f(j1,jj1)*
1267  + f(k1,kk1)*f(l1,ll1)*mat(ii1,jj1,kk1,ll1)
1268  END DO
1269  END DO
1270  END DO
1271  END DO
1272  spatial(i1,j1,k1,l1)=aux
1273  END DO
1274  END DO
1275  END DO
1276  END DO
1277 C
1278  RETURN
1279  END SUBROUTINE push4
1280  SUBROUTINE pk2iso(PKISO,PKFIC,PL,DET,NDI)
1281 C> ISOCHORIC PK2 STRESS TENSOR
1282  IMPLICIT NONE
1283  include 'PARAM_UMAT.INC'
1284 C
1285  INTEGER NDI,I1,J1
1286  DOUBLE PRECISION PKISO(ndi,ndi),
1287  1 pl(ndi,ndi,ndi,ndi),pkfic(ndi,ndi)
1288  DOUBLE PRECISION DET,SCALE2
1289 C
1290  CALL contraction42(pkiso,pl,pkfic,ndi)
1291 C
1292  scale2=det**(-two/three)
1293  DO i1=1,ndi
1294  DO j1=1,ndi
1295  pkiso(i1,j1)=scale2*pkiso(i1,j1)
1296  END DO
1297  END DO
1298 C
1299  RETURN
1300  END SUBROUTINE pk2iso
1301  SUBROUTINE onem(A,AA,AAS,NDI)
1302 C
1303 C> THIS SUBROUTINE GIVES:
1304 C> 2ND ORDER IDENTITY TENSORS - A
1305 C> 4TH ORDER IDENTITY TENSOR - AA
1306 C> 4TH ORDER SYMMETRIC IDENTITY TENSOR - AAS
1307 C
1308  IMPLICIT NONE
1309  include 'PARAM_UMAT.INC'
1310 C
1311  INTEGER I,J,K,L,NDI
1312 C
1313  DOUBLE PRECISION A(ndi,ndi),AA(ndi,ndi,ndi,ndi),
1314  1 aas(ndi,ndi,ndi,ndi)
1315 C
1316  DO i=1,ndi
1317  DO j=1,ndi
1318  IF (i .EQ. j) THEN
1319  a(i,j) = one
1320  ELSE
1321  a(i,j) = zero
1322  END IF
1323  END DO
1324  END DO
1325 C
1326  DO i=1,ndi
1327  DO j=1,ndi
1328  DO k=1,ndi
1329  DO l=1,ndi
1330  aa(i,j,k,l)=a(i,k)*a(j,l)
1331  aas(i,j,k,l)=(one/two)*(a(i,k)*a(j,l)+a(i,l)*a(j,k))
1332  END DO
1333  END DO
1334  END DO
1335  END DO
1336 C
1337  RETURN
1338  END SUBROUTINE onem
1339  SUBROUTINE metiso(CMISO,CMFIC,PL,PKISO,PKFIC,C,UNIT2,DET,NDI)
1340 C> ISOCHORIC MATERIAL ELASTICITY TENSOR
1341  IMPLICIT NONE
1342  include 'PARAM_UMAT.INC'
1343 C
1344  INTEGER NDI,I1,J1,K1,L1
1345  DOUBLE PRECISION UNIT2(ndi,ndi),PL(ndi,ndi,ndi,ndi),
1346  1 cmiso(ndi,ndi,ndi,ndi),pkiso(ndi,ndi),
1347  2 cmfic(ndi,ndi,ndi,ndi),pkfic(ndi,ndi),
1348  3 cisoaux(ndi,ndi,ndi,ndi),
1349  4 cisoaux1(ndi,ndi,ndi,ndi),c(ndi,ndi),
1350  5 plt(ndi,ndi,ndi,ndi),cinv(ndi,ndi),
1351  6 pll(ndi,ndi,ndi,ndi)
1352  DOUBLE PRECISION TRFIC,XX,YY,ZZ,DET,AUX,AUX1
1353 C
1354  CALL matinv3d(c,cinv,ndi)
1355  cisoaux1=zero
1356  cisoaux=zero
1357  CALL contraction44(cisoaux1,pl,cmfic,ndi)
1358 C
1359 C transpose of lagrangian projection tensor
1360  DO i1=1,ndi
1361  DO j1=1,ndi
1362  DO k1=1,ndi
1363  DO l1=1,ndi
1364  plt(i1,j1,k1,l1)=pl(k1,l1,i1,j1)
1365  END DO
1366  END DO
1367  END DO
1368  END DO
1369 C
1370  CALL contraction44(cisoaux,cisoaux1,plt,ndi)
1371 C
1372  trfic=zero
1373  aux=det**(-two/three)
1374  aux1=aux**two
1375  DO i1=1,ndi
1376  trfic=trfic+aux*pkfic(i1,i1)*c(i1,i1)
1377  END DO
1378 C
1379  DO i1=1,ndi
1380  DO j1=1,ndi
1381  DO k1=1,ndi
1382  DO l1=1,ndi
1383  xx=aux1*cisoaux(i1,j1,k1,l1)
1384  pll(i1,j1,k1,l1)=(one/two)*(cinv(i1,k1)*cinv(j1,l1)+
1385  1 cinv(i1,l1)*cinv(j1,k1))-
1386  2 (one/three)*cinv(i1,j1)*cinv(k1,l1)
1387  yy=trfic*pll(i1,j1,k1,l1)
1388  zz=pkiso(i1,j1)*cinv(k1,l1)+cinv(i1,j1)*pkiso(k1,l1)
1389 C
1390  cmiso(i1,j1,k1,l1)=xx+(two/three)*yy-(two/three)*zz
1391  END DO
1392  END DO
1393  END DO
1394  END DO
1395 C
1396  RETURN
1397  END SUBROUTINE metiso
1398  SUBROUTINE setiso(CISO,CFIC,PE,SISO,SFIC,UNIT2,NDI)
1399 C> ISOCHORIC SPATIAL ELASTICITY TENSOR
1400  IMPLICIT NONE
1401  include 'PARAM_UMAT.INC'
1402 C
1403  INTEGER NDI,I1,J1,K1,L1
1404  DOUBLE PRECISION UNIT2(ndi,ndi),PE(ndi,ndi,ndi,ndi),
1405  1 ciso(ndi,ndi,ndi,ndi),siso(ndi,ndi),
1406  2 cfic(ndi,ndi,ndi,ndi),sfic(ndi,ndi),
1407  3 cisoaux(ndi,ndi,ndi,ndi),
1408  4 cisoaux1(ndi,ndi,ndi,ndi)
1409  DOUBLE PRECISION TRFIC,XX,YY,ZZ
1410 C
1411  cisoaux1=zero
1412  cisoaux=zero
1413 
1414  CALL contraction44(cisoaux1,pe,cfic,ndi)
1415  CALL contraction44(cisoaux,cisoaux1,pe,ndi)
1416 C
1417  trfic=zero
1418  DO i1=1,ndi
1419  trfic=trfic+sfic(i1,i1)
1420  END DO
1421 C
1422  DO i1=1,ndi
1423  DO j1=1,ndi
1424  DO k1=1,ndi
1425  DO l1=1,ndi
1426  xx=cisoaux(i1,j1,k1,l1)
1427  yy=trfic*pe(i1,j1,k1,l1)
1428  zz=siso(i1,j1)*unit2(k1,l1)+unit2(i1,j1)*siso(k1,l1)
1429 C
1430  ciso(i1,j1,k1,l1)=xx+(two/three)*yy-(two/three)*zz
1431  END DO
1432  END DO
1433  END DO
1434  END DO
1435 C
1436  RETURN
1437  END SUBROUTINE setiso
1438  SUBROUTINE push2(SIG,PK,F,DET,NDI)
1439 C> PIOLA TRANSFORMATION
1440 C> INPUT:
1441 C> PK - 2ND PIOLA KIRCHOOF STRESS TENSOR
1442 C> F - DEFORMATION GRADIENT
1443 C> DET - DEFORMATION DETERMINANT
1444 C> OUTPUT:
1445 C> SIG - CAUCHY STRESS TENSOR
1446  IMPLICIT NONE
1447  include 'PARAM_UMAT.INC'
1448 C
1449  INTEGER I1,J1,II1,JJ1,NDI
1450  DOUBLE PRECISION PK(ndi,ndi),F(ndi,ndi)
1451  DOUBLE PRECISION SIG(ndi,ndi)
1452  DOUBLE PRECISION AUX,DET
1453 C
1454  DO i1=1,ndi
1455  DO j1=1,ndi
1456  aux=zero
1457  DO ii1=1,ndi
1458  DO jj1=1,ndi
1459  aux=aux+(det**(-one))*f(i1,ii1)*f(j1,jj1)*pk(ii1,jj1)
1460  END DO
1461  END DO
1462  sig(i1,j1)=aux
1463  END DO
1464  END DO
1465 C
1466  RETURN
1467  END SUBROUTINE push2
1468  SUBROUTINE deformation(F,C,B,NDI)
1469 C> RIGHT AND LEFT CAUCHY-GREEN DEFORMATION TENSORS
1470  IMPLICIT NONE
1471  include 'PARAM_UMAT.INC'
1472 C
1473  INTEGER NDI
1474  DOUBLE PRECISION F(ndi,ndi),C(ndi,ndi),B(ndi,ndi)
1475 C RIGHT CAUCHY-GREEN DEFORMATION TENSOR
1476  c=matmul(transpose(f),f)
1477 C LEFT CAUCHY-GREEN DEFORMATION TENSOR
1478  b=matmul(f,transpose(f))
1479  RETURN
1480  END SUBROUTINE deformation
1481  SUBROUTINE fslip(F,FBAR,DET,NDI)
1482 C> DISTORTION GRADIENT
1483  IMPLICIT NONE
1484  include 'PARAM_UMAT.INC'
1485 C
1486  INTEGER NDI,I1,J1
1487  DOUBLE PRECISION F(ndi,ndi),FBAR(ndi,ndi)
1488  DOUBLE PRECISION DET,SCALE1
1489 C
1490 C JACOBIAN DETERMINANT
1491  det = f(1,1) * f(2,2) * f(3,3)
1492  1 - f(1,2) * f(2,1) * f(3,3)
1493 C
1494  IF (ndi .EQ. 3) THEN
1495  det = det + f(1,2) * f(2,3) * f(3,1)
1496  1 + f(1,3) * f(3,2) * f(2,1)
1497  2 - f(1,3) * f(3,1) * f(2,2)
1498  3 - f(2,3) * f(3,2) * f(1,1)
1499  END IF
1500 C
1501  scale1=det**(-one /three)
1502 C
1503  DO i1=1,ndi
1504  DO j1=1,ndi
1505  fbar(i1,j1)=scale1*f(i1,j1)
1506  END DO
1507  END DO
1508 C
1509  RETURN
1510  END SUBROUTINE fslip
1511  SUBROUTINE cmatanisomatfic(CMANISOMATFIC,M0,DANISO,UNIT2,DET,NDI)
1512 C
1513 C> ANISOTROPIC MATRIX: MATERIAL 'FICTICIOUS' ELASTICITY TENSOR
1514  IMPLICIT NONE
1515  include 'PARAM_UMAT.INC'
1516 C
1517  INTEGER NDI,I,J,K,L
1518  DOUBLE PRECISION CMANISOMATFIC(ndi,ndi,ndi,ndi),UNIT2(ndi,ndi),
1519  1 m0(ndi,ndi),daniso(3),det
1520  DOUBLE PRECISION CINV4(ndi,ndi,ndi,ndi),CINV14(ndi,ndi,ndi,ndi)
1521  DOUBLE PRECISION D2UDI4,D2UDI1DI4
1522  DOUBLE PRECISION IMM(ndi,ndi,ndi,ndi),MMI(ndi,ndi,ndi,ndi),
1523  1 mm0(ndi,ndi,ndi,ndi)
1524 C
1525 C-----------------------------------------------------------------------------
1526  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I4
1527  d2udi4=daniso(2)
1528  !2ND DERIVATIVE OF SSEANISO IN ORDER TO I1 AND I4
1529  d2udi1di4=daniso(3)
1530 C
1531  CALL tensorprod2(m0,m0,mm0,ndi)
1532  CALL tensorprod2(unit2,m0,imm,ndi)
1533  CALL tensorprod2(m0,unit2,mmi,ndi)
1534 C
1535  DO i=1,ndi
1536  DO j=1,ndi
1537  DO k=1,ndi
1538  DO l=1,ndi
1539  cinv4(i,j,k,l)=d2udi4*mm0(i,j,k,l)
1540  cinv14(i,j,k,l)=d2udi1di4*(imm(i,j,k,l)+mmi(i,j,k,l))
1541  cmanisomatfic(i,j,k,l)=four*(cinv4(i,j,k,l)+cinv14(i,j,k,l))
1542  END DO
1543  END DO
1544  END DO
1545  END DO
1546 C
1547  RETURN
1548  END SUBROUTINE cmatanisomatfic
1549  SUBROUTINE sdvwrite(DET,LAMBDA,STATEV)
1550 C> VISCOUS DISSIPATION: WRITE STATE VARS
1551  IMPLICIT NONE
1552  include 'PARAM_UMAT.INC'
1553 C
1554  DOUBLE PRECISION STATEV(nsdv),DET,LAMBDA
1555 C write your sdvs here. they should be allocated
1556 C after the viscous terms (check hvwrite)
1557  statev(1)=det
1558  statev(2)=lambda
1559 
1560  RETURN
1561 C
1562  END SUBROUTINE sdvwrite
1563  SUBROUTINE resetdfgrd(DFGRD,NDI)
1564 C
1565  IMPLICIT NONE
1566  include 'PARAM_UMAT.INC'
1567 
1568  INTEGER NDI
1569  DOUBLE PRECISION DFGRD(ndi,ndi)
1570 
1571  dfgrd(1,1)= one
1572  dfgrd(1,2)= zero
1573  dfgrd(1,3)= zero
1574  dfgrd(2,1)= zero
1575  dfgrd(2,2)= one
1576  dfgrd(2,3)= zero
1577  dfgrd(3,1)= zero
1578  dfgrd(3,2)= zero
1579  dfgrd(3,3)= one
1580 
1581  END SUBROUTINE resetdfgrd
1582  SUBROUTINE fibdir(FIB,ST0,ST,NE,NOEL,NDI,VORIF,VD,DISTGR,DFGRD1)
1583 C
1584  IMPLICIT NONE
1585  include 'PARAM_UMAT.INC'
1586 C
1587  INTEGER NDI, NE, NOEL,INOEL,I,J,I1,J1
1588  DOUBLE PRECISION SUM1, DFGRD1(3,3), DNORM
1589  DOUBLE PRECISION VORIF(3),ST(3,3),VD(3),ST0(3,3),DISTGR(3,3)
1590  DOUBLE PRECISION FIB(ne,4)
1591 C
1592  inoel=0
1593  i=0
1594  DO i=1,ne
1595 C ELEMENT IDENTIFICATION
1596  IF(noel.EQ.int(fib(i,1))) THEN
1597  inoel=i
1598  ENDIF
1599  ENDDO
1600 C
1601 C FIB - FIBER ORIENTATION
1602  dnorm=dsqrt(fib(inoel,2)*fib(inoel,2)+
1603  1 fib(inoel,3)*fib(inoel,3)+
1604  2 fib(inoel,4)*fib(inoel,4))
1605 C
1606 C UNDERFORMED FIBER ORIENTATION TENSOR
1607 C
1608  DO i1=1,ndi
1609  j1=i1+1
1610 C FIBER ORIENTATION NORMALIZED VECTOR - FAMILY 1
1611  vorif(i1)=fib(inoel,j1)/dnorm
1612  END DO
1613 C
1614 
1615  DO i=1,ndi
1616  sum1=zero
1617  DO j=1,ndi
1618  sum1=sum1+dfgrd1(i,j)*vorif(j)
1619  ENDDO
1620 C FIBER DIRECTIONS IN THE DEFORMED CONFIGURATION
1621 C -FAMILY 1
1622  vd(i)=sum1
1623  ENDDO
1624  dnorm=dsqrt(vd(1)*vd(1)+
1625  1 vd(2)*vd(2)+
1626  2 vd(3)*vd(3))
1627 C COSINE OF THE ANGLE BETWEEN FIBERS
1628 C
1629 C
1630 C--------------------------------------------------------------------------
1631  DO i=1,ndi
1632  DO j=1,ndi
1633 C STRUCTURAL TENSOR - FAMILY 1
1634  st0(i,j)=vorif(i)*vorif(j)
1635  END DO
1636  END DO
1637 C
1638 C STRUCTURE TENSOR IN THE DEFORMED CONFIGURATION - FAMILY 1
1639  st=matmul(st0,transpose(distgr))
1640  st=matmul(distgr,st)
1641 C
1642 C
1643  RETURN
1644  END SUBROUTINE fibdir
1645  SUBROUTINE indexx(STRESS,DDSDDE,SIG,TNG,NTENS,NDI)
1646 C> INDEXATION: FULL SIMMETRY IN STRESSES AND ELASTICITY TENSORS
1647  IMPLICIT NONE
1648  include 'PARAM_UMAT.INC'
1649 C
1650  INTEGER II1(6),II2(6),NTENS,NDI,I1,J1
1651  DOUBLE PRECISION STRESS(ntens),DDSDDE(ntens,ntens)
1652  DOUBLE PRECISION SIG(ndi,ndi),TNG(ndi,ndi,ndi,ndi)
1653  DOUBLE PRECISION PP1,PP2
1654 C
1655  ii1(1)=1
1656  ii1(2)=2
1657  ii1(3)=3
1658  ii1(4)=1
1659  ii1(5)=1
1660  ii1(6)=2
1661 C
1662  ii2(1)=1
1663  ii2(2)=2
1664  ii2(3)=3
1665  ii2(4)=2
1666  ii2(5)=3
1667  ii2(6)=3
1668 C
1669  DO i1=1,ntens
1670 C STRESS VECTOR
1671  stress(i1)=sig(ii1(i1),ii2(i1))
1672  DO j1=1,ntens
1673 C DDSDDE - FULLY SIMMETRY IMPOSED
1674  pp1=tng(ii1(i1),ii2(i1),ii1(j1),ii2(j1))
1675  pp2=tng(ii1(i1),ii2(i1),ii2(j1),ii1(j1))
1676  ddsdde(i1,j1)=(one/two)*(pp1+pp2)
1677  END DO
1678  END DO
1679 C
1680  RETURN
1681 C
1682  END SUBROUTINE indexx
1683  SUBROUTINE rotation(F,R,U,NDI)
1684 C> COMPUTES ROTATION TENSOR
1685  IMPLICIT NONE
1686  include 'PARAM_UMAT.INC'
1687 C
1688  INTEGER NDI
1689  DOUBLE PRECISION F(ndi,ndi),R(ndi,ndi),U(ndi,ndi),UINV(ndi,ndi)
1690 C
1691  CALL matinv3d(u,uinv,ndi)
1692 C
1693  r = matmul(f,uinv)
1694  RETURN
1695  END SUBROUTINE rotation
1696  SUBROUTINE cmatisomatfic(CMISOMATFIC,CBAR,CBARI1,CBARI2,
1697  1 diso,unit2,unit4,det,ndi)
1698 C> ISOTROPIC MATRIX: MATERIAL 'FICTICIOUS' ELASTICITY TENSOR
1699  IMPLICIT NONE
1700  include 'PARAM_UMAT.INC'
1701 C
1702  INTEGER NDI,I1,J1,K1,L1
1703  DOUBLE PRECISION CMISOMATFIC(ndi,ndi,ndi,ndi),UNIT2(ndi,ndi),
1704  1 cbar(ndi,ndi),diso(5),
1705  2 unit4(ndi,ndi,ndi,ndi)
1706  DOUBLE PRECISION CBARI1,CBARI2
1707  DOUBLE PRECISION DUDI1,DUDI2,D2UD2I1,D2UD2I2,D2UDI1I2
1708  DOUBLE PRECISION AUX,AUX1,AUX2,AUX3,AUX4,DET
1709  DOUBLE PRECISION UIJ,UKL,CIJ,CKL
1710 C
1711  dudi1=diso(1)
1712  dudi2=diso(2)
1713  d2ud2i1=diso(3)
1714  d2ud2i2=diso(4)
1715  d2udi1i2=diso(5)
1716 C
1717  aux1=four*(d2ud2i1+two*cbari1*d2udi1i2+
1718  1 dudi2+cbari1*cbari1*d2ud2i2)
1719  aux2=-four*(d2udi1i2+cbari1*d2ud2i2)
1720  aux3=four*d2ud2i2
1721  aux4=-four*dudi2
1722 
1723  DO i1=1,ndi
1724  DO j1=1,ndi
1725  DO k1=1,ndi
1726  DO l1=1,ndi
1727  uij=unit2(i1,j1)
1728  ukl=unit2(k1,l1)
1729  cij=cbar(i1,j1)
1730  ckl=cbar(k1,l1)
1731  aux=aux1*uij*ukl+
1732  1 aux2*(uij*ckl+cij*ukl)+aux3*cij*ckl+
1733  3 aux4*unit4(i1,j1,k1,l1)
1734  cmisomatfic(i1,j1,k1,l1)=aux
1735  END DO
1736  END DO
1737  END DO
1738  END DO
1739 C
1740  RETURN
1741  END SUBROUTINE cmatisomatfic
1742  SUBROUTINE projlag(C,AA,PL,NDI)
1743 C> LAGRANGIAN PROJECTION TENSOR
1744 C INPUTS:
1745 C IDENTITY TENSORS - A, AA
1746 C ISOCHORIC LEFT CAUCHY GREEN TENSOR - C
1747 C INVERSE OF C - CINV
1748 C OUTPUTS:
1749 C 4TH ORDER SYMMETRIC LAGRANGIAN PROJECTION TENSOR - PL
1750 C
1751  IMPLICIT NONE
1752  include 'PARAM_UMAT.INC'
1753 C
1754  INTEGER I,J,K,L,NDI
1755 C
1756  DOUBLE PRECISION CINV(ndi,ndi),AA(ndi,ndi,ndi,ndi),
1757  1 pl(ndi,ndi,ndi,ndi),c(ndi,ndi)
1758 C
1759  CALL matinv3d(c,cinv,ndi)
1760 C
1761  DO i=1,ndi
1762  DO j=1,ndi
1763  DO k=1,ndi
1764  DO l=1,ndi
1765  pl(i,j,k,l)=aa(i,j,k,l)-(one/three)*(cinv(i,j)*c(k,l))
1766  END DO
1767  END DO
1768  END DO
1769  END DO
1770 C
1771  RETURN
1772  END SUBROUTINE projlag
1773  SUBROUTINE csisomatfic(CISOMATFIC,CMISOMATFIC,DISTGR,DET,NDI)
1774 C> ISOTROPIC MATRIX: SPATIAL 'FICTICIOUS' ELASTICITY TENSOR
1775  IMPLICIT NONE
1776  include 'PARAM_UMAT.INC'
1777 C
1778  INTEGER NDI
1779  DOUBLE PRECISION CMISOMATFIC(ndi,ndi),DISTGR(ndi,ndi),
1780  1 cisomatfic(ndi,ndi,ndi,ndi)
1781  DOUBLE PRECISION DET
1782 C
1783  CALL push4(cisomatfic,cmisomatfic,distgr,det,ndi)
1784 C
1785  RETURN
1786  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