simpeg/notebooks/ThreeMeshDC.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:4ee97c21ba8374164213c87af95095c4694b605af46490529f6e0f5ff4051dea"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from SimPEG import Mesh, Utils, np, SolverLU\n",
      "import matplotlib.pyplot as plt\n",
      "# %pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib\n",
      "from matplotlib.mlab import griddata"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "matplotlib.rcParams.update({'font.size': 16, 'text.usetex': True})"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sz = [15,15]\n",
      "tM = Mesh.TensorMesh(sz)\n",
      "qM = Mesh.TreeMesh(sz)\n",
      "qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.45 else 0)\n",
      "rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 57
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fun = lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.45 else 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fun(np.ones(30))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 59,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def DCfun(mesh, pts):\n",
      "    D = mesh.faceDiv\n",
      "    G = D.T\n",
      "    sigma = 1e-2*np.ones(mesh.nC)\n",
      "    Msigi = mesh.getFaceInnerProduct(1./sigma)\n",
      "    MsigI = Utils.sdInv(Msigi)\n",
      "    A = D*MsigI*G\n",
      "    A[-1,-1] /= mesh.vol[-1] # Remove null space\n",
      "    rhs = np.zeros(mesh.nC)\n",
      "    txind = Utils.meshutils.closestPoints(mesh, pts)\n",
      "    rhs[txind] = np.r_[1,-1]\n",
      "    return A, rhs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tM.vectorCCy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 61,
       "text": [
        "array([ 0.03333333,  0.1       ,  0.16666667,  0.23333333,  0.3       ,\n",
        "        0.36666667,  0.43333333,  0.5       ,  0.56666667,  0.63333333,\n",
        "        0.7       ,  0.76666667,  0.83333333,  0.9       ,  0.96666667])"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "AtM, rhstM = DCfun(tM, pts)\n",
      "AinvtM = SolverLU(AtM)\n",
      "phitM = AinvtM*rhstM"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "AqM, rhsqM = DCfun(qM, pts)\n",
      "AinvqM = SolverLU(AqM)\n",
      "phiqM = AinvqM*rhsqM"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ArM, rhsrM = DCfun(rM, pts)\n",
      "AinvrM = SolverLU(ArM)\n",
      "phirM = AinvrM*rhsrM"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coreind = (qM.gridCC[:,0]-0.5)**2+(qM.gridCC[:,1]-0.5)**2 >0.43**2\n",
      "phiqM[coreind] = 0."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')\n",
      "Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')\n",
      "PHItM = griddata(tM.gridCC[:,0], tM.gridCC[:,1], phitM, Xi, Yi, interp='linear')\n",
      "PHIqM = griddata(qM.gridCC[:,0], qM.gridCC[:,1], phiqM, Xi, Yi, interp='linear')\n",
      "PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(1,3,figsize=(14*1.2,4*1.2))\n",
      "label = [\"(a)\", \"(b)\", \"(c)\"]\n",
      "opts = {}\n",
      "vmin, vmax = PHItM.min(), PHItM.max()\n",
      "dat = axes[0].contourf(Xi, Yi, PHItM, 100)\n",
      "tM.plotGrid(ax=axes[0], **opts)\n",
      "axes[0].set_title('TensorMesh')\n",
      "axes[1].contourf(Xi, Yi, PHIqM, 100)\n",
      "qM.plotGrid(ax=axes[1], **opts)\n",
      "axes[1].set_title('TreeMesh')\n",
      "axes[2].contourf(Xi, Yi, PHIrM, 100)\n",
      "rM.plotGrid(ax=axes[2], **opts)\n",
      "axes[2].set_title('CurvilinearMesh')\n",
      "for i in range(3):\n",
      "    axes[i].set_xlim(0.025, 0.975)\n",
      "    axes[i].set_ylim(0.025, 0.975)\n",
      "    axes[i].text(0., 1.0, label[i], fontsize=24)\n",
      "    if i==0:        \n",
      "        axes[i].set_ylabel(\"y\")\n",
      "    else:\n",
      "        axes[i].set_ylabel(\" \")\n",
      "    axes[i].set_xlabel(\"x\")\n",
      "plt.show()\n",
      "# fig.savefig(\"./ThreeMesh.png\", dpi=100)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}