advection_diffusion_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nutils import *\n",
    "import numpy, unittest\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize\n",
    "ns = function.Namespace()\n",
    "\n",
    "# setup\n",
    "nelems = 20\n",
    "degree = 1\n",
    "tol  = 1e-5\n",
    "h = 1.0 / nelems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Peclet number is:  100\n",
      "The stabilization parameter is:  0.6135673098126083\n"
     ]
    }
   ],
   "source": [
    "# prepare unit domain\n",
    "domain, geom = mesh.rectilinear([numpy.linspace(0, 1.0, nelems+1)])\n",
    "\n",
    "# geometry and basis\n",
    "ns = function.Namespace()\n",
    "ns.x = geom\n",
    "ns.basis = domain.basis('discont', degree=degree)\n",
    "ns.h = h\n",
    "\n",
    "# solution vector\n",
    "ns.u = 'basis_n ?lhs_n'\n",
    "\n",
    "# material parameters\n",
    "ns.Pe = 100\n",
    "\n",
    "# supg stabilization parameter\n",
    "ns.pe = '0.5 Pe h'\n",
    "ns.xi = '(cosh(pe) / sinh(pe)) - (1 / pe)'\n",
    "ns.tau = '(h / (2 Pe)) xi'\n",
    "\n",
    "#  dg penaly parameters\n",
    "ns.C1 = 1       # penalty enforcing C^0 continuity\n",
    "ns.C2 = 1.0       # penalty enforcing C^0 continuity\n",
    "\n",
    "# analytical solution and forcing vector\n",
    "ns.pi  = numpy.pi\n",
    "\n",
    "ns.g   = '(1 / (1 - exp(Pe))) (exp(Pe x_0) - 1) + x_0'\n",
    "ns.dg  = '(Pe / (1 - exp(Pe))) exp(Pe x_0) + 1'\n",
    "ns.ddg  = '((Pe)^2 / (1 - exp(Pe))) exp(Pe x_0)'\n",
    "\n",
    "# ns.g   = '1 + ((exp(Pe x_0) - 1) / ( 1 - exp(Pe)))'\n",
    "# ns.dg  = '-(  Pe    exp(Pe x_0 )) / (exp(Pe) - 1)'\n",
    "# ns.ddg = '-( (Pe)^2 exp(Pe x_0 )) / (exp(Pe) - 1)'\n",
    "\n",
    "# ns.z = 1.5\n",
    "# ns.g  = 'sin( z pi x_0 )'\n",
    "# ns.dg = '(z pi) cos( z pi x_0 )'\n",
    "# ns.ddg = '-(z pi)^2 sin( z pi x_0 )'\n",
    "\n",
    "ns.f   = '- ddg + Pe dg'\n",
    "\n",
    "# print the Peclet number\n",
    "print('The Peclet number is: ',ns.Pe.eval()[0])\n",
    "print('The stabilization parameter is: ',ns.xi.eval()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare Residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = domain.integral('basis_n,0 u_,0 - basis_n,0 Pe u - basis_n f' @ ns, geometry=ns.x, degree=degree*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# res += domain.integral('tau Pe basis_n,0 ( - (u_,0)_,0 + Pe u_,0 - f )' @ ns, geometry=ns.x, degree=degree*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diffusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct residual and inertia vector\n",
    "res -= domain.interfaces.integral('-[basis_n] n_0 {u_,0} - {basis_n,0} n_0 [u]' @ ns, geometry=ns.x, degree=degree*2)\n",
    "\n",
    "# C0 continuity residual\n",
    "res -= domain.interfaces.integral('(1 + 0.5 Pe) C1 [basis_n] n_0 [u]' @ ns, geometry=ns.x, degree=degree*2)\n",
    "\n",
    "# boundary\n",
    "res += domain.boundary.integral('-basis_n u_,0 n_0' @ ns, geometry=ns.x, degree=degree*2)\n",
    "res += domain.boundary.integral(' (1 + 0.5 Pe) C1 basis_n,0 (g - u) n_0' @ ns, geometry=ns.x, degree=degree*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interface C^0 continuity \n",
    "res += domain.interfaces.integral('-[basis_n] n_0 Pe ({u} n_0 - .5 xi [u] n_0)' @ ns, geometry=ns.x, degree=degree*2)\n",
    "\n",
    "# weak boundary imposition\n",
    "res += domain.boundary.integral('basis_n Pe (0.5 (g + u) n_0 - .5 xi (g - u) n_0)' @ ns, geometry=ns.x, degree=degree*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solve > solver returned with residual 5e-14\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e2d6fd0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGgZJREFUeJzt3WtsXOd95/Hff4bD+00SZd0oiZYt3+JEtsPYMlq4WadeuEkTp4XhOotc2gZ1kqZJFy22aZE3QXbRLjZIty0SbCIkTbeNN06bphunTbbbrOPtthAVS5Zs17fUZnwRSVmUxBmR4gw5l2dfnCE1JOcmzuEMzznfDzAwZ+bwnOdwqJ8f/p/nOceccwIAhEes1Q0AAPiLYAeAkCHYASBkCHYACBmCHQBChmAHgJAh2AEgZAh2AAgZgh0AQqatFQcdGhpyIyMjrTg0AATWiRMnzjnnttfariXBPjIyouPHj7fi0AAQWGb2aj3bUYoBgJAh2AEgZAh2AAgZgh0AQoZgB4CQIdgBIGQIdgAImcgE+6WFnP76xOlWNwMANlxkgv3Rpyb123/1lE7PzLe6KQCwoSIT7EuBni9w824A4RaZYJ9KZlrdBABoisgE+0Qy3eomAEBTRCbYJ1MEO4BoiESwFwpOZ1KUYgBEQySC/dzcgrJ5Bk0BREMkgp36OoAoiUSwT1GGARAhkQj2SXrsACIkEsFOKQZAlEQi2FmcBCBKIhHszGEHECXRCPZkWh1tkThVAAh/sGeyeZ2bW9Tuwa5WNwUAmiL0wb604nT3YGeLWwIAzRH6YF+a6rh7gB47gGgIf7AXe+y7KMUAiIjwB3uxx75rgFIMgGiIRLAP9XYwKwZAZIQ+7SZTGQZOAUSKL8FuZoNm9i0ze8HMnjezO/3Yrx8mk2kGTgFEil899j+W9L+cczdIOiTpeZ/22xDnnBfsDJwCiJC2RndgZv2S7pL0y5LknFuUtNjofv2QSmc1v5inFAMgUvzosR+QNC3pa2Z20sy+YmY9Puy3YZPJpcVJ9NgBRIcfwd4m6TZJ/805d6ukS5J+d/VGZvaQmR03s+PT09M+HLa25cVJBDuACPEj2E9LOu2cO1Z8/i15Qb+Cc+6Ic27UOTe6fft2Hw5b29JVHXczhx1AhDQc7M65M5JeN7Priy+9Q9Jzje7XD5PJjBJx01BvR6ubAgBN0/DgadEnJD1sZu2SxiX9ik/7bchkMq1dA12KxazVTQGApvEl2J1zpySN+rEvP3nBThkGQLSEeuXpVCqjPQycAoiY0AZ7Ll/QmYsZZsQAiJzQBvvZ2QXlC067WJwEIGJCG+xTKeawA4im0Ab7RHHVKTV2AFET2mDnBhsAoiq0wT6VTKuvs019nYlWNwUAmiq0wT6RZKojgGgKbbBzHXYAURXaYJ9KseoUQDSFMtjnF3Oamc/SYwcQSaEM9kmmOgKIsFAG+9LiJEoxAKIolMHOnZMARFkog30imZGZtJMeO4AICmWwTybTuqqvQ4l4KE8PAKoKZfJNpZjDDiC6Qhnsk0muww4gukIX7M45b9Up9XUAERW6YL9waVELuQI9dgCRFbpgX1qcRLADiKrQBfvE0hz2AYIdQDSFLtgv3xKPGjuAaApdsE8m0+poi2lrT3urmwIALRHCYPemOppZq5sCAC0RvmBPpSnDAIi08AV7Ms3AKYBIC1WwL+YKOju7oF1MdQQQYaEK9jcuZuSctIdSDIAI8y3YzSxuZifN7G/92ueV4jrsAOBvj/03JT3v4/6u2OTynZMIdgDR5Uuwm9mwpHdJ+oof+1uvy5cToBQDILr86rH/kaTfkVTwaX/rMplMa0t3Qt3tba1sBgC0VMPBbmY/L+msc+5Eje0eMrPjZnZ8enq60cOWNZnkBhsA4EeP/ackvcfMXpH0iKS7zezrqzdyzh1xzo0650a3b9/uw2HXmkplqK8DiLyGg90593vOuWHn3IikByU95px7f8MtW4eJZJqpjgAiLzTz2GczWc1mcpRiAESer6OMzrnHJT3u5z7rNZXyZsSw6hRA1IWmx750gw1KMQCiLjTBzqpTAPCEKtjjMdNVffTYAURbaIJ9KpnRzv5OxWPcYANAtIUm2CeS3GADAKQQBftkKs3iJABQSIK9UHA6k8owcAoACkmwn5tbUDbvmOoIAApJsC/NYacUAwAhCfalVaeUYgAgJME+ubzqlGAHgFAE+0Qyre72uPq7uMEGAIQi2KeS3owYMxYnAUAogn0yxZ2TAGBJOII9mdbuAaY6AoAUgmDPZPM6N7dIjx0AigIf7GeY6ggAKwQ+2Jevw04pBgAkhSHY6bEDwArBD/Zij30nPXYAkBSSYB/qbVdnIt7qpgDAphD8YOdyvQCwQvCDPZnWbq7qCADLAh3szjkv2OmxA8CyQAd7Kp3V/GKee50CQIlAB/tkkqmOALBawIO9uDiJYAeAZcEO9hSrTgFgtWAHezKjRNw01NvR6qYAwKYR8GBPa9dAl2IxbrABAEsaDnYz22tmPzSz583sWTP7TT8aVg8v2CnDAEApP3rsOUm/7Zy7UdJhSR83s5t82G9NU6kMN7AGgFUaDnbn3JRz7sni17OSnpe0p9H91pLLF3TmIpcTAIDVfK2xm9mIpFslHfNzv+WcnV1QvuC0i8VJALCCb8FuZr2S/lrSv3fOXSzz/kNmdtzMjk9PTzd8vKkUc9gBoBxfgt3MEvJC/WHn3LfLbeOcO+KcG3XOjW7fvr3hY04UV51SYweAlfyYFWOSvirpeefcHzbepPosrTplVgwArORHj/2nJH1A0t1mdqr4eKcP+61qKplWX2eb+joTG30oAAiUtkZ34Jz7J0lNXyE0kWSqIwCUE9iVpyxOAoDyAhvsUylusAEA5QQy2OcXc5qZzxLsAFBGIIP98g02KMUAwGqBDPblxUncxBoA1ghksHPnJACoLJDBPpHMyEzayawYAFgjkME+lUzrqr4OJeKBbD4AbKhAJuMkUx0BoKJgBnsyw8ApAFQQuGB3zmkymWaqIwBUELhgv3BpUQu5AqUYAKggcMF+eXESwQ4A5QQu2CeSLE4CgGoavmxvs12+JR41dgDrtLgopVJSMuk9ZmYuf13pcfGi9325nPeQpN5eaWhIuvZaaXRUeu97pd27W3tuCmCwTybT6miLaWtPe6ubAqBVstmVwVzpUSmw5+er77+tTRocXPnYtUtqb5cSCe9956RLl6QzZ6Tvf1/62tekj39c+oVfkL7whZYGfACDPaPdg13y7sgHIJDy+fLBXE/POZmU5uaq7z8W88J4y5aVwbw6rJcepdsNDkrd3dKVZszzz0vf+Ib0uc9Jb3qT9PWvS+961/p/Rg0IXrCnmOoItFyh4JUm6g3i1dvNzlbfv9na8D14sHIQr3709l55MDfqxhulz35W+sAHpAcekB58UPrRj7zXmyx4wZ5M666D21vdDCDYCgUvXK+0hFFab3au+jEGBlYG8TXXVA/j0sDu7fV63UF08KD03e9Kt90m/eIveuHe19fUJgQq2BdzBZ2dXdAupjoi6pzzyhFXWsJY2i6Vqh3M/f0rQ3f/funQofp6zX19UjzenJ/FZjQ8LD3yiHTPPdLv/770B3/Q1MMHKtjfuJiRc9IeSjEIuqWBt/UO/iWTXq+7mt7elWE7PCzdfHN9veb+/mgHsx/uvlt697ulr35V+sxnpI6Oph06UMHOddixaTgnpdPrH/xLJi9Pmaukp2ftrIwbb6w98Dc46JVB2gL1zzucPvYx6Tvfkb79bel972vaYQP1yU8W57DvYnES/JDJrH/wL5n0ptxV09W1Mmy3b5euu652j3kpmNuZ0ht499wjHTggfelLBHsl3OsUKywsrH/wL5n0vr+ajo6VveFt2yoPAK7uNQ8MNPVPb2xSsZj0kY9In/qU9Oyz3jTIJghYsKe1pTuh7vZANRuVlK7+W0/POZOpvv9EYm3g7t9f33S5wUGpkw4EfPChD3nB/t3vEuzlfOSua/SeQ61frouielb/VQvsK1n9txTEw8P1TZdbCmYWsqHVduzwLjlw7FjTDhmoYN+3rVv7tnW3uhnhUW7135X0nC9dqr7/eLz8sux6F5msZ/UfsBndcYf02GPeoHsTfqcDFexYpVCo73oZlQK71uq/WGzlIpPBwfoG/5YCu6eHYAYk6fBh6eGHpdOnpb17N/xwBHsrVVv9V+/V5qotMjFbG8z1Dv4FffUfsJnccYf332PHghPsZnavpD+WFJf0Fefcf/Zjv5vJsfHzumXfoDraShZtlK7+W8/g33pX/91yS31T5vr7CWZgMzh0yJslNTYm3X//hh+u4WA3s7ikL0q6R9JpSU+Y2aPOueca3feGq2f138yMZibO6sLJcc1tiakjO3/lq/9Ke8N790pvfnN9vWZW/wHh0N4u3Xpr0wZQ/eix3y7pJefcuCSZ2SOS7pO08cFebvVfhZ7z7T+Z0l+8NKGdfx+T5kquSpfPVz9GT48S3b064Dqk7h3SnpLVf7UG/1j9B2DJ4cPSl7/szSZLJDb0UH6kzh5Jr5c8Py3pDh/2u9bnPy9985vrWv3X39Wnrlxc+aF90k031DddbmBASiT0yT97Qo+9cFZ/9Eu36L237tmQUwMQPPmC02wmq4vpnC5mskqls7qYzupiyWve85yuzwzpo+m0XvzBUV3/c3dtaLv8CPZy0x7WFI7N7CFJD0nSvn371nekRGLt6r9qveaS1X8/ePK0fusvn9L//Q9v1/5tPXUfMpcv6Ec/ubC+9gLY1AoFp9mFXMUwrhjUxffnFqpf7ydmUn9XQv2dCaV2HNTbrz0ky6Q3/Lz8CPbTkkqHeYclTa7eyDl3RNIRSRodHa0xYljBJz/pPZro2cmLNT88AK3hnNOlxfxy+KbmvcAtF8apMq/NLuRqzl/o62xTf2eiGNBt2re1ezms+7suvzdQfL+/6/K2Pe1tisVK+r6feWBjfyBFfgT7E5IOmtnVkiYkPSjp3/mw303h6Pj5VjcBCC3nnNLZ/KqecklpY35VGBe/vhzSWRVqBHNPe9wL3WIY7x7s1A2dfSsCeHVQL23f29GmeCx4azEaDnbnXM7MfkPS38ub7vinzrlnG27ZJjE2fl6diZgy2RqzX4CIymTzZQN3Tc+5wnvZfPVk7krEV/SMh3rbdc32niq95suv9XW2qS0evSm/vkzZcM59T9L3/NjXZpLLF/TETy5odP9W/dNL51rdHGBDLOYKZQPXqy/nar63mKve6WmPx4qh6/WMB7vbtW9bz3JPeXUYL/WiB7oS6utMqL0tesHcKObiVfHMREqXFvO685ptBDs2rVy+oNlMrkIYl5uxsbInnc5Wn/LbFrOSUoYXxrsHuoqljJVhXC6oOxOsxWg2gr2KsXFvNswdV29tcUsQZvmC01xpMNcoXawueVxarB7MpTMzBophfFVf78qacne5XrP3vCsRl3HNn0Ah2KsYGz+vg1f1alsvN0xAZYWC06XFnDc9br58GFcra8xmqs+6MpP6OlYO8O3f1l29jLEc1An1tBPMUUOwV5DNF/TEKxd0/1uHW90UbDDnnOYX8yvKFsvT46r1lovPZzO1Z2b0drStmAq3Z7BLN+7qW1NTLjcg2NexasocUAPBXsEzEynNL+Z1+MC2VjcFNTjntJArVJyrXHWhSTG0czWSubs9viJwd/R36uBVvRXDuLTs0dsRzZkZaB2CvYKx4vz1O67eqos1/lRG4xZzBc0u15jLlzDW9KBLgnoxX31mRkdbbMUA39aedo1s66kYxqU9aWZmIGgI9gqOvnxe1+/o07beDoK9DmsGACsM9K3uPS8F93pmZuzZ0rXcU75ca1654MSbMtfGzAxECsFeRjZf0PFXZvTAaHTq60t15tTqQK4Q0Kt71nM1lmabaWW5ojOhA0O9axeXrFr5x8wM4MoR7GU8fTqldDY49fWlCxml5rNKpheVSmeVnM8qmc5qYiatPYOddZU38jXqzCsGADvXDgAOlOkt9xcXpfSuvmYGgA1DsJexXF/3KdiPvnxeL0/PqaMtps5EfO1/EzF1tq38b0dbXPGY6eXpOT3+4rRS88XALgntpddSdVwvQ/LqzKXljK097bp6qKfsNTJWlziiujQbCCKCvYyx8fO6YWeftva0+7K/T3zjpM7NLVzx9yXitnwdjaW5zIPd7Rrs9gJ439ZuDXYllp8PFJdrD3YnNFgM6LmF3HIwU2cGooFgX2Ux59XXf+lt/t1wNlco6P63DuuTdx/UQi6vhVxBmWz1/y5kC8rk8lrIFvSW4QG9+9DudV1lbodvZwEgKAj2VZ4+nVxTX48XB+3+0989p785OaEbdvbp+uLjmu29dfWEe9rj2rete8PaDQBLCPZVSuevL9m7tUufve9NevLVGb34xpyOvnx+ed50PGYa2datG3b26/qdfbpuR59u2NmnfVu7GSwE0BIE+ypHi/X1LSX1dTPTB+8c0QfvHJHkTYd85dwlvXBmVj9+Y1YvnJnVMxMp/d0zU8vf05WI67odvbpuR5/SNS7SBAB+IthLLOTyOvHqjN53e/V7sibiMR3c0aeDO/pWvH5pIad/PTunF89c1AtnZvXimVk99sJZLeQKGuJCYgCahGAv8fTplDLZwrrnr/d0tOmWvYO6Ze/gitdT81n1d/GjBtAcpE2Joy+fl5n/118f6E74uj8AqIYVJyXGxs/rxp39Guz2Z/46ALQCwV60VF8PymUEAKASgr3o1GtJLeQKuvMagh1AsBHsRWPjF2Qm3T7C/U0BBBvBXjQ2fl437epnoBNA4BHskjLZvE68NqM7qa8DCAGCXdKp15NazK1//joAbCYEu7wyTMykt/k8fx0AWoFgl7cw6U27BzTQRX0dQPBFPtgz2bxOvp7U4QP01gGEQ+SD/eRr1NcBhEvkg/0o9XUAIdNQsJvZ58zsBTN72sz+xswGa3/X5jI2fl437xlQfyf1dQDh0GiP/R8k3eyce4ukH0v6vcab1DyZbF6nXktShgEQKg0Fu3PufzvncsWnY5KGG29S8zz56owW8wUWJgEIFT9r7L8q6fuV3jSzh8zsuJkdn56e9vGw6zc2fl7xmGl0ZEurmwIAvql5ow0z+4GknWXe+rRz7jvFbT4tKSfp4Ur7cc4dkXREkkZHR926Wuuzo8X6eh/1dQAhUjPYnXM/W+19M/uQpJ+X9A7n3KYI7HqkF/M69XpSv/rTV7e6KQDgq4ZujWdm90r6lKSfcc7N+9Ok5njytRll846BUwCh02iN/QuS+iT9g5mdMrMv+dCmpjj6sldffxvXXwcQMg312J1z1/rVkGYbGz+vN+8ZUG8H9/MGEC6RXHk6v5jTU6eZvw4gnCIZ7Cde9err3N8UQBhFMtiX56/vZ/46gPCJaLBf0FuGB9RDfR1ACEUu2C8t5PXU60kuIwAgtCIX7Cdem1GuwPx1AOEVuWAfGz+vtpjprdTXAYRU5IL92Ph5Hdo7SH0dQGhFLtjPzS1yf1MAoRa5YJdEfR1AqEUu2BNx6usAwi1ywX5oeFDd7dTXAYRX5IKdywgACLvIBTv1dQBhF5lg/+lrh/Trb79Gt1/NjBgA4RaZYvNV/Z36nXtvaHUzAGDDRabHDgBRQbADQMgQ7AAQMgQ7AIQMwQ4AIUOwA0DIEOwAEDIEOwCEjDnnmn9Qs2lJr67z24cknfOxOUHAOUcD5xwNjZzzfufc9lobtSTYG2Fmx51zo61uRzNxztHAOUdDM86ZUgwAhAzBDgAhE8RgP9LqBrQA5xwNnHM0bPg5B67GDgCoLog9dgBAFZs22M3sXjN70cxeMrPfLfN+h5l9s/j+MTMbaX4r/VXHOf+WmT1nZk+b2f8xs/2taKefap1zyXb3m5kzs0DPoKjnfM3sgeLn/KyZ/Y9mt9Fvdfxe7zOzH5rZyeLv9jtb0U4/mdmfmtlZM/uXCu+bmf1J8WfytJnd5msDnHOb7iEpLullSQcktUt6StJNq7b5dUlfKn79oKRvtrrdTTjnfyOpu/j1x6JwzsXt+iT9o6QxSaOtbvcGf8YHJZ2UtKX4/KpWt7sJ53xE0seKX98k6ZVWt9uH875L0m2S/qXC+++U9H1JJumwpGN+Hn+z9thvl/SSc27cObco6RFJ963a5j5J/7349bckvcPMrIlt9FvNc3bO/dA5N198OiZpuMlt9Fs9n7Mk/UdJ/0VSppmN2wD1nO+vSfqic25GkpxzZ5vcRr/Vc85OUn/x6wFJk01s34Zwzv2jpAtVNrlP0p87z5ikQTPb5dfxN2uw75H0esnz08XXym7jnMtJSkkK8p2q6znnUh+W93/8IKt5zmZ2q6S9zrm/bWbDNkg9n/F1kq4zs382szEzu7dprdsY9ZzzZyS938xOS/qepE80p2ktdaX/3q/IZr3nabme9+rpO/VsEyR1n4+ZvV/SqKSf2dAWbbyq52xmMUn/VdIvN6tBG6yez7hNXjnm7fL+Ivt/Znazcy65wW3bKPWc8/sk/Zlz7vNmdqekvyiec2Hjm9cyG5pfm7XHflrS3pLnw1r759nyNmbWJu9PuGp/+mx29ZyzzOxnJX1a0nuccwtNattGqXXOfZJulvS4mb0irxb5aIAHUOv9vf6Ocy7rnPuJpBflBX1Q1XPOH5b0l5LknDsqqVPe9VTCrK5/7+u1WYP9CUkHzexqM2uXNzj66KptHpX0oeLX90t6zBVHJQKq5jkXyxJflhfqQa+9SjXO2TmXcs4NOedGnHMj8sYV3uOcO96a5jasnt/r/ylvkFxmNiSvNDPe1Fb6q55zfk3SOyTJzG6UF+zTTW1l8z0q6YPF2TGHJaWcc1O+7b3Vo8dVRpXfKenH8kbUP1187bPy/mFL3of/V5JekvQjSQda3eYmnPMPJL0h6VTx8Wir27zR57xq28cV4FkxdX7GJukPJT0n6RlJD7a6zU0455sk/bO8GTOnJP3bVrfZh3P+hqQpSVl5vfMPS/qopI+WfM5fLP5MnvH795qVpwAQMpu1FAMAWCeCHQBChmAHgJAh2AEgZAh2AAgZgh0AQoZgB4CQIdgBIGT+P9qsOaef4DxrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b48f898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find lhs such that res == 0 and substitute this lhs in the namespace\n",
    "lhs = solver.solve_linear('lhs', res)\n",
    "ns = ns(lhs=lhs)\n",
    "\n",
    "# plot solution\n",
    "bezier = domain.sample('bezier', 9)\n",
    "x, u = bezier.eval([ns.x, ns.u])\n",
    "x, g = bezier.eval([ns.x, ns.g])\n",
    "plt.plot(x, u)\n",
    "plt.plot(x, g,'r')\n",
    "# plt.ylim(-1.0,1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}