{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Weak Formulation\n", "\n", "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)!\n", "\n", "\n", "\n", "
\n", "**Figure 5**: Discretization using the weak formulation and inner products.\n", "
\n", "\n", "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):\n", "\n", "$$\n", "\\frac{1}{4}\\sum_{i=1}^{4} \\mathbf{P}_i^\\top \\sqrt{v} \\boldsymbol{\\Sigma}^{-1} \\sqrt{v} \\mathbf{P}_i \\mathbf{j} = \\mathbf{D}^\\top v \\phi\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "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 }, "outputs": [], "source": [ "%matplotlib inline\n", "from SimPEG import Mesh\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEw9JREFUeJzt3W2MXOV5xvHrchFKeCkuaeNUNNiW80KhSbbEMVhV5eUD\nBaK2bhvCS01h06hBieBDpEjQKsQ4pFIoKKmAIEKEkhg5Ii6GxCRAcYsPyNRsTIgTQ6BxqNcBAm4l\nWAkH1Fr47ocZm81yZl/OPnueM2f+P2nkOTPHu3eePN6Lue85s44IAQAwWwtyFwAA6E8ECACgEgIE\nAFAJAQIAqIQAAQBUQoAAACrJHiC2b7O9z/ZPejy/yva47ce7t8/WXSMA4M2OyF2ApK9LulHS+inO\neTgi/rymegAAM5D9FUhEbJP08jSnuY5aAAAzlz1AZmil7Z22v2/75NzFAACa0cKazg8lnRgRr9o+\nR9J3JL0nc00AMPAaHyARsX/C/fts32z7+Ih4afK5tvlgLwCYpYioNCZoSgvL6jHnsL1owv0VklwW\nHodEBLcEt7Vr12avoU031pP1bOptLrK/ArH9LUnDkt5m+xeS1ko6UlJExK2SzrX9SUkHJL0m6fxc\ntQ6SsbGx3CW0CuuZFuvZDNkDJCL+eprnvyLpKzWVAwCYoaa0sNAwIyMjuUtoFdYzLdazGTzXHliT\n2I42/e8BgPlmW9HnQ3Q0TFEUuUtoFdYzLdazGQgQAEAltLAAYIDRwgIA1I4AQSl6zGmxnmmxns1A\ngAAAKmEGAgADjBkIAKB2BAhK0WNOi/VMi/VsBgIEAFAJMxAAGGDMQAAAtSNAUIoec1qsZ1qsZzMQ\nIACASpiBAMAAYwYCAKgdAYJS9JjTYj3TYj2bgQABAFTCDAQABhgzEABA7QgQlKLHnBbrmRbr2QwE\nCACgEmYgADDAmIEAAGpHgKAUPea0WM+0WM9mIEAAAJUwAwGAAcYMBABQOwIEpegxp8V6psV6NgMB\nAgCohBkIAAwwZiAAgNoRIChFjzkt1jMt1rMZCBAAQCXMQABggDEDAQDUjgBBKXrMabGeabGezUCA\nAAAqYQYCAAOMGQgAoHYECErRY06L9UyL9WwGAgQAUAkzEAAYYMxAAAC1I0BQih5zWqxnWqxnM2QP\nENu32d5n+ydTnHOD7d22d9oeqrM+AEC57AEi6euSzur1pO1zJC2LiHdLulTSLXUVNsiGh4dzl9Aa\nEaH7H7pfzOfSYX82Q/YAiYhtkl6e4pTVktZ3zx2VdJztRXXUBqSw6Z5NuvnBm3XX9+7KXQqQVPYA\nmYETJD074fj57mOYR/SY04gIXX/79Xpl6Su6bv11vApJhP3ZDEfkLiC1kZERLVmyRJK0cOFCDQ0N\nHX65e2jTcTz98RlnSFLnWBru/slxteNHJRUalbRgwUMNqKcdx1u3do6b8O+ln44P3R8bG9NcNeI6\nENuLJd0TEe8vee4WSVsj4tvd46clrYqIfSXnch1IIrbEUs5NRGjleSs1esqoZEkhnfbkadq+cbvs\nSm+7Rxf7M502XAfi7q3MZkkXS5Lt0yWNl4UH0DSb7tmkXcfuemNnW9p1zC5mIWiN7C0s299S5zXp\n22z/QtJaSUdKioi4NSLutf1h2z+X9CtJH8tX7SAp9EbrAFU88tgjWv76cnmPNf7iuBa+Y6EiQtt2\nbNNH/uwjucvrc4XYn/k1ooWVCi2sdOxCEcO5y2iNoih462lC7M905tLCIkBQih4zmoz9mU4bZiAA\ngD5DgKCHIncBrcJ1C6kVuQuACBAAQEXMQFCKHjOajP2ZDjMQAEDtCBD0UOQuoFWYgaRW5C4AIkAA\nABUxA0EpesxoMvZnOsxAAAC1I0DQQ5G7gFZhBpJakbsAiAABAFTEDASl6DGjydif6TADAQDUjgBB\nD0XuAlqFGUhqRe4CIAIEAFARMxCUoseMJmN/psMMBABQOwIEPRS5C2gVZiCpFbkLgAgQAEBFzEBQ\nih4zmoz9mQ4zEABA7QgQ9FDkLqBVmIGkVuQuACJAAAAVMQNBKXrMaDL2ZzrMQAAAtSNA0EORu4BW\nYQaSWpG7AIgAAQBUxAwEpegxo8nYn+kwAwEA1I4AQQ9F7gJahRlIakXuAiACBABQETMQlKLHjCZj\nf6bDDAQAUDsCBD0UuQtoFWYgqRW5C4AIEABARcxAUIoeM5qM/ZkOMxAAQO0IEPRQ5C6gVZiBpFbk\nLgAiQAAAFTEDQSl6zGgy9mc6zEAAALUjQNBDkbuAVmEGklqRuwCIAAEAVMQMBKXoMaPJ2J/pMAMB\nANSOAEEPRe4CWoUZSGpF7gKgBgSI7bNtP237Z7avKHl+le1x2493b5/NUScA4NdNGSC2f9P2spLH\n35/im9teIOkmSWdJOkXShbZPKjn14Yg4tXv7QorvjTf79Oc+rTs336nOHGk4dzmtsHfPHq276CJt\nXbdO6y66SHv37MldUt9ifzZPzwCxfZ6kpyVtsv2k7Q9NePobib7/Ckm7I2JvRByQdIek1WXlJPp+\nmMLjex/XJXdfopUfXSlJ4g0Jc7N3zx7deOaZ+syGDVpXFPrMhg268cwzCZGK2J/NM9UrkH+Q9MGI\nGJL0MUm32/7L7nOpfqCfIOnZCcfPdR+bbKXtnba/b/vkRN8bk9jWq0te1egfjEoqtPKjKyf8Fx9m\n6xtXXaV1zzyjo9Xp2B8tad0zz+gbV12Vt7A+xf5sniOmeO43IuIFSYqIH9g+Q9L3bL9TUp3/j/1Q\n0okR8artcyR9R9J7ep08MjKiJUuWSJIWLlyooaEhDQ8PS3pjkMlx+fH4i+OSJC3t/DF6zKguvupi\nPfajx/TFtV/MXl+/Hf/XE09oh95othTdPw/+8peNqK/fjtmfaY4P3R8bG9OcRUTpTdJ/SFo26bFj\nJf27pP/t9fdmc5N0uqT7JxxfKemKaf7OHknH93guUN2qS1aFrlZorUKKOO3c0+LOzXfGwYMHc5fW\nl65esyb2dy5XOHzbL8XVa9bkLq0vsT/nR/fnZqWf4VO9AvmkpAW2T46In3Z/Or9i+2xJF8w9uiRJ\nOyS9y/ZiSS90v+6FE0+wvSgi9nXvr1Dn4seXEn1/TBAROmrsKL1v//s0Kmn7xu2yGT9VNXLNNVr7\n6KOH21i/krR22TJdfs01uUvrS+zP5uk5A4mIH0fEbkkbbV/hjrdK+pKkT6X45hHxuqTLJD0g6UlJ\nd0TEU7Yvtf2J7mnn2n7C9o8k/bOk81N8b7zZqYtP1fq/Wq/tG7dLKvjHOUeLly7V5Vu26Po1a3TJ\n0JCuX7NGl2/ZosVLl+YurS+xP5tn2o8ysX20pGslfVCdFtYGSddGxMH5L292+CiTdOxCEcO5y2iN\noigO96Ixd+zPdOb7o0wOSHpN0lslvUXSniaGB1Ibzl1AqxAeqQ3nLgCaWYDsUCdAPiTpj9W52O9f\n5rUqAEDjzSRAPh4Rn4uIAxHxQkSslrR5vgtDbkXuAlpl4lsokUKRuwBoBgESEY+VPHb7/JQDAOgX\n/D4QlOL3LaDJ2J/p8PtAAAC1I0DQQ5G7gFZhBpJakbsAiAABAFTEDASl6DGjydif6TADAQDUjgBB\nD0XuAlqFGUhqRe4CIAIEAFARMxCUoseMJmN/psMMBABQOwIEPRS5C2gVZiCpFbkLgAgQAEBFzEBQ\nih4zmoz9mQ4zEABA7QgQ9FDkLqBVmIGkVuQuACJAAAAVMQNBKXrMaDL2ZzrMQAAAtSNA0EORu4BW\nYQaSWpG7AIgAAQBUxAwEpegxo8nYn+kwAwEA1I4AQQ9F7gJahRlIakXuAiACBABQETMQlKLHjCZj\nf6bDDAQAUDsCBD0UuQtoFWYgqRW5C4AIEABARcxAUIoeM5qM/ZkOMxAAQO0IEPRQ5C6gVZiBpFbk\nLgAiQAAAFTEDQSl6zGgy9mc6zEAAALUjQNBDkbuAVmEGklqRuwCIAAEAVMQMBKXoMaPJ2J/pMAMB\nANSOAEEPRe4CWoUZSGpF7gIgAgQAUBEzEJSix4wmY3+mwwwEAFA7AgQ9FLkLaBVmIKkVuQuAGhAg\nts+2/bTtn9m+osc5N9jebXun7aG6awQAvFnWALG9QNJNks6SdIqkC22fNOmccyQti4h3S7pU0i21\nFzqQhnMX0BoRofsful/M51Iazl0AlP8VyApJuyNib0QckHSHpNWTzlktab0kRcSopONsL6q3TKC6\nTfds0s0P3qy7vndX7lKApHIHyAmSnp1w/Fz3sanOeb7kHCRX5C6gFSJC199+vV5Z+oquW38dr0KS\nKXIXAElH5C4gtZGRES1ZskSStHDhQg0NDWl4eFjSG4NMjmd2bHeO32gXcFzt+FFJhUYlLVjwUAPq\nacdx7n8f/Xp86P7Y2JjmKut1ILZPl3R1RJzdPb5SUkTEtRPOuUXS1oj4dvf4aUmrImJfydfjOhA0\nRkRo5XkrNXrKqGRJIZ325GnavnG77EpvuweS6+frQHZIepftxbaPlHSBpM2Tztks6WLpcOCMl4UH\n0DSb7tmkXcfu6oSHJFnadcwuZiFojawBEhGvS7pM0gOSnpR0R0Q8ZftS25/onnOvpD22fy7pq5I+\nla3gAcJ1C3P3yGOPaPnry7Vqzyp9YPsHtGrPKi0/uFzbdmzLXVrfY382Q/YZSETcL+m9kx776qTj\ny2otCkjgy5//8uH7RVEc7kUDbcFnYQHAAOvnGQgAoE8RIChFjzkt1jMt1rMZCBAAQCXMQABggDED\nAQDUjgBBKXrMabGeabGezUCAAAAqYQYCAAOMGQgAoHYECErRY06L9UyL9WwGAgQAUAkzEAAYYMxA\nAAC1I0BQih5zWqxnWqxnMxAgAIBKmIEAwABjBgIAqB0BglL0mNNiPdNiPZuBAAEAVMIMBAAGGDMQ\nAEDtCBCUosecFuuZFuvZDAQIAKASZiAAMMCYgQAAakeAoBQ95rRYz7RYz2YgQAAAlTADAYABxgwE\nAFA7AgSl6DGnxXqmxXo2AwECAKiEGQgADDBmIACA2hEgKEWPOS3WMy3WsxkIEABAJcxAAGCAMQMB\nANSOAEEpesxpsZ5psZ7NQIAAACphBgIAA4wZCACgdgQIStFjTov1TIv1bAYCBABQCTMQABhgzEAA\nALXLFiC2f8v2A7b/0/a/2j6ux3ljtn9s+0e2f1B3nYOKHnNarGdarGcz5HwFcqWkf4uI90p6UNLf\n9zjvoKThiPjDiFhRW3UDbufOnblLaBXWMy3WsxlyBshqSd/s3v+mpL/ocZ5Fq6124+PjuUtoFdYz\nLdazGXL+YH57ROyTpIh4UdLbe5wXkrbY3mH772qrDgAwpSPm84vb3iJp0cSH1AmEz5ac3uvtU38U\nES/Y/h11guSpiNiWuFRMMjY2lruEVmE902I9myHb23htP6XObGOf7XdI2hoRvz/N31kr6ZWI+FKP\n53kPLwDMUtW38c7rK5BpbJY0IulaSZdI+u7kE2wfJWlBROy3fbSkP5G0rtcXrLoIAIDZy/kK5HhJ\nGyW9U9JeSedFxLjt35X0tYj4U9tLJd2tTnvrCEkbIuKLWQoGAPyaVl2JDgCoT9++PZYLEdOwfbbt\np23/zPYVPc65wfZu2zttD9VdYz+Zbj1tr7I9bvvx7q3sDSWQZPs22/ts/2SKc9ibMzTdelbZm30b\nIOJCxDmzvUDSTZLOknSKpAttnzTpnHMkLYuId0u6VNIttRfaJ2aynl0PR8Sp3dsXai2yv3xdnbUs\nxd6ctSnXs2tWe7OfA4QLEeduhaTdEbE3Ig5IukOddZ1otaT1khQRo5KOs71IKDOT9ZQ6exLT6L5d\n/+UpTmFvzsIM1lOa5d7s5x+sXIg4dydIenbC8XPdx6Y65/mSc9Axk/WUpJXdlsv3bZ9cT2mtxN5M\nb1Z7M+fbeKfFhYhooR9KOjEiXu22YL4j6T2ZawKkCnuz0QESEWf2eq47DFo04ULE/+7xNV7o/vk/\ntu9Wp81AgHQ8L+nECce/131s8jnvnOYcdEy7nhGxf8L9+2zfbPv4iHipphrbhL2ZUJW92c8trEMX\nIkpTXIho+5ju/UMXIj5RV4F9YIekd9lebPtISReos64TbZZ0sSTZPl3S+KHWId5k2vWc2KO3vUKd\nt9ITHr1Zvfvy7M3Z67meVfZmo1+BTONaSRtt/626FyJK0sQLEdVpf93d/YiTQxciPpCr4KaJiNdt\nXybpAXX+Y+K2iHjK9qWdp+PWiLjX9odt/1zSryR9LGfNTTaT9ZR0ru1PSjog6TVJ5+eruNlsf0vS\nsKS32f6FpLWSjhR7s5Lp1lMV9iYXEgIAKunnFhYAICMCBABQCQECAKiEAAEAVEKAAAAqIUAAAJUQ\nIEANbN9n+2Xbky/UBPoWAQLU458kXZS7CCAlAgRIyPby7i8wO9L20bafsH1yRGyVtH/aLwD0kX7+\nKBOgcSLiMdvflfSPkt4q6faI+GnmsoB5QYAA6V2jzgcrvibp8sy1APOGFhaQ3m9LOkbSsZLekrkW\nYN4QIEB6t6jzS882qDM8P2SqjyYH+g4tLCAh238j6f8i4g7bCyQ9YntY0uclvVfSMd2P0v54RGzJ\nWCowZ3ycOwCgElpYAIBKCBAAQCUECACgEgIEAFAJAQIAqIQAAQBUQoAAACohQAAAlfw/RamZT0QY\nxqUAAAAASUVORK5CYII=\n", "text/plain": [ "