From bccae9a7d79b48974ad4ee8787f536bdf6655e6d Mon Sep 17 00:00:00 2001 From: Tero Frondelius Date: Sat, 25 Jul 2015 22:31:23 +0300 Subject: [PATCH] First trials to use Abaqus UMAT interface --- .../2015-07-24-Abaqus-umat-interface.ipynb | 603 ++++++++++++++++++ 1 file changed, 603 insertions(+) create mode 100644 notebooks/2015-07-24-Abaqus-umat-interface.ipynb diff --git a/notebooks/2015-07-24-Abaqus-umat-interface.ipynb b/notebooks/2015-07-24-Abaqus-umat-interface.ipynb new file mode 100644 index 0000000..5395eda --- /dev/null +++ b/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 +}