{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from SimPEG import Mesh, Problem, Utils, np, sp, Tests\n",
      "from simpegPF.Magnetics import BaseMag, MagneticsDiffSecondary\n",
      "from scipy.constants import mu_0\n",
      "from simpegPF.MagAnalytics import spheremodel, CongruousMagBC\n",
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "# hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n",
      "# hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n",
      "# hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3))\n",
      "\n",
      "hxind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n",
      "hyind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n",
      "hzind = ((5,25,1.3),(20, 25.),(5,25,1.3))\n",
      "\n",
      "hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)\n",
      "mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2])\n",
      "\n",
      "chibkg = 0.001\n",
      "chiblk = 0.1\n",
      "chi = np.ones(mesh.nC)*chibkg\n",
      "sph_ind = spheremodel(mesh, 0., 0., 0., 100)\n",
      "chi[sph_ind] = chiblk\n",
      "model = BaseMag.BaseMagModel(mesh)\n",
      "mu = (1.+chi)*mu_0\n",
      "\n",
      "data = BaseMag.BaseMagData()\n",
      "data.setBackgroundField(x=1., y=1., z=0.)\n",
      "xr = np.linspace(-300, 300, 41)\n",
      "yr = np.linspace(-300, 300, 41)\n",
      "X, Y = np.meshgrid(xr, yr)\n",
      "Z = np.ones((xr.size, yr.size))*150\n",
      "rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]\n",
      "data.rxLoc = rxLoc\n",
      "\n",
      "prob = MagneticsDiffSecondary(mesh, model)\n",
      "\n",
      "prob.pair(data)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fields = prob.fields(chi)\n",
      "u = fields['u']\n",
      "B = fields['B']\n",
      "# mesh.plotSlice(B, 'F', view='vec', showIt=True)\n",
      "dpred = data.dpred(chi, fields)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## ?? they are different ?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "MfMui = mesh.getFaceInnerProduct(mu).diagonal()\n",
      "MfMui_test = mesh.getFaceMass(mu).diagonal()\n",
      "print np.linalg.norm(MfMui - MfMui_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9.55511263546\n"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(MfMui/MfMui_test)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 53,
       "text": [
        "[<matplotlib.lines.Line2D at 0x3363350>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD9CAYAAABTJWtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHaxJREFUeJzt3X9Qm/XhB/B3KNz1MC1Kh0AJGwgUCIEkLTTnpBoEiv0h\nSquutiJ2jMsxnXbXr6ec8wbzhmXTm9T6R9fVn93Rzp1r667FwklmWwUOC6srbUcrSPhpqeU3LE3y\n+f7Rbz9fIxioJj6svl93uSPP50Oe95PS502e5wlRCSEEiIiIAPgpHYCIiOYOlgIREUksBSIiklgK\nREQksRSIiEhiKRARkeSxFCYnJ2EymWAwGKDValFSUjJlzp///Gfo9XqkpKTgtttuw8mTJ+VYdXU1\nEhISEBcXh4qKCu+nJyIir1LN9D6F8fFxBAYGwuFwID09HS+88ALS09Pl+EcffQStVougoCBUV1ej\ntLQU9fX1cDqdiI+PR21tLSIiIpCWloaqqiokJib6fKOIiOibmfHwUWBgIADAbrfD6XQiODjYbfzW\nW29FUFAQAMBkMqGrqwsA0NjYiNjYWERFRSEgIAAbNmzAgQMHvJ2fiIi8aMZScLlcMBgMCA0NRUZG\nBrRa7dfO3b17N1avXg0A6O7uRmRkpBzTaDTo7u72QmQiIvIV/5km+Pn5oaWlBUNDQ8jJyYHVaoXZ\nbJ4yr66uDq+++iqOHz8OAFCpVLMKMNt5RETkzhd/pWjWVx8FBQVhzZo1aGpqmjJ28uRJFBUV4eDB\ng7jpppsAABEREbDZbHKOzWaDRqOZ9rGFEHPq9utf/1rxDP8NmeZqLmZipu9DLl/xWAoDAwMYHBwE\nAExMTKCmpgZGo9FtTmdnJ9atW4c9e/YgNjZWLk9NTUVbWxs6Ojpgt9uxb98+5Obm+mATiIjIWzwe\nPurt7UVBQQFcLhdcLhfy8/ORmZmJnTt3AgAsFgt+85vf4NKlSyguLgYABAQEoLGxEf7+/tixYwdy\ncnLgdDpRWFjIK4+IiOY6obA5EGGKuro6pSNMMRczCTE3czHT7DDT7M3FXL7ad874PgVfU6lUPj0+\nRkR0PfLVvpN/5oKIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYC\nERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUkeS2FychImkwkGgwFa\nrRYlJSVT5pw5cwa33nor5s+fjxdffNFtLCoqCikpKTAajVi+fLl3kxMRkdf5exqcP38+6urqEBgY\nCIfDgfT0dBw7dgzp6elyzqJFi/Dyyy9j//79U75fpVLBarUiODjY+8mJiMjrZjx8FBgYCACw2+1w\nOp1TdvAhISFITU1FQEDAtN/viw+WJiIi3/D4SgEAXC4Xli5divPnz6O4uBharXbWD65SqZCVlYV5\n8+bBYrGgqKho2nmlpaXya7PZDLPZPOt1EBF9H1itVlitVp+vRyVm+av80NAQcnJysG3btml32mVl\nZVCr1di6datc1tvbi/DwcFy4cAHZ2dl4+eWXsWLFCvcAKhVfTRARXSNf7TtnffVRUFAQ1qxZg6am\nplk/eHh4OIArh5jy8vLQ2Nh47QmJiOg747EUBgYGMDg4CACYmJhATU0NjEbjtHO/2ljj4+MYGRkB\nAIyNjeHIkSNITk72RmYiIvIRj+cUent7UVBQAJfLBZfLhfz8fGRmZmLnzp0AAIvFgr6+PqSlpWF4\neBh+fn6orKxEa2srPv/8c6xbtw4A4HA4sGnTJqxcudL3W0RERN/YrM8p+CwAzykQEV0zxc8pEBHR\n9Y+lQEREEkuBiIgklgIREUksBSIiklgKREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSWAhER\nSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJI+lMDk5CZPJBIPBAK1Wi5KSkilz\nzpw5g1tvvRXz58/Hiy++6DZWXV2NhIQExMXFoaKiwrvJiYjI61Rihk9+Hh8fR2BgIBwOB9LT0/HC\nCy8gPT1djl+4cAGfffYZ9u/fj5tuuglbt24FADidTsTHx6O2thYRERFIS0tDVVUVEhMT3QP46MOn\niYiuZ77ad854+CgwMBAAYLfb4XQ6ERwc7DYeEhKC1NRUBAQEuC1vbGxEbGwsoqKiEBAQgA0bNuDA\ngQNejE5ERN7mP9MEl8uFpUuX4vz58yguLoZWq53VA3d3dyMyMlLe12g0aGhomHZuaWmp/NpsNsNs\nNs9qHURE3xdWqxVWq9Xn65mxFPz8/NDS0oKhoSHk5OTAarXOaqetUqlmHeLLpUBERFN99RfmsrIy\nn6xn1lcfBQUFYc2aNWhqaprV/IiICNhsNnnfZrNBo9Fce0IiIvrOeCyFgYEBDA4OAgAmJiZQU1MD\no9E47dyvnvBITU1FW1sbOjo6YLfbsW/fPuTm5nopNhER+YLHw0e9vb0oKCiAy+WCy+VCfn4+MjMz\nsXPnTgCAxWJBX18f0tLSMDw8DD8/P1RWVqK1tRVqtRo7duxATk4OnE4nCgsLp1x5REREc8uMl6T6\nPAAvSSUiumaKXZJKRETfHywFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJ\nLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgKREQkeSyFyclJ\nmEwmGAwGaLValJSUTDvv8ccfR1xcHPR6PZqbm+XyqKgopKSkwGg0Yvny5d5NTkREXufvaXD+/Pmo\nq6tDYGAgHA4H0tPTcezYMaSnp8s5hw4dwrlz59DW1oaGhgYUFxejvr4ewJUPlrZarQgODvbtVhAR\nkVfMePgoMDAQAGC32+F0Oqfs4A8ePIiCggIAgMlkwuDgIPr7++W4EMKbeYmIyIc8vlIAAJfLhaVL\nl+L8+fMoLi6GVqt1G+/u7kZkZKS8r9Fo0N3djdDQUKhUKmRlZWHevHmwWCwoKiqadh2lpaXya7PZ\nDLPZ/M22hojoOmW1WmG1Wn2+nhlLwc/PDy0tLRgaGkJOTg6sVuuUnfbXvRo4duwYFi9ejAsXLiA7\nOxsJCQlYsWLFlHlfLgUiIprqq78wl5WV+WQ9s776KCgoCGvWrEFTU5Pb8oiICNhsNnm/q6sLERER\nAIDFixcDAEJCQpCXl4fGxkZvZCYiIh/xWAoDAwMYHBwEAExMTKCmpgZGo9FtTm5uLt58800AQH19\nPW688UaEhoZifHwcIyMjAICxsTEcOXIEycnJvtgGIiLyEo+Hj3p7e1FQUACXywWXy4X8/HxkZmZi\n586dAACLxYLVq1fj0KFDiI2NxQ033IDXXnsNANDX14d169YBABwOBzZt2oSVK1f6eHOIiOjbUAmF\nLw9SqVS8QomI6Br5at/JdzQTEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgKREQksRSIiEhi\nKRARkcRSICIiiaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERJLH\nUpicnITJZILBYIBWq0VJScm08x5//HHExcVBr9ejublZLq+urkZCQgLi4uJQUVHh3eREROR1Hkth\n/vz5qKurQ0tLC06ePIm6ujocO3bMbc6hQ4dw7tw5tLW14Y9//COKi4sBAE6nE4899hiqq6vR2tqK\nqqoqnD592ndbQkRE35r/TBMCAwMBAHa7HU6nE8HBwW7jBw8eREFBAQDAZDJhcHAQfX19aG9vR2xs\nLKKiogAAGzZswIEDB5CYmDhlHe+++203g4jo+rdgAWA2+3YdM5aCy+XC0qVLcf78eRQXF0Or1bqN\nd3d3IzIyUt7XaDTo7u5GT0/PlOUNDQ3TruN//qdUfr1okRmLFpmvcTOIiK5vFy9a4XBYsXq1b9cz\nYyn4+fmhpaUFQ0NDyMnJgdVqhfkrVSWE+FYhzp4t/VbfT0R0/TP/3+2KsrIyn6xl1lcfBQUFYc2a\nNWhqanJbHhERAZvNJu93dXVBo9FMWW6z2aDRaLwQmYiIfMVjKQwMDGBwcBAAMDExgZqaGhiNRrc5\nubm5ePPNNwEA9fX1uPHGGxEaGorU1FS0tbWho6MDdrsd+/btQ25uro82g4iIvMHj4aPe3l4UFBTA\n5XLB5XIhPz8fmZmZ2LlzJwDAYrFg9erVOHToEGJjY3HDDTfgtddeu/LA/v7YsWMHcnJy4HQ6UVhY\nOO1JZiIimjtU4tueEPi2AVSqb31Ogojo+8ZX+06+o5mIiCSWAhERSSwFIiKSWApERCSxFIiISGIp\nEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuB\niIgklgIREUksBSIikjyWgs1mQ0ZGBpKSkqDT6bB9+/Ypcy5duoS8vDzo9XqYTCacOnVKjkVFRSEl\nJQVGoxHLly/3fnoiIvIqlfDwyc99fX3o6+uDwWDA6Ogoli1bhv379yMxMVHOefLJJ7Fw4UI8++yz\nOHv2LB599FHU1tYCAKKjo/Hxxx8jODj46wP46MOniYiuZ77ad3p8pRAWFgaDwQAAUKvVSExMRE9P\nj9uc06dPIyMjAwAQHx+Pjo4OXLhwQY5zh09E9N/Df7YTOzo60NzcDJPJ5LZcr9fjnXfeQXp6Ohob\nG/HZZ5+hq6sLISEhUKlUyMrKwrx582CxWFBUVDTtY5eWlsqvzWYzzGbzN9oYIqLrldVqhdVq9fl6\nPB4+ump0dBRmsxm/+tWvcO+997qNjYyM4IknnkBzczOSk5Nx5swZ/OlPf0JKSgp6enqwePFiXLhw\nAdnZ2Xj55ZexYsUK9wA8fEREdM18te+csRQuX76MtWvXYtWqVdiyZcuMDxgdHY1PPvkEarXabXlZ\nWRnUajW2bt3qHoClQER0zRQ5pyCEQGFhIbRa7dcWwtDQEOx2OwBg165duOOOO6BWqzE+Po6RkREA\nwNjYGI4cOYLk5GQvxyciIm/yeE7h+PHj2LNnj7ysFADKy8vR2dkJALBYLGhtbcUjjzwClUoFnU6H\n3bt3AwD6+/uRl5cHAHA4HNi0aRNWrlzpy20hIqJvaVbnFHwagIePiIiumSKHj4iI6PuFpUBERBJL\ngYiIJJYCERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgK\nREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSPpWCz2ZCRkYGkpCTodDps3759ypxLly4hLy8P\ner0eJpMJp06dkmPV1dVISEhAXFwcKioqvJ+eiIi8SiU8fPJzX18f+vr6YDAYMDo6imXLlmH//v1I\nTEyUc5588kksXLgQzz77LM6ePYtHH30UtbW1cDqdiI+PR21tLSIiIpCWloaqqiq37wV89+HTRETX\nM1/tOz2+UggLC4PBYAAAqNVqJCYmoqenx23O6dOnkZGRAQCIj49HR0cHPv/8czQ2NiI2NhZRUVEI\nCAjAhg0bcODAAa9vABEReY//bCd2dHSgubkZJpPJbbler8c777yD9PR0NDY24rPPPkNXVxe6u7sR\nGRkp52k0GjQ0NEz72KWlpfJrs9kMs9l8bVtBRHSds1qtsFqtPl/PrEphdHQU9913HyorK6FWq93G\nnn76aTzxxBMwGo1ITk6G0WjEvHnzoFKpZh3iy6VARERTffUX5rKyMp+sZ8ZSuHz5MtavX4+HHnoI\n995775TxBQsW4NVXX5X3o6OjERMTg4mJCdhsNrncZrNBo9F4KTYREfmCx3MKQggUFhZCq9Viy5Yt\n084ZGhqC3W4HAOzatQt33HEH1Go1UlNT0dbWho6ODtjtduzbtw+5ubne3wIiIvIaj68Ujh8/jj17\n9iAlJQVGoxEAUF5ejs7OTgCAxWJBa2srHnnkEahUKuh0OuzevfvKA/v7Y8eOHcjJyYHT6URhYeGU\nK4+IiGhu8XhJ6ncSgJekEhFdM0UuSSUiou8XlgIREUksBSIiklgKREQksRSIiEhiKRARkcRSICIi\niaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJ\nHkvBZrMhIyMDSUlJ0Ol02L59+5Q5AwMDuOuuu2AwGKDT6fD666/LsaioKPn5zsuXL/d6eCIi8i6P\nn9Hc19eHvr4+GAwGjI6OYtmyZdi/fz8SExPlnNLSUvznP//B888/j4GBAcTHx6O/vx/+/v6Ijo7G\nxx9/jODg4K8PwM9oJiK6Zop8RnNYWBgMBgMAQK1WIzExET09PW5zwsPDMTw8DAAYHh7GokWL4O/v\nL8e5wyci+u8x63MKHR0daG5uhslkclteVFSEU6dOYfHixdDr9aisrJRjKpUKWVlZSE1Nxa5du7yX\nmoiIfMJ/5inA6Ogo7rvvPlRWVkKtVruNlZeXw2AwwGq14vz588jOzsY///lPLFiwAMePH0d4eDgu\nXLiA7OxsJCQkYMWKFVMev7S0VH5tNpthNpu/1UYREV1vrFYrrFarz9fj8ZwCAFy+fBlr167FqlWr\nsGXLlinjq1evxjPPPIPbbrsNAJCZmYmKigqkpqa6zSsrK4NarcbWrVvdA/CcAhHRNVPknIIQAoWF\nhdBqtdMWAgAkJCSgtrYWANDf34+zZ8/illtuwfj4OEZGRgAAY2NjOHLkCJKTk70cn4iIvMnjK4Vj\nx47h9ttvR0pKClQqFYArh4s6OzsBABaLBQMDA9i8eTM6OzvhcrlQUlKCjRs34tNPP8W6desAAA6H\nA5s2bUJJScnUAHylQER0zXy175zx8JGvsRSIiK6dIoePiIjo+4WlQEREEkuBiIgklgIREUksBSIi\nklgKREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGR\nxFIgIiKJpUBERBJLgYiIJI+lYLPZkJGRgaSkJOh0Omzfvn3KnIGBAdx1110wGAzQ6XR4/fXX5Vh1\ndTUSEhIQFxeHiooKr4f3FavVqnSEKeZiJmBu5mKm2WGm2ZuruXzBYykEBATgD3/4A06dOoX6+nq8\n8sorOH36tNucHTt2wGg0oqWlBVarFVu3boXD4YDT6cRjjz2G6upqtLa2oqqqasr3zlVz8QdgLmYC\n5mYuZpodZpq9uZrLFzyWQlhYGAwGAwBArVYjMTERPT09bnPCw8MxPDwMABgeHsaiRYvg7++PxsZG\nxMbGIioqCgEBAdiwYQMOHDjgo80gIiJv8J/txI6ODjQ3N8NkMrktLyoqwp133onFixdjZGQEf/nL\nXwAA3d3diIyMlPM0Gg0aGhq8FJuIiHxCzMLIyIhYtmyZ+Nvf/jZl7LnnnhNPPPGEEEKIc+fOiejo\naDE8PCzefvtt8bOf/UzOe+utt8Rjjz025fsB8MYbb7zx9g1uvjDjK4XLly9j/fr1eOihh3DvvfdO\nGf/www/xzDPPAABiYmIQHR2Ns2fPQqPRwGazyXk2mw0ajWbK91/pBSIimgs8nlMQQqCwsBBarRZb\ntmyZdk5CQgJqa2sBAP39/Th79ixuueUWpKamoq2tDR0dHbDb7di3bx9yc3O9vwVEROQ1KuHhV/Vj\nx47h9ttvR0pKClQqFQCgvLwcnZ2dAACLxYKBgQFs3rwZnZ2dcLlcKCkpwcaNGwEAhw8fxpYtW+B0\nOlFYWIiSkpLvYJOIiOgb88lBqVk6fPiwiI+PF7GxsWLbtm1efezNmzeLm2++Weh0Orns4sWLIisr\nS8TFxYns7Gxx6dIlOVZeXi5iY2NFfHy8eO+99+TypqYmodPpRGxsrHj88cfl8snJSfHAAw+I2NhY\nYTKZREdHx4yZOjs7hdlsFlqtViQlJYnKykrFc01MTIjly5cLvV4vEhMTxdNPP614pi9zOBzCYDCI\ntWvXzolcP/rRj0RycrIwGAwiLS1tTmS6dOmSWL9+vUhISBCJiYmivr5e0UxnzpwRBoNB3hYuXCgq\nKysVf57Ky8uFVqsVOp1OPPjgg2JyclLxTEII8dJLLwmdTieSkpLESy+9JIRQ9mdKsVJwOBwiJiZG\ntLe3C7vdLvR6vWhtbfXa43/wwQfixIkTbqXw5JNPioqKCiGEENu2bRNPPfWUEEKIU6dOCb1eL+x2\nu2hvbxcxMTHC5XIJIYRIS0sTDQ0NQgghVq1aJQ4fPiyEEOKVV14RxcXFQggh9u7dK37yk5/MmKm3\nt1c0NzcLIa6cvF+yZIlobW1VPNfY2JgQQojLly8Lk8kkjh49qnimq1588UWxceNGcffddwshlP83\njIqKEhcvXnRbpnSmhx9+WOzevVsIceXfcHBwUPFMVzmdThEWFiY6OzsVzdTe3i6io6PF5OSkEEKI\nBx54QLz++uuKP0+ffPKJ0Ol0YmJiQjgcDpGVlSXOnTunaC7FSuHDDz8UOTk58v7zzz8vnn/+ea+u\no7293a0U4uPjRV9fnxDiyg46Pj5eCHGleb/8SiUnJ0d89NFHoqenRyQkJMjlVVVVwmKxyDn19fVC\niCv/EX/wgx9cc7577rlH1NTUzJlcY2NjIjU1VfzrX/+aE5lsNpvIzMwU77//vnyloHSuqKgoMTAw\n4LZMyUyDg4MiOjp6ynKln6er3nvvPZGenq54posXL4olS5aIL774Qly+fFmsXbtWHDlyRPHn6e23\n3xaFhYXy/nPPPScqKioUzaXY3z6a7n0M3d3dPl1nf38/QkNDAQChoaHo7+8HAPT09LhdGXU1y1eX\nR0REyIxfzu/v74+goCB88cUXs87y5fd9KJ3L5XLBYDAgNDRU/lkTpTMBwC9/+Uv8/ve/h5/f//+Y\nKp1LpVIhKysLqamp2LVrl+KZ2tvbERISgs2bN2Pp0qUoKirC2NiY4s/TVXv37sWDDz6o+PMUHByM\nrVu34oc//CEWL16MG2+8EdnZ2Yo/TzqdDkePHsUXX3yB8fFxHDp0CF1dXYrmUqwUrp64VnL9SmUY\nHR3F+vXrUVlZiQULFiiey8/PDy0tLejq6sIHH3yAuro6xTP9/e9/x8033wyj0fi1ly0rkev48eNo\nbm7G4cOH8corr+Do0aOKZnI4HDhx4gR+/vOf48SJE7jhhhuwbds2RTNdZbfb8e677+L++++fMvZd\nZzp//jxeeukldHR0oKenB6Ojo9izZ4+imYArV28+9dRTWLlyJVatWgWDwYB58+YpmkuxUoiIiJjV\n+xi8KTQ0FH19fQCA3t5e3HzzzdNm6erqgkajQUREBLq6uqYsv/o9V6/CcjgcGBoaQnBw8IwZrr7v\nIz8/X77vYy7kAoCgoCCsWbMGH3/8seKZPvzwQxw8eBDR0dF48MEH8f777yM/P1/xXOHh4QCAkJAQ\n5OXlobGxUdFMGo0GGo0GaWlpAID77rsPJ06cQFhYmOI/U4cPH8ayZcsQEhICQNmf86amJvz4xz+W\nf4Zn3bp1+Oijj+bE8/TTn/4UTU1N+Mc//oGbbroJS5YsUfS5UqwUlHgfQ25uLt544w0AwBtvvCF3\nyrm5udi7dy/sdjva29vR1taG5cuXIywsDAsXLkRDQwOEEHjrrbdwzz33THmsv/71r8jMzJxx/eJr\n3vehZK6BgQEMDg4CACYmJlBTUwOj0aj4c1VeXg6bzYb29nbs3bsXd955J9566y1Fc42Pj2NkZAQA\nMDY2hiNHjiA5OVnRTGFhYYiMjMS///1vAEBtbS2SkpJw9913K/rvBwBVVVXy0NFXH+e7zpSQkID6\n+npMTExACIHa2lpotdo58Tx9/vnnAIDOzk6888472Lhxo7L//2Y8E+JDhw4dEkuWLBExMTGivLzc\nq4+9YcMGER4eLgICAoRGoxGvvvqquHjxosjMzJz2Mq/f/va3IiYmRsTHx4vq6mq5/OplXjExMeIX\nv/iFXD45OSnuv/9+eZlXe3v7jJmOHj0qVCqV0Ov18nK9w4cPK5rr5MmTwmg0Cr1eL5KTk8Xvfvc7\nIYRQ/Ln6MqvVKq8+UjLXp59+KvR6vdDr9SIpKUn+zCr9XLW0tIjU1FSRkpIi8vLyxODgoOKZRkdH\nxaJFi8Tw8LBcpnSmiooKeUnqww8/LOx2u+KZhBBixYoVQqvVCr1eL95//33FnyuPb14jIqLvF37y\nGhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiLpfwGexzRJjCWLBQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1fa10d0>"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def MfmuIfun(mu):\n",
      "    \n",
      "#     MfMui = mesh.getFaceMass(1./mu)\n",
      "    \n",
      "    MfMui = mesh.getFaceInnerProduct(1./mu)*1/3\n",
      "    # There are scaling issue for Face inner product!! and bit slower than getFaceMass\n",
      "    \n",
      "    MfMuI = 1/MfMui.diagonal()\n",
      "#     dMfMuI = Utils.sdiag(MfMuI**2)*mesh.aveF2CC.T*Utils.sdiag(mesh.vol*1./mu**2)\n",
      "    dMfMuI = Utils.sdiag(MfMuI**2)*mesh.aveF2CC.T*Utils.sdiag(mesh.vol*1./mu**2)\n",
      "    \n",
      "    return MfMuI, dMfMuI"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d_mu = 0.8*mu\n",
      "# I am not sure why it does not give us second order accuracy\n",
      "Tests.checkDerivative(MfmuIfun, mu, dx=d_mu, num=5, plotIt=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "==================== checkDerivative ====================\n",
        "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n",
        "---------------------------------------------------------\n",
        "0\t1.00e-01\t4.299e-09\t\t1.152e-23\t\tnan"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1\t1.00e-02\t4.299e-10\t\t7.723e-24\t\t0.174"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "2\t1.00e-03\t4.299e-11\t\t1.170e-23\t\t-0.180"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "3\t1.00e-04\t4.299e-12\t\t8.749e-24\t\t0.126"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "4\t1.00e-05\t4.299e-13\t\t1.107e-23\t\t-0.102"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "========================= PASS! =========================\n",
        "Once upon a time, a happy little test passed.\n",
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 55,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def derBbc(chi):\n",
      "    Bbc, const = CongruousMagBC(mesh, np.array([1., 1., 0.]), chi)\n",
      "    bc = Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*Bbc  \n",
      "    print const, np.linalg.norm(Bbc)\n",
      "    dBbc = lambda x: (Utils.sdiag(prob.mesh.vol)*(prob.mesh.faceDiv*prob._Pout.T*const*Bbc))*(np.inner(prob.mesh.vol, x)) \n",
      "    return bc, dBbc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a, b = derBbc(chi)\n",
      "print np.linalg.norm(a), np.linalg.norm(b(chi))\n",
      "print np.linalg.norm(chi)\n",
      "print np.array([1., 1., 0.])\n",
      "print mesh"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "7.65830184557e-07 0.0659506783728\n",
        "55.718105908 55.6914561973\n",
        "1.6089045963\n",
        "[ 1.  1.  0.]\n",
        "  ---- 3-D TensorMesh ----  \n",
        "   x0: -488.58\n",
        "   y0: -488.58\n",
        "   z0: -476.08\n",
        "  nCx: 31\n",
        "  nCy: 31\n",
        "  nCz: 30\n",
        "   hx: 71.40, 54.93, 42.25, 32.50, 23*25.00, 32.50, 42.25, 54.93, 71.40\n",
        "   hy: 71.40, 54.93, 42.25, 32.50, 23*25.00, 32.50, 42.25, 54.93, 71.40\n",
        "   hz: 71.40, 54.93, 42.25, 32.50, 22*25.00, 32.50, 42.25, 54.93, 71.40\n"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d_chi = 0.8*chi\n",
      "# I am not sure why it does not give us second order accuracy\n",
      "Tests.checkDerivative(derBbc, chi, dx=d_chi, num=5, plotIt=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "==================== checkDerivative ====================\n",
        "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n",
        "---------------------------------------------------------\n",
        "7.65830184557e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0659506783728\n",
        "7.0907489096e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0712240073538\n",
        "0\t1.00e-01\t4.455e+00\t\t1.705e-04\t\tnan\n",
        "7.59749260143e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0664780294308\n",
        "1\t1.00e-02\t4.455e-01\t\t1.705e-06\t\t2.000\n",
        "7.65217717349e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0660034136602\n",
        "2\t1.00e-03\t4.455e-02\t\t1.706e-08\t\t2.000\n",
        "7.65768893742e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0659559519034\n",
        "3\t1.00e-04\t4.455e-03\t\t1.749e-10\t\t1.989\n",
        "7.65824055034e-07"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0.0659512057259\n",
        "4\t1.00e-05\t4.455e-04\t\t2.506e-12\t\t1.844\n",
        "========================= PASS! =========================\n",
        "The test be workin!\n",
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 58,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}