{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8.3 Damage function within a fininite element solver\n",
    "@todo - this notebook must be updated - currently not running"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path = ['/home/rch/git/bmcs'] + sys.path\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Define the damage functions\n",
    "\n",
    "\\begin{align}\n",
    "1 - \\frac{f_t}{E \\varepsilon} \\exp(-\\frac{f_t}{G_f} (\\varepsilon - \\varepsilon_0) L_s )\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from ibvpy.mats.mats_damage_fn import \\\n",
    "    IDamageFn, LiDamageFn, JirasekDamageFn, AbaqusDamageFn,\\\n",
    "    PlottableFn, DamageFn, GfDamageFn\n",
    "from traits.api import Float, Property\n",
    "from traitsui.api import View, Item, VGroup, UItem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "f_t = 2.4\n",
    "G_f = 0.090\n",
    "E = 30000.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "omega_fn = GfDamageFn(G_f=G_f, f_t=f_t, L_s=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "eps = np.array([0.00013033], dtype=np.float_)\n",
    "omega_fn.trait_set(L_s = 10)\n",
    "omega = omega_fn(eps)\n",
    "sig = (1 - omega) * E * eps\n",
    "eps, omega, sig, omega_fn.f_t, omega_fn.eps_0, omega_fn.L_s, omega_fn.G_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ax = plt.subplot(111)\n",
    "omega_fn.plot(ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Nonlinear solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "eps_0 = f_t / E\n",
    "eps_ch = G_f / f_t\n",
    "L_ch = E * G_f / f_t**2\n",
    "\n",
    "print('eps_0', eps_0)\n",
    "print( 'L_ch', L_ch)\n",
    "\n",
    "L = 200.0\n",
    "u_max = 0.2\n",
    "eps_max = eps_ch\n",
    "eps = np.linspace(0, eps_max, 4)\n",
    "\n",
    "omega_fn_abaqus = AbaqusDamageFn(s_0=eps_0, s_u=0.03)\n",
    "\n",
    "n_T = 50\n",
    "K_max = 200\n",
    "\n",
    "for N in [2,4,6]:\n",
    "    sig_t = []\n",
    "    eps_t = []\n",
    "    L_el = (N - 1.0) / N * L\n",
    "    L_s = 10#1.0 / N * L\n",
    "    eps_s = 0.0\n",
    "    eps_s_arr = np.array([eps_s], dtype=np.float_)\n",
    "    u_t = np.linspace(0, u_max, n_T)\n",
    "    #omega_fn.L_s = L_s\n",
    "    print('elem size', L_s)\n",
    "    for u in u_t:\n",
    "        print('=== increment %g ' % u,)\n",
    "        for K in range(K_max):\n",
    "            omega = omega_fn(eps_s_arr)\n",
    "            eps_s = eps_s_arr[0]\n",
    "            u_s = eps_s * L_s\n",
    "            u_el = u - u_s\n",
    "            R = 1 / L_el * (u - eps_s_arr * L_s) - \\\n",
    "                (1 - omega) * eps_s_arr\n",
    "            if np.fabs(R) < 1e-8:\n",
    "                print('converged in %d iterations' % K)\n",
    "                break\n",
    "            dR = -L_s / L_el + omega_fn.diff(eps_s_arr) * \\\n",
    "                eps_s_arr - (1 - omega_fn(eps_s_arr))\n",
    "            d_eps_s = -R / dR\n",
    "            eps_s_arr += d_eps_s\n",
    "            if K == K_max - 1:\n",
    "                raise ValueError('No convergence')\n",
    "        sig = ((1.0 - omega_fn(eps_s_arr)) * E * eps_s_arr)[0]\n",
    "        sig_t.append(sig)\n",
    "        eps_t.append(u)\n",
    "    plt.plot(eps_t, sig_t)\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "u_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}