ipy_notebook/convergence_check_multi.v3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(r'../')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "\n",
    "import cell_computation as com\n",
    "import cell_geom as ce\n",
    "import cell_material as ma\n",
    "from copy import deepcopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logging.getLogger('FFC').setLevel(logging.WARNING)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(r\"../m.xml\")\n",
    "# mesh = Mesh(r\"m_fine.xml\")\n",
    "\n",
    "cell = ce.UnitCell(mesh)\n",
    "\n",
    "inc = ce.InclusionCircle(2, (0.5, 0.5), 0.25)\n",
    "inc_di = {'circle_inc': inc}\n",
    "cell.set_append_inclusion(inc_di)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "VFS = VectorFunctionSpace(cell.mesh, \"CG\", 1,\n",
    "                              constrained_domain=ce.PeriodicBoundary_no_corner(2))\n",
    "FS = FunctionSpace(cell.mesh, \"CG\", 1,\n",
    "                   constrained_domain=ce.PeriodicBoundary_no_corner(2))\n",
    "\n",
    "# Set materials\n",
    "E_m, nu_m, Kappa_m = 2e5, 0.4, 7.\n",
    "# n = 1000\n",
    "n = 1000  # 13.Jan\n",
    "E_i, nu_i, Kappa_i = 1000 * E_m, 0.3, n * Kappa_m\n",
    "\n",
    "mat_m = ma.neo_hook_eap(E_m, nu_m, Kappa_m)\n",
    "mat_i = ma.neo_hook_eap(E_i, nu_i, Kappa_i)\n",
    "mat_li = [mat_m, mat_i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = Function(VFS)\n",
    "el_pot_phi = Function(FS)\n",
    "strain_space_w = TensorFunctionSpace(mesh, 'DG', 0)\n",
    "strain_space_E = VectorFunctionSpace(mesh, 'DG', 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deform_grad_with_macro(F_bar, w_component):\n",
    "    return F_bar + grad(w_component)\n",
    "def e_field_with_macro(E_bar, phi):\n",
    "    return E_bar - grad(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "comp = com.MicroComputation(cell, mat_li,\n",
    "                        [deform_grad_with_macro, e_field_with_macro],\n",
    "                        [strain_space_w, strain_space_E])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation with FD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample_num = 8\n",
    "# delta = np.logspace(-2,-4,num=sample_num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_mer_stress(F_bar, E_bar):\n",
    "    comp.input([F_bar, E_bar], [w, el_pot_phi])\n",
    "    comp.comp_fluctuation()\n",
    "    return comp.avg_merge_stress()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conv_check_component(label, F_bar, E_bar, compo, delta):\n",
    "    C_eff_component_FD = np.zeros(shape=(len(delta),6), dtype=float)\n",
    "    if label is 'F':\n",
    "        for i, d in enumerate(delta):\n",
    "            F_minus = deepcopy(F_bar)\n",
    "            F_minus[compo] = F_bar[compo] - d/2\n",
    "            F_plus = deepcopy(F_bar)\n",
    "            F_plus[compo] = F_bar[compo] + d/2\n",
    "\n",
    "            P_minus = avg_mer_stress(F_minus, E_bar)\n",
    "            P_plus  = avg_mer_stress(F_plus, E_bar)\n",
    "            \n",
    "            C_eff_component_FD[i,:] = (P_plus - P_minus)/d\n",
    "    elif label is 'E':\n",
    "        for i, d in enumerate(delta):\n",
    "            E_minus = deepcopy(E_bar)\n",
    "            E_minus[compo] = E_bar[compo] - d/2\n",
    "            E_plus = deepcopy(E_bar)\n",
    "            E_plus[compo] = E_bar[compo] + d/2\n",
    "\n",
    "            P_minus = avg_mer_stress(F_bar, E_minus)\n",
    "            P_plus  = avg_mer_stress(F_bar, E_plus)\n",
    "            \n",
    "            C_eff_component_FD[i,:] = (P_plus - P_minus)/d\n",
    "    else:\n",
    "        raise Exception('no such field label')\n",
    "    \n",
    "    return C_eff_component_FD "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n",
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge stress computation finished\n"
     ]
    }
   ],
   "source": [
    "F_bar = [1.1, 0., 0.1, 1.]\n",
    "\n",
    "E_bar = [0., 0.2]\n",
    "\n",
    "delta = [0.01, 0.01/2, 0.01/4, 0.01/8]\n",
    "\n",
    "C_eff_component_FD = conv_check_component('F', F_bar, E_bar, 3, delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homogenization Method Result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fluctuation computation finished\n",
      "strain computation finished\n",
      "average merge moduli computation finished\n"
     ]
    }
   ],
   "source": [
    "comp = com.MicroComputation(cell, mat_li,\n",
    "                        [deform_grad_with_macro, e_field_with_macro],\n",
    "                        [strain_space_w, strain_space_E])\n",
    "comp.input([F_bar, E_bar], [w, el_pot_phi])\n",
    "comp.comp_fluctuation()\n",
    "C_eff = comp.effective_moduli_2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  3.23014568e+04  -9.13021574e+03  -1.39145501e+03   7.36837080e+04\n",
      "   -1.46915356e-12   1.05645659e-11]\n",
      " [  3.23005307e+04  -9.13012721e+03  -1.39162376e+03   7.36830725e+04\n",
      "   -1.46915153e-12   1.05644829e-11]\n",
      " [  3.23002991e+04  -9.13010507e+03  -1.39166594e+03   7.36829136e+04\n",
      "   -1.46915102e-12   1.05644622e-11]\n",
      " [  3.23002412e+04  -9.13009954e+03  -1.39167649e+03   7.36828739e+04\n",
      "   -1.46915089e-12   1.05644570e-11]]\n"
     ]
    }
   ],
   "source": [
    "print (C_eff_component_FD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  3.23002219e+04,  -9.13009769e+03,  -1.39168000e+03,\n",
       "         7.36828607e+04,  -1.46915085e-12,   1.05644552e-11])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_eff[:,3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convergence Check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "component = C_eff[:,3]\n",
    "\n",
    "tmp = np.outer(np.ones((len(delta),1)),np.transpose(component))\n",
    "\n",
    "error = np.linalg.norm(tmp - C_eff_component_FD, axis=1)/np.linalg.norm(component)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.87577323e-05,   4.68943878e-06,   1.17237868e-06,\n",
       "         2.93114628e-07])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}