weak-formulation.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "title: Weak Formulation\n",
        "description: The weak formulation reduces requirements of differentiability on the original equation allowing us to discretize tensor anisotropy or curvilinear meshes.\n",
        "---"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "XRXnNzNyUsLV40Qxc5IJ",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 2
          }
        }
      },
      "source": [
        "Equation (b) is a vector equation, so really it is two or three equations involving multiple components of $\\vec{j}$. We want to work with a single scalar equation, allow for anisotropic physical properties, and potentially work with non-axis-aligned meshes — how do we do this?! We can use the **weak formulation** where we take the inner product ($\\int \\vec{a} \\cdot \\vec{b} dv$) of the equation with a generic face function, $\\vec{f}$. This reduces requirements of differentiability on the original equation and also allows us to consider tensor anisotropy or curvilinear meshes (Haber 2014).\n",
        "\n",
        "## Scalar equations only, please\n",
        "\n",
        "In Figure 5, we visually walk through the discretization of equation (b). On the left hand side, a dot product requires a *single* cartesian vector, $[\\mathbf{j_x, j_y}]$. However, we have a $j$ defined on each face (2 $j_x$ and 2 $j_y$ in 2D!). There are many different ways to evaluate this inner product: we could approximate the integral using trapezoidal, midpoint or higher order approximations. A simple method is to break the integral into four sections (or 8 in 3D) and apply the midpoint rule for each section using the closest $\\mathbf{j}$ components to compose a cartesian vector. A $\\mathbf{P}_i$ matrix (size 2 × 4) is used to pick out the appropriate faces and compose the corresponding vector (these matrices are shown with colors corresponding to the appropriate face in the figure). On the right hand side, we use a vector identity to integrate by parts. The second term will cancel over the entire mesh (as the normals of adjacent cell faces point in opposite directions) and $\\phi$ on mesh boundary faces are zero by the Dirichlet boundary condition. This leaves us with the divergence, [which we already know how to do](divergence.ipynb)!"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "mdDk7ewGdTdUKlz8g7yg",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "\n",
        "```{figure} ./images/weak-formulation.png\n",
        ":name: b8cf04cf\n",
        "\n",
        "Discretization using the weak formulation and inner products.\n",
        "```"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "E0m5Jh7cTFESjbe0zrGb",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 2
          }
        }
      },
      "source": [
        "The final step is to recognize that, now discretized, we can cancel the general face function $\\mathbf{f}$ and transpose the result (for convention's sake):"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "k3ce9vtzui0ybKCPtcIK",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 3
          }
        }
      },
      "source": [
        "$$\\frac{1}{4} \\sum_{i=1}^4 \\mathbf{P}_i^\\top \\sqrt{v} \\mathbf{\\Sigma}^{-1} \\sqrt{v} \\mathbf{P}_i \\mathbf{j} = \\mathbf{D}^\\top v \\phi$$"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "GAbpJlNVTjUM0Z5vtT50",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "# Implementation\n",
        "\n",
        "We will start by using a SimPEG Mesh that only has a single cell."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": true,
        "iooxa": {
          "id": {
            "block": "YzzyeWqgAuqxx0PAoA4q",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "P9oBPKfsnXyc0cZTOYEe",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "%matplotlib inline\n",
        "import discretize\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "hb8q7AI0MaxLo3uw6PuF",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "jXh403rknR7GhoyOIml3",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEKCAYAAAAmfuNnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZO0lEQVR4nO3df3Bd9X3m8fdj2QSIIUAswLFNTFujXcoYSqhNClPkdsnYniYmm0xib4akKamGTLxtupvOejazJBt2Z5pktmEyAYQ39UCyBIZSSLypY0NZFG/LjxqIbWyIwDHOoNiJZUAGBS/I0mf/uEf0Il/pXknnfs+59vOauaPz63vPo6MjPz7nXl8rIjAzM2u2GUUHMDOzE4MLx8zMknDhmJlZEi4cMzNLwoVjZmZJuHDMzCyJQgtH0gZJByXtGmd9p6TDkrZnjxuq1i2X1Ctpj6R16VKbmdlUqMh/hyPp94FB4DsRcVGN9Z3AFyLij8YsbwOeA64G+oBtwJqIeKbpoc3MbEoKvcKJiK3Ay1MYugTYExF7I+JN4G5gVa7hzMwsVzOLDtCA90vaAeyncrWzG5gHvFi1TR+wtNZgSV1AF8DJJ5/8vvPOO6/JcadvZGSEGTPK//Kac+anFTKCc+atVXI+99xzhyKifbrPU/bCeQp4b0QMSloJfB9YBKjGtjXvDUbEemA9QEdHR/T29jYra256enro7OwsOkZdzpmfVsgIzpm3Vskp6ed5PE+pqzUiXo2IwWx6EzBL0hwqVzQLqjadT+UKyMzMSqrUhSPpXEnKppdQyfsSlTcJLJJ0vqSTgNXAxuKSmplZPYXeUpN0F9AJzJHUB3wJmAUQEd3AR4HPSjoKHAFWR+VtdUclrQW2AG3Ahuy1HTMzK6lCCyci1tRZ/y3gW+Os2wRsakYuMzPLX6lvqZmZ2fHDhWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkkUWjiSNkg6KGnXOOs/IWln9nhE0sVV6/ZJelrSdklPpEttZmZTUfQVzu3A8gnWvwBcFRGLgRuB9WPWL4uISyLisiblMzOznMwscucRsVXSwgnWP1I1+xgwv9mZzMysOYq+wpmM64AfVc0H8ICkJyV1FZTJzMwapIgoNkDlCueHEXHRBNssA24BroyIl7Jl74mI/ZLOBh4E/n1EbK0xtgvoAmhvb3/fPffck/83kbPBwUFmz55ddIy6nDM/rZARnDNvrZJz2bJlT+by0kVEFPoAFgK7Jli/GPgZcMEE23wZ+EK9fV1wwQXRCh5++OGiIzTEOfPTChkjnDNvrZITeCJy+PO+1LfUJJ0H3AdcGxHPVS1/p6TTRqeBDwA13+lmZmblUOibBiTdBXQCcyT1AV8CZgFERDdwA/Bu4BZJAEejcll3DnB/tmwm8L2I2Jz8GzAzs4YV/S61NXXWfwb4TI3le4GLjx1hZmZlVepbamZmdvxw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS6LQwpG0QdJBSbvGWS9J35S0R9JOSZdWrVsuqTdbty5darN89b/ez02/vIlDRw4VHcWsqYq+wrkdWD7B+hXAouzRBdwKIKkNuDlbfyGwRtKFTU1q1iTdO7vZ+8Zeund0Fx3FrKlmFrnziNgqaeEEm6wCvhMRATwm6QxJc4GFwJ6I2Asg6e5s22eam9hGLVz39wAs7X204CT1DQwc4daS5hwaGWJn/5kEf8rtPxdP/WQrs2bMKjrWuMp8LEc9/sLLAOzrLDaHHavQwmnAPODFqvm+bFmt5UtrPYGkLipXR7S3t9PT09OUoHkaHBxsiZwAAwMDRUeoa3h4uLQ5+4/2EwQAQbDv5X20z2wvONX4ynwsx2qF36FW+l3PQ9kLRzWWxQTLj10YsR5YD9DR0RGdnZ25hWuWnp4eyp5zae+jDAwMsOU/rSg6Sl1lPZ79r/ez4r4VnDr8xlvLZrS9gzs/spk5p8wpMNn4ynosq338tsq5Wfac0BrHM09Fv4ZTTx+woGp+PrB/guVmLaN7ZzcjMfK2ZSMx4tdy7LhV9sLZCHwye7fa5cDhiDgAbAMWSTpf0knA6mxbs5ax4+AOhkaG3rZsaGSI7Qe3F5TIrLkKvaUm6S6gE5gjqQ/4EjALICK6gU3ASmAP8Drw6WzdUUlrgS1AG7AhInYn/wbMpuHeD9371vSJdmvFTkxFv0ttTZ31AXxunHWbqBSSmZm1gLLfUjMzs+OEC8fMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJIotHAkLZfUK2mPpHU11v+lpO3ZY5ekYUlnZev2SXo6W/dE+vRmZjYZM4vasaQ24GbgaqAP2CZpY0Q8M7pNRHwd+Hq2/QeBv4iIl6ueZllEHEoY28zMpqjIK5wlwJ6I2BsRbwJ3A6sm2H4NcFeSZGZmlrvCrnCAecCLVfN9wNJaG0o6FVgOrK1aHMADkgK4LSLWjzO2C+gCaG9vp6enZ/rJm2xwcLD0OQcGjjA8PFz6nNAax7MVMkJr5PS5WV5FFo5qLItxtv0g8E9jbqddERH7JZ0NPCjppxGx9ZgnrBTReoCOjo7o7OycZuzm6+npoew5b+19lIGBgdLnhNY4nq2QEVojp8/N8iryllofsKBqfj6wf5xtVzPmdlpE7M++HgTup3KLzszMSqrIwtkGLJJ0vqSTqJTKxrEbSXoXcBXwg6pl75R02ug08AFgV5LUZmY2JYXdUouIo5LWAluANmBDROyWdH22vjvb9MPAAxHx66rh5wD3S4LK9/C9iNicLr2ZmU1Wka/hEBGbgE1jlnWPmb8duH3Msr3AxU2OZ2ZmOfInDZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS2LCwpF0uqTfrLF8cR47l7RcUq+kPZLW1VjfKemwpO3Z44ZGx1pzfHTjR7nxsRsZGhkqOsrxYec98I2LuKrnGvjGRZV5mxKfm+U3c7wVkj4G3AQclDQL+OOI2Jatvh24dDo7ltQG3AxcDfQB2yRtjIhnxmz6fyPij6Y41nLW+0ovew/v5dX+Mzm97XT6X++n/dT2omO1pp33wP/+Mxg6ggAOv1iZB1j8sSKTtSSfm+U30RXOfwbeFxGXAJ8Gvivp32brlMO+lwB7ImJvRLwJ3A2sSjDWpmloZIggODx8mBX3reDGx26k//X+omO1noe+AkNH3r5s6EhluU2Jz81yG/cKB2iLiAMAEfHPkpYBP5Q0H4gc9j0PeLFqvg9YWmO790vaAewHvhARuycxFkldQBdAe3s7PT0900/eZIODgy2RE+CN4Tf4296/5cl9T/L5cz9fdJyayno8rzrcV/NvbnG4jx+XMC+U91jW4nOzfCYqnNck/WZE/AwgIg5I6gS+D/x2Dvuu+bs2Zv4p4L0RMShpZbbvRQ2OrSyMWA+sB+jo6IjOzs4pB06lp6eH0ua84+2z72h7B9f81jVcf/H1zDllTjGZ6ijt8fzJ/MpttDH0rvnlzEuJjyX43GwBExXOZ4EZki4cfW0kIl6TtBxYncO++4AFVfPzqVzFvCUiXq2a3iTpFklzGhlrzTNrxiyEOL3tdDZ/ZHNpf5lL7w9veOs1nLfMOqWy3KbE52a5jVs4EbEDQNIuSd8FvgacnH29DPjuNPe9DVgk6XzgF1RK7N9VbyDpXOBXERGSllB5zeklYKDeWGuOjjM7uOTsS3hqcDG/fvXX/oWejtE3Bjz0FeJwH3rX/ErZ+A0DU+Jzs/wmusIZtRT4KvAIcBpwJ3DFdHccEUclrQW2AG3AhojYLen6bH038FHgs5KOAkeA1RERQM2x081k9d37oXsB+PiORwtOcpxY/DFY/DF+fILdWmkGn5vl10jhDFH5w/4UKlc4L0TESB47j4hNwKYxy7qrpr8FfKvRsWZmVl6NfNLANiqF87vAlcAaSfc2NZWZmR13GrnCuS4insimfwmsknRtEzOZmdlxqO4VTlXZVC+b7hsGzMzsBOMP7zQzsyRcOGZmloQLx8zMknDhmJlZEi4cMzNLwoVjZmZJuHDMzCwJF46ZmSXhwjEzsyRcOGZmloQLx8zMknDhmJlZEi4cMzNLwoVjZmZJuHDMzCyJQgtH0nJJvZL2SFpXY/0nJO3MHo9Iurhq3T5JT0vaLumY/7PHzMzKpZH/8bMpJLUBNwNXA33ANkkbI+KZqs1eAK6KiFckrQDWA0ur1i+LiEPJQpuZ2ZQVeYWzBNgTEXsj4k3gbmBV9QYR8UhEvJLNPgbMT5zRzMxyUtgVDjAPeLFqvo+3X72MdR3wo6r5AB6QFMBtEbG+1iBJXUAXQHt7Oz09PdPJnMTg4GDpcw4MHGF4eLj0OaE1jmcrZITWyOlzs7yKLBzVWBY1N5SWUSmcK6sWXxER+yWdDTwo6acRsfWYJ6wU0XqAjo6O6OzsnHbwZuvp6aHsOW/tfZSBgYHS54TWOJ6tkBFaI6fPzfIq8pZaH7Cgan4+sH/sRpIWA98GVkXES6PLI2J/9vUgcD+VW3RmZlZSRRbONmCRpPMlnQSsBjZWbyDpPOA+4NqIeK5q+TslnTY6DXwA2JUsuZmZTVpht9Qi4qiktcAWoA3YEBG7JV2fre8GbgDeDdwiCeBoRFwGnAPcny2bCXwvIjYX8G2YmVmDinwNh4jYBGwas6y7avozwGdqjNsLXDx2uZmZlZc/acDMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRIuHDMzS8KFY2ZmSbhwzMwsCReOmZkl4cIxM7MkXDhmZpaEC8fMzJJw4ZiZWRKFFo6k5ZJ6Je2RtK7Gekn6ZrZ+p6RLGx1r1ir6X+/npl/exKEjh4qOYtZUhRWOpDbgZmAFcCGwRtKFYzZbASzKHl3ArZMYa9YSund2s/eNvXTv6C46illTzSxw30uAPRGxF0DS3cAq4JmqbVYB34mIAB6TdIakucDCBsZaEz3+wssAfPy2RwtOUt/AwBFu7S1nzqGRIXb2n0nwp9z+c/HUT7Yya8asomONq8zHctTouWnlU2ThzANerJrvA5Y2sM28BscCIKmLytUR7e3t9PT0TCt0CoODgy2RE2BgYKDoCHUNDw+XNmf/0X6CACAI9r28j/aZ7QWnGl+Zj+VYrfA71Eq/63kosnBUY1k0uE0jYysLI9YD6wE6Ojqis7NzEhGL0dPTQ9lz7utsjZxQ3pz9r/ez4r4VnDr8xlvLZrS9gzs/spk5p8wpMNn4ynosx3LOciryTQN9wIKq+fnA/ga3aWSsWal17+xmJEbetmwkRvxajh23iiycbcAiSedLOglYDWwcs81G4JPZu9UuBw5HxIEGx5qV2o6DOxgaGXrbsqGRIbYf3F5QIrPmKuyWWkQclbQW2AK0ARsiYrek67P13cAmYCWwB3gd+PREYwv4Nsym7N4P3fvW9Il2a8VOTEW+hkNEbKJSKtXLuqumA/hco2PNzKy8/EkDZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMySKKRwJJ0l6UFJz2dfz6yxzQJJD0t6VtJuSX9ete7Lkn4haXv2WJn2OzAzs8kq6gpnHfBQRCwCHsrmxzoK/MeI+NfA5cDnJF1Ytf4bEXFJ9tjU/MhmZjYdRRXOKuCObPoO4JqxG0TEgYh4Kpt+DXgWmJcsoZmZ5aqowjknIg5ApViAsyfaWNJC4HeAx6sWr5W0U9KGWrfkzMysXBQRzXli6R+Ac2us+iJwR0ScUbXtKxFRszQkzQZ+DPz3iLgvW3YOcAgI4EZgbkT8yTjju4AugPb29vfdc889U/+mEhkcHGT27NlFx6jLOfPTChnBOfPWKjmXLVv2ZERcNu0niojkD6CXSkkAzAV6x9luFrAF+A8TPNdCYFcj+73ggguiFTz88MNFR2iIc+anFTJGOGfeWiUn8ETk8Gd/UbfUNgKfyqY/Bfxg7AaSBPwN8GxE/PWYdXOrZj8M7GpSTjMzy0lRhfNXwNWSngeuzuaR9B5Jo+84uwK4FviDGm9//pqkpyXtBJYBf5E4v5mZTdLMInYaES8Bf1hj+X5gZTb9j4DGGX9tUwOamVnu/EkDZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMyScOGYmVkSLhwzM0vChWNmZkm4cMzMLAkXjpmZJeHCMTOzJFw4ZmaWhAvHzMyScOGYmVkShRSOpLMkPSjp+ezrmeNst0/S05K2S3pisuPNzKw8irrCWQc8FBGLgIey+fEsi4hLIuKyKY43M7MSKKpwVgF3ZNN3ANckHm9mZokpItLvVBqIiDOq5l+JiGNui0l6AXgFCOC2iFg/mfHZui6gK5u9CNiV33fSNHOAQ0WHaIBz5qcVMoJz5q1VcnZExGnTfZKZeSSpRdI/AOfWWPXFSTzNFRGxX9LZwIOSfhoRWyeTIyup0aJ6YsytuVJyzny1Qs5WyAjOmbdWypnH8zStcCLi34y3TtKvJM2NiAOS5gIHx3mO/dnXg5LuB5YAW4GGxpuZWXkU9RrORuBT2fSngB+M3UDSOyWdNjoNfIB/uR1Wd7yZmZVLUYXzV8DVkp4Hrs7mkfQeSZuybc4B/lHSDuCfgb+PiM0TjW/A+ry+gSZzzny1Qs5WyAjOmbcTKmchbxowM7MTjz9pwMzMknDhmJlZEsdd4bTKx+Y0sh9JCyQ9LOlZSbsl/XnVui9L+kWWf7uklTlmWy6pV9IeScd8ioMqvpmt3ynp0kbH5qmBnJ/I8u2U9Iiki6vW1fz5F5SzU9Lhqp/lDY2OTZzzL6sy7pI0LOmsbF2S4ylpg6SDkmr+e7oSnZv1cpbl3KyXM99zMyKOqwfwNWBdNr0O+Oo42+0D5kx1fIqcwFzg0mz6NOA54MJs/svAF5qQqw34GfAbwEnAjtF9Vm2zEvgRIOBy4PFGxybO+XvAmdn0itGcE/38C8rZCfxwKmNT5hyz/QeB/1PA8fx94FJg1zjrCz83G8xZ+LnZYM5cz83j7gqH1vnYnLr7iYgDEfFUNv0a8Cwwr0l5Ri0B9kTE3oh4E7g7y1ptFfCdqHgMOEOVfw/VyNhkOSPikYh4JZt9DJjfpCwTmc4xKdXxHGMNcFeTsowrKv/w++UJNinDuVk3Z0nOzUaO53imdDyPx8I5JyIOQOUPbODscbYL4AFJT6ry8TeTHZ8qJwCSFgK/AzxetXhtdkm+Icdbf/OAF6vm+zi25MbbppGxeZnsvq6j8jffUeP9/PPWaM73S9oh6UeSfnuSY/PQ8L4knQosB/6uanGq41lPGc7NySrq3GxUbudm0z5poJlUko/NqSennEiaTeWX+/MR8Wq2+FbgRion543A/wD+ZOpp/2V3NZaNfe/8eNs0MjYvDe9L0jIqv9RXVi1u+s9/EjmfAt4bEYPZa3HfBxY1ODYvk9nXB4F/iojqvxmnOp71lOHcbFjB52Yjcj03W7JwokU+NiePnJJmUSmbOyPivqrn/lXVNv8T+OFUc47RByyomp8P7G9wm5MaGJuXRnIiaTHwbWBFRLw0unyCn3/ynFV/iSAiNkm6RdKcRsamzFllNWNupyU8nvWU4dxsSAnOzbryPjePx1tqrfKxOY3kFPA3wLMR8ddj1s2tmv0w+X0K9jZgkaTzJZ1E5Q+XjTWyfzJ7R9DlwOHstmAjY/NSd1+SzgPuA66NiOeqlk/08y8i57nZzxpJS6j8Xr7UyNiUObN87wKuoup8TXw86ynDuVlXSc7NunI/N1O8EyLlA3g3lf+U7fns61nZ8vcAm7Lp36DyroodwG7gi/XGF5TzSiqXqTuB7dljZbbuu8DT2bqNwNwcs62k8o64n40eG+B64PpsWsDN2fqngcsmGtvEn3W9nN+m8t9bjB67J+r9/AvKuTbLsYPKC8i/V8bjmc3/MXD3mHHJjieVK6sDwBCVv2VfV9Jzs17Ospyb9XLmem76o23MzCyJ4/GWmpmZlZALx8zMknDhmJlZEi4cMzNLwoVjZmZJuHDMSkLSZkkDkvL6R7xmpeLCMSuPrwPXFh3CrFlcOGaJSfrd7ENXT87+ZfluSRdFxEPAa0XnM2uWlvwsNbNWFhHbJG0E/htwCvC/IqKwjy8xS8WFY1aMr1D5PKr/B/xZwVnMkvAtNbNinAXMpvI/uZ5ccBazJFw4ZsVYD/wX4E7gqwVnMUvCt9TMEpP0SeBoRHxPUhvwiKQ/AP4r8K+A2ZL6gOsiYkuRWc3y5E+LNjOzJHxLzczMknDhmJlZEi4cMzNLwoVjZmZJuHDMzCwJF46ZmSXhwjEzsyT+P4CcxjbAUZoVAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots(1)\n",
        "mesh = discretize.TensorMesh([1,1])\n",
        "mesh.plot_grid(ax=ax,centers=True,faces=True)\n",
        "ax.set_xlim(-0.5,1.5)\n",
        "ax.set_ylim(-0.5,1.5);"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "bEL0LOMgs2fnWo7TPzdf",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "The figure above is of a single cell, which in 2D has four faces (the green triangles). Two of these are pointed in the x direction and two in the y direction.\n",
        "\n",
        "For this single cell, we need our projection matrices to pick out four cartesian vectors that we will then use in evaluating the full integral.\n",
        "\n",
        "Instead of creating the code from scratch in this notebook, we will leverage some of the internal functions that are in SimPEG. The code is straight forward, but can get a little bit intricate when working with multiple dimensions and indexing in a generalized way (i.e. works for all of the meshes).\n",
        "\n",
        "If you ever want to see the exact code that is executed in a function, you can use two question marks (one for the docs only). For example, execute **`mesh._getFacePxx??`** and Jupyter will show you more than you ever wanted. :)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "iooxa": {
          "id": {
            "block": "WnkAvpsgmPr1AwUiZyzi",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "mTYUldRSnhSmbCTyfMXG",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "# mesh._getFacePxx??"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "NusqKvJ3MlKANE8pOB3b",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "This function creates the P matrices that we need. The matrix just pick out the appropriate faces and are used later on to evaluate our full integral."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "QIfby0zQOGSFgWrebbFK",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "MMDATgoGoAnhFK9c4ZAA",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[1. 0. 0. 0.]\n",
            " [0. 0. 1. 0.]]\n"
          ]
        }
      ],
      "source": [
        "# face-x minus side and the face-y minus side\n",
        "P1 = mesh._getFacePxx()('fXm','fYm')\n",
        "print( P1.todense() )"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "fJFT2oA1PXeXict2NpfQ",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "We can create and plot all of the matrices that were in Figure 5:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "40uVtzW0G0t3yHRCcaaM",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "lrPqqVEm6QxfabDcYLxT",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAACICAYAAACBbYjWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAALjElEQVR4nO3dT6il510H8O+vmQyJHUspM2hoilNoKYQuWhmzsOCiGOg/1IVCC1YClbopRlwUlSwqRtypGzehxiLWFqRdFClIaCOtWDU3MTZJozQI1VAlE0vRASU0/lzkVCaTe+7cZO77vOe55/OBgbl3Luf98pzn+577O+85Z6q7AwAAwPJes3YAAACAfWEAAwAAGMQABgAAMIgBDAAAYBADGAAAwCAGMAAAgEEMYAAAAIMYwAapqrur6hNr59gH1np/ue/Hsdb7zf0/jrXeX+77cUav9U4PYFX1Y1X19aq6papeW1VPVtXb1851mJmyHqaqfquq7rnq69+uql9eM9M2M2WdwUx7d6ash5lp786UdRYz7d+Zsh5mpv07U9YZzLR3Z8p6mJn27q5lre5e69jHUlX3Jbklya1Jnunu31k50lZHZa2qu5Nc7O5PrJPuaFV1Mcnnu/tHq+o1Sb6Z5M7u/o9Vgx3iell3fa13kZ6NoWfo2hi6tt/0bAw9u4E8EwxgZ5M8nOR/kvx4d7+wcqStrs2a5PVJvrT55zckOZvk3zdff7i7Hx8e8ghV9WCSjyf5oSS/2N0/u3Kkra7NmuSXMtFa7xo9G0fP9puujaNr+0vPxtGzV+fMUjd8gt6Q5FySm5PcUlW/keT9SdLd71gz2CFeknUzVb8j2f1nMTY+meTuJD+c5IGq+qMk70zy7e5+35rBDvGSrBOu9a7Rs3H0bL/p2ji6tr/0bBw9exVmuAL2hSSfTfLmJLd198dWjrTVUVlnKNHmWZjH8+JJ4K0TPGN0aNYZ1nrX6Nk4erbfdG0cXdtfejaOnr06O30FrKp+Icn3uvtPq+qmJH9dVe/u7i+vne1aM2Xdprufr6qHknx3lwuUzJV11820d2fKus1Me3emrDOYaf/OlHWbmfbvTFl33Ux7d6as28y0d3cp685fAWOczZsSH03yc939zbXzHGWmrHC1mfbuTFnhWjPt35mywtVm2ru7lHWnP4aecarqjiRPJ/nS2pvyembKClebae/OlBWuNdP+nSkrXG2mvbtrWV0BAwAAGMQVMAAAgEEMYAAAAIMYwAAAAAaZagCrqo+uneHVmDH3jJmTeXPvmhnXccbMyZy5Z8y8i2Zdxxlzy7y/Zl3HGXPPmDlZJ/dUA1iSKe/YzJl7xszJvLl3zYzrOGPmZM7cM2beRbOu44y5Zd5fs67jjLlnzJyskHu2AQwAAGBai3wM/fnz5/vixYsnfruXL1/OhQsXTvx2lzZj7hkzJ8vmfuSRR57r7p1ZlKV6lsx5/8+YOZkzt56djBnv+2TO3DK/1K71LPG747VmzD1j5mSdrp1Z4mAXL17MwcHBEjcNq6mqb62d4Wp6xmmkZ7C8XetZomucTtu65iWIAAAAgxjAAAAABjGAAQAADGIAAwAAGMQABgAAMIgBDAAAYJBFPob+KJfuezDPXXl+9GFz/tzZHNx71/DjAsfn/MBs7FlgG+cHthl+BWyNjbjmcYHjc35gNvYssI3zA9t4CSIAAMAgBjAAAIBBDGAAAACDGMAAAAAGMYABAAAMYgADAAAYxAAGAAAwiAEMAABgEAMYAADAIMcawKrqPVX1T1X1dFX92tKhYB/pGYyha7A8PYPtrjuAVdVNSf4gyXuT3JHkQ1V1x9LBYJ/oGYyha7A8PYOjHecK2J1Jnu7uf+7u55N8NslPLxsL9o6ewRi6BsvTMzjCcQawNyb516u+fmbzPeDk6BmMoWuwPD2DIxxnAKtDvtcv+6Gqj1bVQVUdXL58+caTwX7RMxjjul3TM7hhHtPgCMcZwJ5J8qarvr49ybev/aHuvr+7L3X3pQsXLpxUPtgXegZjXLdregY3zGMaHOE4A9jDSd5aVW+uqrNJPpjkC8vGgr2jZzCGrsHy9AyOcOZ6P9Dd36uqjyX5iyQ3JXmgu59cPBnsET2DMXQNlqdncLTrDmBJ0t1fTPLFhbPAXtMzGEPXYHl6Btsd6z9iBgAA4MYZwAAAAAYxgAEAAAxiAAMAABjEAAYAADCIAQwAAGAQAxgAAMAgBjAAAIBBDGAAAACDGMAAAAAGGT6AnT93dvQhVz0ucHzOD8zGngW2cX5gmzOjD3hw712jDwlMwvmB2dizwDbOD2zjJYgAAACDGMAAAAAGMYABAAAMYgADAAAYxAAGAAAwiAEMAABgEAMYAADAIAYwAACAQQxgAAAAgxjAAAAABjGAAQAADGIAAwAAGMQABgAAMMiZtQOwnEv3PZjnrjw//Ljnz53Nwb13DT8urEXXgMM4N8AYs3XNFbBTbI2NuOZxYS26BhzGuQHGmK1rBjAAAIBBDGAAAACDGMAAAAAGMYABAAAMYgADAAAYxAAGAAAwiAEMAABgEAMYAADAIAYwAACAQa47gFXVA1X1bFU9MSIQ7Ctdg+XpGSxPz+Box7kC9qkk71k4B6BrMMKnomewtE9Fz2Cr6w5g3f2VJN8ZkAX2mq7B8vQMlqdncDTvAQMAABjkxAawqvpoVR1U1cHly5dP6maBq+gZLE/PYAxdY1+d2ADW3fd396XuvnThwoWTulngKnoGy9MzGEPX2FdegggAADDIcT6G/jNJvpbkbVX1TFV9ZPlYsH90DZanZ7A8PYOjnbneD3T3h0YEgX2na7A8PYPl6RkczUsQAQAABjGAAQAADGIAAwAAGMQABgAAMIgBDAAAYBADGAAAwCAGMAAAgEEMYAAAAIMYwAAAAAYxgAEAAAxiADvFzp87u1fHhbXoGnAY5wYYY7aunTnhHOyQg3vvWjsC7AVdAw7j3ABjzNY1V8AAAAAGMYABAAAMYgADAAAYxAAGAAAwiAEMAABgEAMYAADAINXdJ3+jVZeTfOvEbzg5n+S5BW53aTPmnjFzsmzuH+nuCwvd9iu2YM+SOe//GTMnc+bWs5Mx432fzJlb5pfaqZ4lfnc8xIy5Z8ycrNC1RQawpVTVQXdfWjvHKzVj7hkzJ/Pm3jUzruOMmZM5c8+YeRfNuo4z5pZ5f826jjPmnjFzsk5uL0EEAAAYxAAGAAAwyGwD2P1rB3iVbjh3Vb1QVY9V1RNV9WdV9QOb7z9QVc9W1RM3HvMl9natSTLnOp5I5sO6VlVvqqqHquqpqnqyqu45iWNt7O1aM+06LvKYVlW3VNXfVdU/bHr2mycRdmPGtZ4x8y6adR0X+91x8283VdXfV9Wf3+hxrrK3a/1KTfUesH1WVVe6+9zm759O8kh3/25V/USSK0n+uLvfvmpIOAUO61qSzyS5rbsfraof3HzvZ7r7GytGhWlt6dnvJXltd1+pqpuT/FWSe7r7b1aMCtPa9rvj5utfTXIpyeu6+wMrxtxLs10B40VfTfKWJOnuryT5zrpx4NT6apK3dPe/dfejSdLd/5XkqSRvXDUZnB7f71l395XN927e/PEsMZyM///dsapuT/L+JJ9cNdEeM4BNpqrOJHlvksfXzgKn2bauVdXFJO9M8rfjU8Hpcm3PNi+LeizJs0ke7G49gxt0yOPZ7yf5eJL/XS3UnjOAzePWzYPSQZJ/SfKHK+eB02pr16rqXJLPJfmV7v7PlfLBaXBoz7r7he5+R5Lbk9xZVV5aD6/ey3pWVR9I8mx3P7JutP12Zu0AHNt/bx6UgGUd2rXNe1I+l+TT3f358bHgVDnyMa27v1tVf5nkPUlO+kOmYF+8rGdV9a4kP1VV70tyS5LXVdWfdPfPr5JwT7kCBnAdVVV58Rn6p77/BmbgZFXVhap6/ebvtyb5yST/uG4qOF26+9e7+/buvpjkg0m+bPgazwA2uar6TJKvJXlbVT1TVR9ZOxOcQu9K8uEk7958pO9jm2cPgZNzW5KHqurrSR7Oi+8BO8mPyAbYCT6GHgAAYBBXwAAAAAYxgAEAAAxiAAMAABjEAAYAADCIAQwAAGAQAxgAAMAgBjAAAIBBDGAAAACD/B8zAXo72o9XPgAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1080x216 with 4 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "P1 = mesh._getFacePxx()('fXm','fYm')\n",
        "P2 = mesh._getFacePxx()('fXp','fYm')\n",
        "P3 = mesh._getFacePxx()('fXm','fYp')\n",
        "P4 = mesh._getFacePxx()('fXp','fYp')\n",
        "fig, ax = plt.subplots(1, 4, figsize=(15,3))\n",
        "\n",
        "def plot_projection(ii, ax, P):\n",
        "    ax.spy(P, ms=30)\n",
        "    if P.shape[1]==4:\n",
        "        ax.set_xticks(range(4))\n",
        "        ax.set_xticklabels(('x-','x+','y-','y+'))\n",
        "    ax.set_xlabel('P{}'.format(ii+1))\n",
        "\n",
        "list(map(plot_projection, range(4), ax, (P1, P2, P3, P4)));"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "zDZID1FHKMMBwSn6yAau",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "This assumes that we unwrapped our face vectors in certain way and they are ordered [xfaces, yfaces]."
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "7ctcHVZ9rbkU8JOISFLu",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "# All the cells!\n",
        "\n",
        "When we move to more cells the visualization of the matrix is not as clear, so instead we will show a connection from each cell center to the x and y faces that are selected by the P matrix."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": true,
        "iooxa": {
          "id": {
            "block": "WKYZQIAg0wA0So5pZDNO",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "jv4nylw4aMB6cLuAkltO",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "mesh = discretize.TensorMesh([3,4])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "EzFA7goTjbyc3oqloq7T",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "vGpoFc2QY9a8Gtz668VL",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAGoCAYAAACwmRWfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5BlZ3kf+O87I42YWKxGmGgGa2Qkeck6mJWFPTCITawGKmvAlSJJpSTZLntX5RTyYoyr9h+TZCExs5XFy+6SsDEgrZdQ2QTLqsXlxbEcKmHTEEoDQaqRFQERUQRCE80sNqQHCY3oUc+7f5xu9Q91T9++v847tz+fqqnqe/tozqPvvO9zTj9z751Saw0AAABAy/b0XQAAAADAdgwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDMaqlPKNUsrZUsrTpZT/r5Tyj0opl5dSbiml3FdKeaaUMt93nUC/LtAr/pdSyr8vpTxVSvl3pZRf7LtWoD8X6BX/cynliVLKd0spj5dS/nbftQL92qpfrPn+S0opf1JK+XyfdTIaAwwm4S/XWi9P8hNJXpPkf0jynSR/P8n7+ywMaMpmveJ7Sf5ykiuS/DdJ/kEp5fX9lQg0YLNe8X8m+dFa63+W5PVJfq6U8td6rBFow2b9YsVvJvlqL1UxNgYYTEyt9T8m+aMkr6q1/sta6z1Jnuy5LKAxG3rF36m1/rta6/la6xeT/OskN/VbIdCCDb3ikVrr99Z8+3yS/7yfyoDWrO0XSVJKuWn563/UZ12MzgCDiSmlXJPkrUlO9F0L0K6tekUpZX+6vz35ch91AW3Z2CtKKe8upTyd5GSSH0jyiR7LAxqytl+UUvYm+a0k70xSey2MkV3SdwHMpN8vpTyX5EySP0zy93quB2jTdr3io0n+OMmnp10Y0JRNe0Wt9f2llN9McmOSv7L8fWB326xfvCvJF2utD5RS/steq2NkBhhMwl+ptf7LvosAmrdlryilfCDdSz3fUGv1tyWwu23ZK5b7w4lSyk8n+Y0k//1UKwNas65flFJ+KN0A4yf7K4lxMsAAoCmllN9I8pYkN9dav9t3PcBF4ZIkP9J3EUBzXpvkZUm+UkpJkv1J9pdSTie5uta61Gdx7JzPwGAqSil7SykvSneDsaeU8qJSyqV91wW0pZTyN5P8XJK/VGv9dt/1AO0ppewppdxRSrmydF6b5FeSfKbv2oDm/FGSa9O91ezGJO9N9zk6NxpeXJwMMJiWX0hyNslHkvzF5a//j14rAlr095L8cJJ/v/zvuD9dSvlbfRcFNOevJvkPSZ5K8k+S/O/LvwCeV2v9fq319MqvdJ+NcW75ay5CxVuLAQAAgNZ5BQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5l3SdwE79dKXvrRee+212x73ve99Lz/wAz8w+YKGpL7RqG80g9b3wAMP/Gmt9c9OoaSJGKRfzMqfVV/UN5pZqW839Ipkdv68+qK+4bVcW6JXbDQrf159Ud9oZqW+LftFrfWi+vWTP/mTdRD/6l/9q4GO64v6RqO+0QxaX5L7awP7fthfg/SLWfmz6ov6RjMr9e2GXrGTPPqivtG0XF/LtdWqVwybR1/UNxr1jWbUfuEtJAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaZ4ABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaN7EBRinlY6WUb5VSHt7i+6WU8qFSyqOllIdKKT8xqVqAtukXwCB2e6849dSp3Pzxm3P66dN9l9I7WXTksDm9wrpYIYvOLOUwyVdgfDzJmy/w/bckecXyr7cn+cgEawHa9vHoF8D2Pp5d3CuOfe5YPv/Nz+fYZ4/1XUrvZNGRw5Y+Hr3CuogsVsxSDqXWOrnfvJRrk/yzWuurNvnenUnma62/s/z4kSRztdZTF/o9jxw5Uu+///4Lnvfad/9hkuTodS8Zqu5pWFhYyIEDByZ+ni9+/TtJdp7FtOob1k7rGzaHYbWc30oW33j/z2x7bCnlgVrrkUnXtHyuazPlfqFXrNIrVk2zX7Scn16xXt/9YnFpMV944nhqSfaUPTl69dHs27tv3THTWE+j7I9x1TdIFsNo+d5is9omlcNO6RXr6RWr+r632I29Inlhfa30imQ8/eKS8Zc1sKuTPLHm8cnl517QOEopb083Hc3BgwczPz8/0AkWFhZGLnJSlpaWplrfTs817fp2atj6pvX/1Hp+SQbeR42YaL9o+c9KrxjNKPVN4/+r9fwSvWKjvv68Tp49mZW/cqq15mvf+loO7z+87phprqdhzjOu+gbJYhgt31tsVtukchiWXrGeXrGqr3uL3dgrkhfW11qvSEbrF30OMMomz236cpBa611J7kq6yefc3NwFf+OjjxzPwsJCPv3rbxm1xomZn5/Pdv8f43DrnceTJL97x007+u+mVd+wdlrfsDkMq+X8br2z2x+t1reFifQLvWKVXrFqmv2i5fz0ivX67BennjqV6z/01/Jsefb5585csj+f+ZXHcujyQ88/N431NMr+GEd9g2YxjJbvLTbWNskcdkqvWE+vWNXnvcVu7RXJ+vpa6hXJePpFn/8Kyckk16x5fDjJkz3VArRNvwAGMZO94tjnjuV8Pb/uuaW6NBPvZd4pWXTkMDK9YsbJojOLOfQ5wPhUkl9c/hTg1yU5s937zoBdS78ABjGTveL4yeNZXFpc99zi0mLuO3lfTxX1RxYdOYxMr5hxsujMYg4TewtJKeV3kswleWkp5WSSv5Pk0iSptX40yb1J3prk0STPJLl9UrUAbdMvgEHs1l5x4o4TSfL8S24vss8aGCtZdORwYXrFXJLdvS5k0ZnFHCY2wKi1/uw2369JfmVS5wcuHvoFMAi9AhiEXgGzq8+3kAAAAAAMxAADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJp3Sd8FjNOr73x1bjp8UxaXbu23kIfuST7zvuTMyeSKw8mb3pvccMtUS1ibxb69+6Z67nV6zkIONM8e6TSwR2TRaeZa2ohm8njontz9uodz1WXnkg++avf2ikQWK3rOgUY1sC7skY4cVo3zWjpTA4wHTz+Yr/zJV/KDZ6/NlZdemVNPncrLXvyy6Rbx0D3JH7wrOXe2e3zmie5xMtWFsjaLQ5cfyqmnrt2VWciBpjWwNuyRVbLoNHEtbUgTeSyvi0MvOtc93q29IpHFigZyoEGNrAt7pCOHVeO8ls7UACNJFpcWc76ez7e//+0c/sDhHDx9MC//xstz2eJlUzn/3a97eHWBrDh3tvvbtClfUFayePK7T+44i4WFhRw4cGCk87eSxSg5DGttfq3kQIM+877VH1RXnDub0594R25714eTjGcvXtAb1u+R6z90fW6/8fa856feM72L7AA5DGtH+c14FgN7w/praS85NMa9RWfU66l7i+FtzK6VHGjMgNcQ9xZTureY8Rx2ZIz3FrP7GRglOb/3fE790Kl89ce+OrXTXnXZuc2/cebk1Gp4AVl05EBrtlgDW66ZSSvJs889mzsfuDO3ffK26Z23tRwSWazoK4dWuY50esohkcWK5nKgDa1dQxLX0xVyWDViFjP3CowVJSX7L9nfTXZufk8OXX5oOif+4Ku6l+ZsdMXh6Zx/E8NkMT8/n7m5udFO3FgW01wT6/JrLAcacsXhTdfGngPXZH5+PsmY9uIFlN8oq19v2CNTM0AOw9pJfrOexaCayKFR7i06w+bg3mJ4L8iusRxoxIDXEPcW8yP91oPmN+s57MQ4s5i5Aca+vfuyp+zJlZdemRO/9tj0bi5WvOm969+/nCSX7u+en7KVLA5dfignfnX3ZiEHmtXI2rBHVsmi0/u1tDG959HYuuhtfySyWNFIDjSmoXVhj3TksGpc19KZGmDceOjGvP7w6/Pk40fzzFPP9HNhXX7f4elPvCNXXXYuew5c08snva7NYt/efbs2CznQtAbWhj2yShadJq6lDWkij5U//x7/dZom9kciixUN5ECDGlkX9khHDqvGeS2dqQHGiTtOJEluvfN4nskz/RVywy3PfyjKNF+as9baLHrVcxZyoHn2SKeBPSKLTjPX0kY0k8cNt/T6w2kz+yORxYqec6BRDawLe6Qjh1XjvJbO7od4AgAAADPDAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADRvogOMUsqbSymPlFIeLaW8e5PvX1FK+YNSyh+XUr5cSrl9kvUAbdIrgEHoFcCg9AuYTRMbYJRS9ib5rSRvSfLKJD9bSnnlhsN+JclXaq0/nmQuyf9aStk3qZqA9ugVwCD0CmBQ+gXMrkm+AuO1SR6ttT5Wa11McneSt204piZ5cSmlJLk8yXeSPDfBmoD26BXAIPQKYFD6BcyoSyb4e1+d5Ik1j08mObrhmH+Y5FNJnkzy4iS31lrPb/yNSilvT/L2JDl48GDm5+cveOKFhbNZWlra9rhJWlhYSJIta3j66aenUt/CwtkL1rGVcda3XRbD2Gl9w+YwrM3qm0QOw2hhf2wwtl6R7KxftJLFhdaGXjGaYeqbZr/Yqr4W+kUr+2ON3npF0mQeLzCNfjHK/phWPxtWy/cWLWfX6N7Y1T+HbGc33VtMQsu9Imk7v3Hsj0kOMMomz9UNj386yYNJ3pjkR5L8i1LKv661fnfdf1TrXUnuSpIjR47Uubm5C574I48cz8LCQrY7bpIOHDiQJFvWMD8/P5X6PvLI8eU6btrRfzfO+rbLYhg7rW/YHIa1WX2TyGEYLeyPDcbWK5Kd9YtWsrjQ2tArRjNMfdPsF1vV10K/aGV/rNFbr0iazOMFptEvRtkf0+pnw2r53qLl7BrdG7v655Dt7KZ7i0louVckbec3jv0xybeQnExyzZrHh9NNONe6Pcnv1c6jSb6e5EcnWBPQHr0CGIReAQxKv4AZNckBxpeSvKKUct3yB+Lclu5lWmt9M8mbkqSUcjDJf5HksQnWBLRHrwAGoVcAg9IvYEZN7C0ktdbnSinvTPLpJHuTfKzW+uVSyi8vf/+jSY4l+Xgp5d+me6nXr9da/3RSNQHt0SuAQegVwKD0C5hdk/wMjNRa701y74bnPrrm6yeT/NeTrAFon14BDEKvAAalX8BsmuRbSAAAAADGwgADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngDEh39/3/Tz46gdz+unTfZfSO1l05MBWrI2OHFbJgs2ceupUbv74zdZFZLFCDmzGulgli84s5WCAMSGPX/t4zlxxJsc+e6zvUnoni44c2Iq10ZHDKlmwmWOfO5bPf/Pz1kVksUIObMa6WCWLzizlUGqtfdewI0eOHKn333//BY+59t1/mCQ5et1LplHSCywuLeb4E8eTkuwpe3L06qPZt3ffumMWFhZy4MCBidfyxa9/J8nOsxhXfYNkMYyd1jdsDsPaWN+kchjGShbfeP/PbHtsKeWBWuuRSdc0Kdv1i757RbL92tArRjNMfdPsF5vV10q/0CvW67tfLC4t5gtPHE/t+d5ilP0xzn6xXRbDaPneYqteMYkcdkqvWE+vWNXCvcVu6xXJ5j+HtNArkvH0C6/AmIDHzzz+/Ne11nWPdxtZdOTAVqyNjhxWyYLNPH7m8az8ldNuXxey6MiBzVgXq2TRmbUcLum7gEk4et1LsrCwkN+946apn/vUU6dy/YfemGdf9Ozzz505vz/3/vxjOXT5oeefm5+fz9zc5Ou79c7jSbLjLMZR36BZDGOn9Q2bw7DW1jfJHIZx653Hs7CwMPXztqjPXpEMtjb0iun2imS6/WJjfS31C71iPfcWnVH2x6z1C72io1esp1escm/R8XPIqnH0C6/AGLNjnzuW8/X8uueW6tJMvN9op2TRkQNbsTY6clglCzZjXaySRUcObMa6WCWLzizmYIAxZsdPHs/i0uK65xaXFnPfyft6qqg/sujIga1YGx05rJIFm7EuVsmiIwc2Y12skkVnFnOYybeQ9OnEHSf6LqEZsujIga1YGx05rJIFm1lZF3Nzc0m6lwfvVrLoyIHNWBerZNGZxRy8AgMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABo3iV9FzBOr77z1bnp8E1ZXLq130Ieuif5zPuSMyeTKw4nb3pvcsMtUy1hbRb79u6b6rnX6TkLOdA8e6TTwB6RRaeZa2kjmsnjoXty9+sezlWXnUs++Krd2ysSWazoOQca1cC6sEc6clg1zmvpTA0wHjz9YL7yJ1/JD569NldeemVOPXUqL3vxy6ZbxEP3JH/wruTc2e7xmSe6x8lUF8raLA5dfiinnrp2V2YhB5rWwNqwR1bJotPEtbQhTeSxvC4Ovehc93i39opEFisayIEGNbIu7JGOHFaN81o6UwOMJFlcWsz5ej7f/v63c/gDh3Pw9MG8/Bsvz2WLl03l/He/7uHVBbLi3Nnub9OmfEFZyeLJ7z654ywWFhZy4MCBkc7fShaj5DCstfm1kgMN+sz7Vn9QXXHubE5/4h257V0fTjKevXhBb1i/R67/0PW5/cbb856fes/0LrID5DCsHeU341kM7A3rr6W95NAY9xadUa+n7i2GtzG7VnKgMQNeQ9xbTOneYsZz2JEx3lvM7mdglOT83vM59UOn8tUf++rUTnvVZec2/8aZk1Or4QVk0ZEDrdliDWy5ZiatJM8+92zufODO3PbJ26Z33tZySGSxoq8cWuU60ukph0QWK5rLgTa0dg1JXE9XyGHViFnM3CswVpSU7L9kfzfZufk9OXT5oemc+IOv6l6as9EVh6dz/k0Mk8X8/Hzm5uZGO3FjWUxzTazLr7EcaMgVhzddG3sOXJP5+fkkY9qLF1B+o6x+vWGPTM0AOQxrJ/nNehaDaiKHRrm36Aybg3uL4b0gu8ZyoBEDXkPcW8yP9FsPmt+s57AT48xi5gYY+/buy56yJ1deemVO/Npj07u5WPGm965//3KSXLq/e37KVrI4dPmhnPjV3ZuFHGhWI2vDHlkli07v19LG9J5HY+uit/2RyGJFIznQmIbWhT3SkcOqcV1LZ2qAceOhG/P6w6/Pk48fzTNPPdPPhXX5fYenP/GOXHXZuew5cE0vn/S6Not9e/ft2izkQNMaWBv2yCpZdJq4ljakiTxW/vx7/NdpmtgfiSxWNJADDWpkXdgjHTmsGue1dKYGGCfuOJEkufXO43kmz/RXyA23PP+hKNN8ac5aa7PoVc9ZyIHm2SOdBvaILDrNXEsb0UweN9zS6w+nzeyPRBYres6BRjWwLuyRjhxWjfNaOtEP8SylvLmU8kgp5dFSyru3OGaulPJgKeXLpZTPTrIeoE16BTAIvQIYlH4Bs2lir8AopexN8ltJ/lKSk0m+VEr5VK31K2uOOZDkw0neXGv9ZinlqknVA7RJrwAGoVcAg9IvYHZN8hUYr03yaK31sVrrYpK7k7xtwzE/l+T3aq3fTJJa67cmWA/QJr0CGIReAQxKv4AZNcnPwLg6ydp/s+VkkqMbjvlzSS4tpcwneXGSf1Br/ccbf6NSytuTvD1JDh48uO17ghcWzmZpaanXzxpYWFhIsvX7l59++ump1LewcPaCdWxlnPVtl8UwdlrfsDkMa7P6JpHDMFrYHxuMrVckO+sXrWRxobWhV4xmmPqm2S+2qq+FftHK/lijt16RNJnHC0yjX4yyP6bVz4bV8r1Fy9k1ujd29c8h29lN9xaT0HKvSNrObxz7Y5IDjLLJc3WT8/9kkjcl2Z/keCnlC7XWr637j2q9K8ldSXLkyJG63b+7+5FHjmdhYWGi/77xdg4cOJAkW9Yw6X9/ecVHHjm+XMdNO/rvxlnfdlkMY6f1DZvDsDarbxI5DKOF/bHB2HpFsrN+0UoWF1obesVohqlvmv1iq/pa6Bet7I81eusVSZN5vMA0+sUo+2Na/WxYLd9btJxdo3tjV/8csp3ddG8xCS33iqTt/MaxPyY5wDiZ5Jo1jw8neXKTY/601vq9JN8rpXwuyY8necGNBjCz9ApgEHoFMCj9AmbUJD8D40tJXlFKua6Usi/JbUk+teGY/yfJXyylXFJK+TPpXtr11QnWBLRHrwAGoVcAg9IvYEZN7BUYtdbnSinvTPLpJHuTfKzW+uVSyi8vf/+jtdavllL+eZKHkpxP8tu11ocnVRPQHr0CGIReAQxKv4DZNcm3kKTWem+Sezc899ENjz+Q5AOTrANom14BDEKvAAalX8BsmuRbSAAAAADGwgADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwwAAACgeQYYAAAAQPMMMAAAAIDmGWAAAAAAzTPAAAAAAJp3wQFGKeU/K6X8yCbP3zC5kmbD9/d9Pw+++sGcfvp036X0ThadWc9BvxjerK+NQclh1SxnoVcM79RTp3Lzx2+eyXWxU7LozHIOesXwZnld7JQsOrOUw5YDjFLKLUn+XZJPllK+XEp5zZpvf3zShV3sHr/28Zy54kyOffZY36X0ThadWc5BvxjNLK+NnZDDqlnNQq8YzbHPHcvnv/n5mVsXw5BFZ1Zz0CtGM6vrYhiy6MxSDqXWuvk3SnkwyVtqradKKa9N8o+T/K1a6++VUk7UWl89zUJXHDlypN5///0XPObad/9hkuTodS+ZRkkvsLi0mONPHE9KsqfsydGrj2bf3n3rjllYWMiBAwcmXssXv/6dJDvPYlz1DZLFMHZa37A5DGtjfZPKYRgrWXzj/T+z7bGllAdqrUcGOO6i7Bd994pk+7WhV4xmmPqm2S82q6+VfqFXrNd3v1hcWswXnjie2vO9xSj7Y5z9YrsshtHyvcVWvWISOeyUXrGeXrGqhXuL3dYrks1/DmmhVyTj6RcXegvJ3lrrqSSptf6bJG9I8rdLKe9KsvnUgyTJ42cef/7rWuu6x7uNLDq7IAf9Yki7YG0MRA6rZjwLvWJIj595/PmAZnBd7IgsOjOeg14xpBlfFzsii86s5XDJBb73VCnlR2qt/yFJliegc0l+P8mPTaO4YR297iVZWFjI795x09TPfeqpU7n+Q2/Msy969vnnzpzfn3t//rEcuvzQ88/Nz89nbm7y9d165/Ek2XEW46hv0CyGsdP6hs1hWGvrm2QOw7j1zuNZWFgY9297UfaLPntFMtja0Cum2yuS6faLjfW11C/0ivXcW3RG2R+z1i/0io5esZ5escq9RcfPIavG0S8u9AqM/y7JnlLKK1eeqLU+leTNSf7GSGedYcc+dyzn6/l1zy3VpZl4v9FOyaKzS3LQL4awS9bGtuSwahdkoVcMYResi4HJorMLctArhrAL1sXAZNGZxRy2HGDUWv+41vrvk9xTSvn10tmf5H9L8o6pVXiROX7yeBaXFtc9t7i0mPtO3tdTRf2RRWc35KBfDGc3rI1ByGHVrGehVwxn1tfFTsiiM+s56BXDmfV1sROy6MxiDhd6C8mKo0l+M8l9SV6c5J8m+a8mWdTF7MQdJ/ouoRmy6OyyHPSLHdhla2NLcli1i7LQK3ZgZV3Mzc0l6V4evFvJorOLctArdmAXrYttyaIzizlc6C0kK84lOZtkf5IXJfl6rRtehwLQ0S+AQegVwCD0CmCdQQYYX0rXOF6T5C8k+dlSyv890aqAi5V+AQxCrwAGoVcA6wzyFpJfqrWu/IPHp5O8rZTyCxOsCbh46RfAIPQKYBB6BbDOtq/AWNM01j73f02mHOBipl8Ag9ArgEHoFcBGg7yFBAAAAKBXBhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmjfRAUYp5c2llEdKKY+WUt59geNeU0pZKqX89UnWA7RJrwAGoVcAg9IvYDZNbIBRStmb5LeSvCXJK5P8bCnllVsc95tJPj2pWoB26RXAIPQKYFD6BcyuSb4C47VJHq21PlZrXUxyd5K3bXLcryb5ZJJvTbAWoF16BTAIvQIYlH4BM+qSCf7eVyd5Ys3jk0mOrj2glHJ1kr+a5I1JXrPVb1RKeXuStyfJwYMHMz8/f8ETLyyczdLS0rbH9enpp5+eSn0LC2eTZMfnmlZ9w9ppfcPmMKyW82twf4ytVywfO3C/aBjvRnEAABtvSURBVDCLF9ArRjNMfdPsFy3n1+D+6K1XJG3ksbCwkGTrtTmN9TTK/hhnfdtlMYyW7y22qm0SOexUC3tjE7v655AWekVXR//3FrutVySb19dCr+jqGH1/THKAUTZ5rm54/PeT/HqtdamUzQ5f/o9qvSvJXUly5MiROjc3d8ETf+SR41lYWMh2x/Vpfn5+KvV95JHjSZK5uZt29N9Nq75h7bS+YXMYVsv5Nbg/xtYrkp31iwazeAG9YjTD1DfNftFyfg3uj956RdJGHgcOHEiSLWuYxnoaZX+Ms77tshhGy/cWW9U2iRx2qoW9sYld/XNIC70iaePeYrf1imTz+lroFcl49sckBxgnk1yz5vHhJE9uOOZIkruXm8ZLk7y1lPJcrfX3J1gX0Ba9AhiEXgEMSr+AGTXJAcaXkryilHJdkv+Y5LYkP7f2gFrrdStfl1I+nuSfaRqw6+gVwCD0CmBQ+gXMqIkNMGqtz5VS3pnuU333JvlYrfXLpZRfXv7+Ryd1buDioVcAg9ArgEHpFzC7JvkKjNRa701y74bnNm0Ytdb/dtTzvfrOV+emwzdlcenWUX+r0Tx0T/KZ9yVnTiZXHE7e9N7khlumWsLaLPbt3TfVc6/TcxZyuDhMu1c0xR7pNLBHZNFp5lq6iT56RTN5PHRP7n7dw7nqsnPJB1+1e3tFIosVPefQul17b9HAurBHOnJYNc5r6UQHGNP24OkH85U/+Up+8Oy1ufLSK3PqqVN52YtfNt0iHron+YN3Jee6T5vNmSe6x8lUF8raLA5dfiinnrp2V2YhB5rWwNqwR1bJotPEtbQhTeSxvC4Ovehc93i39opEFisayIEGNbIu7JGOHFaN81o6UwOMJFlcWsz5ej7f/v63c/gDh3Pw9MG8/Bsvz2WLl03l/He/7uHVBbLi3Nnub9OmfEFZyeLJ7z654ywWFhae/7TaYbWSxSg5DGttfq3kQIM+877VH1RXnDub0594R25714eTjGcvXtAb1u+R6z90fW6/8fa856feM72L7AA5DGtH+c14FgN7w/praS85NMa9RWfU66l7i+FtzK6VHGjMgNcQ9xZTureY8Rx2ZIz3FnsmVGL/SnJ+7/mc+qFT+eqPfXVqp73qsnObf+PMyanV8AKy6MiB1myxBrZcM5NWkmefezZ3PnBnbvvkbdM7b2s5JLJY0VcOrXId6fSUQyKLFc3lQBtau4Ykrqcr5LBqxCxm7hUYK0pK9l+yv5vs3PyeHLr80HRO/MFXdS/N2eiKw9M5/yaGyWIs//5yY1lMc02sy6+xHGjIFYc3XRt7DlyT+fn5JJP/t9rLb5TVrzfskakZIIdh7SS/Wc9iUE3k0Cj3Fp1hc3BvMbwXZNdYDjRiwGuIe4v5kX7rQfOb9Rx2YpxZzNwAY9/efdlT9uTKS6/MiV97bHo3Fyve9N71719Okkv3d89P2UoWhy4/lBO/unuzkAPNamRt2COrZNHp/VramN7zaGxd9LY/ElmsaCQHGtPQurBHOnJYNa5r6UwNMG48dGNef/j1efLxo3nmqWf6ubAuv+/w9CfekasuO5c9B67p5ZNe12axb+++XZuFHGhaA2vDHlkli04T19KGNJHHyp9/j/86TRP7I5HFigZyoEGNrAt7pCOHVeO8ls7UAOPEHSeSJLfeeTzP5Jn+Crnhluc/FGWaL81Za20Wveo5CznQPHuk08AekUWnmWtpI5rJ44Zbev3htJn9kchiRc850KgG1oU90pHDqnFeS2f3QzwBAACAmWGAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAGNCvr/v+3nw1Q/m9NOn+y6ld7LoyIGtWBsdOaySBZs59dSp3Pzxm62LyGKFHNiMdbFKFp1ZysEAY0Iev/bxnLniTI599ljfpfROFh05sBVroyOHVbJgM8c+dyyf/+bnrYvIYoUc2Ix1sUoWnVnKodRa+65hR44cOVLvv//+Cx5z7bv/MEly9LqXTKOkF1hcWszxJ44nJdlT9uTo1Uezb+++dccsLCzkwIEDE6/li1//TpKdZzGu+gbJYhg7rW/YHIa1sb5J5TCMlSy+8f6f2fbYUsoDtdYjk65pUrbrF333imT7taFXjGaY+qbZLzarr5V+oVes13e/WFxazBeeOJ7a873FKPtjnP1iuyyG0fK9xVa9YhI57JResZ5esaqFe4vd1iuSzX8OaaFXJOPpF16BMQGPn3n8+a9rrese7zay6MiBrVgbHTmskgWbefzM41n5K6fdvi5k0ZEDm7EuVsmiM2s5XNJ3AZNw9LqXZGFhIb97x01TP/epp07l+g+9Mc++6Nnnnztzfn/u/fnHcujyQ88/Nz8/n7m5ydd3653Hk2THWYyjvkGzGMZO6xs2h2GtrW+SOQzj1juPZ2FhYernbVGfvSIZbG3oFdPtFcl0+8XG+lrqF3rFeu4tOqPsj1nrF3pFR69YT69Y5d6i4+eQVePoF16BMWbHPncs5+v5dc8t1aWZeL/RTsmiIwe2Ym105LBKFmzGulgli44c2Ix1sUoWnVnMwQBjzI6fPJ7FpcV1zy0uLea+k/f1VFF/ZNGRA1uxNjpyWCULNmNdrJJFRw5sxrpYJYvOLOYwk28h6dOJO070XUIzZNGRA1uxNjpyWCULNrOyLubm5pJ0Lw/erWTRkQObsS5WyaIzizl4BQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaZ4ABAAAANM8AAwAAAGjeRAcYpZQ3l1IeKaU8Wkp59ybf//lSykPLv+4rpfz4JOsB2qRXAIPQK4BB6RcwmyY2wCil7E3yW0nekuSVSX62lPLKDYd9PcnNtdYbkhxLctek6gHapFcAg9ArgEHpFzC7JvkKjNcmebTW+litdTHJ3UnetvaAWut9tdb/tPzwC0kOT7AeoE16BTAIvQIYlH4BM+qSCf7eVyd5Ys3jk0mOXuD4X0ryR5t9o5Ty9iRvT5KDBw9mfn7+gideWDibpaWlbY/r09NPPz2V+hYWzibJjs81rfqGtdP6hs1hWC3n1+D+GFuvSHbWLxrM4gX0itEMU980+0XL+TW4P3rrFUkbeSwsLCTZem1OYz2Nsj/GWd92WQyj5XuLrWqbRA471cLe2MSu/jmkhV7R1dH/vcVu6xXJ5vW10Cu6OkbfH5McYJRNnqubHljKG9I1jr+w2fdrrXdl+WVdR44cqXNzcxc88UceOZ6FhYVsd1yf5ufnp1LfRx45niSZm7tpR//dtOob1k7rGzaHYbWcX4P7Y2y9ItlZv2gwixfQK0YzTH3T7Bct59fg/uitVyRt5HHgwIEk2bKGaaynUfbHOOvbLothtHxvsVVtk8hhp1rYG5vY1T+HtNArkjbuLXZbr0g2r6+FXpGMZ39McoBxMsk1ax4fTvLkxoNKKTck+e0kb6m1fnuC9QBt0iuAQegVwKD0C5hRk/wMjC8leUUp5bpSyr4ktyX51NoDSik/nOT3kvxCrfVrE6wFaJdeAQxCrwAGpV/AjJrYKzBqrc+VUt6Z5NNJ9ib5WK31y6WUX17+/keTvDfJDyb5cCklSZ6rtR6ZVE1Ae/QKYBB6BTAo/QJm1yTfQpJa671J7t3w3EfXfP03kvyNSdYAtE+vAAahVwCD0i9gNk3yLSQAAAAAY2GAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNu6TvAsbp1Xe+OjcdvimLS7f2W8hD9ySfeV9y5mRyxeHkTe9NbrhlqiWszWLf3n1TPfc6PWchB5pnj3Qa2COy6DRzLW1EM3k8dE/uft3Dueqyc8kHX7V7e0UiixU950CjGlgX9khHDqvGeS2dqQHGg6cfzFf+5Cv5wbPX5spLr8ypp07lZS9+2XSLeOie5A/elZw72z0+80T3OJnqQlmbxaHLD+XUU9fuyizkQNMaWBv2yCpZdJq4ljakiTyW18WhF53rHu/WXpHIYkUDOdCgRtaFPdKRw6pxXktnaoCRJItLizlfz+fb3/92Dn/gcA6ePpiXf+PluWzxsqmc/+7XPby6QFacO9v9bdqULygrWTz53Sd3nMXCwkIOHDgw0vlbyWKUHIa1Nr9WcqBBn3nf6g+qK86dzelPvCO3vevDScazFy/oDev3yPUfuj6333h73vNT75neRXaAHIa1o/xmPIuBvWH9tbSXHBrj3qIz6vXUvcXwNmbXSg40ZsBriHuLKd1bzHgOOzLGe4vZ/QyMkpzfez6nfuhUvvpjX53aaa+67Nzm3zhzcmo1vIAsOnKgNVusgS3XzKSV5Nnnns2dD9yZ2z552/TO21oOiSxW9JVDq1xHOj3lkMhiRXM50IbWriGJ6+kKOawaMYuZewXGipKS/Zfs7yY7N78nhy4/NJ0Tf/BV3UtzNrri8HTOv4lhspifn8/c3NxoJ24si2muiXX5NZYDDbni8KZrY8+BazI/P59kTHvxAspvlNWvN+yRqRkgh2HtJL9Zz2JQTeTQKPcWnWFzcG8xvBdk11gONGLAa4h7i/mRfutB85v1HHZinFnM3ABj39592VP25MpLr8yJX3tsejcXK9703vXvX06SS/d3z0/ZShaHLj+UE7+6e7OQA81qZG3YI6tk0en9WtqY3vNobF30tj8SWaxoJAca09C6sEc6clg1rmvpTA0wbjx0Y15/+PV58vGjeeapZ/q5sC6/7/D0J96Rqy47lz0Hrunlk17XZrFv775dm4UcaFoDa8MeWSWLThPX0oY0kcfKn3+P/zpNE/sjkcWKBnKgQY2sC3ukI4dV47yWztQA48QdJ5Ikt955PM/kmf4KueGW5z8UZZovzVlrbRa96jkLOdA8e6TTwB6RRaeZa2kjmsnjhlt6/eG0mf2RyGJFzznQqAbWhT3SkcOqcV5LZ/dDPAEAAICZYYABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaZ4ABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaZ4ABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaZ4ABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5hlgAAAAAM0zwAAAAACaN9EBRinlzaWUR0opj5ZS3r3J90sp5UPL33+olPITk6xnmr6/7/t58NUP5vTTp/supXey6Mhha7u5VyTWxgo5rJLF5nZ7rwAGp1/AbJrYAKOUsjfJbyV5S5JXJvnZUsorNxz2liSvWP719iQfmVQ90/b4tY/nzBVncuyzx/oupXey6Mhhc7u9VyTWxgo5rJLFC+kVwKD0C5hdpdY6md+4lJuS/N1a608vP/6bSVJr/Z/WHHNnkvla6+8sP34kyVyt9dRWv++RI0fq/ffff8FzX/vuP0ySHL3uJSP+XwxncWkxx584npRkT9mTo1cfzb69+9Yds7CwkAMHDky8li9+/TtJdp7FuOobJIth7LS+YXMY1sb6JpXDMFay+Mb7f2bbY0spD9Raj0yynkn1imT7ftF3r0i2Xxt6xWiGqW+a/WKz+lrpF3rFei30i+1Mo1+Msj+m1c+G1fK9RcvZtdYrls+za38OGcRuubeYlJZ7RdJ2fuPoF5eMv6znXZ3kiTWPTyY5OsAxVydZ1zhKKW9PNxnNwYMHMz8/P1ABCwsLOyp4XE6ePfn817XWfO1bX8vh/YfXHbO0tDTV+nZ6rnHVN0gWwxi2vmllvrG+SeUwikH30RSMrVckw/WLvnpFsv3a0CtGM0p908h9s/pa6xd6xXp99ovtTLNfDHOeafeznWr53qL17JKmekWyi38OGcRuubeYlJZ7RdJ+fslo/WKSA4yyyXMbX+4xyDGptd6V5K6km3zOzc1d8MTfmOtC2e64STj11Klc/6G/lmfLs88/d+aS/fnMrzyWQ5cfev65vuob1DjqGzSLvuqbpLX1TTKHcdTXgLH1imRn/aLPXpEMtjYa+7N6Ab1iNBvra61fNJZfb70i6b9fDEJ9o2m5vpZrS5qsb1f+HDIo9Y1GfaMZtb5JfojnySTXrHl8OMmTQxxzUTn2uWM5X8+ve26pLu3K9zHLoiOHbe3KXpFYGyvksEoWF7RrewWwY/oFzKhJDjC+lOQVpZTrSin7ktyW5FMbjvlUkl9c/hTg1yU5s937VFt3/OTxLC4trntucWkx9528r6eK+iOLjhy2tSt7RWJtrJDDKllc0K7tFcCO6Rcwoyb2FpJa63OllHcm+XSSvUk+Vmv9cinll5e//9Ek9yZ5a5JHkzyT5PZJ1TMtJ+440XcJzZBFRw4Xtlt7RWJtrJDDKllsbTf3CmBn9AuYXZP8DIzUWu9N1xzWPvfRNV/XJL8yyRqA9ukVwCD0CmBQ+gXMpkm+hQQAAABgLAwwAAAAgOYZYAAAAADNM8AAAAAAmmeAAQAAADTPAAMAAABongEGAAAA0DwDDAAAAKB5BhgAAABA8wwwAAAAgOYZYAAAAADNK7XWvmvYkVLKnyR5fIBDX5rkTydczijUNxr1jWbQ+l5ea/2zky5mUgbsF7PyZ9UX9Y1mVurbDb0imZ0/r76ob3gt15boFRvNyp9XX9Q3mlmpb9N+cdENMAZVSrm/1nqk7zq2or7RqG80rdc3Ta1nob7RqG80rdc3ba3nob7RtFxfy7Ul7dc3ba3nob7RqG80o9bnLSQAAABA8wwwAAAAgObN8gDjrr4L2Ib6RqO+0bRe3zS1noX6RqO+0bRe37S1nof6RtNyfS3XlrRf37S1nof6RqO+0YxU38x+BgYAAAAwO2b5FRgAAADAjDDAAAAAAJp30Q8wSilvLqU8Ukp5tJTy7k2+X0opH1r+/kOllJ9orL6fX67roVLKfaWUH2+pvjXHvaaUslRK+eut1VdKmSulPFhK+XIp5bOt1FZKuaKU8gellD9eru32adW2fP6PlVK+VUp5eIvv97o3pk2vmGx9a47TK4aor89+oVesp1dMtr41x+kVQ9SnV7RFv5hsfWuOm3q/0CtGrm9y/aLWetH+SrI3yX9Icn2SfUn+OMkrNxzz1iR/lKQkeV2SLzZW3+uTXLn89Vtaq2/Ncf9vknuT/PWW6ktyIMlXkvzw8uOrGqrtbyX5zeWv/2yS7yTZN8X8firJTyR5eIvv97Y3pv1Lr5h8fWuO0yuGq6+3fqFX7PjPSq8Yob41x+kVw9WnVzTyS7+YfH1rjptqv9ArxlLjxPrFxf4KjNcmebTW+litdTHJ3UnetuGYtyX5x7XzhSQHSikva6W+Wut9tdb/tPzwC0kOT6m2gepb9qtJPpnkW1OsLRmsvp9L8nu11m8mSa11WjUOUltN8uJSSklyebrG8dyU6kut9XPL59xKn3tj2vSKCde3TK8Yvr7e+oVesY5eMeH6lukVw9enV7RDv5hwfcv66Bd6xYgm2S8u9gHG1UmeWPP45PJzOz1mUnZ67l9KN4malm3rK6VcneSvJvnoFOtaMUh+fy7JlaWU+VLKA6WUX2yotn+Y5M8neTLJv03ya7XW89MpbyB97o1p0ytGo1eM5mLvF3rFzo+ZFL1iNHrFZO2mXpHoF6NquV/oFZM39N64ZCLlTE/Z5LmN/y7sIMdMysDnLqW8IV3j+AsTrWjDaTd5bmN9fz/Jr9dal7oB3lQNUt8lSX4yyZuS7E9yvJTyhVrr1xqo7aeTPJjkjUl+JMm/KKX861rrdydc26D63BvTpleMRq8YzcXeL/SKnR8zKXrFaPSKydpNvSLRL0bVcr/QKyZv6L1xsQ8wTia5Zs3jw+mmTDs9ZlIGOncp5YYkv53kLbXWb0+ptmSw+o4kuXu5abw0yVtLKc/VWn+/kfpOJvnTWuv3knyvlPK5JD+eZNLNY5Dabk/y/lprTfJoKeXrSX40yb+ZcG2D6nNvTJteMRq9YvL1tdwv9IqdHzMpesXk69MrhrebekWiX4yq5X6hV0ze8HujTumDPCbxK90A5rEk12X1A0x+bMMxP5P1HxDybxqr74eTPJrk9S3mt+H4j2e6H7Y1SH5/Pslnlo/9M0keTvKqRmr7SJK/u/z1wST/MclLp/xnfG22/vCc3vbGtH/pFZOvb8PxesXO6+u1X+gVO/qz0itGqG/D8XrFzuvTKxr5pV9Mvr4Nx0+tX+gVY6tzIv3ion4FRq31uVLKO5N8Ot2nsX6s1vrlUsovL3//o+k+sfat6TbnM+mmUS3V994kP5jkw8vTxedqrUcaqq83g9RXa/1qKeWfJ3koyfkkv11r3fSf65l2bUmOJfl4KeXfptucv15r/dNJ17ailPI7SeaSvLSUcjLJ30ly6Zr6etsb06ZXTKW+3rTcKwatLz32C71ilV4xlfp6o1eMRq9YT7+YSn290CtGN8l+UZYnIAAAAADNutj/FRIAAABgFzDAAAAAAJpngAEAAAA0zwADAAAAaJ4BBgAAANA8AwymopTyz0spC6WUf9Z3LUC79ApgUPoFMAi9YrYYYDAtH0jyC30XATRPrwAGpV8Ag9ArZogBBmNVSnlNKeWhUsqLSik/UEr5cinlVbXWzyR5qu/6gDboFcCg9AtgEHrF7nBJ3wUwW2qtXyqlfCrJ/5hkf5J/Umt9uOeygMboFcCg9AtgEHrF7mCAwSS8L8mXkjyb5F091wK0S68ABqVfAIPQK2act5AwCS9JcnmSFyd5Uc+1AO3SK4BB6RfAIPSKGWeAwSTcleQ9Sf5pkt/suRagXXoFMCj9AhiEXjHjvIWEsSql/GKS52qtnyil7E1yXynljUl+I8mPJrm8lHIyyS/VWj/dZ61Af/QKYFD6BTAIvWJ3KLXWvmsAAAAAuCBvIQEAAACaZ4ABAAAANM8AAwAAAGieAQYAAADQPAMMAAAAoHkGGAAAAEDzDDAAAACA5v3/MIA21XVXeWYAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 1080x432 with 4 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "P1 = mesh._getFacePxx()('fXm','fYm')\n",
        "P2 = mesh._getFacePxx()('fXp','fYm')\n",
        "P3 = mesh._getFacePxx()('fXm','fYp')\n",
        "P4 = mesh._getFacePxx()('fXp','fYp')\n",
        "fig, ax = plt.subplots(1,4, figsize=(15,6))\n",
        "def plot_projection(ii, ax, P):\n",
        "    x = P * np.r_[mesh.gridFx[:,0], mesh.gridFy[:,0]]\n",
        "    y = P * np.r_[mesh.gridFx[:,1], mesh.gridFy[:,1]]\n",
        "    xx, xy = x[:mesh.nC], x[mesh.nC:]\n",
        "    yx, yy = y[:mesh.nC], y[mesh.nC:]\n",
        "    ax.plot(np.c_[xx, mesh.gridCC[:,0], xx*np.nan].flatten(), np.c_[yx, mesh.gridCC[:,1], yx*np.nan].flatten(), 'k-')\n",
        "    ax.plot(np.c_[xy, mesh.gridCC[:,0], xy*np.nan].flatten(), np.c_[yy, mesh.gridCC[:,1], yy*np.nan].flatten(), 'k-')\n",
        "    ax.plot(xx, yx, 'g>')\n",
        "    ax.plot(xy, yy, 'g^')\n",
        "    mesh.plot_grid(ax=ax, centers=True)\n",
        "    ax.set_title('P{}'.format(ii+1))\n",
        "list(map(plot_projection, range(4), ax.flatten(), (P1, P2, P3, P4)));\n",
        "plt.tight_layout()"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "KPtsBYb1war7X3xTvF4V",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "Each cell (the center is the red dot) has two faces that it uses, one x component and one y component."
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "82G9e2DAu4MWgeggqw79",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "The four $P$ matrices are assembled and summed in the **getFaceInnerProduct** code using this formula:\n",
        "\n",
        "$$\n",
        "\\mathbf{M}_f(\\sigma) = \\frac{1}{4}\\sum_{i=1}^{4} \\mathbf{P}_i^\\top \\sqrt{v} \\boldsymbol{\\Sigma}^{-1} \\sqrt{v} \\mathbf{P}_i\n",
        "$$\n",
        "\n",
        "This function takes a sigma variable which can have any sort of anisotropy:\n",
        "\n",
        "$$\n",
        "\\begin{align}\\begin{aligned}\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & 0 \\\\ 0 & \\mu_{1} \\end{matrix}\\right]\\end{split}\\\\\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & 0 \\\\ 0 & \\mu_{2} \\end{matrix}\\right]\\end{split}\\\\\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & \\mu_{3} \\\\ \\mu_{3} & \\mu_{2} \\end{matrix}\\right]\\end{split}\\end{aligned}\\end{align}\n",
        "$$ (anisotropy)\n",
        "\n",
        "The matrices that are produced by this function can be seen below!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "oG3qra6QUeqwyfIQunfa",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "GzIgLkR3zaDoWst5C0gy",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAADyCAYAAAA89ja4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfbxdZX3n/c83AQaJUMGEB/MwUpuWYWhEm0EtFKGKDdaZqB1vAevTC0R6lyrcdSp1ZgTGscUZK2UUzURMQatwAxLNtEHIWC3DIJgEIxDCQwSUkwRCACu0DiQ5v/ljrRM3J3vtvfY5a6+9rr2/79drv3L22mut69oh58t1rYffUkRgZmZmZmZmzTFj0B0wMzMzMzOzF/JEzczMzMzMrGE8UTMzMzMzM2sYT9TMzMzMzMwaxhM1MzMzMzOzhvFEzczMzMzMrGEGMlGTtETS/ZI2S7pgAO0/IuluSRskretzWyskbZd0T8uyQyStkfRg/ufBNbd/kaQt+fffIOnNfWx/vqTvSNokaaOkD+fL+/530KHtWr6/pP0lfV/SD/P2L86X1/bf33ozStmUtzey+TTIbOrSft+/v7MpPc4mZ9MoZFPejvOphep+jpqkmcADwCnAGLAWOD0i7q2xD48AiyNiRw1tnQg8C3w5Io7Jl/0X4KmIuCQP3IMj4qM1tn8R8GxEfLofbU5q/wjgiIi4U9KBwHrgrcD76PPfQYe2/x9q+P6SBMyKiGcl7QvcCnwYeDs1/fe38kYtm/L2RjafBplNXdrvez45m9LibHI2MSLZlLfvfGoxiDNqxwGbI+KhiHgeuAZYOoB+1CIibgGemrR4KXBV/vNVZL8AdbZfm4jYFhF35j8/A2wC5lLD30GHtmsRmWfzt/vmr6DG//7Wk5HKJhjtfBpkNnVpv++cTclxNmWcTUOeTXmbzqcWg5iozQUebXk/Ro3/AHIB3CxpvaSza24b4LCI2AbZLwRw6AD6cK6ku/LT+7WcPpb0cuBVwB3U/HcwqW2o6ftLmilpA7AdWBMRtX93K83ZlGnCv89a82mQ2dSmfajh+zubkuJsyjTh36ezyWOnWg1ioqY2y+q9/hKOj4hXA6cCf5if4h4lXwBeARwLbAP+ot8NSnox8HXgvIj4Wb/b69J2bd8/InZHxLHAPOA4Scf0qy2bNmdTM9SaT4PMpoL2a/n+zqakOJuawdnksVPtBjFRGwPmt7yfB2ytswMRsTX/czuwkuyygjo9nl8DPHEt8PY6G4+Ix/NfgnHgi/T5++fXGH8d+GpE3JAvruXvoF3bdX//vM2fAt8FljDg//5WyNmUGZl8GmQ2FbVfdz45m5LgbMo4m0Yom/I2Rz6fBjFRWwsslHSkpP2A04BVdTUuaVZ+cySSZgFvAu7pvFXlVgHvzX9+L/DNOhuf+Ieeext9/P75TaFfAjZFxGdaPur730FR23V9f0lzJL0k//lFwBuB+xjwf38r5GzKjEQ+DTKbOrVfx/d3NiXH2ZRxNg15NuXtOJ9a1F71EUBZSc+/BGYCKyLikzW2/ctkR4MA9gG+1s/2JV0NnATMBh4HLgS+AVwLLAB+ArwjIvpy02pB+yeRnboO4BHggxPX/fah/ROA/wXcDYzniz9Gdr1zX/8OOrR9OjV8f0mLyG54nUl2UOTaiPhPkl7KFL/775w8K558anep9tff9dxNEbFkSp0fUaOUTXmbI5tPg8ymLu33PZ/6kU1QPp+cTb1zNjmbGIFsytv32KnFQCZqZjY1v/HK/eOOm+aVWnffI360PiIW97lLZmZA+XxyNplZnVIeO+0z6A6YWS+C3THefTUzs9o5n8ysidLNJk/UzBISwHjtxb7MzLpzPplZE6WcTZ6omSUkCHZGueuszczq5HwysyZKOZs8UTNLTKpHhcxs+DmfzKyJUs0mT9TMEhLA7kTDxsyGm/PJzJoo5WwaxHPUAJB09qDaHvX2R/m7D0P740Spl02dfz/d/qi1XVX7zqb+GoZ/Iym27fbTbz/VsdPAJmrAQP+Dj3j7o/zdk24/gN0RpV42Lf79dPuj1va02y+bTzYtSf8bSbhtt59w+ymPnaY1UZO0RNL9kjZLuqCqTplZe0Gws+Rr1DmfzOpVNp9GnbPJrF4pj52m/MBrSTOBB4BTgDFgLXB6RNxbtM1++mexP7MA2Mlz7Ms/A+BXf+OXp9SH6XjiiSeYM2dO7e02of1R/u5Na3/9+vU7IqJ0ZxYt2jdWrZ5dat0j5z/WqIc21qnXfJp5wKzY9yWH7Hm/+5/+kZkHzOKYlx1WR3dfoEn/Pt3+6LQ9uf1eswnK55OzqbexU2s+TWQTMHL51KTfD4AH7/pJ2/UWLlrQ877L7Ktp33+Q7Y/S2Gk6xUSOAzZHxEMAkq4BlgKFYbM/s3iN3rDX8jXrrptGN8zSJenHvayfPQvESugpn/Z9ySG8/AP/317L1118fj/7aNZYvWYTOJ9K6nns5HxqplPn/lHb5Teu++xA9zUKRmnsNJ2J2lzg0Zb3Y8BrptcdM+tM7EaD7kQKnE9mtXM+leBsMqtdutk0nXvU2n3jva6jlHS2pHWS1u3kuWk0Z2YBjEe5Vxnd7pWQtFTSXZI25L/HJ+TL50v6jqRNkjZK+nClX3T6uuZTazbt/qd/rKlbZsOrbD6NuJ7HTs4ns+mpeuxUp+mcURsD5re8nwdsnbxSRCwHlgMcpEMa+Fdglpaqjgrl90pcTsu9EpJWTbpX4tvAqogISYuAa4GjgF3AH0fEnZIOBNZLWtPpPouadc2n1mza/2XznU1mFUj1qHWNeh47OZ/Mpi/VbJrORG0tsFDSkcAW4DTgjEp6ZWZtBbAzKnuqRtd7JSLi2Zb1Z+VdICK2Advyn5+RtInskp6mTNScT2Y1qzifhpWzyaxmKWfTlCdqEbFL0rnATcBMYEVEbOy0za/+xi+3LRwy/tjCwm1mHP7gVLtoNnSCSo8KlbpXQtLbgD8HDgV+t83nLwdeBdxRVcemq9d8OuZlh7W9Md/ZZFZexfk0lKYydnI+paWwMMiWCouMTGFfoyzlbJrOGTUiYjWwuqK+mFkXgdhd/tbS2ZLWtbxfnl9OM6HUvRIRsRJYKelE4BPAG/fsQHox8HXgvIj4WdmO1cH5ZFavHvNpZDmbzOqVcjZNa6JmZvUbj9JHhXZ0eRZIqXslJkTELZJeIWl2ROyQtC/ZJO2rEXFD2U6Z2fDqIZ/MzGqTajZ5omaWkIpP33e9V0LSrwA/youJvBrYD3hSkoAvAZsi4jNVdcjM0pXy5UVmNrxSziZP1MwSEoidUc2vbdG9EpLOyT9fBvwe8B5JO4GfA+/MJ20nAO8G7pa0Id/lx/JLesxsBFWZT2ZmVUk5m9K8YNNshO3OH9zY7VVGRKyOiF+NiFdExCfzZcvySRoR8amI+JcRcWxEvC4ibs2X3xoRiohF+WfHepJmZlVl01Sf8VhmWzMbPVWOnerMpzSnl2YjKkLsTrTErJkNt6ryaTrPeCy5rZmNkCrHTnXnUyMmap3KyBaVn3XpWRtV44leZ50iZ5NZbyrKpyk/47HMtsPC+TRYRSXyi0rqFy2va1+jrsKxU6351IiJmpmVk90Q6zNqZtY8PeRTt0eHTOcZj6W2NbPR0ePYqVH55ImaWVJ86aOZNVXpfOr26JDpPOOx1LZmNkp6Gjs1Kp88UTNLSAA7Y+agu2FmtpcK82nKz3jsdVszG34Vj51qzSdP1MwSEsiXPppZI1WYT1N+xiPw027bmtloqXjsVGs+eaJmlphxX/poZg1VRT5N5xmPQNttp90pM0taVWOnuvOp8RO1ogpFp8x4R9vla8av62d3zAbKxUSaoyibjrrw0rbL77v4/H52x2zgqsyn/LmMqyctW9by86eAT5XddtQ4nwan1wqOde1rlFU9dqoznxo/UTOzXwjE7nB5fjNrHueTmTVRytnkiZpZYsZ9Rs3MGsr5ZGZNlGo2eaJmlpAIueqjmTWS88nMmijlbPJEzSwhAX6Ompk1kvPJzJoo5WzyRM0sMS4mYmZN5XwysyZKNZuSnagVVXcsqgbZaRuzVARiPNEbYkdFUfW0omprnbYxS4nzqfmcT4NTVMERiqs4Fm1T5b5GQcrZlOxEzWxUpXpUyMyGn/PJzJoo1WzyRM0sIUG6N8Sa2XBzPplZE6WcTZ6omSUkgPFEb4g1s+HmfDKzJko5m9LstdkI241KvcqQtETS/ZI2S7qgzedLJd0laYOkdZJOKLutmY2eqrLJzKxKVY6d6uQzamYJiVBlR4UkzQQuB04BxoC1klZFxL0tq30bWBURIWkRcC1wVMltzWyEVJlPZmZVSTmbpjVRk/QI8AywG9gVEYur6JSZFavwWSDHAZsj4iEASdcAS4E9k62IeLZl/VlkVxCU2nbQnE9m9Uv1WUV1cjaZ1S/VbKrijNrJEbGjgv1UolMJ/vHHFrZdPuPwB/vVHbNKBTBe3an5ucCjLe/HgNdMXknS24A/Bw4FfreXbRugMfnUqcS1s8mGQcX5NOwak03gfGqqKkvtF+1rqvtLScrZ5EsfzRISiJ3jpSsXzZa0ruX98ohY3vK+XWrFXgsiVgIrJZ0IfAJ4Y9ltzWx09JhPZma1SDmbpjtRC+BmSQH890mDQDPrgx6eBbKjyyU1Y8D8lvfzgK1FK0fELZJeIWl2r9sOiPPJrGapPquoZs4ms5qlmk3TnagdHxFbJR0KrJF0X0Tc0rqCpLOBswEWLFgwzebMRlsgxqOy0/drgYWSjgS2AKcBZ7SuIOlXgB/lxUReDewHPAn8tNu2DdAxn5xNZtWqOJ+GmcdOZjVKOZumNb2MiK35n9uBlWQFBiavszwiFkfE4jlz5kynOTMDxplR6tVNROwCzgVuAjYB10bERknnSDonX+33gHskbSCr8vjOyLTdtg9fd8q65ZOzyax6VWTTsPPYyax+VY2d6jblM2qSZgEzIuKZ/Oc3Af+psp6Z2V4iYHeFR4UiYjWwetKyZS0/fwr4VNltm8L5ZFa/qvNpGDmbzOqXcjZN59LHw8gKDEzs52sR8a1KetUnRRWKXNHIUpLq6fuaJZVPziYbFs6nrpLKJnA+1aGo6mJRpcapVHDstE2V1SWbKtVsmvJELX9+0isr7IuZdRGInZFm5aI6OZ/M6ud86s7ZZFa/lLPJ5fnNEhKke1TIzIab88nMmijlbPJEzSwpYjyad7OrmZnzycyaKd1s8kTNLDHjbZ81bWY2eM4nM2uiVLPJEzWzhKRcucjMhpvzycyaKOVs8kSN4gpFp8x4R+E2a8av61d3zAoFYtd4mjfEWu+KsumoCy8t3Oa+i8/vV3fMOnI+jRbnU/9NpYJjr/ua6v5SknI2eaJmlphUT9+b2fBzPplZE6WaTZ6omSUk5cpFZjbcnE9m1kQpZ5MnamaJSbVykZkNP+eTmTVRqtmUZq/NRlWI8ZIvM7NaVZhNkpZIul/SZkkXtPn8XZLuyl+3SXply2ePSLpb0gZJ6yr8hmaWoorHTnXmk8+omSUkSPc6azMbblXlk6SZwOXAKcAYsFbSqoi4t2W1h4HXR8TTkk4FlgOvafn85IjYMe3OmFnyqhw71Z1PnqiZJSSAXeM+EW5mzVNhPh0HbI6IhwAkXQMsBfYMhCLitpb1bwfmVdGwmQ2fisdOteaTJ2oddCrBX1S632X7rd98WaN1KnFdVBrbZbGtDhXl01zg0Zb3Y7zwaPRkZwI3trwP4GZJAfz3iFheRaesHOdT/02l1H6nbXp9DECnfTVVhWOnWvPJEzWzhAS+/8zMmqmHfJo96d6M5ZMGK+12Eu12JOlksoHQCS2Lj4+IrZIOBdZIui8ibinTMTMbPj2OnRqVT56omSXG96iZWVOVzKcdEbG4w+djwPyW9/OArZNXkrQIuAI4NSKenFgeEVvzP7dLWkl2qZInamYjrIexU6PyyTe7mKUkcNVHM2umkvlUwlpgoaQjJe0HnAasal1B0gLgBuDdEfFAy/JZkg6c+Bl4E3BPRd/QzFJU7dip1nzyGTWzhKT80EYzG25V5VNE7JJ0LnATMBNYEREbJZ2Tf74M+DjwUuDzkgB25UfBDwNW5sv2Ab4WEd+adqfMLFlVjp3qzidP1MwSEqjSqo+SlgCXkYXNFRFxyaTP3wV8NH/7LPAHEfHD/LPzgbPIMvBu4P0R8X8q65yZJaXKfIqI1cDqScuWtfx8Fln+TN7uIeCVk5eb2eiqeuxUZz55ojZFRdUdxx9bWLjNjMMf7Fd3bIRERUeFpvMsEElzgQ8BR0fEzyVdS3b6/8pKOmdTVlQ9zdlkdagqn2w4OZ8Gp6iCI/RexbHKfdUl1WzyRM0sMRUWE5nus0D2AV4kaSdwAG1upjWz0eJiR2bWRKlmkydqZgmJ6Ok6624lZqf8LJCI2CLp08BPgJ8DN0fEzWU7ZmbDp8d8MjOrRcrZ5ImaWWJ6OH3frcTslJ8FIulgsrNvRwI/Ba6T9PsR8ddlO2dmwyfVy4vMbLilmk2eqJklReyu7obY6TwL5I3AwxHxRL7ODcBvAp6omY2sSvPJzKwi6WaTJ2pmCam4PP+eZ4EAW8iKgZzRukLRs0DILnl8raQDyC59fAPQepmlmY0YPz7EzJoo5WzqOr2UtELSdkn3tCw7RNIaSQ/mfx7c326aGQCRXWtd5tV1VxG7gIlngWwCrp14FsjE80B44bNANkzc8xYRdwDXA3eSleafQVYRslbOJ7MGqSibhoGzyaxBKhw71a3MGbUrgc8BX25ZdgHw7Yi4RNIF+fuPttl25HQqI1tUftalZ60XVVYumuqzQPLPLgQurKwzU3MlzqdSnE1Wh1Qrq/XBlTibSnM+VaeoPH6nkvpFn9Wxr7qkmk1dz6hFxC3AU5MWLwWuyn++Cnhrxf0yszaC7IbYMq9R4Hwya46y+TQKnE1mzZHy2Gmq96gdFhHbACJim6RDK+yTmRVSstdZ18j5ZDYQzqcunE1mA5FuNvW9mIiks4GzARYsWNDv5syG3vh4mmHTNM4ms+o5n6rhfDKrVqrZNNValY9LOgIg/3N70YoRsTwiFkfE4jlz5kyxOTODiZtd0zx9X6NS+eRsMqtW2XwaYR47mQ1AymOnqU7UVgHvzX9+L/DNarpjZt2Mh0q9RpjzyWxAnE0dOZvMBiTVsVPXSx8lXQ2cBMyWNEZW5e0S4FpJZ5I9T+kd/ezksCiqUHTKjPZ/fWvGr+tndyxRTSwfOyjOp2oUZdNRF17advl9F5/fz+5YwpxPGWdTdZxP1ehUdbFTFcd+76suqWZT14laRJxe8NEbKu6LmZXQxFPzg+J8MmsW51PG2WTWLKlmU9+LiZhZdSLhykVmNtycT2bWRClnkydqZimJdI8KmdmQcz6ZWRMlnE2eqJmlJtHrrM1sBDifzKyJEs0mT9TMEpPqUSEzG37OJzNrolSzyRO1Biiq7lhUDbLTNjb8Uq1cZOkpqp5WVG2t0zY2GpxPVhfnU3WKqjgWVXDsVPWxyn1VKdVs8kTNLCFBukeFzGy4OZ/MrIlSziZP1MxSEhDjaYaNmQ0555OZNVHC2eSJmllqEj19b2YjwPlkZk2UaDZ5omaWFCV7+t7Mhp3zycyaKN1s8kTNLDWJHhUysxHgfDKzJko0m2YMugNm1oP8oY1lXmVIWiLpfkmbJV3Q5vN3Sborf90m6ZUtn71E0vWS7pO0SdLrKvymZpaakvlkZlarisdOdfIZtQbrVIJ//LGFbZfPOPzBfnXHmqKio0KSZgKXA6cAY8BaSasi4t6W1R4GXh8RT0s6FVgOvCb/7DLgWxHxbyXtBxxQTc+s6TqVuHY2jbhEj1rb8HA+9V9RqX3ovdx+lfvqKNFs8kTNLDXVHfE5DtgcEQ8BSLoGWArsmahFxG0t698OzMvXPQg4EXhfvt7zwPNVdczMEtXAI9JmZqlmky99NEtNlHzBbEnrWl5nT9rTXODRlvdj+bIiZwI35j//MvAE8FeSfiDpCkmzpvW9zCx95bKpq2lelt1xWzMbQeXHTl3VmU8+o2aWkqCXo0I7ImJxh8/b7ahtTEk6mWyidkK+aB/g1cAfRcQdki4DLgD+Y9nOmdmQ6S2fCk3nsuyS25rZKKkom6D+fPIZNbPERJR7lTAGzG95Pw/YOnklSYuAK4ClEfFky7ZjEXFH/v56sombmY2wirJpz2XZ+WXVE5dlt7QTt0XE0/nbPZdll9nWzEZPhWOnWvPJEzWz1FR3+n4tsFDSkXkxkNOAVa0rSFoA3AC8OyIe2NOFiMeARyX9Wr7oDbTc22ZmI6qabJrOZdm9bmtmo6C6sVOt+eRLHxNVVKHIFY2Gn8arOX0fEbsknQvcBMwEVkTERknn5J8vAz4OvBT4vCSAXS2XU/4R8NV8kvcQ8P5KOmZJczaNtpL5NFvSupb3yyNieetu2mxT9rLs0tva6HE+9aao6mKnSo1Fn1W5r6noYezUqHzyRM0sJT3c7FpqdxGrgdWTli1r+fks4KyCbTcAne6BM7NRUj6fut0/2+tl2adOuiy767ZmNkJ6Gzs1Kp986aNZUpTdEFvmZWZWq8qyacqXZZfZ1sxGTaVjp1rzyWfUzFLji3jMrKkqyKfpXJZdtO30e2VmSato7FR3PnmiZpYaT9TMrKmqGwxN57LsvbY1sxHXnNtGesonT9TMUuOJmpk1lfPJzJoo0WzqOlGTtAJ4C7A9Io7Jl10EfAB4Il/tY/kM0QasqELRKTPeUbjNmvHr+tUdq1pUV/VxGDif0lGUTUddeGnhNvddfH6/umP94Hzaw9mUFudTb4oqOELvlRqr3FehhLOpTDGRK4ElbZZfGhHH5i8HjVldqnsWyDC4EueTWXM4myZcibPJrDkSHTt1nahFxC3AUzX0xcysJ84nM2siZ5OZVWE65fnPlXSXpBWSDi5aSdLZktZJWvfEE08UrWZmJSnKvUZc13xyNplVz9nUlcdOZgOQ6thpqhO1LwCvAI4FtgF/UbRiRCzPS1IunjNnzhSbM7M9/By1bkrlk7PJrA+cTZ147GQ2KImOnaY0UYuIxyNid0SMA18Ejqu2W2bWVtlrrBt4VKguziezAXE2deRsMhuQhMdOU5qoSTqi5e3bgHuq6Y6ZdaPxcq9R5XwyGxxnUzFnk9ngpDp2KlOe/2rgJGC2pDHgQuAkSceSzT0fAT7Yxz5aBTqV4C8q3e+y/Q3VwCM+g+J8Sl+nEtfOpgQ5n4Bqs+nBu35SaclzK8/51Lsq/+0V7Uv6XO87SzSbuk7UIuL0Nou/1Ie+mFkZiYZNPzifzBrG+QQ4m8waJ9Fs6jpRM7PmaGpVIjMz55OZNVHK2eSJmllqGliVyMwMcD6ZWTMlmk2eqJklpok3u5qZgfPJzJop1WzyRM0sNYmevjezEeB8MrMmSjSbPFGzwgpFrmjUQBVfZy1pCXAZMBO4IiIumfT5u4CP5m+fBf4gIn7Y8vlMYB2wJSLeUl3PzJxNyUn4PpAmW7hoATeu27v6XadKkEWfuRpkdZxPCUk4mzxRM0tNRWGTT7IuB04BxoC1klZFxL0tqz0MvD4inpZ0KrAceE3L5x8GNgEHVdMrM0taooMhMxtyiWbTlB54bWYDFCVf3R0HbI6IhyLieeAaYOkLmoq4LSKezt/eDsyb+EzSPOB3gSum8W3MbJhUk01mZtWqbuxUK59RM0tMD6fvZ0ta1/J+eUQsb3k/F3i05f0YLzxbNtmZwI0t7/8S+BPgwNI9MrOhlurlRWY23FLNJk/UzFJTPmx2RMTiDp+3q1Xbdu+STiabqJ2Qv38LsD0i1ks6qXSPzGy4JToYMrMhl2g2eaJmlpJqb4gdA+a3vJ8HbJ28kqRFZJc3nhoRT+aLjwf+jaQ3A/sDB0n664j4/cp6Z2ZpSfiGfTMbYglnkydqVqjXikadtrEKVRc2a4GFko4EtgCnAWe0riBpAXAD8O6IeGBPFyL+FPjTfJ2TgI94kmZ1cTY1WKKDoRR1quBYVPXR1SD7z/nUUIlmkydqZqmpKGwiYpekc4GbyMrzr4iIjZLOyT9fBnwceCnweUkAu7pcTmlmoyzRwZCZDblEs8kTNbOEiGpP30fEamD1pGXLWn4+Cziryz6+C3y3ul6ZWYqqziczsyqknE2eqJmlJtGwMbMR4HwysyZKNJs8UTNLSYDGB90JM7M2nE9m1kQJZ5MnamapSfSokJmNAOeTmTVRotnkiZpZYlK9ztrMhp/zycyaKNVs8kTNetapjOz4YwvbLp9x+IP96s7oSTRszPqtUzb9zqsvbLv8pjsv7ld3RpPzqXIP3vWTwrL6RXott9/r/qfSxqhzPg1YotnkiZpZSoJkw8bMhpzzycyaKOFsmjHoDphZbzRe7mVmVreqsknSEkn3S9os6YI2nx8l6XuSnpP0kUmfPSLpbkkbJK2r5puZWcqqHDvVmU8+o2aWmFSvszaz4VdFPkmaCVwOnAKMAWslrYqIe1tWewr4EPDWgt2cHBE7pt8bMxsGVY2d6s4nn1EzS02UfJmZ1a2abDoO2BwRD0XE88A1wNIXNBOxPSLWAjsr67uZDa/qxk615lPXiZqk+ZK+I2mTpI2SPpwvP0TSGkkP5n8ePN3OmFkXZYNmBCZqziazhimfTbMlrWt5nT1pT3OBR1vej+XLeunJzZLWt9l3LZxPZg3S29ipUflU5tLHXcAfR8Sdkg4E1ktaA7wP+HZEXJJfn3kB8NEeOmpDqKi6o6tBVkP5ywBnk/WgqHqaq61Vp4d82hERi7vsarJeDj8dHxFbJR0KrJF0X0Tc0sP2VagsnxYuWsCN6/ausNipUmPRZ0WVGjtVcCzaV69tWDHnU3/1OHZqVD51PaMWEdsi4s7852eATWQzx6XAVflqV1F8HaaZVcln1ABnk1kjVZNNY8D8lvfzgK2luxCxNf9zO7CS7FKlWjmfzBqmurFTrfnU0z1qkl4OvAq4AzgsIrbljW0DDu1lX2Y2Na76uDdnk1kzVJRNa4GFko6UtB9wGrCqVPvSrPwMFpJmAW8C7pnat6mG88ls8CocO9WaT6WrPkp6MfB14LyI+JlU7iRifv3l2QALFiwo25yZFRmBs2W9cDaZNUgF+RQRuySdC4D/s1UAABeOSURBVNwEzARWRMRGSefkny+TdDiwDjgIGJd0HnA0MBtYmefAPsDXIuJb0+/V1DifzBqiorFT3flUaqImaV+yoPlqRNyQL35c0hERsU3SEcD2gi+0HFgOsHjxYg8xzaYjXJ6/lbPJrEEqzKeIWA2snrRsWcvPj5FdcjTZz4BXVtOL6XE+mTVExWOnOvOpTNVHAV8CNkXEZ1o+WgW8N//5vcA3e2nYzKbI96gBziazRnI2Ac4ns8ZJdOxU5oza8cC7gbslbciXfQy4BLhW0pnAT4B39KeLNgx6rQbZaZtR5zNqezibbNp6rbbWaRtzPrXoez7VVamx6LNe2+jWju3N+VSdVLOp60QtIm6luKrlG6rtjpl1VWHYSFoCXEZ2nfUVEXHJpM/fxS9KRz8L/EFE/FDSfODLwOHAOLA8Ii6rrmfdOZvMGijRwVDVnE9mDZNoNpUuJmJmDRDVVXSUNBO4HDiFrNzsWkmrIuLeltUeBl4fEU9LOpXsnonXUPCMoEnbmtkoqTCfzMwqk3A29VSe38waoLrrrI8DNkfEQxHxPHAN2TN+ftFUxG0R8XT+9nbym2M7PCPIzEZZgveAmNkIGOJ71MysIURP11nPlrSu5f3yvJLYhLnAoy3vx8jOlhU5E7hxrz698BlBZjaieswnM7NapJxNnqiZpaZ82OyIiMUdPm93/0TbvUs6mWyidsKk5S94RlDpnpnZcEp0MGRmQy7RbPJEzSwxisrSZgyY3/J+HrB1r/akRcAVwKkR8WTL8nbPCDKzEVZhPpmZVSbVbPJEzQaqUwn+otL9I122v9obYtcCCyUdCWwBTgPOaF1B0gLgBuDdEfFAy/KiZwSZDYVOJa6LSmOPfFnshG/YHza9lsHvVFK/qjbqbGfYOZ96lHA2eaJmlpqKDgpFxC5J5wI3kZXnXxERGyWdk3++DPg48FLg89ncjF355ZRtnxEUEaur6Z2ZJSnNg9ZmNuwSzSZP1MwSU+UNsfnEavWkZctafj4LOKvNdp2eEWRmIyrVG/bNbLilmk2eqJmlJtGwMbMR4HwysyZKNJs8UTNLSaR7VMjMhpzzycyaKOFs8kTNLDWJho2ZjQDnk5k1UaLZ5ImaNVZRdceiapCdthkWAjSeaNqYDYmi6mlF1dY6bTNMnE/p6lRZsahSY9HyTvsq+qxTNciptDPKnE97SzmbPFEzS0yqp+/NbPg5n8ysiVLNJk/UzFISJHv63syGnPPJzJoo4WzyRM0sMak+tNHMhp/zycyaKNVs8kTNLDWJHhUysxHgfDKzJko0mzxRM0tJpHtDrJkNOeeTmTVRwtnkiZpZYlK9IdbMhp/zycyaKNVs8kTNktOpBH9R6f6hKtufaNiYDbtOJa6LSmMPXVls59PQ6bWkfqdS+0X7quvxAKNs5PMp0WzyRM0sISLdo0JmNtycT2bWRClnkydqZimJyF5mZk3jfDKzJko4mzxRM0tMqkeFzGz4OZ/MrIlSzSZP1MwSk+qzQMxs+DmfzKyJUs2mGd1WkDRf0nckbZK0UdKH8+UXSdoiaUP+enP/u2s24gIYj3KvIedsMmuYsvk0ApxPZg2S8NipzBm1XcAfR8Sdkg4E1ktak392aUR8un/dM+tNUXXHoaoG2bwcGRRnkyWjqHra0FVbcz5NGPp8mkp1xU4VIatqp442hs1I5FOi2dT1jFpEbIuIO/OfnwE2AXP73TEza09R7lVqX9ISSfdL2izpgjafv0vSXfnrNkmvLLttvzmbzJqnxmw6StL3JD0n6SO9bFsH55NZs9Q8dqosn7pO1Cbt/OXAq4A78kXn5gO4FZIO7mVfZjZFE9WLur26kDQTuBw4FTgaOF3S0ZNWexh4fUQsAj4BLO9h29o4m8waor5segr4EPDpKWxbK+eTWQPUO3aqLJ9KT9QkvRj4OnBeRPwM+ALwCuBYYBvwFwXbnS1pnaR1TzzxRNnmzKxAhUeFjgM2R8RDEfE8cA2wtHWFiLgtIp7O394OzCu7bV2cTWbNUWM2bY+ItcDOXretk/PJrBlqHjtVlk+lJmqS9iULmq9GxA15Jx6PiN0RMQ58MW98LxGxPCIWR8TiOXPmlGnOzAooQONR6lXCXODRlvdjdL4050zgxilu2xfOJrPmKJtPJUwnXxqRTeB8MmuKAY+dprVt12IikgR8CdgUEZ9pWX5ERGzL374NuKdkJ81sOsqXmJ0taV3L++URsbzlvdps0zalJJ1MNlE7oddt+8XZZNZA5fKpsmxqY+DZBM4ns8YZwNipjZ63LVP18Xjg3cDdkjbkyz5Gdl3lsXkDjwAfLNlJs9r1Wg2y0zaDphLXUOd2RMTiDp+PAfNb3s8Dtu7VnrQIuAI4NSKe7GXbPnM2WfJ6rbbWaZsmKJlPlWRTH7atkvOpjaLqip0qNRZ9VrSvThUci/bVaxujYpjyqe6xU1Xbdp2oRcSttJ8Bri7ZKTOrSlDlseG1wEJJRwJbgNOAM1pXkLQAuAF4d0Q80Mu2/eZsMmuY6vJpOvky8GwC55NZo9Q8dqpy2zJn1MysMcpVJSq1p4hdks4FbgJmAisiYqOkc/LPlwEfB14KfD67kodd+X0TbbetpGNmlqhq8qlMNkk6HFgHHASMSzoPODoifuZsMrMXqnfsVGU+eaJmlpiSN7uWEhGrmXSEN5+gTfx8FnBW2W3NbLRVlU8lsukxflGFtuu2Zjbaah47VZZPnqiZpSRA5W+INTOrj/PJzJoo4WzyRM0sNRWdvjczq5zzycyaKNFs8kTNLDVpZo2ZjQLnk5k1UaLZ5ImajbROJfiLSvcPumx/DyVmzSxRnUpcF5XGbkJZbOeTTUVdJfV7fTxAp8cGjHLp/hTzKdVs8kTNLDWJho2ZjQDnk5k1UaLZ5ImaWUIUgXanGTZmNtycT2bWRClnkydqZqlJ9KiQmY0A55OZNVGi2eSJmllqEg0bMxsBzicza6JEs8kTNbOUBJDos0DMbMg5n8ysiRLOJk/UzAoUVXc8ZcY72i5fM35dP7uzR6qVi8ysGkXV0wadTeB8sur1Wl2xU6XGqtqwYk3Np1SzyRM1s9QkGjZmNgKcT2bWRIlmkydqZimJgPFEz9+b2XBzPplZEyWcTZ6omaUmzawxs1HgfDKzJko0mzxRM0tMqtdZm9nwcz6ZWROlmk2eqJmlJtGwMbMR4HwysyZKNJs8UTPrUVGFoqKKRp226VkA42mGjZn110CzCZxP1gidKjgWVYQsWu5qkNXx2GlqPFEzS0q6N8Sa2bBzPplZE6WbTTMG3QEz61FEuVcJkpZIul/SZkkXtPn8KEnfk/ScpI9M+ux8SRsl3SPpakn7V/QNzSxVFWWTmVmlKhw71ckTNbOUTJy+L/PqQtJM4HLgVOBo4HRJR09a7SngQ8CnJ207N1++OCKOAWYCp037+5lZusrmk5lZnSocO9XNlz6aJSUgKjt9fxywOSIeApB0DbAUuHdPaxHbge2SfrfN9vsAL5K0EzgA2FpVx8wsRZXmk5lZRdLNJp9RM0tNdafv5wKPtrwfy5eV6EJsITvL9hNgG/APEXFzj9/EzIZNgpcWmdkIGNZLHyXtL+n7kn6Y349ycb78EElrJD2Y/3lw/7trNuJ6O30/W9K6ltfZk/amgha6yn/flwJHAi8DZkn6/Wl8s545m8waxpc+7uF8MmuQIb/08TngtyPiWUn7ArdKuhF4O/DtiLgkL0JwAfDRPvbVrNE6lZHtVH62Z+UrF+2IiMUdPh8D5re8n0f5yxffCDwcEU8ASLoB+E3gr8t2rgLOJrMSassmSLayWh84nxqoqNx+r2X7O+3LetPAsVOjdD2jFpln87f75q8gO5p+Vb78KuCtfemhmbUoeeq+3On7tcBCSUdK2o+sGMiqkh35CfBaSQdIEvAGYNOUvtIUOZvMmibNS4v6wflk1iSVjp1qVaqYSF4dbj3wK8DlEXGHpMMiYhtARGyTdGgf+2lmkJ++r+aoUETsknQucBNZ1cYVEbFR0jn558skHQ6sAw4CxiWdBxydZ8D1wJ3ALuAHwPJKOtYDZ5NZg1SYT8PA+WTWEAlnU6mJWkTsBo6V9BJgpaRjyjaQ3xdzNsCCBQum1Ekza1HhEZ+IWA2snrRsWcvPj5FdEtlu2wuBCyvrzBQ4m8wapoFHpAfF+WTWIIlmU09VHyPip8B3gSXA45KOAMj/3F6wzfKIWBwRi+fMmTPN7ppZqqfv+8nZZNYQzqa9OJ/MGiDRsVOZqo9z8qNBSHoRWRGB+8juZXlvvtp7gW/2q5NmNqFk1aIGVi6qmrPJrGmcTROcT2ZNku7Yqcylj0cAV+XXWs8Aro2Iv5H0PeBaSWeSFRboWppl/fr1OyT9OH87G9gxxX5XYZTbH+Xv3rT2/3lPWwbE7t2VdyhR/com8O+n2x+9tie331s2QaX5JGkJcBnZ/bNXRMQlkz5X/vmbgX8C3hcRd+afPQI8A+wGdnWpftsvHjsNV9t7tS99bqDtD0CT2h/o2KnOfOo6UYuIu4BXtVn+JFmlt9IiYs/5e0nrBhSeI9/+KH/3oWi/gafmB6Ff2QT+/XT7/m8/ZRXkUz65uRw4hewxImslrYqIe1tWOxVYmL9eA3wh/3PCyRExsEGlx07D1bbbH4L2Kxo71Z1PPd2jZmYDFpFVLirzMjOrU9l86u44YHNEPBQRzwPXkJW1b7UU+HJeBv924CUT936Zmb1AtWOnWvPJEzWz1CR6Q6yZjYBqsmku8GjL+7F8Wdl1ArhZ0vq8eqKZjbrqxk615tMgJ2q1P3Mp1fYlPdt9rbbbfWy6bU/a3+qJm6OnKZm/+ya2H+PjpV42LYP8N9K1bUmHS7pG0o8k3Zv/bv5q1e1LulLSv81/vkLS0VPZoaSTJP1mr+3n214maYukUv+/knRbzx3kBXnZ6P/2TW+/ZDbNlrSu5TV5sKJ2u+5hneMj4tVklx/9oaQTp/WlmiX5fyNlSdotacPEC/h6h3Xfp/ymMUkXSfpIh3V/KOnqkn1YLOm/5W97Gbe9XNIZZdcvaWT+2/ej/R7GTo3Kp1LPUeuHiBjof/ARaf9jwJ+1Wf5FSTMioqfRfES8uYpOjcjffX/aj4DdnoT12yD/jXRrO79JeSVwVUScli87FjgMeGCq7UqaGRG7i9qPiLOmum/gJOBZoOskqrX9fHL2NrIjkyeSlTjvtn3ZCeFkHwP+bPL3z/++1WteTkXS2ZTtoGw+7ehyr8kYML/l/Txga9l1ImLiz+2SVpJdqnRLmY41XfL/Rnrz84g4tsodSvoXZCcpTpQ0KyL+sdP6EbEOWJf/3Mt3fzlwBvC1Nn3YJyJ29bCvib6M0n/7atvvbezUqHzypY8JkXSEpFvyo0v3SPqtfPnpku7Ol30qX3YJ8KJ83a/mR3c2Sfo8cCcwX9J/zbe5W9I78+1OyttYmR+pXzZxJFvSI5Jm5z+/R9Jd+ZGprwzkL2RUxXi5lw2rk4Gdkx5MviEi/pcy7X6vi5afJOk7kr4G3J2v97n8d/9vgUMn2pD0XUmL85+flfTJ/Pf/dkmH5cv/taQ7JP1A0v+UdJiklwPnAOfnefRbykqXf13S2vx1fIfveg/Zjdint/TlIkkr8j49JOlDLZ89m/85sLyUdKakS1v69AFJn+n9P3WCqsmmtcBCSUdK2g84jaysfatVwHvyf7OvBf4hIrZJmiXpQABJs4A3kf0bsiEwaRyyWNJ3e9zFGcBXgJuBf9Oy3+9K+pSk70t6oCUvTpL0N/nPr9cvzvD9QNKBRdkKXAL8Vr7u+crO+F0n6X+QXfZ2iKRv5OOo2yUtytu4SNJXJP2dpAclfSBf/hVJS1v6+1VJe/pvJVQ3dqo1nwZ2Rs2m5Azgpoj4pLKqMwdIehnwKeA3gKfJAuCtEXGBpHMnjkblg6VfA94fEf+vpN8DjgVeSVbydK2kiRn9ccDRwI+BbwFvB66f6ISkfwn8e7LTtzskHdLvL26ZAKKBz/mwWh0DrC/47O20/73+zYLlkP2+HxMRD0t6O1lO/DrZGbp7gRVt2pkF3B4R/17SfwE+APxn4FbgtRERks4C/iQi/ljSMuDZiPg0gLKJ4aURcaukBcBNwL9o087pwNVkz5r6M0n7RsTO/LOjyCZyBwL3S/pCy2cw2Ly8BrhL0p/kfXo/8ME232+oVJVPEbFL0rlk/y5mAisiYqOkc/LPlwGryUpfbyYrf/3+fPPDgJWSIBvjfC0ivjXtTtkgvEjZJY8AD0fE2yrY5zvJqvX9GnAuWb5M2CcijpP0ZuBCsmfftfoI8IcR8b8lvRj4PxRn7gXARyLiLZBdmgm8DlgUEU9J+izwg4h4q6TfBr6c7wdgEfBaspz9gbKDZlcA5wPflPRLZJk+8Tw+66LKsVPd+eSJWlrWAisk7Qt8IyI25L/g342IJyA7ykJ2idA32mz/47z6DMAJwNURsRt4XNLfA/8K+Bnw/Yh4KN/f1fm617fs57eB6ydKi0bEU1V/USsQ4bNl1knR73W33/eH8+1PbFlvq6S/K2jneeBv8p/Xkw18ILu84/9XVt1qP+DhNttCNgA6Ov+fFcBBkg6MiGcmFuRHKt8MnB8Rz0i6g+zo49/mq/xtRDwHPCdpO9n/AMda2hhYXkbE9fnf3VskbQL2jYi7C/4uhkeF+RQRq8kGO63LWs8iB/CHbbZ7iGzQbOmr9NJHSf8KeCIifixpjCwfDo6Ip/NVbsj/XE926eJk/xv4TJ4bN0TEmKRO2TDZmpbx0gnA7wFExN9Jemk+AQP4ZkT8HPi5pO8Ax0XENyRdLulQssnh16dy+eTIqnjsVGc++dLHhETELWSDii3AVyS9h/Y3LBZpvRa703aTDzu0u0nSp3UGJMaj1MuG1kayM0LtFP1ed/p9n3yPRpl/PDvz/xFB9tDOiYN+nwU+FxG/TnYGaf+C7WcAr4uIY/PX3NZJWm4J8Etkl2Q+QjawOb3l8+dafm7tQ/YlBp+XVwDvIzuS+lc9tJs0Z5P12S5+MXYtypcipwNH5XnyI+Ag8slSbiJT9soTgMgeanwW8CLgdklHUW2mxKQ/Jy//CvAuRixTqpLq2Mln1BIi6Z8DWyLii/m1ra8mu4znMmXXbD9NFkSfzTfZOelSoVa3AB+UdBVwCNmA5t+RXU50nKQjyS7leSd7V9r5Ntmp20sj4klJh/isWj2e4emb/uf4tbNLrj6wh71aX/0d2WWAH4iIL8KeI8UHUPx7vU/B8qMm7Xti+y+T3Z92Mm1uhu/gl8gmRvDCy3KeIRsUTbiZ7LKj/5r3/9iI2MALnQ6cFRFX5+vMAh6WdECZjgw6LyPiDknz83YXlelz6nrIJ2eTTdUjZAeqbuSFk6yOlN1r/w6ySw+35MtOBv4D2UGVMvt4RX5m/G5JryP7/S/Khrlkl2UXuYVs0vUJSSeRFbD4WX6VwVJJf0526eNJZJdRAlwJfB94LCI2lvvmBmmPnTxRS8tJwL+TtJOsgtp78psT/xT4DtkRmtUR8c18/eVk90ncSXZPWauVZNdL/5DsaM2fRMRj+RGi75HdCPvrZGGysnXD/FrcTwJ/L2k38AOyI8fWZxGxZNB9sMHK7/96G/CXki4gu0/iEeA8st/Xdr/XnX7fW60ku7T5brIKkn/fY/cuAq6TtAW4HTgyX/4/gOuV3Qz/R8CHgMsl3UX2/6FbyAqOAJBPxn6Hlvu6IuIfJd0K/OuSfTmJwefltcCxLZdWDTXnk9XgYuBLyh6ncUcP251IduBmS8uyW8guwS77IOLz8sndbrL7d28kuwy8XTY8CeyS9EOyCdbkDLgI+Ks8A/+JFx7Y+j7ZJd4LgE+0VAl8PL+Uut2l2tZBytmk8INxrUV+ZGfPDbBmZtZet7xUVi3u0oj4dq0dM7MkSbqIlsJLkz47gOwg2qsj4h/q7psNhu9RMzMzq5Ckl0h6gKwYgidpZjYtkt4I3Ad81pO00eIzamZmZmZmZg3jM2pmZmZmZmYN44mamZmZmZlZw3iiZmZmZmZm1jCeqJmZmZmZmTWMJ2pmZmZmZmYN44mamZmZmZlZw/xf06sfTRQfvPQAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 1080x288 with 6 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "#| label: anisotropy-matrix\n",
        "isotropic       = np.ones(mesh.nC)*4\n",
        "anisotropic_vec = np.r_[np.ones(mesh.nC)*4, np.ones(mesh.nC)]\n",
        "anisotropic     = np.r_[np.ones(mesh.nC)*4, np.ones(mesh.nC), np.ones(mesh.nC)*3]\n",
        "Mf_sig_i = mesh.get_face_inner_product(isotropic)\n",
        "Mf_sig_v = mesh.get_face_inner_product(anisotropic_vec)\n",
        "Mf_sig_a = mesh.get_face_inner_product(anisotropic)\n",
        "clim = (\n",
        "    np.min(map(np.min, (Mf_sig_i.data, Mf_sig_v.data, Mf_sig_a.data))),\n",
        "    np.max(map(np.max, (Mf_sig_i.data, Mf_sig_v.data, Mf_sig_a.data)))\n",
        ")\n",
        "fig, ax = plt.subplots(1,3, figsize=(15,4))\n",
        "def plot_spy(ax, Mf, title):\n",
        "    dense = Mf.toarray()\n",
        "    dense[dense == 0] = np.nan\n",
        "    ms = ax.matshow(dense)#, clim=clim)\n",
        "    plt.colorbar(ms, ax=ax)\n",
        "    ax.set_xlabel(title)\n",
        "list(map(plot_spy, ax, (Mf_sig_i, Mf_sig_v, Mf_sig_a), ('Isotropic', 'Coordinate Anisotropy', 'Full Anisotropy')));"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "gg0HpwxiBxFmCKJlLfvp",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "In the last figure you can see that there is coupling between the cells in x and y due to the non-coordinate aligned anisotropy. This will make the matrix system harder to solve when we try to invert it!"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "mGY80A0iX7A0gH1lrGiQ",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "# Ok, but does it work?\n",
        "\n",
        "Unfortunately the testing of the innerproducts is not that visual. We need to make sure that the integral of some made up function evaluates to the correct value.\n",
        "\n",
        "So to make something up, and integrate it, we will use **sympy**!\n",
        "\n",
        "We will evaluate the integral:\n",
        "$$\n",
        "\\int_0^1 \\int_0^1 \\sigma(x,y) \\vec{j}(x,y) \\cdot \\vec{j}(x,y) \\quad dx dy\n",
        "$$\n",
        "by making up some ridiculous functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "YHg7nx2e6BHty5xLHFmK",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "p5vGrxND1WnUpjXOfnL4",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "It is trivial to see that the answer is 42.\n"
          ]
        }
      ],
      "source": [
        "import sympy\n",
        "from sympy.abc import x, y\n",
        "\n",
        "# Here we will make up some j vectors that vary in space\n",
        "j = sympy.Matrix([\n",
        "    x**2+y*5,\n",
        "    (5**2)*x+y*5\n",
        "])\n",
        "\n",
        "# Create an isotropic sigma vector\n",
        "Sig = sympy.Matrix([\n",
        "    [x*y*432/1163,      0      ],\n",
        "    [     0      , x*y*432/1163]\n",
        "])\n",
        "\n",
        "# Do the inner product!\n",
        "jTSj = j.T*Sig*j\n",
        "ans  = sympy.integrate(sympy.integrate(jTSj, (x,0,1)), (y,0,1))[0] # The `[0]` is to make it an int.\n",
        "\n",
        "print( \"It is trivial to see that the answer is {}.\".format(ans) )"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "Q4LY68w68gHywJXJgboO",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "Wonderful. Everything seems right in the world with that answer. Next step is to sub in the locations of the grid on both the faces in x and y as well as the cell centers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "Wtyf5RJU5fa9fGZZ5yol",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "7Gzd8hQdpwr6jqp2FQFU",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "def get_vectors(mesh):\n",
        "    \"\"\"Gets the vectors sig and [jx, jy] from sympy.\"\"\"\n",
        "    f_jx  = sympy.lambdify((x,y), j[0], 'numpy')\n",
        "    f_jy  = sympy.lambdify((x,y), j[1], 'numpy')\n",
        "    f_sig = sympy.lambdify((x,y), Sig[0], 'numpy')\n",
        "    jx  = f_jx(mesh.gridFx[:,0], mesh.gridFx[:,1])\n",
        "    jy  = f_jy(mesh.gridFy[:,0], mesh.gridFy[:,1])\n",
        "    sig = f_sig(mesh.gridCC[:,0], mesh.gridCC[:,1])\n",
        "    return sig, np.r_[jx, jy]"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "PcA8DQgBBhOOq3iMGuCl",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "Using these vectors ($\\sigma$ and $\\mathbf{j}$) we can evaluate the result on our mesh"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "qddmFAn8izDKO9ynQWyi",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "N5fL7R5F7ohI7NYHfTpZ",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Numerically we get 41.18917558039555.\n"
          ]
        }
      ],
      "source": [
        "n = 5 # get's better if you add cells!\n",
        "mesh = discretize.TensorMesh([n,n])\n",
        "sig, jv = get_vectors(mesh)\n",
        "Msig = mesh.get_face_inner_product(sig)\n",
        "numeric_ans = jv.T.dot(Msig.dot(jv))\n",
        "print( \"Numerically we get {}.\".format(numeric_ans) )"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "Wmu04ObkA2uUXBGHYg32",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "That is pretty close, but to be rigorous we should look at the order of the convergence which should be $\\mathcal{O}(h^2)$.\n",
        "\n",
        "SimPEG has a number of testing functions for \n",
        "[derivatives](http://docs.simpeg.xyz/content/api_core/api_Tests.html#SimPEG.Tests.checkDerivative)\n",
        "and \n",
        "[order of convergence](http://docs.simpeg.xyz/content/api_core/api_Tests.html#SimPEG.Tests.OrderTest) \n",
        "that make this pretty simple!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "collapsed": false,
        "iooxa": {
          "id": {
            "block": "oCa9boN1s3XlZa9SvmdU",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          },
          "outputId": {
            "block": "tGj342IRsaYHRvrqwYBJ",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        },
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/opt/anaconda3/envs/simpeg/lib/python3.8/site-packages/discretize/utils/code_utils.py:247: FutureWarning: setupMesh has been deprecated, please use setup_mesh. It will be removed in version 1.0.0 of discretize.\n",
            "  warnings.warn(\n",
            "."
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "uniformTensorMesh:  Order Test\n",
            "_____________________________________________\n",
            "   h  |    error    | e(i-1)/e(i) |  order\n",
            "~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~\n",
            "   4  |  1.27e+00   |\n",
            "   8  |  3.17e-01   |   3.9942    |  1.9979\n",
            "  16  |  7.93e-02   |   3.9985    |  1.9995\n",
            "  32  |  1.98e-02   |   3.9996    |  1.9999\n",
            "---------------------------------------------\n",
            "Happy little convergence test!\n",
            "\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "\n",
            "----------------------------------------------------------------------\n",
            "Ran 1 test in 0.090s\n",
            "\n",
            "OK\n"
          ]
        }
      ],
      "source": [
        "import sys\n",
        "import unittest\n",
        "from discretize.tests import OrderTest\n",
        "\n",
        "class Testify(OrderTest):\n",
        "    meshDimension = 2\n",
        "    def getError(self):\n",
        "        sig, jv = get_vectors(self.M)\n",
        "        Msig = self.M.get_face_inner_product(sig)\n",
        "        return float(ans) - jv.T.dot(Msig.dot(jv))\n",
        "    def test_order(self):\n",
        "        self.orderTest()\n",
        "\n",
        "# This just runs the unittest:\n",
        "suite = unittest.TestLoader().loadTestsFromTestCase( Testify )\n",
        "unittest.TextTestRunner().run( suite );"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "aTIbN9BFtHIPWlwGCeG6",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "These basic operations on the mesh for inner products and differential operators are core pieces of SimPEG, and they need to be tested continuously to ensure that no one has broken them!\n",
        "\n",
        "On every single change to the SimPEG codebase many hours of testing are completed!"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "iooxa": {
          "id": {
            "block": "XnthZqZ9T4KoIpVw1HqS",
            "project": "5dxSs15hZtvQDKSzFgjP",
            "version": 1
          }
        }
      },
      "source": [
        "# Next up ...\n",
        "\n",
        "In the [next notebook](all-together-now.ipynb) we will use our knowledge of the [mesh](mesh.ipynb) and the [divergence](divergence.ipynb) implementations to bring the DC resistivity equations together and show some pretty figures!"
      ]
    }
  ],
  "metadata": {
    "interpreter": {
      "hash": "28579dcce9b4c93863c7d16f93a67ba69218f1244c98cdb9f71bbf270d0a38fb"
    },
    "iooxa": {
      "id": {
        "block": "W44o6qg1BVf2Vg2rDMof",
        "project": "5dxSs15hZtvQDKSzFgjP",
        "version": 3
      }
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.8.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}