overview.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---\n",
        "title: Getting Started\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "This notebook uses [Python 3.8](https://docs.python.org/3/) and the open source package [SimPEG](https://simpeg.xyz). [SimPEG](https://simpeg.xyz) can be installed using the python package manager PyPi and running:\n",
        "\n",
        "```sh\n",
        "conda install SimPEG --channel conda-forge\n",
        "```\n",
        "\n",
        "This tutorial consists of 3 parts, here, we introduce the problem, in [divergence.ipynb](divergence.ipynb) we build the discrete divergence operator and in [weak-formulation.ipynb](weak-formulation.ipynb), we discretize and solve the DC equations using weak formulation.\n",
        "\n",
        "**Notebooks**\n",
        "- [Mesh](mesh.ipynb)\n",
        "- [Divergence](divergence.ipynb)\n",
        "- [Weak forulation](weak-formulation.ipynb)\n",
        "- [All together now](all-together-now.ipynb)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "# DC Resistivity"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "```{figure} ./images/dc-setup.png\n",
        "Setup of a direct current (DC) resistivity survey.\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "DC resistivity surveys obtain information about subsurface electrical conductivity, $\\sigma$. This physical property is often diagnostic in mineral exploration, geotechnical, environmental and hydrogeologic problems, where the target of interest has a significant electrical conductivity contrast from the background. In a DC resistivity survey, steady state currents are set up in the subsurface by injecting current through a positive electrode and completing the circuit with a return electrode."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "## Deriving the DC equations"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "\n",
        "```{figure} ./images/dc-eqns.png\n",
        "Derivation of the DC resistivity equations.\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "Conservation of charge (which can be derived by taking the divergence of Ampere’s law at steady state) connects the divergence of the current density everywhere in space to the source term which consists of two point sources, one positive and one negative. The flow of current sets up electric fields according to Ohm’s law, which relates current density to electric fields through the electrical conductivity. From Faraday’s law for steady state fields, we can describe the electric field in terms of a scalar potential, $\\phi$, which we sample at potential electrodes to obtain data in the form of potential differences."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "## The finish line\n",
        "*Where are we going??* \n",
        "\n",
        "Here, we are going to do a run through of how to setup and solve the DC resistivity equations for a 2D problem using [SimPEG](http://simpeg.xyz). This is meant to give you a once-over of the whole picture. We will break down the steps to get here in the series of notebooks that follow..."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "# Import numpy, python's n-dimensional array package,\n",
        "# the mesh class with differential operators from SimPEG\n",
        "# matplotlib, the basic python plotting package\n",
        "import numpy as np\n",
        "import discretize\n",
        "from SimPEG import utils\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null,
        "tags": []
      },
      "source": [
        "### Mesh\n",
        "\n",
        "Where we solve things! See [mesh.ipynb](mesh.ipynb) a discussion of how we construct a mesh and the associated properties we need. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The mesh has 8000 cells.\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAW0UlEQVR4nO3df6xk5V3H8fe3u5BiW12QelMX6q6Goo1pTbv0slr1to0uu/5B1GqphkZSs2BbfyT+AZioM/4ImCbaklqXDRLSqoBWtNu6Lmk0V0xYdEssWFq3bqHCBkx/MVWpWhe+/jGD3r079+7s7Jz7nGfO+5WQPGfO2eF5vpD57sznOTORmUiSuusFpScgSSrLRiBJHWcjkKSOsxFIUsfZCCSp4zaXnsCZuvDCC3Pbtm1T/dlnnnmGF73oRbOdUMu55m5wzd1wNmt+8MEHv5iZLx13rrpGsG3bNj7+8Y9P9WeXl5dZWlqa7YRazjV3g2vuhrNZc0T8y1rn/GhIkjrORiBJHWcjkKSOsxFIUsfZCCSp4xprBBFxe0R8PiI+ucb5iIhbIuJYRDwcEa9pai6SpLU1+Y7gDuCKdc7vBi4Z/bMX+L0G5yJJWkNj9xFk5n0RsW2dS64EPpDD78F+ICK2RMTLMvOpJuaz7Ya/AGDx6GH+7rEvD8fbLwA46Xgux6vWPMm4+rpMsebq69LQmltdl4JrLlWXzy0xcyVvKNsKPLHi+PjosVMaQUTsZfiugYWFBZaXl6f+lw4Gg7Hj9c51edyWebRt3JZ5tG3clnm0bTzL5zqb17+1lGwEMeaxsb+Sk5n7gf0AO3bsyGnurFs8epjBYMC91+/mLbceBuDua3cCnHQ8b+Nxa55kXHNdpl1zzXVpcs1trUvpNZeoy2AwaORu6mjyF8pGHw19NDO/c8y5W4HlzLxzdHwUWDrdR0M7duzIab5i4v8+GmrLW9qWj62LdbEu7azL527+IaYREQ9m5o5x50puHz0AvG20e+hy4CtN5QOSpLU19tFQRNwJLAEXRsRx4FeBcwAycx9wENgDHAO+ClzT1FwkSWtrctfQW09zPoF3NvXvlyRNZlOv1ys9hzOyf//+3t69e8/4zx3+7Jd4cfw3H/75N3D4s1/iovPP4+5rd/JjOy4+6XjexuPWPMm45rpMu+aa69Lkmttal9JrLlGXF8d/c/XrL5nqtbPf7z/V6/X2jzvXaFjcBMPiekOueRhbF+tSui7zFhZLklrARiBJHWcjkKSOMyxuafhXOlCruS6GxfWHojWs2bC4IMPiekOueRhbF+tSui6GxZKkmbMRSFLH2QgkqeMMi1sa/pUO1Gqui2Fx/aFoDWs2LC7IsLjekGsextbFupSui2GxJGnmbASS1HE2AknqOBuBJHWcjUCSOs7toy3dDlh6i13NdXH7aP3bJGtYs9tHC3L7aL3b3uZhbF2sS+m6uH1UkjRzNgJJ6jgbgSR1nGFxS8O/0oFazXUxLK4/FK1hzYbFBRkW1xtyzcPYuliX0nUxLJYkzZyNQJI6zkYgSR1nWNzS8K90oFZzXQyL6w9Fa1izYXFBhsX1hlzzMLYu1qV0XQyLJUkzZyOQpI5rtBFExBURcTQijkXEDWPOf0NEfCQiHoqIRyLimibnI0k6VWNhcURsAg4Bu4CbgFv6/f59vV7vC89f0+/3fxH498z80X6//yfAR/r9/u/0er1n13pew2LDYsPi9gWnhsWGxWNFxE6gl5m7Rsc3AmTmTSuuuRG4GHgnsA34GPCKzHxurec1LK435JqHsXWxLqXr0kRYvHmqZ5zMVuCJFcfHgcVV17wPOAA8CbwEeMu4JhARe4G9AAsLCywvL089qcFgMHa83rkuj9syj7aN2zKPto3bMo+2jWf5XGfz+reWJhtBjHls9duPXcAngDcC3wZ8LCL+NjP/7aQ/lLkf2A/DdwRLS0tnPptDw3cEW7Zsgae//P9jOPnYsXWxLtalxXWZ6vXvNJoMi48z/NjneRcx/Jv/StcA9+TQMeAx4NsbnJMkaZUmG8ER4JKI2B4R5wJXMfwYaKXHgTcBRMQCcCnwaBOTWdx+AZee/wLuvnYni9svYHH7Bdx97c5TjudtPG7Nk4xLz7vEmmuuS5NrbmtdSq+5RF0uPb+Zl+xG7yyOiD3Ae4BNwO2Z+ZsRcR1AZu6LiG8G7gBexvCjpJsz8w/We07D4npDrnkYWxfrUroutYXFZOZB4OCqx/atGD8J/GCTc5Akrc87iyWp42wEktRxfg11S+8ULX33Zc118c7i+u+grWHN3llckGFxvSHXPIyti3UpXRe/hlqSNHM2AknqOBuBJHWcYXFLw7/SgVrNdTEsrj8UrWHNhsUFGRbXG3LNw9i6WJfSdTEsliTNnI1AkjrORiBJHWdY3NLwr3SgVnNdDIvrD0VrWLNhcUGGxfWGXPMwti7WpXRdDIslSTNnI5CkjrMRSFLH2QgkqeNsBJLUcW4fbel2wNJb7Gqui9tH698mWcOa3T5akNtH6932Ng9j62JdStfF7aOSpJmzEUhSx9kIJKnjDItbGv6VDtRqrothcf2haA1rNiwuyLC43pBrHsbWxbqUrothsSRp5mwEktRxNgJJ6jjD4paGf6UDtZrrYlhcfyhaw5oNiycUEVcA7wU2Abdl5s1jrlkC3gOcA3wxM79/vec0LK435JqHsXWxLqXr0kRYvHmqZ5zsX7oJ+F3gB4DjwJGIOJCZn1pxzRbg/cAVmfl4RHxTU/ORJI3XZEbwOuBYZj6amV8D7gKuXHXNTwD3ZObjAJn5+QbnI0kao7F3BMBW4IkVx8eBxVXXvAI4JyKWgZcA783MD6x+oojYC+wFWFhYYHl5eepJDQaDseP1znV53JZ5tG3clnm0bdyWebRtPMvnOpvXv7U02QhizGOrA4nNwGuBNwHnAYcj4oHM/MxJfyhzP7AfhhnB0tLSGU9m8ehhBoMB916/m7fcehiAu6/dCXDS8byNx615knHNdZl2zTXXpck1t7Uupddcoi6DwYBpXv9Op7GwOCJ2Ar3M3DU6vhEgM29acc0NwAszszc6/n3gUGb+yVrPa1hcb8g1D2PrYl1K16W2O4uPAJdExPaIOBe4Cjiw6poPA98bEZsj4usYfnT06QbnJElapbGPhjLzRES8C7iX4fbR2zPzkYi4bnR+X2Z+OiIOAQ8DzzHcYvrJpuYkSTpVkxkBmXkQOLjqsX2rjt8NvLvJeUiS1uadxS29U7T03Zc118U7i+u/g7aGNXtncUGGxfWGXPMwti7WpXRdaguLJUkVsBFIUsfZCCSp4wyLWxr+lQ7Uaq6LYXH9oWgNazYsLsiwuN6Qax7G1sW6lK6LYbEkaeZsBJLUcTYCSeo4w+KWhn+lA7Wa62JYXH8oWsOaOxMWR8TXAy/NzM+uevxVmfnwVLM5S4bF9YZc8zC2LtaldF02NCyOiB8H/gn404h4JCIuW3H6jqlmIklqnfUygl8CXpuZ3wVcA3wwIn5kdG7cr49Jkiq03tdQb8rMpwAy8+8j4g3ARyPiIk79yUlJUqXWDIv7/f7V/X5/udfrPQ3Q6/X+o9/v3wn0gB29Xu83NmyWKxgWGxYbFrcvODUsntOwOCJeDXwVOCczP7Xi8XOAqzLzg1PN5iwZFtcbcs3D2LpYl9J12dCwODMfysx/Bv44Iq6PofOA3wbeMdVMJEmtM8kNZYvAxcD9DH+Q/knge5qclCRp40zSCP4H+E/gPOCFwGOZ+Vyjs5IkbZhJGsERho3gMuD1wFsj4kONzkqStGHW2z76vLdn5vPp7L8CV0bE1Q3OSZK0gU77XUO9Xu/JMY8V+XoJcPuo20fdPtrGrZRuH53T7aNt5fbRere9zcPYuliX0nXxh2kkSTNnI5CkjrMRSFLH+cM0LQ3/SgdqNdfFsLj+ULSGNRsWF2RYXG/INQ9j62JdStfFsFiSNHM2AknquEYbQURcERFHI+JYRNywznWXRcSzEfHmJucjSTpVY2FxRGwCDgG7gJuAW/r9/n29Xu8LY677AHAU+Ider/epU55sBcNiw2LD4vYFp4bFhsVjRcROoJeZu0bHNwJk5k2rrvsFht9wehnw0cxc9wvtDIvrDbnmYWxdrEvpujQRFk/ypXPT2go8seL4OLC4amJbgR8G3siwEYwVEXuBvQALCwssLy9PPanBYDB2vN65Lo/bMo+2jdsyj7aN2zKPto1n+Vxn8/q3liYbQYx5bPXbj/cA12fmsxHjLh/9ocz9wH4YviNYWlo689kcGr4j2LJlCzz95f8fw8nHjq2LdbEuLa7LVK9/p9FkIzjO8JfNnncRw183W2kHcNeoCVwI7ImIE5n55w3OS5K0QpO7ho4Al0TE9og4F7gKOLDygszcnpnbMnMb8CHgHU01gcXtF3Dp+S/g7mt3srj9Aha3X8Dd1+485XjexuPWPMm49LxLrLnmujS55rbWpfSaS9Tl0vObeclu9M7iiNjD8OOfTcDtmfmbEXEdQGbuW3XtHRgWt2ZsXayLdWlnXWoLi8nMg8DBVY/tW+Pan2pyLpKk8byzWJI6zkYgSR3n11C39E7R0ndf1lwX7yyu/w7aGtbsncUFGRbXG3LNw9i6WJfSdfFrqCVJM2cjkKSOsxFIUscZFrc0/CsdqNVcF8Pi+kPRGtZsWFyQYXG9Idc8jK2LdSldF8NiSdLM2QgkqeNsBJLUcYbFLQ3/SgdqNdfFsLj+ULSGNRsWF2RYXG/INQ9j62JdStfFsFiSNHM2AknqOBuBJHWcYXFLw7/SgVrNdTEsrj8UrWHNhsUFGRbXG3LNw9i6WJfSdTEsliTNnI1AkjrORiBJHWcjkKSOsxFIUse5fbSl2wFLb7GruS5uH61/m2QNa3b7aEFuH61329s8jK2LdSldF7ePSpJmzkYgSR1nI5CkjjMsbmn4VzpQq7kuhsX1h6I1rNmweEIRcQXwXmATcFtm3rzq/E8C148O/wP4mcx8aL3nNCyuN+Sah7F1sS6l61JVWBwRm4DfBXYDrwTeGhGvXHXZY8D3Z+argF8HxnYrSVJzmswIXgccy8xHM/NrwF3AlSsvyMz7M/Pp0eEDwEUNzkeSNMbmBp97K/DEiuPjwOI6178d+MtxJyJiL7AXYGFhgeXl5aknNRgMxo7XO9flcVvm0bZxW+bRtnFb5tG28Syf62xe/9bSZCOIMY+NDSQi4g0MG8Hrx53PzP2MPjbasWNHLi0tnfFkFo8eZjAYcO/1u3nLrYcBuPvanQAnHc/beNyaJxnXXJdp11xzXZpcc1vrUnrNJeoyGAyY5vXvdBoLiyNiJ9DLzF2j4xsBMvOmVde9CvgzYHdmfuZ0z2tYXG/INQ9j62JdStelqrAYOAJcEhHbI+Jc4CrgwKqJvRy4B7h6kiYgSZq9xj4ayswTEfEu4F6G20dvz8xHIuK60fl9wK8A3wi8PyIATqzVsSRJzWgyIyAzDwIHVz22b8X4p4GfbnIOkqT1eWdxS+8ULX33Zc118c7i+u+grWHN3llckGFxvSHXPIyti3UpXZfawmJJUgVsBJLUcTYCSeo4w+KWhn+lA7Wa62JYXH8oWsOaDYsLMiyuN+Sah7F1sS6l62JYLEmaORuBJHWcjUCSOs6wuKXhX+lArea6GBbXH4rWsGbD4oIMi+sNueZhbF2sS+m6GBZLkmbORiBJHWcjkKSOMyxuafhXOlCruS6GxfWHojWs2bC4IMPiekOueRhbF+tSui6GxZKkmbMRSFLH2QgkqeMMi1sa/pUO1Gqui2Fx/aFoDWs2LC7IsLjekGsextbFupSui2GxJGnmbASS1HE2AknqOMPiloZ/pQO1mutiWFx/KFrDmg2LCzIsrjfkmoexdbEupetiWCxJmjkbgSR1nI1AkjrORiBJHddoI4iIKyLiaEQci4gbxpyPiLhldP7hiHhNk/ORJJ2qse2jEbEJOATsAm4Cbun3+/f1er0vPH9Nv9/fA+wGLgf+AXhfr9e7bb3ndfuo20fdPtq+rZRuH3X76FgRsRPoZeau0fGNAJl504prbgWWM/PO0fFRYCkzn1rred0+Wu+2t3kYWxfrUrouTWwf3TzVM05mK/DEiuPjwOIE12wFTmoEEbEX2AuwsLDA8vLy1JMaDAZjx+ud6/K4LfNo27gt82jbuC3zaNt4ls91Nq9/a2myEcSYx1a//ZjkGjJzP7Afhu8IlpaWzngyn1saFnCaP1sz19wNrrkbmlpzk2HxceDiFccXAU9OcY0kqUFNNoIjwCURsT0izgWuAg6suuYA8LbR7qHLga+slw9IkmavsY+GMvNERLwLuBfYBNyemY9ExHWj8/uAg8Ae4BjwVeCapuYjSRqvyYyAzDzI8MV+5WP7VowTeGeTc5Akrc87iyWp42wEktRxNgJJ6jgbgSR1XHW/UBYRXwD+Zco/fiHwxRlOpwauuRtcczeczZq/JTNfOu5EdY3gbETEx9f6ro155Zq7wTV3Q1Nr9qMhSeo4G4EkdVzXGsHY7+Kec665G1xzNzSy5k5lBJKkU3XtHYEkaRUbgSR13Fw2goi4IiKORsSxiLhhzPmIiFtG5x+OiNeUmOcsTbDmnxyt9eGIuD8iXl1inrN0ujWvuO6yiHg2It68kfNrwiRrjoiliPhERDwSEX+z0XOctQn+3/6GiPhIRDw0WnPV32IcEbdHxOcj4pNrnJ/961dmztU/DL/y+rPAtwLnAg8Br1x1zR7gLxn+QtrlwN+VnvcGrPm7gfNH491dWPOK6/6a4bfgvrn0vDfgv/MW4FPAy0fH31R63huw5l8Cfms0finwZeDc0nM/izV/H/Aa4JNrnJ/569c8viN4HXAsMx/NzK8BdwFXrrrmSuADOfQAsCUiXrbRE52h0645M+/PzKdHhw8w/DW4mk3y3xngZ4E/BT6/kZNryCRr/gngnsx8HCAza1/3JGtO4CUREcCLGTaCExs7zdnJzPsYrmEtM3/9msdGsBV4YsXx8dFjZ3pNTc50PW9n+DeKmp12zRGxFfhhYB/zYZL/zq8Azo+I5Yh4MCLetmGza8Yka34f8B0Mf+b2H4Gfz8znNmZ6Rcz89avRH6YpJMY8tnqP7CTX1GTi9UTEGxg2gtc3OqPmTbLm9wDXZ+azw78sVm+SNW8GXgu8CTgPOBwRD2TmZ5qeXEMmWfMu4BPAG4FvAz4WEX+bmf/W9OQKmfnr1zw2guPAxSuOL2L4N4UzvaYmE60nIl4F3AbszswvbdDcmjLJmncAd42awIXAnog4kZl/vjFTnLlJ/9/+YmY+AzwTEfcBrwZqbQSTrPka4OYcfoB+LCIeA74d+PuNmeKGm/nr1zx+NHQEuCQitkfEucBVwIFV1xwA3jZK3y8HvpKZT230RGfotGuOiJcD9wBXV/y3w5VOu+bM3J6Z2zJzG/Ah4B0VNwGY7P/tDwPfGxGbI+LrgEXg0xs8z1maZM2PM3wHREQsAJcCj27oLDfWzF+/5u4dQWaeiIh3Afcy3HFwe2Y+EhHXjc7vY7iDZA9wDPgqw79RVGvCNf8K8I3A+0d/Qz6RFX9z44RrniuTrDkzPx0Rh4CHgeeA2zJz7DbEGkz43/nXgTsi4h8ZfmxyfWZW+/XUEXEnsARcGBHHgV8FzoHmXr/8iglJ6rh5/GhIknQGbASS1HE2AknqOBuBJHWcjUCSOs5GIM1QRByKiEFEfLT0XKRJ2Qik2Xo3cHXpSUhnwkYgTWH0GwcPR8QLI+JFo+/B/87M/Cvg30vPTzoTc3dnsbQRMvNIRBwAfoPhl7v9Qc138KrbbATS9H6N4Xfh/Bfwc4XnIk3Nj4ak6V3A8IdQXgK8sPBcpKnZCKTp7Qd+GfhD4LcKz0Wamh8NSVMY/fLXicz8o4jYBNwfEW8E+gy/C//Fo2+OfHtm3ltyrtLp+O2jktRxfjQkSR1nI5CkjrMRSFLH2QgkqeNsBJLUcTYCSeo4G4Ekddz/AiYbUx8aDl/QAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Define a unit-cell mesh\n",
        "mesh = discretize.TensorMesh([100, 80])  # setup a mesh on which to solve\n",
        "print(\"The mesh has {nC} cells.\".format(nC=mesh.nC))\n",
        "\n",
        "mesh.plot_grid()\n",
        "plt.axis('tight');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "### Physical Property Model\n",
        "\n",
        "Define an electrical conductivity ($\\sigma$) model, on the cell-centers of the mesh. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEYCAYAAAB82RxTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbaklEQVR4nO3dfZRlVX3m8e9DAzaNIC+NCk2jGAmRYQS1BeJbEDS8jJHgIkvwBWXM9GJFfMmaqCQrS9caTaLJTIY4vrAqSDCJitqgooMgxkHGUV4aBeymRdrmpctGsUUEIUpX1TN/nNN6u6i659xb91SdU/181jqLOvfsu+9v3y5+vXufvfeRbSIiojt2WegAIiJiMEncEREdk8QdEdExSdwRER2TxB0R0TFJ3BERHZPEHRHRMUncEREdk8TdcZLulvSyefqs9ZKOn2Md8xZvnxgukfS+Buqt/f2M4ruMnVcSd9ROprb/g+1r5yGkTpj+vQ3y/fSWbcNfZtEtSdxRSdKuCx1DRPxGEncHSDpI0mWSfiLpLklvHaacpJWSLi+v/1TShyT9C3AI8EVJv5D0zrLs3ZLeJek24BFJu/b2DGeqq+dzzpf0A0kPS7pd0uk129mvzmdJulbSg+Uwwyt7rt0t6c8k3Sbp55I+LWlpz/XnSPp2Gc+ngd5rlvTMnvMdhlFmi2mm723a93O+pDXT2vcPkj7YE/PLZqnnHZIum/be/yXpgprf426S/qr8jG1lGy3p1jrvjw6wnaPFB8VfrjcD7wZ2B54BbAJOKq/fDbysRrklwK3A/wT2pEheL+qtY9rn3g3cAqwE9pj2WbPWVZb7I+CgMqZXA48AB872WTXi2w3YCPxF2bYTgIeBw3vqvLH8zP2ADcC55bXdgXuAPy3rOQPYBryvvG7gmT1xXNJzraqdO7Sl9xx4GvAosHdPXfcBx81Qdno9B5bf2T7l+a7A/cDzav7OfAC4vvyz2xP4KnA58IyF/n3OMZojPe72ez5wgO3/Zvsx25uAfwTOHLDcMRSJ7R22H7H9S9vfqPjsD9rebPvfp73ety7bn7W9xfaU7U8Dd5bv6adfnccBTwTeX7bta8CXgLOmxbrF9gPAF4Gje967G3CB7W221wA3VcRSq5392L4H+Dbwh+VLJwCP2r6+xnvvA66j+AsQ4GRgq+2bq94raS/grcDryz+7R4DLgP3K34lYBJK42+9pwEHlEMGDkh6k6Hk+ZcByK4F7bE8M8NmbZ3m9b12SzpZ0S08cRwLLKz6rX50HAZttT/W8dg+wouf8Rz0/P0qR6Le/94e2e/cvvqciljox1fFJfvOXy2vK87o+Dryu/Pl1wL/UfN9LgE227+x5bV92/H6i45K4228zcJftfXqOvWyfOmC5zcAhs9xonG1T9tlen7UuSU+j6OmfB+xvex9gHaD+zewb3xZgpaTe39dDgB9W1AnF8MQKSb2ff0jPz48Cy3rOn1ozJpj9+9nus8Dxkg4GTmf2xD1TPZ8Hni3pSOAVwCcqPmu7A4CfbT8p2306xb9QYpFI4m6/G4GHyhuFe0haIulISc8fsNyNFEns/ZL2lLRU0gvLaz+mGBMfJKbZ6tqTIhH9BEDSORQ97rnUeQPFmO87yxtvxwN/AFxao95vARPAW8sbrK9ix2GbW4DXlN/XycDv1YwJKr432z8BrgX+ieIv1Q2zFH1cPbZ/CayhSPY32r53+7XyBuols9S1DniupKMl7QH8DcWfx6dnizO6J4m75WxPUiSpo4G7gK3ARcCTBinXc/2ZwL3AOMWNQyj+5/7LcmjjzwaI6XF12b4d+B8UCfPHwH8E/t8c63wMeCVwStmujwBn2/5ejXofA14FvJGiJ/pqiht1272t/NwHgddS9HQrYyrV+d4+SXFDt98wyWz1fJzi+5s+TLKSWb5T22uBvwKupLg5/VTgVNvb+nx+dIx2HPqLiLaQdAjwPeCpth8qX9udYqbLs5OMd15J3BEtVI7n/z3FdML/vNDxRLs0NlQi6WJJ90taN8t1SfqgpI3lwonnNhVLRJdI2hN4CHg58J4FDidGaKa8KGk/SddIurP8775V9TQ5xn0JxfzT2ZwCHFYeq4GPNhhLRGeUc8af6GI/k9mmZEY3XcLj8+L5wL/ZPgz4t/K8r8YSt+3rgAf6FDkN+GcXrgf2kXRgU/FERCy0WfLiaRQ3oin/+4dUWMjNg1aw4wKP8fK1+6YXlLSaolfOEpY8bxl7z0uAEdFtD/OzrbYPmEsdJ710T//0gclaZW++7VfrgV/2vDRme6zibU8pV8ti+z5JT676nIVM3DMtyJjxTmnZ8DGAvbWfj9WJTcYVEYvEV72m7irZWf30gUluvPqQ6oLAkgPv/KXtVXP9zCoLmbjHKeajbncwxQq5iIjWMDDFVGW5OfixpAPL3vaBFBuK9bWQC3CuAM4uZ5ccB/x8+z8XIiLawphtnqx1DOkK4A3lz28AvlD1hsZ63JI+BRwPLJc0TjGtaTcA2xdSrOw6lWK7zkeBc5qKJSJiLkbV454lL74f+IykN1Gs0P2j2WsoNJa4bZ9Vcd3Am5v6/IiIUTBmckQLFfvkxYFu3OWRVBERFaYqN4KcX0ncERF9GJhM4o6I6Jb0uCMiOsTAtpZtxpfEHRHRh3GGSiIiOsUw2a68ncQdEdFPsXKyXZK4IyL6EpOVz7qeX0ncERF9FDcnk7gjIjqjmMedxB0R0SlT6XFHRHRHetwRER1jxOSC7oD9eEncEREVMlQSEdEhRjzmJQsdxg7a1f+PiGiZYgHOLrWOKpLeJmmdpPWS3j5sTOlxR0RUGMXNSUlHAv8FOAZ4DLhK0v+2feegdaXHHRHRhy0mvUuto8KzgOttP2p7Avg6cPowMSVxR0RUmEK1jgrrgJdI2l/SMopn7q4cJp4MlURE9FHcnKydKpdLWttzPmZ7DMD2BkkfAK4BfgHcCkwME1MSd0REH9tvTta01faqWeuyPwZ8DEDSXwPjw8SUxB0RUWFyRPO4JT3Z9v2SDgFeBfzuMPUkcUdE9DHilZOXSdof2Aa82fbPhqkkiTsiosJU9YyRWmy/eBT1JHFHRPRRbDLVrgl4SdwREX0Ysa1lS96TuCMi+rCps7hmXiVxR0T0VWtxzbxK4o6I6MOkxx0R0Tm5ORkR0SFGeZBCRESXGNhWf6+SedGuaCIiWkd5WHBERJeY0a2cHJUk7oiICm3rcTf614ikkyXdIWmjpPNnuP4kSV+UdGv5DLZzmownImJQtpjyLrWO+dJYj1vSEuDDwMsp9py9SdIVtm/vKfZm4HbbfyDpAOAOSZ+w/VhTcUVEDKK4ObnzLHk/BthoexOApEuB04DexG1gL0kCngg8wJBPhIiIaIZ2qgU4K4DNPefjwLHTynwIuALYAuwFvNr21PSKJK0GVgMsZVkjwUZEzKS4ObnzjHHP1FJPOz8JuAU4CDga+JCkvR/3JnvM9irbq3bjCaOPNCKij0l2qXXMlyY/aZwdn2B8MEXPutc5wOUubATuAn6nwZgiIgayfeVknaOKpD8tJ2Ksk/QpSUuHianJxH0TcJikQyXtDpxJMSzS617gRABJTwEOBzY1GFNExMCm2KXW0Y+kFcBbgVW2jwSWUOTFgTU2xm17QtJ5wNUUAV5se72kc8vrFwLvBS6R9F2KoZV32d7aVEwREYOyYdvUyPq4uwJ7SNoGLOPxoxC1K2mM7SuBK6e9dmHPz1uA328yhoiIuSiGSmon7uWS1vacj9keA7D9Q0n/nWKk4d+Br9j+yjAxZeVkRESFAVZObrW9aqYLkvalmBJ9KPAg8FlJr7P9r4PG067JiRERLbN9OuAIbk6+DLjL9k9sbwMuB14wTEzpcUdE9DXQUEk/9wLHSVpGMVRyIrC2/1tmlsQdEVFhFM+ctH2DpDXAtylWiH8HGBumriTuiIg+ilklo9mrxPZ7gPfMtZ4k7oiIPvLosoiIDhrFUMkoJXFHRPTRxk2mkrgjIirk0WURER1ii4kk7oiIbslQSUREh2SMOyKig5K4IyI6JPO4IyI6KPO4IyI6xIaJ0T1IYSSSuCMiKmSoJCKiQzLGHRHRQU7ijojolrbdnGzXiHtERMvYo3l0maTDJd3Sczwk6e3DxJQed0REX2JyBLNKbN8BHA0gaQnwQ+Bzw9SVxB0RUaGBMe4TgR/YvmeYNydxR0T0MeBeJcsl9T4AeMz2TM+VPBP41LAxJXFHRPTjYpy7pq22V/UrIGl34JXAnw8bUhJ3RESFEc8qOQX4tu0fD1tBEndERB8e0c3JHmcxh2ESyHTAiIhKdr2jiqRlwMuBy+cST3rcEREVRjWrxPajwP5zrSeJOyKij6I33a6Vk0ncEREVsslURETHDDAdcF4kcUdE9GHEVB6kEBHRLS3rcDc7HVDSyZLukLRR0vmzlDm+3ClrvaSvNxlPRMTAypuTdY750liPu9z96sMUcxbHgZskXWH79p4y+wAfAU62fa+kJzcVT0TE0FrW5W6yx30MsNH2JtuPAZcCp00r8xrgctv3Ati+v8F4IiKG0rYed5OJewWwued8vHyt128D+0q6VtLNks6eqSJJqyWtlbR2G79qKNyIiMczMDWlWsd8afLm5EytmP4Pjl2B51HsTbsH8C1J19v+/g5vKrZFHAPYW/u17B8tMZurt9y60CGM3EkHHbXQIcR8M7ATzeMeB1b2nB8MbJmhzFbbjwCPSLoOOAr4PhERLdG2edxNDpXcBBwm6dBy/9kzgSumlfkC8GJJu5abrxwLbGgwpoiIwbnmMU8a63HbnpB0HnA1sAS42PZ6SeeW1y+0vUHSVcBtwBRwke11TcUUETG4+b3xWEejC3BsXwlcOe21C6ed/x3wd03GERExJy0bKsnKyYiIfgyexxkjdbRrAX5ERCup5lFRi7SPpDWSvidpg6TfHSaa9LgjIqqMbqjkH4CrbJ9RTtpYNkwlSdwREVVGkLgl7Q28BHgjQLmi/LFh6spQSUREP9sX4NQ5YPn2Vd7lsbqnpmcAPwH+SdJ3JF0kac9hQkrijoioMMDDgrfaXtVzjPVUsyvwXOCjtp8DPALMuGtqlSTuiIgqU6p39DcOjNu+oTxfQ5HIB5bEHRFRQa539GP7R8BmSYeXL50I3N7nLbPKzcmIiH5Gu5z9LcAnyhklm4BzhqkkiTsioq9f33icM9u3AKvmWk8Sd0RElSx5j4jomKmFDmBHSdwREf3sZA9SiIhYFKpmjMy3JO6IiCotS9yV87glnSdp3/kIJiIiqtVZgPNU4CZJn5F0sqR2DfZERDRsFAtwRqkycdv+S+Aw4GMUu1rdKemvJf1Ww7FFRCw8M6ol7yNTa8m7bQM/Ko8JYF9gjaS/bTC2iIh26NrDgiW9FXgDsBW4CHiH7W2SdgHuBN7ZbIgREQuri7NKlgOvsn1P74u2pyS9opmwIiJapGuJ2/a7+1zbMNpwIiJaqGuJOyJiZzbfM0bqSOKOiKgyjzNG6kjijoiokB53RETXjChxS7obeBiYBCZsD7U3dxJ3REQ/ox/jfqntrXOpIIk7IqJKy4ZK8rDgiIgKmqp3AMslre05Vk+rysBXJN08w7Xa0uOOiBidrRXj1i+0vUXSk4FrJH3P9nWDfkh63BERVUa0V4ntLeV/7wc+BxwzTDhJ3BER/dTc0rXqBqakPSXttf1n4PeBdcOElKGSiIgqo7k5+RTgc+UjDXYFPmn7qmEqSuKOiKgygsRtexNw1NxrSuKOiOhL/HrGSGs0OsZdPursDkkbJZ3fp9zzJU1KOqPJeCIiBjaiMe5RaixxS1oCfBg4BTgCOEvSEbOU+wBwdVOxRETMScuegNNkj/sYYKPtTbYfAy4FTpuh3FuAy4D7G4wlImJ4LUvcTY5xrwA295yPA8f2FpC0AjgdOAF4/mwVlSuMVgMsZdnIA41mnHTQSO7DRCy4tu0O2GSPe6YNbKc3/wLgXbYn+1Vke8z2KturduMJIwswIqKWnajHPQ6s7Dk/GNgyrcwq4NJyXuNy4FRJE7Y/32BcERH1uX2zSppM3DcBh0k6FPghcCbwmt4Ctg/d/rOkS4AvJWlHROu0bKikscRte0LSeRSzRZYAF9teL+nc8vqFTX12RMQotW2Mu9EFOLavBK6c9tqMCdv2G5uMJSJiaDtT4o6I6Lx5vvFYRxJ3REQfYicbKomIWAySuCMiuiaJOyKiY1qWuPMEnIiIfka8O6CkJZK+I+lLw4aUxB0RUWW0S97fBmyYSzhJ3BERFTRV76isRzoY+E/ARXOJJ2PcEREVBphVslzS2p7zMdtjPecXAO8E9ppLPEncERH9DDYMstX2qpkuSHoFcL/tmyUdP5eQkrgjIqqMZlbJC4FXSjoVWArsLelfbb9u0Ioyxh0R0cf2lZNznVVi+89tH2z76RS7pX5tmKQN6XFHRFTSVLsmcidxR0T008AmU7avBa4d9v1J3BERFbJXSURE1yRxR0R0S3rcERFdk8QdEdEhO9lT3iMiOi9PwImI6CK3K3MncUdEVEiPOyKiS/KU94iI7snNyYiIjknijojoEpObkxERXZObkxERXZPEHRHRHVmAExHRNfZIHqQgaSlwHfAEity7xvZ7hqkriTsiospoety/Ak6w/QtJuwHfkPRl29cPWlESd0REhVEMldg28IvydLfyGKrmPCw4IqIfA1Oud8BySWt7jtW9VUlaIukW4H7gGts3DBNSetwREVXq94u32l41azX2JHC0pH2Az0k60va6QcNptMct6WRJd0jaKOn8Ga6/VtJt5fFNSUc1GU9ExDDkekddth+keFjwycPE01jilrQE+DBwCnAEcJakI6YVuwv4PdvPBt4LjDUVT0TEsDTlWkffOqQDyp42kvYAXgZ8b5h4mhwqOQbYaHsTgKRLgdOA27cXsP3NnvLXAwc3GE9ExOBGtzvggcDHy07tLsBnbH9pmIqaTNwrgM095+PAsX3Kvwn48kwXygH+1QBLWTaq+CIiKhULcOaeuW3fBjxnzhXRbOLWDK/N2HpJL6VI3C+a6brtMcphlL21X8vWMEXEorcT7Q44DqzsOT8Y2DK9kKRnAxcBp9j+aYPxREQMZRQ97lFqclbJTcBhkg6VtDtwJnBFbwFJhwCXA6+3/f0GY4mIGI4HOOZJYz1u2xOSzgOuBpYAF9teL+nc8vqFwLuB/YGPSAKY6DcHMiJi/o1mr5JRanQBju0rgSunvXZhz89/DPxxkzFERMxZy4ZKsnIyIqIf59FlERHdkx53RETHtCtvJ3FHRFTRVLvGSpK4IyL6MTvVApyIiM4Tbt0CnCTuiIgqSdwRER2TxB0R0SEZ446I6J62zSrJw4IjIvpyMVRS5+hD0kpJ/0fSBknrJb1t2IjS446I6MeMaox7Avivtr8taS/gZknX2L696o3TJXFHRFQZwUiJ7fuA+8qfH5a0geJJYUncERGjNup53JKeTvEYsxuGeX8Sd0RElfqJe7mktT3nY+WjF39N0hOBy4C3235omHCSuCMi+rFhsvZYydZ+D4ORtBtF0v6E7cuHDSmJOyKiygiGSlQ85utjwAbbfz+XujIdMCKiygimAwIvBF4PnCDplvI4dZhw0uOOiOjHwAieOWn7G4DmXBFJ3BERFQxu18rJJO6IiH7MIDcn50USd0RElewOGBHRMUncERFdUmvGyLxK4o6I6MdAy7Z1TeKOiKiSHndERJcMtOR9XiRxR0T0Y3DmcUdEdMwIVk6OUhJ3RESVjHFHRHSInVklERGdkx53RESXGE9OLnQQO0jijojoZ0Tbuo5SEndERJWWTQds9Ak4kk6WdIekjZLOn+G6JH2wvH6bpOc2GU9ExKAMeMq1jiqSLpZ0v6R1c4mpscQtaQnwYeAU4AjgLElHTCt2CnBYeawGPtpUPBERQ3H5IIU6R7VLgJPnGlKTPe5jgI22N9l+DLgUOG1amdOAf3bhemAfSQc2GFNExMA8OVnrqKzHvg54YK7xNDnGvQLY3HM+Dhxbo8wK4L7eQpJWU/TIAX71Va+Z0z8zWmo5sHWhgxixxdgmWJztWoxtAjh8rhU8zM+u/qrXLK9ZfKmktT3nY7bH5hrDdE0m7pkeijl9EKhOGcqGjwFIWmt71dzDa5fF2K7F2CZYnO1ajG2Col1zrcP2nIc2Rq3JoZJxYGXP+cHAliHKREREjyYT903AYZIOlbQ7cCZwxbQyVwBnl7NLjgN+bvu+6RVFRMRvNJa4bU8A5wFXAxuAz9heL+lcSeeWxa4ENgEbgX8E/qRG1SMfL2qJxdiuxdgmWJztWoxtgpa1S9KngG8Bh0sal/Smoepxy9bgR0REf40uwImIiNFL4o6I6JjWJu7FuFy+RpteW7blNknflHTUQsQ5qKp29ZR7vqRJSWfMZ3zDqtMuScdLukXSeklfn+8YB1Xjd/BJkr4o6dayTecsRJyDqFpG3sVcUcl26w5gCfAD4BnA7sCtwBHTypwKfJliLvhxwA0LHfcI2vQCYN/y51Pa3qa67eop9zWKG9JnLHTcI/rz2ge4HTikPH/yQsc9gjb9BfCB8ucDKFb57b7QsVe06yXAc4F1s1zvVK6oc7S1x70Yl8tXtsn2N23/rDy9nmJee9vV+bMCeAtwGXD/fAY3B3Xa9Rrgctv3Athue9vqtMnAXpIEPJEicU/Mb5iDcfUy8q7likptTdyzLYUftEybDBrvmyh6CW1X2S5JK4DTgQvnMa65qvPn9dvAvpKulXSzpLPnLbrh1GnTh4BnUSyE+y7wNrftEeeD61quqNTW/bhHtly+RWrHK+mlFIn7RY1GNBp12nUB8C7bk0VHrhPqtGtX4HnAicAewLckXW/7+00HN6Q6bToJuAU4Afgt4BpJ/9f2Q00H16Cu5YpKbU3ci3G5fK14JT0buAg4xfZP5ym2uajTrlXApWXSXg6cKmnC9ufnJ8Sh1P0d3Gr7EeARSdcBRwFtTdx12nQO8H4Xg8MbJd0F/A5w4/yE2Iiu5YpKbR0qWYzL5SvbJOkQ4HLg9S3utU1X2S7bh9p+uu2nA2uAP2l50oZ6v4NfAF4saVdJyyh2v9wwz3EOok6b7qX4FwSSnkKxu96meY1y9LqWKyq1ssdte0LS9uXyS4CLXS6XL69fSDE74VSK5fKPUvQUWqtmm94N7A98pOydTrjlO7bVbFfn1GmX7Q2SrgJuA6aAi2y3dsvhmn9W7wUukfRdiiGGd9lu9Xav5TLy44HlksaB9wC7QTdzRR1Z8h4R0TFtHSqJiIhZJHFHRHRMEndERMckcUdEdEwSd0RExyRxR0R0TBJ3RETHJHFHZ5X7e98maamkPcv9o49c6LgimpYFONFpkt4HLKXY5Gnc9t8scEgRjUvijk4r99y4Cfgl8ALbkwscUkTjMlQSXbcfxYb/e1H0vCMWvfS4o9MkXUHxJJdDgQNtn7fAIUU0rpW7A0bUUT5xZsL2JyUtAb4p6QTbX1vo2CKalB53RETHZIw7IqJjkrgjIjomiTsiomOSuCMiOiaJOyKiY5K4IyI6Jok7IqJj/j/vJMMk14xV+gAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# model parameters\n",
        "sigma_background = 1.  # Conductivity of the background, S/m\n",
        "sigma_block = 10.  # Conductivity of the block, S/m\n",
        "\n",
        "# add a block to our model\n",
        "x_block = np.r_[0.4, 0.6]\n",
        "y_block = np.r_[0.4, 0.6]\n",
        "\n",
        "# assign them on the mesh\n",
        "sigma = sigma_background * np.ones(mesh.nC)  # create a physical property model \n",
        "\n",
        "block_indices = ((mesh.gridCC[:,0] >= x_block[0]) & # left boundary\n",
        "                 (mesh.gridCC[:,0] <= x_block[1]) & # right boundary\n",
        "                 (mesh.gridCC[:,1] >= y_block[0]) & # bottom boundary\n",
        "                 (mesh.gridCC[:,1] <= y_block[1]))  # top boundary\n",
        "\n",
        "# add the block to the physical property model\n",
        "sigma[block_indices] = sigma_block\n",
        "\n",
        "# plot it!\n",
        "plt.colorbar(mesh.plot_image(sigma)[0])\n",
        "plt.title('electrical conductivity, $\\sigma$');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "### Define a source\n",
        "\n",
        "Define location of the positive and negative electrodes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEYCAYAAAB82RxTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcfElEQVR4nO3dfbRdVX3u8e9DSAxBMECQlxAUK1KRCkoEfC0vWl6uleKlQxClcu3NoBVfOtoq7XXoHbe0VXtvL/X6wkgRsS2KGlDRiyDWItciL0EBCRGJ4S0GxQgIQjHJOc/9Y62jm8M5e629z1777HXyfMZYg7P2mnvu39wn/DIz15xzyTYREdEe2812ABER0Zsk7oiIlknijohomSTuiIiWSeKOiGiZJO6IiJZJ4o6IaJkk7oiIlknibjlJd0t69ZA+a42kI2dYx9Di7RLDhZLOaaDe2t/PIL7L2HYlcUftZGr7BbavHkJIrTD5e+vl++ksOwp/mUW7JHFHJUnbz3YMEfFrSdwtIGlvSZdI+qmkuyS9o59ykpZJurS8/jNJH5H0z8C+wJcl/ULSu8uyd0t6j6Rbgcckbd/ZM5yqro7POVvSDyU9Kul2SSfVbGe3Op8v6WpJD5fDDK/ruHa3pD+TdKukn0v6rKSFHddfJOk7ZTyfBTqvWdJzO86fNIwyXUxTfW+Tvp+zJa2a1L5/kPThjphfPU09fy7pkknv/T+Szq35Pc6X9NflZ2wp22hJt9R5f7SA7RwjfFD85XoT8D5gAfAcYD1wbHn9buDVNcrNA24B/jewI0XyekVnHZM+927gZmAZsMOkz5q2rrLc7wN7lzG9AXgM2Gu6z6oR33xgHfCXZduOBh4FDuio84byM3cF1gJnltcWAPcAf1LWczKwBTinvG7guR1xXNhxraqdT2pL5znwLOBxYOeOuu4Hjpii7OR69iq/s8Xl+fbAA8ChNf/MfBC4rvzd7Qh8HbgUeM5s/3nOMZgjPe7R9xJgd9v/w/Zm2+uBfwRO6bHcYRSJ7c9tP2b7CdvfqvjsD9u+z/Z/THq9a122P297o+1x258F7izf0023Oo8Ang58oGzbN4CvAKdOinWj7QeBLwOHdLx3PnCu7S22VwE3VsRSq53d2L4H+A7we+VLRwOP276uxnvvB66h+AsQ4Dhgk+2bqt4raSfgHcCby9/dY8AlwK7ln4mYA5K4R9+zgL3LIYKHJT1M0fPco8dyy4B7bG/t4bPvm+b1rnVJOl3SzR1xHAQsqfisbnXuDdxne7zjtXuApR3nP+74+XGKRD/x3h/Z7ty/+J6KWOrEVMen+fVfLm8sz+v6FPCm8uc3Af9c832vAtbbvrPjtV148vcTLZfEPfruA+6yvbjj2Mn2CT2Wuw/Yd5objdNtyj7d69PWJelZFD39s4DdbC8GbgPUvZld49sILJPU+ed1X+BHFXVCMTyxVFLn5+/b8fPjwKKO8z1rxgTTfz8TPg8cKWkf4CSmT9xT1fNF4IWSDgJeC1xU8VkTdgcemjgp230Sxb9QYo5I4h59NwCPlDcKd5A0T9JBkl7SY7kbKJLYByTtKGmhpJeX135CMSbeS0zT1bUjRSL6KYCkMyh63DOp83qKMd93lzfejgR+F7i4Rr3fBrYC7yhvsL6eJw/b3Ay8sfy+jgN+u2ZMUPG92f4pcDXwSYq/VNdOU/Qp9dh+AlhFkexvsH3vxLXyBuqF09R1G/BiSYdI2gH4W4rfx2enizPaJ4l7xNkeo0hShwB3AZuA84Fn9FKu4/pzgXuBDRQ3DqH4n/u95dDGn/UQ01Pqsn078L8oEuZPgN8C/n2GdW4GXgccX7brY8Dptr9fo97NwOuBt1D0RN9AcaNuwjvLz30YOI2ip1sZU6nO9/Zpihu63YZJpqvnUxTf3+RhkmVM853aXg38NXA5xc3pPYETbG/p8vnRMnry0F9EjApJ+wLfB/a0/Uj52gKKmS4vTDLediVxR4ygcjz/7ymmE/6X2Y4nRktjQyWSLpD0gKTbprkuSR+WtK5cOPHipmKJaBNJOwKPAK8B3j/L4cQATZUXJe0q6SpJd5b/3aWqnibHuC+kmH86neOB/ctjBfDxBmOJaI1yzvjTXexnMt2UzGinC3lqXjwb+Ffb+wP/Wp531Vjitn0N8GCXIicC/+TCdcBiSXs1FU9ExGybJi+eSHEjmvK/v0eF2dw8aClPXuCxoXzt/skFJa2g6JUzj3mHLmLnoQQYEe32KA9tsr37TOo49qgd/bMHx2qVvenWX64Bnuh4aaXtlRVv26NcLYvt+yU9s+pzZjNxT7UgY8o7pWXDVwLsrF19uI5pMq6ImCO+7lV1V8lO62cPjnHDlftWFwTm7XXnE7aXz/Qzq8xm4t5AMR91wj4UK+QiIkaGgXHGK8vNwE8k7VX2tvei2FCsq9lcgHMZcHo5u+QI4OcT/1yIiBgVxmzxWK2jT5cBf1D+/AfAl6re0FiPW9JngCOBJZI2UExrmg9g+zyKlV0nUGzX+ThwRlOxRETMxKB63NPkxQ8An5P0VooVur8/fQ2FxhK37VMrrht4W1OfHxExCMaMDWihYpe82NONuzySKiKiwnjlRpDDlcQdEdGFgbEk7oiIdkmPOyKiRQxsGbHN+JK4IyK6MM5QSUREqxjGRitvJ3FHRHRTrJwcLUncERFdibHKZ10PVxJ3REQXxc3JJO6IiNYo5nEncUdEtMp4etwREe2RHndERMsYMTarO2A/VRJ3RESFDJVERLSIEZs9b7bDeJLR6v9HRIyYYgHOdrWOKpLeKek2SWskvavfmNLjjoioMIibk5IOAv4rcBiwGbhC0v+1fWevdaXHHRHRhS3GvF2to8LzgetsP257K/BN4KR+YkrijoioMI5qHRVuA14laTdJiyieubusn3gyVBIR0UVxc7J2qlwiaXXH+UrbKwFsr5X0QeAq4BfALcDWfmJK4o6I6GLi5mRNm2wvn7Yu+xPAJwAk/Q2woZ+YkrgjIiqMDWget6Rn2n5A0r7A64GX9lNPEndERBcDXjl5iaTdgC3A22w/1E8lSdwRERXGq2eM1GL7lYOoJ4k7IqKLYpOp0ZqAl8QdEdGFEVtGbMl7EndERBc2dRbXDFUSd0REV7UW1wxVEndERBcmPe6IiNbJzcmIiBYxyoMUIiLaxMCW+nuVDMVoRRMRMXKUhwVHRLSJGdzKyUFJ4o6IqDBqPe5G/xqRdJykOyStk3T2FNefIenLkm4pn8F2RpPxRET0yhbj3q7WMSyN9bglzQM+CryGYs/ZGyVdZvv2jmJvA263/buSdgfukHSR7c1NxRUR0Yvi5uS2s+T9MGCd7fUAki4GTgQ6E7eBnSQJeDrwIH0+ESIiohnaphbgLAXu6zjfABw+qcxHgMuAjcBOwBtsj0+uSNIKYAXAQhY1EmxExFSKm5Pbzhj3VC31pPNjgZuBvYFDgI9I2vkpb7JX2l5ue/l8njb4SCMiuhhju1rHsDT5SRt48hOM96HoWXc6A7jUhXXAXcBvNhhTRERPJlZO1jmqSPqTciLGbZI+I2lhPzE1mbhvBPaXtJ+kBcApFMMine4FjgGQtAdwALC+wZgiIno2zna1jm4kLQXeASy3fRAwjyIv9qyxMW7bWyWdBVxJEeAFttdIOrO8fh7wV8CFkr5HMbTyHtubmoopIqJXNmwZH1gfd3tgB0lbgEU8dRSidiWNsX05cPmk187r+Hkj8DtNxhARMRPFUEntxL1E0uqO85W2VwLY/pGk/0kx0vAfwNdsf62fmLJyMiKiQg8rJzfZXj7VBUm7UEyJ3g94GPi8pDfZ/pde4xmtyYkRESNmYjrgAG5Ovhq4y/ZPbW8BLgVe1k9M6XFHRHTV01BJN/cCR0haRDFUcgywuvtbppbEHRFRYRDPnLR9vaRVwHcoVoh/F1jZT11J3BERXRSzSgazV4nt9wPvn2k9SdwREV3k0WURES00iKGSQUrijojoYhQ3mUrijoiokEeXRUS0iC22JnFHRLRLhkoiIlokY9wRES2UxB0R0SKZxx0R0UKZxx0R0SI2bB3cgxQGIok7IqJChkoiIlokY9wRES3kJO6IiHYZtZuTozXiHhExYuzBPLpM0gGSbu44HpH0rn5iSo87IqIrMTaAWSW27wAOAZA0D/gR8IV+6krijoio0MAY9zHAD23f08+bk7gjIrroca+SJZI6HwC80vZUz5U8BfhMvzElcUdEdONinLumTbaXdysgaQHwOuAv+g0piTsiosKAZ5UcD3zH9k/6rSCJOyKiCw/o5mSHU5nBMAlkOmBERCW73lFF0iLgNcClM4knPe6IiAqDmlVi+3Fgt5nWk8QdEdFF0ZserZWTSdwRERWyyVRERMv0MB1wKJK4IyK6MGI8D1KIiGiXEetwNzsdUNJxku6QtE7S2dOUObLcKWuNpG82GU9ERM/Km5N1jmFprMdd7n71UYo5ixuAGyVdZvv2jjKLgY8Bx9m+V9Izm4onIqJvI9blbrLHfRiwzvZ625uBi4ETJ5V5I3Cp7XsBbD/QYDwREX0ZtR53k4l7KXBfx/mG8rVOzwN2kXS1pJsknT5VRZJWSFotafUWftlQuBERT2VgfFy1jmFp8ubkVK2Y/A+O7YFDKfam3QH4tqTrbP/gSW8qtkVcCbCzdh2xf7TEdK7ceMtshzBwx+598GyHEMNmYBuax70BWNZxvg+wcYoym2w/Bjwm6RrgYOAHRESMiFGbx93kUMmNwP6S9iv3nz0FuGxSmS8Br5S0fbn5yuHA2gZjiojonWseQ9JYj9v2VklnAVcC84ALbK+RdGZ5/TzbayVdAdwKjAPn276tqZgiIno33BuPdTS6AMf25cDlk147b9L53wF/12QcEREzMmJDJVk5GRHRjcFDnDFSx2gtwI+IGEmqeVTUIi2WtErS9yWtlfTSfqJJjzsiosrghkr+AbjC9snlpI1F/VSSxB0RUWUAiVvSzsCrgLcAlCvKN/dTV4ZKIiK6mViAU+eAJROrvMtjRUdNzwF+CnxS0nclnS9px35CSuKOiKjQw8OCN9le3nGs7Khme+DFwMdtvwh4DJhy19QqSdwREVXGVe/obgOwwfb15fkqikTesyTuiIgKcr2jG9s/Bu6TdED50jHA7V3eMq3cnIyI6Gawy9nfDlxUzihZD5zRTyVJ3BERXf3qxuOM2b4ZWD7TepK4IyKqZMl7RETLjM92AE+WxB0R0c029iCFmMLHb/oQz33RftNeX/fdu/ijQ989xIhiW7brnov5b595F+ecci4P/eTh2Q5nZFXNGBm2TAccstu//QM2P7Flymubn9jCmmvvGHJEsS077b0nc9Arns9p7/3Psx3KaBuxBylUJm5JZ0naZRjBbAsuOmcV9tQDZh4f56JzLhlyRLGt2nXPxRx7xpFsN287jjvjKHbZY/FshxQ11elx7wncKOlzko6TNFqDPS3z4I8f5spPXv2UXvfmJ7ZwxSf/Lf9cjaE57b0nIxUpQNttl153F4NYgDNIlYnb9nuB/YFPUOxqdaekv5H0Gw3HNmdN1etObzuGaaK3vWDhfAAWLJyfXvd0zKCWvA9MrTFu2wZ+XB5bgV2AVZI+1GBsc9bkXnd62zFsnb3tCel1d9HCMe53SLoJ+BDw78Bv2f4j4FAgv+U+dfa609uOYZrc256QXvf0WjdUAiwBXm/7WNuft70FwEXWeW2j0c1hE73u8bHx9LZjqKbqbU9Ir3saI9bjrpzHbft9Xa6tHWw425aLzlnFs1+wT3rbMVQHvvR5T+ltT1iwcD4veNkBU17bpo3YPO4swJlFD/74Yf70qP8+22HENiYLvHoz7GGQOpK4IyKqDHHGSB1J3BERFdLjjohomwElbkl3A48CY8BW233tzZ3EHRHRzeDHuI+yvWkmFSRxR0RUGbGhkuwOGBFRQeP1DmCJpNUdx4pJVRn4mqSbprhWW3rcERGDs6li3PrltjdKeiZwlaTv276m1w9JjzsiosqAVk7a3lj+9wHgC8Bh/YSTxB0R0U3NfUqqbmBK2lHSThM/A78D3NZPSBkqiYioMpibk3sAXygfabA98GnbV/RTURJ3RESVASRu2+uBg2deUxJ3RERX4lczRkZGo2Pc5aPO7pC0TtLZXcq9RNKYpJObjCciomcDGuMepMYSt6R5wEeB44EDgVMlHThNuQ8CVzYVS0TEjIzYftxN9rgPA9bZXm97M3AxcOIU5d4OXAI80GAsERH9G7HE3eQY91Lgvo7zDcDhnQUkLQVOAo4GXjJdReUKoxUAC1k08ECjGcfuPZD7MBGzbtR2B2yyxz3VBraTm38u8B7bY90qsr3S9nLby+fztIEFGBFRyzbU494ALOs43wfYOKnMcuDicl7jEuAESVttf7HBuCIi6vPozSppMnHfCOwvaT/gR8ApwBs7C9jeb+JnSRcCX0nSjoiRM2JDJY0lbttbJZ1FMVtkHnCB7TWSziyvn9fUZ0dEDNKojXE3ugDH9uXA5ZNemzJh235Lk7FERPRtW0rcERGtN+Qbj3UkcUdEdCG2saGSiIi5IIk7IqJtkrgjIlpmxBJ3noATEdHNgHcHlDRP0nclfaXfkJK4IyKqDHbJ+zuBtTMJJ4k7IqKCxusdlfVI+wD/CTh/JvFkjDsiokIPs0qWSFrdcb7S9sqO83OBdwM7zSSeJO6IiG56GwbZZHv5VBckvRZ4wPZNko6cSUhJ3BERVQYzq+TlwOsknQAsBHaW9C+239RrRRnjjojoYmLl5Exnldj+C9v72H42xW6p3+gnaUN63BERlTQ+WhO5k7gjIrppYJMp21cDV/f7/iTuiIgK2askIqJtkrgjItolPe6IiLZJ4o6IaJFt7CnvERGtlyfgRES0kUcrcydxR0RUSI87IqJN8pT3iIj2yc3JiIiWSeKOiGgTk5uTERFtk5uTERFtk8QdEdEeWYATEdE29kAepCBpIXAN8DSK3LvK9vv7qSuJOyKiymB63L8Ejrb9C0nzgW9J+qrt63qtKIk7IqLCIIZKbBv4RXk6vzz6qjkPC46I6MbAuOsdsETS6o5jRWdVkuZJuhl4ALjK9vX9hJQed0RElfr94k22l09bjT0GHCJpMfAFSQfZvq3XcBrtcUs6TtIdktZJOnuK66dJurU8rpV0cJPxRET0Q6531GX7YYqHBR/XTzyNJW5J84CPAscDBwKnSjpwUrG7gN+2/ULgr4CVTcUTEdEvjbvW0bUOafeyp42kHYBXA9/vJ54mh0oOA9bZXg8g6WLgROD2iQK2r+0ofx2wT4PxRET0bnC7A+4FfKrs1G4HfM72V/qpqMnEvRS4r+N8A3B4l/JvBb461YVygH8FwEIWDSq+iIhKxQKcmWdu27cCL5pxRTSbuDXFa1O2XtJRFIn7FVNdt72SchhlZ+06YmuYImLO24Z2B9wALOs43wfYOLmQpBcC5wPH2/5Zg/FERPRlED3uQWpyVsmNwP6S9pO0ADgFuKyzgKR9gUuBN9v+QYOxRET0xz0cQ9JYj9v2VklnAVcC84ALbK+RdGZ5/TzgfcBuwMckAWztNgcyImL4BrNXySA1ugDH9uXA5ZNeO6/j5z8E/rDJGCIiZmzEhkqycjIiohvn0WUREe2THndERMuMVt5O4o6IqKLx0RorSeKOiOjGbFMLcCIiWk945BbgJHFHRFRJ4o6IaJkk7oiIFskYd0RE+4zarJI8LDgioisXQyV1ji4kLZP0b5LWSloj6Z39RpQed0REN2ZQY9xbgT+1/R1JOwE3SbrK9u1Vb5wsiTsiosoARkps3w/cX/78qKS1FE8KS+KOiBi0Qc/jlvRsiseYXd/P+5O4IyKq1E/cSySt7jhfWT568VckPR24BHiX7Uf6CSeJOyKiGxvGao+VbOr2MBhJ8ymS9kW2L+03pCTuiIgqAxgqUfGYr08Aa23//UzqynTAiIgqA5gOCLwceDNwtKSby+OEfsJJjzsiohsDA3jmpO1vAZpxRSRxR0RUMHi0Vk4mcUdEdGN6uTk5FEncERFVsjtgRETLJHFHRLRJrRkjQ5XEHRHRjYER29Y1iTsiokp63BERbdLTkvehSOKOiOjG4MzjjohomQGsnBykJO6IiCoZ446IaBE7s0oiIlonPe6IiDYxHhub7SCeJIk7IqKbAW3rOkhJ3BERVUZsOmCjT8CRdJykOyStk3T2FNcl6cPl9VslvbjJeCIiemXA4651VJF0gaQHJN02k5gaS9yS5gEfBY4HDgROlXTgpGLHA/uXxwrg403FExHRF5cPUqhzVLsQOG6mITXZ4z4MWGd7ve3NwMXAiZPKnAj8kwvXAYsl7dVgTBERPfPYWK2jsh77GuDBmcbT5Bj3UuC+jvMNwOE1yiwF7u8sJGkFRY8c4Jdf96oZ/TNjRC0BNs12EAM2F9sEc7Ndc7FNAAfMtIJHeejKr3vVkprFF0pa3XG+0vbKmcYwWZOJe6qHYk4eBKpThrLhKwEkrba9fObhjZa52K652CaYm+2ai22Col0zrcP2jIc2Bq3JoZINwLKO832AjX2UiYiIDk0m7huB/SXtJ2kBcApw2aQylwGnl7NLjgB+bvv+yRVFRMSvNZa4bW8FzgKuBNYCn7O9RtKZks4si10OrAfWAf8I/HGNqgc+XjQi5mK75mKbYG62ay62CUasXZI+A3wbOEDSBklv7asej9ga/IiI6K7RBTgRETF4SdwRES0zsol7Li6Xr9Gm08q23CrpWkkHz0acvapqV0e5l0gak3TyMOPrV512STpS0s2S1kj65rBj7FWNP4PPkPRlSbeUbTpjNuLsRdUy8jbmikq2R+4A5gE/BJ4DLABuAQ6cVOYE4KsUc8GPAK6f7bgH0KaXAbuUPx8/6m2q266Oct+guCF98mzHPaDf12LgdmDf8vyZsx33ANr0l8AHy593p1jlt2C2Y69o16uAFwO3TXO9VbmizjGqPe65uFy+sk22r7X9UHl6HcW89lFX53cF8HbgEuCBYQY3A3Xa9UbgUtv3Atge9bbVaZOBnSQJeDpF4t463DB74+pl5G3LFZVGNXFPtxS+1zKjpNd430rRSxh1le2StBQ4CThviHHNVJ3f1/OAXSRdLekmSacPLbr+1GnTR4DnUyyE+x7wTo/aI85717ZcUWlU9+Me2HL5EVI7XklHUSTuVzQa0WDUade5wHtsjxUduVao067tgUOBY4AdgG9Lus72D5oOrk912nQscDNwNPAbwFWS/p/tR5oOrkFtyxWVRjVxz8Xl8rXilfRC4HzgeNs/G1JsM1GnXcuBi8ukvQQ4QdJW218cToh9qftncJPtx4DHJF0DHAyMauKu06YzgA+4GBxeJ+ku4DeBG4YTYiPalisqjepQyVxcLl/ZJkn7ApcCbx7hXttkle2yvZ/tZ9t+NrAK+OMRT9pQ78/gl4BXStpe0iKK3S/XDjnOXtRp070U/4JA0h4Uu+utH2qUg9e2XFFpJHvctrdKmlguPw+4wOVy+fL6eRSzE06gWC7/OEVPYWTVbNP7gN2Aj5W9060e8R3barardeq0y/ZaSVcAtwLjwPm2R3bL4Zq/q78CLpT0PYohhvfYHuntXstl5EcCSyRtAN4PzId25oo6suQ9IqJlRnWoJCIippHEHRHRMkncEREtk8QdEdEySdwRES2TxB0R0TJJ3BERLZPEHa1V7u99q6SFknYs948+aLbjimhaFuBEq0k6B1hIscnTBtt/O8shRTQuiTtardxz40bgCeBltsdmOaSIxmWoJNpuV4oN/3ei6HlHzHnpcUerSbqM4kku+wF72T5rlkOKaNxI7g4YUUf5xJmttj8taR5wraSjbX9jtmOLaFJ63BERLZMx7oiIlknijohomSTuiIiWSeKOiGiZJO6IiJZJ4o6IaJkk7oiIlvn/MF5esARH7GYAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Define a source\n",
        "a_loc, b_loc = np.r_[0.2, 0.5], np.r_[0.8, 0.5]\n",
        "source_locs = [a_loc, b_loc]\n",
        "\n",
        "# locate it on the mesh\n",
        "source_loc_inds = mesh.closest_points_index(source_locs)\n",
        "a_loc_mesh = mesh.gridCC[source_loc_inds[0],:]\n",
        "b_loc_mesh = mesh.gridCC[source_loc_inds[1],:]\n",
        "\n",
        "# plot it\n",
        "plt.colorbar(mesh.plot_image(sigma)[0])\n",
        "plt.plot(a_loc_mesh[0], a_loc_mesh[1],'wv', markersize=8) # a-electrode\n",
        "plt.plot(b_loc_mesh[0], b_loc_mesh[1],'w^', markersize=8) # b-electrode\n",
        "plt.title('electrical conductivity, $\\sigma$');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "### Assemble and solve the DC system of equations\n",
        "\n",
        "How we construct the divergence operator is discussed in [divergence.ipynb](divergence.ipynb), and the inner product matrix in [weakformulation.ipynb](weakformulation.ipynb). The final system is assembled and discussed in [play.ipynb](play.ipynb) (with widgets!)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [],
      "source": [
        "# Assemble and solve the DC resistivity problem\n",
        "Div = mesh.face_divergence\n",
        "Sigma = mesh.get_face_inner_product(sigma, invert_model=True, invert_matrix=True)\n",
        "Vol = utils.sdiag(mesh.cell_volumes)\n",
        "\n",
        "# assemble the system matrix\n",
        "A = Vol * Div * Sigma * Div.T * Vol\n",
        "\n",
        "# right hand side\n",
        "q = np.zeros(mesh.nC)\n",
        "q[source_loc_inds] = np.r_[+1, -1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "iooxa": null,
        "tags": []
      },
      "outputs": [],
      "source": [
        "from SimPEG import Solver # import the default solver (LU)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "iooxa": null,
        "tags": []
      },
      "outputs": [],
      "source": [
        "# solve the DC resistivity problem\n",
        "Ainv = Solver(A)  # create a matrix that behaves like A inverse\n",
        "phi = Ainv * q"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false,
        "iooxa": null,
        "jupyter": {
          "outputs_hidden": false
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEYCAYAAABSnD3BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2df7AtVXXnP99733uCoIJ5gvB4RJzBHxkjjoI48UdQJMJLlDjlJMSMIqNDOSPGVMYI+WVqJlMTjDNTaqGh3iBRJ07QUUafFspELDQZRYEKEuEFeILCCwh5YDQi8XHvXfNH9zl3n3t2n919Tvc53X3Xp6rr9um9e/fqPn3XWXvttdeWmeE4juP0l6VFC+A4juM0iyt6x3GcnuOK3nEcp+e4onccx+k5rugdx3F6jit6x3GcnuOK3nEcp+e4onccx+k5rug3OZK+Lenlc7rWLZJOm8e15kXZe5rnc3acjbii3wTkSuYRST8Mtktqbj+pxMzsn5nZtVO2P5D/fkl/IunwuuSqKMdIe9PeU4Vr/rGk32uqfWdz4Ip+8/BKMzs82C6Y14UlbamhmVea2eHAc4FTgN+toc0u8ELgLxcthNNtXNE7QyQdK+mTkv5O0l2Sfm1D+U5JV+blD0q6RNL/BI4HPpNb3O/I635b0oWSbgYelrQltIhjbZWR0cz+Fvgc8Ky8nWdKulbS3+dulFflx4vkKrzHXL63S7pZ0vclfUzSIYn2wnu6SNK3JP2DpFslvXrK72FJ0jskfRv4aeAKSb8xTVuOA67onRxJS8BngG8AO4DTgV+X9Iq8fBn4LPAd4Cl5nSvM7HXA3az3GP4oaPZXgJ8HjjCzleBa0bZKyrkT2AX8laStucz/FzgKeCvwUUlPj8mVusecXwLOBE4Ang28ASBxnwO+BbwYeALwH4E/lXRMmfvawG+RPbdfA24HXgm8O793x6mMK/rNw6dyq3ew/dsN5acATzKz/2RmB83sTuB/AOfk5c8HjgV+08weNrN/NLOUS+F9ZnaPmT2y4fg0bX1K0t+TuTG+BPwX4AXA4cDFucxfJPsB+ZWCNlL3OJD5XjN7iOxH4TkJuYaY2f/Oz10zs48Bd+T3WhpJjwN+GzgX2An8lZl9HbgHeEaVthxnQB2+U6cb/KKZfWFC+U8Cx+bKdMAy8Bf5/k7gO6FlXoJ7Co5P09aY/JKOBe4xs7Xg8HfIrPUYqXsE+G6w/yOyH6RSSHo98BtkvRTIfoS2lz0/52XA7Wb2bUknkfVcloAjgQcqtuU4gCt6Z517gLvM7MQJ5cdL2hJR0EWLGhQdn9RWFe4FdkpaCpT98WTujtj1U/eYonDxBkk/SdY7OB34qpmtSroJUMVrPBl4KN9/DvAJ4CXA94GbK0vsOLjrxlnn68AP8gHUQyUtS3qWpFOC8vuAiyUdJukQSS/My+4HnlrxWkVtVeFrwMPAOyRtzePZX8m6v3+jXKl7TDHpPg8j+yH4OwBJ55EPGMeQ9CFJH4oU/Q3wPElPzc9/CLgEeIf5KkHOlLii3zwMokUG2/8JC81slUxJPge4CzgAXEY2sBiW/1OyQcn9wC/np/8h8Lu57//tKUESbZXGzA4CrwLOyuX9APB6M/ubmFypeyxB4X2a2a3AfwO+SvaD8NPA/5vQ1s6C8i8DHwRuBA4BPgq8y8xKDVY7Tgy5keA480XSNrLIn2eb2aMFdc4B3mBmZ85VOKeXuI/eceZM3hN5ZqLa04Db5iCOswlozHUj6XJJD0j6ZkG5JL1P0r58gspzm5LFcTrI01kfVHacmWjMdSPpJcAPgY+Y2diglKRdZBNcdgGnAu81s1MbEcZxHGcT05hFb2ZfZj1MLMbZZD8CZmbXAUdMOYvQcRzHmcAiffQ7GJ1Qsz8/dt/GipLOB84H0LZtz9t69FFzEdBxnG5z8J79B8zsSbO08YqXHmYPPrRaqu6NN//46jYOoC9S0ccmkkT9SGa2G9gN8Jjjd9qOt/96k3I5jtMT7nrb278zaxsPPrTK168+vlTd5WPuqDoTei4sUtHvJ4slHnAc2UxHx3Gc1mDAGmvJem1mkROm9gCvz6NvXgB838zG3DaO4ziLxDAetdVSW1tpzKKX9GfAacB2SfuB3we2ApjZpcBVZBE3+8iSR53XlCyO0ymsanqchpBPphzQdYu+MUVvZkWpYgflBrylqes7juPUgWGsdjyDgM+MbZC2GGYhbqRNSRu/zKZp+p479DKuFScu7QSu6B3HcSZgwKoresdxnH7jFn1P6WtPvc776lDPO01fv/A2k3rmLXnBDHjUffSO4zj9xTB33TiO4/Qag9Vu6/nNq+i9pz47Hel5j+JffHcIv6sFvkzZzNhus2kVveM4TjnEauU13tvFplP0btDNj5YYZBlt/OIX/UxmZZ6PdIEvUzYY28L3pwKbTtE7juNUIYujd0XvOI7Ta9bcom8/Hf+OekHsO2i8B76IL77r7pgqpO61qcc/ZzeOW/SO4zg9xxCrC83oPjuu6B3HcRJ03XXT7Z+pCZjWN8dpBItszjqx59PBZ2SIg7Zcaksh6UxJt0naJ+migjqnSbpJ0i2SvlTHPbhF7ziOM4FswtTsNrGkZeD9wBlkS6leL2mPmd0a1DkC+ABwppndLemomS9Mjy16x3GculjNJ02ltgTPB/aZ2Z1mdhC4Ajh7Q53XAlea2d0AZvZAHfK7Re8sjFZNqEoxT/na+CzqdIHG7m/a9gcvUYMvkJlYtdI28XZJNwSfd5vZ7nx/B3BPULYfOHXD+U8Dtkq6Fngc8F4z+0h1qUdxRe84jpNgrfwv0QEzO7mgLNbIxl+oLcDzgNOBQ4GvSrrOzG4vK0CM3in6Vg6+LkqmNlqGBdRmmNX5AjTx/Dr0nYzQdMx82H7L/oezwdhaVOV+YGfw+Tjg3kidA2b2MPCwpC8DJwEzKXr30TuO40xgMBhbZktwPXCipBMkbQPOAfZsqPNp4MWStkh6LJlrZ++s99A7i95xHKduVmvoKZrZiqQLgKuBZeByM7tF0pvz8kvNbK+kzwM3k2VHvszMvjnrtV3R10HLuppDysrVVXdCE8zsOqpFiloI3WCNuzRrHWCd8fyaqXNmrJldBVy14dilGz6/G3h3LRfMcUXvOI6TYK181E0rcUXvOI4zgSypmSv6hTPXSJuWdCdrpeieWuSGaJwWuWyaCAmvs83S/29F1+zY/5AhHi2R3qDN9ELRO47jNIUZVSZMtRJX9EXUGY69IAtmZiuuzPSOmpnrbNlp2p9RptbPAC5B7B4qveNlB1urxNY3+uKoyoSpVuKK3nEcZwKGW/SO4zi9xwdj+8YUPbRWpl0gLddUPdywzTa5IZr+Eqa411o9CIt41hUe6VTunEqumQp1a8ZQ5xcecUXvOI4zAQMerSfXzcLotvSO4ziNUyrXfKtxRQ+VuoIz9+CaDruoIGCsaiXxBue3yYWTooqs83TXtPEZzpixsuhZRF/RVmev9JmxjuM4vafrFn2jP1OphXAlPUHSZyR9I18I97wm5XEcx6mKmVizpVJbW2nMoi+zEC7wFuBWM3ulpCcBt0n6aL6eYrMkfqCndtEsekZM6vqJG5vKndPWSJw5Uvprr+P5zOsZV4mEqXDe4FkVvoqty16Jp0CYwHAhXABJg4VwQ0VvwOMkCTgceAhYaVAmx3GcilRaM7aVNKnoyyyEewnZCiv3ki2E+8tmtraxIUnnA+cDLB95JDCLxT25uFK701jvTVgplQZQqwc8V5pdXqN1v5DFwytcJylT1wdmpx2MLZmbfup8+XMeuM0GY1vSvZiSJn+mymRKeQVwE3As8BzgEkmPHzvJbLeZnWxmJy8fflj9kjqO40xglaVSW1tpUrIyC+GeB1xpGfuAu4BnNCiT4zhOJQYzY8tsbaVJ181wIVzgb8kWwn3thjp3A6cDfyHpaODpwJ0NyhQl+f1U8Rss4ruedtBseH75PvRC3Cl1UJObpVZ3zZTPbxHPPRn7HjJFnHz0FWxRbH2Jhb9bTWPSm9kKMFgIdy/w8cFCuIPFcIE/AH5G0l8D1wAXmtmBpmRyHMepihk8urZUakuRCjkP6p0iaVXSa+q4h0YnTKUWwjWze4Gfa1IGx3GcWchcN7PbxCVDzgf13kVmJNfC5pgZG+n21eKuSbaRbqJ2SkY8FNctH5UzOFz4qHoUXz+zy6bOaJ5o+zW+bBEBKk3PSL2DLYuTL0NNM2PLhJwDvBX4JHBKHReFhmfGOo7jdJ1BeGUNg7GxkPMdYQVJO4BXA5dSI5vDonccx5maSq6b7ZJuCD7vNrPdw4bG2dj/eQ/ZWOVqNo+0HlzRbyTWR23KRTPNeTNmFExGMqTcORGrpbOROAmadtdUelbzDN2bYsWaZOBW7L0reBejKRKK3ts5uYEqrBl7wMxOLigrE3J+MnBFruS3A7skrZjZpyqIO4YresdxnAlkUTe15LpJhpyb2QmDfUkfAj47q5KHPiv6gh/gqLVRxYpf9ADstFPFY+dPVT7lvPVF5K5vIkXBPKz41HNd9LKCicx3yfH8mBXeopj5jdS1lKCZrUgahJwvA5cPQs7z8lr98iH9VfSO4zg1UcF1M5FUyPmG42+o5aK4onccx5lIH5KabQpFX9pdA/Fu49RunAXPVZ9mybbUoFmVy3dwYLappQCj7c7oopnrqpTJa5V355RNcTB1dssGaPOiImXYFIrecRxnWszEiit6x3GcfuOum66RirCp5LpJ9a1LyjQrI13gRMhDFXdOMjoitSZcj4i6XiafUsldM6ubpgE3TnJKRipqpsCPl4yTTwngC49UZvMpesdxnIq4om8b+feRTlo2ZXmFmPt5vRuFlxla5AXm3jArWeSc8HhbB1WnecAN30uh5R2TtSgeoK7ew4zUuZJgbOprldm0qTD+0fPq/cerK45+kfRP0TuO49RMXXH0i8IVveM4zgTMYKXEoiJtZnMoehX0AWPHBvsV4uwLe3VzCnQuun70cGowdaq84vE+eHQsrqU56kt/VZXcLZPdNcVunskiTOPamYYq+fFS51msRlGcfVkB5oi7bhzHcXqM++gdx3E2AeaKvqMUxc4nompmTqdQJ0PXS7zYIgJotMJoO1CQmTByzmagziX98raquGsquXaa+F5SKQpSpxWkwwgcNvEKVagSfz8DPhjrOI7TY8zcR+84jtNzxKpH3bSAmGulKNJGsWMV3DWJCJ7kMo/TdFET1kRxi+NRNdEAm8IUChovZ7x4U7lzAupcarBSJE0VN09dVMg+mXot4uUK9qzU+SOyNIz76B3HcXqM57rpI7EUCgVWfNR6T8XsTyVTdCQruEz8QhbZGy1Xfn70pAJZKtSdLNSmYvgVVomjL1N30vmzUjQ/I5Y3L3F6cinBNmOZn77LuKJ3HMdJ4FE3juM4PcZ8MHYxlHKXFcbJD/7GXSypwVzFXDqpfPYz96UTKQoKXDsDl87I1SMDrzYyEBarrPFDc6BTSxEOXpw63C2pmPvIObWmRUiMpq7/j8Svb2M7oEBPDpOmFnokBy7FAufjAozrulw3ks4E3gssA5eZ2cUbyn8VuDD/+EPg35nZN2a9bicVveM4zjypI+pG0jLwfuAMYD9wvaQ9ZnZrUO0u4GfN7HuSzgJ2A6fOem1X9I7jOBMwqy288vnAPjO7E0DSFcDZwFDRm9lXgvrXAcfVceH+Kfrk8n7j7pZYhE3Y1UzG5BemNIwcmsIfke42Tp6MHkblFLpxGqBLrpeYO6HW1A8xd0zBfsoNE4vgKao7sc2i0xMZXmOZSJOrVa4F5fn/VvJRl1nZZE5unArhldsl3RB83m1mu/P9HcA9Qdl+JlvrbwQ+V1rICfRP0TuO49RMBR/9ATM7uaAsncF5UFF6KZmif1HpK0+gF4p+mmUDi+Lkh5Z80oqHdTMtKE7OnJ1mZmzB/mD8r3AK4XgFrcWs+3AvMjBbh2XeFut+0QnaZrXiw+OpcgrKUyQs+tJ574oGayMxDNHBVgvfxeAdjbTb5HwmQ6zVE3WzH9gZfD4OuHdjJUnPBi4DzjKzB+u4cLdjhhzHceaAldwSXA+cKOkESduAc4A9YQVJxwNXAq8zs9vrkr9RRS/pTEm3Sdon6aKCOqdJuknSLZK+1KQ8juM4lckHY8tsE5sxWwEuAK4G9gIfN7NbJL1Z0pvzau8EfgL4QK4XbyhorhKNuW7KhBJJOgL4AHCmmd0t6aj6BAj3I+6UQtdMLDa+KKa+nOumMNFZ2X5n6JsZ6QLH4o3DPvZAvngfOhlnHzknJWvKhZAiHTve7OhbdPlDqOBdjZdH8+YlBlO1Nn5s5LwKrp0qaGynYL9oMHbgMRwJaAh218bLYx7RVqWXqcnNZ2ZXAVdtOHZpsP8m4E31XG2dJi36YSiRmR0EBqFEIa8FrjSzuwHM7IEG5XEcx5mKOiz6RdKkoo+FEu3YUOdpwJGSrpV0o6TXxxqSdL6kGyTdsPrDhxsS13EcZxwD1tZUamsrTUbdlOnsbgGeB5wOHAp8VdJ1Gwch8jjU3QCPOX7n5E5UygsSm75dIUWClgLXRyLmflie/P7jtxSNTgithkHXvsBdYoMXL3T9jAgT6SNXcNNUrjeBhcTcJyJwotElBY+ndLqCErHvinyvSpQn4+yneb4p102Ba2bgkgldT1E3TiqQvsychlRbdWC0zI9UnSYVfZlQov1kcacPAw9L+jJwElDbaLPjOM6sdD1NcZOum2QoEfBp4MWStkh6LNkssb0NyuQ4jlOdmuIrF0VjFr2ZrUgahBItA5cPQony8kvNbK+kzwM3kzkiLjOzb5a6QJVJUrEUB8monOBQ6K4Zcd0Mjq2NHQvrFk2SSvm2hm6ayGQTABtknwwnc6wFFQZTzRO+w1DmketPkxagxS97LSTuv2wGjonHI1E3I+VrxcdGjleYkJWeVBjsD143ixwLzw+jaiJunFGX5Pi1RidZVfGT1U27B1rL0OjM2FQoUf753cC7m5TDcRxnJjpuwPQiBUIyTj55/nhbhVZ8sL+0FKkbXHdp2HsY7wWUYWDlhNOvLYyZzy31tdCcI2LdF1jscdOp7Px2J0VyTkFiMLXQYl8b/Tt2fmIwN0kiBcJwDL8oDj7CSN1YnHxE1lJG9OC8Jg1uS/eK204/FL3jOE6juKJ3HMfpNx3vyW4KRR8dVFK8X6uEu2Up4sZZDgZjw/LB8fDYUqIPHea9HkzAsOD81aALuZq7acJxsLgbx8IKQ4YDYMENjgzMziM14CYmmZs+5o6hwHUTGYyNnTN2rahgebVw+b/wJYvNORg/vdDlZ7H7Y3J5kRtqbq+mK3rHcZwe4xOmHMdx+k/XJ0xtDkUfc8MkpneHaQ1G4uQjrpkty+vly+F+Xr51JM5+3I0TumvCeN2Bm+bRteXg/HHLYjVw3oQ97NVo9sqx0xdGl5YanJVUWoKYJzEVVVPkmtHq6N+itgpdIAPXzXJQHrpx8uNhkyHDZov+x1LXbyMedeM4jtNvum6U9FfRV5k5Gx6OBDEXDaYOjodW/LbldTNqa76/bXlleGxLYN0P9kOL/uDquhk1sOS1fjorQVdjYP1b0K9cDaz3gXwjPQYCovHSiemSTpyExT7xnI37sTj4yMDriBUfvCNLq7FjQa92UB5a1JGB17UtwXsTWPeDy44EAcSs98Rgc+HjiVn8RXXngS1agNnpr6J3HMepBXXAtzQZV/SO4zgp3KJvGUmXzfg3Fl/+Lx7kuxS6XnLXzNalcXcNwGO3Hsz+bnl0/diWg8P9bfl5B4PB1h+tbAv2t47Lv7L+la0t5ykQRgZz1+UedqeLUjAM09UX+bHy8naO5XaSUknNBiQGa8PR0KVw4DV32Sw/GrgZg/3l/BVU4M6x5SAIIH8FQ9fQ6tbxuRYj70Xo+km5sWLlCRbuxikaee4I/VP0juM4ddKDOPom89E7juP0Alm5LdmOdKak2yTtk3RRpFyS3peX3yzpuXXI3w+LfqoIm8nfykjoeYHrY7C/HETlhBE2A5fNEdseGR47+jE/GO5v3/pDAA48evjw2P0/fvyYLGsjsfVBzHy+H0YCraWiZmY1TDwqZ0hTIXfR7JWpOPsRN05WIXTXbHlkfX/rw5mfZ/mR9Xd19dB1VfDoYZkrMQja2RBHrzH5LCKfFcgfpcprVZT2sklquIykZeD9wBlkq+tdL2mPmd0aVDsLODHfTgX+OP87E0mLXtIFko6c9UKO4zibnOcD+8zsTjM7CFwBnL2hztnARyzjOuAIScfMeuEyrpsnk/3yfDzvdmxuM85xnE1HBdfNdkk3BNv5QTM7gHuCz/vzY1SsU5mk68bMflfS7wE/B5wHXCLp48AHzexbswrQCiLumDLEFhYJ3SiDCVGP3/KPw2NPO/S7w/1/8/j7Abj8B0cPjz2yGkbdbBtrc2nEjTR+/elwd8wiqPS1ReoWJGAdunSWAt/L8o+DSX0PZe/j2o3rq3Zue96zhvtrWw4FYPUxYQ6EiCxFGSU7Hoo4hlElBcIBMzu5oCy1emjZOpUpNRhrWczed/NtBTgS+ISkP5pVAMdxnNZjJbfJ7Ad2Bp+PA+6dok5lkha9pF8DzgUOAJcBv2lmj0paAu4A3jGrEAtnSmtkbZiCQGPHAFbywdIfrBwyPHb7I08e7v/XlWwQNhyMDesO4uvXCtq3yPWnw634RVBpfkJi8e6R/WEKg/Vjq49Zt+kOPjF7x5Zf9Jz1Y8Fg7KDuSA766JoO8ev38XWqaeD9euBESScAfwucA7x2Q509wAWSriAbhP2+md0364XLRN1sB/6lmX0nPGhma5J+YVYBHMdxWk8Nit7MViRdAFwNLAOXm9ktkt6cl18KXAXsAvYBPyJzl89MGR/9OyeU7a1DCMdxnFZT07iDmV1FpszDY5cG+wa8pZ6rrdOPOPoqMbpDJsfi2og7RwXHs7/h8n4HV9cf6Y9WxudNh+kOBjHzqRQIPw7aDOPoB26conz2UWYet+1hv7xlWMQ1opibJHStBPuD1BgKMmiMxMQvZ++TDlsOjo2nQFgL0h6sBeUDl06RuyYm/7TZZON15zvaW3YyVJvph6J3HMdpEl94xHEcp9+4Rd82km6c8VURQneHcn/MiAskzAgZuE5WBgs4hF3sIIvgj8j6wCvBOWFGyrILjxwMMlaGba2uLuUyTY7KocC1UzZap6mXfCTqpIP/SNNk9Sw6J3o8EdWiEXfN+v7g8Crj7hZYj8aRxctTC48M98NzElE5I6TKIyzcY9jB9zOkf4recRynTtxH32KmTKI0sHLDLza0mEcGyPLjA8saYD3bPKzmbYUDqLMuDh7KMihfG+l9jMtdODcgmjc83pbTAEWDlanBzMFgaNhWuNTfcE2F9WPhYGtti4NHegFhW7HY/rFrxYicv3A6/r/QX0XvOI5TE+r4wiOej95xHKfnbA6LPrLUXlGu7+EK9SO/4Ou/hxb8tK9GfidH0hUs5a6bpSCHfMLZN3L+0PUy7q7J9vOp6mvj5+TC5n+DC3S8C9paEu6WwVcx4o0oGmwdvC4WOVbASNXB1x6cM2KRpt6ByPlR10yBa2coayKFQqtcMyk6/n+zORS94zjOtPhgrOM4zibAFX37iebKLoguia9gH8bRr5+3Hq8cpiUYr1u0FGGKQVNh7P6IWIOom0AmC+uuDdopip0fd+30Lpd4ixmJnY+5fEJ3SaKt8PRh1E34vVYZjYt4/EbcTEujf4FRWSMpElJx9sk4/KK686Lj/xebQtE7juNMi/Com4mkVjwP6p0iaVXSa5qUx3EcpzIllxFssx+/MYu+5Irng3rvIsvRPB2xdAVV5tdHJgxZiV/wteHkqvXK4ZK6tjR5qb9YD3TUYxRzrYRumPxv6K6J1LWRSJzIxYr6wmX7yCPr2PV8wlXB8nlJxjNvxMuJZ39UzI1T4M4YyGUpWRMrSBZmp1yacCzYL1q4JPpaTeuamZcbp+PvcpMWfZkVzwHeCnwSeKBBWRzHcaannqUEF0aTPvrYauanhhUk7QBeDbwMOKWooXwl9fMBlo88MjuYsEZGiFnsgbkjIj2C8JSiFKVDMymMs7exYlUZgQ2vG3lxRiz6tYTFHykfNZ1iFy15rIBpE5UN48zn+c+SuFbSCk41H7Oy49k0Rl/nSBw7QQ9zYN1b0f9AJKBgqkH21GBoQVqDlMUfjaNP9igWq0Xb7JYpQ5MWfcozAfAe4EIzW43UXT/JbLeZnWxmJy8fflhtAjqO45TCLfpCyqxmfjJwRW7xbgd2SVoxs081KJfjOE55rPtRN00q+uSK52Z2wmBf0oeAz86s5GNx8ImBqpGqsZGuMHd9JJPlyKDXSBdXY8dGiPUHIyNRI93umJumwHUTH7Ubv1ZRt37m7mqRb2ICSddPYcrFeii857qeRYnB2KgbI+x7R97x+FyRhEipR5ly3STi4IsGYyceC4+nysu0VRdzsNYlPRH4GPAU4NvAL5nZ9zbU2Ql8BHgymUNvt5m9N9V2Y64bM1sBBiue7wU+PljxfLDqueM4TheYU3jlRcA1ZnYicE3+eSMrwH8ws2cCLwDeIumnUg03OmEqteL5huNvaFIWx3GcqZmP//1s4LR8/8PAtcCFI2KY3Qfcl+//g6S9ZIEvI2HrG+nvzNiIi2KkIHZs5Hg6rGfYRS7qdg+7oGX668XE3TFlyiNRN6kufhV3SMIdEY0wWfSAVRU31RSypl1P8eaj6QpC10zgI054F6NyR72EEfFGyiu4bpLlYbuxqJxFpziYRLWB1u2Sbgg+7zaz3SXPPTpX5JjZfZKOmlRZ0lOAfw58LdVwfxW94zhODYhKbpkDZnZyYVvSF8j86xv5nUoySYeTzT/6dTP7Qap+LxR9clApYpyPnBMuoryWFWipyFwK2xofjbWIGRu2XzRLdhKFi3gnB1tj51Qwl4YWf6L9kSYTX0aFAdpWLR4+Yxx69FEVNK9I5ZHc8rE4+dRznTL+f/2ik49Vi4kvd074pc86eXtW6nr/zOzlhdeQ7pd0TG7NH0PBJFJJW8mU/EfN7Moy1/UVphzHcVLMJ45+D3Buvn8u8OmNFZTFon8Q2Gtm/71sw67oHcdxUsxH0V8MnCHpDrIcYRcDSDpW0g0Zx7AAAA9WSURBVCCo5YXA64CXSbop33alGu6F62aE1Fz6yGDraK83jy1fiwfHj+YNt/GD0S50gWunNCnXTcGxsjH5Qb2Fu0g2KzE3T0E8QMxjF60765SDlLukijslli6hguunjFyNUU/oZPoyZg8Cp0eO3wvsyvf/kimeQP8UveM4Tt103AByRe84jpPAUyAsAFmJrmjyFzgeprDe6y1IgTDSQiROfSQ3e+xaNZKIqommNoi4dhRz58wqhxOdR1A4OyNSUPRaJaNaUiE+ZakSoBWdP1JQnrpWfrOV3DkN03WXZicVveM4ztyoZ6B1obiidxzHSeGKvmcMAmlGDqXdOAMUDY8IK8w4rz5WXNRkwrWTFCWSHdMpIDrLKVYvXjyNJ3LGzBpTUzbTZeFrmypvGRVnxrYSV/SO4zgJtNZtTd8PRR8Z/yyywuNEApJjY6ob240EPBcu7zasMH6wKF996eXfCi36yLWmiL0vda2y53eV2Ghq2XOIvipVjP8oybQHdVI2LQIV4uxjc1IKyyPnzwv30TuO4/Qfd904juP0HVf0HaXwi4umt9xYuqGpMjkJJ4gy60tU0JctnWM9mR2z6LqJ8p4Tza6ZSKVamJEz4oaJ5vOfJykvXmKwtbDu0DVTwwuUcBPVhVv0juM4fccVveM4To8xT4HQDWLhCcmIhfLumOKa8+lwp+Phi47H+tgVzl80yRVnZm0/2I+lFYhcftQdMzkEq0oYfIW4sdqolLYgdV5BhE30nJYFbnkcveM4zmZg5oG0xeKK3nEcJ4Fb9G0j5plJ9YWrlMfCJgpegrn1QGedxFQphUL59rv0z5Far2Zmhu9KQShN0g00uflGxK4x6ibprkm0U8kzV/eX6BOmHMdx+o8PxnaNmLk09a/1zEnx62fmFAZTnFOFRVtG06YNSCWUjxnsqUxkiUD55ABslaxoU1DNio4di99B2aRola7VMK7oHcdx+ozR+cHYpXQVx3GczY2s3DbTNaQnSvpzSXfkf4+cUHdZ0l9J+myZtjeFoo9+Eab1LSQ28GIljkfLNf+tinxl76+IyDXreukXRdGjTCLGXArJtlIPSwVbpP2mX6VCWYbbZK0XbTdyT5We3zyJ/S+V/f8qz0XANWZ2InBN/rmItwF7yza8KRS94zjOtBT9jtVt0QNnAx/O9z8M/GJUHuk44OeBy8o27D56x3GcSZhVWXhku6Qbgs+7zWx3yXOPNrP7skvafZKOKqj3HuAdwOPKCtVfRV8QkRALZ07GxqciNWaN5JiWaSyIadMlDMunFHoRrpxpv5dZ2kzExo9ULUyXEGHhqSwDErJOlTph0feUovz7c8DMTi4qlPQF4MmRot8p07ikXwAeMLMbJZ1WVqj+KnrHcZyaqGvMycxeXngN6X5Jx+TW/DHAA5FqLwReJWkXcAjweEl/amb/etJ13Ue/kbKDmlaiPDmYOeM2TZupc6Llk0fFOjsAO+tgYOL8Ku1XGrhtwDlc2gk9ITa+8rMqeD7JweDY8SYxYM3KbbOxBzg33z8X+PSYKGa/ZWbHmdlTgHOAL6aUPLiidxzHSTOfqJuLgTMk3QGckX9G0rGSrpql4UYVvaQzJd0maZ+ksVAhSb8q6eZ8+4qkk5qUx3EcZxrmEXVjZg+a2elmdmL+96H8+L1mtitS/1oz+4UybTfmo5e0DLyf7JdpP3C9pD1mdmtQ7S7gZ83se5LOAnYDp9YuTGRgNZnKPDmXnfgveLStyeI1zrTXT4yqTZ0Hv4Mkk54NHlWqfKTRydcqIplaoQFmTodQpbyFVIi6aSVNWvTPB/aZ2Z1mdhC4gixOdIiZfcXMvpd/vA44rkF5HMdxqjPruFkLaFLR7wDuCT7vz48V8Ubgc7ECSedLukHSDas/fLhGER3HcSaTTZiyUltbaTK8snRnVdJLyRT9i2Ll+YSD3QCPOX5n7U+zkhun6MRh3dTFykpVgVmfSIV+eeciawaUda0k7i+dnTKsXFKmsYuUl2Gh1OjOSWa0XDSevbKQ/cDO4PNxwL0bK0l6NtlU3rPM7MEG5XEcx5mKNlvrZWjSdXM9cKKkEyRtI4v53BNWkHQ8cCXwOjO7vUFZHMdxpqMHPvrGLHozW5F0AXA1sAxcbma3SHpzXn4p8E7gJ4APSAJYmTR9OCTpbikUbNDA5HZLtZmc1z7h+otiin5/JXdNjfe3EDdRBddLJTdOyTYLz6vaRh1M6zoped7ULqi5u3Qq5bppJY2mQDCzq4CrNhy7NNh/E/CmJmVwHMeZmY67bjZvrptIbH1IYeh80tqaJgh6RmocnZtKvG7/DxQzhXWffH5TDsCWamPeVJCj1jj8eWP4UoKO4zi9xy16x3GcntNtPe+KHki6cUKiofOVBoPb0y+d2YvUxZe/iRz1YZMVplyMVqhykQp1Z2HGV7XT7poNaK3bvhtX9I7jOJMwfMKU4zhOnxHtTm9QBlf0G0nE2ceosgrcPKktwKcL73ilCRDTtB85VuG5TDPlopIsC2aqx94l144resdxnJ7jir6nVMkrnqCzicAWIPdcn9WsA7Mznr+IKRfTUmtHqWxblSz+Bh+W++gdx3H6j0fdOI7j9BrrvOumF4uD17VmYyk6mr0uSp/uZVbEbAN+imwzYmrPVtszWfSg6jQYmaIvs82ApCdK+nNJd+R/jyyod4SkT0j6G0l7Jf2LVNu9UPSO4ziNslZym42LgGvM7ETgmvxzjPcCnzezZwAnAXtTDbuidxzHSTCnpQTPBj6c738Y+MUxOaTHAy8BPghgZgfN7O9TDbuPvm4WkYVws7lZmqTOFAltySvfFE25Ydro3imvxLdLuiH4vDtfCrUMR5vZfdnl7D5JR0XqPBX4O+BPJJ0E3Ai8zcwmLqbtit5xHGcSZrBa2i9zYNLiSZK+ADw5UvQ7JdvfAjwXeKuZfU3Se8lcPL+XOslxHMeZRE1RN2b28qIySfdLOia35o8BHohU2w/sN7Ov5Z8/QbEvf0jvfPRzi76pgyprUfY8Oqa2763OEKymI0aKolIWHakyL5nacK9lmUPUDdma2ufm++cCnx4Xw74L3CPp6fmh04FbUw27Re84jjMJA+azZuzFwMclvRG4G/hXAJKOBS4zs115vbcCH5W0DbgTOC/VsCt6Z2F0pucVUsdSgHVetys0NnA7p8kz1vzMWDN7kMxC33j8XmBX8PkmoHAcIIYresdxnEkYVQZjW4kresdxnBQdT4HQW0Uf9uhatHrfpme+2Snn+BL0PWa+Cn38f3NF7ziO02e6n9TMFb3jOM4kDPA0xe0n5i5wd858WXiETdNLDSavX2NbTTzLLv0/LOJlcovecRynz1RKgdBKNq2iLzIK3NKvzsKt9Sr0YZS+o2LPxCJfMgObQxx9k2xaRe84jlOa+cyMbQxX9I7jOCncR98vpukhdtUDMKBTrpc6Sd1417/YLtHml9DMo24cx3F6j1v0juM4fcaw1dVFCzETruhroM29TmcGmv5iu+4a2iwv/vzSFDeGK3rHcZwUHQ+vbHSFKUlnSrpN0j5JY8tdKeN9efnNkp7bpDyO4zhVMcDWrNTWVhqz6CUtA+8HziBb5/B6SXvMLFz26izgxHw7Ffjj/K/j9J/N4vroOjafhUeapEmL/vnAPjO708wOAlcAZ2+oczbwEcu4DjgiXxTXcRynNdjqaqmtrTTpo98B3BN83s+4tR6rswO4L6wk6Xzg/Pzjj+9629u/Wa+orWA7cGDRQtRMH+8J+nlffbwngKenq0zmH/je1V+wT2wvWb2Vz7BJRR8LKdjYVy1TBzPbDewGkHSDmVVaL7EL9PG++nhP0M/76uM9QXZfs7ZhZmfWIcsiadJ1sx/YGXw+Drh3ijqO4zjODDSp6K8HTpR0gqRtwDnAng119gCvz6NvXgB838zu29iQ4ziOMz2NuW7MbEXSBcDVwDJwuZndIunNefmlwFXALmAf8CPgvBJN725I5EXTx/vq4z1BP++rj/cE/b2vSsg6nsPBcRzHmUyjE6Ycx3GcxeOK3nEcp+e0VtH3MX1CiXv61fxebpb0FUknLULOqqTuK6h3iqRVSa+Zp3zTUua+JJ0m6SZJt0j60rxlrEqJd/AJkj4j6Rv5PZUZN1soki6X9ICk6PyaLuqK2jGz1m1kg7ffAp4KbAO+AfzUhjq7gM+RxeK/APjaouWu4Z5+Bjgy3z+r7fdU9r6Cel8kG4B/zaLlrun7OgK4FTg+/3zUouWu4Z5+G3hXvv8k4CFg26JlT9zXS4DnAt8sKO+Urmhia6tF38f0Ccl7MrOvmNn38o/Xkc0raDtlviuAtwKfBB6Yp3AzUOa+XgtcaWZ3A5hZ2++tzD0Z8DhJAg4nU/Qr8xWzGmb2ZTI5i+iarqidtir6otQIVeu0iaryvpHMCmk7yfuStAN4NXDpHOWalTLf19OAIyVdK+lGSa+fm3TTUeaeLgGeSTZx8a+Bt5l1PKNX93RF7bQ1H31t6RNaRGl5Jb2UTNG/qFGJ6qHMfb0HuNDMVjNDsROUua8twPOA04FDga9Kus7Mbm9auCkpc0+vAG4CXgb8E+DPJf2Fmf2gaeEapGu6onbaquj7mD6hlLySng1cBpxlZg/OSbZZKHNfJwNX5Ep+O7BL0oqZfWo+Ik5F2XfwgJk9DDws6cvASUBbFX2ZezoPuNgy5/Y+SXcBzwC+Ph8RG6FruqJ22uq66WP6hOQ9SToeuBJ4XYutwo0k78vMTjCzp5jZU4BPAP++5Uoeyr2DnwZeLGmLpMeSZWfdO2c5q1Dmnu4m66Eg6Wiy7I93zlXK+umarqidVlr01lz6hIVR8p7eCfwE8IHc+l2xlmcULHlfnaPMfZnZXkmfB24G1oDLzKy1KbRLfld/AHxI0l+TuTwuNLNWpt4dIOnPgNOA7ZL2A78PbIVu6oom8BQIjuM4PaetrhvHcRynJlzRO47j9BxX9I7jOD3HFb3jOE7PcUXvOI7Tc1zRO47j9BxX9I7jOD3HFb3TWfL89jdLOkTSYXn+9GctWi7HaRs+YcrpNJL+M3AIWVKx/Wb2hwsWyXFahyt6p9PkOVuuB/4R+BkzW12wSI7TOtx143SdJ5ItkPE4MsvecZwNuEXvdBpJe8hWSjoBOMbMLliwSI7TOlqZvdJxypCv6LRiZv9L0jLwFUkvM7MvLlo2x2kTbtE7juP0HPfRO47j9BxX9I7jOD3HFb3jOE7PcUXvOI7Tc1zRO47j9BxX9I7jOD3HFb3jOE7P+f/18cLWUNaPeQAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "# look at the results!\n",
        "plt.colorbar(mesh.plot_image(phi)[0])\n",
        "plt.title('Electric Potential, $\\phi$');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "iooxa": null
      },
      "source": [
        "## What just happened!?\n",
        "\n",
        "In the notebooks that follow, we will:\n",
        "\n",
        "- define where variables live on the mesh ([mesh.ipynb](mesh.ipynb))\n",
        "- define the discrete divergence ([divergence.ipynb](divergence.ipynb))\n",
        "- use the weak formulation to define a solveable system of equations ([weak-formulation.ipynb](weak-formulation.ipynb))\n",
        "- solve and play with the DC resistivity equations ([all-together-now.ipynb](all-together-now.ipynb))"
      ]
    }
  ],
  "metadata": {
    "interpreter": {
      "hash": "28579dcce9b4c93863c7d16f93a67ba69218f1244c98cdb9f71bbf270d0a38fb"
    },
    "iooxa": null,
    "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.9.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}