{
 "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.5.4"
  },
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Implementierung einer Finite-Elemente-Funktion\n",
      "\n",
      "In dieser \u00dcbung sollen Sie sich mit DUNE bekannt machen und selbst eine Finite-Elemente-Funktion implementieren.\n",
      "Weiter sollen gegebene Funktionen durch Finite-Elemente-Funktionen interpoliert werden und der Interpolationsfehler berechnet werden.\n",
      "\n",
      "Informationen zu DUNE finden Sie im Einf\u00fchrungs-Notebook.\n",
      "\n",
      "**Hinweis:**\n",
      "Falls Sie die Aufgaben in der bereitgestellen VM bearbeiten, m\u00fcssen Sie zun\u00e4chst die aktuell verf\u00fcgbaren Updates installieren.\n",
      "Geben Sie dazu die Befehle \"`sudo apt update`\" und \"`sudo apt dist-upgrade`\" im Terminal ein; geben Sie als Passwort \"`asdf`\" ein."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import dune.common\n",
      "import dune.geometry\n",
      "import dune.grid"
     ],
     "language": "python",
     "metadata": {
      "scrolled": true
     },
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 1: Integrale berechnen\n",
      "\n",
      "Auf einem $10\\times10$-Gitter f\u00fcr $\\Omega = [0, 1]^2$ soll die Funktion $f(x) = x_0^2 + x_1^2$ integriert werden.\n",
      "Wir benutzen dazu auf jedem Element eine Quadraturformel.\n",
      "\n",
      "Schreiben Sie eine Funktion `integrate(grid, f)`, die das Integral der Funktion `f(x)` auf dem Gitter `grid` mit einer Quadraturregel der Ordnung 2 berechnet.\n",
      "\n",
      "**Hinweis:**\n",
      "Eine Beschreibung, wie Sie Quadraturregeln erhalten, findet sich in der aktuellen Version des Einf\u00fchrungs-Notebooks."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def integrate(grid, f):\n",
      "    pass\n",
      "    # Schreiben Sie hier Ihren Programmcode"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Im Folgenden wird ein Gitter `grid` f\u00fcr $\\Omega = [0,1]^2$ angelegt, die Funktion $f$ definiert und das Integral von $f$ berechnet:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "grid = dune.grid.ugGrid(\n",
      "    (dune.common.reader.structured, \"simplex\", [0.0, 0.0], [1.0, 1.0], [10, 10]),\n",
      "    dimgrid=2\n",
      ")\n",
      "f = lambda x: x[0]**2 + x[1]**2\n",
      "\n",
      "print(\"int(f(x)) = {}\".format(integrate(grid, f)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 2: Lokale Basisfunktionen\n",
      "\n",
      "In dieser Aufgabe sollen Sie die P1-Formfunktionen auf dem Referenzdreieck mit den Eckpunkten $x_0 = (0,0)$, $x_1 = (1,0)$ und $x_2 = (0,1)$ implementieren.\n",
      "\n",
      "Die Basisfunktionen $\\phi_i$ sind affin-linear mit der Eigenschaft $\\phi_i(x_j) = \\delta_{ij}$ f\u00fcr die Eckpunkte $x_j$ des Dreiecks.\n",
      "\n",
      "Schreiben Sie eine Klasse `P1LocalBasis` mit den Methoden `evaluateFunction(self, xlocal)` und `evaluateJacobian(self, xlocal)`.\n",
      "Die Methode `evaluateFunction` soll ein numpy-Array mit den Werten aller P1-Basisfunktionen auf dem Dreieck im Punkt `xlocal` zur\u00fcckgeben.\n",
      "`evaluateJacobian` soll ein zweidimensionales numpy-Array mit den Gradienten aller P1-Basisfunktionen auf dem Dreieck zur\u00fcckgeben, so dass `evaluateJacobian(self, xlocal)[i]` der Gradient der `i`-ten Basisfunktion ist.\n",
      "\n",
      "![Nummerierung der Unterentit\u00e4ten f\u00fcr das Referenzdreieck](https://www.dune-project.org/doxygen/master/gg_triangle.png)\n",
      "\n",
      "**Hinweis:**\n",
      "Im Einf\u00fchrungs-Notebook ist ein kurzer Abschnitt zu objektorientierter Programmierung in Python enthalten, der Ihnen zeigt wie Sie Klassen und Methoden implementieren."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Schreiben Sie hier Ihren Programmcode"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\u00dcberpr\u00fcfen Sie, ob sich die Funktionen wie erwartet verhalten, indem Sie `evaluateFunction` und `evaluateJacobian` in den drei Eckpunkten des Referenzdreiecks aufrufen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Schreiben Sie hier Ihren Programmcode"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 3: Diskrete Gitterfunktionen\n",
      "\n",
      "In dieser Aufgabe sollen Sie eine Klasse `P1Function` implementieren, die eine stetige, st\u00fcckweise affine Funktion auf dem Gitter darstellt.\n",
      "Sei `u` der Koeffizientenvektor einer solchen Gitterfunktion bez\u00fcglich der Knotenbasis.\n",
      "Die Klasse `P1Function` soll den Koeffizientenvektor `u` und das Gitter als Attribute halten und sich in lokalen Koordinaten bez\u00fcglich einzelner Elemente auswerten lassen.\n",
      "\n",
      "Legen Sie eine Klasse `P1Function` an und implementieren Sie die Methoden `__init__(self, grid, u)`, `__call__(self, element, xlocal)` und `gradient(self, element, xlocal)`.\n",
      "Die Methode `__init__` soll hierbei die \u00fcbergebenen Argumente -- das Gitter `grid` und den Koeffizientenvektor `u` -- als Attribute speichern.\n",
      "Weiterhin soll eine Instanz der Klasse `P1LocalBasis` als Attribut gespeichert werden.\n",
      "\n",
      "`__call__` soll die Funktion auf dem Element `element` im Punkt `xlocal` auswerten; verwenden Sie hierzu die `P1LocalBasis`.\n",
      "Die Methode `gradient` soll entsprechend die Ableitung berechnen.\n",
      "Beachten Sie, dass `P1LocalBasis` die Ableitung bez\u00fcglich der Koordinaten auf dem Referenzelement zur\u00fcckgibt; die notwendigen Informationen f\u00fcr die Transformation erhalten Sie aus `element.geometry`.\n",
      "\n",
      "**Hinweis:**\n",
      "Verwenden Sie das `IndexSet` um die f\u00fcr das Element `element` relevanten Eintr\u00e4ge aus `u` zu extrahieren."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Schreiben Sie hier Ihren Programmcode."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 4: Interpolation\n",
      "\n",
      "Sei `f` eine beliebige Pythonfunktion, die eine Funktion $\\Omega \\to \\mathbb{R}$ implementiert.\n",
      "\n",
      "Schreiben Sie eine Funktion `interpolate(grid, f)`, die die Funktion `f` im P1-Funktionenraum interpoliert und den zugeh\u00f6rigen Koeffizientenvektor als numpy-Vektor zur\u00fcckgibt.\n",
      "\n",
      "**Hinweis:**\n",
      "Mit `grid.size(grid.dimension)` erhalten Sie die Anzahl der Knotenpunkte im Gitter.\n",
      "\n",
      "**Hinweis:**\n",
      "Mit einer `for`-Schleife k\u00f6nnen Sie nicht nur \u00fcber Elemente des Gitters iterieren, sondern auch \u00fcber Knotenpunkte."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def interpolate(grid, f):\n",
      "    pass\n",
      "    # Schreiben Sie hier Ihren Programmcode"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "grid = dune.grid.ugGrid(\n",
      "    (dune.common.reader.structured, \"simplex\", [0.0, 0.0], [1.0, 1.0], [10, 10]),\n",
      "    dimgrid=2\n",
      ")\n",
      "\n",
      "f = lambda x: x[0]**2 + x[1]**2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Interpolieren Sie nun $f$ und weisen $u$ den resultierenden Koeffizientenvektor zu.\n",
      "Erzeugen Sie dann eine `P1Function` mit Namen `fh` mit den Koeffizienten $u$.\n",
      "\n",
      "Schreiben Sie eine Funktion `integrateLocal(grid, flocal)` um die interpolierte Funktion \u00fcber das Gitter zu integrieren."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fh = None\n",
      "# Schreiben Sie hier Ihren Programmcode"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Aktuell kann zu jedem Zeitpunkt nur ein UG-Gitter existieren.\n",
      "Bevor ein Gitter in der n\u00e4chsten Aufgabe wieder angelegt werden kann, muss das aktuelle Gitter `grid` daher freigegeben werden, indem `grid` und `fh` explizit gel\u00f6scht werden (`fh` muss gel\u00f6scht werden, da `P1Function` eine weitere Referenz auf das Gitter enth\u00e4lt)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "del grid\n",
      "del fh"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 5: Berechnung des Interpolationsfehlers\n",
      "\n",
      "In dieser Aufgabe soll der Interpolationsfehler in der $L^2$- und $H^1$-Norm f\u00fcr verschiedene Gitterweiten berechnet werden.\n",
      "\n",
      "Interpolieren Sie die gegebene Funktion $f(x) = x_0^2 + x_1^2$ auf einem $n\\times n$-Dreiecksgitter mit $n = 10, 20, 40, 80$.\n",
      "Berechnen Sie jeweils den Interpolationsfehler $\\| f - f_h \\|_{L^2(\\Omega)}$ und $\\| f - f_h \\|_{H^1(\\Omega)}$.\n",
      "\n",
      "Wie verh\u00e4lt sich der Fehler als Funktion von $h = 1/n$?\n",
      "\n",
      "**Hinweis:**\n",
      "Berechnen Sie zun\u00e4chst $\\|f-f_h\\|_{\\Omega}^2 = \\sum_T \\|f-f_h\\|_{T}^2$ f\u00fcr die jeweilige Norm."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "f = lambda x: x[0]**2 + x[1]**2\n",
      "gradf = lambda x: (2*x[0], 2*x[1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Schreiben Sie hier Ihren Programmcode."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Aufgabe 6: Ausgabe als VTK-Datei (optional)\n",
      "\n",
      "Mit dem hier angegebenen Code k\u00f6nnen Sie eine beliebige `P1Function` als VTK-Datei ausgeben.\n",
      "Definieren Sie dazu `fh` als eine `P1Function` Ihrer Wahl und f\u00fchren den Code aus.\n",
      "\n",
      "Die Ausgabedatei k\u00f6nnen Sie mit dem Programm `paraview` betrachten."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "grid = dune.grid.ugGrid(\n",
      "    (dune.common.reader.structured, \"simplex\", [0.0, 0.0], [1.0, 1.0], [10, 10]),\n",
      "    dimgrid=2\n",
      ")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Definieren Sie hier eine `P1Function`\n",
      "fh = None"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "if fh is not None:\n",
      "    gridfunction = grid.localGridFunction(lambda e, x: [fh(e, x)])\n",
      "    grid.writeVTK(\"02-fe-basis\", pointdata={\"f\": gridfunction})\n",
      "\n",
      "del grid\n",
      "del fh"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": null
    }
   ],
   "metadata": {}
  }
 ]
}