Abaqus umat interface

Author(s): Tero Frondelius

Abstract: making the initial version to call Abaqus umat



In [1]:

    
run(`wget http://www.eng.ox.ac.uk/NP/ICP/plasticity_imp/code_imp.f`)









    



--2015-07-25 22:07:34--  http://www.eng.ox.ac.uk/NP/ICP/plasticity_imp/code_imp.f
Resolving www.eng.ox.ac.uk (www.eng.ox.ac.uk)... 163.1.223.199
Connecting to www.eng.ox.ac.uk (www.eng.ox.ac.uk)|163.1.223.199|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 7458 (7,3K) [text/plain]
Saving to: ‘code_imp.f.1’

     0K .......                                               100%  282M=0s

2015-07-25 22:07:34 (282 MB/s) - ‘code_imp.f.1’ saved [7458/7458]



In [2]:

    
run(`head -30 code_imp.f`)









    



***********************************************************************************
**  UMAT, FOR ABAQUS/STANDARD INCORPORATING ELASTIC-PLASTIC LINEAR               **
**  ISOTROPIC HARDENING. LARGE DEFORMATION FORMULATION FOR PLANE STRAIN          **
**  AND AXI-SYMMETRIC ELEMENTS. IMPLICIT INTEGRATION WITH CONSISTENT JACOBIAN    **
***********************************************************************************
***********************************************************************************
**
**
**
*USER SUBROUTINE
      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*80 CMNAME
C
C
      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
      PARAMETER (M=3,N=3,ID=3,ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,
     +          SIX=6.D0, NINE=9.D0, TOLER=1.D-5)
C

some implicit type castings



In [3]:

    
f = open("ABA_PARAM.INC","w")
write(f,"      implicit real*8(a-h,o-z)\n")
write(f,"      parameter (nprecd=2)\n")
close(f)

Let's compile the shared library



In [4]:

    
run(`gfortran -shared -fPIC -o libumat.so code_imp.f`)

Some Abaqus umat interface definitions

variable	explanation
DDSDDE(NTENS,NTENS)	Jacobian matrix of the constitutive model
STRESS(NTENS)	the stress tensor (in vector format)
STATEV(NSTATV)	An array containing the solution-dependent state variables.
SSE	Specific elastic strain energy
SPD	plastic dissipation
SCD	“creep” dissipation
RPL	Volumetric heat generation per unit time
DDSDDT(NTENS)	Variation of the stress increments with respect to the temperature.
DRPLDE(NTENS)	Variation of RPL with respect to the strain increments.
DRPLDT	Variation of RPL with respect to the temperature.
RPL	RPL is used to indicate whether or not a cohesive element is open to the tangential flow of pore fluid.
PNEWDT	Ratio of suggested new time increment to the time increment being used
STRAN(NTENS)	An array containing the total strains at the beginning of the increment.
DSTRAN(NTENS)	Array of strain increments.
TIME(1)	Value of step time at the beginning of the current increment.
TIME(2)	Value of total time at the beginning of the current increment.
DTIME	Time increment.
TEMP	Temperature at the start of the increment.
DTEMP	Increment of temperature.
PREDEF	Array of interpolated values of predefined field variables at this point at the start of the increment, based on the values read in at the nodes.
DPRED	Array of increments of predefined field variables.
CMNAME	User-defined material name, left justified.
NDI	Number of direct stress components at this point.
NSHR	Number of engineering shear stress components at this point.
NTENS	Size of the stress or strain component array (NDI + NSHR).
NSTATV	Number of solution-dependent state variables that are associated with this material type
PROPS(NPROPS)	User-specified array of material constants associated with this user material.
NPROPS	User-defined number of material constants associated with this user material.
COORDS	An array containing the coordinates of this point.
DROT(3,3)	Rotation increment matrix.
CELENT	Characteristic element length
DFGRD0(3,3)	Array containing the deformation gradient at the beginning of the increment.
DFGRD1(3,3)	Array containing the deformation gradient at the end of the increment.
NOEL	Element number.
NPT	Integration point number.
LAYER	Layer number (for composite shells and layered solids).
KSPT	Section point number within the current layer.
KSTEP	Step number.
KINC	Increment number.



In [5]:

    
STRESS = [0. 0. 0. 0.]
p = 0. # EFFECTIVE PLASTIC STRAIN
r = 0. # ISOTROPIC HARDENING VARIABLE
STATEV = [p r]
NTENS = 4 
DDSDDE = zeros(NTENS,NTENS)
SSE = {} # Not used in this example
SPD = {} # Not used in this example
SCD = {} # Not used in this example
RPL = {} # Not used in this example
DDSDDT = {} # Not used in this example
DRPLDE = {} # Not used in this example
DRPLDT = {} # Not used in this example
STRAN = [0. 0. 0. 0.]
DSTRAN = [0. 0. 0. 0.]
TIME = [0. 0.1] # CHECK TIME(2)
DTIME = {} # Not used in this example
TEMP = {} # Not used in this example
DTEMP = {} # Not used in this example
PREDEF = {} # Not used in this example
DPRED = {} # Not used in this example
CMNAME = {} # Not used in this example CHARACTER*80 CMNAME
NDI = {} # Not used in this example
NSHR = {} # Not used in this example
#NTENS correct place
NSTATV = length(STATEV)
PROPS = {} # Not used in this example
NPROPS = {} # Not used in this example
COORDS = {} # Not used in this example
DROT = {} # Not used in this example
PNEWDT = {} # Not used in this example EXPLANATION MISSING
CELENT = {} # Not used in this example
DFGRD0 = {} # Not used in this example
DFGRD1 = {} # Not used in this example
NOEL = {} # Not used in this example
NPT = {} # Not used in this example
LAYER = {} # Not used in this example
KSPT = {} # Not used in this example
KSTEP = {} # Not used in this example
KINC = {} # Not used in this example









    Out[5]:





0-element Array{Any,1}

Finally the ccall of the umat



In [6]:

    
ccall((:umat_, "./libumat"), Int64, 
        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Int64},Ptr{Int64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
        STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,
        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,
        &NTENS,&NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
        NOEL,NPT,LAYER,KSPT,KSTEP,KINC)









    Out[6]:





0

Something happened



In [7]:

    
DDSDDE









    Out[7]:





4x4 Array{Float64,2}:
 282692.0  121154.0  121154.0      0.0
 121154.0  282692.0  121154.0      0.0
 121154.0  121154.0  282692.0      0.0
      0.0       0.0       0.0  80769.2



In [8]:

    
STATEV









    Out[8]:





1x2 Array{Float64,2}:
 0.0  0.0

Let's wrap this to simplified Julia function for testing



In [29]:

    
function isotropichardening!(stress,p,r,jacobian,strain,DSTRAN)
    local STRESS = stress
    local STATEV = [p r]
    local DDSDDE = jacobian
    local STRAN = strain
    o = ccall((:umat_, "./libumat"), Int64, 
        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Void},Ptr{Int64},Ptr{Int64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},
        Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
        STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,
        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,
        NSHR,&NTENS,&NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
        NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
    stress = STRESS
    p = STATEV[1]
    r = STATEV[2]
    jacobian = DDSDDE
    strain = STRAN + DSTRAN
    if o != 0 
        throw("UMAT failed. Return code is $o")
    end
    return o
end









    Out[29]:





isotropichardening! (generic function with 2 methods)

Now somebody should know how to use this function (help needed)



In [30]:

    
p = 0.0
r = 0.0
S = [0. 0. 0. 0.]
strain = [0. 0. 0. 0.]
for i in range(0,0.0001,11)
    DSTRAN = [i 0. 0. 0.]
    out = isotropichardening!(S,p,r,DDSDDE,strain,DSTRAN)
    println(out,S,p,r, DSTRAN)
    #println(i)
end









    



0[0.0 0.0 0.0 0.0]0.00.0[0.0 0.0 0.0 0.0]
0[28.269231713558856 12.11538570784259 12.11538570784259 0.0]0.00.0[0.0001 0.0 0.0 0.0]
0[84.80769514067657 36.34615712352777 36.34615712352777 0.0]0.00.0[0.0002 0.0 0.0 0.0]
0[169.61539028135314 72.69231424705553 72.69231424705553 0.0]0.00.0[0.00030000000000000003 0.0 0.0 0.0]
0[282.69231713558855 121.15385707842589 121.15385707842589 0.0]0.00.0[0.0004 0.0 0.0 0.0]






    



"UMAT failed. Return code is 62421072"
while loading In[30], in expression starting on line 5

 in isotropichardening! at In[29]:21
 in anonymous at no file:7

Let's find out how many different funtions & subroutines are called



In [11]:

    
fil = open("code_imp.f")
sub_list = Set{ASCIIString}()
for line in readlines(fil)
    # Fortran comment = something non whitespce at the firts character
    if ismatch(r"^[^\s]",line)
        #println(line)
        continue
    end
    if ismatch(r"call",lowercase(line))
        call = split(lowercase(line))[2] #divede by white space
        sub = split(call,"(")[1] #divide by "("
        #println(sub)
        push!(sub_list,sub)
    end
end
close(fil)



In [12]:

    
sub_list









    Out[12]:





Set{ASCIIString}({"keffp","kdevia","kmlt1"})

Next let's find out how many functions and subroutines are defined



In [13]:

    
fil = open("code_imp.f")
fun_list = Set{ASCIIString}()
for line in readlines(fil)
    if ismatch(r"^[^\s]",line)
        #println(line)
        continue
    end
    comp = lowercase(line)
    if ismatch(r"subroutine",comp) || ismatch(r"funtion",comp) || ismatch(r"external",comp)
        #println(line)
        call = split(lowercase(line),"(")[1] #divede by white space
        sub = split(call)[end] #divide by "("
        #if length(sub) > 1
        #    sub = sub[2]
        #end
        #println(sub)
        push!(fun_list,sub)
    end
end
close(fil)



In [14]:

    
fun_list









    Out[14]:





Set{ASCIIString}({"keffp","dyadicprod","umat","kdevia","dotprod","kmlt1"})

And Finally are all called subroutines defined



In [16]:

    
setdiff(sub_list,fun_list)









    Out[16]:





Set{ASCIIString}({})