First trials to use Abaqus UMAT interface

notebooks/2015-07-24-Abaqus-umat-interface.ipynb

@@ -0,0 +1,603 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Abaqus umat interface"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Author(s): Tero Frondelius"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Abstract:** making the initial version to call Abaqus umat"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "--2015-07-25 22:07:34--  http://www.eng.ox.ac.uk/NP/ICP/plasticity_imp/code_imp.f\n",
+      "Resolving www.eng.ox.ac.uk (www.eng.ox.ac.uk)... 163.1.223.199\n",
+      "Connecting to www.eng.ox.ac.uk (www.eng.ox.ac.uk)|163.1.223.199|:80... connected.\n",
+      "HTTP request sent, awaiting response... 200 OK\n",
+      "Length: 7458 (7,3K) [text/plain]\n",
+      "Saving to: ‘code_imp.f.1’\n",
+      "\n",
+      "     0K .......                                               100%  282M=0s\n",
+      "\n",
+      "2015-07-25 22:07:34 (282 MB/s) - ‘code_imp.f.1’ saved [7458/7458]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "run(`wget http://www.eng.ox.ac.uk/NP/ICP/plasticity_imp/code_imp.f`)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "***********************************************************************************\n",
+      "**  UMAT, FOR ABAQUS/STANDARD INCORPORATING ELASTIC-PLASTIC LINEAR               **\n",
+      "**  ISOTROPIC HARDENING. LARGE DEFORMATION FORMULATION FOR PLANE STRAIN          **\n",
+      "**  AND AXI-SYMMETRIC ELEMENTS. IMPLICIT INTEGRATION WITH CONSISTENT JACOBIAN    **\n",
+      "***********************************************************************************\n",
+      "***********************************************************************************\n",
+      "**\n",
+      "**\n",
+      "**\n",
+      "*USER SUBROUTINE\n",
+      "      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,\n",
+      "     1 RPL,DDSDDT,DRPLDE,DRPLDT,\n",
+      "     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,\n",
+      "     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,\n",
+      "     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)\n",
+      "C\n",
+      "      INCLUDE 'ABA_PARAM.INC'\n",
+      "C\n",
+      "      CHARACTER*80 CMNAME\n",
+      "C\n",
+      "C\n",
+      "      DIMENSION STRESS(NTENS),STATEV(NSTATV),\n",
+      "     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),\n",
+      "     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),\n",
+      "     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)\n",
+      "C\n",
+      "C\n",
+      "      PARAMETER (M=3,N=3,ID=3,ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,\n",
+      "     +          SIX=6.D0, NINE=9.D0, TOLER=1.D-5)\n",
+      "C\n"
+     ]
+    }
+   ],
+   "source": [
+    "run(`head -30 code_imp.f`)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## some implicit type castings"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "f = open(\"ABA_PARAM.INC\",\"w\")\n",
+    "write(f,\"      implicit real*8(a-h,o-z)\\n\")\n",
+    "write(f,\"      parameter (nprecd=2)\\n\")\n",
+    "close(f)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's compile the shared library"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "run(`gfortran -shared -fPIC -o libumat.so code_imp.f`)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Some Abaqus umat interface definitions\n",
+    "\n",
+    "|variable            | explanation                              |\n",
+    "|--------------------|------------------------------------------|\n",
+    "|DDSDDE(NTENS,NTENS) | Jacobian matrix of the constitutive model|\n",
+    "|STRESS(NTENS)       | the stress tensor (in vector format)     |\n",
+    "|STATEV(NSTATV)      | An array containing the solution-dependent state variables. |\n",
+    "|SSE                 | Specific elastic strain energy |\n",
+    "|SPD                 | plastic dissipation            |\n",
+    "|SCD                 | “creep” dissipation            |\n",
+    "|RPL                 | Volumetric heat generation per unit time |\n",
+    "|DDSDDT(NTENS)       | Variation of the stress increments with respect to the temperature. |\n",
+    "|DRPLDE(NTENS)       | Variation of RPL with respect to the strain increments.|\n",
+    "|DRPLDT              | Variation of RPL with respect to the temperature. |\n",
+    "|RPL          | RPL is used to indicate whether or not a cohesive element is open to the tangential flow of pore fluid.|\n",
+    "|PNEWDT              | Ratio of suggested new time increment to the time increment being used |\n",
+    "|STRAN(NTENS)        | An array containing the total strains at the beginning of the increment. |\n",
+    "|DSTRAN(NTENS)       | Array of strain increments. |\n",
+    "|TIME(1)             | Value of step time at the beginning of the current increment. |\n",
+    "|TIME(2)             | Value of total time at the beginning of the current increment. |\n",
+    "|DTIME               | Time increment.|\n",
+    "|TEMP                | Temperature at the start of the increment. |\n",
+    "|DTEMP               | Increment of temperature. |\n",
+    "|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.|\n",
+    "|DPRED               | Array of increments of predefined field variables. |\n",
+    "|CMNAME              | User-defined material name, left justified. |\n",
+    "|NDI                 | Number of direct stress components at this point. |\n",
+    "|NSHR                | Number of engineering shear stress components at this point. |\n",
+    "|NTENS               | Size of the stress or strain component array (NDI + NSHR). |\n",
+    "|NSTATV              | Number of solution-dependent state variables that are associated with this material type |\n",
+    "|PROPS(NPROPS)       | User-specified array of material constants associated with this user material. |\n",
+    "|NPROPS              | User-defined number of material constants associated with this user material. |\n",
+    "|COORDS              | An array containing the coordinates of this point.  |\n",
+    "|DROT(3,3)           | Rotation increment matrix. |\n",
+    "|CELENT              | Characteristic element length |\n",
+    "|DFGRD0(3,3)         | Array containing the deformation gradient at the beginning of the increment. |\n",
+    "|DFGRD1(3,3)         | Array containing the deformation gradient at the end of the increment. |\n",
+    "|NOEL                | Element number. |\n",
+    "|NPT                 | Integration point number. |\n",
+    "|LAYER               | Layer number (for composite shells and layered solids). |\n",
+    "|KSPT                | Section point number within the current layer. |\n",
+    "|KSTEP               | Step number. |\n",
+    "|KINC                | Increment number. |"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0-element Array{Any,1}"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "STRESS = [0. 0. 0. 0.]\n",
+    "p = 0. # EFFECTIVE PLASTIC STRAIN\n",
+    "r = 0. # ISOTROPIC HARDENING VARIABLE\n",
+    "STATEV = [p r]\n",
+    "NTENS = 4 \n",
+    "DDSDDE = zeros(NTENS,NTENS)\n",
+    "SSE = {} # Not used in this example\n",
+    "SPD = {} # Not used in this example\n",
+    "SCD = {} # Not used in this example\n",
+    "RPL = {} # Not used in this example\n",
+    "DDSDDT = {} # Not used in this example\n",
+    "DRPLDE = {} # Not used in this example\n",
+    "DRPLDT = {} # Not used in this example\n",
+    "STRAN = [0. 0. 0. 0.]\n",
+    "DSTRAN = [0. 0. 0. 0.]\n",
+    "TIME = [0. 0.1] # CHECK TIME(2)\n",
+    "DTIME = {} # Not used in this example\n",
+    "TEMP = {} # Not used in this example\n",
+    "DTEMP = {} # Not used in this example\n",
+    "PREDEF = {} # Not used in this example\n",
+    "DPRED = {} # Not used in this example\n",
+    "CMNAME = {} # Not used in this example CHARACTER*80 CMNAME\n",
+    "NDI = {} # Not used in this example\n",
+    "NSHR = {} # Not used in this example\n",
+    "#NTENS correct place\n",
+    "NSTATV = length(STATEV)\n",
+    "PROPS = {} # Not used in this example\n",
+    "NPROPS = {} # Not used in this example\n",
+    "COORDS = {} # Not used in this example\n",
+    "DROT = {} # Not used in this example\n",
+    "PNEWDT = {} # Not used in this example EXPLANATION MISSING\n",
+    "CELENT = {} # Not used in this example\n",
+    "DFGRD0 = {} # Not used in this example\n",
+    "DFGRD1 = {} # Not used in this example\n",
+    "NOEL = {} # Not used in this example\n",
+    "NPT = {} # Not used in this example\n",
+    "LAYER = {} # Not used in this example\n",
+    "KSPT = {} # Not used in this example\n",
+    "KSTEP = {} # Not used in this example\n",
+    "KINC = {} # Not used in this example"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Finally the ccall of the umat"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ccall((:umat_, \"./libumat\"), Int64, \n",
+    "        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Int64},Ptr{Int64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),\n",
+    "        STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,\n",
+    "        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,\n",
+    "        &NTENS,&NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,\n",
+    "        NOEL,NPT,LAYER,KSPT,KSTEP,KINC)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Something happened"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "4x4 Array{Float64,2}:\n",
+       " 282692.0  121154.0  121154.0      0.0\n",
+       " 121154.0  282692.0  121154.0      0.0\n",
+       " 121154.0  121154.0  282692.0      0.0\n",
+       "      0.0       0.0       0.0  80769.2"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "DDSDDE"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1x2 Array{Float64,2}:\n",
+       " 0.0  0.0"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "STATEV"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's wrap this to simplified Julia function for testing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "isotropichardening! (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 29,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "function isotropichardening!(stress,p,r,jacobian,strain,DSTRAN)\n",
+    "    local STRESS = stress\n",
+    "    local STATEV = [p r]\n",
+    "    local DDSDDE = jacobian\n",
+    "    local STRAN = strain\n",
+    "    o = ccall((:umat_, \"./libumat\"), Int64, \n",
+    "        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Float64},Ptr{Float64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Void},Ptr{Int64},Ptr{Int64},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},\n",
+    "        Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),\n",
+    "        STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,\n",
+    "        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,\n",
+    "        NSHR,&NTENS,&NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,\n",
+    "        NOEL,NPT,LAYER,KSPT,KSTEP,KINC)\n",
+    "    stress = STRESS\n",
+    "    p = STATEV[1]\n",
+    "    r = STATEV[2]\n",
+    "    jacobian = DDSDDE\n",
+    "    strain = STRAN + DSTRAN\n",
+    "    if o != 0 \n",
+    "        throw(\"UMAT failed. Return code is $o\")\n",
+    "    end\n",
+    "    return o\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Now somebody should know how to use this function (help needed)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {
+    "collapsed": false,
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0[0.0 0.0 0.0 0.0]0.00.0[0.0 0.0 0.0 0.0]\n",
+      "0[28.269231713558856 12.11538570784259 12.11538570784259 0.0]0.00.0[0.0001 0.0 0.0 0.0]\n",
+      "0[84.80769514067657 36.34615712352777 36.34615712352777 0.0]0.00.0[0.0002 0.0 0.0 0.0]\n",
+      "0[169.61539028135314 72.69231424705553 72.69231424705553 0.0]0.00.0[0.00030000000000000003 0.0 0.0 0.0]\n",
+      "0[282.69231713558855 121.15385707842589 121.15385707842589 0.0]0.00.0[0.0004 0.0 0.0 0.0]\n"
+     ]
+    },
+    {
+     "ename": "LoadError",
+     "evalue": "\"UMAT failed. Return code is 62421072\"\nwhile loading In[30], in expression starting on line 5",
+     "output_type": "error",
+     "traceback": [
+      "\"UMAT failed. Return code is 62421072\"\nwhile loading In[30], in expression starting on line 5",
+      "",
+      " in isotropichardening! at In[29]:21",
+      " in anonymous at no file:7"
+     ]
+    }
+   ],
+   "source": [
+    "p = 0.0\n",
+    "r = 0.0\n",
+    "S = [0. 0. 0. 0.]\n",
+    "strain = [0. 0. 0. 0.]\n",
+    "for i in range(0,0.0001,11)\n",
+    "    DSTRAN = [i 0. 0. 0.]\n",
+    "    out = isotropichardening!(S,p,r,DDSDDE,strain,DSTRAN)\n",
+    "    println(out,S,p,r, DSTRAN)\n",
+    "    #println(i)\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's find out how many different funtions & subroutines are called"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "fil = open(\"code_imp.f\")\n",
+    "sub_list = Set{ASCIIString}()\n",
+    "for line in readlines(fil)\n",
+    "    # Fortran comment = something non whitespce at the firts character\n",
+    "    if ismatch(r\"^[^\\s]\",line)\n",
+    "        #println(line)\n",
+    "        continue\n",
+    "    end\n",
+    "    if ismatch(r\"call\",lowercase(line))\n",
+    "        call = split(lowercase(line))[2] #divede by white space\n",
+    "        sub = split(call,\"(\")[1] #divide by \"(\"\n",
+    "        #println(sub)\n",
+    "        push!(sub_list,sub)\n",
+    "    end\n",
+    "end\n",
+    "close(fil)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Set{ASCIIString}({\"keffp\",\"kdevia\",\"kmlt1\"})"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sub_list"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Next let's find out how many functions and subroutines are defined"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "fil = open(\"code_imp.f\")\n",
+    "fun_list = Set{ASCIIString}()\n",
+    "for line in readlines(fil)\n",
+    "    if ismatch(r\"^[^\\s]\",line)\n",
+    "        #println(line)\n",
+    "        continue\n",
+    "    end\n",
+    "    comp = lowercase(line)\n",
+    "    if ismatch(r\"subroutine\",comp) || ismatch(r\"funtion\",comp) || ismatch(r\"external\",comp)\n",
+    "        #println(line)\n",
+    "        call = split(lowercase(line),\"(\")[1] #divede by white space\n",
+    "        sub = split(call)[end] #divide by \"(\"\n",
+    "        #if length(sub) > 1\n",
+    "        #    sub = sub[2]\n",
+    "        #end\n",
+    "        #println(sub)\n",
+    "        push!(fun_list,sub)\n",
+    "    end\n",
+    "end\n",
+    "close(fil)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Set{ASCIIString}({\"keffp\",\"dyadicprod\",\"umat\",\"kdevia\",\"dotprod\",\"kmlt1\"})"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fun_list"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## And Finally are all called subroutines defined"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Set{ASCIIString}({})"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "setdiff(sub_list,fun_list)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 0.3.8",
+   "language": "julia",
+   "name": "julia-0.3"
+  },
+  "language_info": {
+   "name": "julia",
+   "version": "0.3.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}