Fixed bugs in meshutils, writing VTR files

Allow fileName to be None, which will output the VTKobject without save a file.

.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

@@ -14,6 +14,34 @@ class BaseDataMisfit(object):
     debug   = False  #: Print debugging information
     counter = None   #: Set this to a SimPEG.Utils.Counter() if you want to count things
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     def __init__(self, survey, **kwargs):
         assert survey.ispaired, 'The survey must be paired to a problem.'
         if isinstance(survey, Survey.BaseSurvey):
@@ -25,7 +53,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 +65,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 +78,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 +87,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

@@ -8,6 +8,34 @@ class InversionDirective(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def inversion(self):
         """This is the inversion of the InversionDirective instance."""
@@ -88,7 +116,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 +164,7 @@ class BetaEstimate_ByEig(InversionDirective):
 class BetaSchedule(InversionDirective):
     """BetaSchedule"""
 
-    coolingFactor = 2.
+    coolingFactor = 8.
     coolingRate = 3
 
     def endIter(self):
@@ -144,6 +172,140 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
+class TargetMisfit(InversionDirective):
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
+
+    def endIter(self):
+        if self.invProb.phi_d < self.target:
+            self.opt.stopNextIteration = True
+
+
+
+class _SaveEveryIteration(InversionDirective):
+    @property
+    def name(self):
+        if getattr(self, '_name', None) is None:
+            self._name = 'InversionModel'
+        return self._name
+    @name.setter
+    def name(self, value):
+        self._name = value
+
+    @property
+    def fileName(self):
+        if getattr(self, '_fileName', None) is None:
+            from datetime import datetime
+            self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
+        return self._fileName
+    @fileName.setter
+    def fileName(self, value):
+        self._fileName = value
+
+
+class SaveModelEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
+
+    def endIter(self):
+        np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
+
+
+class SaveOutputEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
+        f = open(self.fileName+'.txt', 'w')
+        f.write("  #     beta     phi_d     phi_m       f\n")
+        f.close()
+
+    def endIter(self):
+        f = open(self.fileName+'.txt', 'a')
+        f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
+        f.close()
+
+
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel)
+
+
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration
+
+    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
+
+
+
+    """
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+
+
 
 # class UpdateReferenceModel(Parameter):
 

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/Fields.py

@@ -0,0 +1,289 @@
+import Utils, numpy as np, scipy.sparse as sp
+
+class Fields(object):
+    """Fancy Field Storage
+
+        u[:,'phi'] = phi
+        print u[src0,'phi']
+
+    """
+
+    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
+    aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
+    dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    def __init__(self, mesh, survey, **kwargs):
+        self.survey = survey
+        self.mesh = mesh
+        Utils.setKwargs(self, **kwargs)
+        self._fields = {}
+
+        if self.knownFields is None:
+            raise Exception('knownFields cannot be set to None')
+        if self.aliasFields is None:
+            self.aliasFields = {}
+
+        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
+        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
+        self.startup()
+
+    def startup(self):
+        pass
+
+    @property
+    def approxSize(self):
+        """The approximate cost to storing all of the known fields."""
+        sz = 0.0
+        for f in self.knownFields:
+            loc =self.knownFields[f]
+            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
+        return "%e MB"%sz
+
+    def _storageShape(self, loc):
+        nSrc = self.survey.nSrc
+
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+
+        return (nP, nSrc)
+
+    def _initStore(self, name):
+        if name in self._fields:
+            return self._fields[name]
+
+        assert name in self.knownFields, 'field name is not known.'
+
+        loc = self.knownFields[name]
+
+        if type(self.dtype) is dict:
+            dtype = self.dtype[name]
+        else:
+            dtype = self.dtype
+        field = np.zeros(self._storageShape(loc), dtype=dtype)
+
+        self._fields[name] = field
+
+        return field
+
+    def _srcIndex(self, srcTestList):
+        if type(srcTestList) is slice:
+            ind = srcTestList
+        else:
+            ind = self.survey.getSourceIndex(srcTestList)
+        return ind
+
+    def _nameIndex(self, name, accessType):
+
+        if type(name) is slice:
+            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
+            name = None
+
+        if name is None:
+            return
+        if accessType=='set' and name not in self.knownFields:
+            if name in self.aliasFields:
+                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
+            else:
+                raise KeyError('Invalid field name (%s) for setter'%name)
+
+        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
+            raise KeyError('Invalid field name (%s) for getter'%name)
+        return name
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+
+        assert len(key) == 2, 'must be [Src, fieldName]'
+
+        srcTestList, name = key
+        name = self._nameIndex(name, accessType)
+        ind = self._srcIndex(srcTestList)
+        return ind, name
+
+    def __setitem__(self, key, value):
+        ind, name = self._indexAndNameFromKey(key, 'set')
+        if name is None:
+            freq = key
+            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
+            newFields = value
+        elif name in self.knownFields:
+            newFields = {name: value}
+        else:
+            raise Exception('Unknown setter')
+
+        for name in newFields:
+            field = self._initStore(name)
+            self._setField(field, newFields[name], name, ind)
+
+    def __getitem__(self, key):
+        ind, name = self._indexAndNameFromKey(key, 'get')
+        if name is None:
+            out = {}
+            for name in self._fields:
+                out[name] = self._getField(name, ind)
+            return out
+        return self._getField(name, ind)
+
+    def _setField(self, field, val, name, ind):
+        if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1):
+            val = Utils.mkvc(val,2)
+        field[:,ind] = val
+
+    def _getField(self, name, ind):
+        if name in self._fields:
+            out = self._fields[name][:,ind]
+        else:
+            # Aliased fields
+            alias, loc, func = self.aliasFields[name]
+
+            srcII = np.array(self.survey.srcList)[ind]
+            srcII = srcII.tolist()
+
+            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], srcII)
+        if out.shape[0] == out.size or out.ndim == 1:
+            out = Utils.mkvc(out,2)
+        return out
+
+    def __contains__(self, other):
+        if other in self.aliasFields:
+            other = self.aliasFields[other][0]
+        return self._fields.__contains__(other)
+
+
+class TimeFields(Fields):
+    """Fancy Field Storage for time domain problems
+
+        u[:,'phi', timeInd] = phi
+        print u[src0,'phi']
+
+    """
+
+    def _storageShape(self, loc):
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+        nSrc = self.survey.nSrc
+        nT = self.survey.prob.nT + 1
+        return (nP, nSrc, nT)
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+        if len(key) == 2:
+            key += (slice(None,None,None),)
+
+        assert len(key) == 3, 'must be [Src, fieldName, times]'
+
+        srcTestList, name, timeInd = key
+
+        name = self._nameIndex(name, accessType)
+        srcInd = self._srcIndex(srcTestList)
+
+        return (srcInd, timeInd), name
+
+    def _correctShape(self, name, ind, deflate=False):
+        srcInd, timeInd = ind
+        if name in self.knownFields:
+            loc = self.knownFields[name]
+        else:
+            loc = self.aliasFields[name][1]
+        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, 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):
+        srcInd, timeInd = ind
+        shape = self._correctShape(name, ind)
+        if Utils.isScalar(val):
+            field[:,srcInd,timeInd] = val
+            return
+        if val.size != np.array(shape).prod():
+            raise ValueError('Incorrect size for data.')
+        correctShape = field[:,srcInd,timeInd].shape
+        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
+
+    def _getField(self, name, ind):
+        srcInd, timeInd = ind
+
+        if name in self._fields:
+            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][:,srcInd,timeInd]
+            pointerShape = self._correctShape(alias, ind)
+            pointerFields = pointerFields.reshape(pointerShape, order='F')
+
+            timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
+            srcII  = np.array(self.survey.srcList)[srcInd]
+            srcII  = srcII.tolist()
+
+            if timeII.size == 1:
+                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
+                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
+                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.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:
+                        out[i] = out[i][:,:,np.newaxis]
+                out = np.concatenate(out, axis=2)
+
+        shape = self._correctShape(name, ind, deflate=True)
+        return out.reshape(shape, order='F')
+

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

@@ -1,5 +1,6 @@
 import Utils, numpy as np, scipy.sparse as sp
 from Tests import checkDerivative
+from PropMaps import PropMap, Property
 
 
 class IdentityMap(object):
@@ -12,15 +13,42 @@ class IdentityMap(object):
 
     mesh = None      #: A SimPEG Mesh
 
-    def __init__(self, mesh):
+    def __init__(self, mesh, **kwargs):
+        Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+
     @property
     def nP(self):
         """
             :rtype: int
             :return: number of parameters in the model
         """
+        if self.mesh is None:
+            return '*'
         return self.mesh.nC
 
     @property
@@ -31,8 +59,11 @@ class IdentityMap(object):
             :rtype: (int,int)
             :return: shape of the operator as a tuple
         """
+        if self.mesh is None:
+            return ('*', self.nP)
         return (self.mesh.nC, self.nP)
 
+
     def _transform(self, m):
         """
             Changes the model into the physical property.
@@ -85,10 +116,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
@@ -97,17 +128,17 @@ class IdentityMap(object):
 
     def __mul__(self, val):
         if isinstance(val, IdentityMap):
-            if not self.shape[1] == val.shape[0]:
+            if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
             return ComboMap([self, val])
         elif isinstance(val, np.ndarray):
-            if not self.shape[1] == val.shape[0]:
+            if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
         raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
 
     def __str__(self):
-        return "%s(%d,%d)" % (self.__class__.__name__, self.shape[0], self.shape[1])
+        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
 
 class ComboMap(IdentityMap):
     """Combination of various maps."""
@@ -118,10 +149,10 @@ class ComboMap(IdentityMap):
         self.maps = []
         for ii, m in enumerate(maps):
             assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!'
-            if ii > 0 and not self.shape[1] == m.shape[0]:
+            if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
                 prev = self.maps[-1]
-                errArgs = (prev.__name__, prev.shape[0], prev.shape[1], m.__name__, m.shape[0], m.shape[1])
-                raise ValueError('Dimension mismatch in map[%s] (%i, %i) and map[%s] (%i, %i).' % errArgs)
+                errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
+                raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
 
             if isinstance(m, ComboMap):
                 self.maps += m.maps
@@ -154,12 +185,11 @@ class ComboMap(IdentityMap):
         return deriv
 
     def __str__(self):
-        return 'ComboMap[%s]%s' % (' * '.join([m.__str__() for m in self.maps]), str(self.shape))
+        return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
 
 
 class ExpMap(IdentityMap):
     """
-
         Changes the model into the physical property.
 
         A common example of this is to invert for electrical conductivity
@@ -220,6 +250,97 @@ class ExpMap(IdentityMap):
         """
         return Utils.sdiag(np.exp(Utils.mkvc(m)))
 
+class ReciprocalMap(IdentityMap):
+    """
+        Reciprocal mapping. For example, electrical resistivity and conductivity.
+
+        .. math::
+
+            \\rho = \\frac{1}{\sigma}
+
+    """
+    def _transform(self, m):
+        return 1.0 / Utils.mkvc(m)
+
+    def inverse(self, D):
+        return 1.0 / Utils.mkvc(m)
+
+    def deriv(self, m):
+        # TODO: if this is a tensor, you might have a problem.
+        return Utils.sdiag( - Utils.mkvc(m)**(-2) )
+
+
+
+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 FullMap(IdentityMap):
+    """
+    FullMap
+
+    Given a scalar, the FullMap maps the value to the
+    full model space.
+    """
+
+    def __init__(self,mesh,**kwargs):
+        IdentityMap.__init__(self, mesh,**kwargs)
+
+    @property
+    def nP(self):
+        return 1
+
+    def _transform(self, m):
+        """
+            :param m: model (scalar)
+            :rtype: numpy.array
+            :return: transformed model
+        """
+        return np.ones(self.mesh.nC)*m
+
+    def deriv(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: derivative of transformed model
+        """
+        return np.ones([self.mesh.nC,1])
+
+
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap
 
@@ -261,6 +382,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 +519,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

@@ -33,7 +33,7 @@ class InnerProducts(object):
         return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
 
     def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
-        """
+        """r
             :param str projType: 'F' for faces 'E' for edges
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool invProp: inverts the material property
@@ -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 = []
@@ -311,6 +328,7 @@ class TensorView(object):
             v = getattr(np,view)(v) # e.g. np.real(v)
             if clim is None:
                 clim = [v.min(),v.max()]
+            v = np.ma.masked_where(np.isnan(v), v)
             out += (ax.pcolormesh(self.vectorNx, self.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
         elif view in ['vec']:
             U, V = self.r(v.reshape((self.nC,-1), order='F'), 'CC', 'CC', 'M')
@@ -518,13 +536,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 +595,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/Models.py

@@ -31,5 +31,5 @@ class Model(np.ndarray):
     @property
     def transformDeriv(self):
         if getattr(self, '_transformDeriv', None) is None:
-            self.deriv = self.mapping.deriv(self.view(np.ndarray))
-        return self.deriv
+            self._transformDeriv = self.mapping.deriv(self.view(np.ndarray))
+        return self._transformDeriv

SimPEG/Optimization.py

@@ -97,6 +97,8 @@ class Minimize(object):
     tolG = 1e-1        #: Tolerance on gradient norm
     eps = 1e-5         #: Small value
 
+    stopNextIteration = False #: Stops the optimization program nicely.
+
     debug   = False    #: Print debugging information
     debugLS = False    #: Print debugging information for the line-search
 
@@ -113,6 +115,34 @@ class Minimize(object):
 
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def callback(self):
         return getattr(self, '_callback', None)
@@ -186,6 +216,7 @@ class Minimize(object):
                 xt, caught = self.modifySearchDirectionBreak(p)
                 if not caught: return self.xc
             self.doEndIteration(xt)
+            if self.stopNextIteration: break
 
         self.printDone()
         self.finish()
@@ -210,6 +241,7 @@ class Minimize(object):
 
         self.iter = 0
         self.iterLS = 0
+        self.stopNextIteration = False
 
         x0 = self.projection(x0)  # ensure that we start of feasible.
         self.x0 = x0
@@ -444,6 +476,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 +962,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 +982,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
             return p
         else:
 
-            
+
             delx = np.zeros(self.g.size)
             resid = -(1-Active) * self.g
 

SimPEG/Problem.py

@@ -1,279 +1,7 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
 import Maps, Mesh
-
-
-class Fields(object):
-    """Fancy Field Storage
-
-        u[:,'phi'] = phi
-        print u[tx0,'phi']
-
-    """
-
-    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
-    aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
-    dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
-
-    def __init__(self, mesh, survey, **kwargs):
-        self.survey = survey
-        self.mesh = mesh
-        Utils.setKwargs(self, **kwargs)
-        self._fields = {}
-
-        if self.knownFields is None:
-            raise Exception('knownFields cannot be set to None')
-        if self.aliasFields is None:
-            self.aliasFields = {}
-
-        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
-        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
-        self.startup()
-
-    def startup(self):
-        pass
-
-    @property
-    def approxSize(self):
-        """The approximate cost to storing all of the known fields."""
-        sz = 0.0
-        for f in self.knownFields:
-            loc =self.knownFields[f]
-            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
-        return "%e MB"%sz
-
-    def _storageShape(self, loc):
-        nTx = self.survey.nTx
-
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-
-        return (nP, nTx)
-
-    def _initStore(self, name):
-        if name in self._fields:
-            return self._fields[name]
-
-        assert name in self.knownFields, 'field name is not known.'
-
-        loc = self.knownFields[name]
-
-        if type(self.dtype) is dict:
-            dtype = self.dtype[name]
-        else:
-            dtype = self.dtype
-        field = np.zeros(self._storageShape(loc), dtype=dtype)
-
-        self._fields[name] = field
-
-        return field
-
-    def _txIndex(self, txTestList):
-        if type(txTestList) is slice:
-            ind = txTestList
-        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.')
-
-            ind = np.in1d(self.survey.txList, txTestList)
-        return ind
-
-    def _nameIndex(self, name, accessType):
-
-        if type(name) is slice:
-            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
-            name = None
-
-        if name is None:
-            return
-        if accessType=='set' and name not in self.knownFields:
-            if name in self.aliasFields:
-                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
-            else:
-                raise KeyError('Invalid field name (%s) for setter'%name)
-
-        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
-            raise KeyError('Invalid field name (%s) for getter'%name)
-        return name
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-
-        assert len(key) == 2, 'must be [Tx, fieldName]'
-
-        txTestList, name = key
-        name = self._nameIndex(name, accessType)
-        ind = self._txIndex(txTestList)
-        return ind, name
-
-    def __setitem__(self, key, value):
-        ind, name = self._indexAndNameFromKey(key, 'set')
-        if name is None:
-            freq = key
-            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
-            newFields = value
-        elif name in self.knownFields:
-            newFields = {name: value}
-        else:
-            raise Exception('Unknown setter')
-
-        for name in newFields:
-            field = self._initStore(name)
-            self._setField(field, newFields[name], name, ind)
-
-    def __getitem__(self, key):
-        ind, name = self._indexAndNameFromKey(key, 'get')
-        if name is None:
-            out = {}
-            for name in self._fields:
-                out[name] = self._getField(name, ind)
-            return out
-        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):
-            val = Utils.mkvc(val,2)
-        field[:,ind] = val
-
-    def _getField(self, name, ind):
-        if name in self._fields:
-            out = self._fields[name][:,ind]
-        else:
-            # 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]
-
-            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)
-        return out
-
-    def __contains__(self, other):
-        if other in self.aliasFields:
-            other = self.aliasFields[other][0]
-        return self._fields.__contains__(other)
-
-
-class TimeFields(Fields):
-    """Fancy Field Storage for time domain problems
-
-        u[:,'phi', timeInd] = phi
-        print u[tx0,'phi']
-
-    """
-
-    def _storageShape(self, loc):
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-        nTx = self.survey.nTx
-        nT = self.survey.prob.nT + 1
-        return (nP, nTx, nT)
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-        if len(key) == 2:
-            key += (slice(None,None,None),)
-
-        assert len(key) == 3, 'must be [Tx, fieldName, times]'
-
-        txTestList, name, timeInd = key
-
-        name = self._nameIndex(name, accessType)
-        txInd = self._txIndex(txTestList)
-
-        return (txInd, timeInd), name
-
-    def _correctShape(self, name, ind, deflate=False):
-        txInd, 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()
-        nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
-        shape = nP, nTx, nT
-        if deflate:
-             shape = tuple([s for s in shape if s > 1])
-        return shape
-
-    def _setField(self, field, val, name, ind):
-        txInd, timeInd = ind
-        shape = self._correctShape(name, ind)
-        if Utils.isScalar(val):
-            field[:,txInd,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')
-
-    def _getField(self, name, ind):
-        txInd, timeInd = ind
-
-        if name in self._fields:
-            out = self._fields[name][:,txInd,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]
-            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]
-
-            if timeII.size == 1:
-                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
-                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
-                out = func(pointerFields, txII, 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 out[i].ndim == 1:
-                        out[i] = out[i][:,np.newaxis,np.newaxis]
-                    elif out[i].ndim == 2:
-                        out[i] = out[i][:,:,np.newaxis]
-                out = np.concatenate(out, axis=2)
-
-        shape = self._correctShape(name, ind, deflate=True)
-        return out.reshape(shape, order='F')
-
-
+from Fields import Fields, TimeFields
 
 class BaseProblem(object):
     """
@@ -290,15 +18,55 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
 
+    PropMap = None    #: A SimPEG PropertyMap class.
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    @property
+    def mapping(self):
+        "A SimPEG.Map instance or a property map is PropMap is not None"
+        return getattr(self, '_mapping', None)
+    @mapping.setter
+    def mapping(self, val):
+        if self.PropMap is None:
+            val._assertMatchesPair(self.mapPair)
+            self._mapping = val
+        else:
+            self._mapping = self.PropMap(val)
+
     def __init__(self, mesh, mapping=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
         self.mesh = mesh
         self.mapping = mapping or Maps.IdentityMap(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
 
     @property
     def survey(self):
@@ -321,6 +89,7 @@ class BaseProblem(object):
         self.survey._prob = None
         self._survey = None
 
+
     deleteTheseOnModelUpdate = [] # List of strings, e.g. ['_MeSigma', '_MeSigmaI']
 
     @property
@@ -333,7 +102,10 @@ class BaseProblem(object):
     def curModel(self, value):
         if value is self.curModel:
             return # it is the same!
-        self._curModel = Models.Model(value, self.mapping)
+        if self.PropMap is not None:
+            self._curModel = self.mapping(value)
+        else:
+            self._curModel = Models.Model(value, self.mapping)
         for prop in self.deleteTheseOnModelUpdate:
             if hasattr(self, prop):
                 delattr(self, prop)
@@ -437,19 +209,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/PropMaps.py

@@ -0,0 +1,337 @@
+import Utils, Maps, numpy as np, scipy.sparse as sp
+
+class Property(object):
+
+    name           = ''
+    doc            = ''
+
+    defaultVal     = None
+    defaultInvProp = False
+
+    def __init__(self, doc, **kwargs):
+        # Set the default after all other params are set
+        self.doc = doc
+        Utils.setKwargs(self, **kwargs)
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    @property
+    def propertyLink(self):
+        "Can be something like: ('sigma', Maps.ReciprocalMap)"
+        return getattr(self, '_propertyLink', None)
+    @propertyLink.setter
+    def propertyLink(self, value):
+        assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name
+        self._propertyLink = value
+
+    def _getMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sMap'%prop.name, None)
+        def fset(self, val):
+            if prop.propertyLink is not None:
+                linkName, linkMap = prop.propertyLink
+                assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.'
+            # TODO: Check if the mapping can be correct
+            setattr(self, '_%sMap'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getIndexProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sIndex'%prop.name, slice(None))
+        def fset(self, val):
+            setattr(self, '_%sIndex'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return prop.defaultVal
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkMap = linkMapClass(None)
+                if getattr(self, '%sMap'%linkName, None) is None:
+                    return prop.defaultVal
+                m = getattr(self, '%s'%linkName)
+                return linkMap * m
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping * m
+        return property(fget=fget)
+
+    def _getModelDerivProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return None
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkedMap = getattr(self, '%sMap'%linkName)
+                if linkedMap is None:
+                    return None
+                linkMap = linkMapClass(None) * linkedMap
+                m = getattr(self, '%s'%linkName)
+                return linkMap.deriv( m )
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping.deriv( m )
+        return property(fget=fget)
+
+    def _getModelProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            index = getattr(self.propMap, '%sIndex'%prop.name)
+            return self.vector[index]
+        return property(fget=fget)
+
+    def _getModelProjProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            inds = getattr(self.propMap, '%sIndex'%prop.name)
+            if type(inds) is slice:
+                inds = range(*inds.indices(self.nP))
+            nI, nP = len(inds),self.nP
+            return sp.csr_matrix((np.ones(nI), (range(nI), inds) ), shape=(nI, nP))
+        return property(fget=fget)
+
+    def _getModelMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self.propMap, '_%sMap'%prop.name, None)
+        return property(fget=fget)
+
+
+
+class PropModel(object):
+    def __init__(self, propMap, vector):
+        self.propMap = propMap
+        self.vector  = vector
+        assert len(self.vector) == self.nP
+
+    # Pickleing support methods
+    # def __reduce__(self):
+    #     return (dict,{self.propMap,self.vector})
+
+    # def __getstate__(self):
+    #     '''
+    #     Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+    #     Used when doing:
+    #         pickle.dump(pickleFile,object)
+    #     '''
+    #     self.__class__ = ProbModel
+    #     odict = {}
+    #     odict['vec'] = self.__dict__['vector']
+    #     odict['pMap'] = self.__dict__['propMap']
+    #     # Return the dict
+    #     return odict
+
+    # def __setstate__(self,odict):
+    #     '''
+    #     Function that sets a pickle dictionary in to an object.
+
+    #     Used when doing:
+    #         object = pickle.load(pickleFile)
+    #     '''
+    #     # Update the dict
+    #     # Re-hook the methods to the object
+    #     self.propMap = odict['prMap']
+    #     self.vector = odict['vec']
+
+    @property
+    def nP(self):
+        inds = []
+        if getattr(self, '_nP', None) is None:
+            for name in self.propMap._properties:
+                index = getattr(self.propMap, '%sIndex'%name, None)
+                if index is not None:
+                    if type(index) is slice:
+                        inds += range(*index.indices(len(self.vector)))
+                    else:
+                        inds += list(index)
+            self._nP = len(set(inds))
+        return self._nP
+
+    def __contains__(self, val):
+        return val in self.propMap
+
+
+
+_PROPMAPCLASSREGISTRY = {}
+
+class _PropMapMetaClass(type):
+    def __new__(cls, name, bases, attrs):
+        assert name.endswith('PropMap'), 'Please use convention: ___PropMap, e.g. ElectromagneticPropMap'
+        _properties = {}
+        for base in bases:
+            for baseProp in getattr(base, '_properties', {}):
+                _properties[baseProp] = base._properties[baseProp]
+        keys = [key for key in attrs]
+        for attr in keys:
+            if isinstance(attrs[attr], Property):
+                attrs[attr].name = attr
+                attrs[attr + 'Map'  ] = attrs[attr]._getMapProperty()
+                attrs[attr + 'Index'] = attrs[attr]._getIndexProperty()
+                _properties[attr] = attrs[attr]
+                attrs.pop(attr)
+
+        attrs['_properties'] = _properties
+
+        defaultInvProps = []
+        for p in _properties:
+            prop = _properties[p]
+            if prop.defaultInvProp:
+                defaultInvProps += [p]
+            if prop.propertyLink is not None:
+                assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0])
+        if len(defaultInvProps) > 1:
+            raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
+
+        newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
+
+        newClass.PropModel = cls.createPropModelClass(newClass, name, _properties)
+
+        _PROPMAPCLASSREGISTRY[name] = newClass
+        return newClass
+
+    def createPropModelClass(self, name, _properties):
+
+        attrs = dict()
+
+        for attr in _properties:
+            prop = _properties[attr]
+
+            attrs[attr          ] = prop._getProperty()
+            attrs[attr + 'Map'  ] = prop._getModelMapProperty()
+            attrs[attr + 'Proj' ] = prop._getModelProjProperty()
+            attrs[attr + 'Model'] = prop._getModelProperty()
+            attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
+
+        return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
+
+
+class PropMap(object):
+    __metaclass__ = _PropMapMetaClass
+
+    def __init__(self, mappings):
+        """
+            PropMap takes a multi parameter model and maps it to the equivalent PropModel
+        """
+        if type(mappings) is dict:
+            assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"'
+            self.setup(mappings['maps'], slices=mappings['slices'])
+        elif type(mappings) is list:
+            self.setup(mappings)
+        elif isinstance(mappings, Maps.IdentityMap):
+            self.setup([(self.defaultInvProp, mappings)])
+        else:
+            raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
+    def setup(self, maps, slices=None):
+        """
+            Sets up the maps and slices for the PropertyMap
+
+
+            :param list maps: [('sigma', sigmaMap), ('mu', muMap), ...]
+            :param list slices: [('sigma', slice(0,nP)), ('mu', [1,2,5,6]), ...]
+
+        """
+        assert np.all([
+                type(m) is tuple and
+                len(m)==2 and
+                type(m[0]) is str and
+                m[0] in self._properties and
+                isinstance(m[1], Maps.IdentityMap)
+                for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties]))
+        if slices is None:
+            slices = dict()
+        else:
+            assert np.all([
+                s in self._properties and
+                (type(slices[s]) in [slice, list] or isinstance(slices[s], np.ndarray))
+                for s in slices]), 'Slices must be for each property'
+
+        self.clearMaps()
+
+        nP = 0
+        for name, mapping in maps:
+            setattr(self, '%sMap'%name, mapping)
+            setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
+            nP += mapping.nP
+        self.nP = nP
+
+    @property
+    def defaultInvProp(self):
+        for name in self._properties:
+            p = self._properties[name]
+            if p.defaultInvProp:
+                return p.name
+
+    def clearMaps(self):
+        for name in self._properties:
+            setattr(self, '%sMap'%name, None)
+            setattr(self, '%sIndex'%name, None)
+
+    def __call__(self, vec):
+        return self.PropModel(self, vec)
+
+    def __contains__(self, val):
+        activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None]
+        return val in activeMaps

SimPEG/Regularization.py

@@ -27,6 +27,34 @@ class BaseRegularization(object):
         self.mapping = mapping or Maps.IdentityMap(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -40,7 +68,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 +109,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 +224,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 +240,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
@@ -261,14 +289,64 @@ class Tikhonov(BaseRegularization):
         return self._Wzz
 
     @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wx, self.Wxx)
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx, self.Wxx)
             if self.mesh.dim > 1:
                 wlist += (self.Wy, self.Wyy)
             if self.mesh.dim > 2:
                 wlist += (self.Wz, self.Wzz)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wsmooth)
             self._W = sp.vstack(wlist)
         return self._W
 
+    @Utils.timeIt
+    def eval(self, m):
+        if self.smoothModel == True:
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - 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(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.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

@@ -1,4 +1,4 @@
-import Utils, numpy as np, scipy.sparse as sp
+import Utils, numpy as np, scipy.sparse as sp, uuid
 
 
 class BaseRx(object):
@@ -6,18 +6,48 @@ 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
 
     storeProjections = True #: Store calls to getP (organized by mesh)
 
     def __init__(self, locs, rxType, **kwargs):
+        self.uid = str(uuid.uuid4())
         self.locs = locs
         self.rxType = rxType
         self._Ps = {}
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+
     @property
     def rxType(self):
         """Receiver Type"""
@@ -111,37 +141,50 @@ 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.uid = str(uuid.uuid4())
         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
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
 
     @property
     def nD(self):
@@ -155,59 +198,79 @@ 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.uid = str(uuid.uuid4())
         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)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
     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
 
 
@@ -222,23 +285,64 @@ class BaseSurvey(object):
     mtrue = None     #: True model, if data is synthetic
 
     counter = None   #: A SimPEG.Utils.Counter object
+    srcPair = BaseSrc  #: Source Pair
 
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    txPair = BaseTx  #: Transmitter Pair
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
 
     @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
+        self._sourceOrder = dict()
+        [self._sourceOrder.setdefault(src.uid, ii) for ii, src in enumerate(self._srcList)]
+
+    def getSourceIndex(self, sources):
+        if type(sources) is not list:
+            sources = [sources]
+        for src in sources:
+            if getattr(src,'uid',None) is None:
+                raise KeyError('Source does not have a uid: %s'%str(src))
+        inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
+        if None in inds:
+            raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds))
+        return inds
 
     @property
     def prob(self):
@@ -282,12 +386,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,19 +459,21 @@ 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)
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
+        return self.dobs

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_Fields.py

@@ -1,71 +1,53 @@
 import unittest
 from SimPEG import *
 
-class DataAndFieldsTest(unittest.TestCase):
+
+class FieldsTest(unittest.TestCase):
 
     def setUp(self):
         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
 
+    def test_contains(self):
+        F = self.F
+        nSrc = F.survey.nSrc
+        self.assertTrue('b' not in F)
+        self.assertTrue('e' not in F)
+        e = np.random.rand(F.mesh.nE, nSrc)
+        F[:, 'e'] = e
+        self.assertTrue('b' not in F)
+        self.assertTrue('e' in F)
+
     def test_overlappingFields(self):
         self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
                             knownFields={'b':'F'},
                             aliasFields={'b':['b',(lambda F, b, ind: b)]})
 
-    def test_data(self):
-        V = []
-        for tx in self.D.survey.txList:
-            for rx in tx.rxList:
-                v = np.random.rand(rx.nD)
-                V += [v]
-                self.D[tx, rx] = v
-                self.assertTrue(np.all(v == self.D[tx, rx]))
-        V = np.concatenate(V)
-        self.assertTrue(np.all(V == Utils.mkvc(self.D)))
-
-        D2 = Survey.Data(self.D.survey, V)
-        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
-
-    def test_contains(self):
-        F = self.F
-        nTx = F.survey.nTx
-        self.assertTrue('b' not in F)
-        self.assertTrue('e' not in F)
-        e = np.random.rand(F.mesh.nE, nTx)
-        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_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 +61,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 +102,68 @@ 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)
-        self.D = Survey.Data(survey)
+        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.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[0] 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 +172,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 +207,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 +225,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 +274,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 +310,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 +322,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[0] 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_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_PropMaps.py

@@ -0,0 +1,193 @@
+import unittest
+from SimPEG import *
+from scipy.constants import mu_0
+
+
+class MyPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True)
+    mu    = Maps.Property("Mu", defaultVal=mu_0)
+
+class MyReciprocalPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True, propertyLink=('rho',   Maps.ReciprocalMap))
+    rho   = Maps.Property("Electrical Resistivity",                       propertyLink=('sigma', Maps.ReciprocalMap))
+    mu    = Maps.Property("Mu", defaultVal=mu_0, propertyLink=('mui',   Maps.ReciprocalMap))
+    mui    = Maps.Property("Mu", defaultVal=1./mu_0, propertyLink=('mu',   Maps.ReciprocalMap))
+
+
+class TestPropMaps(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def test_setup(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        assert expMap.nP == 3
+
+        PM1 = MyPropMap(expMap)
+        PM2 = MyPropMap([('sigma', expMap)])
+        PM3 = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':slice(0,3)}})
+
+        for PM in [PM1,PM2,PM3]:
+            assert PM.defaultInvProp == 'sigma'
+            assert PM.sigmaMap is not None
+            assert PM.sigmaMap is expMap
+            assert PM.sigmaIndex == slice(0,3)
+            assert getattr(PM, 'sigma', None) is None
+            assert PM.muMap is None
+            assert PM.muIndex is None
+
+            assert 'sigma' in PM
+            assert 'mu' not in PM
+            assert 'mui' not in PM
+
+            m = PM(np.r_[1.,2,3])
+
+            assert 'sigma' in m
+            assert 'mu' not in m
+            assert 'mui' not in m
+
+            assert m.mu == mu_0
+            assert m.muModel is None
+            assert m.muMap is None
+            assert m.muDeriv is None
+
+            assert np.all(m.sigmaModel == np.r_[1.,2,3])
+            assert m.sigmaMap is expMap
+            assert np.all(m.sigma == np.exp(np.r_[1.,2,3]))
+            assert m.sigmaDeriv is not None
+
+            assert m.nP == 3
+
+    def test_slices(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}})
+        assert PM.sigmaIndex == [2,1,0]
+        m = PM(np.r_[1.,2,3])
+        assert np.all(m.sigmaModel == np.r_[3,2,1])
+        assert np.all(m.sigma == np.exp(np.r_[3,2,1]))
+
+    def test_multiMap(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap([('sigma', expMap), ('mu', iMap)])
+
+        pm = PM(np.r_[1.,2,3,4,5,6])
+
+        assert pm.nP == 6
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1.,2,3])
+        assert np.all(pm.sigma == np.exp([1.,2,3]))
+        assert np.all(pm.muModel == [4.,5,6])
+        assert np.all(pm.mu == [4.,5,6])
+
+
+    def test_multiMapCompressed(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap({'maps':[('sigma', expMap), ('mu', iMap)],'slices':{'mu':[0,1,2]}})
+
+        pm = PM(np.r_[1,2.,3])
+
+        assert pm.nP == 3
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1,2,3])
+        assert np.all(pm.sigma == np.exp([1,2,3]))
+        assert np.all(pm.muModel == [1,2,3])
+        assert np.all(pm.mu == [1,2,3])
+
+    def test_Projections(self):
+        m = Mesh.TensorMesh((3,))
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        v = np.r_[1,2.,3]
+        pm = PM(v)
+
+        assert pm.sigmaProj is not None
+        assert pm.rhoProj   is None
+        assert pm.muProj    is None
+        assert pm.muiProj   is None
+
+        assert np.all(pm.sigmaProj * v == pm.sigmaModel)
+
+    def test_Links(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == [1,2,3])
+        assert np.all(pm.rho == 1./np.r_[1,2,3])
+        assert pm.sigmaMap is iMap
+        assert pm.rhoMap is None
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' in PM
+        assert 'rho' not in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' in pm
+        assert 'rho' not in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        assert pm.mu == mu_0
+        assert pm.mui == 1.0/mu_0
+        assert pm.muMap is None
+        assert pm.muDeriv is None
+        assert pm.muiMap is None
+        assert pm.muiDeriv is None
+
+        PM = MyReciprocalPropMap([('rho', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == 1./np.r_[1,2,3])
+        assert np.all(pm.rho == [1,2,3])
+        assert pm.sigmaMap is None
+        assert pm.rhoMap is iMap
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' not in PM
+        assert 'rho' in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' not in pm
+        assert 'rho' in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)])
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)])
+
+        MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine
+
+if __name__ == '__main__':
+    unittest.main()
+
+

SimPEG/Tests/test_SurveyAndData.py

@@ -0,0 +1,53 @@
+import unittest
+from SimPEG import *
+
+class TestData(unittest.TestCase):
+
+    def setUp(self):
+        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.])
+        srcLoc = np.r_[0,0,0.]
+        rxList0 = Survey.BaseRx(XYZ, 'exi')
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
+        rxList1 = Survey.BaseRx(XYZ, 'bxi')
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
+        rxList2 = Survey.BaseRx(XYZ, 'bxi')
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
+        rxList3 = Survey.BaseRx(XYZ, 'bxi')
+        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)
+
+    def test_data(self):
+        V = []
+        for src in self.D.survey.srcList:
+            for rx in src.rxList:
+                v = np.random.rand(rx.nD)
+                V += [v]
+                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)))
+
+        D2 = Survey.Data(self.D.survey, V)
+        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
+
+    def test_uniqueSrcs(self):
+        srcs = self.D.survey.srcList
+        srcs += [srcs[0]]
+        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
+
+    def test_sourceIndex(self):
+        survey = self.D.survey
+        srcs = survey.srcList
+        assert survey.getSourceIndex([srcs[1],srcs[0]]) == [1,0]
+        assert survey.getSourceIndex([srcs[1],srcs[2],srcs[2]]) == [1,2,2]
+        SrcNotThere = Survey.BaseSrc(srcs[0].rxList, loc=np.r_[0,0,0])
+        self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere])
+        self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere])
+
+if __name__ == '__main__':
+    unittest.main()

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", "FullMap"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
+
 class MapTests(unittest.TestCase):
 
     def setUp(self):
@@ -13,17 +16,36 @@ 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_reciprocalMap(self):
+        # Note that log/reciprocal 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))
+
+        maps = Maps.ReciprocalMap(self.mesh2)
+        self.assertTrue(maps.test(v2, dx=dv2))
+        maps = Maps.ReciprocalMap(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 +112,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_regularization.py

@@ -22,7 +22,7 @@ class RegularizationTests(unittest.TestCase):
             mapping = r.mapPair(self.mesh2)
             reg = r(self.mesh2, mapping=mapping)
             m = np.random.rand(mapping.nP)
-            reg.mref = m[:]*0
+            reg.mref = m[:]*np.mean(m)
             passed = checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
             self.assertTrue(passed)
 

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/ModelBuilder.py

@@ -3,10 +3,27 @@ import scipy.ndimage as ndi
 import scipy.sparse as sp
 from matutils import mkvc
 
-
-def getIndecesBlock(p0,p1,ccMesh):
+def addBlock(gridCC, modelCC, p0, p1, blockProp):
     """
-        Creates a vector containing the block indexes in the cell centerd mesh.
+        Add a block to an exsisting cell centered model, modelCC
+
+        :param numpy.array, gridCC: mesh.gridCC is the cell centered grid
+        :param numpy.array, modelCC: cell centered model
+        :param numpy.array, p0: bottom, southwest corner of block
+        :param numpy.array, p1: top, northeast corner of block
+        :blockProp float, blockProp: property to assign to the model 
+
+        :return numpy.array, modelBlock: model with block
+    """
+    ind = getIndicesBlock(p0, p1, gridCC)
+    modelBlock = modelCC.copy()
+    modelBlock[ind] = blockProp
+    return modelBlock
+
+
+def getIndicesBlock(p0,p1,ccMesh):
+    """
+        Creates a vector containing the block indices in the cell centers mesh.
         Returns a tuple
 
         The block is defined by the points
@@ -78,7 +95,7 @@ def defineBlock(ccMesh,p0,p1,vals=[0,1]):
         vals[1]  conductivity of the ground
     """
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0]
 
@@ -132,7 +149,7 @@ def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     # The depth is always defined on the last one.
     p1[len(p1)-1] -= depth
 
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0];
 
@@ -153,16 +170,58 @@ def scalarConductivity(ccMesh,pFunction):
 
     return mkvc(sigma)
 
+def layeredModel(ccMesh, layerTops, layerValues):
+    """
+        Define a layered model from layerTops (z-positive up)
+
+        :param numpy.array ccMesh: cell-centered mesh
+        :param numpy.array layerTops: z-locations of the tops of each layer
+        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top) 
+        :rtype: numpy.array
+        :return: M, layered model on the mesh 
+    """
+
+    descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
+
+    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts 
+    # assert ascending or descending, "Layers must be listed in either ascending or descending order"
+
+    # start from bottom up
+    if not descending:
+        zprop = np.hstack([mkvc(layerTops,2),mkvc(layerValues,2)])
+        zprop.sort(axis=0)
+        layerTops, layerValues  = zprop[::-1,0], zprop[::-1,1]
+
+    # put in vector form
+    layerTops, layerValues = mkvc(layerTops), mkvc(layerValues)
+
+    # initialize with bottom layer
+    dim = ccMesh.shape[1]
+    if dim == 3:
+        z = ccMesh[:,2]
+    elif dim == 2:
+        z = ccMesh[:,1]
+    elif dim == 1:
+        z = ccMesh[:,0]
+
+    model = np.zeros(ccMesh.shape[0])
+
+    for i, top in enumerate(layerTops):
+        zind = z <= top 
+        model[zind] = layerValues[i]
+
+    return model 
+
 
 
 def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
-        Create a random model by convolving a kernal with a
+        Create a random model by convolving a kernel with a
         uniformly distributed model.
 
         :param int,tuple shape: shape of the model.
         :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
+        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
         :param int its: number of smoothing iterations
         :param list bounds: bounds on the model, len(list) == 2
         :rtype: numpy.ndarray

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/codeutils.py

@@ -58,9 +58,11 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
         print 'Method '+name+' was not overwritten.'
 
 
-def setKwargs(obj, **kwargs):
+def setKwargs(obj, ignore=[],  **kwargs):
     """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
     for attr in kwargs:
+        if attr in ignore:
+            continue
         if hasattr(obj, attr):
             setattr(obj, attr, kwargs[attr])
         else:

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
@@ -198,6 +205,200 @@ 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
+    if model is not None:
+        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])
+
+
+    # Check the extension of the fileName
+    if fileName is not None:
+        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.SetInputData(vtkObj)
+        vtrWriteFilter.SetFileName(fileName)
+        vtrWriteFilter.Update()
+    else:
+        return vtkObj
+
+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

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"],