{ "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": [ "# Poisson-Problem mit adaptiver Gitterverfeinerung\n", "\n", "In dieser \u00dcbung soll das Poisson-Problem auf einem Gebiet mit einspringender Ecke (\"L-shaped domain\") mit adaptiver Gitterverfeinerung gel\u00f6st werden.\n", "Dazu wird der Residuenfehlersch\u00e4tzer aus der Vorlesung implementiert, der eine Sch\u00e4tzung $\\eta_T$ f\u00fcr den Fehler eines einzelnen Elements $T$ berechnet.\n", "Die Elemente mit dem gr\u00f6\u00dften gesch\u00e4tzten Fehler werden verfeinert und der L\u00f6ser ein weiteres Mal aufgerufen.\n", "\n", "Wir betrachten\n", "\\begin{equation}\n", " -\\Delta u = 0\n", "\\end{equation}\n", "auf dem Gebiet $\\Omega = [-1,1] \\times [-1,1] \\setminus (0,1] \\times [-1,0)$ (ein Quadrat ohne die untere rechte Ecke).\n", "Als Dirichlet-Randbedingung w\u00e4hlen wir\n", "\\begin{equation}\n", " u(r, \\theta) = r^{\\frac23} \\sin(\\frac23 \\theta)\n", "\\end{equation}\n", "auf $\\partial\\Omega$ analog zu Aufgabe 2 der zweiten \u00dcbung.\n", "(Wir verwenden die Polarkoordinaten $(r, \\theta)$ nur f\u00fcr die Problembeschreibung, die Rechnung selbst erfolgt in kartesischen Koordinaten $(x, y)$.)\n", "\n", "*Hinweis:*\n", "Laden Sie die Datei `06-ldomain.msh` herunter und speichern Sie die Datei im gleichen Verzeichnis wie dieses Notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vorspann, Problembeschreibung und Anlegen des Gitters\n", "\n", "Zun\u00e4chst werden die ben\u00f6tigen Module geladen." ] }, { "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 rechter Seite $f$, Dirichlet-Randwerten $u_0$ und Boolscher Funktion $\\gamma$ f\u00fcr den Teil des Gebietes, auf dem Dirichlet-Randbedingung gefordert wird:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "f = lambda x: 0.0\n", "def u0(x):\n", " phi = np.arctan2(x[1], x[0])\n", " if phi < 0.0:\n", " phi += 2*np.pi\n", " return (x[0]**2 + x[1]**2)**(1./3) * np.sin(2./3 * phi)\n", "gamma = lambda x: x[0] < -0.9999 or x[0] > 0.9999 or x[1] < -0.9999 or x[1] > 0.9999 or (x[0] > -0.00001 and x[1] < 0.00001)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das Gitter wird diesmal aus der Datei `06-ldomain.msh` geladen.\n", "Dies erledigt der folgende Code:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "if 'grid' not in globals():\n", " grid = dune.grid.ugGrid(\n", " (dune.common.reader.gmsh, \"06-ldomain.msh\"),\n", " dimgrid=2\n", " )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Au\u00dferdem werden die Quadraturregeln unter einem einfacheren Namen verf\u00fcgbar gemacht:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "quadratureRule = grid._module.quadratureRule" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## L\u00f6ser f\u00fcr das Poisson-Problem\n", "\n", "Die Klassen `LocalAssembler` und `GlobalAssembler` sind im Wesentlichen aus dem vorherigen Notebook `05-Poisson-with-dune-functions.ipynb` \u00fcbernommen.\n", "Sie k\u00f6nnen alternativ auch Ihren Programmcode aus der vierten \u00dcbung verwenden." ] }, { "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 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": "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(self.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 = self.basis.localView()\n", " localIndexSet = self.basis.localIndexSet()\n", "\n", " for element in self.basis.gridView.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": [ "Die Funktion `solve(basis)` l\u00f6st das Poisson-Problem mit Ansatzfunktionen aus der gegebenen Basis und liefert den Koeffizientenvektor der L\u00f6sung zur\u00fcck." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def solve(basis):\n", " localAssembler = LocalAssembler()\n", " globalAssembler = GlobalAssembler(basis, localAssembler)\n", "\n", " A, b = globalAssembler.assemble(u0, gamma)\n", " return spsolve(A.tocsr(), b)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `P1Function` f\u00fcr Auswertung der Ableitung einer P1-Funktion\n", "\n", "Der Python-Anbindung von `dune-functions` fehlt noch die Auswertung der Ableitung einer Gitterfunktion.\n", "Daher wird in dieser \u00dcbung die Klasse `P1Function` bereitgestellt, deren Methode `gradient(self, element, xlocal)` die Auswertung der Ableitung erm\u00f6glicht:\n", "\n", "```\n", " u = P1Function(basis, x)\n", " du = u.gradient(element, xlocal)\n", "```\n", "Dabei ist `basis` eine `dune-functions`-Basis, `x` ein Koeffizientenvektor, `element` ein Gitterelement und `xlocal` lokale Koordinaten in diesem Element." ] }, { "cell_type": "code", "collapsed": false, "input": [ "class P1Function:\n", " def __init__(self, basis, u):\n", " self._basis = basis\n", " self._localView = basis.localView()\n", " self._grid = basis.gridView\n", " self._u = u\n", " def _index(self, element, i):\n", " return self._grid.indexSet.subIndex(element, i, self._grid.dimension)\n", " def __call__(self, element, xlocal):\n", " idx = [self._index(element, i) for i in range(3)]\n", " self._localView.bind(element)\n", " phi = self._localView.tree().finiteElement.localBasis.evaluateFunction(xlocal)\n", " self._localView.unbind()\n", " return np.dot(phi, self._u[idx])\n", " def gradient(self, element, xlocal):\n", " idx = [self._index(element, i) for i in range(3)]\n", " JinvT = element.geometry.jacobianInverseTransposed(xlocal)\n", " self._localView.bind(element)\n", " refGradPhi = self._localView.tree().finiteElement.localBasis.evaluateJacobian(xlocal)\n", " self._localView.unbind()\n", " refGrad = np.tensordot(refGradPhi, self._u[idx], axes=[(0,), (0,)])\n", " grad = [ np.dot(JinvT, g) for g in refGrad ]\n", " return np.array(grad)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Intersections\n", "\n", "Um den Residuensch\u00e4tzer aus der Vorlesung implementieren zu k\u00f6nnen, brauchen Sie Zugriff auf die Elementkanten eines Gitters.\n", "In DUNE hei\u00dfen diese Kanten `Intersection`s.\n", "\n", "F\u00fcr ein Element `element` k\u00f6nnen Sie mit `for intersection in grid.intersections(element)` \u00fcber die Kanten des Elements iterieren.\n", "\n", "Wie f\u00fcr Elemente k\u00f6nnen Sie mit `intersection.type` den Geometrietyp der Intersection erhalten, wie etwa f\u00fcr Quadraturregeln ben\u00f6tigt wird.\n", "\n", "Weiterhin k\u00f6nnen Sie \u00fcber `intersection.outside` auf das Nachbarelement zugreifen, sofern dieses existiert.\n", "F\u00fcr eine Kante am Gebietsrand liefert `intersection.outside` den Wert `None` zur\u00fcck.\n", "\n", "\u00dcber `intersection.integrationOuterNormal(xlocal)` erhalten Sie die \u00e4u\u00dfere Normale am Punkt `xlocal`, die bereits mit dem Integrationselement skaliert ist.\n", "\n", "Das `intersection`-Objekt verf\u00fcgt wie die Elemente \u00fcber ein Geometrieobjekt mit den Methoden `position(xlocal)` und `positionLocal(x)`, die lokale Koordinaten auf der Kante in globale Koordinaten umwandelt und umgekehrt.\n", "Die Eigenschaft `intersection.geometry.volume` liefert das \"Volumen\" der Intersection zur\u00fcck, d.h. die Kantenl\u00e4nge $h_e$.\n", "\n", "Um lokale Koordinaten `xlocal` der Kante in lokale Koordinaten des inneren bzw. \u00e4u\u00dferen Elements umzurechnen, verketten Sie `intersection.position` mit `element.geometry.localPosition` bzw. `intersection.outside.geometry.localPosition`.\n", "\n", "### Beispiel\n", "\n", "Das folgende Beispiel z\u00e4hlt die Kanten am Gebietsrand und im Inneren." ] }, { "cell_type": "code", "collapsed": false, "input": [ "nBoundaryEdges = 0\n", "nInteriorEdges = 0\n", "\n", "for element in grid.elements():\n", " for intersection in grid.intersections(element):\n", " if intersection.outside is None:\n", " nBoundaryEdges += 1\n", " else:\n", " nInteriorEdges += 1\n", " \n", "print(\"Das Gitter hat {} Kanten am Rand und {} Kanten im Inneren.\".format(nBoundaryEdges, nInteriorEdges))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Beispiel\n", "\n", "Das folgende Beispiel gibt f\u00fcr die erste Kante des ersten Elements die Koordinaten des Kantenmittelpunkts in globalen Koordinaten und lokalen Koordinaten des inneren Elements und, falls die Kante nicht zum Rand geh\u00f6rt, des \u00e4u\u00dferen Elements zur\u00fcck." ] }, { "cell_type": "code", "collapsed": false, "input": [ "firstElement = next(iter(grid.elements()))\n", "firstIntersection = next(iter(grid.intersections(firstElement)))\n", "intersectionGeometry = firstIntersection.geometry\n", "\n", "centerGlobal = intersectionGeometry.center\n", "print(\"Globale Koordinaten des Mittelpunkts: {}\".format(centerGlobal))\n", "\n", "centerInside = firstElement.geometry.localPosition(centerGlobal)\n", "print(\"In lokalen Koordinaten des Elements: {}\".format(centerInside))\n", "\n", "outside = firstIntersection.outside\n", "if outside is not None:\n", " centerOutside = outside.geometry.localPosition(centerGlobal)\n", " print(\"In lokalen Koordinaten des Nachbarelements: {}\".format(centerOutside))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 1 (Fehlersch\u00e4tzer)\n", "\n", "Implementieren Sie die Funktion `computeError(basis, x, f)`, die den Fehlersch\u00e4tzer bereitstellt.\n", "Hierbei ist `x` der Koeffizientenvektor der diskreten L\u00f6sung und `f(x)` die rechte Seite.\n", "Die Funktion soll ein Paar $(\\eta, (\\eta_T)_{T\\in\\mathcal{T}})$ mit dem globalen Fehler $\\eta$ und den Elementfehlern $\\eta_T$ zur\u00fcckgeben.\n", "Die lokalen Fehler $\\eta_T$ sollen hierbei als numpy-Array gespeichert werden; verwenden Sie `grid.indexSet.index(T)` als Index f\u00fcr das Element `T`.\n", "\n", "Der Fehler $\\eta$ setzt sich aus den lokalen Fehlern zusammen:\n", "\\begin{equation}\n", "\\eta = \\left( \\sum_T \\eta_T^2 \\right)^{1/2}\n", "\\end{equation}\n", "\n", "Der lokale Fehler $\\eta_T$ wird aus dem Residuum und Spr\u00fcngen der Ableitung von $u$ auf den Kanten gesch\u00e4tzt:\n", "\\begin{equation}\n", "\\eta_T = \\left( h_T^2 \\|R_T\\|_{0,T}^2 + \\frac12 \\sum_e h_e \\|R_e\\|_{0,e}^2 \\right)^{1/2}\n", "\\end{equation}\n", "mit dem Residuumsterm\n", "\\begin{equation}\n", "R_T = \\Delta u_h + f.\n", "\\end{equation}\n", "Da wir ausschlie\u00dflich Lagrange-Elemente erster Ordnung betrachten f\u00e4llt der erste Term weg.\n", "\n", "Weiterhin beinhaltet der Sch\u00e4tzer die Sprungterme der Ableitung in Normalenrichtung\n", "\\begin{equation}\n", "R_e = \\left[ \\frac{\\partial u_h}{\\partial n} \\right] = \\left[ \\nabla u_h \\cdot n \\right].\n", "\\end{equation}\n", "\n", "Hierbei ist $h_T = \\mathrm{diam}(T) \\approx \\sqrt{2 \\mathrm{vol}(T)}$ der Elementdurchmesser von $T$ und $h_e$ die Kantenl\u00e4nge von $e$.\n", "In DUNE erhalten Sie mit `element.geometry.volume` den Term $\\mathrm{vol}(T)$ und mit `intersection.geometry.volume` die Kantenl\u00e4nge $h_e$.\n", "\n", "Verwenden Sie f\u00fcr alle Integrale eine Quadraturformel zweiter Ordnung.\n", "\n", "*Hinweis:*\n", "Aus der Basis erhalten Sie mit `grid = basis.gridView` wieder das Gitter." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def computeErrors(basis, x, f):\n", " pass" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 2 (Aufruf des L\u00f6sers und Fehlersch\u00e4tzers)\n", "\n", "Als weitere Hilfestellungen finden Sie hier zwei Methoden `addErrorsToVtkWriter` und `addFunctionToVtkWriter`.\n", "\n", "Mit der Funktion `addErrorsToVtkWriter(basis, etas, vtkWriter)` k\u00f6nnen Sie den gesch\u00e4tzten Fehler $\\eta_T$ zu einem `VTKWriter` hinzuf\u00fcgen.\n", "Hierbei ist `basis` die globale `dune-functions`-Basis, die Sie verwendet haben.\n", "\n", "Mit `addFunctionToVtkWriter(basis, x, vktWriter)` k\u00f6nnen Sie die durch die Basis `basis` und Koeffizientenvektor `x` gegebenen Funktion zu einem `VTKWriter` hinzuf\u00fcgen." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def addErrorsToVtkWriter(basis, etas, vtkWriter):\n", " etaFunction = grid.localGridFunction(lambda element, _: [etas[ basis.gridView.indexSet.index(element) ]])\n", " etaFunction.addToVTKWriter(\"eta\", vtkWriter, dune.common.DataType.CellData)\n", "\n", "def addFunctionToVtkWriter(basis, x, vtkWriter):\n", " u = basis.asFunction(x)\n", " u.addToVTKWriter(\"u\", vtkWriter, dune.common.DataType.PointData)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implementieren Sie eine Funktion `solveAndEstimateError(grid, filename)`, die das Poisson-Problem auf dem gegebenen Gitter l\u00f6st und den Fehler mit dem Fehlersch\u00e4tzer aus Aufgabe 1 sch\u00e4tzt.\n", "Die Funktion soll ein Tupel $(\\eta, (\\eta_T)_T)$ wie in Aufgabe 1 zur\u00fcckgeben.\n", "\n", "Geben Sie mit einem `vtkWriter` auch eine VTK-Datei mit der L\u00f6sung und dem gesch\u00e4tzten lokalen Fehler als Skalarfeld auf den Elementen aus." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def solveAndEstimateError(grid, filename=\"Adaptive-Refinement\"):\n", " pass" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 3 (Implementierung einer Verfeinerungsstrategie)\n", "\n", "Die folgende Funktion `refine(grid, etas, etaMax)` verfeinert alle Elemente $T$ des Gitters `grid`, deren gesch\u00e4tzter Fehler $\\eta_T$ oberhalb der Schranke `etaMax` liegt:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def refine(grid, etas, etaMax):\n", " from dune.common import Marker\n", " hierarchicalGrid = grid.hierarchicalGrid\n", " mark = lambda element: Marker.keep if etas[ grid.indexSet.index(element) ] < etaMax else Marker.refine\n", " \n", " hierarchicalGrid.mark(mark)\n", " hierarchicalGrid.adapt()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implementieren Sie eine Funktion `adaptGrid(grid, etas, refineFraction)`, die das Gitter anhand der im Parameter `etas` gegebenen lokalen Fehlerindikaturen $\\eta_T$ adaptiv verfeinert.\n", "Hierbei soll ein Anteil von Elementen mit dem gr\u00f6\u00dften Fehler $\\eta_T$ verfeinert werden.\n", "Dieser Anteil wird mittels des skalaren Parameters `refineFraction` bestimmt.\n", "\n", "*Hinweis:*\n", "Einen Vektor k\u00f6nnen Sie mit `x.copy()` kopieren und mit `x.sort()` sortieren." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def adaptGrid(grid, etas, refineFraction=0.3):\n", " pass" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 4 (Adaptiver L\u00f6ser)\n", "\n", "Implementieren Sie eine Funktion `adaptiveSolve(grid)`, die das Poisson-Problem auf dem Gebiet mit der einspringenden Ecke adaptiv l\u00f6st.\n", "Hierbei soll die Iteration abgebrochen werden, wenn der gesch\u00e4tzte Fehler den Wert $10^{-2}$ unterschreitet oder wenn 10 Verfeinerungsschritte erreicht wurden.\n", "\n", "Geben Sie bei jeder Iteration die Anzahl der Elemente und den gesch\u00e4tzten Fehler $\\eta$ aus.\n", "\n", "Verwenden Sie die folgende Funktion `filenameForLevel(level)` als Dateiname f\u00fcr die Ausgabe beim Aufruf von `solveAndEstimateError`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def filenameForLevel(level):\n", " return \"Adaptive-Refinement-lvl={:04d}\".format(level)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "def adaptiveSolve(grid):\n", " pass" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rufen Sie den L\u00f6ser auf.\n", "\n", "*Hinweis:* Wenn das Gitter bereits in einem vorherigen Durchlauf ver\u00e4ndert wurde, sollten Sie den Python-Prozess neu starten (Kernel -> Restart), da sonst auf einem bereits verfeinerten Gitter gestartet wird." ] }, { "cell_type": "code", "collapsed": false, "input": [ "adaptiveSolve(grid)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 5\n", "\n", "Betrachten Sie die Ergebnisse mit ParaView.\n", "In der Werkzeugleiste k\u00f6nnen Sie zwischen der Funktion `u` und gesch\u00e4tztem Fehler `eta` ausw\u00e4hlen.\n", "Das Gitter k\u00f6nnen Sie sich anzeigen lassen, indem Sie \"Surface With Edges\" an Stelle von \"Surface\" ausw\u00e4hlen.\n", "\n", "Was f\u00e4llt Ihnen bei der Verfeinerung auf?" ] } ], "metadata": {} } ] }