Modify Gudni's notebook:

- use SolverLU in simpeg - visualize fields Add gitignore
2026-10-28 13:30:34 +08:00 · 2015-01-29 10:35:39 -08:00
1 parent 47e7132f26
commit 8e14db7dba
@@ -0,0 +1,44 @@
+*.py[cod]
+
+# C extensions
+*.so
+
+# Folders
+docs/build/*
+
+# Packages
+*.egg
+*.egg-info
+dist
+build
+eggs
+parts
+bin
+var
+sdist
+develop-eggs
+.installed.cfg
+libgit
+lib64
+__pycache__
+
+# Installer logs
+pip-log.txt
+
+# Unit test / coverage reports
+.coverage
+.tox
+nosetests.xml
+
+# Translations
+*.mo
+
+# Mr Developer
+.mr.developer.cfg
+.project
+.pydevproject
+
+*.sublime-project
+*.sublime-workspace
+docs/_build/
+*.ipynb_checkpoints
@@ -0,0 +1,356 @@
+{
+ "metadata": {
+  "name": "",
+  "signature": "sha256:746cad20da01888f7eb04de3a80de5697a42a8f6e49a6b2f2c5c73200761f12a"
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+  {
+   "cells": [
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from SimPEG import *\n",
+      "from scipy.constants import mu_0\n",
+      "def omega(freq):\n",
+      "    \"\"\"Change frequency to angular frequency, omega\"\"\"\n",
+      "    return 2.*np.pi*freq"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 3
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "%pylab inline"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "Populating the interactive namespace from numpy and matplotlib\n"
+       ]
+      }
+     ],
+     "prompt_number": 4
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# M = Mesh.TensorMesh([[(100.,32)],[(100.,34)],[(100.,18)]], x0='CCC')\n",
+      "M = Mesh.TensorMesh([[(100.,20)],[(100.,20)],[(100.,18)]], x0='CCC')"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 5
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "print M"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "  ---- 3-D TensorMesh ----  \n",
+        "   x0: -1000.00\n",
+        "   y0: -1000.00\n",
+        "   z0: -900.00\n",
+        "  nCx: 20\n",
+        "  nCy: 20\n",
+        "  nCz: 18\n",
+        "   hx: 20*100.00\n",
+        "   hy: 20*100.00\n",
+        "   hz: 18*100.00\n"
+       ]
+      }
+     ],
+     "prompt_number": 6
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# Setup the model\n",
+      "conds = [1e-2,1]\n",
+      "sig = Utils.ModelBuilder.defineBlock(M.gridCC,[0,0,-300],[500,500,-100],conds)\n",
+      "sig[M.gridCC[:,2]>0] = 1e-8\n",
+      "sigBG = np.zeros(M.nC) + conds[0]\n",
+      "sigBG[M.gridCC[:,2]>0] = 1e-8\n",
+      "colorbar(M.plotImage(log10(sig)))"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "pyout",
+       "prompt_number": 7,
+       "text": [
+        "<matplotlib.colorbar.Colorbar instance at 0x000000000A9B6088>"
+       ]
+      },
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAEKCAYAAADZ8ATAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFX28PFvd2cRZEcMIWwawWELhFVlnEQHMTIjyrih\nIzCCC6KC2yCIksTfOw7ighuiqCgoiNFRBBE3JCAqiwjIFhZJIIkJSkJYZAkJ/f5xq5NOp5fqpLfq\nPp/nqSfV1VXV9xg8qT51614QQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIUQ4uwX4ETgK/Ap8\nBgzU3usMfAD8DpQBm4EHAHPgmymEEOJB4ABwLdAAsAB/A54CEoFDwDNAnLZ/Z+BdoGnAWyqEEBGu\nKerq+joX778LLAlcc4QQQriTBpzGdamjCBgVuOYIEZqkFihCRUvgIHDGzftFgWuOEKFJkrYIFSXA\nObj+N1kCtAlcc4QQQrjTFDiG65r2O0hNWwghQsqDQDFwDdAQiAauQvUeOR91tT2d6t4jF6CSufQe\nEUKIILkFWI+66i5CXV1fpL3XGchC1b7LgE3AeKTMJ0QNk4FtwBZgARALtAC+AnYBXwLNHPbfDeQA\ng+2299HOsRt4we+tFkII/0tD5brdwCMu9nlRe38zkOzvBnUE9qISNcD7qG5X04GJ2rZHgGnaelfU\n1U+0duwewKS9tw7or61/hgpWCCGMyoLKcR1ROW8T0MVhnyGofAcwAFhT3w/19LXyCKrvbEMgSvv5\nKzAUmKvtMxf1BBuoWuR72jF5qIAGAPFAY1TiBphnd4wQQhhRf1SOy0PlvIWoHGjPPleuRVUl4qgH\nT0m7FHgW2I9K1mWoskgc6nFjtJ+2RrQBCuyOLwASnGwv1LYLIYRRJQD5dq9t+c7TPm3r86GeknYi\ncD/q8r8N0Ai41WEfq7YIIUQk0Zv3TA6v65Uvozy83xf4HtXVCuAj4GJUt6zW2s944Dft/UKgnd3x\nbVF/WQqp+delrbatlsTEROsvv/yiPwIhRKT6BdXts87OAutJ/bsfBZrYvXbMd+2oWVFwto/L3KeX\np6SdAzyOGnHtJDAIVZf+A3VD8int5yJt/8WoHibPob4WdNL2t6Lq4wO01yNQd1RrUQk7o47hhB6r\nNZ2MjAwyMx3/2BpbOMZliykjIyPYTfG5cIzLZDIl1vccJ4H/p3Pfx9R9OXs/onJcR1T5+CbgZod9\nFgP3ourdF6FKzAeoB09JezPqpuGPqDEhfgJmoxqfBYxBFeFv1Pbfrm3fDlQA46j+KjAOeBv1B+Az\n4PP6NFwIIXwhuu6HVqAS8heoniRvAjuAu7T3X0PluiGoG5Z/ALfV/eMUT0kbVPe+6Q7bSlFX3c48\nqS2ONgA99DdNCCH8T08SdGOZtth7zeH1vfX7iJrkSbIASE1NDXYT/CIc4wrHmCB84/KFBjqXUCFJ\nOwDC9X+YcIwrHGOC8I3LF6J1LqGint8MhBDC2IyWBI3WXiGE8KlQuorWQ5K2ECKiGS0JGq29Qgjh\nU3KlLYQQBiJJWwghDCSUuvPpIUlbCBHRjJYEjdZeIYTwKSmPCCGEgRgtCRqtvUII4VNypS2EEAZi\ntCRotPYKIYRPGe1KWwaMEkJEND+N8tcCNZ/uLuBL1IS+zkwGtgFbUBPIxHo6sSRtIURE89Mof5NQ\nSbszsFx77agjcAfQGzXXgAUY7unEkrSFEBEtSufipaHAXG19LnCtk32OAKeBhtpHNETH/JFS0xZC\nRLRovVmwwqvTxlE9F+QB7bWjUuBZYD9wAjVt2deeTixJWwgR0aJcZMFvK2H1GbeHfgW0drJ9isNr\nK9Vz5dpLBO5HlUkOAx8A/wTmu/tQPUn7QtRMwjbno2Zofxd4H+hA9eS+Zdo+k4HRQCUwHlWIB+iD\nmtz3LNSElxN0fL4QQvhNtMX59sstcLnd62l/1NrlCjenPYBK6MVAPPCbk336At8DJdrrj4BL8JC0\n9dS0dwLJ2tIHOA58jOtCe1fUVPJdgTTgFcCkvTcLNYN7J21J0/H5QgjhN1FR+hYvLQZGaeujgEVO\n9skBLkJ1TjGhJkvf7unE3t6IHISaCj4f14X2a4D3UAX2PG3/Aai/No2Bddp+83BenBdCiICJjtW3\neGka6kp8F+qCfZq2vQ2wVFvfjMqDPwI/a9tmezqxt38/hqMSMrgutLcB1tgdUwAkoJJ4gd32Qm27\nEEIEj3/u7JWiLnId/Qr8ze71dG3RzZvmxgBXA484ec9VoV0IIUKbwbpjeNPcq4ANwO/aa1eF9kKg\nnd1xbVFX2IXauv12F30SV9itdwTO86KZQohwlJ2dTXZ2tu9PHMZJ+2aqSyNQXWh/ipqF9sWoxzGf\nQ5U/OqHq2FZUZ/IB2usRwIvOP+oyL5olhIgEqamppKamVr3OzMz0zYld9B4JVXqT9tmo+swddtum\nAVmo3iB5qC5/oO5+Zmk/K4BxVJdOxqG6/DVAdfn7vM4tF0IIXwjTK+0/gHMctrkqtAM8qS2ONqCe\nsRdCiNDgfc+QoDLY3xghhPAxg2VBgzVXCCF8zGBZ0GDNFUIIHwvTG5FCCBGeDJYFDdZcIYTwMYNl\nQYM1VwghfMxgWdBgzRVCCB+TLn9CCGEgBsuCMkekECKyWXQu3rkBNct6JWriXleaAR8CO1BPkV/k\n6cQRnbSLix+id+94AFau/BfDh3eveq9r11ZkZV3Pzp33UlHxOLNnX13r+JSUDlRWTq213HZbr4DF\n4Ki+MQGYzSYeeWQgOTn3cOLEFIqLH2LmzCEBab8r9Y3rrbeucfq7qqh4nJYtGwQsDhGC/DOz7xZg\nGLDKw34voIb06AIkoZK3x+ZGpMTE5jRsGM3GjUVER5vp27cNq1fvr3q/QYMo8vIO88knO3nwwYux\nWl2PPJuc/BpFRUerXh85csqvbXfFVzG9/fY1DBjQlokTv2LTpmIaN46lY8dmgQqjFl/ENX78MiZO\n/KrqtclkYtGimzh2rJySkhMBiUOEKP9kwRwd+zQFLqV6hpsK1FyRbkVs0h44sD1r1xZitUK/fgmU\nlBynoOBI1fsbNhSxYUMRAGPGJLs918GDx/n99+N+ba8evogpNbUjw4d3JynpVXJyDlZt37rV2RR3\ngeGLuI4eLefo0fKq1506tWDAgLbccMMH/m28CH3Be7jmPNRQ128BPVFjM01ATenoUsQl7UOHHsFq\ntRIbG4XZbKK0dCLR0RZiYy2Ulk7EaoWWLb2aSILVq2+jYcNo9uwp5bXXNvDOOz97PsiHfBnTddd1\nYe/eQwwenMiSJTcTE2Phhx/yefjhr2okykDwx+/K5q67+lJcfIxFi/RcEImw5iILZv8K2UVuj3Q1\nG/ujwBKdn9wbuBdYDzyPmmt3qqeDIkpS0ixMJhNr1oxh7NilbNpUzMKF17FgwVY++cS7/4F//fUo\n48Yt5ccff+XMGStDhnRi9uyrueCCFqSnZ/snACd8GVNiYnPat2/Krbf2YMyYxZSXV/Kf/1zON9+M\npHv3WZSXV/opitp8GZe9mBgLo0b15LXXNnDmjEy4FPHOcr459Xy12GRurLWLu9nY9SjQlvXa6w+p\nniDdpYhL2vn5R+jR41yioy0sWbKTRo1i6NWrNUOHLuTgQe9KHLt3l7J7d2nV640bi7FYzDz88MVk\nZq4MWELwZUxms4nY2ChGjlxUVR656aYPKSp6iCFDOgX0ytSXcdm7/vquNG9+FrNnb/Bha4Vh+b88\nYnKxvRg1SXpn1ATAg1A9TtyKqKS9devdtG/flKgoM9HRFg4fnlSVpPbuHQ9Aly4zKSw86uFMrq1d\nW8DZZ8fQqlVDDhz4w1dNd8nXMRUVHcNqtdaoZx88eJyDB4/Tvn1Tv8TgjD9/V2PH9uGLL35h/36P\n93xEJPBPFhyGmpnrHNTs6xtRUza2AV6nenLf+4D5qDl4fwFu83TiiEraaWnziYmxMGfOUJYt20NW\n1jbS01M4daqSadNWAypp1Ufv3vEcP366XleC3vB1TKtW7WPkyJ507tySXbtKAGjRogHnnNOQvLwy\nv8TgjL9+V126nMPAge0ZNux9XzdZGJV/suDH2uLIcTb2zUA/b04cUUm7oOAIZrOJpKQ47rzzU3Jz\ny+jRI46MjGxyc2smpKgoM926tQKgceNYWrZsQM+ecZSXV7Jjh7oKvf/+i9i3r4zt23/HaoUrr0xk\nypRLefnl9VRWBqY04uuY3ntvK1OmXMqcOUMZP/5zTp+u5KmnBrF7dwnLlu0OSEz+iMvmrrv68uuv\nR1myZGfAYhEhToZmDW3Jya05daqSXbtKaNIklm7dWrFq1b5a+yUkNOann+4CwGq10rt3PMOGdSEv\nr4zERDUfscVi4skn/0q7dk04ffoMu3eXMH7858yZU/uOhVFiOnmygkGD3uH5568kO3sUx4+fJjs7\nj0GD3uH06TOGjQvgrLOiGDEiiZdeWoebbvci0hgsC7oqkDtqBrwBdENN0nsbsBt4H+hA9cS+tkug\nycBo1COc44Evte19UBP7noV6CmiCk8+yQoZXQYQyqzUdAJPJRzNHh4hwjMsWkzAGk8kE+nOYK1br\nEzo/T3XEq+/n1Zvex9gdH7XMQXVN+Qp153M51V1VugI3aT/TgFeoDnQWavb2TtqSVu8IhBCiPmJ1\nLiFCT9K2PWo5R3tte9RyKDBX2zYXuFZbvwZ4DziNugLfAwwA4oHGwDptv3l2xwghRHD4Z+wRv9HT\nFGePWt4PxAEHtH0OaK9BdWlZY3d8AZCASuIFdtsLte21hOPX1HCMCcI3LhFBQigh66Gnua4etbRn\n1RafyMjIqFpPTU0lNTXVV6cWQhhUdnY22dnZvj9xGPYecfao5WTU0zyttZ/xgG1EoUKgndLine truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x3eba8d0>"
+       ]
+      }
+     ],
+     "prompt_number": 7
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# Get the mass matrix \n",
+      "# The model\n",
+      "Msig = M.getEdgeInnerProduct(sig)\n",
+      "MsigBG = M.getEdgeInnerProduct(sigBG)\n",
+      "Mmu = M.getFaceInnerProduct(mu_0, invProp=True)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 8
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "freq = 1e1\n",
+      "C = M.edgeCurl\n",
+      "A = C.T*Mmu*C - 1j*omega(freq)*Msig\n",
+      "ABG = C.T*Mmu*C - 1j*omega(freq)*MsigBG"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 9
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# Need to solve x and y polarizations of the source.\n",
+      "from simpegMT.Utils import get1DEfields\n",
+      "# Get a 1d solution for a halfspace background\n",
+      "mesh1d = Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))\n",
+      "e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
+      "# Setup x (east) polarization (_x)\n",
+      "ex_x = np.zeros(M.vnEx,dtype=complex)\n",
+      "ey_x = np.zeros((M.nEy,1),dtype=complex)\n",
+      "ez_x = np.zeros((M.nEz,1),dtype=complex)\n",
+      "# Assign the source to ex_x\n",
+      "for i in arange(M.vnEx[0]):\n",
+      "    for j in arange(M.vnEx[2]):\n",
+      "        ex_x[i,j,:] = e0_1d\n",
+      "eBG_x = np.vstack((Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))\n",
+      "rhs_x = ABG.dot(eBG_x)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 10
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "M.vnEy"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "pyout",
+       "prompt_number": 11,
+       "text": [
+        "array([21, 20, 19])"
+       ]
+      }
+     ],
+     "prompt_number": 11
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# Setup y (north) polarization (_y)\n",
+      "ex_y = np.zeros(M.nEx, dtype='complex128')\n",
+      "ey_y = np.zeros((M.vnEy), dtype='complex128')\n",
+      "ez_y = np.zeros(M.nEz, dtype='complex128')\n",
+      "# Assign the source to ex_x\n",
+      "for i in arange(M.vnEy[0]):\n",
+      "    for j in arange(M.vnEy[2]):\n",
+      "        ey_y[i,j,:] = e0_1d        \n",
+      "# eBG_y = np.vstack((ex_y,Utils.mkvc(M.r(ey_y,'Ey','Ey','V'),2),ez_y))\n",
+      "# rhs_y = ABG.dot(eBG_y)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 12
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "eBG_y = np.r_[ex_y,Utils.mkvc(ey_y),ez_y]\n",
+      "rhs_y = ABG.dot(eBG_y)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 13
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "We are using splu in scipy package. This is bit slow, but on the cluster you can use mumps, which might a lot faster. We can think about having better iterative solver. "
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "%%time\n",
+      "# Solve the systems for each polarization\n",
+      "Ainv = SolverLU(A)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "Wall time: 37.7 s\n"
+       ]
+      }
+     ],
+     "prompt_number": 14
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "%%time\n",
+      "e_x = Ainv*rhs_x"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "Wall time: 202 ms\n"
+       ]
+      }
+     ],
+     "prompt_number": 23
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "I want to visualize electrical field, which is a vector, so I average them on to cell center. Also I want to see current density ($\\vec{j} = \\sigma \\vec{e}$)."
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "e_x_CC = M.aveE2CCV*e_x\n",
+      "j_x_CC = Utils.sdiag(np.r_[sig, sig, sig])*e_x_CC"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 26
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Then use \"plotSlice\" function, to visualize 2D sections"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "fig, ax = plt.subplots(1,2, figsize = (12, 5))\n",
+      "dat0 = M.plotSlice(abs(e_x_CC), vType='CCv', view='vec', streamOpts={'color': 'k'}, normal='X', ax = ax[0])\n",
+      "cb0 = plt.colorbar(dat0[0], ax = ax[0])\n",
+      "dat1 = M.plotSlice(abs(j_x_CC), vType='CCv', view='vec', streamOpts={'color': 'k'}, normal='X', ax = ax[1])\n",
+      "cb1 = plt.colorbar(dat1[0], ax = ax[1])"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [
+      {
+       "metadata": {},
+       "output_type": "display_data",
+       "png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAFRCAYAAABDmqrPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXeYZEd19//p7sm7K62kVVqlVUJIgCQECJkobIxlbIOx\nMcE4gP2+tvEjzIvBYIyDMDhgko2xMD8wGJNBYJBAApRWAi0SyjmtpNVGbd6d1OGm3x/VvdOh6py6\n3Xd6Zqfv93nmmZmuc6vqhq77ved+zzmQI0eOHDly5MiRI0eOHDly5MiRI0eOHDly5MiRI0eOHDly\n5MiRI0eOHDly5MiRI0eOHDly5MiRI0eOHDly5MiRI0eOHDly5MiRI0eOHDly5MiRI0eOHDly5MiR\nI4cX3gz8uOn/KWDNgswkR44cOXL44M3k63aOHD2huNATyLGk8SJgHbAP2A38BHiuw3YFsCHj8S8G\nbgMqwOct7b8APATMANcBJ2Y8fo4cOXIcbMjX7Rw55hE58c4xXzgE+B7wb8BhwHHA+4FqH+ewBfgA\n8DlL2yrgW8D7MPO7Dfh6/6aWI0eOHIsO+bqdI0eOHAcpngvsFdrfTOsryxg4pf73OPBRjCdlX91u\nrN52AcYbsxe4C3ipx1w+QKfn5I8wnpwGJoBZ4Gke/eXIkSPHUkS+bufIMc/IPd455gsPAxHw38BF\nGO+ELz4CPBv4OeBw4C8wC/xxGG/M39f7exfG+7FK6a9g+ewZwN1N/88C64Fnpphnjhw5ciwl5Ot2\njhzzjJx455gvTGG0ggnwGWAH8F3gKGW7IvAW4O3ANszCfTNQA34HuBL4Qd32GsyrxlcqfSaWz5YB\nk22fTQLLlb5y5MiRY6kiX7dz5Jhn5MQ7x3ziIcxifALGI7Ea+Fdlm1WY15OPWdpOAn4L87qy8fNC\n4BilT5vnZBqjZ2zGoZgbT44cOXIMKvJ1O0eOeUROvHP0Cw8DX0B/JbgLE81+mqVtI/BFzOvKxs8K\n4F+UPm2ek/uBc5r+XwacWv88R44cOXLk63aOHJkjJ9455gtnAH+O0feB8Z68Efipsl2MiWb/GHAs\nUMJoBkeALwG/Bryi/vkYcGHTGO1o2AzV/x6t/wb4X8zN5DfqNn+HCfp5xHsPc+TIkWNpIV+3c+TI\nkeMgxWpMmqfNmNeDm4FPMafF+33gxib7iLno+DHg4/Vt9gFrmYuOP7/+/26M/vAKzM3BhkswN4Tm\nn79tav8F4EFMgE6eDzZHjhyDjnzdzpEjR44cOXLkyJEjR44cOXLkyJEjR44cOXLkyJEjR44cOXLk\nyJEjR44cOXLkyJEjR44cOXLkyJEjR44cOXLkyJEjR3awJag/yHFGYlKP5siR4yDH3cC53W48Bkkl\n/WZ7MeWuc/QNRyewfaEnkSNHjmzQ9brd5ZoNB9m6vQSJNwk8sNBz6DM+CVy80JNoQtCHMT4FvLUP\n4/ginOf+Pw38cQr7fpyD+cYF0NsalXww5QZ/bX4txXVxMSMxGeQGCdcDL1voSfQZ+T4PBi6B7tfQ\n1Gs2HHzr9tBCTyBHjhw55gvDCz2BHDly5MjhjUFYs/PKlTly5MiRI0eOHDly9AG5x3tJ4PyFnsAC\n4LkLPYE+4zkLPYGDEvkCl2NxYs1CT2ABsGahJ7AAWLPQEzjoMAhr9iDs4wBgEIn38xZ6An3GoD1o\nZINBeG2Z42DEyQs9gQVAvs85dAzCmp0T7xw5cixZ5Atcjhw5chw8GIQ1e0A03rdhMtxoeBS40dF2\nJfCUsO164IYUc7oZuL/ts1uBe4VtJoHLlH6vArYI7RuA65Q+vgqUhXbteM4A31DGALiDzmPQwGbg\nGmX7bwLTHuM08AUgbvq/itnXdtwO3KP0dQ9m/r7YhPu436f0VQO+pvS/FtgotG+u20i4DJCSOd0O\nPCi078F8T9qxDbhWGXt+MJzyJ8diwR707FS3IGfveRJz3btQwaxlEh6oz8WF3cjfCTDrelVo34j8\n3W1GBPzU8vmD9bm4sBf3WtvA7cjf/82Y+4cLIea+JuFhYKfl8y1K3za4jttW4HFhu9uQ99MGn/Pc\njjsx90IXtgBPCO0h5hqXoJ13n3nfjF9WrjTXafdIu2YfjOv2gBDvG9EvYDAE6CpH21eQCe0DwPdT\nzOkaOsnWdZjFz4XdwH8p/X4N+WbzIPo8/xVD9ly4AfiZ0L4X+LwyBtiPQQOPAN9Ttr8U/wU0AD5B\n6yU/CXzWYnsdZtGUsA79JtOMh4AfCH1J1+d+4HNK/9/GEA0XHhHGb+BSzI3dhbXIDyTbga9bPt8A\nXKGMPT8YSvmTY7FgF/p38Hpk4v0oMrEpAz9WxrgdmdjsQHfqaPNcj0wUmxHW+2vHHZhj5sIO4C6l\n77XIDwiP1X9cqKI7n+6qz6UdT5CeeLvO7+OYY+rCtcjrnA3bkJ1iNtyIfH/SjmcFtyOwAe369Jn3\nNbQ6pFzQvk/ZIO2afTCu2wfjnOcRCe5UkFIbmAu3lGKsmM7DHyp9hJZt2hEpfQQefQTAiNBe67G9\ngQow5mgrA+PK9j42DcwCE22fzVg+a/Rr+7wZ08CRnmM3xlou9LVKGcu1ra+Nbf+bkWDOx6hgU1Xa\nXdeNzzU3PzgYvSE5wFyPPjbSmqz1oW0PZp2W/FNau49NhP+VGmL/jmnrfoT+HdTuL9p9IUDfj5rD\npqL0bUMVWOYYQ+pLW8dsqAErUm4j3d8afWrrrXZMtGOundPEw6aBiPTHLT0GYc0eIOLtk1tdWoh9\nFuA0+dvH6byIfRZoH9IsXbq+i2MvxNp3YaviXpg0EtggitLC1gwb8XQR7BnsC3q7zRrPsbU+p5W+\nfOajEWvtYaKKuXFL15cP8bZt73PNzQ8GaIHL0QEfYq6t2T4Ol36Q92ZbGzTirTl1GjbSNyZEXkN8\nvufSw7nmXGiH6z5UBVY6tmlIKtKuDGXSCQQS9PWyikzmfRxYmo12ThrXr8+++fCP3jEIa/Yg7GMK\njOCuOnos8sU5BByWYqz9dC6iy5A9uDFwotLvMcintQAcIbSHwDOQF+kJ5EUywI+USgtTjOxRrgDP\nwv8tg4t4n2qxdXnC2200MuxrP4r7RgFmX09R+l+FPGef4/lMZYxlyPscAcdZPl844j0I3pPBxVFK\n+wgyKUmQ3zSBIUbSGlNAJ4xHId87hvH3JLoI9nJ0z+YhSt9HOvpuYAj5/hSin5MJ7OfE9y2pzzYl\n3GthN97ubrYLgdXI50SaJ5j9096qrkBf5aTrM8S+ZtugPdxlg0FYswdE4z2Bn2e0AuxztEm6aTBf\nzKkUc7J5OfYia60C5ABPMDo56cs+gyGh0hiPKmPYNHrNKCPrzhqQvNqTyLrICrKmuR2zdN40puvj\ntOMI9JtUWuItSUE2IS/qk+hBpA8I/YMJEJOurQpyDAMYvaB0bVUw13A7CsDRSt/zg6UepLN04UNo\ntyN7myvIQWMJ9uu1Ga77QQMBepzJVmTCMoufxhbc5Ge3MkYVOW4nQv9+a/e3Gvo6tRP7N80lQZHg\nIsP7cB+LfhHvKnJQLph5Sse7hpzkAMx3QHpgKSPr2SPk2IB22/4Q76W+bg+Ix3sGv13VNN6aDCSN\n1MR2EUfKGCH6Zeaj0+uHFMVnkdI03tLDUhp9N6TTeD8EvF7pz0d33YyKYK/1NY2f9EWymcX9Ngf8\njqemOXRdn7OkeyjNDgOywC1BJMgZIRo22prbD6mJtGYn6HKUNITGZeuj8e41fiiLe4fLJkuPt3T/\nqWI80WnRDfHW7DUbn/uo9sCindcsrr1sMQhr9iDsI36BOg07SeOtLeJpgyvbF2xNQ+Vz4fdKvH11\nZb0E4fnYVZBfs2VBvF26aE0vDUbSk8bjvV/oUyPNUmAmzHmyNM+HtE8+enmfG4XrRpgHV+ZIA981\nu5c+sgiu1BwyjTGy0nh3S7x99Nvz7dQBN1Hshni7HAGSg6BKuhS0DWgxR7ZxfIi3FiuVRXClJg/0\n5S0+cQK9YxDW7AGRmoB/cKXU1ovnox02ku3j8dYufC2DhE+7dulrC4Kvx1tanKTAS8iOeNsWJR8Z\nyT0eNs2QvNq9Em/f4EvpePl6vKXz6rox+9zUFzUuwrwGeRR4j8PmE/X2u4FnN33+Ocz7YFtOr7dh\n8nveB3woq8kODnoNnvQJjOy1D21Nb9j0Sn40p40P8e7V4+3zPZeCK9OuES69uHZf6VZq4hvI7ztO\nt46MBmL0Y67d79Nce8NKXzl8sTSP4lvPbP1/3eEwfiQ8+0y7fQP3HQO7t8NLLXbfGIGfPw1WOfq4\n5yiYnIQXKWM08IMJOP1EOLXJ/kfL4KST4BRHH5uehPsPhYuEMb4M/PpZsMzxOu3WlTA04T4W+4fh\nh8vgdU3t7RKxK4fgvNPhGEcfj/wUdhwDLzrbPU+A/x2CC8+GFZZ+bhqHY06HUx19PLUH7j4CfkkZ\nA8za9MRdsP9sOLfJ/r4VEI/B2W19fKMMv/Z8GBLI7uUV+KXzYbQenKWlhV0LPOs5cMR5c3Nq4KpZ\nePmLYPhQ+7YPLIehVXDKefb22cfhZ4fDhU3t7fL4O0fghPNg1fPtfezaDZuOhmc72gHWFeCc82HZ\nGfb93XIXTO+AM9r62HAtRMvhVKFvG2zpilMigwWuBHwSeDlGBH8rcDmtVSleCZwGnA48H/gUcEG9\n7fPAvwP/09bvy4BXAWdjzlaa3JQ5vJBFVpNeM5JkkW6wGS6ipDllInSZWBbpBrUMGrY0utCdx3sH\nduKqSU36pfGeb+LdON7SNewjNfFdJWfoh692aZLSVuQe7xYIC3ESQ0HoI/HxbLTZF9oWyjjs/Kyl\nPYCCclnGARSFxS+uQUn4MsdVKCoLYFSFkrBgRBUoeXgHgkkYdnhqo1koCdKIcAaGUnicq7Line truncated
+       "text": [
+        "<matplotlib.figure.Figure at 0x12a396d8>"
+       ]
+      }
+     ],
+     "prompt_number": 32
+    },
+    {
+     "cell_type": "markdown",
+     "metadata": {},
+     "source": [
+      "Is it reasonable?: Based on that you put resistive target that makes sense to me; current does not want to flow on resistive target so they just do roundabout:). And see air interface. It is continuous on current but not on electric field, which looks reasonable. "
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "# # Get a 1d solution for a halfspace background\n",
+      "# mesh1d = simpeg.Mesh.TensorMesh(M.hz,M.x0[2])\n",
+      "# e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
+      "# # Setup x (east) polarization (_x)\n",
+      "# ex_x = np.zeros(M.vnEx)\n",
+      "# ey_x = np.zeros(M.vnEy)\n",
+      "# ez_x = np.zeros(M.vnEz)\n",
+      "# # Assign the source to ex_x\n",
+      "# for i in arange(M.vnEx[0]):\n",
+      "#     for j in arange(M.vnEx[2]):\n",
+      "#         ex_x[i,j,:] = e0_1d\n",
+      "# eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 33
+    }
+   ],
+   "metadata": {}
+  }
+ ]
+}