{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# L\u00f6sen der Poisson-Gleichung mit der Methode der Finiten Elemente\n",
      "\n",
      "Es soll die Gleichung\n",
      "\n",
      "        -\u0394u = f(x,y) in   \u03a9 = (0,1)\u00b2\n",
      "     u(x,y) = g(x,y) auf \u2202\u03a9\n",
      "\n",
      "f\u00fcr gegebene Funktionen $f$ und $g$ mit der Methode der Finite Elemente auf einem 128x128-Dreiecksgitter gel\u00f6st werden.\n",
      "Hierbei sollen Lagrange-Elemente erster Ordnung (P1-Funktionen) f\u00fcr Ansatz- und Testraum verwendet werden."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import dune.common\n",
      "import dune.grid\n",
      "import numpy as np\n",
      "from scipy.sparse import lil_matrix\n",
      "from scipy.sparse.linalg import spsolve"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 1\n",
      "\n",
      "\u00dcbernehmen Sie die Klassen `P1LocalBasis` und `P1Function` von Ihrer L\u00f6sung des zweiten \u00dcbungsblatts oder dem L\u00f6sungsvorschlag."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 2\n",
      "\n",
      "Implementieren Sie eine Klasse `P1Index`, deren Methode `__call__` dem lokalen Freiheitsgrad `i` auf dem Element `element` die Nummer des zugeh\u00f6rigen globalen Freiheitsgrads zuordnet."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class P1Index:\n",
      "    def __init__(self, grid):\n",
      "        self._grid = grid\n",
      "    def __call__(self, element, i):\n",
      "        pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 3\n",
      "\n",
      "Implementieren Sie eine neue Klasse `LocalAssembler`.\n",
      "Diese soll die Werte\n",
      "$$\n",
      "  A_{T;i,j} = a_T(\\phi_{T;i}, \\phi_{T;j}) = \\int_T \\nabla\\phi_{T;i} \\nabla\\phi_{T;j}\n",
      "$$\n",
      "und\n",
      "$$\n",
      "  b_{T;i} = l_T(\\phi_{T;i}) = \\int_T f \\phi_{T;i}\n",
      "$$\n",
      "f\u00fcr die lokalen Formfunktionen $\\phi_{T;i}$ ($i\\in\\{0,1,2\\}$) und rechte Seite $f$ berechnen.\n",
      "\n",
      "Implementieren Sie dazu die Methode `assembleLocal`, die die Integrale auf dem Element `element` mit einer Quadraturformel zweiter Ordnung berechnet und ein Paar aus der lokalen Matrix `Alocal` und der lokalen rechten Seite `blocal` zur\u00fcckgibt.\n",
      "\n",
      "Verwenden Sie `self._rhs` als Funktion $f: \\Omega\\to\\mathbb{R}$ und `self._basis` mit den Methoden `evaluateFunction` und `evaluateJacobian` wie in den vorherigen \u00dcbungen implementiert als Formfunktionen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class LocalAssembler:\n",
      "    def __init__(self, grid, basis, f, nDofsPerElement):\n",
      "        self._grid = grid\n",
      "        self._basis = basis\n",
      "        self._f = f\n",
      "\n",
      "        # Anzahl der Freiheitsgrade auf einem einzelnen Element\n",
      "        self._nDofsPerElement = nDofsPerElement\n",
      "    \n",
      "    def assembleLocal(self, element):\n",
      "        pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 4\n",
      "\n",
      "Implementieren Sie eine neue Klasse `GlobalAssembler`.\n",
      "Diese soll die Matrix $A$ mit $$A_{i,j} = a(\\phi_i, \\phi_j)$$ und die rechte Seite $b$ mit $$b_i = l(\\phi_i)$$ aufstellen.\n",
      "\n",
      "Implementieren Sie dazu die Methode `assemble`.\n",
      "Diese soll zun\u00e4chst die Matrix `A` als leere `lil_matrix` und `b` als Nullvektor geeigneter Gr\u00f6\u00dfe anlegen.\n",
      "Anschlie\u00dfend soll auf jedem Element `element` mit dem `LocalAssembler` aus Aufgabe 2 der lokale Anteil an `A` bzw. `b` berechnet werden und zu `A` bzw. `b` hinzuaddiert werden.\n",
      "Schlie\u00dflich soll das Paar `A`, `b` zur\u00fcckgegeben werden.\n",
      "\n",
      "Verwenden Sie `self._index` als Funktion $g_T: \\{0, ..., n\\} \\to \\mathbb{N}$ um zu den lokalen Freiheitsgraden $i$ auf dem Element $T$ die zugeh\u00f6rige Nummer des globalen Freiheitsgrades $g_T(i)$ zu berechnen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class GlobalAssembler:\n",
      "    def __init__(self, grid, localAssembler, index, nDofs):\n",
      "        self._grid = grid\n",
      "        self._localAssembler = localAssembler\n",
      "        self._index = index\n",
      "\n",
      "        # Anzahl aller Freiheitsgrade\n",
      "        self._nDofs = nDofs\n",
      "    \n",
      "    def assemble(self):\n",
      "        pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 5\n",
      "\n",
      "Implementieren Sie die Funktion `dirichletP1(grid, g, A, b)`.\n",
      "Die Funktion soll f\u00fcr die Freiheitsgrade $i$ mit zugeh\u00f6rigen Koordinaten $x_i$, die auf $\\partial\\Omega$ liegen, die Dirichlet-Randbedingung implementieren.\n",
      "Setzen Sie dazu $A_{i,j} = \\delta_{i,j}$ und $b_i = g(x_i)$ f\u00fcr alle diese $i$.\n",
      "\n",
      "Sie k\u00f6nnen die Funktion `clearRow(A, i)` verwenden, um vorhandene Nicht-Null-Eintr\u00e4ge aus der `i`-te Zeile von `A` zu entfernen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def clearRow(A, i):\n",
      "    \"remove entries from row `i` from lil_matrix `A`\"\n",
      "    A.data[i] = list()\n",
      "    A.rows[i] = list()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def isDirichlet(x):\n",
      "    return x[0] < 1e-3 or x[0] > 0.999 or x[1] < 1e-3 or x[1] > 0.999"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def dirichletP1(grid, g, A, b):\n",
      "    pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 6\n",
      "\n",
      "Implementieren Sie die Funktion `poisson(n, f, g)`, die das Poissonproblem auf einem strukturierten $n\\times n$-Dreiecksgitter l\u00f6st und als VTK-Datei ausgibt.\n",
      "Dabei ist $f: \\Omega \\to \\mathbb{R}$ die rechte Seite und $g: \\Omega \\to \\mathbb{R}$ die Funktion, die Dirichlet-Randwerte auf $\\partial\\Omega$ implementiert.\n",
      "\n",
      "Verwenden Sie den `GlobalAssembler` mit der `P1Index`-Klasse und dem `LocalAssembler` mit der `P1LocalFunction`-Klasse um `A` und `b` aufzustellen und `dirichletP1` um die Dirichlet-Randbedingung zu ber\u00fccksichtigen.\n",
      "L\u00f6sen Sie das Gleichungssystem mit `spsolve` wie auf dem ersten \u00dcbungsblatt.\n",
      "\n",
      "Die Ausgabe als VTK-Datei kann mit der Funktion `output` erfolgen; der Parameter `flocal` sollte hierbei eine Instanz der Klasse `P1Function` sein."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def output(grid, flocal, filename=\"04-FE-Poisson\"):\n",
      "    gridfunction = grid.localGridFunction(lambda e, x: [flocal(e, x)])\n",
      "    grid.writeVTK(filename, pointdata={\"u\": gridfunction})"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def poisson(n, f, g):\n",
      "    pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 7\n",
      "\n",
      "Rufen Sie die Funktion `poisson` f\u00fcr die Funktionen\n",
      "$$f(x, y) = 20 \\pi^2 \\sin(10\\pi x) * \\cos(20\\pi y)$$\n",
      "und\n",
      "$$g(x, y) = \\sin(\\pi x)$$\n",
      "auf, und betrachten Sie die L\u00f6sung mit dem Programm `paraview`.\n",
      "Reichen Sie bei Ihrer Abgabe auch einen Screenshot der L\u00f6sung ein (in ParaView: File -> Save Screenshot...).\n",
      "\n",
      "**Hinweis:**\n",
      "Nach dem \u00d6ffnen der Datei in ParaView m\u00fcssen Sie links auf die \"Apply\"-Schaltfl\u00e4che klicken, damit die Funktion angezeigt wird.\n",
      "Sie k\u00f6nnen sich auch das Gitter anzeigen lassen, indem Sie in der Werkzeugleiste oben \"Surface With Edges\" statt dem vorgegebenen \"Surface\" ausw\u00e4hlen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 128\n",
      "f = lambda x: 20 * np.pi**2 * np.sin(10*np.pi*x[0]) * np.cos(20*np.pi*x[1])\n",
      "g = lambda x: np.sin(np.pi * x[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "poisson(n, f, g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}