Fenics demo.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Let us solve a simple Poisson equation\n",
      "The tutorial on FEniCS can be found [here](http://fenicsproject.org/documentation/tutorial/)\n",
      "\n",
      "$$-\\frac{\\partial^2 u}{\\partial x^2} - \\frac{\\partial^2 u}{\\partial y^2} = f(x, y),$$\n",
      "\n",
      "\n",
      "with (some boundary conditions)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Variational form\n",
      "$$-\\int_{\\Omega} (\\nabla^2 u)dx = \\int_{\\Omega} fv dx,$$ for all $v$,  or by integrating by parts: \n",
      "\n",
      "\n",
      "$$\n",
      "   \\int_{\\Omega} (\\nabla^2 u)dx = -\\int_{\\Omega} (\\nabla u, \\nabla v) dx  - \\int_{\\partial \\Omega} \\frac{\\partial u}{\\partial n} dx\n",
      "$$\n",
      "\n",
      "In the **weak form** we get\n",
      "\n",
      "$$a(u, v) = L(v).$$\n",
      "\n",
      "By taking **test** and **trial** spaces for $u$ and $v$ we get the linear system  \n",
      "\n",
      "$$\n",
      "    \\sum_{ij} A_{ij} u_j = f_j.\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Simple Poisson problem in the FEniCS package"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\"\"\"\n",
      "FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.\n",
      "Simplest example of computation and visualization with FEniCS.\n",
      "\n",
      "-Laplace(u) = f on the unit square.\n",
      "u = u0 on the boundary.\n",
      "u0 = u = 1 + x^2 + 2y^2, f = -6.\n",
      "\"\"\"\n",
      "\n",
      "from dolfin import *\n",
      "\n",
      "# Create mesh and define function space\n",
      "mesh = UnitSquareMesh(6, 4)\n",
      "#mesh = UnitCube(6, 4, 5)\n",
      "V = FunctionSpace(mesh, 'Lagrange', 1)\n",
      "\n",
      "# Define boundary conditions\n",
      "u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')\n",
      "\n",
      "def u0_boundary(x, on_boundary):\n",
      "    return on_boundary\n",
      "\n",
      "bc = DirichletBC(V, u0, u0_boundary)\n",
      "\n",
      "# Define variational problem\n",
      "u = TrialFunction(V)\n",
      "v = TestFunction(V)\n",
      "f = Constant(-6.0)\n",
      "a = inner(nabla_grad(u), nabla_grad(v))*dx\n",
      "L = f*v*dx"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## This is the solve!\n",
      "The solve is quite simple"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute solution\n",
      "u = Function(V)\n",
      "solve(a == L, u, bc)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot solution and mesh\n",
      "plot(u)\n",
      "plot(mesh)\n",
      "\n",
      "# Dump solution to file in VTK format\n",
      "file = File('poisson.pvd')\n",
      "file << u\n",
      "\n",
      "# Hold plot\n",
      "interactive()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now more complicated example\n",
      "The Stokes equation for the fluid flow."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\"\"\"This demo solves the Stokes equations, using quadratic elements for\n",
      "the velocity and first degree elements for the pressure (Taylor-Hood\n",
      "elements). The sub domains for the different boundary conditions\n",
      "used in this simulation are computed by the demo program in\n",
      "src/demo/mesh/subdomains.\"\"\"\n",
      "\n",
      "# Copyright (C) 2007 Kristian B. Oelgaard\n",
      "#\n",
      "# This file is part of DOLFIN.\n",
      "#\n",
      "# DOLFIN is free software: you can redistribute it and/or modify\n",
      "# it under the terms of the GNU Lesser General Public License as published by\n",
      "# the Free Software Foundation, either version 3 of the License, or\n",
      "# (at your option) any later version.\n",
      "#\n",
      "# DOLFIN is distributed in the hope that it will be useful,\n",
      "# but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
      "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n",
      "# GNU Lesser General Public License for more details.\n",
      "#\n",
      "# You should have received a copy of the GNU Lesser General Public License\n",
      "# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.\n",
      "#\n",
      "# Modified by Anders Logg, 2008-2009.\n",
      "#\n",
      "# First added:  2007-11-16\n",
      "# Last changed: 2009-11-26\n",
      "# Begin demo\n",
      "\n",
      "from __future__ import print_function\n",
      "from dolfin import *\n",
      "\n",
      "# Load mesh and subdomains\n",
      "mesh = Mesh(\"../dolfin_fine.xml.gz\")\n",
      "sub_domains = MeshFunction(\"size_t\", mesh, \"../dolfin_fine_subdomains.xml.gz\")\n",
      "\n",
      "#plot(mesh)\n",
      "#plot(sub_domains)\n",
      "\n",
      "# Define function spaces\n",
      "V = VectorFunctionSpace(mesh, \"CG\", 2)\n",
      "Q = FunctionSpace(mesh, \"CG\", 1)\n",
      "W = V * Q\n",
      "\n",
      "# No-slip boundary condition for velocity \n",
      "# x1 = 0, x1 = 1 and around the dolphin\n",
      "noslip = Constant((0, 0))\n",
      "bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0)\n",
      "\n",
      "# Inflow boundary condition for velocity\n",
      "# x0 = 1\n",
      "inflow = Expression((\"-sin(x[1]*pi)\", \"0.0\"))\n",
      "bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1)\n",
      "\n",
      "# Boundary condition for pressure at outflow\n",
      "# x0 = 0\n",
      "zero = Constant(0)\n",
      "bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2)\n",
      "\n",
      "# Collect boundary conditions\n",
      "bcs = [bc0, bc1, bc2]\n",
      "\n",
      "# Define variational problem\n",
      "(u, p) = TrialFunctions(W)\n",
      "(v, q) = TestFunctions(W)\n",
      "f = Constant((0, 0))\n",
      "a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx\n",
      "L = inner(f, v)*dx"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute solution\n",
      "w = Function(W)\n",
      "solve(a == L, w, bcs)\n",
      "\n",
      "# # Split the mixed solution using a shallow copy\n",
      "(u, p) = w.split()\n",
      "\n",
      "# Save solution in VTK format\n",
      "ufile_pvd = File(\"velocity.pvd\")\n",
      "ufile_pvd << u\n",
      "pfile_pvd = File(\"pressure.pvd\")\n",
      "pfile_pvd << p\n",
      "\n",
      "# Plot solution\n",
      "plot(u)\n",
      "plot(p)\n",
      "interactive()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:FFC:Reusing form from cache.\n"
       ]
      }
     ],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "def css_styling():\n",
      "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
      "    return HTML(styles)\n",
      "css_styling()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
        "<style>\n",
        "    @font-face {\n",
        "        font-family: \"Computer Modern\";\n",
        "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
        "    }\n",
        "    div.cell{\n",
        "        /*width:80%;*/\n",
        "        /*margin-left:auto !important;\n",
        "        margin-right:auto;*/\n",
        "    }\n",
        "    h1 {\n",
        "        font-family: 'Alegreya Sans', sans-serif;\n",
        "    }\n",
        "    h2 {\n",
        "        font-family: 'Fenix', serif;\n",
        "    }\n",
        "    h3{\n",
        "\t\tfont-family: 'Fenix', serif;\n",
        "        margin-top:12px;\n",
        "        margin-bottom: 3px;\n",
        "       }\n",
        "\th4{\n",
        "\t\tfont-family: 'Fenix', serif;\n",
        "       }\n",
        "    h5 {\n",
        "        font-family: 'Alegreya Sans', sans-serif;\n",
        "    }\t   \n",
        "    div.text_cell_render{\n",
        "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
        "        line-height: 1.2;\n",
        "        font-size: 160%;\n",
        "        /*width:70%;*/\n",
        "        /*margin-left:auto;*/\n",
        "        margin-right:auto;\n",
        "    }\n",
        "    .CodeMirror{\n",
        "            font-family: \"Source Code Pro\";\n",
        "\t\t\tfont-size: 90%;\n",
        "    }\n",
        "/*    .prompt{\n",
        "        display: None;\n",
        "    }*/\n",
        "    .text_cell_render h1 {\n",
        "        font-weight: 200;\n",
        "        font-size: 50pt;\n",
        "\t\tline-height: 110%;\n",
        "        color:#CD2305;\n",
        "        margin-bottom: 0.5em;\n",
        "        margin-top: 0.5em;\n",
        "        display: block;\n",
        "    }\t\n",
        "    .text_cell_render h5 {\n",
        "        font-weight: 300;\n",
        "        font-size: 16pt;\n",
        "        color: #CD2305;\n",
        "        font-style: italic;\n",
        "        margin-bottom: .5em;\n",
        "        margin-top: 0.5em;\n",
        "        display: block;\n",
        "    }\n",
        "    \n",
        "    li {\n",
        "        line-height: 110%;\n",
        "    }\n",
        "    .warning{\n",
        "        color: rgb( 240, 20, 20 )\n",
        "        }  \n",
        "</style>\n",
        "<script>\n",
        "    MathJax.Hub.Config({\n",
        "                        TeX: {\n",
        "                           extensions: [\"AMSmath.js\"]\n",
        "                           },\n",
        "                tex2jax: {\n",
        "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
        "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
        "                },\n",
        "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
        "                \"HTML-CSS\": {\n",
        "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
        "                }\n",
        "        });\n",
        "</script>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<IPython.core.display.HTML at 0x10e232bd0>"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}