{
 "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": [
      "# L\u00f6sen der Poisson-Gleichung mit dune-functions\n",
      "\n",
      "Bisher wurden alle Finite-Elemente-Funktionen von Hand implementiert.\n",
      "Dies hatte didaktische Gr\u00fcnde --- man versteht solche Funktionen besser wenn man sie einmal selbst implementiert hat.\n",
      "\n",
      "Das DUNE-System enth\u00e4lt eine Sammlung von fertigen FE-Funktionen im Modul `dune-functions`.\n",
      "Dieser Text zeigt, wie diese Funktionen verwendet werden.\n",
      "Als Beispiel wird der L\u00f6ser f\u00fcr die Poisson-Gleichung vom letzten \u00dcbungszettel damit neu implementiert.\n",
      "Das Programm folgt dabei der gleichen Struktur.\n",
      "Die Hauptarbeit wird wieder von zwei Klassen `LocalAssembler` und `GlobalAssembler` erledigt."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Vorspann, Problembeschreibung und Anlegen des Gitters\n",
      "\n",
      "Zun\u00e4chst werden die ben\u00f6tigen Module geladen.\n",
      "Neu ist hier das Modul `dune.functions`."
     ]
    },
    {
     "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": [
      "Es folgt die Problembeschreibung mit Gitteraufl\u00f6sung, rechter Seite $f$, Dirichlet-Randwerten $u_0$ und Boolscher Funktion $\\gamma$ f\u00fcr den Teil des Randes, auf dem Dirichlet-Randbedingung gefordert wird:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 32\n",
      "f = lambda x: 20 * np.pi**2 * np.sin(2*np.pi*x[0]) * np.cos(4*np.pi*x[1])\n",
      "u0 = lambda x: np.sin(2 * np.pi * x[0])\n",
      "gamma = lambda x: x[0] < 1e-5 or x[0] > 0.9999 or x[1] < 1e-5 or x[1] > 0.9999"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Lege ein Gitter f\u00fcr das Gebiet $\\Omega = (0,1)^2$ an:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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": "markdown",
     "metadata": {},
     "source": [
      "## Anlegen einer dune-functions-Basis\n",
      "\n",
      "Nun kann eine globale Basis aus dune-functions angelegt werden.\n",
      "Hierzu wird das Gitter und der gew\u00fcnschte Funktionenraum ben\u00f6tigt, hier zum Beispiel eine Lagrange-Basis zweiter Ordnung:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "basis = dune.functions.defaultGlobalBasis(grid, dune.functions.Lagrange(order=2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Im `GlobalAssembler` wird die globale Basis auf einzelne Elemente eingeschr\u00e4nkt werden (`LocalView`) und die Formfunktionen der Einschr\u00e4nkung an den `LocalAssembler` \u00fcbergeben werden; diese Formfunktionen entsprechen der `P1LocalBasis` der \u00dcbungen und besitzen die gleichen Methoden `evaluateFunction` und `evaluateJacobian`.\n",
      "\n",
      "Im `GlobalAssembler` wird weiterhin die Zuordnung von lokalen Freiheitsgraden zu globalen Freiheitsgraden ben\u00f6tigt.\n",
      "Auch diese wird von der Einschr\u00e4nkung bereit gestellt; siehe die Beschreibung des `GlobalAssembler` unten.\n",
      "\n",
      "F\u00fcr die Ausgabe kann eine Basis zusammen mit einem Koeffizientenvektor als Funktion interpretiert werden.\n",
      "Dies entspricht der `P1Function`-Klasse aus den \u00dcbungen.\n",
      "Der Programmcode hierzu findet sich im letzten Abschnitt."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## `LocalAssembler`\n",
      "\n",
      "Im `LocalAssembler` wird der Anteil eines einzelnen Elements an der Systemmatrix und rechten Seite berechnet.\n",
      "Dies erledigt die Methode `assembleLocal` mit den Argumenten `self`, `element` und `localBasis`.\n",
      "\n",
      "Die `localBasis` verh\u00e4lt sich mit den Methoden `evaluateFunction` und `evaluateJacobian` wie die Klasse `P1LocalBasis` der \u00dcbungen.\n",
      "Die Anzahl der Freiheitsgrade wird jetzt aus der `localBasis` mit `len(localBasis)` ermittelt; diese Funktion fehlte in der `P1LocalBasis` noch."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "class LocalAssembler:\n",
      "    def assembleLocal(self, element, localBasis):\n",
      "        n = len(localBasis)\n",
      "\n",
      "        localA = np.zeros((n, n))\n",
      "        localb = np.zeros(n)\n",
      "\n",
      "        geometry = element.geometry\n",
      "        for pt in grid._module.quadratureRule(element.type, 2):\n",
      "            x = geometry.position(pt.position)\n",
      "            jit = geometry.jacobianInverseTransposed(pt.position)\n",
      "\n",
      "            phi = localBasis.evaluateFunction(pt.position)\n",
      "            gradPhiRef = localBasis.evaluateJacobian(pt.position)\n",
      "            gradPhi = [ np.dot(jit, np.array(g)[0]) for g in gradPhiRef ]\n",
      "\n",
      "            w = pt.weight * geometry.integrationElement(pt.position)\n",
      "        \n",
      "            for i in range(n):\n",
      "                for j in range(n):\n",
      "                    localA[i, j] += w*np.dot(gradPhi[i], gradPhi[j])\n",
      "                localb[i] += w*phi[i][0]*f(x)\n",
      "        \n",
      "        return localA, localb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## `GlobalAssembler`\n",
      "\n",
      "Im `GlobalAssembler` wird der `LocalAssembler` f\u00fcr jedes Element aufgerufen und die Systemmatrix aus den einzelnen Teilen zusammengesetzt.\n",
      "Hier werden auch die Dirichlet-Daten direkt ber\u00fccksichtigt statt nach der Assemblierung von einer separaten Funktion gesetzt zu werden.\n",
      "\n",
      "Um die Dirichlet-Werte zu setzen, muss eine gegebene Funktion interpoliert werden.\n",
      "In der \u00dcbung gab es dazu die Funktion `interpolate`.\n",
      "Dies ist nun eine Methode der globalen Basis.\n",
      "\n",
      "Mit der neuen `interpolate`-Methode ist es auch m\u00f6glich eine Boolsche Funktion $\\gamma: \\Omega \\to \\{\\text{True}, \\text{False}\\}$ zu interpolieren.\n",
      "Dies gibt einen Vektor von Wahrheitswerten zur\u00fcck, wobei jeder Eintrag besagt, ob $\\gamma$ an der dem Freiheitsgrad zugeordneten Position wahr oder falsch ist.\n",
      "Dies wird verwendet, um den Dirichlet-Rand zu markieren, der nicht immer dem gesamten Gebietrand $\\partial\\Omega$ entspricht.\n",
      "\n",
      "Der `LocalView` einer globalen Basis ist die Einschr\u00e4nkung der Basis auf ein einzelnes Element.\n",
      "Im Programm wird dies daf\u00fcr benutzt, die lokalen Formfunktionen zu erhalten, die an den `LocalAssembler` \u00fcbergeben wird.\n",
      "\n",
      "Der `LocalView` wird mit `localView = basis.localView()` angelegt.\n",
      "Bevor man ihn verwendet, muss man allerdings noch das Element angeben, auf das eingeschr\u00e4nkt werden soll.\n",
      "Dies erledigt der Aufruf `localView.bind(element)`.\n",
      "Dann erh\u00e4lt man die Formfunktionen mit `localView.tree().finiteElement.localBasis`.\n",
      "\n",
      "Das `LocalIndexSet` einer Basis liefert die Zuordung von lokalen Freiheitsgraden auf einem Element zu den zugeh\u00f6rigen globalen Freiheitsgraden.\n",
      "In den \u00dcbungen wurde dies bisher fest einkodiert; durch die zus\u00e4tzliche Abstraktion funktioniert das neue Programm z.B. auch mit Lagrange-Basen beliebiger Ordnung.\n",
      "Um die globale Nummer des Freiheitsgrads `i` auf dem Element zu erhalten, wird `localIndexSet.index(i)[0]` aufgerufen.\n",
      "Die Indizierung `[0]` ist notwendig, da in dune-functions Freiheitsgrade im Allgemeinen mit Tupeln von Zahlen nummeriert werden, was hier aber nicht verwendet werden soll.\n",
      "\n",
      "Auch dem `LocalIndexSet` muss man nach dem Anlegen noch mitteilen, welche Einschr\u00e4nkung betrachtet wird.\n",
      "Das `LocalIndexSet` wird daher vor der Verwendung mit `localIndexSet.bind(localView)` an eine Einschr\u00e4nkung gebunden; der `localView` muss hierbei bereits an ein Element gebunden sein.\n",
      "\n",
      "Nach Verwendung wird die Bindung des `LocalIndexSet` und `LocalView` mit der Methode `unbind()` wieder aufgehoben."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class GlobalAssembler:\n",
      "    def __init__(self, basis, localAssembler):\n",
      "        self.basis = basis\n",
      "        self.localAssembler = localAssembler\n",
      "\n",
      "    def assemble(self, u0, gamma):\n",
      "        N = len(basis)\n",
      "        A = lil_matrix((N, N))\n",
      "        # `constraints` wird als Vektor von Wahrheitswerten angelegt\n",
      "        constraints = np.ndarray(N, dtype=np.bool_)\n",
      "        b = np.zeros(N)\n",
      "        x = np.ndarray(N)\n",
      "\n",
      "        self.basis.interpolate(x, u0)\n",
      "        self.basis.interpolate(constraints, gamma)\n",
      "\n",
      "        localView = basis.localView()\n",
      "        localIndexSet = basis.localIndexSet()\n",
      "\n",
      "        for element in grid.elements():\n",
      "            localView.bind(element)\n",
      "            localIndexSet.bind(localView)\n",
      "\n",
      "            # localBasis: Formfunktionen auf einem einzelnen Element\n",
      "            localBasis = localView.tree().finiteElement.localBasis\n",
      "\n",
      "            # localN: Anzahl der Freiheitsgrade auf dem aktuellen Element\n",
      "            localN = len(localBasis)\n",
      "\n",
      "            localA, localb = self.localAssembler.assembleLocal(element, localBasis)\n",
      "\n",
      "            for i in range(localN):\n",
      "                # gi: Globaler Freiheitsgrad\n",
      "                gi = localIndexSet.index(i)[0]\n",
      "                # Dirichlet-Randbedingung, falls \\gamma am Freiheitsgrad wahr ist:\n",
      "                if constraints[gi]:\n",
      "                    A[gi, gi] = 1\n",
      "                    b[gi] = x[gi]\n",
      "                else:\n",
      "                    b[gi] += localb[i]\n",
      "                    for j in range(localN):\n",
      "                        gj = localIndexSet.index(j)[0]\n",
      "                        A[gi, gj] += localA[i, j]\n",
      "\n",
      "            localIndexSet.unbind()\n",
      "            localView.unbind()\n",
      "            \n",
      "        return A, b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Assemblierung, L\u00f6sen und Ausgabe als VTK-Datei\n",
      "\n",
      "Der globale Assemblierer erh\u00e4lt jetzt die Basis als Argument."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "localAssembler = LocalAssembler()\n",
      "globalAssembler = GlobalAssembler(basis, localAssembler)\n",
      "\n",
      "A, b = globalAssembler.assemble(u0, gamma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "Acsr = A.tocsr()\n",
      "x = spsolve(Acsr, b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Die Methode `basis.asFunction(x)` liefert die FE-Funktion zum Koeffizientenvektor `x`.\n",
      "Diese kann direkt an den `VTKWriter` \u00fcbergeben werden."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u = basis.asFunction(x)\n",
      "vtkWriter = grid.vtkWriter()\n",
      "u.addToVTKWriter(\"u\", vtkWriter, dune.common.DataType.PointData)\n",
      "vtkWriter.write(\"05-Poisson\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}