Bump version: 0.1.2 → 0.1.3

Bug/fields object

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.1
-files = setup.py SimPEG/__init__.py
+current_version = 0.1.3
+files = setup.py SimPEG/__init__.py docs/conf.py
 

.travis.yml

@@ -4,7 +4,7 @@ python:
 
 # Setup anaconda
 before_install:
-  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.3.0-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.3.0-Linux-x86_64.sh -O miniconda.sh; fi
+  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi
   - chmod +x miniconda.sh
   - ./miniconda.sh -b
   - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
@@ -22,7 +22,7 @@ install:
 
 # Run test
 script:
-  - nosetests --with-cov --cov SimPEG --cov-config .coveragerc
+  - nosetests --with-cov --cov SimPEG --cov-config .coveragerc -v -s
 
 # Calculate coverage
 after_success:
@@ -31,3 +31,4 @@ after_success:
 notifications:
   email:
     - rowanc1@gmail.com
+    - lindseyheagy@gmail.com

LICENSE

@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2013-2014 SimPEG Developers
+Copyright (c) 2013-2015 SimPEG Developers
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

PROJECTS.rst

@@ -1,6 +1,36 @@
 - Electromagnetics (`simpegEM <http://simpegem.rtfd.org/>`_)
-- Magnetotellurics (`simpegMT <http://simpegmt.rtfd.org/>`_)
-- Electromagnetics 1D (`simpegEM1D <http://simpegem1d.rtfd.org/>`_)
-- Direct Current Resistivity (`simpegDC <http://simpeg-dc.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegem.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegem
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegem/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegem?branch=master
 - Potential Fields (`simpegPF <http://simpegpf.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegpf.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegpf
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegpf/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegpf?branch=master
 - Ground Water Flow (`simpegFLOW <http://simpegflow.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegflow.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegflow
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegflow/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegflow?branch=master
+- Direct Current Resistivity (`simpegDC <http://simpeg-dc.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegdc.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegdc
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegdc/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegdc?branch=master
+- Electromagnetics 1D (`simpegEM1D <http://simpegem1d.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegEM1D.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegEM1D
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegEM1D/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegEM1D?branch=master
+- Magnetotellurics (`simpegMT <http://simpegmt.rtfd.org/>`_)
+    .. image:: https://travis-ci.org/simpeg/simpegmt.svg?branch=master
+        :target: https://travis-ci.org/simpeg/simpegmt
+        :alt: Master Branch
+    .. image:: https://coveralls.io/repos/simpeg/simpegmt/badge.png?branch=master
+        :target: https://coveralls.io/r/simpeg/simpegmt?branch=master

README.md

@@ -1,32 +0,0 @@
-![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png)
-
-Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
-
-The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:
-
-* modular with respect to the spacial discretization, optimization routine, and geophysical problem
-* built with the inverse problem in mind
-* provides a framework for geophysical and hydrogeologic problems
-* supports 1D, 2D and 3D problems
-* designed for large-scale inversions
-
-Documentation:
-[http://simpeg.readthedocs.org/en/latest/](http://simpeg.readthedocs.org/en/latest/)
-
-Code:
-[https://github.com/simpeg/simpeg](https://github.com/simpeg/simpeg)
-
-Tests:
-[https://travis-ci.org/simpeg/simpeg](https://travis-ci.org/simpeg/simpeg)
-
-Build Status:
-[![Build Status](https://travis-ci.org/simpeg/simpeg.svg?branch=master)](https://travis-ci.org/simpeg/simpeg)
-
-Coverage Status:
-[![Coverage Status](https://coveralls.io/repos/simpeg/simpeg/badge.png?branch=master)](https://coveralls.io/r/simpeg/simpeg?branch=master)
-
-Bugs & Issues:
-[https://github.com/simpeg/simpeg/issues](https://github.com/simpeg/simpeg/issues)
-
-Code Snippets & Tutorials:
-[http://www.row1.ca/simpeg](http://www.row1.ca/simpeg)

README.rst

@@ -0,0 +1,60 @@
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+    :alt: SimPEG Logo
+
+======
+SimPEG
+======
+
+.. image:: https://img.shields.io/pypi/v/SimPEG.svg
+    :target: https://crate.io/packages/SimPEG/
+    :alt: Latest PyPI version
+
+.. image:: https://img.shields.io/pypi/dm/SimPEG.svg
+    :target: https://crate.io/packages/SimPEG/
+    :alt: Number of PyPI downloads
+
+.. image:: https://img.shields.io/badge/license-MIT-blue.svg
+    :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
+    :alt: BSD 3 clause license.
+
+.. image:: https://img.shields.io/travis/simpeg/simpeg.svg
+    :target: https://travis-ci.org/simpeg/simpeg
+    :alt: Travis CI build status
+
+.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
+    :target: https://coveralls.io/r/simpeg/simpeg?branch=master
+    :alt: Coverage status
+
+Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+
+The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:
+
+* modular with respect to the spacial discretization, optimization routine, and geophysical problem
+* built with the inverse problem in mind
+* provides a framework for geophysical and hydrogeologic problems
+* supports 1D, 2D and 3D problems
+* designed for large-scale inversions
+
+
+Website:
+http://simpeg.xyz
+
+
+Documentation:
+http://docs.simpeg.xyz
+
+
+Code:
+https://github.com/simpeg/simpeg
+
+
+Tests:
+https://travis-ci.org/simpeg/simpeg
+
+
+Bugs & Issues:
+https://github.com/simpeg/simpeg/issues
+
+
+Code Snippets & Tutorials:
+http://www.row1.ca/simpeg

SimPEG/DataMisfit.py

@@ -25,7 +25,6 @@ class BaseDataMisfit(object):
     def eval(self, m, u=None):
         """eval(m, u=None)
 
-            :param Problem,Survey forward: forward simulation
             :param numpy.array m: geophysical model
             :param numpy.array u: fields
             :rtype: float
@@ -38,7 +37,6 @@ class BaseDataMisfit(object):
     def evalDeriv(self, m, u=None):
         """evalDeriv(m, u=None)
 
-            :param Problem,Survey forward: forward simulation
             :param numpy.array m: geophysical model
             :param numpy.array u: fields
             :rtype: numpy.array
@@ -52,7 +50,6 @@ class BaseDataMisfit(object):
     def eval2Deriv(self, m, v, u=None):
         """eval2Deriv(m, v, u=None)
 
-            :param Problem,Survey forward: forward simulation
             :param numpy.array m: geophysical model
             :param numpy.array v: vector to multiply
             :param numpy.array u: fields
@@ -62,19 +59,20 @@ class BaseDataMisfit(object):
         """
         raise NotImplementedError('This method should be overwritten.')
 
-    def target(self, forward):
-        """target(forward)
+    # TODO: implement target misfit as a property, or possibly as an inversion directive.
 
-            Target for data misfit. By default this is the number of data,
-            which satisfies the Discrepancy Principle.
+    # def target(self, forward):
+    #     """target(forward)
 
-            :param Problem,Survey forward: forward simulation
-            :rtype: float
-            :return: data misfit target
+    #         Target for data misfit. By default this is the number of data,
+    #         which satisfies the Discrepancy Principle.
 
-        """
-        prob, survey = self.splitForward(forward)
-        return survey.nD
+    #         :rtype: float
+    #         :return: data misfit target
+
+    #     """
+    #     prob, survey = self.splitForward(forward)
+    #     return survey.nD
 
 
 class l2_DataMisfit(BaseDataMisfit):

SimPEG/Directives.py

@@ -88,7 +88,7 @@ class BetaEstimate_ByEig(InversionDirective):
     """BetaEstimate"""
 
     beta0 = None       #: The initial Beta (regularization parameter)
-    beta0_ratio = 0.1  #: estimateBeta0 is used with this ratio
+    beta0_ratio = 1e2  #: estimateBeta0 is used with this ratio
 
     def initialize(self):
         """
@@ -136,7 +136,7 @@ class BetaEstimate_ByEig(InversionDirective):
 class BetaSchedule(InversionDirective):
     """BetaSchedule"""
 
-    coolingFactor = 2.
+    coolingFactor = 8.
     coolingRate = 3
 
     def endIter(self):
@@ -145,6 +145,24 @@ class BetaSchedule(InversionDirective):
             self.invProb.beta /= self.coolingFactor
 
 
+
+class SaveModelEveryIteration(InversionDirective):
+    """SaveModelEveryIteration"""
+
+    @property
+    def modelName(self):
+        if getattr(self, '_modelName', None) is None:
+            from datetime import datetime
+            self._modelName = 'inversionModel-%s'%datetime.now().strftime('%Y-%m-%d')
+        return self._modelName
+    @modelName.setter
+    def modelName(self, value):
+        self._modelName = value
+
+    def endIter(self):
+        np.save('%03d-%s' % (self.opt.iter, self.modelName), self.opt.xc)
+
+
 # class UpdateReferenceModel(Parameter):
 
 #     mref0 = None

SimPEG/Examples/DCfwd.py

@@ -0,0 +1,71 @@
+from SimPEG import Mesh, Utils, np, SolverLU
+import matplotlib.pyplot as plt
+import matplotlib
+from matplotlib.mlab import griddata
+
+## 2D DC forward modeling example with Tensor and Curvilinear Meshes
+
+# Step1: Generate Tensor and Curvilinear Mesh
+sz = [40,40]
+# Tensor Mesh
+tM = Mesh.TensorMesh(sz)
+# Curvilinear Mesh
+rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
+# Step2: Direct Current (DC) operator
+def DCfun(mesh, pts):
+    D = mesh.faceDiv
+    G = D.T
+    sigma = 1e-2*np.ones(mesh.nC)
+    Msigi = mesh.getFaceInnerProduct(1./sigma)
+    MsigI = Utils.sdInv(Msigi)
+    A = D*MsigI*G
+    A[-1,-1] /= mesh.vol[-1] # Remove null space
+    rhs = np.zeros(mesh.nC)
+    txind = Utils.meshutils.closestPoints(mesh, pts)
+    rhs[txind] = np.r_[1,-1]
+    return A, rhs
+
+pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5]))
+
+#Step3: Solve DC problem (LU solver)
+AtM, rhstM = DCfun(tM, pts)
+AinvtM = SolverLU(AtM)
+phitM = AinvtM*rhstM
+
+ArM, rhsrM = DCfun(rM, pts)
+AinvrM = SolverLU(ArM)
+phirM = AinvrM*rhsrM
+
+#Step4: Making Figure
+fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+label = ["(a)", "(b)"]
+opts = {}
+vmin, vmax = phitM.min(), phitM.max()
+dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
+
+#TODO: At the moment Curvilinear Mesh do not have plotimage 
+
+Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
+Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
+PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
+axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
+
+cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
+cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
+
+tM.plotGrid(ax=axes[0], **opts)
+axes[0].set_title('TensorMesh')
+rM.plotGrid(ax=axes[1], **opts)
+axes[1].set_title('CurvilinearMesh')
+for i in range(2):
+    axes[i].set_xlim(0.025, 0.975)
+    axes[i].set_ylim(0.025, 0.975)
+    axes[i].text(0., 1.0, label[i], fontsize=20)
+    if i==0:        
+        axes[i].set_ylabel("y")
+    else:
+        axes[i].set_ylabel(" ")
+    axes[i].set_xlabel("x")
+
+plt.show()

SimPEG/Examples/Linear.py

@@ -9,8 +9,8 @@ class LinearProblem(Problem.BaseProblem):
 
     surveyPair = LinearSurvey
 
-    def __init__(self, model, G, **kwargs):
-        Problem.BaseProblem.__init__(self, model, **kwargs)
+    def __init__(self, mesh, G, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, **kwargs)
         self.G = G
 
     def fields(self, m, u=None):

SimPEG/InvProblem.py

@@ -45,6 +45,8 @@ class BaseInvProblem(object):
         self.prob, self.survey = dmisfit.prob, dmisfit.survey
         #TODO: Remove: (and make iteration printers better!)
         self.opt.parent = self
+        self.reg.parent = self
+        self.dmisfit.parent = self
 
     @Utils.callHooks('startup')
     def startup(self, m0):

SimPEG/Maps.py

@@ -12,7 +12,8 @@ class IdentityMap(object):
 
     mesh = None      #: A SimPEG Mesh
 
-    def __init__(self, mesh):
+    def __init__(self, mesh, **kwargs):
+        Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
     @property
@@ -85,10 +86,10 @@ class IdentityMap(object):
         """
         print 'Testing %s' % str(self)
         if m is None:
-            m = np.random.rand(self.nP)
+            m = abs(np.random.rand(self.nP))
         if 'plotIt' not in kwargs:
             kwargs['plotIt'] = False
-        return checkDerivative(lambda m : [self * m, self.deriv(m)], m, **kwargs)
+        return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4, **kwargs)
 
     def _assertMatchesPair(self, pair):
         assert (isinstance(self, pair) or
@@ -159,7 +160,6 @@ class ComboMap(IdentityMap):
 
 class ExpMap(IdentityMap):
     """
-
         Changes the model into the physical property.
 
         A common example of this is to invert for electrical conductivity
@@ -220,6 +220,46 @@ class ExpMap(IdentityMap):
         """
         return Utils.sdiag(np.exp(Utils.mkvc(m)))
 
+
+class LogMap(IdentityMap):
+    """
+        Changes the model into the physical property.
+
+        If \\(p\\) is the physical property and \\(m\\) is the model, then
+
+        ..math::
+
+            p = \\log(m)
+
+        and
+
+        ..math::
+
+            m = \\exp(p)
+
+        NOTE: If you have a model which is log conductivity (ie. \\(m = \\log(\\sigma)\\)),
+        you should be using an ExpMap
+
+    """
+
+    def __init__(self, mesh, **kwargs):
+        IdentityMap.__init__(self, mesh, **kwargs)
+
+    def _transform(self, m):
+        return np.log(Utils.mkvc(m))
+
+    def deriv(self, m):
+        mod = Utils.mkvc(m)
+        deriv = np.zeros(mod.shape)
+        tol = 1e-16 # zero
+        ind = np.greater_equal(np.abs(mod),tol)
+        deriv[ind] = 1.0/mod[ind]
+        return Utils.sdiag(deriv)
+
+    def inverse(self, m):
+        return np.exp(Utils.mkvc(m))
+
+
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap
 
@@ -261,6 +301,62 @@ class Vertical1DMap(IdentityMap):
                     ), shape=(repNum, 1))
         return sp.kron(sp.identity(self.nP), repVec)
 
+
+class Map2Dto3D(IdentityMap):
+    """Map2Dto3D
+
+        Given a 2D vector, this will extend to the full
+        3D model space.
+    """
+
+    normal = 'Y' #: The normal
+
+    def __init__(self, mesh, **kwargs):
+        assert mesh.dim == 3, 'Only works for a 3D Mesh'
+        IdentityMap.__init__(self, mesh, **kwargs)
+        assert self.normal in ['X','Y','Z'], 'For now, only "Y" normal is supported'
+
+    @property
+    def nP(self):
+        """Number of model properties.
+
+           The number of cells in the
+           last dimension of the mesh."""
+        if self.normal == 'Z':
+            return self.mesh.nCx * self.mesh.nCy
+        elif self.normal == 'Y':
+            return self.mesh.nCx * self.mesh.nCz
+        elif self.normal == 'X':
+            return self.mesh.nCy * self.mesh.nCz
+
+    def _transform(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: transformed model
+        """
+        m = Utils.mkvc(m)
+        if self.normal == 'Z':
+            return Utils.mkvc(m.reshape(self.mesh.vnC[[0,1]], order='F')[:,:,np.newaxis].repeat(self.mesh.nCz,axis=2))
+        elif self.normal == 'Y':
+            return Utils.mkvc(m.reshape(self.mesh.vnC[[0,2]], order='F')[:,np.newaxis,:].repeat(self.mesh.nCy,axis=1))
+        elif self.normal == 'X':
+            return Utils.mkvc(m.reshape(self.mesh.vnC[[1,2]], order='F')[np.newaxis,:,:].repeat(self.mesh.nCx,axis=0))
+
+    def deriv(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: scipy.csr_matrix
+            :return: derivative of transformed model
+        """
+        inds = self * np.arange(self.nP)
+        nC, nP = self.mesh.nC, self.nP
+        P = sp.csr_matrix(
+                    (np.ones(nC),
+                    (range(nC), inds)
+                ), shape=(nC, nP))
+        return P
+
 class Mesh2Mesh(IdentityMap):
     """
         Takes a model on one mesh are translates it to another mesh.
@@ -342,6 +438,118 @@ class ActiveCells(IdentityMap):
     def deriv(self, m):
         return self.P
 
+class ActiveCellsTopo(IdentityMap):
+    """
+        Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
+
+    """
+
+    indActive   = None #: Active Cells
+    valInactive = None #: Values of inactive Cells
+    nC          = None #: Number of cells in the full model
+
+    def __init__(self, mesh, indActive, nC=None):
+        self.mesh  = mesh
+
+        self.nC = nC or mesh.nC
+
+        if indActive.dtype is not bool:
+            z = np.zeros(self.nC,dtype=bool)
+            z[indActive] = True
+            indActive = z
+        self.indActive = indActive
+
+        self.indInactive = np.logical_not(indActive)
+        inds = np.nonzero(self.indActive)[0]
+        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
+
+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.indActive.sum()
+
+    def _transform(self, m):
+        val_temp = np.zeros(self.mesh.nC)
+        val_temp[self.indActive] = m
+        valInactive = np.zeros(self.mesh.nC)
+        #1D
+        if self.mesh.dim == 1:
+            z_temp = self.mesh.gridCC
+            val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])]
+        #2D
+        elif self.mesh.dim == 2:
+            act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            for i in range(self.mesh.nCx):
+                act_tempx = act_temp[i,:] == 1
+                val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])]
+            valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
+        #3D
+        elif self.mesh.dim == 3:
+            act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            for i in range(self.mesh.nCx*self.mesh.nCy):
+                act_tempxy = act_temp[i,:] == 1
+                val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])]
+            valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
+
+        self.valInactive = valInactive
+
+        return self.P*m + self.valInactive
+
+    def inverse(self, D):
+        return self.P.T*D
+
+    def deriv(self, m):
+        return self.P
+
+
+class Weighting(IdentityMap):
+    """
+        Model weight parameters.
+
+    """
+
+    weights     = None #: Active Cells
+    nC          = None #: Number of cells in the full model
+
+    def __init__(self, mesh, weights=None, nC=None):
+        self.mesh  = mesh
+
+        self.nC = nC or mesh.nC
+
+        if weights is None:
+            weights = np.ones(self.nC)
+
+        self.weights = np.array(weights, dtype=float)
+
+        self.P = Utils.sdiag(self.weights)
+
+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.nC
+
+    def _transform(self, m):
+        return self.P*m
+
+    def inverse(self, D):
+        Pinv = Utils.sdiag(self.weights**(-1.))
+        return Pinv*D
+
+    def deriv(self, m):
+        return self.P
+
 
 class ComplexMap(IdentityMap):
     """ComplexMap

SimPEG/Mesh/CurvilinearMesh.py

@@ -10,24 +10,24 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
     """
-    LogicallyRectMesh is a mesh class that deals with logically rectangular meshes.
+    CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
-    Example of a logically rectangular mesh:
+    Example of a curvilinear mesh:
 
     .. plot::
             :include-source:
 
             from SimPEG import Mesh, Utils
             X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.LogicallyRectMesh([X, Y])
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
     """
 
     __metaclass__ = Utils.SimPEGMetaClass
 
-    _meshType = 'LRM'
+    _meshType = 'Curv'
 
     def __init__(self, nodes):
         assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
@@ -38,7 +38,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
             assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
 
         assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
-        assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh."
+        assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
 
         BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None)
 
@@ -343,7 +343,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
 
             from SimPEG import Mesh, Utils
             X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.LogicallyRectMesh([X, Y])
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """
@@ -435,9 +435,9 @@ if __name__ == '__main__':
     dee3 = True
     if dee3:
         X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False)
-        M = LogicallyRectMesh([X, Y, Z])
+        M = CurvilinearMesh([X, Y, Z])
     else:
         X, Y = Utils.ndgrid(h1, h2, vector=False)
-        M = LogicallyRectMesh([X, Y])
+        M = CurvilinearMesh([X, Y])
 
     print M.r(M.normals, 'F', 'Fx', 'V')

SimPEG/Mesh/CylMesh.py

@@ -1,15 +1,20 @@
 import numpy as np
 import scipy.sparse as sp
 from scipy.constants import pi
-from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
+from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
 from TensorMesh import BaseTensorMesh
 from InnerProducts import InnerProducts
+from View import CylView
 
 
-class CylMesh(BaseTensorMesh, InnerProducts):
+class CylMesh(BaseTensorMesh, InnerProducts, CylView):
     """
         CylMesh is a mesh class for cylindrical problems
 
+        .. note::
+
+            for a cylindrically symmetric mesh use [hx, 1, hz]
+
         ::
 
             cs, nc, npad = 20., 30, 8
@@ -22,10 +27,15 @@ class CylMesh(BaseTensorMesh, InnerProducts):
 
     _unitDimensions = [1, 2*np.pi, 1]
 
-    def __init__(self, h, x0=None):
+    def __init__(self, h, x0=None, cartesianOrigin=None):
         BaseTensorMesh.__init__(self, h, x0)
-        assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]"
         assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
+        if self.dim == 2:
+            print 'Warning, a disk mesh has not been tested thoroughly.'
+        cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin
+        assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh."
+        self.cartesianOrigin = np.array(cartesianOrigin, dtype=float)
+
 
     @property
     def isSymmetric(self):
@@ -277,6 +287,18 @@ class CylMesh(BaseTensorMesh, InnerProducts):
                 #                                   kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
         return self._aveE2CC
 
+    @property
+    def aveE2CCV(self):
+        "Construct the averaging operator on cell edges to cell centers."
+        if getattr(self, '_aveE2CCV', None) is None:
+            # The number of cell centers in each direction
+            n = self.vnC
+            if self.isSymmetric:
+                return self.aveE2CC
+            else:
+                raise NotImplementedError('wrapping in the averaging is not yet implemented')
+        return self._aveE2CCV
+
 
     @property
     def aveF2CC(self):
@@ -293,3 +315,84 @@ class CylMesh(BaseTensorMesh, InnerProducts):
                 #                                    kron3(speye(n[2]), av(n[1]), speye(n[0])),
                 #                                    kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
         return self._aveF2CC
+
+    @property
+    def aveF2CCV(self):
+        "Construct the averaging operator on cell faces to cell centers."
+        if getattr(self, '_aveF2CCV', None) is None:
+            n = self.vnC
+            if self.isSymmetric:
+                avR = av(n[0])[:,1:]
+                avR[0,0] = 1.
+                self._aveF2CCV = sp.block_diag((sp.kron(speye(n[2]), avR),
+                                                sp.kron(av(n[2]), speye(n[0]))), format="csr")
+            else:
+                raise NotImplementedError('wrapping in the averaging is not yet implemented')
+        return self._aveF2CCV
+
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+        """
+            Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
+        """
+
+        assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
+
+
+        if locType == 'F':
+            # do this three times for each component
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            return sp.vstack((X,Y,Z))
+        if locType == 'E':
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
+            Z = spzeros(Mrect.nEz, self.nE)
+            return sp.vstack((X,Y,Z))
+
+        grid = getattr(Mrect, 'grid' + locType)
+        # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
+        theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
+        theta[theta < 0] += np.pi*2.0
+        r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
+
+        if locType in ['CC', 'N', 'Fz', 'Ez']:
+            G, proj = np.c_[r, theta, grid[:,2]], np.ones(r.size)
+        else:
+            dotMe = {
+                        'Fx': Mrect.normals[:Mrect.nFx,:],
+                        'Fy': Mrect.normals[Mrect.nFx:(Mrect.nFx+Mrect.nFy),:],
+                        'Fz': Mrect.normals[-Mrect.nFz:,:],
+                        'Ex': Mrect.tangents[:Mrect.nEx,:],
+                        'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
+                        'Ez': Mrect.tangents[-Mrect.nEz:,:],
+                    }[locType]
+            if 'F' in locType:
+                normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
+                proj = ( normals * dotMe ).sum(axis=1)
+            if 'E' in locType:
+                tangents = np.c_[-np.sin(theta), np.cos(theta), np.zeros(theta.size)]
+                proj = ( tangents * dotMe ).sum(axis=1)
+            G = np.c_[r, theta, grid[:,2]]
+
+        interpType = locType
+        if interpType == 'Fy':
+            interpType = 'Fx'
+        elif interpType == 'Ex':
+            interpType = 'Ey'
+
+        Pc2r = self.getInterpolationMat(G, interpType)
+        Proj = sdiag(proj)
+        return Proj * Pc2r
+
+
+
+if __name__ == '__main__':
+
+    from SimPEG import *
+    hx = np.r_[1,1,0.5]
+    hz = np.r_[2,1]
+    M = Mesh.CylMesh([hx, 1,hz], x0='00N')
+
+    M.plotImage(np.random.rand(M.nC), showIt=False)
+    M.plotGrid(centers=True, showIt=True)

SimPEG/Mesh/DiffOperators.py

@@ -422,11 +422,9 @@ class DiffOperators(object):
 
         def fget(self):
             if(self._edgeCurl is None):
-                assert self.dim > 2, "Edge Curl only programed for 3D."
+                assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
                 # The number of cell centers in each direction
-                n1 = self.nCx
-                n2 = self.nCy
-                n3 = self.nCz
+                n = self.vnC
 
                 # Compute lengths of cell edges
                 L = self.edge
@@ -435,26 +433,32 @@ class DiffOperators(object):
                 S = self.area
 
                 # Compute divergence operator on faces
-                d1 = ddx(n1)
-                d2 = ddx(n2)
-                d3 = ddx(n3)
+                if self.dim == 2:
 
-                D32 = kron3(d3, speye(n2), speye(n1+1))
-                D23 = kron3(speye(n3), d2, speye(n1+1))
-                D31 = kron3(d3, speye(n2+1), speye(n1))
-                D13 = kron3(speye(n3), speye(n2+1), d1)
-                D21 = kron3(speye(n3+1), d2, speye(n1))
-                D12 = kron3(speye(n3+1), speye(n2), d1)
+                    D21 = sp.kron(ddx(n[1]), speye(n[0]))
+                    D12 = sp.kron(speye(n[1]), ddx(n[0]))
+                    C = sp.hstack((-D21, D12), format="csr")
+                    self._edgeCurl = C*sdiag(1/S)
 
-                O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
-                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
-                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+                elif self.dim == 3:
 
-                C = sp.vstack((sp.hstack((O1, -D32, D23)),
-                               sp.hstack((D31, O2, -D13)),
-                               sp.hstack((-D21, D12, O3))), format="csr")
+                    D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
+                    D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
+                    D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
+                    D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0]))
+                    D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]))
+                    D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0]))
+
+                    O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
+                    O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
+                    O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+
+                    C = sp.vstack((sp.hstack((O1, -D32, D23)),
+                                   sp.hstack((D31, O2, -D13)),
+                                   sp.hstack((-D21, D12, O3))), format="csr")
+
+                    self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
 
-                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
             return self._edgeCurl
         return locals()
     _edgeCurl = None
@@ -567,35 +571,58 @@ class DiffOperators(object):
     @property
     def aveF2CC(self):
         "Construct the averaging operator on cell faces to cell centers."
-        if getattr(self, '_aveF2CC', None) is None:
-            n = self.vnC
-            if(self.dim == 1):
-                self._aveF2CC = av(n[0])
-            elif(self.dim == 2):
-                self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
-                                                 sp.kron(av(n[1]), speye(n[0]))), format="csr")
-            elif(self.dim == 3):
-                self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
-                                                   kron3(speye(n[2]), av(n[1]), speye(n[0])),
-                                                   kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
-        return self._aveF2CC
-
+        if(self.dim == 1):
+            return self.aveFx2CC
+        elif(self.dim == 2):
+            return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
+        elif(self.dim == 3):
+            return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
 
     @property
     def aveF2CCV(self):
         "Construct the averaging operator on cell faces to cell centers."
-        if getattr(self, '_aveF2CCV', None) is None:
+        if(self.dim == 1):
+            return self.aveFx2CC
+        elif(self.dim == 2):
+            return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr")
+        elif(self.dim == 3):
+            return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
+
+    @property
+    def aveFx2CC(self):
+        "Construct the averaging operator on cell faces in the x direction to cell centers."
+        if getattr(self, '_aveFx2CC', None) is None:
             n = self.vnC
             if(self.dim == 1):
-                self._aveF2CCV = av(n[0])
+                self._aveFx2CC = av(n[0])
             elif(self.dim == 2):
-                self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
-                                               sp.kron(av(n[1]), speye(n[0]))), format="csr")
+                self._aveFx2CC = sp.kron(speye(n[1]), av(n[0]))
             elif(self.dim == 3):
-                self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])),
-                                               kron3(speye(n[2]), av(n[1]), speye(n[0])),
-                                               kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
-        return self._aveF2CCV
+                self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0]))
+        return self._aveFx2CC
+
+    @property
+    def aveFy2CC(self):
+        "Construct the averaging operator on cell faces in the y direction to cell centers."
+        if self.dim < 2: return None
+        if getattr(self, '_aveFy2CC', None) is None:
+            n = self.vnC
+            if(self.dim == 2):
+                self._aveFy2CC = sp.kron(av(n[1]), speye(n[0]))
+            elif(self.dim == 3):
+                self._aveFy2CC = kron3(speye(n[2]), av(n[1]), speye(n[0]))
+        return self._aveFy2CC
+
+    @property
+    def aveFz2CC(self):
+        "Construct the averaging operator on cell faces in the z direction to cell centers."
+        if self.dim < 3: return None
+        if getattr(self, '_aveFz2CC', None) is None:
+            n = self.vnC
+            if(self.dim == 3):
+                self._aveFz2CC = kron3(av(n[2]), speye(n[1]), speye(n[0]))
+        return self._aveFz2CC
+
 
     @property
     def aveCC2F(self):
@@ -616,36 +643,60 @@ class DiffOperators(object):
     @property
     def aveE2CC(self):
         "Construct the averaging operator on cell edges to cell centers."
-        if getattr(self, '_aveE2CC', None) is None:
-            # The number of cell centers in each direction
-            n = self.vnC
-            if(self.dim == 1):
-                self._aveE2CC = speye(n[0])
-            elif(self.dim == 2):
-                self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
-                                               sp.kron(speye(n[1]), av(n[0]))), format="csr")
-            elif(self.dim == 3):
-                self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
-                                                  kron3(av(n[2]), speye(n[1]), av(n[0])),
-                                                  kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
-        return self._aveE2CC
+        if(self.dim == 1):
+            return self.aveEx2CC
+        elif(self.dim == 2):
+            return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
+        elif(self.dim == 3):
+            return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
 
     @property
     def aveE2CCV(self):
         "Construct the averaging operator on cell edges to cell centers."
-        if getattr(self, '_aveE2CCV', None) is None:
+        if(self.dim == 1):
+            return self.aveEx2CC
+        elif(self.dim == 2):
+            return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr")
+        elif(self.dim == 3):
+            return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
+
+    @property
+    def aveEx2CC(self):
+        "Construct the averaging operator on cell edges in the x direction to cell centers."
+        if getattr(self, '_aveEx2CC', None) is None:
             # The number of cell centers in each direction
             n = self.vnC
             if(self.dim == 1):
-                raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
+                self._aveEx2CC = speye(n[0])
             elif(self.dim == 2):
-                self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])),
-                                                sp.kron(speye(n[1]), av(n[0]))), format="csr")
+                self._aveEx2CC = sp.kron(av(n[1]), speye(n[0]))
             elif(self.dim == 3):
-                self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])),
-                                                kron3(av(n[2]), speye(n[1]), av(n[0])),
-                                                kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
-        return self._aveE2CCV
+                self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0]))
+        return self._aveEx2CC
+
+    @property
+    def aveEy2CC(self):
+        "Construct the averaging operator on cell edges in the y direction to cell centers."
+        if self.dim < 2: return None
+        if getattr(self, '_aveEy2CC', None) is None:
+            # The number of cell centers in each direction
+            n = self.vnC
+            if(self.dim == 2):
+                self._aveEy2CC = sp.kron(speye(n[1]), av(n[0]))
+            elif(self.dim == 3):
+                self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0]))
+        return self._aveEy2CC
+
+    @property
+    def aveEz2CC(self):
+        "Construct the averaging operator on cell edges in the z direction to cell centers."
+        if self.dim < 3: return None
+        if getattr(self, '_aveEz2CC', None) is None:
+            # The number of cell centers in each direction
+            n = self.vnC
+            if(self.dim == 3):
+                self._aveEz2CC = kron3(speye(n[2]), av(n[1]), av(n[0]))
+        return self._aveEz2CC
 
     @property
     def aveN2CC(self):

SimPEG/Mesh/InnerProducts.py

@@ -328,7 +328,7 @@ class InnerProducts(object):
         iijj = ndgrid(i, j)
         ii, jj = iijj[:, 0], iijj[:, 1]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             fN1 = M.r(M.normals, 'F', 'Fx', 'M')
             fN2 = M.r(M.normals, 'F', 'Fy', 'M')
 
@@ -353,7 +353,7 @@ class InnerProducts(object):
 
             PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
                                            getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
                 PXX = I2x2 * PXX
@@ -376,7 +376,7 @@ class InnerProducts(object):
         iijjkk = ndgrid(i, j, k)
         ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             fN1 = M.r(M.normals, 'F', 'Fx', 'M')
             fN2 = M.r(M.normals, 'F', 'Fy', 'M')
             fN3 = M.r(M.normals, 'F', 'Fz', 'M')
@@ -410,7 +410,7 @@ class InnerProducts(object):
 
             PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
                                            getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
                                            getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
@@ -432,7 +432,7 @@ class InnerProducts(object):
         iijj = ndgrid(i, j)
         ii, jj = iijj[:, 0], iijj[:, 1]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
             eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
 
@@ -452,7 +452,7 @@ class InnerProducts(object):
 
             PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
                                            getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
                 PXX = I2x2 * PXX
@@ -466,7 +466,7 @@ class InnerProducts(object):
         iijjkk = ndgrid(i, j, k)
         ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
             eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
             eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
@@ -495,7 +495,7 @@ class InnerProducts(object):
 
             PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
                                            getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
                                            getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))

SimPEG/Mesh/TensorMesh.py

@@ -13,7 +13,7 @@ class BaseTensorMesh(BaseRectangularMesh):
     _unitDimensions = [1, 1, 1]
 
     def __init__(self, h_in, x0_in=None):
-        assert type(h_in) is list, 'h_in must be a list'
+        assert type(h_in) in [list, tuple], 'h_in must be a list'
         assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3'
         h = range(len(h_in))
         for i, h_i in enumerate(h_in):
@@ -168,21 +168,21 @@ class BaseTensorMesh(BaseRectangularMesh):
 
         """
 
-        if   key is 'Fx':
+        if   key == 'Fx':
             ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
-        elif key is 'Fy':
+        elif key == 'Fy':
             ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
-        elif key is 'Fz':
+        elif key == 'Fz':
             ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
-        elif key is 'Ex':
+        elif key == 'Ex':
             ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
-        elif key is 'Ey':
+        elif key == 'Ey':
             ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
-        elif key is 'Ez':
+        elif key == 'Ez':
             ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
-        elif key is 'CC':
+        elif key == 'CC':
             ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
-        elif key is 'N':
+        elif key == 'N':
             ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
 
         return [t for t in ten if t is not None]
@@ -204,10 +204,12 @@ class BaseTensorMesh(BaseRectangularMesh):
         if locType == 'N' and self._meshType == 'CYL':
             #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction!
             tensors[0] = np.r_[0.,tensors[0]]
+            tensors[1] = np.r_[tensors[1], 2.0*np.pi]
 
         inside = np.ones(pts.shape[0],dtype=bool)
         for i, tensor in enumerate(tensors):
-            inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max())
+            TOL = np.diff(tensor).min() * 1.0e-10
+            inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL)
         return inside
 
     def getInterpolationMat(self, loc, locType, zerosOutside=False):

SimPEG/Mesh/View.py

@@ -109,6 +109,10 @@ class TensorView(object):
                 vc = (self.aveN2CC*v).reshape(self.vnC, order='F')
             elif vType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
                 aveOp = 'ave' + vType[0] + '2CCV'
+                # n = getattr(self,'vn'+vType[0])
+                # if 'x' in vType: v = np.r_[v,np.zeros(n[1]),np.zeros(n[2])]
+                # if 'y' in vType: v = np.r_[np.zeros(n[0]),v,np.zeros(n[2])]
+                # if 'z' in vType: v = np.r_[np.zeros(n[0]),np.zeros(n[1]),v]
                 v = getattr(self,aveOp)*v # average to cell centers
                 ind_xyz = {'x':0,'y':1,'z':2}[vType[1]]
                 vc = self.r(v.reshape((self.nC,-1),order='F'), 'CC','CC','M')[ind_xyz]
@@ -190,9 +194,16 @@ class TensorView(object):
             M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
 
         """
+        if type(vType) in [list, tuple]:
+            assert ax is None, "cannot specify an axis to plot on with this function."
+            fig, axs = plt.subplots(1,len(vType))
+            out = []
+            for vTypeI, ax in zip(vType, axs):
+                out += [self.plotSlice(v,vType=vTypeI, normal=normal, ind=ind, grid=grid, view=view, ax=ax, clim=clim, showIt=False, pcolorOpts=pcolorOpts, streamOpts=streamOpts, gridOpts=gridOpts)]
+            return out
         viewOpts = ['real','imag','abs','vec']
         normalOpts = ['X', 'Y', 'Z']
-        vTypeOpts = ['CC', 'CCv','F','E']
+        vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
 
         # Some user error checking
         assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
@@ -219,9 +230,15 @@ class TensorView(object):
             elif vType == 'CCv':
                 assert view == 'vec', 'Other types for CCv not supported'
             else:
-                # Now just deal with 'F' and 'E'
+                # Now just deal with 'F' and 'E' (x,y,z, maybe...)
                 aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
-                v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
+                Av = getattr(self,aveOp)
+                if v.size == Av.shape[1]:
+                    v = Av * v
+                else:
+                    v = self.r(v,vType[0],vType) # get specific component
+                    v = Av * v
+                # we should now be averaged to cell centers (might be a vector)
             v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
             if view == 'vec':
                 outSlice = []
@@ -518,13 +535,52 @@ class TensorView(object):
 
         return animate(fig, animateFrame, frames=len(frames))
 
+class CylView(object):
+
+    def _plotCylTensorMesh(self, plotType, *args, **kwargs):
+
+        if not self.isSymmetric:
+            raise Exception('We have not yet implemented this type of view.')
+        assert plotType in ['plotImage', 'plotGrid']
+        # Hackity Hack:
+        # Just create a TM and use its view.
+        from SimPEG.Mesh import TensorMesh
+        M = TensorMesh([self.hx, self.hz], x0=[self.x0[0], self.x0[2]])
+
+        ax = kwargs.get('ax', None)
+        if ax is None:
+            fig = plt.figure()
+            ax = plt.subplot(111)
+            kwargs['ax'] = ax
+        else:
+            assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
+            fig = ax.figure
+
+        # Don't show things in the TM.plotImage
+        showIt = kwargs.get('showIt', False)
+        kwargs['showIt'] = False
+
+        out = getattr(M, plotType)(*args, **kwargs)
+
+        ax.set_xlabel('x')
+        ax.set_ylabel('z')
+
+        if showIt: plt.show()
+
+        return out
 
 
-class LomView(object):
+    def plotGrid(self, *args, **kwargs):
+        return self._plotCylTensorMesh('plotGrid', *args, **kwargs)
+
+    def plotImage(self, *args, **kwargs):
+        return self._plotCylTensorMesh('plotImage', *args, **kwargs)
+
+class CurvView(object):
     """
-    Provides viewing functions for LogicallyOrthogonalMesh
+    Provides viewing functions for CurvilinearMesh
 
-    This class is inherited by LogicallyOrthogonalMesh
+    This class is inherited by CurvilinearMesh
 
     """
     def __init__(self):
@@ -538,8 +594,8 @@ class LomView(object):
             :include-source:
 
             from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLomGird([3,3],'rotate')
-            M = Mesh.LogicallyOrthogonalMesh([X, Y])
+            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """

SimPEG/Mesh/__init__.py

@@ -1,5 +1,5 @@
 from TensorMesh import TensorMesh
 from CylMesh import CylMesh
-from LogicallyRectMesh import LogicallyRectMesh
+from CurvilinearMesh import CurvilinearMesh
 from TreeMesh import TreeMesh
 from BaseMesh import BaseMesh

SimPEG/Optimization.py

@@ -444,6 +444,7 @@ class Minimize(object):
         else: # Assume inversion is the parent
             group.attrs['phi_d'] = self.parent.phi_d
             group.attrs['phi_m'] = self.parent.phi_m
+            group.attrs['beta'] = self.parent.beta
             group.setArray('m', self.xc)
             group.setArray('dpred', self.parent.dpred)
 
@@ -929,7 +930,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
 
     @approxHinv.setter
     def approxHinv(self, value):
-        self._approxHinv = value        
+        self._approxHinv = value
 
     @Utils.timeIt
     def findSearchDirection(self):
@@ -949,7 +950,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
             return p
         else:
 
-            
+
             delx = np.zeros(self.g.size)
             resid = -(1-Active) * self.g
 

SimPEG/Problem.py

@@ -7,7 +7,7 @@ class Fields(object):
     """Fancy Field Storage
 
         u[:,'phi'] = phi
-        print u[tx0,'phi']
+        print u[src0,'phi']
 
     """
 
@@ -43,14 +43,14 @@ class Fields(object):
         return "%e MB"%sz
 
     def _storageShape(self, loc):
-        nTx = self.survey.nTx
+        nSrc = self.survey.nSrc
 
         nP = {'CC': self.mesh.nC,
               'N':  self.mesh.nN,
               'F':  self.mesh.nF,
               'E':  self.mesh.nE}[loc]
 
-        return (nP, nTx)
+        return (nP, nSrc)
 
     def _initStore(self, name):
         if name in self._fields:
@@ -70,17 +70,17 @@ class Fields(object):
 
         return field
 
-    def _txIndex(self, txTestList):
-        if type(txTestList) is slice:
-            ind = txTestList
+    def _srcIndex(self, srcTestList):
+        if type(srcTestList) is slice:
+            ind = srcTestList
         else:
-            if type(txTestList) is not list:
-                txTestList = [txTestList]
-            for txTest in txTestList:
-                if txTest not in self.survey.txList:
-                    raise KeyError('Invalid Transmitter, not in survey list.')
+            if type(srcTestList) is not list:
+                srcTestList = [srcTestList]
+            for srcTest in srcTestList:
+                if srcTest not in self.survey.srcList:
+                    raise KeyError('Invalid Source, not in survey list.')
 
-            ind = np.in1d(self.survey.txList, txTestList)
+            ind = np.in1d(self.survey.srcList, srcTestList)
         return ind
 
     def _nameIndex(self, name, accessType):
@@ -107,11 +107,11 @@ class Fields(object):
         if len(key) == 1:
             key += (None,)
 
-        assert len(key) == 2, 'must be [Tx, fieldName]'
+        assert len(key) == 2, 'must be [Src, fieldName]'
 
-        txTestList, name = key
+        srcTestList, name = key
         name = self._nameIndex(name, accessType)
-        ind = self._txIndex(txTestList)
+        ind = self._srcIndex(srcTestList)
         return ind, name
 
     def __setitem__(self, key, value):
@@ -139,7 +139,7 @@ class Fields(object):
         return self._getField(name, ind)
 
     def _setField(self, field, val, name, ind):
-        if isinstance(val, np.ndarray) and (field.shape[1] == 1 or val.ndim == 1):
+        if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1):
             val = Utils.mkvc(val,2)
         field[:,ind] = val
 
@@ -150,19 +150,18 @@ class Fields(object):
             # Aliased fields
             alias, loc, func = self.aliasFields[name]
 
-            txII   = np.array(self.survey.txList)[ind]
-            if isinstance(txII, np.ndarray):
-                txII = txII.tolist()
-            if len(txII) == 1:
-                txII = txII[0]
+            srcII   = np.array(self.survey.srcList)[ind]
+            if isinstance(srcII, np.ndarray):
+                srcII = srcII.tolist()
+            if len(srcII) == 1:
+                srcII = srcII[0]
 
             if type(func) is str:
                 assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
                 func = getattr(self, func)
-            out = func(self._fields[alias][:,ind], txII)
-
-        if out.shape[1] == 1:
-            out = Utils.mkvc(out)
+            out = func(self._fields[alias][:,ind], srcII)
+        if isinstance(out, np.ndarray) and (out.shape[0] == out.size or out.ndim == 1):
+            out = Utils.mkvc(out,2)
         return out
 
     def __contains__(self, other):
@@ -175,7 +174,7 @@ class TimeFields(Fields):
     """Fancy Field Storage for time domain problems
 
         u[:,'phi', timeInd] = phi
-        print u[tx0,'phi']
+        print u[src0,'phi']
 
     """
 
@@ -184,9 +183,9 @@ class TimeFields(Fields):
               'N':  self.mesh.nN,
               'F':  self.mesh.nF,
               'E':  self.mesh.nE}[loc]
-        nTx = self.survey.nTx
+        nSrc = self.survey.nSrc
         nT = self.survey.prob.nT + 1
-        return (nP, nTx, nT)
+        return (nP, nSrc, nT)
 
     def _indexAndNameFromKey(self, key, accessType):
         if type(key) is not tuple:
@@ -196,74 +195,76 @@ class TimeFields(Fields):
         if len(key) == 2:
             key += (slice(None,None,None),)
 
-        assert len(key) == 3, 'must be [Tx, fieldName, times]'
+        assert len(key) == 3, 'must be [Src, fieldName, times]'
 
-        txTestList, name, timeInd = key
+        srcTestList, name, timeInd = key
 
         name = self._nameIndex(name, accessType)
-        txInd = self._txIndex(txTestList)
+        srcInd = self._srcIndex(srcTestList)
 
-        return (txInd, timeInd), name
+        return (srcInd, timeInd), name
 
     def _correctShape(self, name, ind, deflate=False):
-        txInd, timeInd = ind
+        srcInd, timeInd = ind
         if name in self.knownFields:
             loc = self.knownFields[name]
         else:
             loc = self.aliasFields[name][1]
-        nP, total_nTx, total_nT = self._storageShape(loc)
-        nTx = np.ones(total_nTx, dtype=bool)[txInd].sum()
+        nP, total_nSrc, total_nT = self._storageShape(loc)
+        nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
         nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
-        shape = nP, nTx, nT
+        shape = nP, nSrc, nT
         if deflate:
              shape = tuple([s for s in shape if s > 1])
+        if len(shape) == 1:
+            shape = shape + (1,)
         return shape
 
     def _setField(self, field, val, name, ind):
-        txInd, timeInd = ind
+        srcInd, timeInd = ind
         shape = self._correctShape(name, ind)
         if Utils.isScalar(val):
-            field[:,txInd,timeInd] = val
+            field[:,srcInd,timeInd] = val
             return
         if val.size != np.array(shape).prod():
             raise ValueError('Incorrect size for data.')
-        correctShape = field[:,txInd,timeInd].shape
-        field[:,txInd,timeInd] = val.reshape(correctShape, order='F')
+        correctShape = field[:,srcInd,timeInd].shape
+        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
 
     def _getField(self, name, ind):
-        txInd, timeInd = ind
+        srcInd, timeInd = ind
 
         if name in self._fields:
-            out = self._fields[name][:,txInd,timeInd]
+            out = self._fields[name][:,srcInd,timeInd]
         else:
             # Aliased fields
             alias, loc, func = self.aliasFields[name]
             if type(func) is str:
                 assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
                 func = getattr(self, func)
-            pointerFields = self._fields[alias][:,txInd,timeInd]
+            pointerFields = self._fields[alias][:,srcInd,timeInd]
             pointerShape = self._correctShape(alias, ind)
             pointerFields = pointerFields.reshape(pointerShape, order='F')
 
             timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
-            txII   = np.array(self.survey.txList)[txInd]
-            if isinstance(txII, np.ndarray):
-                txII = txII.tolist()
-            if len(txII) == 1:
-                txII = txII[0]
+            srcII   = np.array(self.survey.srcList)[srcInd]
+            if isinstance(srcII, np.ndarray):
+                srcII = srcII.tolist()
+            if len(srcII) == 1:
+                srcII = srcII[0]
 
             if timeII.size == 1:
                 pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
                 pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
-                out = func(pointerFields, txII, timeII)
+                out = func(pointerFields, srcII, timeII)
             else: #loop over the time steps
                 nT = pointerShape[2]
                 out = range(nT)
                 for i, TIND_i in enumerate(timeII):
                     fieldI = pointerFields[:,:,i]
-                    if fieldI.ndim == 2 and fieldI.shape[1] == 1:
-                        fieldI = Utils.mkvc(fieldI)
-                    out[i] = func(fieldI, txII, TIND_i)
+                    if fieldI.shape[0] == fieldI.size:
+                        fieldI = Utils.mkvc(fieldI,2)
+                    out[i] = func(fieldI, srcII, TIND_i)
                     if out[i].ndim == 1:
                         out[i] = out[i][:,np.newaxis,np.newaxis]
                     elif out[i].ndim == 2:
@@ -321,6 +322,7 @@ class BaseProblem(object):
         self.survey._prob = None
         self._survey = None
 
+
     deleteTheseOnModelUpdate = [] # List of strings, e.g. ['_MeSigma', '_MeSigmaI']
 
     @property
@@ -437,19 +439,7 @@ class BaseTimeProblem(BaseProblem):
             del self.timeMesh
             return
 
-        if type(value) is not list:
-            raise Exception('timeSteps must be a np.ndarray or a list of scalars and tuples.')
-
-        proposed = []
-        for v in value:
-            if Utils.isScalar(v):
-                proposed += [float(v)]
-            elif type(v) is tuple and len(v) == 2:
-                proposed += [float(v[0])]*int(v[1])
-            else:
-                raise Exception('timeSteps list must contain only scalars and len(2) tuples.')
-
-        self._timeSteps = np.array(proposed)
+        self._timeSteps = Utils.meshTensor(value)
         del self.timeMesh
 
     @property

SimPEG/Regularization.py

@@ -40,7 +40,7 @@ class BaseRegularization(object):
     @property
     def inv(self): return self.parent.inv
     @property
-    def objFunc(self): return self.parent
+    def invProb(self): return self.parent
     @property
     def reg(self): return self
     @property
@@ -81,7 +81,7 @@ class BaseRegularization(object):
         """
         mD = self.mapping.deriv(m - self.mref)
         r = self.W * ( self.mapping * (m - self.mref) )
-        return mD.T * ( self.W.T * r )
+        return  mD.T * ( self.W.T * r )
 
     @Utils.timeIt
     def eval2Deriv(self, m, v=None):
@@ -196,7 +196,7 @@ class Tikhonov(BaseRegularization):
 
 
     """
-
+    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
     alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
@@ -212,7 +212,7 @@ class Tikhonov(BaseRegularization):
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self,'_Ws', None) is None:
-            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+                self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
         return self._Ws
 
     @property
@@ -272,3 +272,45 @@ class Tikhonov(BaseRegularization):
             self._W = sp.vstack(wlist)
         return self._W
 
+    @Utils.timeIt
+    def eval(self, m):
+        if self.smoothModel == True:
+            r1 = self.W * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.smoothModel == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        if self.smoothModel == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(self.mref)
+            r1 = self.W * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (self.mref) )
+            out1 = mD1.T * ( self.W.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1-out2
+        elif self.smoothModel == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
+

SimPEG/Survey.py

@@ -6,7 +6,7 @@ class BaseRx(object):
 
     locs = None   #: Locations (nRx x nDim)
 
-    knownRxTypes = None  #: Set this to a list of strings to ensure that txType is known
+    knownRxTypes = None  #: Set this to a list of strings to ensure that srcType is known
 
     projGLoc = 'CC'  #: Projection grid location, default is CC
 
@@ -111,37 +111,23 @@ class BaseTimeRx(BaseRx):
         return P
 
 
-class BaseTx(object):
-    """SimPEG Transmitter Object"""
+class BaseSrc(object):
+    """SimPEG Source Object"""
 
     loc    = None #: Location [x,y,z]
 
     rxList = None #: SimPEG Receiver List
     rxPair = BaseRx
 
-    knownTxTypes = None #: Set this to a list of strings to ensure that txType is known
-
-    def __init__(self, loc, txType, rxList, **kwargs):
+    def __init__(self, rxList, **kwargs):
         assert type(rxList) is list, 'rxList must be a list'
         for rx in rxList:
             assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
         assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
 
-        self.loc    = loc
-        self.txType = txType
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
-    @property
-    def txType(self):
-        """Transmitter Type"""
-        return getattr(self, '_txType', None)
-    @txType.setter
-    def txType(self, value):
-        known = self.knownTxTypes
-        if known is not None:
-            assert value in known, "txType must be in ['%s']" % ("', '".join(known))
-        self._txType = value
 
     @property
     def nD(self):
@@ -155,59 +141,59 @@ class BaseTx(object):
 
 
 class Data(object):
-    """Fancy data storage by Tx and Rx"""
+    """Fancy data storage by Src and Rx"""
 
     def __init__(self, survey, v=None):
         self.survey = survey
         self._dataDict = {}
-        for tx in self.survey.txList:
-            self._dataDict[tx] = {}
+        for src in self.survey.srcList:
+            self._dataDict[src] = {}
         if v is not None:
             self.fromvec(v)
 
     def _ensureCorrectKey(self, key):
         if type(key) is tuple:
             if len(key) is not 2:
-                raise KeyError('Key must be [Tx, Rx]')
-            if key[0] not in self.survey.txList:
-                raise KeyError('Tx Key must be a transmitter in the survey.')
+                raise KeyError('Key must be [Src, Rx]')
+            if key[0] not in self.survey.srcList:
+                raise KeyError('Src Key must be a source in the survey.')
             if key[1] not in key[0].rxList:
-                raise KeyError('Rx Key must be a receiver for the transmitter.')
+                raise KeyError('Rx Key must be a receiver for the source.')
             return key
-        elif isinstance(key, self.survey.txPair):
-            if key not in self.survey.txList:
-                raise KeyError('Key must be a transmitter in the survey.')
+        elif isinstance(key, self.survey.srcPair):
+            if key not in self.survey.srcList:
+                raise KeyError('Key must be a source in the survey.')
             return key, None
         else:
-            raise KeyError('Key must be [Tx] or [Tx,Rx]')
+            raise KeyError('Key must be [Src] or [Src,Rx]')
 
     def __setitem__(self, key, value):
-        tx, rx = self._ensureCorrectKey(key)
-        assert rx is not None, 'set data using [Tx, Rx]'
+        src, rx = self._ensureCorrectKey(key)
+        assert rx is not None, 'set data using [Src, Rx]'
         assert isinstance(value, np.ndarray), 'value must by ndarray'
-        assert value.size == rx.nD, "value must have the same number of data as the transmitter."
-        self._dataDict[tx][rx] = Utils.mkvc(value)
+        assert value.size == rx.nD, "value must have the same number of data as the source."
+        self._dataDict[src][rx] = Utils.mkvc(value)
 
     def __getitem__(self, key):
-        tx, rx = self._ensureCorrectKey(key)
+        src, rx = self._ensureCorrectKey(key)
         if rx is not None:
-            if rx not in self._dataDict[tx]:
+            if rx not in self._dataDict[src]:
                 raise Exception('Data for receiver has not yet been set.')
-            return self._dataDict[tx][rx]
+            return self._dataDict[src][rx]
 
-        return np.concatenate([self[tx,rx] for rx in tx.rxList])
+        return np.concatenate([self[src,rx] for rx in src.rxList])
 
     def tovec(self):
-        return np.concatenate([self[tx] for tx in self.survey.txList])
+        return np.concatenate([self[src] for src in self.survey.srcList])
 
     def fromvec(self, v):
         v = Utils.mkvc(v)
         assert v.size == self.survey.nD, 'v must have the correct number of data.'
         indBot, indTop = 0, 0
-        for tx in self.survey.txList:
-            for rx in tx.rxList:
+        for src in self.survey.srcList:
+            for rx in src.rxList:
                 indTop += rx.nD
-                self[tx, rx] = v[indBot:indTop]
+                self[src, rx] = v[indBot:indTop]
                 indBot += rx.nD
 
 
@@ -226,19 +212,19 @@ class BaseSurvey(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    txPair = BaseTx  #: Transmitter Pair
+    srcPair = BaseSrc  #: Source Pair
 
     @property
-    def txList(self):
-        """Transmitter List"""
-        return getattr(self, '_txList', None)
+    def srcList(self):
+        """Source List"""
+        return getattr(self, '_srcList', None)
 
-    @txList.setter
-    def txList(self, value):
-        assert type(value) is list, 'txList must be a list'
-        assert np.all([isinstance(tx, self.txPair) for tx in value]), 'All transmitters must be instances of %s' % self.txPair.__name__
-        assert len(set(value)) == len(value), 'The txList must be unique'
-        self._txList = value
+    @srcList.setter
+    def srcList(self, value):
+        assert type(value) is list, 'srcList must be a list'
+        assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
+        assert len(set(value)) == len(value), 'The srcList must be unique'
+        self._srcList = value
 
     @property
     def prob(self):
@@ -282,12 +268,12 @@ class BaseSurvey(object):
     @property
     def vnD(self):
         """Vector number of data"""
-        return np.array([tx.nD for tx in self.txList])
+        return np.array([src.nD for src in self.srcList])
 
     @property
-    def nTx(self):
-        """Number of Transmitters"""
-        return len(self.txList)
+    def nSrc(self):
+        """Number of Sources"""
+        return len(self.srcList)
 
     @Utils.count
     @Utils.requires('prob')
@@ -355,16 +341,17 @@ class BaseSurvey(object):
         "Check if the data is synthetic."
         return self.mtrue is not None
 
-    def makeSyntheticData(self, m, std=0.05, u=None):
+    def makeSyntheticData(self, m, std=0.05, u=None, force=False):
         """
             Make synthetic data given a model, and a standard deviation.
 
             :param numpy.array m: geophysical model
             :param numpy.array std: standard deviation
             :param numpy.array u: fields for the given model (if pre-calculated)
+            :param bool force: force overwriting of dobs
 
         """
-        if getattr(self, 'dobs', None) is not None:
+        if getattr(self, 'dobs', None) is not None and not force:
             raise Exception('Survey already has dobs.')
         self.mtrue = m
         self.dtrue = self.dpred(m, u=u)

SimPEG/Tests/TestUtils.py

@@ -1,9 +1,9 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from numpy.linalg import norm
-from SimPEG.Utils import mkvc, sdiag
+from SimPEG.Utils import mkvc, sdiag, diagEst
 from SimPEG import Utils
-from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
+from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
 import numpy as np
 import scipy.sparse as sp
 import unittest
@@ -115,7 +115,7 @@ class OrderTest(unittest.TestCase):
                 max_h = max([np.max(hi) for hi in self.M.h])
             return max_h
 
-        elif 'LRM' in self._meshType:
+        elif 'Curv' in self._meshType:
             if 'uniform' in self._meshType:
                 kwrd = 'rect'
             elif 'rotate' in self._meshType:
@@ -126,10 +126,10 @@ class OrderTest(unittest.TestCase):
                 raise Exception('Lom not supported for 1D')
             elif self.meshDimension == 2:
                 X, Y = Utils.exampleLrmGrid([nc, nc], kwrd)
-                self.M = LogicallyRectMesh([X, Y])
+                self.M = CurvilinearMesh([X, Y])
             elif self.meshDimension == 3:
                 X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
-                self.M = LogicallyRectMesh([X, Y, Z])
+                self.M = CurvilinearMesh([X, Y, Z])
             return 1./nc
 
     def getError(self):
@@ -267,7 +267,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
     # Ensure we are about precision
     order0 = order0[E0[1:] > eps]
     order1 = order1[E1[1:] > eps]
-    belowTol = order1.size == 0 and order0.size > 0
+    belowTol = (order1.size == 0 and order0.size >= 0)
     # Make sure we get the correct order
     correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder
 

SimPEG/Tests/test_CurvilinearMesh.py

@@ -0,0 +1,104 @@
+import numpy as np
+import unittest
+from SimPEG.Mesh import TensorMesh, CurvilinearMesh
+from SimPEG.Utils import ndgrid
+
+
+class BasicCurvTests(unittest.TestCase):
+
+    def setUp(self):
+        a = np.array([1, 1, 1])
+        b = np.array([1, 2])
+        c = np.array([1, 4])
+        gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
+        X, Y = ndgrid(gridIt([a, b]), vector=False)
+        self.TM2 = TensorMesh([a, b])
+        self.Curv2 = CurvilinearMesh([X, Y])
+        X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
+        self.TM3 = TensorMesh([a, b, c])
+        self.Curv3 = CurvilinearMesh([X, Y, Z])
+
+    def test_area_3D(self):
+        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
+        self.assertTrue(np.all(self.Curv3.area == test_area))
+
+    def test_vol_3D(self):
+        test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
+        np.testing.assert_almost_equal(self.Curv3.vol, test_vol)
+        self.assertTrue(True)  # Pass if you get past the assertion.
+
+    def test_vol_2D(self):
+        test_vol = np.array([1, 1, 1, 2, 2, 2])
+        t1 = np.all(self.Curv2.vol == test_vol)
+        self.assertTrue(t1)
+
+    def test_edge_3D(self):
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
+        t1 = np.all(self.Curv3.edge == test_edge)
+        self.assertTrue(t1)
+
+    def test_edge_2D(self):
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
+        t1 = np.all(self.Curv2.edge == test_edge)
+        self.assertTrue(t1)
+
+    def test_tangents(self):
+        T = self.Curv2.tangents
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy)))
+
+        T = self.Curv3.tangents
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx)))
+
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy)))
+
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz)))
+
+    def test_normals(self):
+        N = self.Curv2.normals
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy)))
+
+        N = self.Curv3.normals
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx)))
+
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy)))
+
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz)))
+
+    def test_grid(self):
+        self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
+        self.assertTrue(np.all(self.Curv2.gridN == self.TM2.gridN))
+        self.assertTrue(np.all(self.Curv2.gridFx == self.TM2.gridFx))
+        self.assertTrue(np.all(self.Curv2.gridFy == self.TM2.gridFy))
+        self.assertTrue(np.all(self.Curv2.gridEx == self.TM2.gridEx))
+        self.assertTrue(np.all(self.Curv2.gridEy == self.TM2.gridEy))
+
+        self.assertTrue(np.all(self.Curv3.gridCC == self.TM3.gridCC))
+        self.assertTrue(np.all(self.Curv3.gridN == self.TM3.gridN))
+        self.assertTrue(np.all(self.Curv3.gridFx == self.TM3.gridFx))
+        self.assertTrue(np.all(self.Curv3.gridFy == self.TM3.gridFy))
+        self.assertTrue(np.all(self.Curv3.gridFz == self.TM3.gridFz))
+        self.assertTrue(np.all(self.Curv3.gridEx == self.TM3.gridEx))
+        self.assertTrue(np.all(self.Curv3.gridEy == self.TM3.gridEy))
+        self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz))
+
+
+if __name__ == '__main__':
+    unittest.main()

SimPEG/Tests/test_LogicallyRectMesh.py

@@ -1,104 +0,0 @@
-import numpy as np
-import unittest
-from SimPEG.Mesh import TensorMesh, LogicallyRectMesh
-from SimPEG.Utils import ndgrid
-
-
-class BasicLRMTests(unittest.TestCase):
-
-    def setUp(self):
-        a = np.array([1, 1, 1])
-        b = np.array([1, 2])
-        c = np.array([1, 4])
-        gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
-        X, Y = ndgrid(gridIt([a, b]), vector=False)
-        self.TM2 = TensorMesh([a, b])
-        self.LRM2 = LogicallyRectMesh([X, Y])
-        X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
-        self.TM3 = TensorMesh([a, b, c])
-        self.LRM3 = LogicallyRectMesh([X, Y, Z])
-
-    def test_area_3D(self):
-        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
-        self.assertTrue(np.all(self.LRM3.area == test_area))
-
-    def test_vol_3D(self):
-        test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
-        np.testing.assert_almost_equal(self.LRM3.vol, test_vol)
-        self.assertTrue(True)  # Pass if you get past the assertion.
-
-    def test_vol_2D(self):
-        test_vol = np.array([1, 1, 1, 2, 2, 2])
-        t1 = np.all(self.LRM2.vol == test_vol)
-        self.assertTrue(t1)
-
-    def test_edge_3D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
-        t1 = np.all(self.LRM3.edge == test_edge)
-        self.assertTrue(t1)
-
-    def test_edge_2D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
-        t1 = np.all(self.LRM2.edge == test_edge)
-        self.assertTrue(t1)
-
-    def test_tangents(self):
-        T = self.LRM2.tangents
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy)))
-
-        T = self.LRM3.tangents
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx)))
-
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy)))
-
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz)))
-
-    def test_normals(self):
-        N = self.LRM2.normals
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy)))
-
-        N = self.LRM3.normals
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx)))
-
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy)))
-
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz)))
-
-    def test_grid(self):
-        self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC))
-        self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN))
-        self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx))
-        self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy))
-        self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx))
-        self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy))
-
-        self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC))
-        self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN))
-        self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx))
-        self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy))
-        self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz))
-        self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx))
-        self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy))
-        self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz))
-
-
-if __name__ == '__main__':
-    unittest.main()

SimPEG/Tests/test_Survey.py

@@ -7,22 +7,22 @@ class DataAndFieldsTest(unittest.TestCase):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        txLoc = np.r_[0,0,0.]
+        srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
-        Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
-        txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
-        survey = Survey.BaseSurvey(txList=txList)
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
         self.D = Survey.Data(survey)
         self.F = Problem.Fields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}, dtype={"phi":float,"e":complex,"b":complex})
-        self.Tx0 = Tx0
-        self.Tx1 = Tx1
+        self.Src0 = Src0
+        self.Src1 = Src1
         self.mesh = mesh
         self.XYZ = XYZ
 
@@ -33,12 +33,12 @@ class DataAndFieldsTest(unittest.TestCase):
 
     def test_data(self):
         V = []
-        for tx in self.D.survey.txList:
-            for rx in tx.rxList:
+        for src in self.D.survey.srcList:
+            for rx in src.rxList:
                 v = np.random.rand(rx.nD)
                 V += [v]
-                self.D[tx, rx] = v
-                self.assertTrue(np.all(v == self.D[tx, rx]))
+                self.D[src, rx] = v
+                self.assertTrue(np.all(v == self.D[src, rx]))
         V = np.concatenate(V)
         self.assertTrue(np.all(V == Utils.mkvc(self.D)))
 
@@ -47,25 +47,25 @@ class DataAndFieldsTest(unittest.TestCase):
 
     def test_contains(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         self.assertTrue('b' not in F)
         self.assertTrue('e' not in F)
-        e = np.random.rand(F.mesh.nE, nTx)
+        e = np.random.rand(F.mesh.nE, nSrc)
         F[:, 'e'] = e
         self.assertTrue('b' not in F)
         self.assertTrue('e' in F)
 
-    def test_uniqueTxs(self):
-        txs = self.D.survey.txList
-        txs += [txs[0]]
-        self.assertRaises(AssertionError, Survey.BaseSurvey, txList=txs)
+    def test_uniqueSrcs(self):
+        srcs = self.D.survey.srcList
+        srcs += [srcs[0]]
+        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
 
     def test_SetGet(self):
         F = self.F
-        nTx = F.survey.nTx
-        e = np.random.rand(F.mesh.nE, nTx) + np.random.rand(F.mesh.nE, nTx)*1j
+        nSrc = F.survey.nSrc
+        e = np.random.rand(F.mesh.nE, nSrc) + np.random.rand(F.mesh.nE, nSrc)*1j
         F[:, 'e'] = e
-        b = np.random.rand(F.mesh.nF, nTx) + np.random.rand(F.mesh.nF, nTx)*1j
+        b = np.random.rand(F.mesh.nF, nSrc) + np.random.rand(F.mesh.nF, nSrc)*1j
         F[:, 'b'] = b
 
         self.assertTrue(np.all(F[:, 'e'] == e))
@@ -79,31 +79,31 @@ class DataAndFieldsTest(unittest.TestCase):
             self.assertTrue(np.all(F[:, 'b'] == b*0))
 
         b = np.random.rand(F.mesh.nF,1)
-        F[self.Tx0, 'b'] = b
-        self.assertTrue(np.all(F[self.Tx0, 'b'] == Utils.mkvc(b)))
+        F[self.Src0, 'b'] = b
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
-        b = np.random.rand(F.mesh.nF)
-        F[self.Tx0, 'b'] = b
-        self.assertTrue(np.all(F[self.Tx0, 'b'] == b))
+        b = np.random.rand(F.mesh.nF,1)
+        F[self.Src0, 'b'] = b
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
         phi = np.random.rand(F.mesh.nC,2)
-        F[[self.Tx0,self.Tx1], 'phi'] = phi
-        self.assertTrue(np.all(F[[self.Tx0,self.Tx1], 'phi'] == phi))
+        F[[self.Src0,self.Src1], 'phi'] = phi
+        self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
 
         fdict = F[:,:]
         self.assertTrue(type(fdict) is dict)
         self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
 
         b = np.random.rand(F.mesh.nF, 2)
-        F[[self.Tx0, self.Tx1],'b'] = b
-        self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,))
-        self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,))
-        self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0]))
-        self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1]))
+        F[[self.Src0, self.Src1],'b'] = b
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1))
+        self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2)))
+        self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2)))
 
     def test_assertions(self):
-        freq = [self.Tx0, self.Tx1]
-        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTx)
+        freq = [self.Src0, self.Src1]
+        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
         def fun(): self.F[freq, 'b'] = bWrongSize
         self.assertRaises(ValueError, fun)
         def fun(): self.F[-999.]
@@ -120,69 +120,69 @@ class FieldsTest_Alias(unittest.TestCase):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        txLoc = np.r_[0,0,0.]
+        srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
-        Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
-        txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
-        survey = Survey.BaseSurvey(txList=txList)
+        Src3 = Survey.BaseSrc([rxList3],loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
         self.D = Survey.Data(survey)
         self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
-        self.Tx0 = Tx0
-        self.Tx1 = Tx1
+        self.Src0 = Src0
+        self.Src1 = Src1
         self.mesh = mesh
         self.XYZ = XYZ
 
     def test_contains(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         self.assertTrue('b' not in F)
         self.assertTrue('e' not in F)
-        e = np.random.rand(F.mesh.nE, nTx)
+        e = np.random.rand(F.mesh.nE, nSrc)
         F[:, 'e'] = e
         self.assertTrue('b' in F)
         self.assertTrue('e' in F)
 
     def test_simpleAlias(self):
         F = self.F
-        nTx = F.survey.nTx
-        e = np.random.rand(F.mesh.nE, nTx)
+        nSrc = F.survey.nSrc
+        e = np.random.rand(F.mesh.nE, nSrc)
         F[:, 'e'] = e
         self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e ))
 
         e = np.random.rand(F.mesh.nE,1)
-        F[self.Tx0, 'e'] = e
-        self.assertTrue(np.all(F[self.Tx0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e)))
+        F[self.Src0, 'e'] = e
+        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e))
 
         def f():
-            F[self.Tx0, 'b'] = F[self.Tx0, 'b']
+            F[self.Src0, 'b'] = F[self.Src0, 'b']
         self.assertRaises(KeyError, f) # can't set a alias attr.
 
     def test_aliasFunction(self):
         def alias(e, ind):
-            self.assertTrue(ind is self.Tx0)
+            self.assertTrue(ind is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1)
-        F[self.Tx0, 'e'] = e
-        F[self.Tx0, 'b']
+        F[self.Src0, 'e'] = e
+        F[self.Src0, 'b']
 
 
         def alias(e, ind):
             self.assertTrue(type(ind) is list)
-            self.assertTrue(ind[0] is self.Tx0)
-            self.assertTrue(ind[1] is self.Tx1)
+            self.assertTrue(ind[0] is self.Src0)
+            self.assertTrue(ind[1] is self.Src1)
             return self.F.mesh.edgeCurl * e
         F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,2)
-        F[[self.Tx0, self.Tx1], 'e'] = e
-        F[[self.Tx0, self.Tx1], 'b']
+        F[[self.Src0, self.Src1], 'e'] = e
+        F[[self.Src0, self.Src1], 'b']
 
 
 class FieldsTest_Time(unittest.TestCase):
@@ -191,34 +191,34 @@ class FieldsTest_Time(unittest.TestCase):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        txLoc = np.r_[0,0,0.]
+        srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
-        Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
-        txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
-        survey = Survey.BaseSurvey(txList=txList)
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
         survey.pair(prob)
         self.F = Problem.TimeFields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'})
-        self.Tx0 = Tx0
-        self.Tx1 = Tx1
+        self.Src0 = Src0
+        self.Src1 = Src1
         self.mesh = mesh
         self.XYZ = XYZ
 
     def test_contains(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         nT = F.survey.prob.nT + 1
         self.assertTrue('b' not in F)
         self.assertTrue('e' not in F)
         self.assertTrue('phi' not in F)
-        e = np.random.rand(F.mesh.nE, nTx, nT)
+        e = np.random.rand(F.mesh.nE, nSrc, nT)
         F[:, 'e', :] = e
         self.assertTrue('e' in F)
         self.assertTrue('b' not in F)
@@ -226,11 +226,11 @@ class FieldsTest_Time(unittest.TestCase):
 
     def test_SetGet(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         nT = F.survey.prob.nT + 1
-        e = np.random.rand(F.mesh.nE, nTx, nT)
+        e = np.random.rand(F.mesh.nE, nSrc, nT)
         F[:, 'e'] = e
-        b = np.random.rand(F.mesh.nF, nTx, nT)
+        b = np.random.rand(F.mesh.nF, nSrc, nT)
         F[:, 'b'] = b
 
         self.assertTrue(np.all(F[:, 'e'] == e))
@@ -244,39 +244,39 @@ class FieldsTest_Time(unittest.TestCase):
             self.assertTrue(np.all(F[:, 'b'] == b*0))
 
         b = np.random.rand(F.mesh.nF,1,nT)
-        F[self.Tx0, 'b'] = b
-        self.assertTrue(np.all(F[self.Tx0, 'b'] == b[:,0,:]))
+        F[self.Src0, 'b'] = b
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b[:,0,:]))
 
         b = np.random.rand(F.mesh.nF,1,nT)
-        F[self.Tx0, 'b', 0] = b[:,:,0]
-        self.assertTrue(np.all(F[self.Tx0, 'b', 0] == b[:,0,0]))
+        F[self.Src0, 'b', 0] = b[:,:,0]
+        self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2)))
 
         phi = np.random.rand(F.mesh.nC,2,nT)
-        F[[self.Tx0,self.Tx1], 'phi'] = phi
-        self.assertTrue(np.all(F[[self.Tx0,self.Tx1], 'phi'] == phi))
+        F[[self.Src0,self.Src1], 'phi'] = phi
+        self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
 
         fdict = F[:]
         self.assertTrue(type(fdict) is dict)
         self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
 
         b = np.random.rand(F.mesh.nF, 2, nT)
-        F[[self.Tx0, self.Tx1],'b'] = b
-        self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,nT))
-        self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,nT))
-        self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0,:]))
-        self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1,:]))
-        self.assertTrue(np.all(F[self.Tx0,'b',1] == b[:,0,1]))
-        self.assertTrue(np.all(F[self.Tx1,'b',1] == b[:,1,1]))
-        self.assertTrue(np.all(F[self.Tx0,'b',4] == b[:,0,4]))
-        self.assertTrue(np.all(F[self.Tx1,'b',4] == b[:,1,4]))
+        F[[self.Src0, self.Src1],'b'] = b
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,nT))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
+        self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
+        self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
+        self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2)))
+        self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2)))
 
 
         b = np.random.rand(F.mesh.nF, 2, nT)
-        F[[self.Tx0, self.Tx1],'b', 0] = b[:,:,0]
+        F[[self.Src0, self.Src1],'b', 0] = b[:,:,0]
 
     def test_assertions(self):
-        freq = [self.Tx0, self.Tx1]
-        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTx)
+        freq = [self.Src0, self.Src1]
+        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
         def fun(): self.F[freq, 'b'] = bWrongSize
         self.assertRaises(ValueError, fun)
         def fun(): self.F[-999.]
@@ -293,35 +293,35 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        txLoc = np.r_[0,0,0.]
+        srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Tx0 = Survey.BaseTx(txLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Tx1 = Survey.BaseTx(txLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Tx2 = Survey.BaseTx(txLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Tx3 = Survey.BaseTx(txLoc, 'VMD', [rxList3])
-        Tx4 = Survey.BaseTx(txLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
-        txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
-        survey = Survey.BaseSurvey(txList=txList)
+        Src3 = Survey.BaseSrc( [rxList3],loc=srcLoc)
+        Src4 = Survey.BaseSrc( [rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
         survey.pair(prob)
-        def alias(b, txInd, timeInd):
+        def alias(b, srcInd, timeInd):
             return self.F.mesh.edgeCurl.T * b + timeInd
         self.F = Problem.TimeFields(mesh, survey, knownFields={'b':'F'}, aliasFields={'e':['b','E',alias]})
-        self.Tx0 = Tx0
-        self.Tx1 = Tx1
+        self.Src0 = Src0
+        self.Src1 = Src1
         self.mesh = mesh
         self.XYZ = XYZ
 
     def test_contains(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         nT = F.survey.prob.nT + 1
         self.assertTrue('b' not in F)
         self.assertTrue('e' not in F)
-        b = np.random.rand(F.mesh.nF, nTx, nT)
+        b = np.random.rand(F.mesh.nF, nSrc, nT)
         F[:, 'b', :] = b
         self.assertTrue('e' in F)
         self.assertTrue('b' in F)
@@ -329,9 +329,9 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
 
     def test_simpleAlias(self):
         F = self.F
-        nTx = F.survey.nTx
+        nSrc = F.survey.nSrc
         nT = F.survey.prob.nT + 1
-        b = np.random.rand(F.mesh.nF, nTx, nT)
+        b = np.random.rand(F.mesh.nF, nSrc, nT)
         F[:, 'b', :] = b
         self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:,:,0] ))
 
@@ -341,57 +341,57 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
             e[i] = e[i][:,:,np.newaxis]
         e = np.concatenate(e, axis=2)
         self.assertTrue(np.all(F[:, 'e', :] == e ))
-        self.assertTrue(np.all(F[self.Tx0, 'e', :] == e[:,0,:] ))
-        self.assertTrue(np.all(F[self.Tx1, 'e', :] == e[:,1,:] ))
+        self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
+        self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
         for t in range(nT):
-            self.assertTrue(np.all(F[self.Tx1, 'e', t] == e[:,1,t] ))
+            self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) ))
 
         b = np.random.rand(F.mesh.nF,nT)
-        F[self.Tx0, 'b',:] = b
+        F[self.Src0, 'b',:] = b
         Cb = F.mesh.edgeCurl.T * b
         for i in range(Cb.shape[1]):
             Cb[:,i] += i
-        self.assertTrue(np.all(F[self.Tx0, 'e',:] == Cb))
+        self.assertTrue(np.all(F[self.Src0, 'e',:] == Cb))
 
         def f():
-            F[self.Tx0, 'e'] = F[self.Tx0, 'e']
+            F[self.Src0, 'e'] = F[self.Src0, 'e']
         self.assertRaises(KeyError, f) # can't set a alias attr.
 
     def test_aliasFunction(self):
         nT = self.F.survey.prob.nT + 1
         count = [0]
-        def alias(e, txInd, timeInd):
+        def alias(e, srcInd, timeInd):
             count[0] += 1
-            self.assertTrue(txInd is self.Tx0)
+            self.assertTrue(srcInd is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1,nT)
-        F[self.Tx0, 'e', :] = e
-        F[self.Tx0, 'b', :]
+        F[self.Src0, 'e', :] = e
+        F[self.Src0, 'b', :]
         self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
         e = np.random.rand(F.mesh.nE,1,1)
-        F[self.Tx0, 'e', 1] = e
+        F[self.Src0, 'e', 1] = e
         count[0] = 0
-        F[self.Tx0, 'b', 1]
+        F[self.Src0, 'b', 1]
         self.assertTrue(count[0] == 1) # ensure that this is called only once.
 
 
-        def alias(e, txInd, timeInd):
+        def alias(e, srcInd, timeInd):
             count[0] += 1
-            self.assertTrue(type(txInd) is list)
-            self.assertTrue(txInd[0] is self.Tx0)
-            self.assertTrue(txInd[1] is self.Tx1)
+            self.assertTrue(type(srcInd) is list)
+            self.assertTrue(srcInd[0] is self.Src0)
+            self.assertTrue(srcInd[1] is self.Src1)
             return self.F.mesh.edgeCurl * e
         F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,2, nT)
-        F[[self.Tx0, self.Tx1], 'e', :] = e
+        F[[self.Src0, self.Src1], 'e', :] = e
         count[0] = 0
-        F[[self.Tx0, self.Tx1], 'b', :]
+        F[[self.Src0, self.Src1], 'b', :]
         self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
         e = np.random.rand(F.mesh.nE,2, 1)
-        F[[self.Tx0, self.Tx1], 'e', 1] = e
+        F[[self.Src0, self.Src1], 'e', 1] = e
         count[0] = 0
-        F[[self.Tx0, self.Tx1], 'b', 1]
+        F[[self.Src0, self.Src1], 'b', 1]
         self.assertTrue(count[0] == 1) # ensure that this is called only once.
 
 

SimPEG/Tests/test_TreeMesh.py

@@ -483,22 +483,22 @@ class SimpleOctreeOperatorTests(unittest.TestCase):
         self.oM2 = TreeMesh([h1,h2])
 
     def test_faceDiv(self):
-        self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0)
-        self.assertTrue((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum() == 0)
+        self.assertAlmostEqual((self.tM.faceDiv - self.oM.faceDiv).toarray().sum(), 0)
+        self.assertAlmostEqual((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum(), 0)
 
     def test_nodalGrad(self):
-        self.assertTrue((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum() == 0)
-        self.assertTrue((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum() == 0)
+        self.assertAlmostEqual((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum(), 0)
+        self.assertAlmostEqual((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum(), 0)
 
     def test_edgeCurl(self):
-        self.assertTrue((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum() == 0)
-        # self.assertTrue((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum() == 0)
+        self.assertAlmostEqual((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum(), 0)
+        # self.assertAlmostEqual((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum(), 0)
 
     def test_InnerProducts(self):
-        self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() < TOL)
-        self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() < TOL)
-        self.assertTrue((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum() < TOL)
-        self.assertTrue((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum() < TOL)
+        self.assertAlmostEqual((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum(), 0)
+        self.assertAlmostEqual((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum(), 0)
+        self.assertAlmostEqual((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum(), 0)
+        self.assertAlmostEqual((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum(), 0)
 
 
 if __name__ == '__main__':

SimPEG/Tests/test_boundaryPoisson.py

@@ -20,9 +20,9 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
         j_fun = lambda x: -np.pi*np.sin(np.pi*x)
         q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        j_anal = j_fun(self.M.gridFx)
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
 
         #TODO: Check where our boundary conditions are CCx or Nx
         # vec = self.M.vectorNx
@@ -38,32 +38,32 @@ class Test1D_InhomogeneousDirichlet(OrderTest):
         V = Utils.sdiag(self.M.vol)
         G = -Pin.T*Pin*self.M.faceDiv.T * V
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*phi_bc)
+        j = McI*(G*xc_ana + P*phi_bc)
         q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
 
         # Rearrange if we know q to solve for x
         A = V*D*Pin.T*Pin*McI*G
-        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
         # A = D*McI*G
-        # rhs = q_anal - D*McI*P*phi_bc
+        # rhs = q_ana - D*McI*P*phi_bc
 
 
         if self.myTest == 'j':
-            err = np.linalg.norm((j-j_anal), np.inf)
+            err = np.linalg.norm((j-j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-V*q_anal), np.inf)
+            err = np.linalg.norm((q-V*q_ana), np.inf)
         elif self.myTest == 'xc':
             #TODO: fix the null space
             solver = SolverCG(A, maxiter=1000)
             xc = solver * (rhs)
             print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
         elif self.myTest == 'xcJ':
             #TODO: fix the null space
             xc = Solver(A) * (rhs)
             print np.linalg.norm(Utils.mkvc(A*xc) - rhs)
             j = McI*(G*xc + P*phi_bc)
-            err = np.linalg.norm((j-j_anal), np.inf)
+            err = np.linalg.norm((j-j_ana), np.inf)
 
         return err
 
@@ -102,11 +102,11 @@ class Test2D_InhomogeneousDirichlet(OrderTest):
         j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
         q_fun = lambda x: -2*(np.pi**2)*phi(x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        jX_anal = j_funX(self.M.gridFx)
-        jY_anal = j_funY(self.M.gridFy)
-        j_anal = np.r_[jX_anal,jY_anal]
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
 
         #TODO: Check where our boundary conditions are CCx or Nx
         # fxm,fxp,fym,fyp = self.M.faceBoundaryInd
@@ -126,26 +126,26 @@ class Test2D_InhomogeneousDirichlet(OrderTest):
         McI = Utils.sdInv(self.M.getFaceInnerProduct())
         G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*bc)
+        j = McI*(G*xc_ana + P*bc)
         q = D*j
 
         # self.M.plotImage(j, 'FxFy', showIt=True)
 
         # Rearrange if we know q to solve for x
         A = D*McI*G
-        rhs = q_anal - D*McI*P*bc
+        rhs = q_ana - D*McI*P*bc
 
         if self.myTest == 'j':
-            err = np.linalg.norm((j-j_anal), np.inf)
+            err = np.linalg.norm((j-j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-q_anal), np.inf)
+            err = np.linalg.norm((q-q_ana), np.inf)
         elif self.myTest == 'xc':
             xc = Solver(A) * (rhs)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
         elif self.myTest == 'xcJ':
             xc = Solver(A) * (rhs)
             j = McI*(G*xc + P*bc)
-            err = np.linalg.norm((j-j_anal), np.inf)
+            err = np.linalg.norm((j-j_ana), np.inf)
 
         return err
 
@@ -182,9 +182,9 @@ class Test1D_InhomogeneousNeumann(OrderTest):
         j_fun = lambda x: np.pi*np.cos(np.pi*x)
         q_fun = lambda x: -(np.pi**2)*np.sin(np.pi*x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        j_anal = j_fun(self.M.gridFx)
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
 
         #TODO: Check where our boundary conditions are CCx or Nx
         vecN = self.M.vectorNx
@@ -200,24 +200,24 @@ class Test1D_InhomogeneousNeumann(OrderTest):
         V = Utils.sdiag(self.M.vol)
         G = -Pin.T*Pin*self.M.faceDiv.T * V
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*phi_bc)
+        j = McI*(G*xc_ana + P*phi_bc)
         q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
 
         # Rearrange if we know q to solve for x
         A = V*D*Pin.T*Pin*McI*G
-        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
         # A = D*McI*G
-        # rhs = q_anal - D*McI*P*phi_bc
+        # rhs = q_ana - D*McI*P*phi_bc
 
 
         if self.myTest == 'j':
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-V*q_anal), np.inf)
+            err = np.linalg.norm((q-V*q_ana), np.inf)
         elif self.myTest == 'xc':
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -225,7 +225,7 @@ class Test1D_InhomogeneousNeumann(OrderTest):
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
             j = McI*(G*xc + P*phi_bc)
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -261,11 +261,11 @@ class Test2D_InhomogeneousNeumann(OrderTest):
         j_funY = lambda x: np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
         q_fun = lambda x: -2*(np.pi**2)*phi(x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        jX_anal = j_funX(self.M.gridFx)
-        jY_anal = j_funY(self.M.gridFy)
-        j_anal = np.r_[jX_anal,jY_anal]
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
 
         #TODO: Check where our boundary conditions are CCx or Nx
 
@@ -290,21 +290,21 @@ class Test2D_InhomogeneousNeumann(OrderTest):
         V = Utils.sdiag(self.M.vol)
         G = -Pin.T*Pin*self.M.faceDiv.T * V
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*phi_bc)
+        j = McI*(G*xc_ana + P*phi_bc)
         q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
 
         # Rearrange if we know q to solve for x
         A = V*D*Pin.T*Pin*McI*G
-        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
 
         if self.myTest == 'j':
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-V*q_anal), np.inf)
+            err = np.linalg.norm((q-V*q_ana), np.inf)
         elif self.myTest == 'xc':
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -312,7 +312,7 @@ class Test2D_InhomogeneousNeumann(OrderTest):
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
             j = McI*(G*xc + P*phi_bc)
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -346,9 +346,9 @@ class Test1D_InhomogeneousMixed(OrderTest):
         j_fun = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x)
         q_fun = lambda x: -0.25*(np.pi**2)*np.cos(0.5*np.pi*x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        j_anal = j_fun(self.M.gridFx)
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
 
         #TODO: Check where our boundary conditions are CCx or Nx
         vecN = self.M.vectorNx
@@ -364,24 +364,24 @@ class Test1D_InhomogeneousMixed(OrderTest):
         V = Utils.sdiag(self.M.vol)
         G = -Pin.T*Pin*self.M.faceDiv.T * V
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*phi_bc)
+        j = McI*(G*xc_ana + P*phi_bc)
         q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
 
         # Rearrange if we know q to solve for x
         A = V*D*Pin.T*Pin*McI*G
-        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
         # A = D*McI*G
-        # rhs = q_anal - D*McI*P*phi_bc
+        # rhs = q_ana - D*McI*P*phi_bc
 
 
         if self.myTest == 'j':
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-V*q_anal), np.inf)
+            err = np.linalg.norm((q-V*q_ana), np.inf)
         elif self.myTest == 'xc':
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -389,7 +389,7 @@ class Test1D_InhomogeneousMixed(OrderTest):
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
             j = McI*(G*xc + P*phi_bc)
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -425,11 +425,11 @@ class Test2D_InhomogeneousMixed(OrderTest):
         j_funY = lambda x: -0.5*np.pi*np.cos(0.5*np.pi*x[:,0])*np.sin(0.5*np.pi*x[:,1])
         q_fun = lambda x: -2*((0.5*np.pi)**2)*phi(x)
 
-        xc_anal = phi(self.M.gridCC)
-        q_anal = q_fun(self.M.gridCC)
-        jX_anal = j_funX(self.M.gridFx)
-        jY_anal = j_funY(self.M.gridFy)
-        j_anal = np.r_[jX_anal,jY_anal]
+        xc_ana = phi(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
 
         #TODO: Check where our boundary conditions are CCx or Nx
 
@@ -454,21 +454,21 @@ class Test2D_InhomogeneousMixed(OrderTest):
         V = Utils.sdiag(self.M.vol)
         G = -Pin.T*Pin*self.M.faceDiv.T * V
         D = self.M.faceDiv
-        j = McI*(G*xc_anal + P*phi_bc)
+        j = McI*(G*xc_ana + P*phi_bc)
         q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
 
         # Rearrange if we know q to solve for x
         A = V*D*Pin.T*Pin*McI*G
-        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
 
         if self.myTest == 'j':
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
         elif self.myTest == 'q':
-            err = np.linalg.norm((q-V*q_anal), np.inf)
+            err = np.linalg.norm((q-V*q_ana), np.inf)
         elif self.myTest == 'xc':
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
-            err = np.linalg.norm((xc-xc_anal), np.inf)
+            err = np.linalg.norm((xc-xc_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
@@ -476,7 +476,7 @@ class Test2D_InhomogeneousMixed(OrderTest):
             #TODO: fix the null space
             xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
             j = McI*(G*xc + P*phi_bc)
-            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
             if info > 0:
                 print 'Solve does not work well'
                 print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)

SimPEG/Tests/test_cylMesh.py

@@ -133,6 +133,80 @@ class TestCyl2DMesh(unittest.TestCase):
     def test_lightOperators(self):
         self.assertTrue(self.mesh.nodalGrad is None)
 
+    def test_getInterpMatCartMesh_Cells(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'CC')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
+    def test_getInterpMatCartMesh_Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf = Mc.getInterpolationMatCartMesh(Mr, 'F')
+        mf = np.ones(Mc.nF)
+
+        frect = Pf * mf
+
+        fxcc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        fycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        fzcc = Mr.r(frect, 'F', 'Fz')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(fxcc[indX]) - 1) < TOL
+        assert np.abs(float(fxcc[indY]) - 0) < TOL
+        assert np.abs(float(fycc[indX]) - 0) < TOL
+        assert np.abs(float(fycc[indY]) - 1) < TOL
+        assert np.abs((fzcc - 1).sum()) < TOL
+
+        mag = (fxcc**2 + fycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
+        me = np.ones(Mc.nE)
+
+        erect = Pe * me
+
+        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
+        ezcc = Mr.r(erect, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
 MESHTYPES = ['uniformCylMesh']
@@ -160,9 +234,9 @@ class TestFaceDiv2D(OrderTest):
         F = self.M.projectFaceVector(Fc)
 
         divF = self.M.faceDiv.dot(F)
-        divF_anal = call3(sol, self.M.gridCC)
+        divF_ana = call3(sol, self.M.gridCC)
 
-        err = np.linalg.norm((divF-divF_anal), np.inf)
+        err = np.linalg.norm((divF-divF_ana), np.inf)
         return err
 
     def test_order(self):
@@ -198,9 +272,9 @@ class TestEdgeCurl2D(OrderTest):
 
         Fc = cylF2(self.M, solR, solZ)
         Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
-        curlE_anal = self.M.projectFaceVector(Fc)
+        curlE_ana = self.M.projectFaceVector(Fc)
 
-        err = np.linalg.norm((curlE-curlE_anal), np.inf)
+        err = np.linalg.norm((curlE-curlE_ana), np.inf)
         return err
 
     def test_order(self):

SimPEG/Tests/test_innerProduct.py

@@ -7,7 +7,7 @@ from SimPEG import Utils
 class TestInnerProducts(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
-    meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+    meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
     meshDimension = 3
     meshSizes = [16, 32]
 
@@ -154,7 +154,7 @@ class TestInnerProducts(OrderTest):
 class TestInnerProducts2D(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
-    meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+    meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
     meshDimension = 2
     meshSizes = [4, 8, 16, 32, 64, 128]
 

SimPEG/Tests/test_innerProductDerivs.py

@@ -7,9 +7,9 @@ from TestUtils import checkDerivative
 class TestInnerProductsDerivs(unittest.TestCase):
 
     def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
-        if meshType == 'LRM':
+        if meshType == 'Curv':
             hRect = Utils.exampleLrmGrid(h,'rotate')
-            mesh = Mesh.LogicallyRectMesh(hRect)
+            mesh = Mesh.CurvilinearMesh(hRect)
         elif meshType == 'Tree':
             mesh = Mesh.TreeMesh(h)
         elif meshType == 'Tensor':
@@ -24,9 +24,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
-        if meshType == 'LRM':
+        if meshType == 'Curv':
             hRect = Utils.exampleLrmGrid(h,'rotate')
-            mesh = Mesh.LogicallyRectMesh(hRect)
+            mesh = Mesh.CurvilinearMesh(hRect)
         elif meshType == 'Tree':
             mesh = Mesh.TreeMesh(h)
         elif meshType == 'Tensor':
@@ -137,65 +137,65 @@ class TestInnerProductsDerivs(unittest.TestCase):
 
 
 
-    def test_FaceIP_2D_float_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
-    def test_FaceIP_3D_float_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM'))
-    def test_FaceIP_2D_isotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM'))
-    def test_FaceIP_3D_isotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM'))
-    def test_FaceIP_2D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM'))
-    def test_FaceIP_3D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM'))
-    def test_FaceIP_2D_tensor_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM'))
-    def test_FaceIP_3D_tensor_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM'))
+    def test_FaceIP_2D_float_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],0, False, 'Curv'))
+    def test_FaceIP_3D_float_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Curv'))
+    def test_FaceIP_2D_isotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],1, False, 'Curv'))
+    def test_FaceIP_3D_isotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Curv'))
+    def test_FaceIP_2D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],2, False, 'Curv'))
+    def test_FaceIP_3D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Curv'))
+    def test_FaceIP_2D_tensor_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],3, False, 'Curv'))
+    def test_FaceIP_3D_tensor_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Curv'))
 
-    def test_FaceIP_2D_float_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM'))
-    def test_FaceIP_3D_float_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM'))
-    def test_FaceIP_2D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM'))
-    def test_FaceIP_3D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM'))
-    def test_FaceIP_2D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM'))
-    def test_FaceIP_3D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM'))
+    def test_FaceIP_2D_float_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'Curv'))
+    def test_FaceIP_3D_float_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Curv'))
+    def test_FaceIP_2D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'Curv'))
+    def test_FaceIP_3D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Curv'))
+    def test_FaceIP_2D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'Curv'))
+    def test_FaceIP_3D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Curv'))
 
-    def test_EdgeIP_2D_float_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM'))
-    def test_EdgeIP_3D_float_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM'))
-    def test_EdgeIP_2D_isotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM'))
-    def test_EdgeIP_3D_isotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM'))
-    def test_EdgeIP_2D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM'))
-    def test_EdgeIP_3D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM'))
-    def test_EdgeIP_2D_tensor_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM'))
-    def test_EdgeIP_3D_tensor_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM'))
+    def test_EdgeIP_2D_float_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, False, 'Curv'))
+    def test_EdgeIP_3D_float_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Curv'))
+    def test_EdgeIP_2D_isotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, False, 'Curv'))
+    def test_EdgeIP_3D_isotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Curv'))
+    def test_EdgeIP_2D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, False, 'Curv'))
+    def test_EdgeIP_3D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Curv'))
+    def test_EdgeIP_2D_tensor_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],3, False, 'Curv'))
+    def test_EdgeIP_3D_tensor_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Curv'))
 
-    def test_EdgeIP_2D_float_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM'))
-    def test_EdgeIP_3D_float_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM'))
-    def test_EdgeIP_2D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM'))
-    def test_EdgeIP_3D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM'))
-    def test_EdgeIP_2D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM'))
-    def test_EdgeIP_3D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM'))
+    def test_EdgeIP_2D_float_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, True, 'Curv'))
+    def test_EdgeIP_3D_float_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Curv'))
+    def test_EdgeIP_2D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, True, 'Curv'))
+    def test_EdgeIP_3D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Curv'))
+    def test_EdgeIP_2D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, True, 'Curv'))
+    def test_EdgeIP_3D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Curv'))
 
 
 

SimPEG/Tests/test_interpolation.py

@@ -34,7 +34,7 @@ class TestInterpolation1D(OrderTest):
     def getError(self):
         funX = lambda x: np.cos(2*np.pi*x)
 
-        anal = call1(funX, self.LOCS)
+        ana = call1(funX, self.LOCS)
 
         if 'CC' == self.type:
             grid = call1(funX, self.M.gridCC)
@@ -43,7 +43,7 @@ class TestInterpolation1D(OrderTest):
 
         comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
 
-        err = np.linalg.norm((comp - anal), 2)
+        err = np.linalg.norm((comp - ana), 2)
         return err
 
     def test_orderCC(self):
@@ -82,11 +82,11 @@ class TestInterpolation2d(OrderTest):
         funY = lambda x, y: np.cos(2*np.pi*x)
 
         if 'x' in self.type:
-            anal = call2(funX, self.LOCS)
+            ana = call2(funX, self.LOCS)
         elif 'y' in self.type:
-            anal = call2(funY, self.LOCS)
+            ana = call2(funY, self.LOCS)
         else:
-            anal = call2(funX, self.LOCS)
+            ana = call2(funX, self.LOCS)
 
         if 'F' in self.type:
             Fc = cartF2(self.M, funX, funY)
@@ -101,7 +101,7 @@ class TestInterpolation2d(OrderTest):
 
         comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
 
-        err = np.linalg.norm((comp - anal), np.inf)
+        err = np.linalg.norm((comp - ana), np.inf)
         return err
 
     def test_orderCC(self):
@@ -165,13 +165,13 @@ class TestInterpolation2dCyl(OrderTest):
         funY = lambda x, y: np.cos(2*np.pi*x)
 
         if 'x' in self.type:
-            anal = call2(funX, self.LOCS)
+            ana = call2(funX, self.LOCS)
         elif 'y' in self.type:
-            anal = call2(funY, self.LOCS)
+            ana = call2(funY, self.LOCS)
         elif 'z' in self.type:
-            anal = call2(funY, self.LOCS)
+            ana = call2(funY, self.LOCS)
         else:
-            anal = call2(funX, self.LOCS)
+            ana = call2(funX, self.LOCS)
 
         if 'Fx' == self.type:
             Fc = cartF2Cyl(self.M, funX, funY)
@@ -192,7 +192,7 @@ class TestInterpolation2dCyl(OrderTest):
 
         comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
 
-        err = np.linalg.norm((comp - anal), np.inf)
+        err = np.linalg.norm((comp - ana), np.inf)
         return err
 
     def test_orderCC(self):
@@ -234,13 +234,13 @@ class TestInterpolation3D(OrderTest):
         funZ = lambda x, y, z: np.cos(2*np.pi*x)
 
         if 'x' in self.type:
-            anal = call3(funX, self.LOCS)
+            ana = call3(funX, self.LOCS)
         elif 'y' in self.type:
-            anal = call3(funY, self.LOCS)
+            ana = call3(funY, self.LOCS)
         elif 'z' in self.type:
-            anal = call3(funZ, self.LOCS)
+            ana = call3(funZ, self.LOCS)
         else:
-            anal = call3(funX, self.LOCS)
+            ana = call3(funX, self.LOCS)
 
         if 'F' in self.type:
             Fc = cartF3(self.M, funX, funY, funZ)
@@ -255,7 +255,7 @@ class TestInterpolation3D(OrderTest):
 
         comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
 
-        err = np.linalg.norm((comp - anal), np.inf)
+        err = np.linalg.norm((comp - ana), np.inf)
         return err
 
     def test_orderCC(self):

SimPEG/Tests/test_maps.py

@@ -6,6 +6,9 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
+MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
+
 class MapTests(unittest.TestCase):
 
     def setUp(self):
@@ -13,17 +16,31 @@ class MapTests(unittest.TestCase):
         a = np.array([1, 1, 1])
         b = np.array([1, 2])
         self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
+        self.mesh3 = Mesh.TensorMesh([a, b, [3,4]], x0=np.array([3, 5, 2]))
         self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
 
-    def test_transforms(self):
-        for M in dir(Maps):
-            try:
-                maps = getattr(Maps, M)(self.mesh2)
-                assert isinstance(maps, Maps.IdentityMap)
-            except Exception, e:
-                continue
+    def test_transforms2D(self):
+        for M in MAPS_TO_TEST_2D:
+            maps = getattr(Maps, M)(self.mesh2)
             self.assertTrue(maps.test())
 
+    def test_transforms3D(self):
+        for M in MAPS_TO_TEST_3D:
+            maps = getattr(Maps, M)(self.mesh3)
+            self.assertTrue(maps.test())
+
+
+    def test_transforms_logMap(self):
+        # Note that log maps can be kinda finicky, so we are being explicit about the random seed.
+        v2 = np.r_[ 0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519, 0.79754415]
+        dv2 = np.r_[ 0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762, 0.44252488]
+        v3 = np.r_[ 0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109, 0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864, 0.84130763, 0.22123854]
+        dv3 = np.r_[ 0.96827838, 0.26072111, 0.45090749, 0.10573893, 0.65276365, 0.15646586, 0.51679682, 0.23071984, 0.95106218, 0.14201845, 0.25093564, 0.3732866 ]
+        maps = Maps.LogMap(self.mesh2)
+        self.assertTrue(maps.test(v2, dx=dv2))
+        maps = Maps.LogMap(self.mesh3)
+        self.assertTrue(maps.test(v3, dx=dv3))
+
     def test_Mesh2MeshMap(self):
         maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
         self.assertTrue(maps.test())
@@ -90,5 +107,34 @@ class MapTests(unittest.TestCase):
         self.assertRaises(ValueError, lambda: actMap * vertMap * expMap )
 
 
+    def test_map2Dto3D_x(self):
+        M2 = Mesh.TensorMesh([2,4])
+        M3 = Mesh.TensorMesh([3,2,4])
+        m = np.random.rand(M2.nC)
+        m2to3 = Maps.Map2Dto3D(M3, normal='X')
+        m = np.arange(m2to3.nP)
+        self.assertTrue(m2to3.test())
+        self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[0,:,:] ) == m))
+
+
+    def test_map2Dto3D_y(self):
+        M2 = Mesh.TensorMesh([3,4])
+        M3 = Mesh.TensorMesh([3,2,4])
+        m = np.random.rand(M2.nC)
+        m2to3 = Maps.Map2Dto3D(M3, normal='Y')
+        m = np.arange(m2to3.nP)
+        self.assertTrue(m2to3.test())
+        self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,0,:] ) == m))
+
+    def test_map2Dto3D_z(self):
+        M2 = Mesh.TensorMesh([3,2])
+        M3 = Mesh.TensorMesh([3,2,4])
+        m = np.random.rand(M2.nC)
+        m2to3 = Maps.Map2Dto3D(M3, normal='Z')
+        m = np.arange(m2to3.nP)
+        self.assertTrue(m2to3.test())
+        self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,:,0] ) == m))
+
+
 if __name__ == '__main__':
     unittest.main()

SimPEG/Tests/test_operators.py

@@ -4,7 +4,7 @@ from TestUtils import OrderTest
 import matplotlib.pyplot as plt
 
 #TODO: 'randomTensorMesh'
-MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
 call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
 cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
@@ -35,20 +35,41 @@ class TestCurl(OrderTest):
         E = self.M.projectEdgeVector(Ec)
 
         Fc = cartF3(self.M, solX, solY, solZ)
-        curlE_anal = self.M.projectFaceVector(Fc)
+        curlE_ana = self.M.projectFaceVector(Fc)
 
         curlE = self.M.edgeCurl.dot(E)
-        if self._meshType == 'rotateLRM':
+        if self._meshType == 'rotateCurv':
             # Really it is the integration we should be caring about:
             # So, let us look at the l2 norm.
-            err = np.linalg.norm(self.M.area*(curlE - curlE_anal), 2)
+            err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2)
         else:
-            err = np.linalg.norm((curlE - curlE_anal), np.inf)
+            err = np.linalg.norm((curlE - curlE_ana), np.inf)
         return err
 
     def test_order(self):
         self.orderTest()
 
+class TestCurl2D(OrderTest):
+    name = "Cell Grad 2D - Dirichlet"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 2
+    meshSizes = [8, 16, 32, 64]
+
+    def getError(self):
+        #Test function
+        ex = lambda x, y: np.cos(y)
+        ey = lambda x, y: np.cos(x)
+        sol = lambda x, y: -np.sin(x)+np.sin(y)
+
+        sol_curl2d = call2(sol, self.M.gridCC)
+        Ec = cartE2(self.M, ex, ey)
+        sol_ana = self.M.edgeCurl*self.M.projectFaceVector(Ec)
+        err = np.linalg.norm((sol_curl2d-sol_ana), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
 
 class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
     name = "Cell Grad 1D - Dirichlet"
@@ -65,13 +86,13 @@ class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
 
         xc = sol(self.M.gridCC)
 
-        gradX_anal = fx(self.M.gridFx)
+        gradX_ana = fx(self.M.gridFx)
 
         bc = np.array([1,1])
         self.M.setCellGradBC('dirichlet')
         gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc
 
-        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+        err = np.linalg.norm((gradX-gradX_ana), np.inf)
 
         return err
 
@@ -93,12 +114,12 @@ class TestCellGrad2D_Dirichlet(OrderTest):
         xc = call2(sol, self.M.gridCC)
 
         Fc = cartF2(self.M, fx, fy)
-        gradX_anal = self.M.projectFaceVector(Fc)
+        gradX_ana = self.M.projectFaceVector(Fc)
 
         self.M.setCellGradBC('dirichlet')
         gradX = self.M.cellGrad.dot(xc)
 
-        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+        err = np.linalg.norm((gradX-gradX_ana), np.inf)
 
         return err
 
@@ -122,12 +143,12 @@ class TestCellGrad3D_Dirichlet(OrderTest):
         xc = call3(sol, self.M.gridCC)
 
         Fc = cartF3(self.M, fx, fy, fz)
-        gradX_anal = self.M.projectFaceVector(Fc)
+        gradX_ana = self.M.projectFaceVector(Fc)
 
         self.M.setCellGradBC('dirichlet')
         gradX = self.M.cellGrad.dot(xc)
 
-        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+        err = np.linalg.norm((gradX-gradX_ana), np.inf)
 
         return err
 
@@ -149,12 +170,12 @@ class TestCellGrad2D_Neumann(OrderTest):
         xc = call2(sol, self.M.gridCC)
 
         Fc = cartF2(self.M, fx, fy)
-        gradX_anal = self.M.projectFaceVector(Fc)
+        gradX_ana = self.M.projectFaceVector(Fc)
 
         self.M.setCellGradBC('neumann')
         gradX = self.M.cellGrad.dot(xc)
 
-        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+        err = np.linalg.norm((gradX-gradX_ana), np.inf)
 
         return err
 
@@ -178,12 +199,12 @@ class TestCellGrad3D_Neumann(OrderTest):
         xc = call3(sol, self.M.gridCC)
 
         Fc = cartF3(self.M, fx, fy, fz)
-        gradX_anal = self.M.projectFaceVector(Fc)
+        gradX_ana = self.M.projectFaceVector(Fc)
 
         self.M.setCellGradBC('neumann')
         gradX = self.M.cellGrad.dot(xc)
 
-        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+        err = np.linalg.norm((gradX-gradX_ana), np.inf)
 
         return err
 
@@ -206,14 +227,14 @@ class TestFaceDiv3D(OrderTest):
         F = self.M.projectFaceVector(Fc)
 
         divF = self.M.faceDiv.dot(F)
-        divF_anal = call3(sol, self.M.gridCC)
+        divF_ana = call3(sol, self.M.gridCC)
 
-        if self._meshType == 'rotateLRM':
+        if self._meshType == 'rotateCurv':
             # Really it is the integration we should be caring about:
             # So, let us look at the l2 norm.
-            err = np.linalg.norm(self.M.vol*(divF-divF_anal), 2)
+            err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2)
         else:
-            err = np.linalg.norm((divF-divF_anal), np.inf)
+            err = np.linalg.norm((divF-divF_ana), np.inf)
         return err
 
     def test_order(self):
@@ -236,9 +257,9 @@ class TestFaceDiv2D(OrderTest):
         F = self.M.projectFaceVector(Fc)
 
         divF = self.M.faceDiv.dot(F)
-        divF_anal = call2(sol, self.M.gridCC)
+        divF_ana = call2(sol, self.M.gridCC)
 
-        err = np.linalg.norm((divF-divF_anal), np.inf)
+        err = np.linalg.norm((divF-divF_ana), np.inf)
 
         return err
 
@@ -262,9 +283,9 @@ class TestNodalGrad(OrderTest):
         gradE = self.M.nodalGrad.dot(phi)
 
         Ec = cartE3(self.M, solX, solY, solZ)
-        gradE_anal = self.M.projectEdgeVector(Ec)
+        gradE_ana = self.M.projectEdgeVector(Ec)
 
-        err = np.linalg.norm((gradE-gradE_anal), np.inf)
+        err = np.linalg.norm((gradE-gradE_ana), np.inf)
 
         return err
 
@@ -288,9 +309,9 @@ class TestNodalGrad2D(OrderTest):
         gradE = self.M.nodalGrad.dot(phi)
 
         Ec = cartE2(self.M, solX, solY)
-        gradE_anal = self.M.projectEdgeVector(Ec)
+        gradE_ana = self.M.projectEdgeVector(Ec)
 
-        err = np.linalg.norm((gradE-gradE_anal), np.inf)
+        err = np.linalg.norm((gradE-gradE_ana), np.inf)
 
         return err
 

SimPEG/Tests/test_utils.py

@@ -246,5 +246,18 @@ class TestSequenceFunctions(unittest.TestCase):
         self.assertTrue(np.all(true == listArray))
         self.assertTrue(true.shape == listArray.shape)
 
+class TestDiagEst(unittest.TestCase):
+
+    def setUp(self):
+        self.n = 10
+        self.A = np.random.rand(self.n,self.n)
+        self.Adiag = np.diagonal(self.A)
+
+    def testOnes(self):
+        Adiagtest = diagEst(self.A,self.n,self.n)
+        r = np.abs(Adiagtest-self.Adiag)
+        self.assertTrue(r.dot(r) < TOL)
+
+
 if __name__ == '__main__':
     unittest.main()

SimPEG/Utils/__init__.py

@@ -1,7 +1,7 @@
 from matutils import *
 from codeutils import *
-from meshutils import exampleLrmGrid, meshTensor, closestPoints, writeUBCTensorMesh, writeUBCTensorModel
-from lrmutils import volTetra, faceInfo, indexCube
+from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
+from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
 from CounterUtils import *

SimPEG/Utils/matutils.py

@@ -341,6 +341,61 @@ def invPropertyTensor(M, tensor, returnMatrix=False):
     return T
 
 
+def diagEst(matFun, n, k=None, approach='Probing'):
+    """ 
+    Estimate the diagonal of a matrix, A. Note that the matrix may be a function which returns A times a vector.
+
+    Three different approaches have been implemented, 
+        1. Probing : uses cyclic permutations of vectors with ones and zeros (default)
+        2. Ones : random +/- 1 entries
+        3. Random : random vectors
+
+    :param lambda (numpy.array) matFun: matrix to estimate the diagonal of
+    :param int64 n: size of the vector that should be used to compute matFun(v)
+    :param int64 k: number of vectors to be used to estimate the diagonal
+    :param str approach: approach to be used for getting vectors
+    :rtype: numpy.array
+    :return: est_diag(A)
+
+    Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf
+    """
+
+    if type(matFun).__name__=='ndarray':
+        A = matFun
+        matFun = lambda v: A.dot(v) 
+
+    if k is None:
+        k = np.floor(n/10.)
+
+    if approach =='Ones':
+        def getv(n,i=None):
+            v = np.random.randn(n)
+            v[v<0]  = -1.
+            v[v>=0] = 1.
+            return v
+
+    elif approach == 'Random':
+        def getv(n,i=None):
+            return np.random.randn(n)
+
+    else: #if approach == 'Probing':
+        def getv(n,i):
+            v = np.zeros(n)
+            v[i:n:k] = 1.
+            return v
+
+    Mv = np.zeros(n)
+    vv = np.zeros(n)
+
+    for i in range(0,k):
+        vk = getv(n,i)
+        Mv += matFun(vk)*vk
+        vv += vk*vk
+
+    d   = Mv/vv
+
+    return d
+
 
 from scipy.sparse.linalg import LinearOperator
 

SimPEG/Utils/meshutils.py

@@ -3,6 +3,7 @@ from scipy import sparse as sp
 from matutils import mkvc, ndgrid, sub2ind, sdiag
 from codeutils import asArray_N_x_Dim
 from codeutils import isScalar
+import os
 
 def exampleLrmGrid(nC, exType):
     assert type(nC) == list, "nC must be a list containing the number of nodes"
@@ -142,11 +143,16 @@ def readUBCTensorMesh(fileName):
     tensMsh = Mesh.TensorMesh([h1,h2,h3],x0)
     return tensMsh
 
-
 def readUBCTensorModel(fileName, mesh):
     """
-        ReadUBC 3DTensor mesh model and generate 3D Tensor mesh model in simpegTD
+        Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
 
+        Input:
+        :param fileName, path to the UBC GIF mesh file to read
+        :param mesh, TensorMesh object, mesh that coresponds to the model 
+
+        Output:
+        :return numpy array, model with TensorMesh ordered
     """
     f = open(fileName, 'r')
     model = np.array(map(float, f.readlines()))
@@ -158,19 +164,20 @@ def readUBCTensorModel(fileName, mesh):
 
     return model
 
-def writeUBCTensorMesh(mesh, fileName):
+def writeUBCTensorMesh(fileName, mesh):
     """
         Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
 
-        :param simpeg.Mesh.TensorMesh mesh: The mesh
         :param str fileName: File to write to
+        :param simpeg.Mesh.TensorMesh mesh: The mesh
+        
     """
     assert mesh.dim == 3
     s = ''
     s += '%i %i %i\n' %tuple(mesh.vnC)
-    origin = mesh.x0
+    origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
     origin.dtype = float
-    origin[2] = origin[2]+mesh.hz.sum()
+
     s += '%.2f %.2f %.2f\n' %tuple(origin)
     s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
     s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
@@ -179,14 +186,14 @@ def writeUBCTensorMesh(mesh, fileName):
     f.write(s)
     f.close()
 
-def writeUBCTensorModel(mesh, model, fileName):
+def writeUBCTensorModel(fileName, mesh, model):
     """
         Writes a model associated with a SimPEG TensorMesh
         to a UBC-GIF format model file.
 
+        :param str fileName: File to write to
         :param simpeg.Mesh.TensorMesh mesh: The mesh
         :param numpy.ndarray model: The model
-        :param str fileName: File to write to
     """
 
     # Reshape model to a matrix
@@ -199,6 +206,198 @@ def writeUBCTensorModel(mesh, model, fileName):
     np.savetxt(fileName, modelMatTR.ravel())
 
 
+def readVTRFile(fileName):
+    """
+        Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
+
+        Input:
+        :param vtrFileName, path to the vtr model file to write to
+
+        Output:
+        :return SimPEG TensorMesh object
+        :return SimPEG model dictionary
+        
+    """
+    # Import
+    from vtk import vtkXMLRectilinearGridReader as vtrFileReader
+    from vtk.util.numpy_support import vtk_to_numpy
+
+    # Read the file
+    vtrReader = vtrFileReader()
+    vtrReader.SetFileName(fileName)
+    vtrReader.Update()
+    vtrGrid = vtrReader.GetOutput()
+    # Sort information
+    hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates())))
+    xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0]
+    hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates())))
+    yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0]
+    zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates()))
+    # Check the direction of hz
+    if np.all(zD < 0):
+        hz = np.abs(zD[::-1])
+        zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1]
+    else:
+        hz = np.abs(zD)
+        zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0]
+    x0 = np.array([xR,yR,zR])
+
+    # Make the SimPEG object
+    from SimPEG import Mesh
+    tensMsh = Mesh.TensorMesh([hx,hy,hz],x0)
+
+    # Grap the models
+    modelDict = {}
+    for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()):
+        modelName = vtrGrid.GetCellData().GetArrayName(i)
+        if np.all(zD < 0):
+            modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
+            tM = tensMsh.r(modFlip,'CC','CC','M')
+            modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V')
+        else:
+            modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i))
+        modelDict[modelName] = modArr
+
+    # Return the data
+    return tensMsh, modelDict
+
+def writeVTRFile(fileName,mesh,model=None):
+    """
+    Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
+
+    Input:
+    :param str, path to the output vtk file
+    :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+    :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells
+
+    """
+    # Import
+    from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter
+    from vtk.util.numpy_support import numpy_to_vtk
+
+    # Deal with dimensionalities
+    if mesh.dim >= 1:
+        vX = mesh.vectorNx
+        xD = mesh.nNx
+        yD,zD = 1,1
+        vY, vZ = np.array([0,0])
+    if mesh.dim >= 2:
+        vY = mesh.vectorNy
+        yD = mesh.nNy
+    if mesh.dim == 3:
+        vZ = mesh.vectorNz
+        zD = mesh.nNz
+    # Use rectilinear VTK grid.
+    # Assign the spatial information.
+    vtkObj = rectGrid()
+    vtkObj.SetDimensions(xD,yD,zD)
+    vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
+    vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
+    vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
+
+    # Assign the model('s) to the object
+    for item in model.iteritems():
+        # Convert numpy array
+        vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
+        vtkDoubleArr.SetName(item[0])
+        vtkObj.GetCellData().AddArray(vtkDoubleArr)
+    # Set the active scalar
+    vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+    vtkObj.Update()
+
+
+    # Check the extension of the fileName
+    ext = os.path.splitext(fileName)[1]
+    if ext is '':
+        fileName = fileName + '.vtr'
+    elif ext not in '.vtr':
+        raise IOError('{:s} is an incorrect extension, has to be .vtr')
+    # Write the file.
+    vtrWriteFilter = rectWriter()
+    vtrWriteFilter.SetInput(vtkObj)
+    vtrWriteFilter.SetFileName(fileName)
+    vtrWriteFilter.Update()
+
+
+def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
+    """
+        Extracts Core Mesh from Global mesh
+        xyzlim: 2D array [ndim x 2]
+        mesh: SimPEG mesh
+        This function ouputs: 
+            - actind: corresponding boolean index from global to core
+            - meshcore: core SimPEG mesh  
+        Warning: 1D and 2D has not been tested
+    """
+    from SimPEG import Mesh
+    if mesh.dim ==1:
+        xyzlim = xyzlim.flatten()
+        xmin, xmax = xyzlim[0], xyzlim[1]
+        
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)        
+        
+        xc = mesh.vectorCCx[xind]
+
+        hx = mesh.hx[xind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
+        
+    elif mesh.dim ==2:
+        xmin, xmax = xyzlim[0,0], xyzlim[0,1]
+        ymin, ymax = xyzlim[1,0], xyzlim[1,1]
+
+        yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+
+        xc = mesh.vectorCCx[xind]
+        yc = mesh.vectorCCy[yind]
+
+        hx = mesh.hx[xind]
+        hy = mesh.hy[yind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
+               & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
+        
+    elif mesh.dim==3:
+        xmin, xmax = xyzlim[0,0], xyzlim[0,1]
+        ymin, ymax = xyzlim[1,0], xyzlim[1,1]
+        zmin, zmax = xyzlim[2,0], xyzlim[2,1]
+        
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
+        yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+
+        xc = mesh.vectorCCx[xind]
+        yc = mesh.vectorCCy[yind]
+        zc = mesh.vectorCCz[zind]
+
+        hx = mesh.hx[xind]
+        hy = mesh.hy[yind]
+        hz = mesh.hz[zind]
+        
+        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
+        
+        meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
+        
+        actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
+               & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
+               & (mesh.gridCC[:,2]>zmin) & (mesh.gridCC[:,2]<zmax)
+                
+    else:
+        raise(Exception("Not implemented!"))
+    
+    
+    return actind, meshCore
+
+
 if __name__ == '__main__':
     from SimPEG import Mesh
     import matplotlib.pyplot as plt

SimPEG/__init__.py

@@ -15,4 +15,7 @@ import Directives
 import Inversion
 import Tests
 
-__version__ = '0.1.1'
+__version__   = '0.1.3'
+__author__    = 'Rowan Cockett'
+__license__   = 'MIT'
+__copyright__ = 'Copyright 2014 Rowan Cockett'

docs/api_Examples.rst

@@ -0,0 +1,11 @@
+.. _api_Examples:
+
+Examples
+********
+
+Forward problem
+===============
+
+* `Example 1: Direct Current <http://www.seogi.me/s/notebooks/DCEx.html>`_
+* `Example 2: Seismic-Acoustic <http://www.seogi.me/s/notebooks/SeismicEx.html>`_
+

docs/api_InnerProducts.rst

@@ -75,7 +75,7 @@ We multiply by  square-root of volume on each side of the tensor conductivity to
 
 .. math::
     \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
-    \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM}
+    \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
 
 Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
 We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
@@ -114,7 +114,7 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of
 
 .. math::
 
-    \mathbf{P}_{(i)} =  \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)}
+    \mathbf{P}_{(i)} =  \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)}
 
 .. note::
 

docs/api_Maps.rst

@@ -119,7 +119,6 @@ When these are used in the inverse problem, this is extremely important!!
     expMap.test(m, plotIt=True)
 
 
-
 The API
 =======
 
@@ -154,6 +153,14 @@ Vertical 1D Map
     :undoc-members:
 
 
+Map 2D Cross-Section to 3D Model
+--------------------------------
+
+.. autoclass:: SimPEG.Maps.Map2Dto3D
+    :members:
+    :undoc-members:
+
+
 Mesh to Mesh Map
 ----------------
 
@@ -170,8 +177,8 @@ Mesh to Mesh Map
     v = Utils.mkvc(V)
     modh = Maps.Mesh2Mesh([M,M2])
     modH = Maps.Mesh2Mesh([M2,M])
-    H = modH.transform(v)
-    h = modh.transform(H)
+    H = modH * v
+    h = modh * H
     ax = plt.subplot(131)
     M.plotImage(v, ax=ax)
     ax.set_title('Fine Mesh (Original)')
@@ -196,7 +203,7 @@ Combo Map
 ---------
 
 The ComboMap holds the information for multiplying and combining
-maps. It also uses the chain rule create the derivative.
+maps. It also uses the chain rule to create the derivative.
 Remember, any time that you make your own combination of mappings
 be sure to test that the derivative is correct.
 
@@ -204,10 +211,3 @@ be sure to test that the derivative is correct.
     :members:
     :undoc-members:
 
-
-Non Linear Map
---------------
-
-.. autoclass:: SimPEG.Maps.NonLinearMap
-    :members:
-    :undoc-members:

docs/api_Mesh.rst

@@ -29,7 +29,7 @@ the implementations.
     tM = Mesh.TensorMesh(sz)
     qM = Mesh.TreeMesh(sz)
     qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0)
-    rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
 
     fig, axes = plt.subplots(1,3,figsize=(14,5))
     opts = {}
@@ -38,7 +38,7 @@ the implementations.
     qM.plotGrid(ax=axes[1], **opts)
     axes[1].set_title('TreeMesh')
     rM.plotGrid(ax=axes[2], **opts)
-    axes[2].set_title('LogicallyRectMesh')
+    axes[2].set_title('CurvilinearMesh')
     plt.show()
 
 

docs/api_MeshCode.rst

@@ -18,10 +18,10 @@ Tree Mesh
     :undoc-members:
 
 
-Logically Rectangular Mesh
-==========================
+Curvilinear Mesh
+================
 
-.. automodule:: SimPEG.Mesh.LogicallyRectMesh
+.. automodule:: SimPEG.Mesh.Curvilinear
     :show-inheritance:
     :members:
     :undoc-members:

docs/api_Utils.rst

@@ -23,10 +23,10 @@ Solver Utilities
     :members:
     :undoc-members:
 
-LRM Utilities
+Curv Utilities
 =============
 
-.. automodule:: SimPEG.Utils.lrmutils
+.. automodule:: SimPEG.Utils.curvutils
     :members:
     :undoc-members:
 

docs/api_bigPicture.rst

@@ -7,7 +7,7 @@ Why SimPEG?
 The Big Picture
 ===============
 
-.. image:: simpeg-framework.png
+.. image:: InversionWorkflow.png
    :width: 400 px
    :alt: Framework
    :align: center

docs/api_installing.rst

@@ -25,7 +25,10 @@ You can download the package manager and use it to install the dependencies abov
 Installing SimPEG
 =================
 
-SimPEG is not **yet** on pip.
+SimPEG is on pip!!::
+
+    pip install SimPEG
+
 
 Installing from Source
 ----------------------

docs/conf.py

@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.0.1'
+version = '0.1.3'
 # The full version, including alpha/beta/rc tags.
-release = '0.0.1'
+release = '0.1.3'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -94,7 +94,13 @@ pygments_style = 'sphinx'
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'default'
+try:
+    import sphinx_rtd_theme
+    html_theme = 'sphinx_rtd_theme'
+    html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
+    pass
+except Exception, e:
+    html_theme = 'default'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the

docs/index.rst

@@ -102,3 +102,10 @@ Project Index & Search
 * :ref:`modindex`
 * :ref:`search`
 
+Examples
+********
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Examples

ez_setup.py

@@ -1,485 +0,0 @@
-#!python
-"""Bootstrap distribute installation
-
-If you want to use setuptools in your package's setup.py, just include this
-file in the same directory with it, and add this to the top of your setup.py::
-
-    from distribute_setup import use_setuptools
-    use_setuptools()
-
-If you want to require a specific version of setuptools, set a download
-mirror, or use an alternate download directory, you can do so by supplying
-the appropriate options to ``use_setuptools()``.
-
-This file can also be run as a script to install or upgrade setuptools.
-"""
-import os
-import sys
-import time
-import fnmatch
-import tempfile
-import tarfile
-from distutils import log
-
-try:
-    from site import USER_SITE
-except ImportError:
-    USER_SITE = None
-
-try:
-    import subprocess
-
-    def _python_cmd(*args):
-        args = (sys.executable,) + args
-        return subprocess.call(args) == 0
-
-except ImportError:
-    # will be used for python 2.3
-    def _python_cmd(*args):
-        args = (sys.executable,) + args
-        # quoting arguments if windows
-        if sys.platform == 'win32':
-            def quote(arg):
-                if ' ' in arg:
-                    return '"%s"' % arg
-                return arg
-            args = [quote(arg) for arg in args]
-        return os.spawnl(os.P_WAIT, sys.executable, *args) == 0
-
-DEFAULT_VERSION = "0.6.14"
-DEFAULT_URL = "http://pypi.python.org/packages/source/d/distribute/"
-SETUPTOOLS_FAKED_VERSION = "0.6c11"
-
-SETUPTOOLS_PKG_INFO = """\
-Metadata-Version: 1.0
-Name: setuptools
-Version: %s
-Summary: xxxx
-Home-page: xxx
-Author: xxx
-Author-email: xxx
-License: xxx
-Description: xxx
-""" % SETUPTOOLS_FAKED_VERSION
-
-
-def _install(tarball):
-    # extracting the tarball
-    tmpdir = tempfile.mkdtemp()
-    log.warn('Extracting in %s', tmpdir)
-    old_wd = os.getcwd()
-    try:
-        os.chdir(tmpdir)
-        tar = tarfile.open(tarball)
-        _extractall(tar)
-        tar.close()
-
-        # going in the directory
-        subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
-        os.chdir(subdir)
-        log.warn('Now working in %s', subdir)
-
-        # installing
-        log.warn('Installing Distribute')
-        if not _python_cmd('setup.py', 'install'):
-            log.warn('Something went wrong during the installation.')
-            log.warn('See the error message above.')
-    finally:
-        os.chdir(old_wd)
-
-
-def _build_egg(egg, tarball, to_dir):
-    # extracting the tarball
-    tmpdir = tempfile.mkdtemp()
-    log.warn('Extracting in %s', tmpdir)
-    old_wd = os.getcwd()
-    try:
-        os.chdir(tmpdir)
-        tar = tarfile.open(tarball)
-        _extractall(tar)
-        tar.close()
-
-        # going in the directory
-        subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0])
-        os.chdir(subdir)
-        log.warn('Now working in %s', subdir)
-
-        # building an egg
-        log.warn('Building a Distribute egg in %s', to_dir)
-        _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir)
-
-    finally:
-        os.chdir(old_wd)
-    # returning the result
-    log.warn(egg)
-    if not os.path.exists(egg):
-        raise IOError('Could not build the egg.')
-
-
-def _do_download(version, download_base, to_dir, download_delay):
-    egg = os.path.join(to_dir, 'distribute-%s-py%d.%d.egg'
-                       % (version, sys.version_info[0], sys.version_info[1]))
-    if not os.path.exists(egg):
-        tarball = download_setuptools(version, download_base,
-                                      to_dir, download_delay)
-        _build_egg(egg, tarball, to_dir)
-    sys.path.insert(0, egg)
-    import setuptools
-    setuptools.bootstrap_install_from = egg
-
-
-def use_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
-                   to_dir=os.curdir, download_delay=15, no_fake=True):
-    # making sure we use the absolute path
-    to_dir = os.path.abspath(to_dir)
-    was_imported = 'pkg_resources' in sys.modules or \
-        'setuptools' in sys.modules
-    try:
-        try:
-            import pkg_resources
-            if not hasattr(pkg_resources, '_distribute'):
-                if not no_fake:
-                    _fake_setuptools()
-                raise ImportError
-        except ImportError:
-            return _do_download(version, download_base, to_dir, download_delay)
-        try:
-            pkg_resources.require("distribute>="+version)
-            return
-        except pkg_resources.VersionConflict:
-            e = sys.exc_info()[1]
-            if was_imported:
-                sys.stderr.write(
-                "The required version of distribute (>=%s) is not available,\n"
-                "and can't be installed while this script is running. Please\n"
-                "install a more recent version first, using\n"
-                "'easy_install -U distribute'."
-                "\n\n(Currently using %r)\n" % (version, e.args[0]))
-                sys.exit(2)
-            else:
-                del pkg_resources, sys.modules['pkg_resources']    # reload ok
-                return _do_download(version, download_base, to_dir,
-                                    download_delay)
-        except pkg_resources.DistributionNotFound:
-            return _do_download(version, download_base, to_dir,
-                                download_delay)
-    finally:
-        if not no_fake:
-            _create_fake_setuptools_pkg_info(to_dir)
-
-def download_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL,
-                        to_dir=os.curdir, delay=15):
-    """Download distribute from a specified location and return its filename
-
-    `version` should be a valid distribute version number that is available
-    as an egg for download under the `download_base` URL (which should end
-    with a '/'). `to_dir` is the directory where the egg will be downloaded.
-    `delay` is the number of seconds to pause before an actual download
-    attempt.
-    """
-    # making sure we use the absolute path
-    to_dir = os.path.abspath(to_dir)
-    try:
-        from urllib.request import urlopen
-    except ImportError:
-        from urllib2 import urlopen
-    tgz_name = "distribute-%s.tar.gz" % version
-    url = download_base + tgz_name
-    saveto = os.path.join(to_dir, tgz_name)
-    src = dst = None
-    if not os.path.exists(saveto):  # Avoid repeated downloads
-        try:
-            log.warn("Downloading %s", url)
-            src = urlopen(url)
-            # Read/write all in one block, so we don't create a corrupt file
-            # if the download is interrupted.
-            data = src.read()
-            dst = open(saveto, "wb")
-            dst.write(data)
-        finally:
-            if src:
-                src.close()
-            if dst:
-                dst.close()
-    return os.path.realpath(saveto)
-
-def _no_sandbox(function):
-    def __no_sandbox(*args, **kw):
-        try:
-            from setuptools.sandbox import DirectorySandbox
-            if not hasattr(DirectorySandbox, '_old'):
-                def violation(*args):
-                    pass
-                DirectorySandbox._old = DirectorySandbox._violation
-                DirectorySandbox._violation = violation
-                patched = True
-            else:
-                patched = False
-        except ImportError:
-            patched = False
-
-        try:
-            return function(*args, **kw)
-        finally:
-            if patched:
-                DirectorySandbox._violation = DirectorySandbox._old
-                del DirectorySandbox._old
-
-    return __no_sandbox
-
-def _patch_file(path, content):
-    """Will backup the file then patch it"""
-    existing_content = open(path).read()
-    if existing_content == content:
-        # already patched
-        log.warn('Already patched.')
-        return False
-    log.warn('Patching...')
-    _rename_path(path)
-    f = open(path, 'w')
-    try:
-        f.write(content)
-    finally:
-        f.close()
-    return True
-
-_patch_file = _no_sandbox(_patch_file)
-
-def _same_content(path, content):
-    return open(path).read() == content
-
-def _rename_path(path):
-    new_name = path + '.OLD.%s' % time.time()
-    log.warn('Renaming %s into %s', path, new_name)
-    os.rename(path, new_name)
-    return new_name
-
-def _remove_flat_installation(placeholder):
-    if not os.path.isdir(placeholder):
-        log.warn('Unkown installation at %s', placeholder)
-        return False
-    found = False
-    for file in os.listdir(placeholder):
-        if fnmatch.fnmatch(file, 'setuptools*.egg-info'):
-            found = True
-            break
-    if not found:
-        log.warn('Could not locate setuptools*.egg-info')
-        return
-
-    log.warn('Removing elements out of the way...')
-    pkg_info = os.path.join(placeholder, file)
-    if os.path.isdir(pkg_info):
-        patched = _patch_egg_dir(pkg_info)
-    else:
-        patched = _patch_file(pkg_info, SETUPTOOLS_PKG_INFO)
-
-    if not patched:
-        log.warn('%s already patched.', pkg_info)
-        return False
-    # now let's move the files out of the way
-    for element in ('setuptools', 'pkg_resources.py', 'site.py'):
-        element = os.path.join(placeholder, element)
-        if os.path.exists(element):
-            _rename_path(element)
-        else:
-            log.warn('Could not find the %s element of the '
-                     'Setuptools distribution', element)
-    return True
-
-_remove_flat_installation = _no_sandbox(_remove_flat_installation)
-
-def _after_install(dist):
-    log.warn('After install bootstrap.')
-    placeholder = dist.get_command_obj('install').install_purelib
-    _create_fake_setuptools_pkg_info(placeholder)
-
-def _create_fake_setuptools_pkg_info(placeholder):
-    if not placeholder or not os.path.exists(placeholder):
-        log.warn('Could not find the install location')
-        return
-    pyver = '%s.%s' % (sys.version_info[0], sys.version_info[1])
-    setuptools_file = 'setuptools-%s-py%s.egg-info' % \
-            (SETUPTOOLS_FAKED_VERSION, pyver)
-    pkg_info = os.path.join(placeholder, setuptools_file)
-    if os.path.exists(pkg_info):
-        log.warn('%s already exists', pkg_info)
-        return
-
-    log.warn('Creating %s', pkg_info)
-    f = open(pkg_info, 'w')
-    try:
-        f.write(SETUPTOOLS_PKG_INFO)
-    finally:
-        f.close()
-
-    pth_file = os.path.join(placeholder, 'setuptools.pth')
-    log.warn('Creating %s', pth_file)
-    f = open(pth_file, 'w')
-    try:
-        f.write(os.path.join(os.curdir, setuptools_file))
-    finally:
-        f.close()
-
-_create_fake_setuptools_pkg_info = _no_sandbox(_create_fake_setuptools_pkg_info)
-
-def _patch_egg_dir(path):
-    # let's check if it's already patched
-    pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO')
-    if os.path.exists(pkg_info):
-        if _same_content(pkg_info, SETUPTOOLS_PKG_INFO):
-            log.warn('%s already patched.', pkg_info)
-            return False
-    _rename_path(path)
-    os.mkdir(path)
-    os.mkdir(os.path.join(path, 'EGG-INFO'))
-    pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO')
-    f = open(pkg_info, 'w')
-    try:
-        f.write(SETUPTOOLS_PKG_INFO)
-    finally:
-        f.close()
-    return True
-
-_patch_egg_dir = _no_sandbox(_patch_egg_dir)
-
-def _before_install():
-    log.warn('Before install bootstrap.')
-    _fake_setuptools()
-
-
-def _under_prefix(location):
-    if 'install' not in sys.argv:
-        return True
-    args = sys.argv[sys.argv.index('install')+1:]
-    for index, arg in enumerate(args):
-        for option in ('--root', '--prefix'):
-            if arg.startswith('%s=' % option):
-                top_dir = arg.split('root=')[-1]
-                return location.startswith(top_dir)
-            elif arg == option:
-                if len(args) > index:
-                    top_dir = args[index+1]
-                    return location.startswith(top_dir)
-        if arg == '--user' and USER_SITE is not None:
-            return location.startswith(USER_SITE)
-    return True
-
-
-def _fake_setuptools():
-    log.warn('Scanning installed packages')
-    try:
-        import pkg_resources
-    except ImportError:
-        # we're cool
-        log.warn('Setuptools or Distribute does not seem to be installed.')
-        return
-    ws = pkg_resources.working_set
-    try:
-        setuptools_dist = ws.find(pkg_resources.Requirement.parse('setuptools',
-                                  replacement=False))
-    except TypeError:
-        # old distribute API
-        setuptools_dist = ws.find(pkg_resources.Requirement.parse('setuptools'))
-
-    if setuptools_dist is None:
-        log.warn('No setuptools distribution found')
-        return
-    # detecting if it was already faked
-    setuptools_location = setuptools_dist.location
-    log.warn('Setuptools installation detected at %s', setuptools_location)
-
-    # if --root or --preix was provided, and if
-    # setuptools is not located in them, we don't patch it
-    if not _under_prefix(setuptools_location):
-        log.warn('Not patching, --root or --prefix is installing Distribute'
-                 ' in another location')
-        return
-
-    # let's see if its an egg
-    if not setuptools_location.endswith('.egg'):
-        log.warn('Non-egg installation')
-        res = _remove_flat_installation(setuptools_location)
-        if not res:
-            return
-    else:
-        log.warn('Egg installation')
-        pkg_info = os.path.join(setuptools_location, 'EGG-INFO', 'PKG-INFO')
-        if (os.path.exists(pkg_info) and
-            _same_content(pkg_info, SETUPTOOLS_PKG_INFO)):
-            log.warn('Already patched.')
-            return
-        log.warn('Patching...')
-        # let's create a fake egg replacing setuptools one
-        res = _patch_egg_dir(setuptools_location)
-        if not res:
-            return
-    log.warn('Patched done.')
-    _relaunch()
-
-
-def _relaunch():
-    log.warn('Relaunching...')
-    # we have to relaunch the process
-    # pip marker to avoid a relaunch bug
-    if sys.argv[:3] == ['-c', 'install', '--single-version-externally-managed']:
-        sys.argv[0] = 'setup.py'
-    args = [sys.executable] + sys.argv
-    sys.exit(subprocess.call(args))
-
-
-def _extractall(self, path=".", members=None):
-    """Extract all members from the archive to the current working
-       directory and set owner, modification time and permissions on
-       directories afterwards. `path' specifies a different directory
-       to extract to. `members' is optional and must be a subset of the
-       list returned by getmembers().
-    """
-    import copy
-    import operator
-    from tarfile import ExtractError
-    directories = []
-
-    if members is None:
-        members = self
-
-    for tarinfo in members:
-        if tarinfo.isdir():
-            # Extract directories with a safe mode.
-            directories.append(tarinfo)
-            tarinfo = copy.copy(tarinfo)
-            tarinfo.mode = 448 # decimal for oct 0700
-        self.extract(tarinfo, path)
-
-    # Reverse sort directories.
-    if sys.version_info < (2, 4):
-        def sorter(dir1, dir2):
-            return cmp(dir1.name, dir2.name)
-        directories.sort(sorter)
-        directories.reverse()
-    else:
-        directories.sort(key=operator.attrgetter('name'), reverse=True)
-
-    # Set correct owner, mtime and filemode on directories.
-    for tarinfo in directories:
-        dirpath = os.path.join(path, tarinfo.name)
-        try:
-            self.chown(tarinfo, dirpath)
-            self.utime(tarinfo, dirpath)
-            self.chmod(tarinfo, dirpath)
-        except ExtractError:
-            e = sys.exc_info()[1]
-            if self.errorlevel > 1:
-                raise
-            else:
-                self._dbg(1, "tarfile: %s" % e)
-
-
-def main(argv, version=DEFAULT_VERSION):
-    """Install or upgrade setuptools and EasyInstall"""
-    tarball = download_setuptools()
-    _install(tarball)
-
-
-if __name__ == '__main__':
-    main(sys.argv[1:])

setup.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-"""SimPEG: Simulation and Parameter Estimation for Geophysics
+"""SimPEG: Simulation and Parameter Estimation in Geophysics
 
 SimPEG is a python package for simulation and gradient based
 parameter estimation in the context of geophysical applications.
@@ -11,33 +11,40 @@ from Cython.Build import cythonize
 import numpy as np
 
 CLASSIFIERS = [
+'Development Status :: 4 - Beta',
 'Intended Audience :: Developers',
+'Intended Audience :: Science/Research',
 'License :: OSI Approved :: MIT License',
 'Programming Language :: Python',
 'Topic :: Scientific/Engineering',
 'Topic :: Scientific/Engineering :: Mathematics',
+'Topic :: Scientific/Engineering :: Physics',
 'Operating System :: Microsoft :: Windows',
 'Operating System :: POSIX',
 'Operating System :: Unix',
 'Operating System :: MacOS',
+'Natural Language :: English',
 ]
 
 import os, os.path
 
+with open("README.rst") as f:
+    LONG_DESCRIPTION = ''.join(f.readlines())
+
 setup(
     name = "SimPEG",
-    version = "0.1.1",
+    version = "0.1.3",
     packages = find_packages(),
     install_requires = ['numpy>=1.7',
-                        'scipy>=0.13',
-                        'matplotlib>=1.3',
+                        'scipy>=0.13'
                        ],
     author = "Rowan Cockett",
-    author_email = "rowanc1@gmail.com",
-    description = "SimPEG: Simulation and Parameter Estimation for Geophysics",
+    author_email = "rowan@3ptscience.com",
+    description = "SimPEG: Simulation and Parameter Estimation in Geophysics",
+    long_description = LONG_DESCRIPTION,
     license = "MIT",
     keywords = "geophysics inverse problem",
-    url = "http://simpeg.3ptscience.com/",
+    url = "http://simpeg.xyz/",
     download_url = "http://github.com/simpeg/simpeg",
     classifiers=CLASSIFIERS,
     platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],