C **********************************************************
C **
C ** ABAQUS User Subroutines
C **
C **********************************************************
C ** Created by James Innes
C **********************************************************
C **
C **  DLOAD - define a non-uniform load ********************
C****************************************************************
C   *** ABAQUS UMAT SUBROUTINE ***
C       for 3D elastic analysis
C       output the elastic stress solutions and the temperature field
C       with element number and integration point number 
C
CC      STATEV(1) - Von Mises effective stress	
CC      STATEV(2) - The temperature field
C************************************************************************************

       SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD,
     1 RPL, DDSDDT, DRPLDE, DRPLDT, 
     2 STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, CMNAME,
     3 NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT,
     4 CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, KSTEP, KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*8 CMNAME
      DIMENSION STRESS(NTENS), STATEV(NSTATV),
     1 DDSDDE(NTENS, NTENS), DDSDDT(NTENS), DRPLDE(NTENS),
     2 STRAN(NTENS), DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1),
     3 PROPS(NPROPS), COORDS(3), DROT(3,3), DFGRD0(3,3), DFGRD1(3,3)
C
** 
C     NEL - The total number of elements
C     NUPT - The maximun number of integration points in the element 
C     You can change the values of NEL and NUPT for your case.
 
      PARAMETER(NEL=5422, NUPT=9)
**
C
      DOUBLE PRECISION YM
      DIMENSION DSTRESS(NTENS), KMARK(NEL,NUPT), IMARK(NEL,NUPT),
     1          E(10),TVALS(10)
      LOGICAL FLAG
C
      SAVE
C
**
C     You can change the following three file names for your case as well as the working directory.   
C     The first file is for checking the maximum elastic stress
C     The second file is for saving the elastic solutions, which used for further computations.      
C     The third file is for saving the temperature field if applied, which used for further computations.      
 
      OPEN(UNIT=16,FILE='/elastic3d_check_holed.dat')
      OPEN(UNIT=17,FILE='/elastic3dP_holed.dat',
     1 STATUS='UNKNOWN',ACCESS='DIRECT',FORM='FORMATTED',RECL=17)
      OPEN(UNIT=19,FILE='/em3dP_holed.dat',
     1 STATUS='UNKNOWN',ACCESS='DIRECT',FORM='FORMATTED',RECL=17)
**
C
C INITIALISING FLAG
C
C     YM=PROPS(1) - Constant Young's modulus. 
C          For temperature-depedent Young's modulus, let PROPS(1)=0.
C     PNU - Poisson's ratio
C
      YM=PROPS(1) 
      PNU=PROPS(2)

      IF(KINC.NE.KINI) THEN
        KINI=KINC
        DO I1=1, NEL
           DO J1=1, NUPT
              KMARK(I1,J1)=0
           END DO
        END DO
      END IF      
C
C INITIALISING STRESSES
C
      DO K1=1, NTENS
         DSTRESS(K1)=0
      END DO
C

      IF(YM.EQ.0) THEN 
**
C     Define temperature-dependent Young's Modulus
      
      LENGTH=7 			
      TMIN=0.0
      TMAX=750.0
      E(1)=1.98298E+5
      E(2)=1.9229E+5
      E(3)=1.7727E+5
      E(4)=1.6976E+5
      E(5)=1.6225E+5
      E(6)=1.45E+5
      E(7)=1.4723E+5
      TVALS(1)=20
      TVALS(2)=100
      TVALS(3)=300
      TVALS(4)=400
      TVALS(5)=500
      TVALS(6)=600
      TVALS(7)=700
**
C
CC    Interpolate linearly to determine Young's modulus at TVAL
C
      CALL ZZPRIN(LENGTH,E,TVALS,TMAX,TMIN,DTEMP,FLAG,YM)

      END IF

C
C DEFINING TERMS FOR JACOBIAN
C
      TERM4 =(1.0+PNU)*(1.0-2*PNU) 
C
      TERM1 = YM*(1.0-PNU)/TERM4
C
      TERM2 = YM*PNU/TERM4
C
      TERM3 = YM/2./(1.0+PNU)
C
C FILLING IN JACOBIAN
C
      DO K1=1,NDI
         DDSDDE(K1,K1)=TERM1
      END DO
C
      DO K1=2,NDI
         N2=K1-1
         DO K2=1,N2
            DDSDDE(K2,K1) = TERM2 
            DDSDDE(K1,K2) = TERM2
         END DO
      END DO
C
      DO K1=NDI+1,NTENS 	
         DDSDDE(K1,K1)=TERM3
      END DO
C COMPUTING STRESS FIELD
C
        DO K=1,NTENS
         DO I=1,NTENS
            DSTRESS(K) = DSTRESS(K)+DDSDDE(K,I)*DSTRAN(I)
         END DO
      END DO
C
C UPDATING STRESS TENSOR
C
      DO K=1, NTENS
         STRESS(K) = STRESS(K)+DSTRESS(K)
      END DO
C
C WRITING SOLUTION TO FILE
C
      IF (KMARK(NOEL,NPT).EQ.0) GO TO 100
         IF (IMARK(NOEL,NPT).EQ.1) GO TO 100
            KCOUNTT=KCOUNTT+2
	    WRITE(19,170,REC=KCOUNTT-1)NOEL, NPT
	    KCOUNTT=KCOUNTT+1
            WRITE(19,180,REC=KCOUNTT)DTEMP
 	    
            KCOUNTS=KCOUNTS+2
            WRITE(17,170,REC=KCOUNTS-1)NOEL, NPT
	    DO K=1,NTENS     
	    KCOUNTS=KCOUNTS+1
            WRITE(17,180,REC=KCOUNTS)DSTRESS(K)
            END DO
            IMARK(NOEL,NPT)=1
 100  CONTINUE
C CHF chf
      HP=(DSTRESS(1)+DSTRESS(2)+DSTRESS(3))/3.0
      V=1.5*((DSTRESS(1)-HP)**2+(DSTRESS(2)-HP)**2+(DSTRESS(3)-HP)**2
     1   +2*DSTRESS(4)**2+2*DSTRESS(5)**2+2*DSTRESS(6)**2)
      VMSTRESS=DSQRT(V)
      STATEV(1)=DSQRT(V)
      STATEV(2)=DTEMP 
C
C SETTING FLAG
C
      KMARK(NOEL,NPT)=1
C
170    FORMAT(I17)
180    FORMAT(E17.10)
 
C
      RETURN
      END
C
C
      SUBROUTINE URDFIL(LSTOP,LOVRWRT,KSTEP,KINC,DTIME,TIME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION ARRAY(513),JRRAY(NPRECD,513),TIME(2)
      EQUIVALENCE (ARRAY(1),JRRAY(1,1))
C
 
C
            TSTRESS=0.0
            XSTRESS=0.0
            YSTRESS=0.0
            ZSTRESS=0.0
C
      CALL POSFIL(KSTEP, KINC, ARRAY, JRCD)
C
      DO 1000 K=1,999999
         CALL DBFILE(0, ARRAY, JRCD)
         IF(JRCD.NE.0) GO TO 10
         KEY=JRRAY(1,2)
C
C
         IF(KEY.EQ.1) THEN
            JEL=JRRAY(1,3)
            JPNT=JRRAY(1,4)
C
         ELSE IF(KEY.EQ.11) THEN
            IF (ARRAY(3).GT.XSTRESS) THEN
               XSTRESS=ARRAY(3)
               IELE=JEL
               IINT=JPNT
            ELSE IF (ARRAY(4).GT.YSTRESS) THEN
               YSTRESS=ARRAY(4)
               JELE=JEL
               JINT=JPNT
            ELSE IF (ARRAY(5).GT.ZSTRESS) THEN
               ZSTRESS=ARRAY(5)
               KELE=JEL
               KINT=JPNT
            END IF 
C
         ELSE IF(KEY.EQ.12) THEN
            IF (ARRAY(3).GT.TSTRESS) THEN
               TSTRESS=ARRAY(3)
               LELE=JEL
               LINT=JPNT
            END IF
         END IF
 1000 CONTINUE
 10   CONTINUE

C
120   FORMAT('ELEMENT: ',I5,' POINT: ',I3)
130   FORMAT('--------------------------------')
   
      WRITE(16,*)'KINC=',KINC
      WRITE(16,*)'URDFIL'
      WRITE(16,*)'HIGEST EFFECTIVE STRESS:'
      WRITE(16,*)TSTRESS
      WRITE(16,*)'LOCATED AT:'
      WRITE(16,120)LELE,LINT
      WRITE(16,130)
      WRITE(16,*)'HIGHEST STRESS IN X-DIRECTION:'
      WRITE(16,*)XSTRESS
      WRITE(16,*)'LOCATED AT:'
      WRITE(16,120)IELE,IINT
      WRITE(16,130)
      WRITE(16,*)'HIGHEST STRESS IN Y-DIRECTION:'
      WRITE(16,*)YSTRESS
      WRITE(16,*)'LOCATED AT:'
      WRITE(16,120)JELE,JINT
      WRITE(16,130)
      WRITE(16,*)'HIGHEST STRESS IN Z-DIRECTION:'
      WRITE(16,*)ZSTRESS
      WRITE(16,*)'LOCATED AT:'
      WRITE(16,120)KELE,KINT
     
      RETURN
      END

CC======================================================================
C
      SUBROUTINE ZZPRIN(LENGTH,PVALS,TVALS,TMAX,TMIN,VALT,FLAG,PROP)
C
C     Returns the PRoperty value PROP corresponding to the input
C     temperature TVAL by linear INterpolation from property values in
C     PVALS and corresponding temperature values in TVALS.   
C     Extrapolation is allowed but FLAG is set if TVAL lies outside
C     the range TMIN to TMAX.
C
C     Input Arguments.
C
C     LENGTH - length of arrays TVALS and PVALS, Integer
C     TMAX - maximum temperature for valid extrapolation, DP
C     TMIN - minimum temperature for valid extrapolation, DP
C     TVALS(LENGTH) - temperature value list, DP
C     PVALS(LENGTH) - property value list, DP
C     VALT - value of temperature input, DP
C
C     Output Arguments.
C
C     FLAG - flag set if allowed temperature range violated, DP
C     PROP - value of property output, DP
C
C     Local Variables.
C
C     FRACT - interpolation fraction, DP
C     IFLOW - low index for interpolation, Integer
C
C     Declarations.
C
      DOUBLE PRECISION FRACT, PVALS, PROP, TMAX, TMIN, TVALS, VALT
C
      INTEGER LENGTH, ILOW
      
C
      LOGICAL FLAG
      
C
      DIMENSION PVALS(LENGTH), TVALS(LENGTH)
C
CC     Check for out-of-bounds temperature.
C
      FLAG = .FALSE.
C
      IF ((VALT.LT.TMIN).OR.(VALT.GT.TMAX)) THEN
         FLAG = .TRUE.
      ENDIF
C
C     If LENGTH is 1 return the single value
C
      IF (LENGTH .EQ. 1) THEN
        PROP = PVALS(1)
        RETURN
      ENDIF
C
C     Find low entry for inter(extra)polation.
C
C     Set as first.
C
      ILOW=1
C
C     Check for value below or on bottom of interpolation range.
C
      IF (VALT.LE.TVALS(1)) THEN
         GO TO 200
      ENDIF
C
C     Check for value above interpolation range.
C
      IF (VALT.GT.TVALS(LENGTH)) THEN
         ILOW = LENGTH - 1
         GO TO 200
      ENDIF
C
C     Loop to find ILOW=1,......LENGTH-1 such that 
C     (comment) TVALS(ILOW).LE.VALT.LT.TVALS(ILOW+1)
C
100   CONTINUE
C       
      IF (VALT.GT.TVALS(ILOW + 1)) THEN
         ILOW = ILOW + 1
         GO TO 100
      ENDIF
C
200   CONTINUE
C
CC     Perform inter(extra)polation.
C
      FRACT = (VALT - TVALS(ILOW))/(TVALS(ILOW + 1) - TVALS(ILOW))
      PROP = PVALS(ILOW) + FRACT*(PVALS(ILOW + 1) - PVALS(ILOW))
C
      RETURN
      END
C
C ************************************************

C
      SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,
     1 COORDS,JLTYP,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION TIME(2), COORDS (3)
      CHARACTER*80 SNAME
C    
C **********************************************************
C Force increases along the x-axis.
C Max Force is excerted at the radius of the hole
C Force distribution is non-linear
C A=Constant, Xr=Max X-Coord of radius, Fmax=Max Force 
C A * Xr^2 = Fmax
C In our case, givens are: Xr=15, Fmax=50
C Therefore: A=50/(15^2)
C Giving our force equation as:
C F = A * X^2
C When A = 50/(15^2)
C **********************************************************
C
	Fmax = 50             !** Max force
	R = 15                !** Radius of circle
	A = Fmax/(R**2)       !** Constant
	F = A*(COORDS(1)**2)  !** Force applied
C
      RETURN
      END                    
C
C **********************************************************
C **
C **  UTEMP - define a non-uniform temp' distro' ***********
C
      SUBROUTINE UTEMP(TEMP,NSECPT,KSTEP,KINC,TIME,NODE,COORDS)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION TEMP(NSECPT), TIME(2), COORDS(3)
C
	TEMP(1) = SQRT( COORDS(1)**2 + COORDS(2)**2 )
C
      RETURN
      END