{
 "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.3"
  },
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Diskretisierung der Stokes-Gleichung mit dem Taylor-Hood-Element\n",
      "\n",
      "In dieser Aufgabe soll die Stokes-Gleichung\n",
      "\\begin{align}\n",
      "-\\mu \\Delta \\mathbf{v} + \\nabla p & = \\mathbf{f} \\\\\n",
      "\\nabla \\cdot \\mathbf{v} & = 0\n",
      "\\end{align}\n",
      "auf $\\Omega = (0,1)^2$ mit dem Taylor-Hood-Element gel\u00f6st werden.\n",
      "Hierbei steht $\\mathbf{v}$ f\u00fcr den Geschwindigkeitsvektor, $p$ f\u00fcr den Druck und $\\mathbf{f}$ f\u00fcr eine Kraftdichte.\n",
      "Der Parameter $\\mu$ beschreibt die Viskosit\u00e4t der Fl\u00fcssigkeit.\n",
      "\n",
      "Als Randbedingung wird\n",
      "\\begin{equation}\n",
      "\\mathbf{v}(x) = \\mathbf{v}_0(x) \\begin{cases}\n",
      "  (1, 0)^T & \\text{ falls } x_1 = 1 \\\\\n",
      "  (0, 0)^T & \\text{sonst}\n",
      "  \\end{cases}\n",
      "\\end{equation}\n",
      "f\u00fcr $x \\in \\partial \\Omega$\n",
      "und\n",
      "\\begin{equation}\n",
      "p(x) = 0\n",
      "\\end{equation}\n",
      "f\u00fcr $x \\in \\partial \\Omega$ mit $x_1 = 0$ gew\u00e4hlt.\n",
      "\n",
      "In der schwachen Formulierung werden Funktionen $(\\mathbf{v}, p) \\in V^2 \\times W$ gesucht, so dass\n",
      "\\begin{align}\n",
      "\\int \\mu \\nabla \\psi \\cdot \\nabla \\mathbf{v} - \\nabla \\cdot \\psi p & = \\int \\psi \\mathbf{f} && \\text{ f\u00fcr alle } \\psi\\in V_0^2 \\\\\n",
      "\\int \\phi \\nabla \\cdot \\mathbf{v} & = 0 && \\text{ f\u00fcr alle } \\phi \\in W_0.\n",
      "\\end{align}\n",
      "\n",
      "F\u00fcr das Taylor-Hood-Element sind $V$ die Lagrange-Elemente zweiter Ordnung ($P_2$) und $W$ Lagrange-Elemente erster Ordnung ($P_1$).\n",
      "Die R\u00e4ume der Testfunktionen $V_0$ und $W_0$ sind die Teilmenge von $V$ bzw. $W$, die auf dem jeweiligen Dirichlet-Rand verschwinden."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ihr Programm sollte mit den folgenden Import-Anweisungen auskommen:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy.sparse import lil_matrix\n",
      "from scipy.sparse.linalg import spsolve\n",
      "\n",
      "import dune.common\n",
      "import dune.functions\n",
      "import dune.grid"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 1\n",
      "\n",
      "Studieren Sie das Notebook `05-Poisson-with-dune-functions.ipynb`.\n",
      "\n",
      "In diesem wird das Poisson-Problem der letzten \u00dcbung mit FE-Funktionen aus dem `dune-functions`-Modul diskretisiert.\n",
      "Diese Implementierung sollen Sie als Grundlage f\u00fcr die Bearbeitung der folgenden Aufgaben verwenden."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 2\n",
      "\n",
      "Implementieren Sie eine Klasse `LocalAssembler`, die den Anteil eines einzelnen Dreiecks an der Systemmatrix und der rechten Seite berechnet.\n",
      "Verwenden Sie hierbei eine Quadraturformel zweiter Ordnung.\n",
      "\n",
      "Implementieren Sie zun\u00e4chst die Methoden `indexv(self, i, k)` und `indexp(self, i)`, die die Nummerierung der Freiheitsgrade f\u00fcr Geschwindigkeit und Druck bereitstellen.\n",
      "\n",
      "In der Methode `indexv` bezeichnet `i` den Freiheitsgrad auf dem Dreieck und `k` die Komponente des Geschwindigkeitsvektors $\\mathbf{v}$.\n",
      "Da der Druck eine skalare Gr\u00f6\u00dfe ist, gibt es in `indexp` nur die Nummer `i` des Freiheitsgrads auf dem Dreieck.\n",
      "\n",
      "Die Freiheitsgrade sollen hierbei wie folgt nummeriert werden:\n",
      "zun\u00e4chst die zur Geschwindigkeit geh\u00f6rigen Freiheitsgrade, dann die zum Druck geh\u00f6rigen Freiheitsgrade.\n",
      "Erstere sollen zuerst nach den zugeh\u00f6rigen Freiheitsgraden in $V$, dann nach der Komponente von $\\mathbf{v}$ sortiert werden.\n",
      "F\u00fcr Dreiecke mit Lagrange-Elementen zweiter Ordnung f\u00fcr $\\mathbf{v}$ und erster Ordnung f\u00fcr $p$ sollen die Freiheitsgrade also in der Reihenfolge\n",
      "\n",
      "  $v_0^0, v_1^0, v_0^1, v_1^1, \\dots, v_0^5, v_1^5, p^0, p^1, p^2$\n",
      "  \n",
      "nummeriert werden.\n",
      "\n",
      "Implementieren Sie dann die Methode `assembleLocal(self, element, localV, localW)`.\n",
      "Die Argumente `localV` und `localW` sind hierbei die lokalen Formfunktionen der Funktionsr\u00e4ume V bzw. W.\n",
      "\n",
      "Die Methode soll ein Paar bestehend aus lokaler Matrix `localA` und lokaler rechten Seite `localb` zur\u00fcckgeben.\n",
      "Verwenden Sie `self.mu` als Koeffizienten $\\mu$ und `self.f` als Kraftdichte $f: \\Omega \\to \\mathbb{R}^2$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "class LocalAssembler:\n",
      "    def __init__(self, dimension, mu, f):\n",
      "        self.dimension = dimension\n",
      "        self.mu = mu\n",
      "        self.f = f\n",
      "\n",
      "    def indexv(self, i, k):\n",
      "        pass\n",
      "    \n",
      "    def indexp(self, i):\n",
      "        pass\n",
      "\n",
      "    def assembleLocal(self, element, localV, localW):\n",
      "        pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 3\n",
      "\n",
      "Implementieren Sie eine Klasse `GlobalAssembler`.\n",
      "Diese soll die Systemmatrix `A` und rechte Seite `b` aufstellen.\n",
      "\n",
      "Implementieren Sie dazu zun\u00e4chst die Methoden `indexv(self, i, k)` und `indexp(self, i)`.\n",
      "Diese sollen die globalen Freiheitsgrade analog zur Nummerierung im `LocalAssembler` nummerieren, allerdings hier mit der *globalen* Nummer `i` des Freiheitsgrads.\n",
      "\n",
      "Implementieren Sie dann die Methode `assemble(self, v0, gammav, gammap)`.\n",
      "Die Funktion $\\mathbf{v}_0: \\Omega \\to \\mathbb{R^d}$ stellt die vorgegebenen Werte auf dem Dirichlet-Rand dar und ist auf das ganze Gebiet $\\Omega$ fortgesetzt.\n",
      "Der Dirichlet-Rand selbst ist durch Boolsche Funktionen `gammav` bzw. `gammap` gegeben.\n",
      "\n",
      "Die Methode soll ein Tupel $(A, b)$ bestehend aus Systemmatrix $A$ und rechter Seite $b$ zur\u00fcckgeben."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class GlobalAssembler:\n",
      "    def __init__(self, grid, V, W, localAssembler):\n",
      "        self.grid = grid\n",
      "        self.V = V\n",
      "        self.W = W\n",
      "        self.localAssembler = localAssembler\n",
      "        \n",
      "    def indexv(self, i, k):\n",
      "        pass\n",
      "    \n",
      "    def indexp(self, i):\n",
      "        pass\n",
      "\n",
      "    def assemble(self, v0, gammav, p0, gammap):\n",
      "        pass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 4\n",
      "\n",
      "Rufen Sie den `GlobalAssembler` f\u00fcr das am Anfang beschriebene Stokes-Problem auf und l\u00f6sen Sie das lineare Gleichungssystem $Ax = b$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 64\n",
      "f = lambda x: np.array((0.0, 0.0))\n",
      "\n",
      "v0 = lambda x: np.array((1.0 if x[1] > 0.9999 else 0.0, 0.0))\n",
      "gammav = lambda x: x[0] < 1e-5 or x[0] > 0.9999 or x[1] < 1e-5 or x[1] > 0.9999\n",
      "\n",
      "p0 = lambda x: 0.0\n",
      "gammap = lambda x: x[1] < 1e-5\n",
      "\n",
      "if 'grid' not in globals():\n",
      "    grid = dune.grid.ugGrid(\n",
      "        (dune.common.reader.structured, \"simplex\", [0.0, 0.0], [1.0, 1.0], [n, n]),\n",
      "        dimgrid=2\n",
      "    )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# ...\n",
      "# x = ..."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 4\n",
      "\n",
      "Betrachten Sie die Ausgabe in ParaView.\n",
      "\n",
      "Achtung: der `VTKWriter` schreibt das Ergebnis als drei separate skalare Felder.\n",
      "Um das Geschwindigkeitsfeld zu betrachten, wandeln Sie daher die skalaren Funktionen `v0` und `v1` zu einer vektorwertigen Funktion um, indem Sie in der Werkzeugleiste den \"Calculator\" ausw\u00e4hlen und als Ausdruck `v0 * iHat + v1 * jHat` eingeben.\n",
      "F\u00fcgen Sie dann einen Glyph-Filter hinzu.\n",
      "\n",
      "Sie k\u00f6nnen sich auch die Str\u00f6mungslinien anzeigen lassen.\n",
      "F\u00fcgen Sie dazu einen \"Stream Tracer\"-Filter hinzu und w\u00e4hlen als Seed Type die \"High Resolution Line Source\" aus."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vtkWriter = grid.vtkWriter()\n",
      "\n",
      "if 'x' in globals():\n",
      "    # Limitierung in dune-python; schreibe alle Komponenten von v als skalare Funktion\n",
      "    for k in range(grid.dimension):\n",
      "        xk = x[k : grid.dimension*len(V)+k : grid.dimension]\n",
      "        vk = V.asFunction(xk)\n",
      "        vk.addToVTKWriter(\"v{}\".format(k), vtkWriter, dune.common.DataType.PointData)\n",
      "\n",
      "    xp = x[grid.dimension*len(V) :]\n",
      "    p = W.asFunction(xp)\n",
      "    p.addToVTKWriter(\"p\", vtkWriter, dune.common.DataType.PointData)\n",
      "\n",
      "    vtkWriter.write(\"05-Stokes\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}