{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SolKx analytic solution with error analysis\n",
    "\n",
    "The SolKx analytic solution is run over a range of resolutions with the error (in velocity and pressure) measured at each configuration. The error convergence is analysed against the expected element type order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "import glucifer\n",
    "from underworld import function as fn\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# SolKx object contains the analytic solution and parameter (ie, viscosity, density) profiles\n",
    "sol = fn.analytic.SolKx(B=10.0, nx=2.0, nz=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "elTypeToUse = \"Q1/dQ0\"\n",
    "\n",
    "# expected order of convergence results for Q1/dQ0\n",
    "# velocity for Q1 -> 2nd order convergence\n",
    "# pressure is dQ0 -> 1st order convergence \n",
    "expected_order = np.array([1.9,1.9,0.9])\n",
    "\n",
    "# Resolutions to test against\n",
    "res = [(32,32), (64,64), (128,128), (256,256)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def analyticMeasure( numeric, analytic, mesh ):\n",
    "    '''\n",
    "    Measure the integrated error between 2 MeshVariables (or uw.Functions) over the mesh\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    rel, err: ndarray, ndarray\n",
    "        Array representing the relative and absolute errors between the MeshVariables\n",
    "        The lengths are the number of components of the field\n",
    "    '''\n",
    "    \n",
    "    # functional description of L2 norm sums; error and magnitude\n",
    "    fn_err = fn.math.pow(numeric - analytic, 2.)\n",
    "    fn_norm = fn.math.pow(analytic, 2.)\n",
    "    \n",
    "    # perform L2 norm evaluation        \n",
    "    err = fn.math.sqrt(uw.utils.Integral( fn_err, mesh ).evaluate()).evaluate()[0]\n",
    "    norm = fn.math.sqrt(uw.utils.Integral( fn_norm, mesh ).evaluate()).evaluate()[0]\n",
    "    \n",
    "    # calculate relative error measures\n",
    "    rel = err / norm\n",
    "    \n",
    "    return rel, err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def analyticSol_SingleResolutionTest( sol, elType, resolution ):\n",
    "    \"\"\"\n",
    "    Test a FEM model against an analytic solution\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    sol : fn.analyic._SolBase\n",
    "        The python wrapped analyic solution.\n",
    "        \n",
    "    elType: string\n",
    "        The element type desired, as these equation solve stokes this must\n",
    "        describe a mesh/subMesh pair, eg 'Q1/dQ0'.\n",
    "    \n",
    "    resolution: integer tuple\n",
    "        The resolution of the mesh.\n",
    "    \"\"\"\n",
    "\n",
    "    # 'global' python vars are visible outside this function's scope\n",
    "    global mesh\n",
    "    global velocityField\n",
    "    global pressureField\n",
    "    global solver\n",
    "    \n",
    "    mesh = uw.mesh.FeMesh_Cartesian(elType, resolution, (0.,0.), (1.,1.))\n",
    "\n",
    "    velocityField = mesh.add_variable(2)\n",
    "    velocityField.data[:] = (0.,0.)\n",
    "    pressureField = mesh.subMesh.add_variable(1)\n",
    "    pressureField.data[:] = 0.\n",
    "\n",
    "    # currently not checking strain rate\n",
    "    # strainRate = fn.tensor.symmetric( velocityField.fn_gradient )\n",
    "    \n",
    "    # freeslip\n",
    "    IWalls = mesh.specialSets[\"MinI_VertexSet\"] + mesh.specialSets[\"MaxI_VertexSet\"]\n",
    "    JWalls = mesh.specialSets[\"MinJ_VertexSet\"] + mesh.specialSets[\"MaxJ_VertexSet\"]\n",
    "    freeslip = uw.conditions.DirichletCondition(velocityField, (IWalls, JWalls))\n",
    "\n",
    "    stokesSystem = uw.systems.Stokes(velocityField,pressureField,\n",
    "                                     sol.fn_viscosity,sol.fn_bodyforce,\n",
    "                                     conditions=[freeslip,])\n",
    "    #Run the BSSCR Solver\n",
    "    # # can optionally set penalty this way\n",
    "    solver=uw.systems.Solver(stokesSystem)\n",
    "    solver.options.scr.ksp_rtol=1e-9  # tolerance for the scr solve, default is 1e-5.\n",
    "    solver.options.A11.ksp_rtol=1e-5  # tolerance for the inv K approximation\n",
    "    solver.solve()\n",
    "    \n",
    "    v_err, abs_v_err = analyticMeasure( velocityField, sol.fn_velocity, mesh )\n",
    "    p_err, abs_p_err = analyticMeasure( pressureField, sol.fn_pressure, mesh )\n",
    "    \n",
    "    return v_err, p_err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# store the error measures\n",
    "vx = []\n",
    "vy = []\n",
    "p  = []\n",
    "\n",
    "ssize = len(res) # the number of samples\n",
    "\n",
    "for i in xrange(ssize):\n",
    "    test_err = analyticSol_SingleResolutionTest( sol, elTypeToUse, res[i])\n",
    "    \n",
    "    vx.append(test_err[0][0]) # get vx error\n",
    "    vy.append(test_err[0][1]) # get vy\n",
    "    p.append(test_err[1][0])  # get pressure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res_x = np.asarray( res )[:,0] # assume square resolutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h = np.log(np.reciprocal(res_x.astype(float)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = np.vstack([h, np.ones(len(h))]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# take the log of all dependent variables\n",
    "yvx = np.log(vx)\n",
    "yvy = np.log(vy)\n",
    "yp  = np.log(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vxm, c = np.linalg.lstsq(A, yvx)[0]\n",
    "vym, c = np.linalg.lstsq(A, yvy)[0]\n",
    "pm, c = np.linalg.lstsq(A, yp)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test = np.array([vxm, vym, pm])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if uw.rank()==0:\n",
    "    np.all(np.less(expected_order, test))\n",
    "\n",
    "    if not np.all(np.less(expected_order, test)):\n",
    "        raise RuntimeError(\"The numerical order of convergence is outside the error threshold.\"+\n",
    "        \"The Relative error found were \", test,\" the threshold is \", expected_order )  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VNX5wPHvi0ECKAQXEERj2KxAiFgqbmhwAUVFiooF\nQUkQBFraqFQF/AFSUASLWpGwVGi1BTdAgoIkLFFA2QRCAJVFQEEFQUJMIECS8/vj3qRDmMlMkknu\nLO/nefJk5t4z577n3jvvnDlz5o4YY1BKKRX6qjkdgFJKqaqhCV8ppcKEJnyllAoTmvCVUipMaMJX\nSqkwoQlfKaXChCZ8PxGRFSKS6Ke6bhCRHSKSLSJd/VDfKBF5uwKP3yoiN1c0jkAlIp1FZJ7TcfiT\niMSKyOoAiKNQRJo4HUdlEZFfReQKD+seFZGVpTz2FhH5vrJic8eRhC8ie0TkVie2HSTGAP8wxtQx\nxqT4qU6fvnAhIrNEZMwZDzSmtTHmMz/FEYjGAi8W3fGUpESki4isFJGjIvKDiEwXkdpVGqmPjDGZ\nwFERudvpUBze/lnsDtBb/qjLGHO+MWZvaUVctuvuvKrS/aM9fEBEznE6hhKige2+FAzA2CuFiIgv\ny7zUcda+EpF2QB1jzHqXxZ6ehHWAvwENgauAxsDEssRQHhU4xrOBgf6MpRzKdIxCnPMvfsaYKv8D\n9gC3eljXH9gJHAY+BBq6rOsEfA0cBd4A0oFED/VEAv8GfgG2AX8Fvi8Rw9NABnAC68XvGWAXkA1s\nBbq5lH8UWAW8DmRhJeRbXdavwOqZr7If/wlwQSn7oGQ7L7GX7wLygeN2PdU97L+SsTcEPgAOAbuB\nIS7lRwFvudx/D/jR3o/pwFUuMZ0C8uxtL3A9XvY2jgNRLnW1BX4GzrHvJ9r75giwGLi8lH1wHbDa\njmMTcEuJ/TnW3p+5QBMPyxoCC+zt7QAeK9Hu94G37WN21rkC/B8wvcSyQqCJD+fx74EML+f5U/Zx\nOgrMAc718VwvBAbbbdrtsmyQveyYfb41sfdhFvAOEOFSRyP7eLk7h3oA60ssewL40L7dBet5kw18\nDzxZSjs9HnPXfQmcC7wM7LPPvylADXvdLfZ2/gocBA4A9wF3Ad/Y+2iYS70CPIv1fPnZbnuUvS7a\n3u4j9rYOAcPtdZ2Bk/bfr8AmN+3pC6S43N8JvOty/zugjZv2XQCk2MdmjX18PrPXfWqXzbH36YMu\nbX7Spc19/ZVn3R6ryqzcyxPhrISPlVR+BuKA6sA/gE/tdRfZO/I+rAT3Z/ugeUr447ESRB37xM8A\nvisRw0Z7XdFJdz/QwL79oH1wiu4/Cpy2t3sO1hMmy+UkW2GfGE2BGvb9FzzE5rGdLrF19LL/imO3\nT/4NwAg7tivsJ8IddvmSCb8vUMve9iTXkx6YBYzxdLyApUA/l3UTgCn27fuwklEL+xgNB1Z7aEMj\nrCdxZ/v+bfb9C132517gN3ZdER6WfYr1Ilzd3p+HgHiXdp8E7rXv13ATx3vAUyWW+ZrwXwVmezlO\na4AGQBRWUhzg4zlQCCyxH1fDZdl8oDbWO4w8IA0rwZ2PlaD7lIjhGNDaTWw17XVNXZatAx60b/8A\n3GDfrgtc7aGNpR5zzkyIr2C9sNW127AAGGevuwXr+VV0Dj9mH8v/YJ2rLbFevKLt8n8BPsd6wa8O\nJBcdC/6X8Kdhvci0sffVle6eD27aFAP8Yt9uaJ9z39n3mwBHXMoWuLTvHfsvEmgF7MdO+C77Isbl\nflGbR9ltvgurI1O3vLnV6zlbWRV7eaJ4Svj/BMa73K+N9YS9HOhDieSB9UrrKeHvBm53ud+PsxP+\no17i3MT/ksWjwP4S69cCD9u3V2D3Iuz7g4BFHup1185T2D0jT/vHU+zAtcDeEmWeBd70doJjJZRC\n4Hz7vreE3w9YVuIY3GjfXgQkuKyrZp/Al7nZ7tPAv0ss+wQ7Ydn7c3SJ9WcswxpSOQ3Ucln2AjDT\npd3pXo5xKnYSdlnmNeEDd2D1aJuWUmYP0NPl/kv878XR2zlQiMs7Hpdl17nc3wD81eX+y8CkEo/Z\nD9zkIb63gOfs282xXgCKXlz2Yr0DOd/Lfij1mHNmws/hzIR3PfCtffsW+3Fi3z/Pfmy7Eu3tat/e\njkunCCsxn7K3H42ViF3fMa0Fenh7PriU3wdcDTyE9cKxButFrS/2uyDX9tnbPQU0d1k3jrMTfhOX\n+0Vtruay7CBwbWmxVeQv0MbwG2HtaACMMblYQzKX2utKfqK930tdruvdfRp+xuNF5BER2WR/KHcU\n61X6IpciB0o8fp+9nSI/udw+jnXSeoqtZDuPYLXTV66xRwOXisgv9t9RYBhQv+SDRKSaiIwXkV0i\nkoWVlAxntrM0c4HrRKSBiNwCFBhjimaDRAOvFcVht8l4aFc00KNEzDcCl7iUcXfMXJc1wuqJHXdZ\ntq/E9rzNgjiK1Tv2mYhcB/wXuN8Ys9tL8YMut13PCV/OAXfn9yGX2ydK1H+Cs8+587HeibozB+hp\n3+6FlchO2vfvB+4G9tkz0K7zUIdPx1xELsbqqX/pUnYxcKFLsSPGznp2W9y1t6h90cB8l7q2Y734\nN3Ap72nf++JToCNwM9awZzoQj5WkP3VT/mKsXrrrMdvnplxJR4wxhRWIs0wiKqvicvoB60ACYM+A\nuBAr0f4IlJyi2NhLXY2xxvzBepdQUtHJhYhcDkzH6jV8YS/bxJkfOpVMXJdjvS0tK0/tLO0FrCTj\ncvt7rJ7SlT487mHgXqwe+3ciUhcr6RW103h8JGCMyRKRVOAPWMMK77is/g4Ya4yZ40Mc32P1sh4v\nbXNelv0AXCAite2ECdYxOeChvDtbsHpuPhGRtljDEn2NMem+Ps4NX84Bb7GXSkQaYQ13fOOhSBpw\nsYjEYR3PpOING/Ml0M3+wHgI1tCXu+eQr8f8MFYya2WM+bFMDXGv6N39FyVXiEi0m/KufNmvn2E9\nT67A6qkfw3ruXIc1hFjSz1ifvV2GNcQF7veXo5zs4Z8rIjVc/s7B6nEkiEgbEamB9fZ8jTHmO+Bj\noLWIdBWRc0TkT5z5al7S+8AwEYkSkUuBP3qJpzbWW67Ddi84AWhdokx9ERkiIhEi8iDWWPLHZW65\n53aWd07uOuBXEXlaRCLt/dPKnoFS0nlYw2RH7STzImc+AQ5ivUX1Fv8jWL3A2S7LpwHDRaQlgIjU\nFZEHPNTxH+BeEelk7+9Ie15yIw/lz2KM2Y81jvuifQ61wRpyKst3DhZh9dxKcj03a9gxtsbqlQ4x\nxiwqwzbc8fc54M4twHJjzGl3K40x+VjPk4lAPawXAESkuoj0EpE6xpgCrA83Czxsw6djbvfcZwCv\n2r19RORSEelUzrZNA16wO2qIyMVy5ndWSpsddBC4wsssr6Iefk1jzA/ASuBOrBflTSUL2730ecBo\nEalp749HSxT7Ce/PrUrlZML/GOsV/4T9f5QxZhnWrIl5WL20GKyeB8aYI1gfpE7E6i38BmtM7+RZ\nNVvG2HXswRqnfb9E2TNe5Y0xXwF/xxqr+wlrOGdViTrXYo11Hsaanne/Mabo7bLPvbHS2uljXSVj\nLwTuwRpz3IP1NngG1gfWJb2F1Ts6gDUT6fMS698EWtlvlYu+jFQynhSs/fCjseZ7F8XxIdaH5e/Y\nw0VbsJ4kZzfAStb3YX3I9zPW29+h/O+c9Na7L9ITa//9gDXc9H/GmBXutukhjk1Aloj8rsR2tnLm\n+dkXazbFRcCbYn3h5lcRycQzj8exnOdAyWXezpOHgaleyszB+sD8vRJDC32APfZxHIA15HN2QN6P\nuWuMRbPg1thlUyn93VVp7X0N6911qogcwzqPr/Xxse9jvSAcEZENHtq1E+uF7jP7/q9Ynwuuchl2\nKlnvEKwhtB+Bmfafq9HAW/Zzy1NHqELv6ryRM2Mv44NFJmC97TmJtTMSjDHZforN27YF6+1vL2OM\nuzG1kuUHAg8ZYzqWc3uPYs1OCdlvnIYrEbkDGGSM6e50LP4iIrHAVGPMjU7HogJHRXv4qVhjcldj\nTUkcVvGQPLPf/te13wKPsBev8VD2ErEuUSAiciXWfOiQ+vq88g9jTFooJXuwvmmryV6VVKGEb4xZ\n6vI2cA2lf4jqD9djvZM4hDWD4D6XWQUlnYs1zpeNNXd8PtZcXaWUCksVGtI5oyKRFOAdY8xsr4WV\nUkpVOa/TMkUkjTNnwwjWBwsjjDEL7TIjgNOa7JVSKnBVuIcvIn2xvpF3aynDK4hIpX76rJRSocoY\n45eL0FVoDF9E7sS62FHX0pJ9kcr6unAg/I0aNcrxGLR92rZwa9/I+HhGYQ05lPY3smNHx2Mt758/\nVXSWzutYX+RJE5GNIjLFDzEppZRP8mvV8ulLK/m1alVFOAGvQpdWMMY091cgSilVVjclJPDOkiVQ\n4OmLwLAiMpIOiX75MbqgF2gXTwta8fHxTodQqUK5faHcNgjt9nXu3p1jLVqQ62F9LjA3Lo5O3bpV\nZVgBy2/TMr1uSMRU1baUUuHj0KFDPN+1K/dnZNAxL694GuGKyEjmxsUxKiWF+vXPunBs0BARjJ8+\ntHU84V9xxRXs2+fLVUSDV3R0NHv37nU6DKVCVmFhIUvmz2fVrFlEHD9Ofq1adEhMpFO3blSrFtwD\nGSGV8O3GVEkMTgmHNiqlKoc/E35wv/QppZTymSZ8pZQKE5rwlVIqTAR8wn/llX8GRB1KKRXsAvpD\n2/3799O6dRe2bVvMpZeW5fe9/VtHRemHtkqp8gqbD22nTJnPsWP/ZMqU+Y7WoZRSoSCgevhjxrzB\nokW7iIysC8CBA7Br12iaNRtNUec8L+8YXbo0Y+RI979J7o86ikyYMIH169fz/vvvFy9LSkrCGMO8\nefNITk7mnnvuITc3l6uvvppRo0bRu3fvUtuolFJl4c8eflVe8c2447r86NGjpl27wQaOGjBu/n4x\n7doNMkePHnVbl7/qKLJv3z5Tu3Ztk5OTY4wxpqCgwDRs2NCsW7fOpKammoYNG5pDhw6Zxx57zPTo\n0cNjPZ7arpRS3tj5wz952F8Ved2QDwnfmNIStu+J2h91FOnQoYN5++23jTHGpKammmbNmhWv+/Of\n/2xiY2NN48aNzS+//OKxDk34Sqny8mfCD7gx/KioKNLSxtGixdgzlrdoMY60tBeIioqqkjqK9OzZ\nkzlz5gAwZ84cevXqVbyuf//+bN26lb59+1KvXj2f61RKKScEXMIHyMnJITs7moiIbTRvnkRExHay\ns6PJzfV0TbzKqQPgwQcfJD09nQMHDjB//vzihF9YWMiAAQN49NFHmTJlCt9++22Z6lVKqaoWkAl/\n8uS5FBTsIykpjczMl0hKSqWgYB9vvDG3SusAuOiii7jllltISEigSZMmXHnllQCMGzeOatWqMXPm\nTIYOHUqfPn30g1mlVEALyIRfp05N0tISmTgxiRo1ajBxYhJpaYnUqVOzSuso0qtXL5YtW8bDDz8M\nwMaNG3n11Vd5++23ERGeeeYZqlWrxvjx48tct1JKVZWAmpYZqsKhjUqpyuHPaZkV+olDpZR/FRQU\nsGTePFb/61/F13W/KSGBzt27B/113ZXztIdfBcKhjariin656YGMDOJdfrkpPTKSD0Lgl5tU+egP\noASZcGijqpjCwkKG3HADE9aupbab9bnA0+3b8/rnn2tPP8yEzbV0lAoXS+bN44GMDLfJHqA2cH9G\nBqkffliVYakQowlfqQCwatYs4vPySi3TMS+PlTNnVlFEKhRpwlcqAEQcP4639+xil1OqvDThKxUA\n8mvVwtunPMYup1R5acJXKgDclJBAemRkqWVWREbSITGxiiJSoUgTvlIBoHP37nwQF4enKz3lAnPj\n4ujUrVtVhqVCjCZ8pQJAtWrVGJWSwtPt27M8MrJ4eMcAyyMjebp9e0alpOiUTB/ob1h7VqF5+CIy\nBrgPKAQOAn2NMT95KKvz8JXyorCwkCXz57Nq1qzib9p2SEykU7dumux9EAi/Ye1vAfPFKxE5zxiT\nY98eArQ0xgzyUDboEr67nzj8y1/+QlpaGrVq1WLDhg3FyydNmsTKlSuZP//s384N5DYqFUqGDX+N\n8S9ez/Dh6xg37k9Oh+MXAZPwz6hI5FngMmOM2x+KDcaE/91339GyZUsOHjxI7dq1KSwspHHjxrz7\n7rv8/ve/Z/Xq1cWXS77mmmsYOXIk3dyMsQZyG5UKZq6/YZ1T+ye2NlnIyZR3aHbusjL/hnWgCqiE\nLyJjgUeALKCjMeaIh3LlSvjyvJ9+u3dU+dp58803M2DAAHr37k1aWhqDBw9m586dDB48mAsvvJC/\n/e1vbNu2jZtvvpmffvqJ6tWrn1WHJnylKkdWVha3dX6GjXVqQtvZsHQ8bEqA4m81HKVduxFl/qW7\nQFKlCV9E0oAGrouwPksaYYxZ6FLuGaCmMWa0h3qCrocPkJyczEcffcTHH39MYmIil112Gc8//zxr\n166lV69e7N69m2HDhpGVlUVycrLbOgK9jUoFq6XfLmVAygCytp/D0dkfQ24Ll7XBn+whwHr4xRWJ\nXAYsMsbEelhvRo0aVXw/Pj6e+Pj4gE+Ghw8fJjo6mh07dtC6dWvWrFlTPIxz1VVXMX36dHr37s2c\nOXO44YYb3NYR6G1UKtgcOX6Ep1KfYsXeFUzpMoUb699I+/Zj2bHj5eIyLVoMZe3a54Iu2aenp5Oe\nnl58//nnnw+MhC8izYwxu+zbQ4AOxpgeHsoGZQ8foEuXLuTn53PkyBG+/PLL4uUvvPAC7777LseP\nH2fnzp0eHx8MbVQqGBhjmJ05m6dSn+IPrf/A2FvHct6557F//35+97v5HD58KzExM9izZwAXXbSM\nDRu6B/1snUC6WuZ4EdkiIpuB24G/+CGmgFPyJw6L9OnTh61bt9KnTx+HIlMqfOzN2kuX2V2Y8PkE\nUnqm8Oqdr3LeuecB/vsN61Cn18OvgLy8PBo0aMDGjRtp2rSpx3LB3EalnJZfmM9ra17jxVUvMvSG\noTx1/VNUP+fMyREvvDCdu+++ibi4lsXLMjK2s3jxap59tn9Vh+xXATmG73VDIZjwJ02axKJFi1i6\ndGmp5YK5jUo5adOPm+i/sD91I+sy7Z5pNLugmdMhVTn9TdsAEBMTA8CH+oMUSvnd8dPHGbViFG9t\neYuXbn+JR+MeRcQ/U7TDmSb8ctqzZ4/TISgVklJ3pzLwo4Fcf9n1ZA7KpH5t/R1ff9GEr5QKCD/n\n/syTqU+y6rtVTOkyhbua3+V0SCFHr8aklHKUMYa3Mt6idXJr6teqz9ZBWzXZVxLt4SulHPPt0W8Z\n+NFADh8/zKJei/hto986HVJI04SvKk1BQQFL5s1j9b/+VXyp35sSEujcvbte6jfM5RfmM+mLSUxY\nPYFnbnyGJ65/gohqmo4qm07LrALh0MaSDh06xPNdu/JARgbxeXnFF2BKj4zkg7g4RqWkUL++fhgX\njjb8sIH+C/tzca2LmXrPVJrUa+J0SAEt5Ofh+6NnGEi9y3BL+IWFhQy54QYmrF1LbTfrc4Gn27fn\n9c8/155+GMk5lcPIFSOZnTmbiXdMpHeb3jrV0gf+TPgYY6rkz9rU2UouP3jwoBncvr1ZHhlpCsEY\nMIVglkdGmsHt25uDBw+6rcffdfiTp7aHqkXvv2+WR0YaY+97d3/LIiPN4rlznQ5VVZFFOxaZ6Fei\nTZ95fczPuT87HU5QsfOHf/KwvyryuiEfEn5BQYEZ3L69yfGQJHLADG7f3hQUFHjcOf6ow9/CLeEP\n79Kl+IXW018hmOF33+10qKqSHcw5aHp+0NPEvBpjluxa4nQ4QcmfCT+g3k8vmTePBzIy3A4DANQG\n7s/IILWUb7f6o44iMTExjB8/nlatWnHhhRfSr18/Tp065fVx4S7i+HG8vf8Uu5wKTcYYZm2aRWxy\nLI3rNCZzUCadmnZyOqywF1AJf9WsWcTn5ZVapmNeHitnzqzUOlzNnj2btLQ0du/ezTfffMPYsWN9\nelw4y69VC2+fWBi7nAo9u37Zxe1v387k9ZP55OFPmHDHBGqf66kLpqpSQCV8f/QM/d27HDJkCI0a\nNSIqKooRI0YwZ84cnx4Xzm5KSCA9MrLUMisiI+mQmFhFEamqcLrgNONXjee6f17H3c3vZu1ja2nb\nsK3TYSkXAZXw/dEz9HfvsnHjxsW3o6Oj+eGHH3x6XDjr3L07H8TFkethfS4wNy6OTm5+8F0Fp3UH\n1tFuRjvS96azvv96nrz+SZ1XH4ACKuH7o2fo797l999/X3x73759NGrUyKfHhbNq1aoxKiWFp9u3\nZ3lkZPELsAGWR0bydPv2jEpJ0SmZIeDXk7+S9EkS971zH8/c+AyLH15MTL0Yp8NSHgTUPHx/zN/2\n5xzwmJgY6tSpw6JFi6hZsyb33Xcf8fHx/O1vfyu9saW0MZwUFhayZP58Vs2aVfxdiA6JiXTq1k2T\nfQj4eMfHDF40mFtjbuXlO17mwloXOh1SSAqLefjLSsyhX1aOefgVqcMYY6644gozfvx407JlS1Ov\nXj2TkJBgTpw44dNjS2ujUsHsx19/ND3e72GavtbULN291OlwQh5+nJYZUD38Iv7oGfqjjpiYGN58\n801uvfVWn8p7Eq49fBVajDHM3DSTYcuG0a9tP0beMpKa1Ws6HVbIC/lLKwQKTfhKWXYc2cGAhQM4\nfvo4M+6dQdwlcU6HFDb8mfB1ILUUep0PFe5OFZxi3GfjuOHNG/j9b37PF/2+0GQfxHTeVCm+/fZb\np0NQyjFr9q+h/8L+XF73cr4c8CXRUdFOh6QqSBO+UuoM2SezGb5sOPO+mscrnV+hR6se+m43ROiQ\njlKqWMo3KbSa0ooTp0+wdfBWHmr9kCb7EKI9fKUUP/76I0MWD2HLwS281e0tOsZ0dDokVQkc7+FH\nR0cjIiH9Fx2tY58qMBWaQqZ/OZ02U9tw5YVXkjEwQ5N9CHN8WqZSyhlfH/6aAQsHcKrgFDPunUFs\ng1inQ1Ju6LRMpVS5nSo4xZhPx9BhVgd6tOrB6sTVmuzDhI7hKxVGVn+3mgEfDaBpvaZsHLCRy+pe\n5nRIqgr5JeGLyFPAROAiY8wv/qhTKeU/x/KOMWzZMBZ8s4DX7nyN+6+6X2ffhKEKD+mISGPgDmBf\nxcNRSvnb/K/m02pKK/IL89k6aCsPtHxAk32Y8kcP/xXgr0CKH+pSSvnJgewDDFk8hO0/b2f2/bO5\nOfpmp0NSDqtQwheRrsD3xpjMUO8xFBQUsGTePFb/61/FV9+8KSGBzt2767XdVUApNIVM2zCNkekj\nGdxuMLPvn01kROk/CqTCg9eELyJpQAPXRVg/XvQcMBxrOMd1nUejR48uvh0fH098fLzvkTro0KFD\nPN+1Kw9kZDA2L694B6QvX86Ql19mVEoK9evXdzpMpdj+83b6L+wPQPqj6bSq38rhiFRZpaenk56e\nXil1l3sevoi0BpYCx7ESfWPgAHCtMeaQm/JBOQ/fn7+gpVRlOZl/khdWvsCUDVMYEz+Gx9s9TjXR\n8zEUBMQ8fGPMVmPMJcaYJsaYGGA/0NZdsg9mS+bN44GMDLfJHqA2cH9GBqkffliVYSlVbOW+lcRN\njWPLoS1sfnwzg343SJO9csufZ4XBy5BOMFo1axbxeXmllumYl8fKmTOrKCKlLFl5WTy+8HF6zu3J\nC7e9wPyH5nNpnUudDksFML8lfLunH3Jz8COOH/f6KiZ2OaWqgjGGD7Z/QKspragm1dg2eBvdr+ru\ndFgqCOg3bb3Ir1XL61sXY5dTqrJ9f+x7/rT4T+w8spN3H3iXmy6/yemQVBDRgT4vbkpIID2y9Clt\nKyIj6ZCYWEURqXBUUFjA5HWTaTutLddccg2bHt+kyV6VmV4t0wudpaOctvXQVvov7E9EtQim3zOd\nqy6+yumQVBXy5ywdTfg+KJqHf39GBh1d5uGviIxkblyczsNXlWLCpGSyrz7AtC+nMbbjWPr/tr/O\nvglDmvAdUFhYyJL581k1a1bxN207JCbSqVs37dkrv5u7YS493u7FndfczozuM2h0fiOnQ1IO0YSv\nVIjKPZXLsGXDmLX2LXLeHc7w7rUYN+5PToelHKQJX6kQM2bMG8xZs4K9bZZRNyuaWis7s2f7SzRr\nNppL7an1eXnH6NKlGSNH/tHRWFXV8mfC12mZSjks+2Q2e1ptYPfxTzg9703ydj5UvG7XrtHs2gVw\nlHbtRvDnPz/sWJwq+Ongs1IOWrJrCbHJsURUj2DnE9toV/czIKtEKSvZp6W9QFRUlBNhqhChQzpK\nOeDoiaM8lfoUy/csZ8a9M7ijqXXR2aysLNq3H8uOHS8Xl23RYihr1z6nyT5MBcTF05RS5bPwm4XE\nJsdSM6ImmYMyi5M9QE5ODtnZ0UREbKN58yQiIraTnR1Nbm6ugxGrUKEJX6kqcuT4EXrP603SkiT+\n0/0/vHH3G5xf4/wzykyePJeCgn0kJaWRmfkSSUmpFBTs44035joUtQolmvCVqgJzt88lNjmWi2td\nzJaBW4i/It5tuTp1apKWlsjEiUnUqFGDiROTSEtLpE6dmlUbsApJOoavVCU6lHuIPy36ExkHM5jZ\ndSY3Xn6j0yGpIKNj+EoFOGMMczLn0Ca5DTFRMWx+fLMme+U4nYevlJ/9+OuPDPp4EDt/2UlKzxSu\nvfRap0NSCtAevlJ+Y4zh35v/TdzUOGLrx7JxwEZN9iqgaA9fKT/Yn72fxz96nAPZB1jSewltG7Z1\nOiSlzqI9fKUqwBjDjC9n0HZaW6679DrW9V+nyV4FLO3hK1VOe7P20n9hf46eOMryR5YT2yDW6ZCU\nKpX28JUqo0JTyBvr3qDd9HbcFnMbax5bo8leBQXt4StVBrt+2cVjKY9xsuAkKxNW6s8NqqCiPXyl\nfFBQWMCra17lun9eR9cru7IqYZUmexV0tIevlBdfH/6afin9OEfO4Yt+X9D8wuZOh6RUuWgPXykP\n8gvzmbB6AjfNvImerXuS3jddk70KatrDV8qNrYe2krggkfNrnM/6/uuJqRfjdEhKVZj28JVycbrg\nNGM/G0vUhVywAAARLElEQVTHf3fksWseY2mfpZrsVcjQHr5Sts0/bSZhQQKXnHcJGwds5LK6lzkd\nklJ+VaEevoiMEpH9IrLR/rvTX4EpVVVOFZxi5IqRdHq7E39p/xcW9VqkyV6FJH/08CcZYyb5oR6l\nqtz6A+tJTEm0LmE8cDONzm/kdEhKVRp/JHy/XJhfqaqUl5/H6PTRzNo8i1c6v0LP1j0R0VNZhTZ/\nfGj7RxHZLCL/FJG6fqhPqUr1+fefc/XUq9n1yy62DNxCr9hemuxVWPDawxeRNKCB6yLAACOAKcAY\nY4wRkbHAJKCfp7pGjx5dfDs+Pp74+PhyBa1UeRw/fZznlj/HnK1zeP2u13mg5QNOh6TUWdLT00lP\nT6+Uuv32m7YiEg0sNMa08bBef9NWOebTvZ/SL6Uf7Ru357U7X+OiWhc5HZJSPvHnb9pWaAxfRC4x\nxvxk3+0ObK14SEr5T86pHJ5d+izzv55P8t3JdL2yq9MhKeWYin5oO0FErgYKgb3A4xWOSCk/Wfrt\nUvov7E/8FfFsHbSVejXrOR2SUo7y25CO1w3pkI6qIsfyjvHXtL/yya5PmHbPNO5qfpfTISlVbv4c\n0tFLK6iQsnjnYmKTYxGEzEGZmuyVcqGXVlAh4eiJozyx5Ak+3fcps+6bxW1NbnM6JKUCjvbwVdBb\n8PUCWie35vxzzydzUKYme6U80B6+ClqHjx/mz4v/zPof1jPn/jncHH2z0yEpFdC0h6+C0vvb3ic2\nOZaG5zUkY2CGJnulfKA9fBVUDuYc5I+L/si2n7cxr8c8rr/seqdDUipoaA9fBQVjDP/d8l/aTG1D\n8wuas+nxTZrslSoj7eGrgPfDrz8w8KOB7Mnaw8e9PqZdo3ZOh6RUUNIevgpYxhhmbZrF1VOvpu0l\nbdnQf4Mme6UqQHv4KiB9d+w7Hv/ocX7K+YnUPqlcfcnVToekVNDTHr4KKMYYpm2Yxm+n/5abLruJ\ndY+t02SvlJ9oD18FjD1H9/DYwsf49eSvpD+aTqv6rZwOSamQoj185bhCU8jkdZP53Yzf0blpZz7v\n97kme6UqgfbwlaN2HtlJv5R+FJgCVieu5sqLrnQ6JKVClvbwlSMKCguY9MUkrn/zerpf1Z3P+n6m\nyV6pSqY9fFXlvvr5KxJTEqlxTg3WPraWphc0dTokpcKC9vBVlckvzOfFlS/SYVYH+rTpw/JHl2uy\nV6oKaQ9fVYm/vjyaFXU+ol7NemwYsIEroq5wOiSlwo728FWlOl1wmqELh/L3n8fyUNOHSO2dqsle\nKYdoD19VmoyfMui7oC/ZP5zETJ1PVsQ+5Da//DSnUqoc9EfMld+NfP41Zu58h0PRW2iy+3byv7ya\n3buep1mz0Vx6qVUmL+8YXbo0Y+TIPzoaq1KBzp8/Yq4JX/nV5p8202duH77fls2x/y6CbHdfoDpK\nu3YjSEt7gaioqCqPUalg4s+Er2P4yi9OFZxidPpoOr3diaE3DmXPuM20azEFyCpRUpO9Uk7RHr6q\nsE0/bqLvgr5cVucypt0zjUvrWOM2WVlZtG8/lh07Xi4u26LFUNaufU6TvVI+0h6+CginCk4xcsVI\nOv+nM09d/xQLey4sTvYAOTk5ZGdHExGxjebNk4iI2E52djS5ubkORq1U+NKEr8pl448baTe9HZt+\n2sTmgZt5JO4RRM7shEyePJeCgn0kJaWRmfkSSUmpFBTs44035joUtVLhTadlqjI5mX+SsZ+NZdqX\n0/h7p7/Tu03vsxJ9kTp1apKWlkhcXEsAJk5Monfv7SxevLoqQ1ZK2XQMX/nsyx++pO+CvsRExTDt\nnmk0PL+h0yEpFfICagxfRIaIyFcikiki4/0RlAosJ/NP8tzy5+gyuwvP3vgsC/6wQJO9UkGoQkM6\nIhIP3AvEGmPyReQiv0SlAsaGHzbQ98O+NL2gKZsf36yJXqkgVtEx/EHAeGNMPoAx5nDFQ1KB4GT+\nSZ7/9Hne3PQmr3R+hZ6te3ocq1dKBYeKDum0AG4WkTUiskJE2vkjKOWs9QfWc830a/jq8FdkDMyg\nV2wvTfZKhQCvPXwRSQMauC4CDPCc/fh6xpjrROR3wHtAE091jR49uvh2fHw88fHx5QpaVY68/Dye\nT3+emZtn8mrnV/lD6z9ooleqiqWnp5Oenl4pdVdolo6ILAJeMsZ8at/fBbQ3xhxxU1Zn6QSwdQfW\nkbAggSsvvJLku5NpcF4D7w9SSlU6f87SqegY/ofArcCnItICqO4u2avAlZefx+j00fxr87947c7X\n6NGqh/bqlQpRFU34s4CZIpIJnAQeqXhIqqqs3b+Wvgv60vLilmQMzNBevVIhTr94FYby8vMYuWIk\nb2W8xT/u+gcPtnxQe/VKBahAGtJRQWbN/jUkLEigdf3WbBm0hfq16zsdklKqimjCDxMnTp9g5IqR\nvL3lbV6/63UebPWg0yEppaqYJvww8MX3X5CwIIG4S+LIHJTJxbUvdjokpZQDNOGHsBOnT/B/K/6P\n/2b+l9fvep0HWj7gdEhKKQdpwg9Rq79bTWJKIm0vacuWgVu0V6+U0oQfao6fPs5zy59jztY5TL5r\nMve3vN/pkJRSAUITfghZ9d0qEhck8ttGvyVzUCYX1dKLlyql/kcTfgg4fvo4I5aN4N1t7zK5y2S6\nX9Xd6ZCUUgFIf9M2yK3ct5K4qXEcOn6IzEGZmuyVUh5pDz9I5Z7KZcTyEby37T2m3D2Fbr/p5nRI\nSqkApz38IPTZvs+ImxrHkRNHyByUqcleKeUT7eEHkdxTuQxbNoy5X80l+e5kul7Z1emQlFJBRHv4\nQeLTvZ/SZmobsvKyyByUqcleKVVm2sMPcDmnchi2dBjzv55P8t3J3HvlvU6HpJQKUtrDD2Dpe9OJ\nmxpH9qlsMgdlarJXSlWI9vADUM6pHJ5d+iwffv0hU++Zyj0t7nE6JKVUCNAefoBZsWcFbZLbkHs6\nl8xBmZrslVJ+oz38AJFzKoen055m4Y6FTLtnGl2ad3E6JKVUiNEefgBYvmc5scmx5OXnkTkoU5O9\nUqpSaA/fQb+e/JWn057mo50fMf2e6dzV/C6nQ1JKhTDt4Ttk2bfLiE2O5VTBKbYO2qrJXilV6bSH\nX8WyT2bzdNrTLNq5iOn3TufOZnc6HZJSKkxoD78Kpe1Oo01yGwoKC8gclKnJXilVpbSHXwWyT2Yz\nNHUoS3YvYfo90+ncrLPTISmlwpD28CtZ6u5UYpNjAdgycIsme6WUY7SHX0mO5R1jaOpQUr9NZca9\nM+jUtJPTISmlwpz28CvBkl1LiE2OpZpUI3NQpiZ7pVRA0B6+Hx3LO8ZTqU+x9NulzLxvJrc3ud3p\nkJRSqliFevgi8o6IbLT/9ojIRn8FFmwW71xMbHIs1atVJ3NQpiZ7pVTAqVAP3xjzh6LbIvIykFXh\niIJMVl4WTy55khV7VzDrvlnc1uQ2p0NSSim3/DmG3wOY48f6At6inYuITY4lMiKSLQO3aLJXSgU0\nv4zhi0gH4CdjzG5/1BfosvKyeGLJE6TvTeff3f7NrTG3Oh2SUkp55bWHLyJpIrLF5S/T/u/680s9\nCZPe/cc7PiZmQhNqRdQic1CmJnulVNDw2sM3xtxR2noROQfoDlzjra7Ro0cX346Pjyc+Pt5rgIFm\n/pb5nH4viuH9hnPeuec5HY5SKsSkp6eTnp5eKXWLMaZiFYjcCTxjjOnopZyp6LYCwfDhr/Pii+0Z\nPnwd48b9yelwlFIhTkQwxog/6vLHGP5DhPBwzpgxb7Bo0S4iI+sCcOAAwLW8994iVq8eDUBe3jG6\ndGnGyJF/dCxOpZTypsI9fJ83FKQ9/KysLO64YwQbNowDotyUOEq7diNIS3uBqCh365VSqvz82cPX\nSyt4ERUVRVraONq1G8HZXzPQZK+UCh7aw/dRVlYW7duPZceOl4uXtWgxlLVrn9Nkr5SqNNrDd0BO\nTg7Z2dFERGyjefMkIiK2k50dTW5urtOhKaWUTzTh+2jy5LkUFOwjKSmNzMyXSEpKpaBgH2+8Mdfp\n0JRSyid6tUwf1alTk7S0ROLiWgIwcWISvXtvZ/Hi1Q5HppRSvtExfKWUCmA6hq+UUqrMNOErpVSY\n0ISvlFJhQhO+UkqFCU34SikVJjThK6VUmNCEr5RSYUITvlJKhQlN+EopFSY04SulVJjQhK+UUmFC\nE75SSoUJTfhKKRUmNOErpVSY0ISvlFJhQhO+UkqFCU34SikVJjThK6VUmNCEr5RSYUITvlJKhQlN\n+EopFSY04SulVJioUMIXkTgR+UJENonIOhFp56/AlFJK+VdFe/gTgFHGmLbAKGBixUMKTunp6U6H\nUKlCuX2h3DbQ9qn/qWjCLwTq2rejgAMVrC9ohfpJF8rtC+W2gbZP/U9EBR//BLBERP4OCHBDxUNS\nSilVGbwmfBFJAxq4LgIMMAK4HfiLMeZDEXkAmAncURmBKqWUqhgxxpT/wSJZxpgol/vHjDF1PZQt\n/4aUUiqMGWPEH/VUdEjngIjcYoz5VERuA3Z4KuivgJVSSpVPRRN+f+AfInIOkAcMqHhISimlKkOF\nhnSUUkoFj0r9pq2IDBGRr0QkU0TGeyhzp4h8LSI7ROSZyozHn0RklIjsF5GN9t+dHsrtFZGMoi+n\nVXWc5VWG9gXl8QMQkadEpFBELvCwvsBu+yYR+bCq46soH9r3qH3cvhGRR6o6vvISkTEuz6lPROQS\nD+WC8viVoX1lP37GmEr5A+KBVCDCvn+RmzLVgF1ANFAd2Az8prJi8nP7RgFP+lDuW6Ce0/FWRvuC\n/Pg1Bj4B9gAXeCiT7XScldU+oB6wG+t7NFFFt52O28e2nedyewiQHErHz5f2lff4VWYPfxAw3hiT\nD2CMOeymzLXATmPMPmPMaeAd4L5KjMnffPkgWgjeaxZ5a18wH79XgL96KRPMEw28ta8zkGqMOWaM\nycLqnLl9FxdojDE5LndrY30B1J2gPH4+tq9cx68yE1EL4GYRWSMiKzxcZ+dS4HuX+/vtZcHijyKy\nWUT+KSJup6NifWdhiYisF5H+VRmcH3hrX1AePxHpCnxvjMn0UrSGfY2oz0UkWF7IfG1fyWN3gCA4\ndkVEZKyIfAf0AkZ6KBaUxw98al+5jl+FZumU8qWs5+y66xljrhOR3wHvAU0qsr2q5uVLZ1OAMcYY\nIyJjgUlAPzfV3GiM+VFELgbSROQrY8yqyo7dF35qX0Dycm4O58wvCHrqCUbbxy4GWC4iW4wxeyol\n4DLyU/sCVmnnpjFmoTHmOeA5+3OjIcBoN9UE4/ErS/vKrEIJ3xjj8Vu1IjIQmGeXW29/eHShMeaI\nS7EDwOUu9xsTQNfjKa19JcwAFnqo40f7/88iMh9rGCQgEr4f2hewx89T20SkNXAFkCEighXzlyJy\nrTHmUIk6io7dHhFJB9pijYk7zg/tO4D1OVuRxsCKyom27Mpwbs4GFuEmIQbj8XPDU/vKdfwqc0jn\nQ+BWABFpAVQvkewB1gPNRCRaRM4F/gCkVGJMflPik/PuwFY3ZWqJyHn27dpAJ3flApEv7SMIj58x\nZqsx5hJjTBNjTAzWMFTbksleRKLsNiEiF2FdJ2p71UdcNr62D1gC3CEidUWkHtY7giVVHW95iEgz\nl7vdgK/clAnK4we+tY/yHr9K/KS5OvA2kAlsAG6xlzcEPnIpdyfwDbATeNaJT8XL2b63gC1YM1M+\nBBqUbB8QY6/fZO+HkGpfMB8/l/i/xZ7FAvwWmG7fvt5u/yYgA+jrdKz+bJ99v6993HYAjzgdaxna\n9IHLubkAaBhKx8+X9pX3+OkXr5RSKkwE63RBpZRSZaQJXymlwoQmfKWUChOa8JVSKkxowldKqTCh\nCV8ppcKEJnyllAoTmvCVUipM/D/Aqx9yOHhpBQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5dd9704390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if uw.rank()==0:\n",
    "    uw.matplotlib_inline()\n",
    "    import matplotlib.pyplot as pyplot\n",
    "    import matplotlib.pylab as pylab\n",
    "    pyplot.plot(h, yvx, '*', label='vx', markersize=10)\n",
    "    pyplot.plot(h, yvy, '-', label='vy', markersize=10)\n",
    "    pyplot.plot(h, yp, 'o', label='p', markersize=10)\n",
    "    pyplot.legend(loc=2)\n",
    "    pyplot.title(\"Log graph of relative error (L2 norm) vs element width\")\n",
    "    pyplot.savefig(\"solCx\") if uw.nProcs() > 1 else pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# figViscosity = glucifer.Figure()\n",
    "# figForce = glucifer.Figure(**figViscosity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# figViscosity.append( glucifer.objects.Surface(mesh, fn=sol.fn_viscosity, onMesh=True) )\n",
    "# figViscosity.append( glucifer.objects.Mesh(mesh))\n",
    "# figViscosity.show()\n",
    "# figForce.append( glucifer.objects.Surface(mesh, fn=sol.fn_bodyforce[1], onMesh=True) )\n",
    "# figForce.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# figSR = glucifer.Figure(**figViscosity)\n",
    "# figSR.append( glucifer.objects.Surface(mesh, fn=pressureField ) )\n",
    "# figSR.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}