example/Manual on cell_computation.py.v3.ipynb

{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  },
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Table of Contents\n",
      "* [Overview](#Overview)\n",
      "* [Example on 2D uni field unit cell computation](#Example-on-2D-uni-field-unit-cell-computation)\n",
      "* [Concepts in implementation](#Concepts-in-implementation)\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Overview"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the main part of the module consisting of class `MicroComputation` and several assisting functions which cover setup of the solver parameters and relation between extended strain and displacement."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Example on 2D uni field unit cell computation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Module import__"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "from dolfin import *\n",
      "\n",
      "import numpy as np\n",
      "\n",
      "import sys\n",
      "sys.path.append('../')\n",
      "\n",
      "import cell_geom as geom\n",
      "import cell_material as mat\n",
      "import cell_computation as comp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Choose linear backend__\n",
      "\n",
      "The embedded backends can be viewed through `linear_algebra_backend()`, namely _Eigen_, _PETSc_ (default), and _STL_. Notice that the backend should be specified at the beginning of all the calculation and initialization, as the corresponding matrix and vector objects are casted in a style appropriate for the chose backend."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Linear Backend\n",
      "parameters['linear_algebra_backend'] = 'PETSc'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Define unit cell geometry (mesh, inclusions)__\n",
      "\n",
      "The geometrical properties are set according to the _Manual on cell_geom.py_. Inclusions such as circle and rectangular in 2D and 3D are available. Extensions are possible in this context. Mesh can be imported the same as in _FEniCS_."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Define Geometry\n",
      "mesh = Mesh(r'../m_fine.xml')\n",
      "cell = geom.UnitCell(mesh)\n",
      "\n",
      "# Add inclusion\n",
      "inc = geom.InclusionCircle(2, (0.5, 0.5), 0.25)\n",
      "inc_di = {'circle_inc': inc}\n",
      "cell.set_append_inclusion(inc_di)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Define materials in composites__\n",
      "\n",
      "The material properties are set according to the [manual on cell_material.py][mat help]. Three types of materials are provided in the material libraries. Users can also specify their own materials. The definition of a new material follows the steps, \n",
      " 1. free energy function `psi`, \n",
      " 2. invariant relations with physical variables, \n",
      " 3. initialize an instance of `Material` with `psi` and required parameters, \n",
      " 4. append invariant relations to the `Material` instance, \n",
      " 5. call this instance with the physical field variables and complete instantiation\n",
      "\n",
      "Details and restrictions are referred in [manual on cell_material.py][mat help].\n",
      "\n",
      "[mat help]:http://localhost:8888/notebooks/Manual%20on%20cell_material.py.ipynb"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Define Material\n",
      "E_m, nu_m, E_i, nu_i = 10.0, 0.3, 1000.0, 0.3\n",
      "mat_m = mat.st_venant_kirchhoff(E_m, nu_m)\n",
      "mat_i = mat.st_venant_kirchhoff(E_i, nu_i)\n",
      "\n",
      "mat_li = [mat_m, mat_i]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Create MicroComputation instance (pre-processing)__\n",
      "\n",
      "Here a general relation between the physical field variable and strain is given. This is realized in the function `deform_grad_with_macro`. `w` stands for the fluctuation that needs to be solved, while `F` is a general strain measure, which are considered in the material total energy. `F_bar` is the general strain from macro computation serving as input for `MicroComputation`.\n",
      "\n",
      "The instantiation steps are as follows\n",
      " 1. Function spaces for physical field variables and generate Functions for field variables. These field variables are often regarded as general displacements, and in `MicroComputation` they are fluctuation to solve.\n",
      " 2. Define relation between general strains and general displacements.\n",
      " 3. Give function space for general strain, these are for the post processing.\n",
      " 4. Initialize a `MicroComputation` instance using geometry (`cell`), materials (`mat_li`), general strain-displacement relation (`deform_grad_with_macro`), and the strain function space for post-processing (`strain_space`).\n",
      " 5. Specify general strain from macro field and input them into instance of `MicroComputation`. It completes the creation a `MicroComputation`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Define Computation\n",
      "# Step 1: Field variables\n",
      "VFS = VectorFunctionSpace(cell.mesh, \"CG\", 1, \n",
      "                          constrained_domain=geom.PeriodicBoundary_no_corner(2))\n",
      "w = Function(VFS)\n",
      "\n",
      "# Step 2: Strain and field variable relations\n",
      "def deform_grad_with_macro(F_bar, w_component):\n",
      "    return F_bar + grad(w_component)\n",
      "\n",
      "# Step 3: Strain space for post processing\n",
      "strain_space = TensorFunctionSpace(mesh, 'DG', 0)\n",
      "\n",
      "# Step 4: Initialization\n",
      "compute = comp.MicroComputation(cell, mat_li, \n",
      "                                [deform_grad_with_macro],\n",
      "                                [strain_space])\n",
      "# Step 5: Complete instantiation\n",
      "F_bar = [0.9, 0., 0., 1.]\n",
      "compute.input([F_bar], [w])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Solve fluctuation__\n",
      "\n",
      "It begins with setting of solver parameters. The solving step is just call the member function of this instance `comp_fluctuation()`\n",
      "\n",
      "There are two classes of methods, _direct_ (lu solver) and _iterative_ (krylov solver). To list the solvers in direct methods, command `lu_solver_methods()`, while for listing iterative method solvers command is `krylov_solver_methods()`. For iterative solvers a preconditioner is needed, which can be viewed using command `krylov_solver_preconditioners()`. A complete summary of solvers can be referred in the [website of _PETSc_][PETSc solver].\n",
      "\n",
      "[PETSc solver]:http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Parameters for solving of fluctuation\n",
      "comp.set_solver_parameters('non_lin_newton', lin_method='direct', linear_solver='mumps')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".-------------------.\n",
        "| Solver Parameters |\n",
        ".-------------------.\n",
        "direct method is used\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "compute.comp_fluctuation(print_progress=True, print_solver_info=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "fluctuation computation finished\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Post Processing and view the results__\n",
      "\n",
      "Calculation of homogenized quantity is the central part of homogenization method. Hence various post processing of the `MicroComputation` result are implemented in this class. These are `comp_strain()`, `comp_stress()` in terms of calculating general strain and stress in the unit cell, `avg_merge_strain()`, `avg_merge_stress()`, `avg_merge_moduli()` in terms of calculating the averaged quantities in a unit cell, and `effective_moduli_2()` for calculating the homogenized tangent moduli.\n",
      "\n",
      "`effective_moduli_2()` is the most consuming part. Specifying a good solver for it can speed up this process. This is involved with using function `set_post_solver_parameters()`. \n",
      "\n",
      "Visualizing the results are also realized in the current module. Wrapping the visualization tools in _FEniCS_ is included. The member functions are `view_fluctuation()`, `view_displacement()`, and `view_post_processing()`. When multiple fields are considered, specifying the component of visualization is needed."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "compute.comp_strain()\n",
      "compute.comp_stress()\n",
      "compute.avg_merge_strain()\n",
      "compute.avg_merge_stress()\n",
      "compute.avg_merge_moduli()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "strain computation finished\n",
        "strain computation finished\n",
        "stress computation finished\n",
        "average merge strain computation finished\n",
        "average merge stress computation finished\n",
        "average merge moduli computation finished\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "array([[  2.68263384e+02,  -1.47306216e-02,   1.45226354e-02,\n",
        "          1.16307090e+02],\n",
        "       [ -1.47306216e-02,   7.64588022e+01,   7.75364924e+01,\n",
        "         -1.46396985e-02],\n",
        "       [  1.45226354e-02,   7.75364924e+01,   7.64302939e+01,\n",
        "          1.44729548e-02],\n",
        "       [  1.16307090e+02,  -1.46396985e-02,   1.44729548e-02,\n",
        "          2.72444929e+02]])"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Post processing solver parameters\n",
      "comp.set_post_solver_parameters(lin_method='iterative',)\n",
      "\n",
      "# Homogenized tangent moduli\n",
      "compute.effective_moduli_2()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "+----------------------------+\n",
        "| Post Processing Parameters |\n",
        "+----------------------------+\n",
        "iterative method is used\n",
        "a valid preconditioner should be provided\n",
        "average merge moduli computation finished\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([[  1.05215604e+01,   6.81100019e-04,   8.10922941e-04,\n",
        "          6.09319740e+00],\n",
        "       [  6.81100019e-04,   3.03957695e+00,   4.12022593e+00,\n",
        "         -2.43801203e-04],\n",
        "       [  8.10922941e-04,   4.12022593e+00,   2.98311033e+00,\n",
        "         -3.19622254e-04],\n",
        "       [  6.09319740e+00,  -2.43801203e-04,  -3.19622254e-04,\n",
        "          1.74423266e+01]])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# View results\n",
      "compute.view_fluctuation()\n",
      "compute.view_displacement()\n",
      "compute.view_post_processing('strain', (0,1))\n",
      "compute.view_post_processing('stress', (0,1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- __Call help for each functions__\n",
      "\n",
      "Note that in-code docstrings are given to explain functions arguments and outputs. Calling docstrings in IPython is simply add a question mark behind a function, and two question marks will show the detailed implementation the method. In Python context, `help()` is used to list the docstring."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "help(compute.comp_fluctuation)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Help on method comp_fluctuation in module cell_computation:\n",
        "\n",
        "comp_fluctuation(self, print_progress=False, print_solver_info=False) method of cell_computation.MicroComputation instance\n",
        "    Solve fluctuation, solver parameters are set before solving\n",
        "    \n",
        "    :param print_progress: (bool) print detailed solving progress\n",
        "    :param print_solver_info: (bool) print detailed solver info\n",
        "    \n",
        "    :return: updated self.w_merge\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "collapsed": true
     },
     "source": [
      "__This example can work as a [simulation template][template].__\n",
      "\n",
      "[template]:http://localhost:8888/notebooks/simulation_template.ipynb"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Concepts in implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The design concepts of this module are elaborated in this part. The main idea of the `MicroComputation` is to establish a unit cell. When the deformation or field variables from macro state are input, micro state computation can be realized. This computation will fall into several parts such as solving fluctuation, and in the end calculating effective tangent moduli.\n",
      "\n",
      "If the geometry and materials of this unit cell does not change. This instance of `MicroComputation` can be reused for another computation. This will make `MicroComputation` act like a factory that produce the result of computation in the micro scale.\n",
      "\n",
      "Various functions and member methods are targeted to achieve this computation. \n",
      "\n",
      "- __Initializtion__\n",
      "\n",
      "When an instance is initialized, several methods are called behind scene. Once `input()` is called, the instance is complete. Arguments for `input()` are `F_bar_li` and `w_li` representing the macro field input and to be solved fluctuation respectively. In `input()` `F_bar_li` is transformed to a list of Functions suitable for the later formulation done by `_F_bar_init()`. `set_field()` is called to merge and split fields for the variational formulaiton later on. Then general strain measures are created through `extend_strain()` also for the later formulation, where the input is the splitted _FEniCS_ Functions.\n",
      "\n",
      "- __Pre-processing__\n",
      "\n",
      "Pre-processing is the problem formulation stage. `_total_energy()` is to insert general strains into materials to construct its dependency within the physical fields. The summation of all the free enerygies of every material components in composite make the total energy complete.\n",
      "\n",
      "Then boundary condition is provided with `_bc_fixed_corner()`, which fix all the corners of the unit cell.\n",
      "\n",
      "`_fem_formulation_composite` follows with derivation of the nonlinear problem using powerful functions defined in _FEniCS_, `derivative()`\n",
      "\n",
      "- __Solving__\n",
      "\n",
      "Solving step is accomplished by `comp_fluctuation()` with previously specified solver parameters.\n",
      "\n",
      "- __Post-processing__\n",
      "\n",
      "Post-processing plays the key role in homogenization problem. It starts with `comp_strain`, where convenient _FEniCS_ function `project()` is used. Then material energy is updated with the calculated general strains in `_energy_update()`. It is for the purpose that, general stresses are conjugated part in the total energy. Formulating the total energy in general strain will lead to direct calculation of general stresses with the help of `derivative()`. It is done by `comp_stress()`. \n",
      "\n",
      "`avg_merge_strain()`, `avg_merge_stress()`, and `avg_merge_moduli()` are implemented using the trick of multiplying trial and test functions defining on a constant Function Space. Detailed elaboration is given in the report.\n",
      "\n",
      "`effective_moduli_2()` is based on `avg_merge_moduli()`, a LTKL term is subtracted from the averaged merged moduli. The calculation of this term is realized in `sensitivity()`. L matrix is assembled with boundary conditions. Some techniques are taken in imposing the boundary conditions. Using the _FEniCS_ `solve()` on `K_a` and `L[:, i]` gives a intermediate result. Left multiplying with the transpose of L gives the required term.\n",
      "\n",
      "- __Things to notice__\n",
      "\n",
      "In the module only the relation between mechanical displacement and deformation gradient is given in the function `deform_grad_with_macro`, other kinds of relation between general strain and displacement should be specified by the user.\n",
      "\n",
      "Another thing to notice is that a good solver needs to be chosen in complicated cases such as multi field or 3D. Direct solvers are rather slow in this circumstances, while iterative solvers will not always converge and requires a lot of try-outs with the right preconditioners."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}