JuliaFEM.jl/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constitutive modelling using JuliaFEM\n",
    "\n",
    "Author(s): Jukka Aho\n",
    "\n",
    "**Abstract**: This is reproduction of the results from notebook [*Ideal plastic Von Mises material*](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-09-24-Ideal%20plastic%20Von%20Mises%20material.ipynb) made by Olli Väinölä. Small strain theory is used, see https://en.wikipedia.org/wiki/Flow_plasticity_theory. The purpose of this notebook is to show how one can easily design and simulate material models using JuliaFEM. This is 2d version. For 3d version see Olli's notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuliaFEM: IntegrationPoint, Field, FieldSet, TimeStep, Increment\n",
    "using ForwardDiff\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The atomic structure here is `IntegrationPoint`. It has identical `Field`-structure like elements and can store multidimensional variables which can be time-dependent also. That way one can easily store, for example, measured strain and run material simulation for real measured data and fit material parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaFEM.Increment[[0.0]])])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ip = IntegrationPoint([0.0, 0.0], 1.0)\n",
    "\n",
    "ip.fields[\"total strain\"] = Field()\n",
    "ip.fields[\"plastic strain\"] = Field()\n",
    "ip.fields[\"plastic potential\"] = Field()\n",
    "ip.fields[\"elastic strain\"] = Field()\n",
    "ip.fields[\"effective plastic strain\"] = Field()\n",
    "ip.fields[\"stress\"] = Field()\n",
    "ip.fields[\"young\"] = Field(200.0e9)\n",
    "ip.fields[\"poisson\"] = Field(0.3)\n",
    "ip.fields[\"yield stress\"] = Field(200.0e6)\n",
    "ip.fields[\"plastic rate parameter\"] = Field(0.0) # material parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accessing parameters is done just like with `FieldSet`s in general:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: [2.0e11]\n",
      "INFO: [2.0e8]\n"
     ]
    }
   ],
   "source": [
    "info(last(ip.fields[\"young\"]))\n",
    "info(last(ip.fields[\"yield stress\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next to material model: ideal plastic material model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "calculate_stress! (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\" Ideal plastic material model. \"\"\"\n",
    "function calculate_stress!(ip::IntegrationPoint, strain::Matrix, time::Number)\n",
    "\n",
    "    dim = size(strain, 1)\n",
    "    \n",
    "    poisson = last(ip.fields[\"poisson\"])[1]\n",
    "    young = last(ip.fields[\"young\"])[1]\n",
    "    mu = young/(2*(1+poisson))\n",
    "    lambda = young*poisson/((1+poisson)*(1-2*poisson))\n",
    "    if dim == 2\n",
    "        lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d\n",
    "    end\n",
    "\n",
    "    # yield function\n",
    "    function f(stress)\n",
    "        # https://en.wikipedia.org/wiki/Yield_surface\n",
    "        # https://en.wikipedia.org/wiki/Von_Mises_yield_criterion\n",
    "        stress_y = last(ip.fields[\"yield stress\"])[1]\n",
    "        s1 = stress[1,1]\n",
    "        s2 = stress[2,2]\n",
    "        s12 = stress[1,2]\n",
    "        stress_v = sqrt(s1^2 - s1*s2 + s2^2 + 3*s12^2)\n",
    "        return stress_v - stress_y\n",
    "    end\n",
    "\n",
    "    strain_plastic_prev = last(ip.fields[\"plastic strain\"])[1]\n",
    "    strain_elastic = strain - strain_plastic_prev\n",
    "    stress_trial = lambda*trace(strain_elastic)*I + 2*mu*(strain_elastic)\n",
    "\n",
    "    if f(stress_trial) <= 0.0\n",
    "        #info(\"time=$time: no yield\")\n",
    "        push!(ip.fields[\"stress\"], TimeStep(time, Increment(Matrix[stress_trial])))\n",
    "        push!(ip.fields[\"total strain\"], TimeStep(time, Increment(Matrix[strain])))\n",
    "        return true\n",
    "    else\n",
    "        #info(\"time=$time: yield, f(stress_trial) = $(f(stress_trial))\")\n",
    "        dt = time - ip.fields[\"total strain\"][end].time\n",
    "\n",
    "        # associated flow rule, plastic potential ψ(σ) = f\n",
    "        psi = f\n",
    "\n",
    "        \"\"\" Calculate equations\n",
    "\n",
    "            dσ - C (dϵₜ - dγ*dΨ/dσ) = 0\n",
    "                         σₑ(σ) - σy = 0\n",
    "\n",
    "            dϵₜ = total strain\n",
    "            dϵₑ = elastic strain\n",
    "            dϵₚ = plastic strain\n",
    "        \n",
    "            dϵₜ = dϵₑ + dϵₚ\n",
    "            => dϵₑ = dϵₜ - dϵₚ = dϵₜ - dγ*dψ/dσ\n",
    "            dσ = λ⋅tr(dϵₑ)I + 2μ⋅dϵₑ\n",
    "            => dσ - λ⋅tr(dϵₑ)I - 2μ⋅dϵₑ = 0\n",
    "        \"\"\"\n",
    "        function residual(params::Vector)\n",
    "            dstress = reshape(params[1:prod(size(strain))], size(strain))\n",
    "            gamma = params[end]\n",
    "            strain_prev = last(ip.fields[\"total strain\"])[1]\n",
    "            stress_prev = last(ip.fields[\"stress\"])[1]\n",
    "            dstrain_total = 1/dt*(strain - strain_prev)\n",
    "            stress_tot = stress_prev + dstress\n",
    "            # derivative of plastic potential ψ(σ) w.r.t 2nd order tensor σ(ϵ)\n",
    "            # https://en.wikipedia.org/wiki/Tensor_derivative_%28continuum_mechanics%29\n",
    "            # \"Derivatives of scalar valued functions of second-order tensors\"\n",
    "            dpsi_dstress = derivative(psi, stress_tot)\n",
    "            dstrain_plastic = gamma*dpsi_dstress\n",
    "            dstrain_elastic = dstrain_total - dstrain_plastic\n",
    "            dstress_elastic = lambda*trace(dstrain_elastic)*I + 2*mu*dstrain_elastic\n",
    "            stress_delta = dstress - dstress_elastic\n",
    "            return [vec(stress_delta); psi(stress_tot)]\n",
    "        end\n",
    "\n",
    "        # solve equations using Newton iterations. Jacobian is calcualated\n",
    "        # using automatic differentiation as usual.\n",
    "        stress_prev = last(ip.fields[\"stress\"])[1]\n",
    "        initial_gamma = last(ip.fields[\"plastic rate parameter\"])[1]\n",
    "        params = [vec(stress_prev); initial_gamma]\n",
    "        dparams = zeros(5)\n",
    "        l = [1, 4, 3, 5] # <-- reorder to voigt\n",
    "        for iterations=1:10\n",
    "            A = ForwardDiff.jacobian(residual, params)[l,l]\n",
    "            b = -residual(params)[l]\n",
    "            dparams = A \\ b\n",
    "            params[l] += dparams\n",
    "            norm(dparams) < 1.0e-7 && break\n",
    "        end\n",
    "        params[2] = params[3]\n",
    "\n",
    "        # save all kind of stuff to integration point\n",
    "\n",
    "        dstress = reshape(params[1:prod(size(strain))], size(strain))\n",
    "        stress = stress_prev + dstress\n",
    "        dpsi_dstress = derivative(psi, stress)\n",
    "        gamma = params[end]\n",
    "        dstrain_plastic = gamma*dpsi_dstress\n",
    "        strain_plastic_prev = last(ip.fields[\"plastic strain\"])[1]\n",
    "        strain_plastic  = strain_plastic_prev + dstrain_plastic\n",
    "\n",
    "        push!(ip.fields[\"stress\"], TimeStep(time, Increment(Matrix[stress])))\n",
    "        push!(ip.fields[\"total strain\"], TimeStep(time, Increment(Matrix[strain])))\n",
    "        push!(ip.fields[\"elastic strain\"], TimeStep(time, Increment(Matrix[strain_elastic])))\n",
    "        push!(ip.fields[\"plastic strain\"], TimeStep(time, Increment(Matrix[strain_plastic])))\n",
    "        push!(ip.fields[\"plastic rate parameter\"], TimeStep(time, Increment(params[end])))\n",
    "        push!(ip.fields[\"plastic potential\"], TimeStep(time, Increment(psi(stress))))\n",
    "        push!(ip.fields[\"derivative of plastic potential\"], TimeStep(time, Increment(Matrix[dpsi_dstress])))\n",
    "    end\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small \"simulator\" to study the behavior of material model in `IntegrationPoint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elapsed time: 0.031678742 seconds\n"
     ]
    }
   ],
   "source": [
    "function run(steps=11)\n",
    "\n",
    "    # initialization\n",
    "    ip = IntegrationPoint([0.0, 0.0], 1.0)\n",
    "    ip.fields[\"total strain\"] = Field()\n",
    "    ip.fields[\"plastic strain\"] = Field(Matrix[zeros(2,2)])\n",
    "    ip.fields[\"plastic potential\"] = Field()\n",
    "    ip.fields[\"derivative of plastic potential\"] = Field()\n",
    "    ip.fields[\"elastic strain\"] = Field()\n",
    "    ip.fields[\"effective plastic strain\"] = Field()\n",
    "    ip.fields[\"stress\"] = Field()\n",
    "    ip.fields[\"young\"] = Field(200.0e9)\n",
    "    ip.fields[\"poisson\"] = Field(0.3)\n",
    "    ip.fields[\"yield stress\"] = Field(200.0e6)\n",
    "    ip.fields[\"plastic rate parameter\"] = Field(0.0) # material parameter\n",
    "\n",
    "    strain = zeros(2, 2)\n",
    "    strain[1,1] = 2.0e-3\n",
    "    strain[2,2] = -0.3*2.0e-3\n",
    "    #strain[1,2] = 1.0e-3\n",
    "    #strain[2,1] = 1.0e-3\n",
    "\n",
    "    # calculate stress in integration point\n",
    "    #for (time, omega) in enumerate(linspace(0, 4*pi, steps))\n",
    "    #    calculate_stress!(ip, sin(omega)*strain, Float64(time))\n",
    "    #end\n",
    "\n",
    "    for (time, k) in enumerate(linspace(0, 1, steps))\n",
    "        calculate_stress!(ip, k*strain, Float64(time))\n",
    "    end\n",
    "\n",
    "    return ip\n",
    "end\n",
    "\n",
    "tic()\n",
    "ip = run()\n",
    "toc();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAHqCAYAAACQvy+/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8TXf+x/HXTSKJNZIIYo1BEV3s69Rata+1r6FGi6JoY9rOTBeqGrSoMrT2pbWlxBpKSkoq1DL9Waat0ao9hNgSkeT8/rhy68oiy01ulvfz8biPuGf7fs7tbbyd8/1+j8kwDAMRERERyXcc7F2AiIiIiNiHgqCIiIhIPqUgKCIiIpJPKQiKiIiI5FMKgiIiIiL5lIKgiIiISD6lICgiIiKSTykIioiIiORTCoIiIiIi+ZSCoIhIDuPj40OlSpXsXYaI5AMKgiJ5WHx8PF988QXNmzfHw8MDZ2dnSpUqxXPPPcff/vY3Nm/ebLX90qVLcXBwYNmyZXaqOOvFxMQwY8YMGjZsiJubGy4uLpQpU4Z69eoxZswY9u3bZ7X9e++9h4ODQ5LlWclkMmEymbKtPRHJv5zsXYCIZI34+Hg6depEcHAw7u7udOrUiXLlyhEbG8v//d//sXr1av773//SuXPnJPvm1RBy584dmjdvztGjR/H29qZXr16ULl2aO3fucOzYMRYuXEhUVBTNmjWza5179uyxa/sikn8oCIrkUV999RXBwcHUqlWLvXv3UrRoUav10dHRhIeHJ7uvYRjZUWK2mzVrFkePHqVt27Zs3rwZJyfrX4E3b97k9OnTye6bnZ+JbguLSHbRrWGRPOrAgQMA+Pn5JQmBAAULFqR58+aW9y1atGDYsGEADB06FAcHB8vr3LlzwJ+3Sffu3cvq1atp2LAhRYoUsQou9+7d46OPPqJWrVoUKVKEokWL0qRJE77++utk61y2bBlNmjTBy8uLggULUqFCBdq1a8fatWuttvvPf/5Dv3798PHxwdXVlZIlS1K3bl3Gjx9PXFxcuj6TkSNHJgmBAMWLF6dRo0aW9z4+PnzwwQcAtGzZ0uozSeTn54eDgwNnz57ls88+49lnn6VQoUK0bNkSgAcPHjB37lw6dOhAxYoVcXV1xdPTkzZt2rBjx45k60yuj+Cjt+1DQkJo0aIFxYoVw83NjU6dOqUYYFOzc+dOOnfuTMmSJXF1daVChQp069aN3bt3J9tuchwcHCznmii178nBgwdxcHCgR48eKdZVo0YNXF1duXnzptXy4OBgOnToQIkSJXB1daVKlSr4+/sTFRWV7nMXETNdERTJo0qUKAHAf//73zRtP3ToUNzd3dm0aRPdunWjVq1alnVubm5W286cOZNdu3bRpUsXWrdubfmL+ObNm7Rq1Ypjx45Rt25dXn75ZRISEtixYwf9+/fnxIkTTJ482XKct99+m2nTpvGXv/yFvn374ubmxsWLFzl06BDr16+nd+/egDkENmzYEEdHR7p06UKlSpW4desWv/zyC/Pnz+fDDz9MNthl9jMZP348GzduZO/evfj5+eHj45PituPGjSM0NJROnTrRqVMnHB0dAbh+/Tqvv/46TZs2pW3btnh5eXHx4kU2b95Mhw4d+OKLL3j55ZeTHC+l2/Nbtmxh06ZNdOjQgZEjR3LixAm2bdvGoUOHOHnyJJ6enmk6t3fffZfJkydTtGhRunXrRvny5blw4QIHDhxg1apVtG7dOk31pLYuue9Jw4YNqVatGtu2bSMyMhIPDw+rfcLDw/nvf/9Lz549KV68uGX5+++/z/vvv4+np6clvB4/fpwZM2awbds2wsLCkv0Hj4g8gSEiedLRo0cNZ2dnw8HBwRg0aJARGBho/Pbbb6nus2TJEsNkMhnLli1Ldv27775rmEwmo0iRIsaxY8eSrB8yZIhhMpmM6dOnWy2PiYkx2rVrZzg4OFjt5+HhYZQvX96Ijo5Ocqxr165Z/jxhwgTDZDIZQUFBSba7efOmkZCQkOp5JdqyZYthMpkMFxcXY9SoUcbWrVuNixcvprpP4jnv3bs32fWJ51yuXLlkP9/79+8bFy5cSLI8KirKePrppw0PD48k51+xYkWjUqVKVssS/9sUKFDA2LNnj9W6t956yzCZTEZAQECq55IoODjYMJlMRuXKlZM9//PnzydpN6XvhMlkMlq2bGm17Enfk48++sgwmUzG3Llzk6wbNWqUYTKZjC1btliW7dmzxzCZTEbTpk2NqKgoq+2XLl1qmEwmY/z48amftIgkS7eGRfKoWrVqsXLlSkqVKsXKlSt56aWXqFSpEp6envTo0YMtW7Zk+NgjRozgueees1p2/fp1Vq5cSf369XnjjTes1rm4uDBt2jQMw2D16tWW5SaTiQIFCljdak2U3JUtV1fXJMvc3NzSPLilY8eOzJ49m4IFCzJ//nw6depE2bJl8fb2ZuDAgYSGhqbpOMnx9/enYsWKSZY7OztTpkyZJMuLFSvG0KFDuXHjBocOHUpzO3379k1yK3bEiBEAaT7OZ599Bpiv2Hl7eydZX7Zs2TTXk5rkvicAgwYNSvZ2c2xsLF9//TWlSpWiffv2luVz5swB4IsvvqBYsWJW+wwZMoTnnnuOVatW2aRmkfxGt4ZF8rBevXrRvXt3QkJC2L9/P0ePHuX7779n48aNbNy4kcGDB7N06dJ0H7dBgwZJlh06dIiEhATA3EfscQ8ePADg1KlTlmUDBgzgs88+w9fXl969e9O8eXMaNWqU5FZ03759mTNnDt26daNnz560bt2apk2bUrlyZavtjh07xsaNG62Wubu7M27cOMv7MWPGMHz4cHbt2kVYWBhHjx7lwIEDrF69mtWrV/PPf/6T999/P30fCMl/JolOnDjB9OnT2bdvH5cvXyYmJsZq/cWLF9PcTr169ZIsK1euHAA3btxI0zF++OEHHBwcaNeuXZrbzYiUPpOyZcvSunVrdu3axalTp6hRowYAmzdv5saNG0yYMMHqHwdhYWEUKFCAtWvXJjtoJzY2loiICG7cuIG7u3vWnIxIHqUgKJLHOTk50aZNG9q0aQNAQkICGzZsYNiwYSxfvpzu3bvTtWvXdB2zdOnSSZZdv34dMAfClK5MmUwm7t69a3n/6aef8pe//IUlS5Ywbdo0pk2bhpOTEx06dGDmzJmWoFe/fn1CQ0P58MMPWb9+PStWrACgWrVqvPvuu/Tt2xeA48ePWwZ3JPLx8bEKgmAeKNOlSxe6dOkCmEPqF198wbhx45g8eTI9evRI9kpWej8TMIeuVq1akZCQQOvWrenWrRvFihXDwcGBo0ePsmnTJu7fv5/mdh7tN5cosX9kfHx8mo5x8+ZN3N3dcXFxSXO7GZHSZwLmQTa7du1i2bJlTJs2DcByhXDIkCFW216/fp34+PhUA7rJZOLOnTsKgiLppFvDIvmMg4MDvXr1Yvz48QCEhISk+xjJ3YpNvIo3YcIEEhISkn3Fx8dbjUh1cHBg3LhxHDt2jCtXrrBhwwa6d+9OUFAQ7dq1IzY21rJto0aN2Lx5Mzdv3mT//v3885//5MqVK/Tv399yzCFDhiRp83//+98Tz6dAgQKMGjWKfv36ARmbxy+l29NTpkwhJiaGnTt3snXrVj755BPee+89/vWvf6V6FTErFS9enBs3biS5MpmcxCtzyY3MfnxU7+NSu2XfvXt3ihUrxsqVKzEMg6tXr7J9+3Zq1arFM888Y7Wtm5sbHh4eKX6vEr9b5cuXf+L5iIg1BUGRfKpIkSKA9fx4iSNd03pl6VENGzbM1BM4vLy86N69O2vWrKFly5acOXOGEydOJNmuQIECNG7cmPfff9/SdywoKChDbT4u8TN5VGY+E4Bff/0VT0/PZCep3rt3b4aOmVmNGze2jOZ+ksQrbIlTCD3q8OHDGa7B1dWV3r17c/HiRXbt2sXq1auJj49PcjUwsd7IyEhOnjyZ4fZEJHkKgiJ51FdffcW3336bbJ+qy5cv88UXXwBYBZTEARq///57utvz8vJiwIABHD58mClTplj6Cz7qzJkz/Pbbb4C5X9f+/fuTbPPgwQMiIyMxmUwUKlQIMM//l9zVq8uXLwNYtnuSf//73xw8eDDZdadPn2bdunWYTCabfSZgnhz6+vXr/PT
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7fa8b80dd450>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.legend.Legend object at 0x7fa8a16155d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "steps = length(ip.fields[\"total strain\"])\n",
    "eps11 = Float64[]\n",
    "sig11 = Float64[]\n",
    "sig22 = Float64[]\n",
    "sig12 = Float64[]\n",
    "principals = Vector{Float64}[]\n",
    "for i=1:steps\n",
    "    # extract from integration points\n",
    "    # field -> timestep -> increment -> (vector of tensors, take first) -> (first component)\n",
    "    strain = ip.fields[\"total strain\"][i][end][1]*1.0e6\n",
    "    stress = ip.fields[\"stress\"][i][end][1]*1.0e-6\n",
    "    push!(eps11, strain[1,1])\n",
    "    push!(sig11, stress[1,1])\n",
    "    push!(sig22, stress[2,2])\n",
    "    push!(sig12, stress[1,2])\n",
    "    push!(principals, eigvals(stress))\n",
    "end\n",
    "\n",
    "PyPlot.figure(figsize=(7, 5))\n",
    "PyPlot.plot(eps11, sig11, \"-bo\", label=\"s11\")\n",
    "PyPlot.plot(eps11, sig22, \"-ro\", label=\"s22\")\n",
    "PyPlot.plot(eps11, sig12, \"-go\", label=\"s12\")\n",
    "PyPlot.title(\"Stress-Strain curve\")\n",
    "PyPlot.xlabel(\"Strain [ustr]\")\n",
    "PyPlot.ylabel(\"Stress [MPa]\")\n",
    "#PyPlot.ylim([-250, 250])\n",
    "#PyPlot.xlim([-2100, 2100])\n",
    "PyPlot.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAI6CAYAAADMlPm5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcTekfB/DPvXVblKiESHYlyyRZasZSRJYsg2RsZZuxMwYzZcY2YwnDGDOGkW1sGUIYQpYak3ZrsiRLiiI7abnn9wfdn1S0n+69n/frdV+v6TzPPedzLtP9ep7nnCMRBEEAERERkYqTih2AiIiIqCyw6CEiIiK1wKKHiIiI1AKLHiIiIlILLHqIiIhILbDoISIiIrXAooeIiIjUAoseIiIiUgsseoiIiEgtsOghojxJpVI4ODiU6TE7duwIqbTsfy3dvHkTUqkUHh4eZX5sIio7LHqIKF8SiaTMj1fWx3z/+CWtTp06qFu3bonvl4gKT1PsAERUPsXGxqJChQpleszNmzfj1atXZXrMsiBmIUdE/8eih4jy1KhRozI/Zq1atcr8mESkPji9RaTk3l2PEhsbiz59+sDIyAj6+vpo164djh49mus9GzduhFQqxaZNm3D48GF07NgRlSpVyrGeJq81PXPmzIFUKsWpU6ewa9cutG7dGnp6ejA2NsagQYOQmJiYZ8bU1FR4eXmhadOm0NPTQ+XKlWFtbY3vvvsOL1++VPTLa03PyZMnIZVKMXfuXISEhKBz586oXLkyDAwM4OzsjMjIyFzHS0xMxLx58/Dpp5+ievXq0NbWRs2aNTF48GBcvny5UJ9vfjZt2gR7e3uYmJhAV1cX5ubmcHZ2xs6dO3Pkvn37tuLPKPv17tqh7M/5/v37GDVqFGrWrAlNTU1s2rRJ0Sc0NBT9+/dXnIu5uTm++uorJCUl5cp148YNjBkzBg0aNECFChVgbGyM5s2bY+zYsUhNTVX0S09Px8qVK2FjYwMjIyPo6emhbt266NOnDwIDA0vkMyIqbzjSQ6Qi4uPjYW9vr/iCS0xMhK+vL7p164Zt27bB1dU113t27dqFw4cPo3v37hg3bhxu3bqVoz2/aZnff/8d/v7+6N27NxwcHHDmzBn4+vri3LlzOHv2LLS0tHLkcnBwwO3bt2Fra4tx48ZBLpfjypUrWLFiBcaOHQtzc/OPHjM0NBQLFy6Ek5MTJkyYgGvXrsHPzw9BQUE4cuQIPvvsM0XfoKAgLF68GI6OjrCxsYG+vj6uXr2KXbt2wd/fH6dPn0bz5s0L9fm+y9PTE4sWLUK9evXg5uaGSpUqITExEeHh4di1axdcXV1Rt25dzJ49GytWrAAATJ06VfF+a2vrHPtLTU1F27ZtUbFiRfTv3x9SqRTVq1cHAKxfvx5jxoyBrq4uevXqhVq1auHq1atYt24d9u/fjzNnzihGyJKSktCqVSs8e/YMPXr0wIABA5CWloYbN25gy5YtmDhxIoyMjAAA7u7u2LFjB5o1a4bhw4dDV1cXd+/exenTpxEQEIBOnToV+fMhKrcEIlJq8fHxgkQiESQSiTBjxowcbREREYJMJhMMDQ2Fp0+fKrZv2LBBkEgkgoaGhhAQEJDnfiUSieDg4JBj2+zZswWJRCJUqlRJuHjxYo62L774QpBIJMLOnTtzbLezsxMkEomwaNGiXMd4+PChkJaWpvi5Q4cOglQqzdHnxIkTivP77bffcrTt27dPkEgkQsOGDQW5XK7YnpycLDx//jzX8c6dOyfo6+sL3bp1y7E9+zP08PDI66PIxcjISKhVq5bw6tWrXG0PHjzI8XPt2rWFunXr5ruv7HMbPny4kJWVlaPtypUrgkwmExo2bCgkJibmaAsMDBQ0NDSEvn37KratXLlSkEgkwsqVK3Md5+XLl4q8jx8/FiQSidCqVascn1u2hw8f5puXSJlxeotIRVSuXBk//PBDjm0tW7bE4MGD8fjxY+zZsyfXe3r37o0uXboU+liTJk1CkyZNcmwbPXo0ACA8PFyxLTIyEmfOnEGLFi0wc+bMXPsxMjKCtrZ2gY7ZsGFDjBs3Lse2Xr16oUOHDrh+/TqCg4MV201MTKCnp5drH82bN4eDgwNOnDiBrKysAh03LxKJBDKZLM/L642NjQu9P21tbSxdujTX/lavXo3MzEz88ssvMDU1zdHm6OgIFxcX7N+/Hy9evMjRpqOjk+sYurq6iu3Zo2na2tp5jqxljwYRqRoWPUQqwsbGJs8v+g4dOgAAzp49m6utdevWRTqWra1trm1mZmYAgEePHim2nTlzBgDQtWvXIh3nXe3atctze37nd/DgQbi4uMDU1BRaWlqK9TQHDhxAeno6Hjx4UOQsgwcPRnx8PKysrODp6YmAgAA8efKkyPurU6cOqlSpkmt7SEgIgDfrg+bMmZPrlZycjKysLFy5cgXAmyJWX18f48ePR//+/bF27VrExMTk2q+BgQFcXFxw+vRpWFtbY/78+Th58mSO9VVEqohreohURLVq1fLcnr02JK8v5ey2wqpcuXKubZqab36dvDuC8vjxYwBAzZo1i3ScdxXm/H755RdMnToVRkZGcHJygrm5OSpUqACJRII9e/bg3LlzeP36dZGzLF++HPXq1cOGDRuwaNEiLFq0CJqamujevTuWLVuG+vXrF2p/+f05PHz4EACwZMmSfN8rkUgUIz3m5uYICwvDnDlzcPjwYfj5+QF4c1XcN998g4kTJyre5+vri8WLF2Pbtm2YPXs2gDcjRP3798fSpUtRtWrVQp0DkTJg0UOkIu7fv5/n9nv37gEAKlWqlKuttO8fk10cJSQkFHtfBT2/zMxMzJkzB6ampoiKispVLJ0+fbrYWaRSKSZPnozJkycjJSUF//77L3bs2IG///4bly5dwqVLl3Is5v6Y/P4cKlWqBIlEgidPnkBfX79A+7K0tMSOHTuQlZWFc+fO4dixY/j1118xefJk6OnpYcSIEQDeFDizZ8/G7NmzkZCQgKCgIGzcuBFbtmzBzZs3ERQUVOD8RMqC01tEKiIqKgrPnz/Ptf3kyZMAgBYtWpRxIsDOzg4AEBAQAEEQirWv4ODgPPfx/vk9ePAAT548gb29fa6C5/nz54iKiirRYs/ExAR9+/aFr68vHBwcEBcXh0uXLinaNTQ0irx+yM7ODoIgFKkA0dDQgI2NDWbMmIHt27cDAPbt25dnXzMzM3zxxRcICAhA/fr18e+//+aYpiRSFSx6iFTE48ePMW/evBzbIiIisHXrVlSuXBl9+/Yt80w2Njawt7fH2bNnsXjx4lztDx8+LPA007Vr1/D777/n2LZv3z4EBQWhYcOGijU/VatWRYUKFRAREZFjgW9GRgYmT56smDIqqvT09DxHizIyMpCamgqJRJLjTtbGxsZITk5GWlpaoY81YcIEyGQyTJ06FdeuXcszy7sLuKOiovKcxsweDcvO9eDBA1y4cCFXv+fPn+P58+eQyWSFGqkiUhac3iJSEe3bt8e6desQGhoKe3t7JCUlwdfXFwCwZs2aAk+PlLQtW7agY8eO8PT0xO7du9GhQwcIgoBr167h6NGjuHLlSo779OQ3IuTs7Ixp06bh0KFDaN68Oa5fvw4/Pz/o6upi/fr1in5SqRSTJk3CokWL0KxZM/Tq1Qvp6ek4ceIEHj9+rLh6q6hevnyJdu3aoUGDBrCxsUHt2rWRlpaGo0ePIjY2Fr1794aFhYWif+fOnREREQFnZ2e0a9cO2trasLa2Rs+ePT96LAsLC6xfvx4jRoxAkyZN4OzsjIYNGyIjIwO3b99GcHAwqlWrplisvHnzZqxduxafffYZ6tWrB0NDQ8TFxWH//v3Q0dHBlClTALyZbrSxsUGzZs3QrFkz1KpVC0+fPsWBAwdw//59xVQYkcoR9YJ5Iiq2d+8xExsbK/Tu3VswNDQU9PT0hM8++0w4cuRIrvds3LhRkEqlwqZNm/Ldb1736ZkzZ44glUqFU6dOfTDH+x4+fCjMnDlTsLCwEHR0dARDQ0OhRYsWwqxZs4SXL18q+nXs2DHf+/TMnTtXCAkJETp37iwYGBgIBgYGQteuXYWIiIhcx8vMzBR+/vlnwcrKStDV1RVMTU2FYcOGCbdv3xbc3d0FqVQq3Lp1q0D
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7fa8a1931150>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{Any,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x7fa8a1513050>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function plot_principal()\n",
    "    PyPlot.figure(figsize=(6, 6))\n",
    "    n = 100\n",
    "    s = linspace(-300, 300, n)\n",
    "    s1 = repmat(s', n, 1)\n",
    "    s2 = repmat(s, 1, n)\n",
    "    sigma = sqrt(s1.^2 + s2.^2 - s1.*s2)\n",
    "    contour(s1, s2, sigma, [200], colors=\"k\")\n",
    "    p1 = [p[1] for p in principals]\n",
    "    p2 = [p[2] for p in principals]\n",
    "    PyPlot.plot(p1, p2, \"-bo\")\n",
    "    axis(\"equal\")\n",
    "    xlabel(\"sigma 1\")\n",
    "    ylabel(\"sigma 2\")\n",
    "    title(\"principal stress\")\n",
    "    dir = last(ip.fields[\"derivative of plastic potential\"])[1]\n",
    "    PyPlot.plot([p1[end], p1[end]+dir[2,2]*50],\n",
    "                [p2[end], p2[end]+dir[1,1]*50], \"-r\")\n",
    "end\n",
    "plot_principal()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.0",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}