{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append(\"build-clang-5.0/language_bindings/python\")\n",
    "import muSpectre as msp\n",
    "import matplotlib.pyplot as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "ename": "RuntimeError",
     "evalue": "Failure at load step 3 of 30. In Newton-Raphson step 4:\nConjugate gradient has not converged. After 1000 steps, the solver  FAILED with  |r|/|b| =            -nan, cg_tol = 1e-06\n\nThe applied boundary condition is F =\n    0 0.003\n    0     0\nand the load increment is ΔF =\n    0 0.001\n    0     0\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-995a33fb918b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     71\u001b[0m \u001b[0;31m## de_geus: is described in de Geus et al. see Ref above\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     72\u001b[0m \u001b[0;31m## newton_cg: classical Newton-Conjugate Gradient solver. Recommended\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 73\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmsp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolvers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnewton_cg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrve\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDel0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meq_tol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m: Failure at load step 3 of 30. In Newton-Raphson step 4:\nConjugate gradient has not converged. After 1000 steps, the solver  FAILED with  |r|/|b| =            -nan, cg_tol = 1e-06\n\nThe applied boundary condition is F =\n    0 0.003\n    0     0\nand the load increment is ΔF =\n    0 0.001\n    0     0\n"
     ]
    }
   ],
   "source": [
    "## currently, muSpectre is restricted to odd-numbered resolutions for\n",
    "## reasons explained in T.W.J. de Geus, J. Vondřejc, J. Zeman,\n",
    "## R.H.J. Peerlings, M.G.D. Geers, Finite strain FFT-based non-linear\n",
    "## solvers made simple, Computer Methods in Applied Mechanics and\n",
    "## Engineering, Volume 318, 2017\n",
    "## https://doi.org/10.1016/j.cma.2016.12.032\n",
    "resolution = [51, 51]\n",
    "center = np.array([r//2 for r in resolution])\n",
    "incl = resolution[0]//5\n",
    "\n",
    "\n",
    "## Domain dimensions\n",
    "lengths = [7., 5.]\n",
    "## formulation (small_strain or finite_strain)\n",
    "formulation = msp.Formulation.finite_strain\n",
    "\n",
    "## build a computational domain\n",
    "rve = msp.Cell(resolution, lengths, formulation)\n",
    "\n",
    "## define the material properties of the matrix and inclusion\n",
    "hard = msp.material.MaterialHyperElastoPlastic1_2d.make(\n",
    "    rve, \"hard\", 10e9, .33, .006e9,  h=.01e9)\n",
    "#hard = msp.material.MaterialLinearElastic1_2d.make(\n",
    "#    rve, \"hard\", 10e9, .33)\n",
    "soft = msp.material.MaterialLinearElastic1_2d.make(\n",
    "    rve, \"soft\",  70e9, .33)\n",
    "\n",
    "## assign each pixel to exactly one material\n",
    "for i, pixel in enumerate(rve):\n",
    "    if np.linalg.norm(center - np.array(pixel),2)<incl:\n",
    "        hard.add_pixel(pixel)\n",
    "    else:\n",
    "        soft.add_pixel(pixel)\n",
    "\n",
    "## define the convergence tolerance for the Newton-Raphson increment\n",
    "tol = 1e-5\n",
    "## tolerance for the solver of the linear cell\n",
    "cg_tol = 1e-6\n",
    "eq_tol = 1e-8\n",
    "\n",
    "\n",
    "## Macroscopic strain\n",
    "def get_Del0(gammas):\n",
    "    def Del0(gamma):\n",
    "        Del0 = np.array([[.0, gamma],\n",
    "                         [0,  0.]])\n",
    "        return Del0#.5*(Del0 + Del0.T)\n",
    "    return [Del0(gamma) for gamma in gammas]\n",
    "\n",
    "del_gamma = 1e-3 # increments\n",
    "Del_gamma = 3e-2 # total strain\n",
    "nb_gamma = int(Del_gamma/del_gamma)\n",
    "Del0 = get_Del0((i*del_gamma for i in range(1, nb_gamma+1)))\n",
    "\n",
    "\n",
    "maxiter = 1000 ## for linear cell solver\n",
    "\n",
    "## Choose a solver for the linear cells. Currently avaliable:\n",
    "## SolverCG, SolverCGEigen, SolverBiCGSTABEigen, SolverGMRESEigen,\n",
    "## SolverDGMRESEigen, SolverMINRESEigen.\n",
    "## See Reference for explanations\n",
    "solver = msp.solvers.SolverCG(rve, cg_tol, maxiter, verbose=False)\n",
    "\n",
    "\n",
    "## Verbosity levels:\n",
    "## 0: silent,\n",
    "## 1: info about Newton-Raphson loop,\n",
    "verbose = 1\n",
    "\n",
    "## Choose a solution strategy. Currently available:\n",
    "## de_geus: is described in de Geus et al. see Ref above\n",
    "## newton_cg: classical Newton-Conjugate Gradient solver. Recommended\n",
    "result = msp.solvers.newton_cg(rve, Del0, solver, tol, eq_tol, verbose=verbose)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mean_stress = [res.stress.reshape(-1, 4).mean(axis=0) for res in result]\n",
    "\n",
    "mean_tau = np.array([stress[2] for stress in mean_stress])\n",
    "gammas = np.array([del0[0,1] for del0 in Del0])\n",
    "plt.plot(gammas, mean_tau)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "collapsed": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "name": "plastic_check.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}