copy docs badges from readme to docs index

Feat/codecov

.bumpversion.cfg

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

.gitignore

@@ -39,3 +39,5 @@ nosetests.xml
 *.sublime-workspace
 docs/_build/
 Makefile
+docs/warnings.txt
+.DS_Store

.travis.yml

@@ -24,18 +24,25 @@ env:
     - TEST_DIR=tests/examples
     - TEST_DIR=tests/em/fdem/inverse/adjoint
     - TEST_DIR=tests/em/fdem/forward
+    - TEST_DIR=tests/docs;
+      GAE_PYTHONPATH=${HOME}/.cache/google_appengine;
+      PATH=$PATH:${HOME}/google-cloud-sdk/bin;
+      PYTHONPATH=${PYTHONPATH}:${GAE_PYTHONPATH};
+      CLOUDSDK_CORE_DISABLE_PROMPTS=1
 
 # Setup anaconda
 before_install:
-  - 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
+# Install packages
+  - 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
   - conda update --yes conda
 
-# Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git
@@ -46,11 +53,26 @@ install:
 
 # Run test
 script:
+  # test docs
   - nosetests $TEST_DIR --with-cov --cov SimPEG --cov-config .coveragerc -v -s
 
 # Calculate coverage
 after_success:
-  - coveralls --config_file .coveragerc
+  - bash <(curl -s https://codecov.io/bash)
+  - if [ "$TRAVIS_BRANCH" = "master" -a "$TRAVIS_PULL_REQUEST" = "false" ]; then
+      if [ ${TEST_DIR} == "tests/docs" ]; then
+        python scripts/fetch_gae_sdk.py $(dirname "${GAE_PYTHONPATH}");
+        openssl aes-256-cbc -K $encrypted_93066031461c_key -iv $encrypted_93066031461c_iv
+          -in docs/credentials.tar.gz.enc -out credentials.tar.gz -d ;
+        if [ ! -d ${HOME}/google-cloud-sdk ]; then curl https://sdk.cloud.google.com | bash; fi ;
+        tar -xzf credentials.tar.gz ;
+        gcloud auth activate-service-account --key-file client-secret.json ;
+        gcloud config set project simpegdocs;
+        gcloud -q components update gae-python;
+        gcloud -q preview app deploy ./docs/app.yaml --version ${TRAVIS_COMMIT} --promote;
+      fi;
+    fi
+
 
 notifications:
   email:

README.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
     :alt: SimPEG Logo
 
 ======
@@ -21,14 +21,21 @@ SimPEG
     :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
-
 .. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
     :alt: gitter chat room at https://gitter.im/simpeg/simpeg
     :target: https://gitter.im/simpeg/simpeg
-    
+
+.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
+    :target: https://codecov.io/gh/simpeg/simpeg
+
+.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
+    :target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
+    :alt: Code issues
+
+.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
+    :target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&utm_medium=referral&utm_content=simpeg/simpeg&utm_campaign=Badge_Grade
+    :alt: codacy
+
 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:

SimPEG/DCIP/BaseDC.py

@@ -162,8 +162,8 @@ class ProblemDC_CC(Problem.BaseProblem):
         """
             Makes the matrix A(m) for the DC resistivity problem.
 
-            :param numpy.array m: model
+            :param numpy.ndarray m: model
-            :rtype: scipy.csc_matrix
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/BaseIP.py

@@ -71,7 +71,7 @@ class ProblemIP(Problem.BaseProblem):
             Makes the matrix A(m) for the DC resistivity problem.
 
             :param numpy.array m: model
-            :rtype: scipy.csc_matrix
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/DCIPUtils.py

@@ -476,7 +476,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
     fid.write('! ' + surveyType + ' FORMAT\n')
 
     if iptype!=0:
-        fid.write('IPTYPE=%i\n'%iptype)
+        fid.write('IPTYPE={0:d}\n'.format(iptype))
 
     else:
         fid.write('! ' + stype + ' FORMAT\n')
@@ -512,7 +512,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
 
                 if surveyType == 'SURFACE':
 
-                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("{0:f} ".format(ii) for ii in mkvc(tx[0,:]))
                     M = M[:,0]
                     N = N[:,0]
 
@@ -521,7 +521,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
                     # Flip sign for z-elevation to depth
                     tx[2::2,:] = -tx[2::2,:]
 
-                    fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
+                    fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[::2,:]))
                     M = M[:,0::2]
                     N = N[:,0::2]
 
@@ -529,22 +529,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
                     M[:,1::2] = -M[:,1::2]
                     N[:,1::2] = -N[:,1::2]
 
-                fid.write('%i\n'% nD)
+                fid.write('{0:d}\n'.format(nD))
                 np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
 
         if dim=='3D':
 
             if surveyType == 'SURFACE':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
+                fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:2,:]))
                 M = M[:,0:2]
                 N = N[:,0:2]
 
             if surveyType == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
+                fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:3,:]))
 
-            fid.write('%i\n'% nD)
+            fid.write('{0:d}\n'.format(nD))
             np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
             fid.write('\n')
 

SimPEG/Directives.py

@@ -15,7 +15,7 @@ class InversionDirective(object):
     @inversion.setter
     def inversion(self, i):
         if getattr(self,'_inversion',None) is not None:
-            print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__
+            print 'Warning: InversionDirective {0!s} has switched to a new inversion.'.format(self.__name__)
         self._inversion = i
 
     @property
@@ -47,7 +47,7 @@ class DirectiveList(object):
     def __init__(self, *directives, **kwargs):
         self.dList = []
         for d in directives:
-            assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not %s' % d.__name__
+            assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not {0!s}'.format(d.__name__)
             self.dList.append(d)
         Utils.setKwargs(self, **kwargs)
 
@@ -68,7 +68,7 @@ class DirectiveList(object):
     def inversion(self, i):
         if self.inversion is i: return
         if getattr(self,'_inversion',None) is not None:
-            print 'Warning: %s has switched to a new inversion.' % self.__name__
+            print 'Warning: {0!s} has switched to a new inversion.'.format(self.__name__)
         for d in self.dList:
             d.inversion = i
         self._inversion = i
@@ -79,7 +79,7 @@ class DirectiveList(object):
             return
 
         directives = ['initialize', 'endIter', 'finish']
-        assert ruleType in directives, 'Directive type must be in ["%s"]' % '", "'.join(directives)
+        assert ruleType in directives, 'Directive type must be in ["{0!s}"]'.format('", "'.join(directives))
         for r in self.dList:
             getattr(r, ruleType)()
 
@@ -141,7 +141,7 @@ class BetaSchedule(InversionDirective):
 
     def endIter(self):
         if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
-            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
             self.invProb.beta /= self.coolingFactor
 
 
@@ -167,7 +167,7 @@ class TargetMisfit(InversionDirective):
 
 
 
-class _SaveEveryIteration(InversionDirective):
+class SaveEveryIteration(InversionDirective):
     @property
     def name(self):
         if getattr(self, '_name', None) is None:
@@ -181,42 +181,42 @@ class _SaveEveryIteration(InversionDirective):
     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'))
+            self._fileName = '{0!s}-{1!s}'.format(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):
+class SaveModelEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
+        print "SimPEG.SaveModelEveryIteration will save your models as: '###-{0!s}.npy'".format(self.fileName)
 
     def endIter(self):
-        np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
+        np.save('{0:03d}-{1!s}'.format(self.opt.iter, self.fileName), self.opt.xc)
 
 
-class SaveOutputEveryIteration(_SaveEveryIteration):
+class SaveOutputEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
+        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-{0!s}.txt'".format(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.write(' {0:3d} {1:1.4e} {2:1.4e} {3:1.4e} {4:1.4e}\n'.format(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
         f.close()
 
-class SaveOutputDictEveryIteration(_SaveEveryIteration):
+class SaveOutputDictEveryIteration(SaveEveryIteration):
     """SaveOutputDictEveryIteration"""
 
     def initialize(self):
-        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-{0!s}.npz'".format(self.fileName)
 
     def endIter(self):
         # Save the data.
@@ -328,7 +328,7 @@ class Update_IRLS(InversionDirective):
 
         # Beta Schedule
         if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
-            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
             self.invProb.beta /= self.coolingFactor
 
 
@@ -340,11 +340,11 @@ class Update_IRLS(InversionDirective):
             phim_new = self.reg.eval(self.invProb.curModel)
             self.f_change = np.abs(self.f_old - phim_new) / self.f_old
 
-            print "Regularization decrease: %6.3e" % (self.f_change)
+            print "Regularization decrease: {0:6.3e}".format((self.f_change))
 
             # Check for maximum number of IRLS cycles
             if self.IRLSiter == self.maxIRLSiter:
-                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
+                print "Reach maximum number of IRLS cycles: {0:d}".format(self.maxIRLSiter)
                 self.opt.stopNextIteration = True
                 return
 

SimPEG/EM/Analytics/FDEMDipolarfields.py

@@ -0,0 +1,302 @@
+from __future__ import division
+import numpy as np
+from scipy.constants import mu_0, pi, epsilon_0
+from scipy.special import erf
+from SimPEG import Utils
+
+omega = lambda f: 2.*np.pi*f
+# TODO:
+# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.)
+# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.):
+
+    """
+        Computing Analytic Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ey = front*(dx*dy  / r**2)*mid
+        Ez = front*(dx*dz  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ez = front*(dy*dz  / r**2)*mid
+        Ex = front*(dy*dx  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ex = front*(dz*dx  / r**2)*mid
+        Ey = front*(dz*dy  / r**2)*mid
+        return Ex, Ey, Ez
+
+
+def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.))
+        Ey_galvanic = front*(dx*dy  / r**2)*mid
+        Ez_galvanic = front*(dx*dz  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.))
+        Ez_galvanic = front*(dy*dz  / r**2)*mid
+        Ex_galvanic = front*(dy*dx  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.))
+        Ex_galvanic = front*(dz*dx  / r**2)*mid
+        Ey_galvanic = front*(dz*dy  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+
+def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Ex_inductive = front*(k**2 * r**2)
+        Ey_inductive = np.zeros_like(Ex_inductive)
+        Ez_inductive = np.zeros_like(Ex_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_inductive = front*(k**2 * r**2)
+        Ez_inductive = np.zeros_like(Ey_inductive)
+        Ex_inductive = np.zeros_like(Ey_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_inductive = front*(k**2 * r**2)
+        Ex_inductive = np.zeros_like(Ez_inductive)
+        Ey_inductive = np.zeros_like(Ez_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+
+def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx = sig*Ex
+    Jy = sig*Ey
+    Jz = sig*Ez
+    return Jx, Jy, Jz
+
+
+def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx_galvanic = sig*Ex_galvanic
+    Jy_galvanic = sig*Ey_galvanic
+    Jz_galvanic = sig*Ez_galvanic
+    return Jx_galvanic, Jy_galvanic, Jz_galvanic
+
+
+def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Jx_inductive = sig*Ex_inductive
+    Jy_inductive = sig*Ey_inductive
+    Jz_inductive = sig*Ez_inductive
+    return Jx_inductive, Jy_inductive, Jz_inductive
+
+
+def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Hy = front*(-dz  / r)
+        Hz = front*(dy  / r)
+        Hx = np.zeros_like(Hy)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Y':
+        Hx = front*(dz  / r)
+        Hz = front*(-dx  / r)
+        Hy = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Z':
+        Hx = front*(-dy / r)
+        Hy = front*(dx  / r)
+        Hz = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+
+def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic flux densites from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Bx = mu*Hx
+    By = mu*Hy
+    Bz = mu*Hz
+    return Bx, By, Bz
+
+
+def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Electric vector potentials from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*r)
+
+    if orientation.upper() == 'X':
+        Ax = front*np.exp(-1j*k*r)
+        Ay = np.zeros_like(Ax)
+        Az = np.zeros_like(Ax)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Y':
+        Ay = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Az = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Z':
+        Az = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Ay = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+
+
+
+

SimPEG/EM/Analytics/__init__.py

@@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
 from DC import DCAnalyticHalf, DCAnalyticSphere
+from FDEMDipolarfields import *

SimPEG/EM/Base.py

@@ -20,10 +20,10 @@ class BaseEMProblem(Problem.BaseProblem):
         Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
 
-    surveyPair = Survey.BaseSurvey
+    surveyPair = Survey.BaseSurvey #: The survey to pair with.
-    dataPair = Survey.Data
+    dataPair = Survey.Data #: The data to pair with.
 
-    PropMap = EMPropMap
+    PropMap = EMPropMap #: The property mapping
 
     Solver = SimpegSolver
     solverOpts = {}
@@ -217,6 +217,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
     def eval(self, f):
         """
         Project fields to receiver locations
+
         :param Fields u: fields object
         :rtype: numpy.ndarray
         :return: data

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -6,11 +6,11 @@ from SimPEG.EM.Utils import omega
 from SimPEG.Utils import Zero, Identity, sdiag
 
 
-class Fields(SimPEG.Problem.Fields):
+class FieldsFDEM(SimPEG.Problem.Fields):
     """
 
     Fancy Field Storage for a FDEM survey. Only one field type is stored for
-    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+    each problem, the rest are computed. The fields object acts like an array and is indexed by
 
     .. code-block:: python
 
@@ -42,7 +42,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total electric field
         """
         if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
-            raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting e from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
 
@@ -56,7 +56,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total magnetic flux density
         """
         if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
-            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting b from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
@@ -70,7 +70,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total magnetic field
         """
         if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
-            raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting h from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
 
@@ -84,7 +84,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: total current density
         """
         if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
-            raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting j from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
 
@@ -92,7 +92,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -100,7 +100,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -110,7 +110,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -118,7 +118,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
-            raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting bDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -128,7 +128,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -136,7 +136,7 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
-            raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting hDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -146,7 +146,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -154,18 +154,18 @@ class Fields(SimPEG.Problem.Fields):
         :return: derivative times a vector (or tuple for adjoint)
         """
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields3D_e(Fields):
+class Fields3D_e(FieldsFDEM):
     """
     Fields object for Problem3D_e.
 
-    :param Mesh mesh: mesh
+    :param BaseMesh mesh: mesh
-    :param Survey survey: survey
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -180,9 +180,6 @@ class Fields3D_e(Fields):
                     'h' : ['eSolution','CCV','_h'],
                   }
 
-    def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -426,12 +423,12 @@ class Fields3D_e(Fields):
 
 
 
-class Fields3D_b(Fields):
+class Fields3D_b(FieldsFDEM):
     """
     Fields object for Problem3D_b.
 
-    :param Mesh mesh: mesh
+    :param BaseMesh mesh: mesh
-    :param Survey survey: survey
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -446,9 +443,6 @@ class Fields3D_b(Fields):
                     'h' : ['bSolution','CCV','_h'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -693,12 +687,12 @@ class Fields3D_b(Fields):
         return Zero()
 
 
-class Fields3D_j(Fields):
+class Fields3D_j(FieldsFDEM):
     """
     Fields object for Problem3D_j.
 
-    :param Mesh mesh: mesh
+    :param BaseMesh mesh: mesh
-    :param Survey survey: survey
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -713,9 +707,6 @@ class Fields3D_j(Fields):
                     'b' : ['jSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -988,12 +979,12 @@ class Fields3D_j(Fields):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields3D_h(Fields):
+class Fields3D_h(FieldsFDEM):
     """
     Fields object for Problem3D_h.
 
-    :param Mesh mesh: mesh
+    :param BaseMesh mesh: mesh
-    :param Survey survey: survey
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -1008,9 +999,6 @@ class Fields3D_h(Fields):
                     'b' : ['hSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -1,7 +1,7 @@
 from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyFDEM import Survey as SurveyFDEM
-from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
+from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -31,10 +31,11 @@ class BaseFDEMProblem(BaseEMProblem):
         if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+
     """
 
     surveyPair = SurveyFDEM
-    fieldsPair = Fields
+    fieldsPair = FieldsFDEM
 
     def fields(self, m):
         """
@@ -64,7 +65,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take sensitivity product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -99,7 +100,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take adjoint product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -153,8 +154,8 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :rtype: tuple
-        :return: s_m, s_e (nE or nF, nSrc)
+        :return: (s_m, s_e) (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
@@ -194,7 +195,7 @@ class Problem3D_e(BaseFDEMProblem):
 
     which we solve for :math:`\mathbf{e}`.
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'eSolution'
@@ -269,7 +270,7 @@ class Problem3D_e(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -305,7 +306,7 @@ class Problem3D_b(BaseFDEMProblem):
     .. note ::
         The inverse problem will not work with full anisotropy
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'bSolution'
@@ -400,7 +401,7 @@ class Problem3D_b(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -444,6 +445,7 @@ class Problem3D_j(BaseFDEMProblem):
 
         \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
 
+
     and solve for \\\(\\\mathbf{j}\\\) using
 
     .. math ::
@@ -453,7 +455,7 @@ class Problem3D_j(BaseFDEMProblem):
     .. note::
         This implementation does not yet work with full anisotropy!!
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'jSolution'
@@ -529,8 +531,8 @@ class Problem3D_j(BaseFDEMProblem):
             \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray (nE, nSrc)
+        :rtype: numpy.ndarray
-        :return: RHS
+        :return: RHS (nE, nSrc)
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -549,7 +551,7 @@ class Problem3D_j(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -591,7 +593,7 @@ class Problem3D_h(BaseFDEMProblem):
 
         \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'hSolution'
@@ -608,9 +610,11 @@ class Problem3D_h(BaseFDEMProblem):
         .. math::
             \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
 
+
         :param float freq: Frequency
         :rtype: scipy.sparse.csr_matrix
         :return: A
+
         """
 
         MeMu = self.MeMu
@@ -653,6 +657,7 @@ class Problem3D_h(BaseFDEMProblem):
         :param float freq: Frequency
         :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
+
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -666,7 +671,7 @@ class Problem3D_h(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/RxFDEM.py

@@ -11,8 +11,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
     """
 
     def __init__(self, locs, orientation=None, component=None):
-        assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation
+        assert(orientation in ['x','y','z']), "Orientation {0!s} not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented.".format(orientation)
-        assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
+        assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not {0!s}".format(component)
 
         self.projComp = orientation
         self.component = component
@@ -25,10 +25,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
 
     def eval(self, src, mesh, f):
         """
-        Project fields to recievers to get data.
+        Project fields to receivers to get data.
 
-        :param Source src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
-        :param Mesh mesh: mesh used
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :rtype: numpy.ndarray
         :return: fields projected to recievers
@@ -44,8 +44,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
-        :param Source src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
-        :param Mesh mesh: mesh used
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/SrcFDEM.py

@@ -23,8 +23,8 @@ class BaseSrc(Survey.BaseSrc):
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term
         """
         s_m = self.s_m(prob)
@@ -37,10 +37,10 @@ class BaseSrc(Survey.BaseSrc):
         - :code:`s_mDeriv` : derivative of the magnetic source term
         - :code:`s_eDeriv` : derivative of the electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term derivatives times a vector
         """
         if v is not None:
@@ -52,7 +52,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic flux density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic flux density
         """
@@ -64,7 +64,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -76,7 +76,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary electric field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary electric field
         """
@@ -88,7 +88,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary current density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary current density
         """
@@ -100,7 +100,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -110,7 +110,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -120,7 +120,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of magnetic source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -133,7 +133,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of electric source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -162,7 +162,7 @@ class RawVec_e(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -191,7 +191,7 @@ class RawVec_m(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -220,7 +220,7 @@ class RawVec(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -232,7 +232,7 @@ class RawVec(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -301,7 +301,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -339,7 +339,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -350,7 +350,7 @@ class MagDipole(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -364,7 +364,7 @@ class MagDipole(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -416,7 +416,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -455,7 +455,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -479,7 +479,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -530,7 +530,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -567,7 +567,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -578,7 +578,7 @@ class CircularLoop(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -591,7 +591,7 @@ class CircularLoop(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """

SimPEG/EM/Static/DC/FieldsDC.py

@@ -9,7 +9,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -18,7 +18,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _eDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -26,7 +26,7 @@ class Fields(SimPEG.Problem.Fields):
 
     def _jDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)

SimPEG/EM/Static/DC/FieldsDC_2D.py

@@ -32,7 +32,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -41,7 +41,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -49,7 +49,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
 
     def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+            raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
 
         if adjoint:
             return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)

SimPEG/EM/Static/DC/ProblemDC.py

@@ -46,7 +46,7 @@ class BaseDCProblem(BaseEMProblem):
             du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
 
             for rx in src.rxList:
-                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                 df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
                 Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
         return Utils.mkvc(Jv)
@@ -69,7 +69,7 @@ class BaseDCProblem(BaseEMProblem):
             u_src = f[src, self._solutionType]
             for rx in src.rxList:
                 PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                 df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
                 ATinvdf_duT = self.Ainv * df_duT

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -60,7 +60,7 @@ class BaseDCProblem_2D(BaseEMProblem):
                 dRHS_dm_v = self.getRHSDeriv(ky, src, v)
                 du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
                 for rx in src.rxList:
-                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                     df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
                     # Trapezoidal intergration
                     Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
@@ -101,7 +101,7 @@ class BaseDCProblem_2D(BaseEMProblem):
                     ky = self.kys[iky]
                     AT = self.getA(ky)
                     PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                     df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
 
                     ATinvdf_duT = self.Ainv[iky] * df_duT

SimPEG/EM/Static/IP/ProblemIP.py

@@ -56,7 +56,7 @@ class BaseIPProblem(BaseEMProblem):
             du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
 
             for rx in src.rxList:
-                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                 df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
                 Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
         # Conductivity (d u / d log sigma)
@@ -83,7 +83,7 @@ class BaseIPProblem(BaseEMProblem):
             u_src = f[src, self._solutionType]
             for rx in src.rxList:
                 PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                 df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
                 ATinvdf_duT = self.Ainv * df_duT
                 dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -83,7 +83,7 @@ class BaseSIPProblem(BaseEMProblem):
                 for rx in src.rxList:
                     timeindex = rx.getTimeP(self.survey.times)
                     if timeindex[tind]:
-                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                         df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
                         Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
 
@@ -122,7 +122,7 @@ class BaseSIPProblem(BaseEMProblem):
                 for rx in src.rxList:
                     timeindex = rx.getTimeP(self.survey.times)
                     if timeindex[tind]:
-                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                         df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
                         df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
                         Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
@@ -153,7 +153,7 @@ class BaseSIPProblem(BaseEMProblem):
                     timeindex = rx.getTimeP(self.survey.times)
                     if timeindex[tind]:
                         PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
-                        df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
                         df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
                         ATinvdf_duT = self.Ainv * df_duT
                         dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)

SimPEG/EM/TDEM/BaseTDEM.py

@@ -47,7 +47,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         self.waveformType = "GENERAL"
 
     def fields(self, m):
-        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nCalculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
         self.curModel = m
         # Create a fields storage object
         F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -55,7 +55,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             # Set the initial conditions
             F[src,:,0] = src.getInitialFields(self.mesh)
         F = self.forward(m, self.getRHS, F=F)
-        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nDone calculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
         return F
 
     def forward(self, m, RHS, F=None):
@@ -70,11 +70,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
                 if Ainv is not None:
                     Ainv.clean()
                 A = self.getA(tInd)
-                if self.verbose: print 'Factoring...   (dt = %e)'%dt
+                if self.verbose: print 'Factoring...   (dt = {0:e})'.format(dt)
                 Ainv = self.Solver(A, **self.solverOpts)
                 if self.verbose: print 'Done'
             rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving...   (tInd = %d)'%tInd
+            if self.verbose: print '    Solving...   (tInd = {0:d})'.format(tInd)
             sol = Ainv * rhs
             if self.verbose: print '    Done...'
             if sol.ndim == 1:
@@ -95,11 +95,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
                 if Ainv is not None:
                     Ainv.clean()
                 A = self.getA(tInd)
-                if self.verbose: print 'Factoring (Adjoint)...   (dt = %e)'%dt
+                if self.verbose: print 'Factoring (Adjoint)...   (dt = {0:e})'.format(dt)
                 Ainv = self.Solver(A, **self.solverOpts)
                 if self.verbose: print 'Done'
             rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving (Adjoint)...   (tInd = %d)'%tInd
+            if self.verbose: print '    Solving (Adjoint)...   (tInd = {0:d})'.format(tInd)
             sol = Ainv * rhs
             if self.verbose: print '    Done...'
             if sol.ndim == 1:
@@ -112,7 +112,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
+            :param FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -123,21 +123,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
 
         """
-        if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nCalculating J(v)\n{1!s}'.format('*'*50, '*'*50)
         self.curModel = m
         if f is None:
             f = self.fields(m)
         p = self.Gvec(m, v, f)
         y = self.solveAh(m, p)
         Jv = self.survey.evalDeriv(f, v=y)
-        if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nDone calculating J(v)\n{1!s}'.format('*'*50, '*'*50)
         return - mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
-            :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
+            :param numpy.ndarray v: vector (or a :class:`SimPEG.Survey.Data` object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param FieldsTDEM u: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (model object)
 
@@ -148,7 +148,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
 
         """
-        if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nCalculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
         self.curModel = m
         if f is None:
             f = self.fields(m)
@@ -159,6 +159,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         p = self.survey.evalDeriv(f, v=v, adjoint=True)
         y = self.solveAht(m, p)
         w = self.Gtvec(m, y, f)
-        if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
+        if self.verbose: print '{0!s}\nDone calculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
         return - mkvc(w)
 

SimPEG/EM/TDEM/TDEM_b.py

@@ -87,8 +87,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a model)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param FieldsTDEM u: Fields resulting from m
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply G by a vector
@@ -125,9 +125,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a fields)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param FieldsTDEM u: Fields resulting from m
-            :rtype: np.ndarray (like a model)
+            :rtype: numpy.ndarray
-            :return: p
+            :return: p (like a model)
 
             Multiply G.T by a vector
         """
@@ -153,8 +153,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAh(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
+            :param FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -200,8 +200,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAht(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
+            :param FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -270,8 +270,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
+            :param FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -315,8 +315,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhtVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
+            :param FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where

SimPEG/EM/Utils/testingUtils.py

@@ -58,7 +58,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
                 Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
 
     if verbose:
-        print '  Fetching %s problem' % (fdemType)
+        print '  Fetching {0!s} problem'.format((fdemType))
 
     if fdemType == 'e':
         survey = EM.FDEM.Survey(Src)
@@ -94,7 +94,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
 
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
-    print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
+    print 'Cross Checking Forward: {0!s}, {1!s} formulations - {2!s}'.format(fdemType1, fdemType2, comp)
 
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -1,7 +1,7 @@
 from SimPEG import *
-import SimPEG.DCIP as DC
+import SimPEG.EM.Static.DC as DC
 
-def run(plotIt=False):
+def run(plotIt=True):
     cs = 25.
     hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
     hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -21,10 +21,10 @@ def run(plotIt=False):
     #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
     #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
 
-    rx = DC.RxDipole(xyz_rxP, xyz_rxN)
+    rx = DC.Rx.Dipole(xyz_rxP, xyz_rxN)
-    src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
+    src = DC.Src.Dipole([rx], np.r_[-200, 0, -12.5], np.r_[+200, 0, -12.5])
-    survey = DC.SurveyDC([src])
+    survey = DC.Survey([src])
-    problem = DC.ProblemDC_CC(mesh)
+    problem = DC.Problem3D_CC(mesh)
     problem.pair(survey)
     try:
         from pymatsolver import MumpsSolver
@@ -65,4 +65,4 @@ def run(plotIt=False):
 
 
 if __name__ == '__main__':
-    print run(plotIt=True)
+    print run()

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -19,10 +19,13 @@ def run(plotIt=True):
         Morrison Casing Model, and the results are used in a 2016 SEG abstract by
         Yang et al.
 
-        - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
+        .. code-block:: text
+
+            Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
 
 
         The model consists of:
+
         - Air: Conductivity 1e-8 S/m, above z = 0
         - Background: conductivity 1e-2 S/m, below z = 0
         - Casing: conductivity 1e6 S/m
@@ -107,7 +110,7 @@ def run(plotIt=True):
     # Mesh
     mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
 
-    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
+    print 'Mesh Extent xmax: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
     print 'Number of cells', mesh.nC
 
     if plotIt is True:
@@ -215,7 +218,7 @@ def run(plotIt=True):
     # ------------ Problem and Survey ---------------
     survey = FDEM.Survey(sg_p + dg_p)
     mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
     problem.pair(survey)
 
     # ------------- Solve ---------------------------

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
 def run(plotIt=True):
     """
         MT: 1D: Inversion
-        =======================
+        =================
 
         Forward model 1D MT data.
         Setup and run a MT 1D inversion.
@@ -50,7 +50,7 @@ def run(plotIt=True):
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
-    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
     mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
 
     ## Setup the layout of the survey, set the sources and the connected receivers
@@ -76,7 +76,7 @@ def run(plotIt=True):
     survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
         fig.suptitle('Target - smooth true')
 
 

SimPEG/Examples/MT_3D_Foward.py

@@ -12,7 +12,7 @@ except:
 def run(plotIt=True, nFreq=1):
     """
         MT: 3D: Forward
-        =======================
+        ===============
 
         Forward model 3D MT data.
 
@@ -46,16 +46,15 @@ def run(plotIt=True, nFreq=1):
     survey = MT.Survey(srcList)
 
     ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
+    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
     problem.pair(survey)
-    problem.Solver = Solver
 
     # Calculate the data
     fields = problem.fields(sig)
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = MT.Data(survey,dataVec)
+    mtData = MT.Data(survey, dataVec)
     # Add plots
     if plotIt:
         pass

SimPEG/Examples/Maps_ComboMaps.py

@@ -0,0 +1,62 @@
+from SimPEG import Mesh, Maps, np
+
+def run(plotIt=True):
+    """
+
+        Maps: ComboMaps
+        ===============
+
+        We will use an example where we want a 1D layered earth as
+        our model, but we want to map this to a 2D discretization to do our forward
+        modeling. We will also assume that we are working in log conductivity still,
+        so after the transformation we want to map to conductivity space.
+        To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
+        which does the first part of what we just described. The second part will be
+        done by the :class:`SimPEG.Maps.ExpMap` described above.
+
+        .. code-block:: python
+            :linenos:
+
+            M = Mesh.TensorMesh([7,5])
+            v1dMap = Maps.SurjectVertical1D(M)
+            expMap = Maps.ExpMap(M)
+            myMap = expMap * v1dMap
+            m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+            sig = myMap * m
+
+        If you noticed, it was pretty easy to combine maps. What is even cooler is
+        that the derivatives also are made for you (if everything goes right).
+        Just to be sure that the derivative is correct, you should always run the test
+        on the mapping that you create.
+
+    """
+
+
+    M = Mesh.TensorMesh([7,5])
+    v1dMap = Maps.SurjectVertical1D(M)
+    expMap = Maps.ExpMap(M)
+    myMap = expMap * v1dMap
+    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+    sig = myMap * m
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    figs, axs = plt.subplots(1,2)
+    axs[0].plot(m, M.vectorCCy, 'b-o')
+    axs[0].set_title('Model')
+    axs[0].set_ylabel('Depth, y')
+    axs[0].set_xlabel('Value, $m_i$')
+    axs[0].set_xlim(0,3)
+    axs[0].set_ylim(0,1)
+    clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
+    axs[1].set_title('Physical Property')
+    axs[1].set_ylabel('Depth, y')
+    clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
+    plt.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Maps_Mesh2Mesh.py

@@ -0,0 +1,41 @@
+from SimPEG import Mesh, Maps, Utils
+
+def run(plotIt=True):
+    """
+
+        Maps: Mesh2Mesh
+        ===============
+
+        This mapping allows you to go from one mesh to another.
+
+    """
+
+    M = Mesh.TensorMesh([100,100])
+    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
+    h1 = h1/h1.sum()
+    M2 = Mesh.TensorMesh([h1,h1])
+    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+    v = Utils.mkvc(V)
+    modh = Maps.Mesh2Mesh([M,M2])
+    modH = Maps.Mesh2Mesh([M2,M])
+    H = modH * v
+    h = modh * H
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    ax = plt.subplot(131)
+    M.plotImage(v, ax=ax)
+    ax.set_title('Fine Mesh (Original)')
+    ax = plt.subplot(132)
+    M2.plotImage(H,clim=[0,1],ax=ax)
+    ax.set_title('Course Mesh')
+    ax = plt.subplot(133)
+    M.plotImage(h,clim=[0,1],ax=ax)
+    ax.set_title('Fine Mesh (Interpolated)')
+    plt.show()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Mesh_Operators_CahnHilliard.py

@@ -98,7 +98,7 @@ def run(plotIt=True, n=60):
             ii = int(ii)
             out = M.plotImage(PHIS[ii][1],ax=ax)
             ax.axis('off')
-            ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
+            ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii][0]))
         plt.show()
 
 if __name__ == '__main__':

SimPEG/Examples/Mesh_QuadTree_FaceDiv.py

@@ -29,15 +29,15 @@ def run(plotIt=True, n=60):
         axes[0].set_ylim([-1,17])
 
         for ii, loc in zip(range(M.nC),M.gridCC):
-            axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
+            axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='r')
 
         axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
         for ii, loc in zip(range(M.nFx),M.gridFx):
-            axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
+            axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='g')
 
         axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
         for ii, loc in zip(range(M.nFy),M.gridFy):
-            axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
+            axes[0].text(loc[0]+0.2,loc[1]+0.2,'{0:d}'.format((ii+M.nFx)), color='m')
 
         axes[1].spy(M.faceDiv)
         axes[1].set_title('Face Divergence')

SimPEG/Examples/Utils_surface2ind_topo.py

@@ -2,8 +2,12 @@ from SimPEG import *
 from SimPEG.Utils import surface2ind_topo
 
 
-def run(plotIt=False, nx = 5, ny = 5):
+def run(plotIt=True, nx=5, ny=5):
     """
+
+        Utils: surface2ind_topo
+        =======================
+
         Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
         a topographic surface.
 
@@ -13,27 +17,25 @@ def run(plotIt=False, nx = 5, ny = 5):
     xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max())
     topo = 0.4*np.sin(xtopo*5) # define a topographic surface
 
-    Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array
+    Topo = np.hstack([Utils.mkvc(xtopo,2), Utils.mkvc(topo,2)]) #make it an array
 
-    indcc = surface2ind_topo(mesh, Topo,'CC')
+    indcc = surface2ind_topo(mesh, Topo, 'CC')
 
     if plotIt:
         from matplotlib.pylab import plt
         from scipy.interpolate import interp1d
-        fig, ax = plt.subplots(1,1,figsize=(6,6))
+        fig, ax = plt.subplots(1,1, figsize=(6,6))
         mesh.plotGrid(ax=ax, nodes=True, centers=True)
         ax.plot(xtopo,topo,'k',linewidth=1)
-        # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3)
         ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3)
 
-
         aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T
         a = aveN2CC * indcc
         a[a > 0] = 1.
         a[a < 0.25] = np.nan
         a = a.reshape(mesh.vnN, order='F')
         masked_array = np.ma.array(a, mask=np.isnan(a))
-        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2)
+        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap=plt.cm.gray, alpha=0.2)
         plt.show()
 
 

SimPEG/Examples/__init__.py

@@ -10,6 +10,8 @@ import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Inversion_IRLS
 import Inversion_Linear
+import Maps_ComboMaps
+import Maps_Mesh2Mesh
 import Mesh_Basic_ForwardDC
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
@@ -22,7 +24,7 @@ import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
 
 ##### AUTOIMPORTS #####
 
@@ -38,7 +40,7 @@ if __name__ == '__main__':
 
     # Create the examples dir in the docs folder.
     fName = os.path.realpath(__file__)
-    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'content', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
@@ -57,7 +59,7 @@ if __name__ == '__main__':
         if line == "##### AUTOIMPORTS #####\n":
             inimports = not inimports
             if inimports:
-                out += '\n'.join(["import %s"%_ for _ in exfiles])
+                out += '\n'.join(["import {0!s}".format(_) for _ in exfiles])
                 out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
                 out += '\n##### AUTOIMPORTS #####\n'
     f.close()
@@ -74,11 +76,11 @@ if __name__ == '__main__':
 
         docstr = runFunction.__doc__
         if docstr is None:
-            doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
+            doc = '{0!s}\n{1!s}'.format(name.replace('_',' '), '='*len(name))
         else:
             doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
 
-        out = """.. _examples_%s:
+        out = """.. _examples_{0!s}:
 
 .. --------------------------------- ..
 ..                                   ..
@@ -88,21 +90,21 @@ if __name__ == '__main__':
 ..                                   ..
 .. --------------------------------- ..
 
-%s
+{1!s}
 
 .. plot::
 
     from SimPEG import Examples
-    Examples.%s.run()
+    Examples.{2!s}.run()
 
-.. literalinclude:: ../../SimPEG/Examples/%s.py
+.. literalinclude:: ../../../SimPEG/Examples/{3!s}.py
     :language: python
     :linenos:
-"""%(name,doc,name,name)
+""".format(name, doc, name, name)
 
-        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
+        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
 
-        print 'Creating: %s.rst'%name
+        print 'Creating: {0!s}.rst'.format(name)
         f = open(rst, 'w')
         f.write(out)
         f.close()

SimPEG/FLOW/Richards/Empirical.py

@@ -31,7 +31,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -44,7 +44,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -116,7 +116,7 @@ class RichardsMap(object):
         ax.semilogx(self.k(h, m), h)
 
     def _assertMatchesPair(self, pair):
-        assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__)
+        assert isinstance(self, pair), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
 
 
 

SimPEG/FLOW/Richards/RichardsProblem.py

@@ -140,7 +140,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         for ii, dt in enumerate(self.timeSteps):
             bc = self.getBoundaryConditions(ii, u[ii])
             u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
-            if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
+            if self.debug: print "Solving Fields ({0:4d}/{1:d} - {2:3.1f}% Done) {3:d} Iterations, {4:4.2f} seconds".format(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
         return u
 
     @Utils.timeIt

SimPEG/Fields.py

@@ -37,7 +37,7 @@ class Fields(object):
         for f in self.knownFields:
             loc =self.knownFields[f]
             sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
-        return "%e MB"%sz
+        return "{0:e} MB".format(sz)
 
     def _storageShape(self, loc):
         nSrc = self.survey.nSrc
@@ -84,12 +84,12 @@ class Fields(object):
             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)
+                raise KeyError("Invalid field name ({0!s}) for setter, you can't set an aliased property".format(name))
             else:
-                raise KeyError('Invalid field name (%s) for setter'%name)
+                raise KeyError('Invalid field name ({0!s}) for setter'.format(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)
+            raise KeyError('Invalid field name ({0!s}) for getter'.format(name))
         return name
 
     def _indexAndNameFromKey(self, key, accessType):

SimPEG/MT/SrcMT.py

@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)

SimPEG/MT/Utils/dataUtils.py

@@ -19,7 +19,7 @@ def getAppRes(MTdata):
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
-def rotateData(MTdata,rotAngle):
+def rotateData(MTdata, rotAngle):
     '''
     Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
     '''
@@ -44,19 +44,19 @@ def rotateData(MTdata,rotAngle):
     return MT.Data.fromRecArray(outRec)
 
 
-def appResPhs(freq,z):
+def appResPhs(freq, z):
     app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
     app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
     return app_res, app_phs
 
-def skindepth(rho,freq):
+def skindepth(rho, freq):
     ''' Function to calculate the skindepth of EM waves'''
     return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
 
-def rec2ndarr(x,dt=float):
+def rec2ndarr(x, dt=float):
     return x.view((dt, len(x.dtype.names)))
 
-def makeAnalyticSolution(mesh,model,elev,freqs):
+def makeAnalyticSolution(mesh, model, elev, freqs):
     from SimPEG import MT
     data1D = []
     for freq in freqs:
@@ -70,7 +70,7 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
     dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
     return dataRec
 
-def plotMT1DModelData(problem,models,symList=None):
+def plotMT1DModelData(problem, models, symList=None):
     from SimPEG import MT
     # Setup the figure
     fontSize = 15

SimPEG/Maps.py

@@ -41,8 +41,8 @@ class IdentityMap(object):
             If this is a meshless mapping (i.e. nP is defined independently)
             the shape will be the the shape (nP,nP).
 
-            :rtype: (int,int)
+            :rtype: tuple
-            :return: shape of the operator as a tuple
+            :return: shape of the operator as a tuple (int,int)
         """
         if self._nP is not None:
             return (self.nP, self.nP)
@@ -86,7 +86,7 @@ class IdentityMap(object):
             The derivative of the transformation.
 
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
         """
@@ -101,7 +101,7 @@ class IdentityMap(object):
             :return: passed the test?
 
         """
-        print 'Testing %s' % str(self)
+        print 'Testing {0!s}'.format(str(self))
         if m is None:
             m = abs(np.random.rand(self.nP))
         if 'plotIt' not in kwargs:
@@ -111,21 +111,21 @@ class IdentityMap(object):
     def _assertMatchesPair(self, pair):
         assert (isinstance(self, pair) or
             isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
-            ), "Mapping object must be an instance of a %s class."%(pair.__name__)
+            ), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
 
     def __mul__(self, val):
         if isinstance(val, IdentityMap):
             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)))
+                raise ValueError('Dimension mismatch in {0!s} and {1!s}.'.format(str(self), str(val)))
             return ComboMap([self, val])
         elif isinstance(val, np.ndarray):
             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)))
+                raise ValueError('Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(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(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
+        return "{0!s}({1!s},{2!s})".format(self.__class__.__name__, self.shape[0], self.shape[1])
 
 
 class ComboMap(IdentityMap):
@@ -140,7 +140,7 @@ class ComboMap(IdentityMap):
             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.__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)
+                raise ValueError('Dimension mismatch in map[{0!s}] ({1!s}, {2!s}) and map[{3!s}] ({4!s}, {5!s}).'.format(*errArgs))
 
             if isinstance(m, ComboMap):
                 self.maps += m.maps
@@ -173,7 +173,7 @@ class ComboMap(IdentityMap):
         return deriv
 
     def __str__(self):
-        return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
+        return 'ComboMap[{0!s}]({1!s},{2!s})'.format(' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
 
 
 class ExpMap(IdentityMap):
@@ -216,7 +216,7 @@ class ExpMap(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -366,7 +366,7 @@ class SurjectVertical1D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
         """
         repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
@@ -427,7 +427,7 @@ class Surject2Dto3D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
         """
         inds = self * np.arange(self.nP)

SimPEG/Mesh/BaseMesh.py

@@ -7,8 +7,8 @@ class BaseMesh(object):
     BaseMesh does all the counting you don't want to do.
     BaseMesh should be inherited by meshes with a regular structure.
 
-    :param numpy.array,list n: number of cells in each direction (dim, )
+    :param numpy.array n: (or list) number of cells in each direction (dim, )
-    :param numpy.array,list x0: Origin of the mesh (dim, )
+    :param numpy.array x0: (or list) Origin of the mesh (dim, )
 
     """
 
@@ -34,8 +34,8 @@ class BaseMesh(object):
         """
         Origin of the mesh
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
-        :return: x0
+        :return: x0, (dim, )
         """
         return self._x0
 
@@ -116,8 +116,8 @@ class BaseMesh(object):
         """
         Total number of edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
-        :return: [nEx, nEy, nEz]
+        :return: [nEx, nEy, nEz], (dim, )
 
         .. plot::
             :include-source:
@@ -173,8 +173,8 @@ class BaseMesh(object):
         """
         Total number of faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
-        :return: [nFx, nFy, nFz]
+        :return: [nFx, nFy, nFz], (dim, )
 
         .. plot::
             :include-source:
@@ -200,8 +200,8 @@ class BaseMesh(object):
         """
         Face Normals
 
-        :rtype: numpy.array (sum(nF), dim)
+        :rtype: numpy.array
-        :return: normals
+        :return: normals, (sum(nF), dim)
         """
         if self.dim == 2:
             nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)]
@@ -218,8 +218,8 @@ class BaseMesh(object):
         """
         Edge Tangents
 
-        :rtype: numpy.array (sum(nE), dim)
+        :rtype: numpy.array
-        :return: normals
+        :return: normals, (sum(nE), dim)
         """
         if self.dim == 2:
             tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)]
@@ -236,8 +236,9 @@ class BaseMesh(object):
         Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
 
         :param numpy.array fV: face vector with shape (nF, dim)
-        :rtype: numpy.array with shape (nF, )
+        :rtype: numpy.array
-        :return: projected face vector
+        :return: projected face vector, (nF, )
+
         """
         assert isinstance(fV, np.ndarray), 'fV must be an ndarray'
         assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
@@ -248,8 +249,9 @@ class BaseMesh(object):
         Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
 
         :param numpy.array eV: edge vector with shape (nE, dim)
-        :rtype: numpy.array with shape (nE, )
+        :rtype: numpy.array
-        :return: projected edge vector
+        :return: projected edge vector, (nE, )
+
         """
         assert isinstance(eV, np.ndarray), 'eV must be an ndarray'
         assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
@@ -295,7 +297,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of cells in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: [nCx, nCy, nCz]
         """
         return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
@@ -335,7 +337,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of nodes in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: [nNx, nNy, nNz]
         """
         return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
@@ -345,7 +347,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEx
         """
         return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
@@ -355,7 +357,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
@@ -365,7 +367,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
@@ -375,7 +377,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFx
         """
         return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
@@ -385,7 +387,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
@@ -395,7 +397,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
@@ -520,7 +522,7 @@ class BaseRectangularMesh(BaseMesh):
         assert xType in outType, 'You cannot change type of components.'
         if type(x) == list:
             for i, xi in enumerate(x):
-                assert isinstance(x, np.ndarray), "x[%i] must be a numpy array" % i
+                assert isinstance(x, np.ndarray), "x[{0:d}] must be a numpy array".format(i)
                 assert xi.size == x[0].size, "Number of elements in list must not change."
 
             x_array = np.ones((x.size, len(x)))

SimPEG/Mesh/CurvilinearMesh.py

@@ -4,14 +4,28 @@ from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 from View import CurvView
 
+
 # Some helper functions.
-length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
+def length2D(x):
-length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
+    return (x[:, 0]**2 + x[:, 1]**2)**0.5
-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 CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
+def length3D(x):
+    return (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
+
+
+def normalize2D(x):
+    return x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
+
+
+def normalize3D(x):
+    return x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
+
+
+# Curvi Mesh
+
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
+                      CurvView):
     """
     CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
@@ -31,12 +45,16 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
     _meshType = 'Curv'
 
     def __init__(self, nodes):
-        assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
+        assert type(nodes) == list, ("'nodes' variable must be a list of "
+                                     "np.ndarray")
         assert len(nodes) > 1, "len(node) must be greater than 1"
 
         for i, nodes_i in enumerate(nodes):
-            assert isinstance(nodes_i, np.ndarray), ("nodes[%i] is not a numpy array." % i)
+            assert isinstance(nodes_i, np.ndarray), ("nodes[{0:d}] is not a"
-            assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
+                                                     "numpy array.".format(i))
+            assert nodes_i.shape == nodes[0].shape, ("nodes[{0:d}] is not the "
+                                                     "same shape as nodes[0]"
+                                                     .format(i))
 
         assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
         assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
@@ -48,121 +66,113 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
         for i, node_i in enumerate(nodes):
             self._gridN[:, i] = Utils.mkvc(node_i.astype(float))
 
-    def gridCC():
+    @property
-        doc = "Cell-centered grid."
+    def gridCC(self):
+        """
+        Cell-centered grid
+        """
+        if getattr(self, '_gridCC', None) is None:
+            self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:, i]
+                                           for i in range(self.dim)]).reshape(
+                                           (-1, self.dim), order='F')
+        return self._gridCC
 
-        def fget(self):
+    @property
-            if self._gridCC is None:
+    def gridN(self):
-                self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F')
+        """
-            return self._gridCC
+        Nodal grid.
-        return locals()
+        """
-    _gridCC = None  # Store grid by default
+        if getattr(self, '_gridN', None) is None:
-    gridCC = property(**gridCC())
+            raise Exception("Someone deleted this. I blame you.")
+        return self._gridN
 
-    def gridN():
+    @property
-        doc = "Nodal grid."
+    def gridFx(self):
+        """
+        Face staggered grid in the x direction.
+        """
 
-        def fget(self):
+        if getattr(self, '_gridFx', None) is None:
-            if self._gridN is None:
+            N = self.r(self.gridN, 'N', 'N', 'M')
-                raise Exception("Someone deleted this. I blame you.")
+            if self.dim == 2:
-            return self._gridN
+                XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
-        return locals()
+                self._gridFx = np.c_[XY[0], XY[1]]
-    _gridN = None  # Store grid by default
+            elif self.dim == 3:
-    gridN = property(**gridN())
+                XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] +
+                       n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
+                self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+        return self._gridFx
 
-    def gridFx():
+    @property
-        doc = "Face staggered grid in the x direction."
+    def gridFy(self):
+        """
+        Face staggered grid in the y direction.
+        """
 
-        def fget(self):
+        if getattr(self, '_gridFy', None) is None:
-            if self._gridFx is None:
+            N = self.r(self.gridN, 'N', 'N', 'M')
-                N = self.r(self.gridN, 'N', 'N', 'M')
+            if self.dim == 2:
-                if self.dim == 2:
+                XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
-                    XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+                self._gridFy = np.c_[XY[0], XY[1]]
-                    self._gridFx = np.c_[XY[0], XY[1]]
+            elif self.dim == 3:
-                elif self.dim == 3:
+                XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] +
-                    XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
+                       n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
-                    self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+                self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-            return self._gridFx
+        return self._gridFy
-        return locals()
-    _gridFx = None  # Store grid by default
-    gridFx = property(**gridFx())
 
-    def gridFy():
+    @property
-        doc = "Face staggered grid in the y direction."
+    def gridFz(self):
+        """
+        Face staggered grid in the y direction.
+        """
 
-        def fget(self):
+        if getattr(self, '_gridFz', None) is None:
-            if self._gridFy is None:
+            N = self.r(self.gridN, 'N', 'N', 'M')
-                N = self.r(self.gridN, 'N', 'N', 'M')
+            XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] +
-                if self.dim == 2:
+                   n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
-                    XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
+            self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-                    self._gridFy = np.c_[XY[0], XY[1]]
+        return self._gridFz
-                elif self.dim == 3:
-                    XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
-                    self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-            return self._gridFy
-        return locals()
-    _gridFy = None  # Store grid by default
-    gridFy = property(**gridFy())
 
-    def gridFz():
+    @property
-        doc = "Face staggered grid in the z direction."
+    def gridEx(self):
+        """
+        Edge staggered grid in the x direction.
+        """
+        if getattr(self, '_gridEx', None) is None:
+            N = self.r(self.gridN, 'N', 'N', 'M')
+            if self.dim == 2:
+                XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
+                self._gridEx = np.c_[XY[0], XY[1]]
+            elif self.dim == 3:
+                XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
+                self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+        return self._gridEx
 
-        def fget(self):
+    @property
-            if self._gridFz is None and self.dim == 3:
+    def gridEy(self):
-                N = self.r(self.gridN, 'N', 'N', 'M')
+        """
-                XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
+        Edge staggered grid in the y direction.
-                self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+        """
-            return self._gridFz
+        if getattr(self, '_gridEy', None) is None:
-        return locals()
+            N = self.r(self.gridN, 'N', 'N', 'M')
-    _gridFz = None  # Store grid by default
+            if self.dim == 2:
-    gridFz = property(**gridFz())
+                XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+                self._gridEy = np.c_[XY[0], XY[1]]
+            elif self.dim == 3:
+                XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
+                self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+        return self._gridEy
 
-    def gridEx():
+    @property
-        doc = "Edge staggered grid in the x direction."
+    def gridEz(self):
-
+        """
-        def fget(self):
+        Edge staggered grid in the z direction.
-            if self._gridEx is None:
+        """
-                N = self.r(self.gridN, 'N', 'N', 'M')
+        if getattr(self, '_gridEz', None) is None and self.dim == 3:
-                if self.dim == 2:
+            N = self.r(self.gridN, 'N', 'N', 'M')
-                    XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
+            XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
-                    self._gridEx = np.c_[XY[0], XY[1]]
+            self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-                elif self.dim == 3:
+        return self._gridEz
-                    XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
-                    self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-            return self._gridEx
-        return locals()
-    _gridEx = None  # Store grid by default
-    gridEx = property(**gridEx())
-
-    def gridEy():
-        doc = "Edge staggered grid in the y direction."
-
-        def fget(self):
-            if self._gridEy is None:
-                N = self.r(self.gridN, 'N', 'N', 'M')
-                if self.dim == 2:
-                    XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
-                    self._gridEy = np.c_[XY[0], XY[1]]
-                elif self.dim == 3:
-                    XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
-                    self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-            return self._gridEy
-        return locals()
-    _gridEy = None  # Store grid by default
-    gridEy = property(**gridEy())
-
-    def gridEz():
-        doc = "Edge staggered grid in the z direction."
-
-        def fget(self):
-            if self._gridEz is None and self.dim == 3:
-                N = self.r(self.gridN, 'N', 'N', 'M')
-                XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
-                self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-            return self._gridEz
-        return locals()
-    _gridEz = None  # Store grid by default
-    gridEz = property(**gridEz())
 
     # --------------- Geometries ---------------------
     #
@@ -194,78 +204,94 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
     #            | /              | /
     #            D -------------- C
     #       node(i+1,j,k)      node(i+1,j+1,k)
-    def vol():
-        doc = "Construct cell volumes of the 3D model as 1d array."
 
-        def fget(self):
+    @property
-            if(self._vol is None):
+    def vol(self):
-                if self.dim == 2:
+        """
-                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
+        Construct cell volumes of the 3D model as 1d array
-                    normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
+        """
-                    self._vol = area
-                elif self.dim == 3:
-                    # Each polyhedron can be decomposed into 5 tetrahedrons
-                    # However, this presents a choice so we may as well divide in two ways and average.
-                    A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1)
 
-                    vol1 = (Utils.volTetra(self.gridN, A, B, D, E) +  # cutted edge top
+        if getattr(self, '_vol', None) is None:
-                            Utils.volTetra(self.gridN, B, E, F, G) +  # cutted edge top
+            if self.dim == 2:
-                            Utils.volTetra(self.gridN, B, D, E, G) +  # middle
+                A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
-                            Utils.volTetra(self.gridN, B, C, D, G) +  # cutted edge bottom
+                normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros(
-                            Utils.volTetra(self.gridN, D, E, G, H))   # cutted edge bottom
+                                              (self.nN, 1))], A, B, C, D)
+                self._vol = area
+            elif self.dim == 3:
+                # Each polyhedron can be decomposed into 5 tetrahedrons
+                # However, this presents a choice so we may as well divide in
+                # two ways and average.
+                A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC +
+                                                         1)
 
-                    vol2 = (Utils.volTetra(self.gridN, A, F, B, C) +  # cutted edge top
+                vol1 = (Utils.volTetra(self.gridN, A, B, D, E) +  # cutted edge top
-                            Utils.volTetra(self.gridN, A, E, F, H) +  # cutted edge top
+                        Utils.volTetra(self.gridN, B, E, F, G) +  # cutted edge top
-                            Utils.volTetra(self.gridN, A, H, F, C) +  # middle
+                        Utils.volTetra(self.gridN, B, D, E, G) +  # middle
-                            Utils.volTetra(self.gridN, C, H, D, A) +  # cutted edge bottom
+                        Utils.volTetra(self.gridN, B, C, D, G) +  # cutted edge bottom
-                            Utils.volTetra(self.gridN, C, G, H, F))   # cutted edge bottom
+                        Utils.volTetra(self.gridN, D, E, G, H))   # cutted edge bottom
 
-                    self._vol = (vol1 + vol2)/2
+                vol2 = (Utils.volTetra(self.gridN, A, F, B, C) +  # cutted edge top
-            return self._vol
+                        Utils.volTetra(self.gridN, A, E, F, H) +  # cutted edge top
-        return locals()
+                        Utils.volTetra(self.gridN, A, H, F, C) +  # middle
-    _vol = None
+                        Utils.volTetra(self.gridN, C, H, D, A) +  # cutted edge bottom
-    vol = property(**vol())
+                        Utils.volTetra(self.gridN, C, G, H, F))   # cutted edge bottom
 
-    def area():
+                self._vol = (vol1 + vol2)/2
-        doc = "Face areas."
+        return self._vol
 
-        def fget(self):
+    @property
-            if(self._area is None or self._normals is None):
+    def area(self):
-                # Compute areas of cell faces
+        if (getattr(self, '_area', None) is None or
-                if(self.dim == 2):
+            getattr(self, '_normals', None) is None):
-                    xy = self.gridN
+            # Compute areas of cell faces
-                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
+            if(self.dim == 2):
-                    edge1 = xy[B, :] - xy[A, :]
+                xy = self.gridN
-                    normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
+                A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
-                    area1 = length2D(edge1)
+                                       self.nCy]))
-                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
+                edge1 = xy[B, :] - xy[A, :]
-                    # Note that we are doing A-D to make sure the normal points the right way.
+                normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
-                    # Think about it. Look at the picture. Normal points towards C iff you do this.
+                area1 = length2D(edge1)
-                    edge2 = xy[A, :] - xy[D, :]
+                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
-                    normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
+                                       self.nNy]))
-                    area2 = length2D(edge2)
+                # Note that we are doing A-D to make sure the normal points the
-                    self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
+                # right way.
-                    self._normals = [normalize2D(normal1), normalize2D(normal2)]
+                # Think about it. Look at the picture. Normal points towards C
-                elif(self.dim == 3):
+                # iff you do this.
+                edge2 = xy[A, :] - xy[D, :]
+                normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
+                area2 = length2D(edge2)
+                self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
+                self._normals = [normalize2D(normal1), normalize2D(normal2)]
 
-                    A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz]))
+            elif(self.dim == 3):
-                    normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
 
-                    A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz]))
+                A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array(
-                    normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
+                                             [self.nNx, self.nCy, self.nCz]))
+                normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B,
+                                                average=False,
+                                                normalizeNormals=False)
 
-                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz]))
+                A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array(
-                    normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
+                                             [self.nCx, self.nNy, self.nCz]))
+                normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E,
+                                                average=False,
+                                                normalizeNormals=False)
 
-                    self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
+                A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array(
-                    self._normals = [normal1, normal2, normal3]
+                                             [self.nCx, self.nCy, self.nNz]))
-            return self._area
+                normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D,
-        return locals()
+                                                average=False,
-    _area = None
+                                                normalizeNormals=False)
-    area = property(**area())
 
-    def normals():
+                self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2),
-        doc = """Face normals: calling this will average
+                                   Utils.mkvc(area3)]
+                self._normals = [normal1, normal2, normal3]
+        return self._area
+
+    @property
+    def normals(self):
+        """
+        Face normals: calling this will average
         the computed normals so that there is one
         per face. This is especially relevant in
         3D, as there are up to 4 different normals
@@ -276,58 +302,64 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
             NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
         """
 
-        def fget(self):
+        if getattr(self, '_normals', None) is None:
-            if(self._normals is None):
+            self.area  # calling .area will create the face normals
-                self.area  # calling .area will create the face normals
+        if self.dim == 2:
-            if self.dim == 2:
+            return normalize2D(np.r_[self._normals[0], self._normals[1]])
-                return normalize2D(np.r_[self._normals[0], self._normals[1]])
+        elif self.dim == 3:
-            elif self.dim == 3:
+            normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
-                normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
+            normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
-                normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
+            normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
-                normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
+            return normalize3D(np.r_[normal1, normal2, normal3])
-                return normalize3D(np.r_[normal1, normal2, normal3])
-        return locals()
-    _normals = None
-    normals = property(**normals())
 
-    def edge():
+    @property
-        doc = "Edge legnths."
+    def edge(self):
-
+        """
-        def fget(self):
+        Edge lengths
-            if(self._edge is None or self._tangents is None):
+        """
-                if(self.dim == 2):
+        if getattr(self, '_edge', None) is None:
-                    xy = self.gridN
+            if(self.dim == 2):
-                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
+                xy = self.gridN
-                    edge1 = xy[D, :] - xy[A, :]
+                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
-                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
+                                                                  self.nNy]))
-                    edge2 = xy[B, :] - xy[A, :]
+                edge1 = xy[D, :] - xy[A, :]
-                    self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
+                A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
-                    self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
+                                                                   self.nCy]))
-                elif(self.dim == 3):
+                edge2 = xy[B, :] - xy[A, :]
-                    xyz = self.gridN
+                self._edge = np.r_[Utils.mkvc(length2D(edge1)),
-                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
+                                   Utils.mkvc(length2D(edge2))]
-                    edge1 = xyz[D, :] - xyz[A, :]
+                self._tangents = np.r_[edge1, edge2]/np.c_[self._edge,
-                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz]))
+                                                           self._edge]
-                    edge2 = xyz[B, :] - xyz[A, :]
+            elif(self.dim == 3):
-                    A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
+                xyz = self.gridN
-                    edge3 = xyz[E, :] - xyz[A, :]
+                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
-                    self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
+                                                                   self.nNy,
-                    self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
+                                                                   self.nNz]))
+                edge1 = xyz[D, :] - xyz[A, :]
+                A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
+                                                                   self.nCy,
+                                                                   self.nNz]))
+                edge2 = xyz[B, :] - xyz[A, :]
+                A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx,
+                                                                   self.nNy,
+                                                                   self.nCz]))
+                edge3 = xyz[E, :] - xyz[A, :]
+                self._edge = np.r_[Utils.mkvc(length3D(edge1)),
+                                   Utils.mkvc(length3D(edge2)),
+                                   Utils.mkvc(length3D(edge3))]
+                self._tangents = (np.r_[edge1, edge2, edge3] /
+                                  np.c_[self._edge, self._edge, self._edge])
             return self._edge
-        return locals()
+        return self._edge
-    _edge = None
-    edge = property(**edge())
 
-    def tangents():
+    @property
-        doc = "Edge tangents."
+    def tangents(self):
-
+        """
-        def fget(self):
+        Edge tangents
-            if(self._tangents is None):
+        """
-                self.edge  # calling .edge will create the tangents
+        if getattr(self, '_tangents', None) is None:
-            return self._tangents
+            self.edge  # calling .edge will create the tangents
-        return locals()
+        return self._tangents
-    _tangents = None
-    tangents = property(**tangents())
 
 
 

SimPEG/Mesh/CylMesh.py

@@ -68,8 +68,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array 
-        :return: vnFx
+        :return: vnFx, (dim, )
         """
         return self.vnC
 
@@ -78,8 +78,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array 
-        :return: vnEy or None if dim < 2
+        :return: vnEy or None if dim < 2, (dim, )
         """
         nNx = self.nNx if self.isSymmetric else self.nNx - 1
         return np.r_[nNx, self.nCy, self.nNz]
@@ -89,8 +89,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array 
-        :return: vnEz or None if nCy > 1
+        :return: vnEz or None if nCy > 1, (dim, )
         """
         if self.isSymmetric:
             return np.r_[self.nNx, self.nNy, self.nCz]

SimPEG/Mesh/DiffOperators.py

@@ -18,13 +18,15 @@ def checkBC(bc):
 
     for bc_i in bc:
         assert type(bc_i) is str, "each bc must be a string"
-        assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
+        assert bc_i in ['dirichlet', 'neumann'], ("each bc must be either,"
+                                                  "'dirichlet' or 'neumann'")
     return bc
 
 
 def ddxCellGrad(n, bc):
     """
-        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
+        Create 1D derivative operator from cell-centers to nodes this means we
+        go from n to n+1
 
         For Cell-Centered **Dirichlet**, use a ghost point::
 
@@ -52,7 +54,8 @@ def ddxCellGrad(n, bc):
     """
     bc = checkBC(bc)
 
-    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
+    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n,
+                   format="csr")
     # Set the first side
     if(bc[0] == 'dirichlet'):
         D[0, 0] = 2
@@ -65,10 +68,11 @@ def ddxCellGrad(n, bc):
         D[-1, -1] = 0
     return D
 
+
 def ddxCellGradBC(n, bc):
     """
-
+        Create 1D derivative operator from cell-centers to nodes this means we
-        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
+        go from n to n+1
 
         For Cell-Centered **Dirichlet**, use a ghost point::
 
@@ -99,7 +103,7 @@ def ddxCellGradBC(n, bc):
     """
     bc = checkBC(bc)
 
-    ij   = (np.array([0, n]),np.array([0, 1]))
+    ij   = (np.array([0, n]), np.array([0, 1]))
     vals = np.zeros(2)
 
     # Set the first side
@@ -112,7 +116,7 @@ def ddxCellGradBC(n, bc):
         vals[1] = 2
     elif(bc[1] == 'neumann'):
         vals[1] = 0
-    D = sp.csr_matrix((vals, ij), shape=(n+1,2))
+    D = sp.csr_matrix((vals, ij), shape=(n+1, 2))
     return D
 
 
@@ -121,175 +125,166 @@ class DiffOperators(object):
         Class creates the differential operators that you need!
     """
     def __init__(self):
-        raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
+        raise Exception('DiffOperators is a base class providing differential'
+                        'operators on meshes and cannot run on its own.'
+                        'Inherit to your favorite Mesh class.')
 
-    def faceDiv():
+    @property
-        doc = "Construct divergence operator (face-stg to cell-centres)."
+    def faceDiv(self):
-
+        """
-        def fget(self):
+        Construct divergence operator (face-stg to cell-centres).
-            if(self._faceDiv is None):
+        """
-                # The number of cell centers in each direction
+        if getattr(self, '_faceDiv', None) is None:
-                n = self.vnC
+            n = self.vnC
-                # Compute faceDivergence operator on faces
+            # Compute faceDivergence operator on faces
-                if(self.dim == 1):
+            if(self.dim == 1):
-                    D = ddx(n[0])
+                D = ddx(n[0])
-                elif(self.dim == 2):
+            elif(self.dim == 2):
-                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
+                D1 = sp.kron(speye(n[1]), ddx(n[0]))
-                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
+                D2 = sp.kron(ddx(n[1]), speye(n[0]))
-                    D = sp.hstack((D1, D2), format="csr")
+                D = sp.hstack((D1, D2), format="csr")
-                elif(self.dim == 3):
+            elif(self.dim == 3):
-                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
+                D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
-                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
+                D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
-                    D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
-                    D = sp.hstack((D1, D2, D3), format="csr")
-                # Compute areas of cell faces & volumes
-                S = self.area
-                V = self.vol
-                self._faceDiv = sdiag(1/V)*D*sdiag(S)
-
-            return self._faceDiv
-        return locals()
-    _faceDiv = None
-    faceDiv = property(**faceDiv())
-
-    def faceDivx():
-        doc = "Construct divergence operator in the x component (face-stg to cell-centres)."
-
-        def fget(self):
-            if(self._faceDivx is None):
-                # The number of cell centers in each direction
-                n = self.vnC
-                # Compute faceDivergence operator on faces
-                if(self.dim == 1):
-                    D1 = ddx(n[0])
-                elif(self.dim == 2):
-                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
-                elif(self.dim == 3):
-                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
-                # Compute areas of cell faces & volumes
-                S = self.r(self.area, 'F', 'Fx', 'V')
-                V = self.vol
-                self._faceDivx = sdiag(1/V)*D1*sdiag(S)
-
-            return self._faceDivx
-        return locals()
-    _faceDivx = None
-    faceDivx = property(**faceDivx())
-
-    def faceDivy():
-        doc = "Construct divergence operator in the y component (face-stg to cell-centres)."
-
-        def fget(self):
-            if(self.dim < 2): return None
-            if(self._faceDivy is None):
-                # The number of cell centers in each direction
-                n = self.vnC
-                # Compute faceDivergence operator on faces
-                if(self.dim == 2):
-                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
-                elif(self.dim == 3):
-                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
-                # Compute areas of cell faces & volumes
-                S = self.r(self.area, 'F', 'Fy', 'V')
-                V = self.vol
-                self._faceDivy = sdiag(1/V)*D2*sdiag(S)
-
-            return self._faceDivy
-        return locals()
-    _faceDivy = None
-    faceDivy = property(**faceDivy())
-
-    def faceDivz():
-        doc = "Construct divergence operator in the z component (face-stg to cell-centres)."
-
-        def fget(self):
-            if(self.dim < 3): return None
-            if(self._faceDivz is None):
-                # The number of cell centers in each direction
-                n = self.vnC
-                # Compute faceDivergence operator on faces
                 D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
-                # Compute areas of cell faces & volumes
+                D = sp.hstack((D1, D2, D3), format="csr")
-                S = self.r(self.area, 'F', 'Fz', 'V')
+            # Compute areas of cell faces & volumes
-                V = self.vol
+            S = self.area
-                self._faceDivz = sdiag(1/V)*D3*sdiag(S)
+            V = self.vol
+            self._faceDiv = sdiag(1/V)*D*sdiag(S)
+        return self._faceDiv
 
-            return self._faceDivz
+    @property
-        return locals()
+    def faceDivx(self):
-    _faceDivz = None
+        """
-    faceDivz = property(**faceDivz())
+        Construct divergence operator in the x component (face-stg to
+        cell-centres).
+        """
+        if getattr(self, '_faceDivx', None) is None:
+            # The number of cell centers in each direction
+            n = self.vnC
+            # Compute faceDivergence operator on faces
+            if(self.dim == 1):
+                D1 = ddx(n[0])
+            elif(self.dim == 2):
+                D1 = sp.kron(speye(n[1]), ddx(n[0]))
+            elif(self.dim == 3):
+                D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
+            # Compute areas of cell faces & volumes
+            S = self.r(self.area, 'F', 'Fx', 'V')
+            V = self.vol
+            self._faceDivx = sdiag(1/V)*D1*sdiag(S)
 
-    def nodalGrad():
+        return self._faceDivx
-        doc = "Construct gradient operator (nodes to edges)."
 
-        def fget(self):
+    @property
-            if(self._nodalGrad is None):
+    def faceDivy(self):
-                # The number of cell centers in each direction
+        if(self.dim < 2):
-                n = self.vnC
+            return None
-                # Compute divergence operator on faces
+        if getattr(self, '_faceDivy', None) is None:
-                if(self.dim == 1):
+            # The number of cell centers in each direction
-                    G = ddx(n[0])
+            n = self.vnC
-                elif(self.dim == 2):
+            # Compute faceDivergence operator on faces
-                    D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
+            if(self.dim == 2):
-                    D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
+                D2 = sp.kron(ddx(n[1]), speye(n[0]))
-                    G = sp.vstack((D1, D2), format="csr")
+            elif(self.dim == 3):
-                elif(self.dim == 3):
+                D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
-                    D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
+            # Compute areas of cell faces & volumes
-                    D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
+            S = self.r(self.area, 'F', 'Fy', 'V')
-                    D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
+            V = self.vol
-                    G = sp.vstack((D1, D2, D3), format="csr")
+            self._faceDivy = sdiag(1/V)*D2*sdiag(S)
-                # Compute lengths of cell edges
+        return self._faceDivy
-                L = self.edge
-                self._nodalGrad = sdiag(1/L)*G
-            return self._nodalGrad
-        return locals()
-    _nodalGrad = None
-    nodalGrad = property(**nodalGrad())
 
-    def nodalLaplacian():
+    @property
-        doc = "Construct laplacian operator (nodes to edges)."
+    def faceDivz(self):
+        """
+        Construct divergence operator in the z component (face-stg to
+        cell-centres).
+        """
+        if(self.dim < 3):
+            return None
+        if getattr(self, '_faceDivz', None) is None:
+            # The number of cell centers in each direction
+            n = self.vnC
+            # Compute faceDivergence operator on faces
+            D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
+            # Compute areas of cell faces & volumes
+            S = self.r(self.area, 'F', 'Fz', 'V')
+            V = self.vol
+            self._faceDivz = sdiag(1/V)*D3*sdiag(S)
+        return self._faceDivz
 
-        def fget(self):
+    @property
-            if(self._nodalLaplacian is None):
+    def nodalGrad(self):
-                print 'Warning: Laplacian has not been tested rigorously.'
+        """
-                # The number of cell centers in each direction
+        Construct gradient operator (nodes to edges).
-                n = self.vnC
+        """
-                # Compute divergence operator on faces
+        if getattr(self, '_nodalGrad', None) is None:
-                if(self.dim == 1):
+            # The number of cell centers in each direction
-                    D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
+            n = self.vnC
-                    L  = - D1.T*D1
+            # Compute divergence operator on faces
-                elif(self.dim == 2):
+            if(self.dim == 1):
-                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                G = ddx(n[0])
-                    D2 = sdiag(1./self.hy) * ddx(n[1])
+            elif(self.dim == 2):
-                    L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
+                D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
-                    L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
+                D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
-                    L  = L1 + L2
+                G = sp.vstack((D1, D2), format="csr")
-                elif(self.dim == 3):
+            elif(self.dim == 3):
-                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
-                    D2 = sdiag(1./self.hy) * ddx(n[1])
+                D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
-                    D3 = sdiag(1./self.hz) * ddx(n[2])
+                D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
-                    L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
+                G = sp.vstack((D1, D2, D3), format="csr")
-                    L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
+            # Compute lengths of cell edges
-                    L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
+            L = self.edge
-                    L  = L1 + L2 + L3
+            self._nodalGrad = sdiag(1/L)*G
-                self._nodalLaplacian = L
+        return self._nodalGrad
-            return self._nodalLaplacian
+
-        return locals()
+    @property
-    _nodalLaplacian = None
+    def nodalLaplacian(self):
-    nodalLaplacian = property(**nodalLaplacian())
+        """
+        Construct laplacian operator (nodes to edges).
+        """
+        if getattr(self, '_nodalLaplacian', None) is None:
+            print 'Warning: Laplacian has not been tested rigorously.'
+            # The number of cell centers in each direction
+            n = self.vnC
+            # Compute divergence operator on faces
+            if(self.dim == 1):
+                D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
+                L = - D1.T*D1
+            elif(self.dim == 2):
+                D1 = sdiag(1./self.hx) * ddx(n[0])
+                D2 = sdiag(1./self.hy) * ddx(n[1])
+                L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
+                L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
+                L = L1 + L2
+            elif(self.dim == 3):
+                D1 = sdiag(1./self.hx) * ddx(n[0])
+                D2 = sdiag(1./self.hy) * ddx(n[1])
+                D3 = sdiag(1./self.hz) * ddx(n[2])
+                L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
+                L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
+                L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
+                L = L1 + L2 + L3
+            self._nodalLaplacian = L
+        return self._nodalLaplacian
 
     def setCellGradBC(self, BC):
         """
-        Function that sets the boundary conditions for cell-centred derivative operators.
+        Function that sets the boundary conditions for cell-centred derivative
+        operators.
 
         Examples::
+            # Neumann in all directions
+            BC = 'neumann'
 
-            BC = 'neumann'                                            # Neumann in all directions
+            # 3D, Dirichlet in y Neumann else
-            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
+            BC = ['neumann', 'dirichlet', 'neumann']
-            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
-                                                                      #     Dirichlet else
 
+            # 3D, Neumann in x on bottom of domain,  Dirichlet else
+            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
         """
+
         if(type(BC) is str):
             BC = [BC]*self.dim
         if(type(BC) is list):
@@ -323,47 +318,69 @@ class DiffOperators(object):
             G = sp.vstack((G1, G2, G3), format="csr")
         return G
 
-    def cellGrad():
+    @property
-        doc = "The cell centered Gradient, takes you to cell faces."
+    def cellGrad(self):
+        """
+        The cell centered Gradient, takes you to cell faces.
+        """
+        if getattr(self, '_cellGrad', None) is None:
+            G = self._cellGradStencil()
+            S = self.area  # Compute areas of cell faces & volumes
+            V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+            self._cellGrad = sdiag(S/V)*G
+        return self._cellGrad
 
-        def fget(self):
+    @property
-            if(self._cellGrad is None):
+    def cellGradBC(self):
-                G = self._cellGradStencil()
+        """
-                # Compute areas of cell faces & volumes
+        The cell centered Gradient boundary condition matrix
-                S = self.area
+        """
-                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+        if getattr(self, '_cellGradBC', None) is None:
-                self._cellGrad = sdiag(S/V)*G
+            BC = self.setCellGradBC(self._cellGradBC_list)
-            return self._cellGrad
+            n = self.vnC
-        return locals()
+            if(self.dim == 1):
-    _cellGrad = None
+                G = ddxCellGradBC(n[0], BC[0])
-    cellGrad = property(**cellGrad())
+            elif(self.dim == 2):
+                G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+                G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+                G = sp.block_diag((G1, G2), format="csr")
+            elif(self.dim == 3):
+                G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+                G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+                G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
+                G = sp.block_diag((G1, G2, G3), format="csr")
+            # Compute areas of cell faces & volumes
+            S = self.area
+            V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+            self._cellGradBC = sdiag(S/V)*G
+        return self._cellGradBC
 
-    def cellGradBC():
+    # def cellGradBC():
-        doc = "The cell centered Gradient boundary condition matrix"
+    #     doc = "The cell centered Gradient boundary condition matrix"
 
-        def fget(self):
+    #     def fget(self):
-            if(self._cellGradBC is None):
+    #         if(self._cellGradBC is None):
-                BC = self.setCellGradBC(self._cellGradBC_list)
+    #             BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.vnC
+    #             n = self.vnC
-                if(self.dim == 1):
+    #             if(self.dim == 1):
-                    G = ddxCellGradBC(n[0], BC[0])
+    #                 G = ddxCellGradBC(n[0], BC[0])
-                elif(self.dim == 2):
+    #             elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+    #                 G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
-                    G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+    #                 G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
-                    G = sp.block_diag((G1, G2), format="csr")
+    #                 G = sp.block_diag((G1, G2), format="csr")
-                elif(self.dim == 3):
+    #             elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+    #                 G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
-                    G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+    #                 G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
-                    G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
+    #                 G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
-                    G = sp.block_diag((G1, G2, G3), format="csr")
+    #                 G = sp.block_diag((G1, G2, G3), format="csr")
-                # Compute areas of cell faces & volumes
+    #             # Compute areas of cell faces & volumes
-                S = self.area
+    #             S = self.area
-                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+    #             V = self.aveCC2F*self.vol  # Average volume between adjacent cells
-                self._cellGradBC = sdiag(S/V)*G
+    #             self._cellGradBC = sdiag(S/V)*G
-            return self._cellGradBC
+    #         return self._cellGradBC
-        return locals()
+    #     return locals()
-    _cellGradBC = None
+    # _cellGradBC = None
-    cellGradBC = property(**cellGradBC())
+    # cellGradBC = property(**cellGradBC())
 
     def _cellGradxStencil(self):
         BC = ['neumann', 'neumann']
@@ -376,20 +393,19 @@ class DiffOperators(object):
             G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
         return G1
 
-
+    @property
-    def cellGradx():
+    def cellGradx(self):
-        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+        """
-
+        Cell centered Gradient in the x dimension. Has neumann boundary
-        def fget(self):
+        conditions.
-            if getattr(self, '_cellGradx', None) is None:
+        """
-                G1 = self._cellGradxStencil()
+        if getattr(self, '_cellGradx', None) is None:
-                # Compute areas of cell faces & volumes
+            G1 = self._cellGradxStencil()
-                V = self.aveCC2F*self.vol
+            # Compute areas of cell faces & volumes
-                L = self.r(self.area/V, 'F','Fx', 'V')
+            V = self.aveCC2F*self.vol
-                self._cellGradx = sdiag(L)*G1
+            L = self.r(self.area/V, 'F','Fx', 'V')
-            return self._cellGradx
+            self._cellGradx = sdiag(L)*G1
-        return locals()
+        return self._cellGradx
-    cellGradx = property(**cellGradx())
 
     def _cellGradyStencil(self):
         if self.dim < 2: return None
@@ -401,19 +417,17 @@ class DiffOperators(object):
             G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
         return G2
 
-    def cellGrady():
+    @property
-        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+    def cellGrady(self):
-        def fget(self):
+        if self.dim < 2:
-            if self.dim < 2: return None
+            return None
-            if getattr(self, '_cellGrady', None) is None:
+        if getattr(self, '_cellGrady', None) is None:
-                G2 = self._cellGradyStencil()
+            G2 = self._cellGradyStencil()
-                # Compute areas of cell faces & volumes
+            # Compute areas of cell faces & volumes
-                V = self.aveCC2F*self.vol
+            V = self.aveCC2F*self.vol
-                L = self.r(self.area/V, 'F','Fy', 'V')
+            L = self.r(self.area/V, 'F', 'Fy', 'V')
-                self._cellGrady = sdiag(L)*G2
+            self._cellGrady = sdiag(L)*G2
-            return self._cellGrady
+        return self._cellGrady
-        return locals()
-    cellGrady = property(**cellGrady())
 
     def _cellGradzStencil(self):
         if self.dim < 3: return None
@@ -422,66 +436,61 @@ class DiffOperators(object):
         G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
         return G3
 
-    def cellGradz():
+    @property
-        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+    def cellGradz(self):
-        def fget(self):
+        """
-            if self.dim < 3: return None
+        Cell centered Gradient in the x dimension. Has neumann boundary
-            if getattr(self, '_cellGradz', None) is None:
+        conditions.
-                G3 = self._cellGradzStencil()
+        """
-                # Compute areas of cell faces & volumes
+        if self.dim < 3:
-                V = self.aveCC2F*self.vol
+            return None
-                L = self.r(self.area/V, 'F','Fz', 'V')
+        if getattr(self, '_cellGradz', None) is None:
-                self._cellGradz = sdiag(L)*G3
+            G3 = self._cellGradzStencil()
-            return self._cellGradz
+            # Compute areas of cell faces & volumes
-        return locals()
+            V = self.aveCC2F*self.vol
-    cellGradz = property(**cellGradz())
+            L = self.r(self.area/V, 'F', 'Fz', 'V')
+            self._cellGradz = sdiag(L)*G3
+        return self._cellGradz
 
-    def edgeCurl():
+    @property
-        doc = "Construct the 3D curl operator."
+    def edgeCurl(self):
+        """
+        Construct the 3D curl operator.
+        """
+        if getattr(self, '_edgeCurl', None) is None:
+            assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
 
-        def fget(self):
+            n = self.vnC  # The number of cell centers in each direction
-            if(self._edgeCurl is None):
+            L = self.edge  # Compute lengths of cell edges
-                assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
+            S = self.area # Compute areas of cell faces
-                # The number of cell centers in each direction
-                n = self.vnC
 
-                # Compute lengths of cell edges
+            # Compute divergence operator on faces
-                L = self.edge
+            if self.dim == 2:
 
-                # Compute areas of cell faces
+                D21 = sp.kron(ddx(n[1]), speye(n[0]))
-                S = self.area
+                D12 = sp.kron(speye(n[1]), ddx(n[0]))
+                C = sp.hstack((-D21, D12), format="csr")
+                self._edgeCurl = C*sdiag(1/S)
 
-                # Compute divergence operator on faces
+            elif self.dim == 3:
-                if self.dim == 2:
 
-                    D21 = sp.kron(ddx(n[1]), speye(n[0]))
+                D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
-                    D12 = sp.kron(speye(n[1]), ddx(n[0]))
+                D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
-                    C = sp.hstack((-D21, D12), format="csr")
+                D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
-                    self._edgeCurl = C*sdiag(1/S)
+                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]))
 
-                elif self.dim == 3:
+                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])
 
-                    D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
+                C = sp.vstack((sp.hstack((O1, -D32, D23)),
-                    D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
+                               sp.hstack((D31, O2, -D13)),
-                    D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
+                               sp.hstack((-D21, D12, O3))), format="csr")
-                    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])
+                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
-                    O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
+        return self._edgeCurl
-                    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))
-
-            return self._edgeCurl
-        return locals()
-    _edgeCurl = None
-    edgeCurl = property(**edgeCurl())
 
     def getBCProjWF(self, BC, discretization='CC'):
         """
@@ -489,16 +498,19 @@ class DiffOperators(object):
         The weak form boundary condition projection matrices.
 
         Examples::
+            # Neumann in all directions
+            BC = 'neumann'
 
-            BC = 'neumann'                                            # Neumann in all directions
+            # 3D, Dirichlet in y Neumann else
-            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
+            BC = ['neumann', 'dirichlet', 'neumann']
-            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
-                                                                      #     Dirichlet else
 
+            # 3D, Neumann in x on bottom of domain, Dirichlet else
+            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
         """
 
         if discretization is not 'CC':
-            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+            raise NotImplementedError('Boundary conditions only implemented'
+                                      'for CC discretization.')
 
         if(type(BC) is str):
             BC = [BC for _ in self.vnC]  # Repeat the str self.dim times
@@ -510,35 +522,34 @@ class DiffOperators(object):
         for i, bc_i in enumerate(BC):
             BC[i] = checkBC(bc_i)
 
-
         def projDirichlet(n, bc):
             bc = checkBC(bc)
-            ij = ([0,n], [0,1])
+            ij = ([0, n], [0, 1])
-            vals = [0,0]
+            vals = [0, 0]
             if(bc[0] == 'dirichlet'):
                 vals[0] = -1
             if(bc[1] == 'dirichlet'):
                 vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+            return sp.csr_matrix((vals, ij), shape=(n+1, 2))
 
         def projNeumannIn(n, bc):
             bc = checkBC(bc)
             P = sp.identity(n+1).tocsr()
             if(bc[0] == 'neumann'):
-                P = P[1:,:]
+                P = P[1:, :]
             if(bc[1] == 'neumann'):
-                P = P[:-1,:]
+                P = P[:-1, :]
             return P
 
         def projNeumannOut(n, bc):
             bc = checkBC(bc)
-            ij   = ([0, 1],[0, n])
+            ij   = ([0, 1], [0, n])
             vals = [0,0]
             if(bc[0] == 'neumann'):
                 vals[0] = 1
             if(bc[1] == 'neumann'):
                 vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(2,n+1))
+            return sp.csr_matrix((vals, ij), shape=(2, n+1))
 
         n = self.vnC
         indF = self.faceBoundaryInd
@@ -550,6 +561,7 @@ class DiffOperators(object):
             Pin = projNeumannIn(n[0], BC[0])
 
             Pout = projNeumannOut(n[0], BC[0])
+
         elif(self.dim == 2):
             Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
             Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
@@ -564,12 +576,14 @@ class DiffOperators(object):
             P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
             P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
             Pout = sp.block_diag((P1, P2), format="csr")
+
         elif(self.dim == 3):
             Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
             Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
             Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
             Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
-            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] |
+                          indF[5])]
             Pbc = Pbc*sdiag(self.area[indF])
 
             P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
@@ -586,36 +600,36 @@ class DiffOperators(object):
 
     def getBCProjWF_simple(self, discretization='CC'):
         """
-
         The weak form boundary condition projection matrices
         when mixed boundary condition is used
-
-
         """
 
         if discretization is not 'CC':
-            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+            raise NotImplementedError('Boundary conditions only implemented'
+                                      'for CC discretization.')
 
         def projBC(n):
-            ij = ([0,n], [0,1])
+            ij = ([0, n], [0, 1])
-            vals = [0,0]
+            vals = [0, 0]
             vals[0] = 1
             vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+            return sp.csr_matrix((vals, ij), shape=(n+1, 2))
 
         def projDirichlet(n, bc):
             bc = checkBC(bc)
-            ij = ([0,n], [0,1])
+            ij = ([0, n], [0, 1])
-            vals = [0,0]
+            vals = [0, 0]
             if(bc[0] == 'dirichlet'):
                 vals[0] = -1
             if(bc[1] == 'dirichlet'):
                 vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+            return sp.csr_matrix((vals, ij), shape=(n+1, 2))
 
-        BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
+        BC = [['dirichlet', 'dirichlet'], ['dirichlet', 'dirichlet'],
+              ['dirichlet', 'dirichlet']]
         n = self.vnC
         indF = self.faceBoundaryInd
+
         if(self.dim == 1):
             Pbc = projDirichlet(n[0], BC[0])
             B = projBC(n[0])
@@ -653,9 +667,11 @@ class DiffOperators(object):
         if(self.dim == 1):
             return self.aveFx2CC
         elif(self.dim == 2):
-            return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
+            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")
+            return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC,
+                                      self.aveFz2CC), format="csr")
 
     @property
     def aveF2CCV(self):
@@ -665,11 +681,16 @@ class DiffOperators(object):
         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")
+            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."
+        """
+        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):
@@ -682,8 +703,12 @@ class DiffOperators(object):
 
     @property
     def aveFy2CC(self):
-        "Construct the averaging operator on cell faces in the y direction to cell centers."
+        """
-        if self.dim < 2: return None
+        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):
@@ -694,7 +719,10 @@ class DiffOperators(object):
 
     @property
     def aveFz2CC(self):
-        "Construct the averaging operator on cell faces in the z direction to cell centers."
+        """
+        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
@@ -711,12 +739,18 @@ class DiffOperators(object):
             if(self.dim == 1):
                 self._aveCC2F = avExtrap(n[0])
             elif(self.dim == 2):
-                self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
+                self._aveCC2F = sp.vstack((sp.kron(speye(n[1]),
-                                           sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
+                                           avExtrap(n[0])),
+                                           sp.kron(avExtrap(n[1]),
+                                           speye(n[0]))), format="csr")
             elif(self.dim == 3):
-                self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
+                self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]),
-                                           kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
+                                                 avExtrap(n[0])),
-                                           kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
+                                           kron3(speye(n[2]), avExtrap(n[1]),
+                                                 speye(n[0])),
+                                           kron3(avExtrap(n[2]), speye(n[1]),
+                                                 speye(n[0]))),
+                                          format="csr")
         return self._aveCC2F
 
     @property
@@ -727,7 +761,8 @@ class DiffOperators(object):
         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")
+            return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC,
+                                     self.aveEz2CC), format="csr")
 
     @property
     def aveE2CCV(self):
@@ -737,11 +772,15 @@ class DiffOperators(object):
         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")
+            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."
+        """
+        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
@@ -755,8 +794,12 @@ class DiffOperators(object):
 
     @property
     def aveEy2CC(self):
-        "Construct the averaging operator on cell edges in the y direction to cell centers."
+        """
-        if self.dim < 2: return None
+        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
@@ -768,8 +811,12 @@ class DiffOperators(object):
 
     @property
     def aveEz2CC(self):
-        "Construct the averaging operator on cell edges in the z direction to cell centers."
+        """
-        if self.dim < 3: return None
+        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
@@ -793,7 +840,10 @@ class DiffOperators(object):
 
     @property
     def aveN2E(self):
-        "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
+        """
+        Construct the averaging operator on cell nodes to cell edges, keeping
+        each dimension separate.
+        """
 
         if getattr(self, '_aveN2E', None) is None:
             # The number of cell centers in each direction
@@ -802,16 +852,24 @@ class DiffOperators(object):
                 self._aveN2E = av(n[0])
             elif(self.dim == 2):
                 self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
-                                          sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
+                                          sp.kron(av(n[1]), speye(n[0]+1))),
+                                         format="csr")
             elif(self.dim == 3):
-                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
+                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1),
-                                          kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
+                                                av(n[0])),
-                                          kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
+                                          kron3(speye(n[2]+1), av(n[1]),
+                                                speye(n[0]+1)),
+                                          kron3(av(n[2]), speye(n[1]+1),
+                                                speye(n[0]+1))),
+                                         format="csr")
         return self._aveN2E
 
     @property
     def aveN2F(self):
-        "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
+        """
+        Construct the averaging operator on cell nodes to cell faces, keeping
+        each dimension separate.
+        """
         if getattr(self, '_aveN2F', None) is None:
             # The number of cell centers in each direction
             n = self.vnC
@@ -819,9 +877,14 @@ class DiffOperators(object):
                 self._aveN2F = av(n[0])
             elif(self.dim == 2):
                 self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
-                                          sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
+                                          sp.kron(speye(n[1]+1), av(n[0]))),
+                                         format="csr")
             elif(self.dim == 3):
-                self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
+                self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]),
-                                          kron3(av(n[2]), speye(n[1]+1), av(n[0])),
+                                                speye(n[0]+1)),
-                                          kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
+                                          kron3(av(n[2]), speye(n[1]+1),
+                                                av(n[0])),
+                                          kron3(speye(n[2]+1), av(n[1]),
+                                                av(n[0]))),
+                                         format="csr")
         return self._aveN2F

SimPEG/Mesh/InnerProducts.py

@@ -16,7 +16,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
         return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -27,7 +27,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
         return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -39,7 +39,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
@@ -115,13 +115,12 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: function
             :return: dMdmu(u), the derivative of the inner product matrix (u)
 
             Given u, dMdmu returns (nF, nC*nA)
 
-            :param np.ndarray u: vector that multiplies dMdmu
+            :param numpy.ndarray u: vector that multiplies dMdmu
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
         return self._getInnerProductDeriv(prop, 'F', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -133,7 +132,7 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
         return self._getInnerProductDeriv(prop, 'E', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -145,7 +144,7 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
         fast = None
@@ -169,7 +168,7 @@ class InnerProducts(object):
             :param numpy.array v: vector to multiply (required in the general implementation)
             :param list P: list of projection matrices
             :param str projType: 'F' for faces 'E' for edges
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (n, nC*nA)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
@@ -422,7 +421,7 @@ class InnerProducts(object):
     def _getEdgePx(M):
         """Returns a function for creating projection matrices"""
         def Px(xEdge):
-            assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
+            assert xEdge == 'eX0', 'xEdge = {0!s}, not eX0'.format(xEdge)
             return sp.identity(M.nC)
         return Px
 

SimPEG/Mesh/MeshIO.py

@@ -6,13 +6,11 @@ class TensorMeshIO(object):
     @classmethod
     def readUBC(TensorMesh, fileName):
         """
-            Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
+            Read UBC GIF 3D tensor mesh and generate 3D TensorMesh in SimPEG.
 
-            Input:
+            :param string fileName: path to the UBC GIF mesh file
-            :param fileName, path to the UBC GIF mesh file
+            :rtype: TensorMesh
-
+            :return: The tensor mesh for the fileName.
-            Output:
-            :param SimPEG TensorMesh object
         """
 
         # Interal function to read cell size lines for the UBC mesh files.
@@ -48,11 +46,9 @@ class TensorMeshIO(object):
             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
+            :param string fileName: path to the vtr model file to read
-
+            :rtype: tuple
-            Output:
+            :return: (TensorMesh, modelDictionary)
-            :return SimPEG TensorMesh object
-            :return SimPEG model dictionary
 
         """
         # Import
@@ -102,9 +98,8 @@ class TensorMeshIO(object):
         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 string fileName: path to the output vtk file
-        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+        :param dict models: dictionary of numpy.array - Name('s) and array('s). Match number of cells
-        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
 
         """
         # Import
@@ -162,12 +157,9 @@ class TensorMeshIO(object):
         """
             Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
 
-            Input:
+            :param string fileName: path to the UBC GIF mesh file to read
-            :param fileName, path to the UBC GIF mesh file to read
+            :rtype: numpy.ndarray
-            :param mesh, TensorMesh object, mesh that coresponds to the model
+            :return: model with TensorMesh ordered
-
-            Output:
-            :return numpy array, model with TensorMesh ordered
         """
         f = open(fileName, 'r')
         model = np.array(map(float, f.readlines()))
@@ -183,8 +175,7 @@ class TensorMeshIO(object):
             Writes a model associated with a SimPEG TensorMesh
             to a UBC-GIF format model file.
 
-            :param str fileName: File to write to
+            :param string fileName: File to write to
-            :param simpeg.Mesh.TensorMesh mesh: The mesh
             :param numpy.ndarray model: The model
         """
 
@@ -201,17 +192,17 @@ class TensorMeshIO(object):
         """
             Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
 
-            :param str fileName: File to write to
+            :param string fileName: File to write to
-            :param simpeg.Mesh.TensorMesh mesh: The mesh
+            :param dict models: A dictionary of the models
 
         """
         assert mesh.dim == 3
         s = ''
-        s += '%i %i %i\n' %tuple(mesh.vnC)
+        s += '{0:d} {1:d} {2:d}\n'.format(*tuple(mesh.vnC))
         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
 
-        s += '%.2f %.2f %.2f\n' %tuple(origin)
+        s += '{0:.2f} {1:.2f} {2:.2f}\n'.format(*tuple(origin))
         s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
         s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
         s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
@@ -231,9 +222,8 @@ class TreeMeshIO(object):
         """
             Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
 
-            :param str fileName: File to write to
+            :param string fileName: File to write to
-            :param simpeg.Mesh.TreeMesh mesh: The mesh
+            :param dict models: The models in a dictionary, where the keys is the name of the of the model file
-            :param dictionary models: The models in a dictionary, where the keys is the name of the of the model file
         """
 
         # Calculate information to write in the file.
@@ -286,10 +276,9 @@ class TreeMeshIO(object):
 
             Input:
             :param str meshFile: path to the UBC GIF OcTree mesh file to read
+            :rtype: SimPEG.Mesh.TreeMesh
+            :return: The octree mesh
 
-            Output:
-            :return SimPEG.Mesh.TreeMesh mesh: The octree mesh
-            :return list of ndarray's: models as a list of numpy array's
         """
 
         ## Read the file lines
@@ -335,11 +324,9 @@ class TreeMeshIO(object):
         """
             Read UBC OcTree model and get vector
 
-            Input:
+            :param string fileName: path to the UBC GIF model file to read
-            :param fileName, path to the UBC GIF model file to read
+            :rtype: numpy.ndarray
-
+            :return: OcTree model
-            Output:
-            :return numpy array, OcTree model
         """
 
         if type(fileName) is list:

SimPEG/Mesh/TensorMesh.py

@@ -23,8 +23,8 @@ class BaseTensorMesh(BaseMesh):
                 h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i)
             elif type(h_i) is list:
                 h_i = Utils.meshTensor(h_i)
-            assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i)
+            assert isinstance(h_i, np.ndarray), ("h[{0:d}] is not a numpy array.".format(i))
-            assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
+            assert len(h_i.shape) == 1, ("h[{0:d}] must be a 1D numpy array.".format(i))
             h[i] = h_i[:] # make a copy.
 
         x0 = np.zeros(len(h))
@@ -41,7 +41,7 @@ class BaseTensorMesh(BaseMesh):
                 elif x_i == 'N':
                     x0[i] = -h_i.sum()
                 else:
-                    raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i)
+                    raise Exception("x0[{0:d}] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative.".format(i))
 
         if isinstance(self, BaseRectangularMesh):
             BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
@@ -198,8 +198,8 @@ class BaseTensorMesh(BaseMesh):
         Determines if a set of points are inside a mesh.
 
         :param numpy.ndarray pts: Location of points to test
-        :rtype numpy.ndarray
+        :rtype numpy.ndarray:
-        :return inside, numpy array of booleans
+        :return: inside, numpy array of booleans
         """
         pts = Utils.asArray_N_x_Dim(pts, self.dim)
 
@@ -221,7 +221,7 @@ class BaseTensorMesh(BaseMesh):
 
         :param numpy.ndarray loc: Location of points to interpolate to
         :param str locType: What to interpolate (see below)
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: M, the interpolation matrix
 
         locType can be::
@@ -239,7 +239,7 @@ class BaseTensorMesh(BaseMesh):
             'CCVz'  -> z-component of vector field defined on cell centers
         """
         if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
-            raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
+            raise Exception('Symmetric CylMesh does not support {0!s} interpolation, as this variable does not exist.'.format(locType))
 
         loc = Utils.asArray_N_x_Dim(loc, self.dim)
 
@@ -289,7 +289,7 @@ class BaseTensorMesh(BaseMesh):
             :param bool returnP: returns the projection matrices
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"

SimPEG/Mesh/TreeMesh.py

@@ -177,7 +177,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
         return l
 
     def __str__(self):
-        outStr = '  ---- %sTreeMesh ----  '%('Oc' if self.dim == 3 else 'Quad')
+        outStr = '  ---- {0!s}TreeMesh ----  '.format(('Oc' if self.dim == 3 else 'Quad'))
         def printH(hx, outStr=''):
             i = -1
             while True:
@@ -213,7 +213,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
             outStr += printH(self.hy, outStr='\n   hy:')
             outStr += printH(self.hz, outStr='\n   hz:')
         outStr += '\n  nC: {0:d}'.format(self.nC)
-        outStr += '\n  Fill: %2.2f%%'%(self.fill*100)
+        outStr += '\n  Fill: {0:2.2f}%'.format((self.fill*100))
         return outStr
 
     @property
@@ -1875,7 +1875,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
 
         :param numpy.ndarray locs: Location of points to interpolate to
         :param str locType: What to interpolate (see below)
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: M, the interpolation matrix
 
         locType can be::
@@ -2210,7 +2210,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
 
         ax.set_xlabel('y' if normal == 'X' else 'x')
         ax.set_ylabel('y' if normal == 'Z' else 'z')
-        ax.set_title('Slice %d, %s = %4.2f' % (ind,normal,indLoc))
+        ax.set_title('Slice {0:d}, {1!s} = {2:4.2f}'.format(ind, normal, indLoc))
 
         if grid:
             _ = antiNormalInd
@@ -2240,7 +2240,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
             if key < 0 : #Handle negative indices
                 key += len( self )
             if key >= len( self ) :
-                raise IndexError, "The index (%d) is out of range."%key
+                raise IndexError, "The index ({0:d}) is out of range.".format(key)
 
             self._numberCells() # no-op if numbered
             index   = self._i2cc[key]

SimPEG/Mesh/View.py

@@ -171,7 +171,7 @@ class TensorView(object):
                         iz = ix + iy*nX
                         if iz < self.nCz:
                             ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
-                                     '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
+                                     '#{0:.0f}'.format(iz),color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
 
         ax.set_title(vType)
         if showIt: plt.show()
@@ -221,10 +221,10 @@ class TensorView(object):
         vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
 
         # Some user error checking
-        assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
+        assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
         assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.'
-        assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
+        assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
-        assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts)
+        assert normal in normalOpts, "normal must be in ['{0!s}']".format("','".join(normalOpts))
         assert type(grid) is bool, 'grid must be a boolean'
 
         szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
@@ -295,7 +295,7 @@ class TensorView(object):
 
         ax.set_xlabel('y' if normal == 'X' else 'x')
         ax.set_ylabel('y' if normal == 'Z' else 'z')
-        ax.set_title('Slice %d' % ind)
+        ax.set_title('Slice {0:.0f}'.format(ind))
         return out
 
 
@@ -316,11 +316,11 @@ class TensorView(object):
         vTypeOptsV = ['CCv','F','E']
         vTypeOpts = vTypeOptsCC + vTypeOptsV
         if view == 'vec':
-            assert vType in vTypeOptsV, "vType must be in ['%s'] when view='vec'" % "','".join(vTypeOptsV)
+            assert vType in vTypeOptsV, "vType must be in ['{0!s}'] when view='vec'".format("','".join(vTypeOptsV))
-        assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
+        assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
 
         viewOpts = ['real','imag','abs','vec']
-        assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
+        assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
 
 
         if ax is None:

SimPEG/Optimization.py

@@ -121,7 +121,7 @@ class Minimize(object):
     @callback.setter
     def callback(self, value):
         if self.callback is not None:
-            print 'The callback on the %s Optimization was replaced.' % self.__name__
+            print 'The callback on the {0!s} Optimization was replaced.'.format(self.__name__)
         self._callback = value
 
 
@@ -131,7 +131,7 @@ class Minimize(object):
 
         Minimizes the function (evalFunction) starting at the location x0.
 
-        :param def evalFunction: function handle that evaluates: f, g, H = F(x)
+        :param callable evalFunction: function handle that evaluates: f, g, H = F(x)
         :param numpy.ndarray x0: starting location
         :rtype: numpy.ndarray
         :return: x, the last iterate of the optimization algorithm
@@ -372,8 +372,8 @@ class Minimize(object):
             Else, a modifySearchDirectionBreak call is preformed.
 
             :param numpy.ndarray p: searchDirection
-            :rtype: numpy.ndarray,bool
+            :rtype: tuple
-            :return: (xt, passLS)
+            :return: (xt, passLS) numpy.ndarray, bool
         """
         # Projected Armijo linesearch
         self._LS_t = 1
@@ -408,8 +408,8 @@ class Minimize(object):
             evalFunction returns a False indicating the break was not caught.
 
             :param numpy.ndarray p: searchDirection
-            :rtype: numpy.ndarray,bool
+            :rtype: tuple
-            :return: (xt, breakCaught)
+            :return: (xt, breakCaught) numpy.ndarray, bool
         """
         self.printDone(inLS=True)
         print 'The linesearch got broken. Boo.'
@@ -855,7 +855,7 @@ class NewtonRoot(object):
             if self.comments and self.doLS: print '\tLinesearch:\n'
             # Enter Linesearch
             while True and self.doLS:
-                if self.comments: print '\t\tResid: %e\n'%norm(rt)
+                if self.comments: print '\t\tResid: {0:e}\n'.format(norm(rt))
                 if norm(rt) <= norm(r) or norm(rt) < self.tol:
                     break
 
@@ -873,7 +873,7 @@ class NewtonRoot(object):
             if norm(rt) < self.tol:
                 break
             if self.iter > self.maxIter:
-                print 'NewtonRoot stopped by maxIters (%d). norm: %4.4e' % (self.maxIter, norm(rt))
+                print 'NewtonRoot stopped by maxIters ({0:d}). norm: {1:4.4e}'.format(self.maxIter, norm(rt))
                 break
 
         return x

SimPEG/Problem.py

@@ -49,7 +49,7 @@ class BaseProblem(object):
 
     def pair(self, d):
         """Bind a survey to this problem instance using pointers."""
-        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
+        assert isinstance(d, self.surveyPair), "Data object must be an instance of a {0!s} class.".format((self.surveyPair.__name__))
         if d.ispaired:
             raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
         self._survey = d

SimPEG/PropMaps.py

@@ -19,85 +19,85 @@ class Property(object):
         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
+        assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("{0!s}", Maps.ReciprocalMap)'.format(self.name)
         self._propertyLink = value
 
     def _getMapProperty(self):
         prop = self
         def fget(self):
-            return getattr(self, '_%sMap'%prop.name, None)
+            return getattr(self, '_{0!s}Map'.format(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.'
+                assert getattr(self, '{0!s}Map'.format(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)
+            setattr(self, '_{0!s}Map'.format(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))
+            return getattr(self, '_{0!s}Index'.format(prop.name), slice(None))
         def fset(self, val):
-            setattr(self, '_%sIndex'%prop.name, val)
+            setattr(self, '_{0!s}Index'.format(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)
+            mapping = getattr(self, '{0!s}Map'.format(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:
+                if getattr(self, '{0!s}Map'.format(linkName), None) is None:
                     return prop.defaultVal
-                m = getattr(self, '%s'%linkName)
+                m = getattr(self, '{0!s}'.format(linkName))
                 return linkMap * m
 
-            m = getattr(self, '%sModel'%prop.name)
+            m = getattr(self, '{0!s}Model'.format(prop.name))
             return mapping * m
         return property(fget=fget)
 
     def _getModelDerivProperty(self):
         prop = self
         def fget(self):
-            mapping = getattr(self, '%sMap'%prop.name)
+            mapping = getattr(self, '{0!s}Map'.format(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)
+                linkedMap = getattr(self, '{0!s}Map'.format(linkName))
                 if linkedMap is None:
                     return None
                 linkMap = linkMapClass(None) * linkedMap
-                m = getattr(self, '%sModel'%linkName)
+                m = getattr(self, '{0!s}Model'.format(linkName))
                 return linkMap.deriv( m )
 
-            m = getattr(self, '%sModel'%prop.name)
+            m = getattr(self, '{0!s}Model'.format(prop.name))
             return mapping.deriv( m )
         return property(fget=fget)
 
     def _getModelProperty(self):
         prop = self
         def fget(self):
-            mapping = getattr(self, '%sMap'%prop.name)
+            mapping = getattr(self, '{0!s}Map'.format(prop.name))
             if mapping is None:
                 return None
-            index = getattr(self.propMap, '%sIndex'%prop.name)
+            index = getattr(self.propMap, '{0!s}Index'.format(prop.name))
             return self.vector[index]
         return property(fget=fget)
 
     def _getModelProjProperty(self):
         prop = self
         def fget(self):
-            mapping = getattr(self, '%sMap'%prop.name)
+            mapping = getattr(self, '{0!s}Map'.format(prop.name))
             if mapping is None:
                 return None
-            inds = getattr(self.propMap, '%sIndex'%prop.name)
+            inds = getattr(self.propMap, '{0!s}Index'.format(prop.name))
             if type(inds) is slice:
                 inds = range(*inds.indices(self.nP))
             nI, nP = len(inds),self.nP
@@ -107,7 +107,7 @@ class Property(object):
     def _getModelMapProperty(self):
         prop = self
         def fget(self):
-            return getattr(self.propMap, '_%sMap'%prop.name, None)
+            return getattr(self.propMap, '_{0!s}Map'.format(prop.name), None)
         return property(fget=fget)
 
 
@@ -123,7 +123,7 @@ class PropModel(object):
         inds = []
         if getattr(self, '_nP', None) is None:
             for name in self.propMap._properties:
-                index = getattr(self.propMap, '%sIndex'%name, None)
+                index = getattr(self.propMap, '{0!s}Index'.format(name), None)
                 if index is not None:
                     if type(index) is slice:
                         inds += range(*index.indices(len(self.vector)))
@@ -163,9 +163,9 @@ class _PropMapMetaClass(type):
             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])
+                assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '{0!s}' is trying to link to '{1!s}'".format(prop.name, prop.propertyLink[0])
         if len(defaultInvProps) > 1:
-            raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
+            raise Exception('You have more than one default inversion property: {0!s}'.format(defaultInvProps))
 
         newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
 
@@ -187,7 +187,7 @@ class _PropMapMetaClass(type):
             attrs[attr + 'Model'] = prop._getModelProperty()
             attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
 
-        return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
+        return type('PropModel', (PropModel, ), attrs)
 
 
 class PropMap(object):
@@ -223,7 +223,7 @@ class PropMap(object):
                 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]))
+                for m in maps]), "Use signature: [{0!s}]".format((', '.join(["('{0!s}', {1!s}Map)".format(p, p) for p in self._properties])))
         if slices is None:
             slices = dict()
         else:
@@ -236,8 +236,8 @@ class PropMap(object):
 
         nP = 0
         for name, mapping in maps:
-            setattr(self, '%sMap'%name, mapping)
+            setattr(self, '{0!s}Map'.format(name), mapping)
-            setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
+            setattr(self, '{0!s}Index'.format(name), slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
         self.nP = nP
 
@@ -250,12 +250,12 @@ class PropMap(object):
 
     def clearMaps(self):
         for name in self._properties:
-            setattr(self, '%sMap'%name, None)
+            setattr(self, '{0!s}Map'.format(name), None)
-            setattr(self, '%sIndex'%name, None)
+            setattr(self, '{0!s}Index'.format(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]
+        activeMaps = [name for name in self._properties if getattr(self, '{0!s}Map'.format(name)) is not None]
         return val in activeMaps

SimPEG/Regularization.py

@@ -10,7 +10,7 @@ class RegularizationMesh(object):
     are not necessarily true differential operators, but are constructed from
     a SimPEG Mesh.
 
-    :param Mesh mesh: problem mesh
+    :param BaseMesh mesh: problem mesh
     :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
     """
 
@@ -383,8 +383,8 @@ class BaseRegularization(object):
 
         :param numpy.array m: geophysical model
         :param numpy.array v: vector to multiply
-        :rtype: scipy.sparse.csr_matrix or numpy.ndarray
+        :rtype: scipy.sparse.csr_matrix
-        :return: WtW or WtW*v
+        :return: WtW, or if v is supplied WtW*v (numpy.ndarray)
 
         The regularization is:
 
@@ -650,8 +650,8 @@ class Tikhonov(Simple):
     Note if the key word argument `mrefInSmooth` is False, then mref is not
     included in the smoothness contribution.
 
-    :param Mesh mesh: SimPEG mesh
+    :param BaseMesh mesh: SimPEG mesh
-    :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
+    :param IdentityMap mapping: regularization mapping, takes the model from model space to the thing you want to regularize
     :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
     :param bool mrefInSmooth: (default = False) put mref in the smoothness component?
     :param float alpha_s: (default 1e-6) smallness weight
@@ -671,7 +671,7 @@ class Tikhonov(Simple):
     alpha_yy     = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'],    "Weight for the second derivative in the y direction")
     alpha_zz     = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'],    "Weight for the second derivative in the z direction")
 
-    def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property

SimPEG/Survey.py

@@ -26,7 +26,7 @@ class BaseRx(object):
     def rxType(self, value):
         known = self.knownRxTypes
         if known is not None:
-            assert value in known, "rxType must be in ['%s']" % ("', '".join(known))
+            assert value in known, "rxType must be in ['{0!s}']".format(("', '".join(known)))
         self._rxType = value
 
     @property
@@ -125,7 +125,7 @@ class BaseSrc(object):
     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 isinstance(rx, self.rxPair), 'rxList must be a {0!s}'.format(self.rxPair.__name__)
         assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
         self.uid = str(uuid.uuid4())
         self.rxList = rxList
@@ -227,7 +227,7 @@ class BaseSurvey(object):
     @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 np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of {0!s}'.format(self.srcPair.__name__)
         assert len(set(value)) == len(value), 'The srcList must be unique'
         self._srcList = value
         self._sourceOrder = dict()
@@ -238,10 +238,10 @@ class BaseSurvey(object):
             sources = [sources]
         for src in sources:
             if getattr(src,'uid',None) is None:
-                raise KeyError('Source does not have a uid: %s'%str(src))
+                raise KeyError('Source does not have a uid: {0!s}'.format(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))
+            raise KeyError('Some of the sources specified are not in this survey. {0!s}'.format(str(inds)))
         return inds
 
     @property
@@ -263,7 +263,7 @@ class BaseSurvey(object):
     def pair(self, p):
         """Bind a problem to this survey instance using pointers"""
         assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'."
-        assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__)
+        assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a {0!s} class.".format((p.surveyPair.__name__))
         if p.ispaired:
             raise Exception("The problem object is already paired to a survey. Use prob.unpair()")
         self._prob = p
@@ -311,7 +311,6 @@ class BaseSurvey(object):
         if f is None: f = self.prob.fields(m)
         return Utils.mkvc(self.eval(f))
 
-
     @Utils.count
     def eval(self, f):
         """eval(f)
@@ -322,7 +321,7 @@ class BaseSurvey(object):
 
                 d_\\text{pred} = \mathbf{P} f(m)
         """
-        raise NotImplemented('eval is not yet implemented.')
+        raise NotImplementedError('eval is not yet implemented.')
 
     @Utils.count
     def evalDeriv(self, f):
@@ -334,7 +333,7 @@ class BaseSurvey(object):
 
                 \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
         """
-        raise NotImplemented('eval is not yet implemented.')
+        raise NotImplementedError('eval is not yet implemented.')
 
     @Utils.count
     def residual(self, m, f=None):

SimPEG/Tests.py

@@ -4,7 +4,6 @@ from SimPEG.Utils import mkvc, sdiag, diagEst
 from SimPEG import Utils
 from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
 from SimPEG.Mesh.TreeMesh import TreeMesh as Tree
-import numpy as np
 import scipy.sparse as sp
 import unittest
 import inspect
@@ -200,10 +199,10 @@ class OrderTest(unittest.TestCase):
                     print '_____________________________________________'
                     print '   h  |    error    | e(i-1)/e(i) |  order'
                     print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
-                    print '%4i  |  %8.2e   |' % (nc, err)
+                    print '{0:4d}  |  {1:8.2e}   |'.format(nc, err)
                 else:
                     order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
-                    print '%4i  |  %8.2e   |   %6.4f    |  %6.4f' % (nc, err, err_old/err, order[-1])
+                    print '{0:4d}  |  {1:8.2e}   |   {2:6.4f}    |  {3:6.4f}'.format(nc, err, err_old/err, order[-1])
                 err_old = err
                 max_h_old = max_h
             print '---------------------------------------------'
@@ -237,7 +236,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
         Compares error decay of 0th and 1st order Taylor approximation at point
         x0 for a randomized search direction.
 
-        :param lambda fctn: function handle
+        :param callable fctn: function handle
         :param numpy.array x0: point at which to check derivative
         :param int num: number of times to reduce step length, h
         :param bool plotIt: if you would like to plot
@@ -258,8 +257,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
             Tests.checkDerivative(simplePass, np.random.randn(5))
     """
 
-    print "%s checkDerivative %s" % ('='*20, '='*20)
+    print "{0!s} checkDerivative {1!s}".format('='*20, '='*20)
-    print "iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order\n%s" % ('-'*57)
+    print "iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order\n{0!s}".format(('-'*57))
 
     f0, J0 = fctn(x0)
 
@@ -290,7 +289,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
 
         order0 = np.log10(E0[:-1]/E0[1:])
         order1 = np.log10(E1[:-1]/E1[1:])
-        print " %d   %1.2e    %1.3e     %1.3e      %1.3f" % (i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
+        print " {0:d}   {1:1.2e}    {2:1.3e}     {3:1.3e}      {4:1.3f}".format(i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
 
     # Ensure we are about precision
     order0 = order0[E0[1:] > eps]
@@ -302,10 +301,10 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
     passTest = belowTol or correctOrder
 
     if passTest:
-        print "%s PASS! %s" % ('='*25, '='*25)
+        print "{0!s} PASS! {1!s}".format('='*25, '='*25)
         print happiness[np.random.randint(len(happiness))]+'\n'
     else:
-        print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57)
+        print "{0!s}\n{1!s} FAIL! {2!s}\n{3!s}".format('*'*57, '<'*25, '>'*25, '*'*57)
         print sadness[np.random.randint(len(sadness))]+'\n'
 
 
@@ -314,7 +313,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
         ax = ax or plt.subplot(111)
         ax.loglog(h, E0, 'b')
         ax.loglog(h, E1, 'g--')
-        ax.set_title('Check Derivative - %s' % ('PASSED :)' if passTest else 'FAILED :('))
+        ax.set_title('Check Derivative - {0!s}'.format(('PASSED :)' if passTest else 'FAILED :(')))
         ax.set_xlabel('h')
         ax.set_ylabel('Error')
         leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best',

SimPEG/Utils/ModelBuilder.py

@@ -7,11 +7,11 @@ def addBlock(gridCC, modelCC, p0, p1, blockProp):
     """
         Add a block to an exsisting cell centered model, modelCC
 
-        :param numpy.array, gridCC: mesh.gridCC is the cell centered grid
+        :param numpy.array gridCC: mesh.gridCC is the cell centered grid
-        :param numpy.array, modelCC: cell centered model
+        :param numpy.array modelCC: cell centered model
-        :param numpy.array, p0: bottom, southwest corner of block
+        :param numpy.array p0: bottom, southwest corner of block
-        :param numpy.array, p1: top, northeast corner of block
+        :param numpy.array p1: top, northeast corner of block
-        :blockProp float, blockProp: property to assign to the model 
+        :blockProp float blockProp: property to assign to the model
 
         :return numpy.array, modelBlock: model with block
     """
@@ -147,7 +147,7 @@ def getIndicesSphere(center,radius,ccMesh):
 
     if dimMesh == 1:
        # Define the reference points
-        
+
         ind  = np.abs(center[0] - ccMesh[:,0]) < radius
 
     elif dimMesh == 2:
@@ -222,14 +222,14 @@ def layeredModel(ccMesh, layerTops, layerValues):
 
         :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) 
+        :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 
+        :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 
+    # 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
@@ -253,10 +253,10 @@ def layeredModel(ccMesh, layerTops, layerValues):
     model = np.zeros(ccMesh.shape[0])
 
     for i, top in enumerate(layerTops):
-        zind = z <= top 
+        zind = z <= top
         model[zind] = layerValues[i]
 
-    return model 
+    return model
 
 
 
@@ -265,9 +265,9 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
         Create a random model by convolving a kernel with a
         uniformly distributed model.
 
-        :param int,tuple shape: shape of the model.
+        :param 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 kernel that is used.
+        :param numpy.ndarray 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/SolverUtils.py

@@ -8,12 +8,12 @@ def _checkAccuracy(A, b, X, accuracyTol):
     if nrm_b > 0:
         nrm /= nrm_b
     if nrm > accuracyTol:
-        msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
+        msg = '### SolverWarning ###: Accuracy on solve is above tolerance: {0:e} > {1:e}'.format(nrm, accuracyTol)
         print msg
         warnings.warn(msg, RuntimeWarning)
 
 
-def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
+def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6, name=None):
     """
     Wraps a direct Solver.
 
@@ -72,11 +72,11 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
         if factorize and hasattr(self.solver, 'clean'):
             return self.solver.clean()
 
-    return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
+    return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
 
 
-def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
+def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5, name=None):
     """
     Wraps an iterative Solver.
 
@@ -128,13 +128,13 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
     def clean(self):
         pass
 
-    return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
+    return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
 
 from scipy.sparse import linalg
-Solver   = SolverWrapD(linalg.spsolve, factorize=False)
+Solver   = SolverWrapD(linalg.spsolve, factorize=False, name="Solver")
-SolverLU = SolverWrapD(linalg.splu, factorize=True)
+SolverLU = SolverWrapD(linalg.splu, factorize=True, name="SolverLU")
-SolverCG = SolverWrapI(linalg.cg)
+SolverCG = SolverWrapI(linalg.cg, name="SolverCG")
 
 
 class SolverDiag(object):

SimPEG/Utils/codeutils.py

@@ -32,7 +32,7 @@ def memProfileWrapper(towrap, *funNames):
         if hasattr(towrap,f):
             attrs[f] = profile(getattr(towrap,f))
         else:
-            print '%s not found in %s Class' % (f, towrap.__name__)
+            print '{0!s} not found in {1!s} Class'.format(f, towrap.__name__)
 
     return type(towrap.__name__ + 'MemProfileWrap', (towrap,), attrs)
 
@@ -65,7 +65,7 @@ def setKwargs(obj, ignore=None,  **kwargs):
         if hasattr(obj, attr):
             setattr(obj, attr, kwargs[attr])
         else:
-            raise Exception('%s attr is not recognized' % attr)
+            raise Exception('{0!s} attr is not recognized'.format(attr))
 
     hook(obj,hook, silent=True)
     hook(obj,setKwargs, silent=True)
@@ -74,7 +74,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
     titles = ''
     widths = 0
     for printer in printers:
-        titles += ('{:^%i}'%printer['width']).format(printer['title']) + ''
+        titles += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['title']) + ''
         widths += printer['width']
     print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
     print pad + titles
@@ -83,7 +83,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
 def printLine(obj, printers, pad=''):
     values = ''
     for printer in printers:
-        values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))
+        values += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['format'] % printer['value'](obj))
     print pad + values
 
 def checkStoppers(obj, stoppers):
@@ -104,12 +104,12 @@ def checkStoppers(obj, stoppers):
     return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
 
 def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
-    print pad + "%s%s%s" % ('-'*25,stop,'-'*25)
+    print pad + "{0!s}{1!s}{2!s}".format('-'*25, stop, '-'*25)
     for stopper in stoppers:
         l = stopper['left'](obj)
         r = stopper['right'](obj)
         print pad + stopper['str'] % (l<=r,l,r)
-    print pad + "%s%s%s" % ('-'*25,done,'-'*25)
+    print pad + "{0!s}{1!s}{2!s}".format('-'*25, done, '-'*25)
 
 def callHooks(match, mainFirst=False):
     """
@@ -144,14 +144,14 @@ def callHooks(match, mainFirst=False):
 
 
         extra = """
-            If you have things that also need to run in the method %s, you can create a method::
+            If you have things that also need to run in the method {0!s}, you can create a method::
 
-                def _%s*(self, ... ):
+                def _{1!s}*(self, ... ):
                     pass
 
-            Where the * can be any string. If present, _%s* will be called at the start of the default %s call.
+            Where the * can be any string. If present, _{2!s}* will be called at the start of the default {3!s} call.
             You may also completely overwrite this function.
-        """ % (match, match, match, match)
+        """.format(match, match, match, match)
         doc = wrapper.__doc__
         wrapper.__doc__ = ('' if doc is None else doc) + extra
         return wrapper
@@ -186,7 +186,7 @@ def asArray_N_x_Dim(pts, dim):
         elif len(pts.shape) == 1:
             pts = pts[:,np.newaxis]
 
-        assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, %d) not (%d, %d)" % ((dim,)+pts.shape)
+        assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, {0:d}) not ({1:d}, {2:d})".format(*((dim,)+pts.shape))
 
         return pts
 
@@ -207,17 +207,17 @@ def requires(var):
 
         .. note::
 
-            To use survey.%s(), SimPEG requires that a problem be bound to the survey.
+            To use survey.{0!s}(), SimPEG requires that a problem be bound to the survey.
             If a problem has not been bound, an Exception will be raised.
             To bind a problem to the Data object::
 
                 survey.pair(myProblem)
 
-            """ % f.__name__
+            """.format(f.__name__)
         else:
             extra = """
-                To use *%s* method, SimPEG requires that the %s be specified.
+                To use *{0!s}* method, SimPEG requires that the {1!s} be specified.
-            """ % (f.__name__, var)
+            """.format(f.__name__, var)
         @wraps(f)
         def requiresVarWrapper(self,*args,**kwargs):
             if getattr(self, var, None) is None:

SimPEG/Utils/curvutils.py

@@ -80,7 +80,7 @@ def indexCube(nodes, gridSize, n=None):
     # Make sure that we choose from the possible nodes.
     possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
     for node in nodes:
-        assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
+        assert node in possibleNodes, "Nodes must be chosen from: '{0!s}'".format(possibleNodes)
     dim = gridSize.size
     if n is None:
         n = gridSize - 1

SimPEG/Utils/interputils.py

@@ -25,7 +25,7 @@ def interpmat(locs, x, y=None, z=None):
         :param numpy.ndarray x: Tensor vector of 1st dimension of grid.
         :param numpy.ndarray y: Tensor vector of 2nd dimension of grid. None by default.
         :param numpy.ndarray z: Tensor vector of 3rd dimension of grid. None by default.
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: Interpolation matrix
 
         .. plot::

SimPEG/Utils/matutils.py

@@ -27,7 +27,7 @@ def mkvc(x, numDims=1):
 
     if isinstance(x, Zero):
         return x
-    
+
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -278,7 +278,7 @@ class TensorType(object):
         else:
             raise Exception('Unexpected shape of tensor')
     def __str__(self):
-        return 'TensorType[%i]: %s' % (self._tt, self._tts)
+        return 'TensorType[{0:d}]: {1!s}'.format(self._tt, self._tts)
     def __eq__(self, v): return self._tt == v
     def __le__(self, v): return self._tt <= v
     def __ge__(self, v): return self._tt >= v
@@ -355,9 +355,9 @@ def diagEst(matFun, n, k=None, approach='Probing'):
         2. Ones : random +/- 1 entries
         3. Random : random vectors
 
-    :param lambda (numpy.array) matFun: matrix to estimate the diagonal of
+    :param callable matFun: takes a (numpy.array) and multiplies it by a matrix to estimate the diagonal
-    :param int64 n: size of the vector that should be used to compute matFun(v)
+    :param int 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 int 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)
@@ -422,9 +422,9 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
-    @property 
+    @property
     def transpose(self): return Zero()
-    
+
     @property
     def T(self): return Zero()
 

SimPEG/Utils/meshutils.py

@@ -83,7 +83,7 @@ def closestPoints(mesh, pts, gridLoc='CC'):
     """
         Move a list of points to the closest points on a grid.
 
-        :param simpeg.Mesh.BaseMesh mesh: The mesh
+        :param BaseMesh mesh: The mesh
         :param numpy.ndarray pts: Points to move
         :param string gridLoc: ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ex', 'Ey', 'Ez']
         :rtype: numpy.ndarray
@@ -104,16 +104,20 @@ def closestPoints(mesh, pts, gridLoc='CC'):
 
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
     """
-        Extracts Core Mesh from Global mesh
+    Extracts Core Mesh from Global mesh
-        xyzlim: 2D array [ndim x 2]
+
-        mesh: SimPEG mesh
+    :param numpy.ndarray xyzlim: 2D array [ndim x 2]
-        This function ouputs:
+    :param BaseMesh mesh: The mesh
-            - actind: corresponding boolean index from global to core
+
-            - meshcore: core SimPEG mesh
+    This function ouputs::
-        Warning: 1D and 2D has not been tested
+
+        - 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:
+    if mesh.dim == 1:
         xyzlim = xyzlim.flatten()
         xmin, xmax = xyzlim[0], xyzlim[1]
 
@@ -125,11 +129,11 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
 
-        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        meshCore = Mesh.TensorMesh([hx, hy], x0=x0)
 
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
 
-    elif mesh.dim ==2:
+    elif mesh.dim == 2:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
 
@@ -144,12 +148,12 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
 
-        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        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:
+    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]
@@ -168,7 +172,7 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         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)
+        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) \

SimPEG/Utils/modelutils.py

@@ -26,7 +26,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
             gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F')
 
             if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
-                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
 
             Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes
             actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
@@ -47,7 +47,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
 
             gridTopo = Ftopo(mesh.vectorNx)
             if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
-                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
 
             Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes
             actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')

SimPEG/__init__.py

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

docs/Makefile

@@ -2,7 +2,7 @@
 #
 
 # You can set these variables from the command line.
-SPHINXOPTS    =
+SPHINXOPTS    = -n -w warnings.txt
 SPHINXBUILD   = sphinx-build
 PAPER         =
 BUILDDIR      = _build

docs/_templates/layout.html

@@ -0,0 +1,22 @@
+{# Import the theme's layout. #}
+{% extends "!layout.html" %}
+
+{% block extrahead %}
+{{ super() }}
+
+<meta name="description" content="Simulation and Parameter Estimation in Geophysics">
+<meta name="author" content="SimPEG Developers">
+<meta name="keywords" content="python, geophysics, inversion, electromagnetics, magnetotellurics, magnetics, gravity, DC, flow inverse problems, open source, finite volume">
+
+
+<script>
+  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
+  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
+  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
+  })(window,document,'script','https://www.google-analytics.com/analytics.js','ga');
+
+  ga('create', 'UA-45185336-1', 'auto');
+  ga('send', 'pageview');
+
+</script>
+{% endblock %}

docs/api_FiniteVolume.rst

@@ -1,19 +0,0 @@
-.. _api_FiniteVolume:
-
-Finite Volume
-*************
-
-Any numerical implementation requires the discretization of continuous functions into discrete approximations. These approximations are typically organized in a mesh, which defines boundaries, locations, and connectivity. Of specific interest to geophysical simulations, we require that averaging, interpolation and differential operators be defined for any mesh. In SimPEG, we have implemented a staggered mimetic finite volume approach (`Hyman and Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This approach requires the definitions of variables at either cell-centers, nodes, faces, or edges as seen in the figure below.
-
-.. image:: images/finitevolrealestate.png
-   :width: 400 px
-   :alt: FiniteVolume
-   :align: center
-
-
-.. toctree::
-   :maxdepth: 2
-
-   api_Mesh
-   api_DiffOps
-   api_InnerProducts

docs/api_MeshCode.rst

@@ -1,36 +0,0 @@
-.. _api_MeshCode:
-
-Tensor Mesh
-===========
-
-.. automodule:: SimPEG.Mesh.TensorMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Cylindrical Mesh
-================
-
-.. automodule:: SimPEG.Mesh.CylMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Tree Mesh
-=========
-
-.. autoclass:: SimPEG.Mesh.TreeMesh.TreeMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Curvilinear Mesh
-================
-
-.. automodule:: SimPEG.Mesh.CurvilinearMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:

docs/app.yaml

@@ -0,0 +1,95 @@
+# application: simpegdocs
+# version: 1
+runtime: python27
+api_version: 1
+threadsafe: yes
+
+handlers:
+
+# favicon
+- url: /images/logo-block\.ico
+  static_files: /images/logo-block.ico
+  upload: /images/logo-block\.ico
+
+# all css
+- url: /(.*\.css)
+  mime_type: text/css
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.css)
+
+# webfonts
+- url: /(.*\.(eot|svg|ttf|woff|woff2|otf))
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.(eot|svg|ttf|woff|woff2|otf))
+
+# javascript
+- url: /(.*\.js)
+  mime_type: text/javascript
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.js)
+
+# plain text source
+- url: /(.*\.txt)
+  mime_type: text/plain
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.txt)
+
+# images
+- url: /_images/(.*\.(gif|png|jpg|ico))
+  static_files: _build/html/_images/\1
+  upload: _build/html/_images/(.*\.(gif|png|jpg|ico))
+
+# redirect en/latest traffic
+- url: /en/latest/(.*\.html)
+  script: simpegdocs.app
+
+# raw html
+- url: /(.*\.html)
+  mime_type: text/html
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.html)
+
+# serve index files
+- url: /(.+)/
+  static_files: _build/html/\1/index.html
+  upload: _build/html/(.+)/index.html
+
+- url: /(.+)
+  static_files: _build/html/\1/index.html
+  upload: _build/html/(.+)/index.html
+
+- url: /
+  static_files: _build/html/index.html
+  upload: _build/html/index.html
+
+- url: .*
+  script: simpegdocs.app
+
+# Recommended file skipping declaration from the GAE tutorials
+skip_files:
+  - ^(.*/)?app\.yaml
+  - ^(.*/)?app\.yml
+  - ^(.*/)?#.*#
+  - ^(.*/)?.*~
+  - ^(.*/)?.*\.py[co]
+  - ^(.*/)?.*/RCS/.*
+  - ^(.*/)?\..*
+  - ^(.*/)?tests$
+  - ^(.*/)?test$
+  - ^test/(.*/)?
+  - ^COPYING.LESSER
+  - ^README\..*
+  - \.gitignore
+  - ^\.git/.*
+  - \.*\.lint$
+  - ^(.*/)?.*\.doctree$
+
+libraries:
+- name: webapp2
+  version: "2.5.2"
+- name: PIL
+  version: "1.1.7"
+- name: numpy
+  version: "latest"
+- name: jinja2
+  version: "latest"

docs/conf.py

@@ -28,7 +28,7 @@ sys.path.append('../')
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'matplotlib.sphinxext.plot_directive']
+extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'matplotlib.sphinxext.plot_directive']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -44,16 +44,16 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'SimPEG'
-copyright = u'2013, SimPEG Developers'
+copyright = u'2013 - 2016, SimPEG Developers'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.10'
+version = '0.1.12'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.10'
+release = '0.1.12'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -124,12 +124,12 @@ except Exception, e:
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
 # pixels large.
-#html_favicon = None
+html_favicon = './images/logo-block.ico'
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+html_static_path = []
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
@@ -229,6 +229,12 @@ man_pages = [
 # If true, show URL addresses after external links.
 #man_show_urls = False
 
+# Intersphinx
+intersphinx_mapping = {'python': ('http://docs.python.org/2', None),
+                       'numpy': ('http://docs.scipy.org/doc/numpy/', None),
+                       'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None),
+                       'matplotlib': ('http://matplotlib.sourceforge.net/', None)}
+
 
 # -- Options for Texinfo output ------------------------------------------------
 
@@ -251,3 +257,36 @@ texinfo_documents = [
 #texinfo_show_urls = 'footnote'
 
 autodoc_member_order = 'bysource'
+
+def supress_nonlocal_image_warn():
+    import sphinx.environment
+    sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn
+
+def _supress_nonlocal_image_warn(self, msg, node):
+    from docutils.utils import get_source_line
+
+    if not msg.startswith('nonlocal image URI found:'):
+        self._warnfunc(msg, '{0!s}:{1!s}'.format(*get_source_line(node)))
+
+supress_nonlocal_image_warn()
+
+
+nitpick_ignore = [
+    ('py:class', 'IdentityMap'),
+    ('py:class', 'BaseSurvey'),
+    ('py:class', 'BaseSrc'),
+    ('py:class', 'BaseRx'),
+    ('py:class', 'Survey'),
+    ('py:class', 'FieldsFDEM'),
+    ('py:class', 'Fields3D_e'),
+    ('py:class', 'Fields3D_b'),
+    ('py:class', 'Fields3D_j'),
+    ('py:class', 'Fields3D_h'),
+    ('py:class', 'SurveyTDEM'),
+    ('py:class', 'SrcTDEM'),
+    ('py:class', 'EMPropMap'),
+    ('py:class', 'Data'),
+    ('py:class', 'SurveyDC'),
+    ('py:class', 'BaseMTFields'),
+    ('py:class', 'SolverLU'),
+]

docs/content/api_core/api_Examples.rst

@@ -7,7 +7,7 @@ Examples
     :maxdepth: 1
     :glob:
 
-    examples/*
+    ../examples/*
 
 
 External Notebooks

docs/content/api_core/api_FiniteVolume.rst

@@ -0,0 +1,27 @@
+.. _api_FiniteVolume:
+
+Finite Volume
+*************
+
+Any numerical implementation requires the discretization of continuous
+functions into discrete approximations. These approximations are typically
+organized in a mesh, which defines boundaries, locations, and connectivity. Of
+specific interest to geophysical simulations, we require that averaging,
+interpolation and differential operators be defined for any mesh. In SimPEG,
+we have implemented a staggered mimetic finite volume approach (`Hyman and
+Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This
+approach requires the definitions of variables at either cell-centers, nodes,
+faces, or edges as seen in the figure below.
+
+.. image:: ../../images/finitevolrealestate.png
+   :width: 400 px
+   :alt: FiniteVolume
+   :align: center
+
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Mesh
+   api_DiffOps
+   api_InnerProducts

docs/content/api_core/api_ForwardProblem.rst

@@ -52,13 +52,15 @@ We can take the derivative of the PDE:
 
     \nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0
 
-If the forward problem is invertible, then we can rearrange for \\(\\frac{\\partial u}{\\partial m}\\):
+If the forward problem is invertible, then we can rearrange for
+\\(\\frac{\\partial u}{\\partial m}\\):
 
 .. math::
 
     J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)
 
-This can often be computed given a vector (i.e. \\(J(v)\\)) rather than stored, as \\(J\\) is a large dense matrix.
+This can often be computed given a vector (i.e. \\(J(v)\\)) rather than
+stored, as \\(J\\) is a large dense matrix.
 
 
 
@@ -67,13 +69,45 @@ The API
 
 Problem
 -------
-.. automodule:: SimPEG.Problem
+
+.. autoclass:: SimPEG.Problem.BaseProblem
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Problem.BaseTimeProblem
+    :members:
+    :undoc-members:
+
+Fields
+------
+
+.. autoclass:: SimPEG.Fields.Fields
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Fields.TimeFields
     :members:
     :undoc-members:
 
 Survey
 ------
-.. automodule:: SimPEG.Survey
+
+.. autoclass:: SimPEG.Survey.BaseSurvey
     :members:
     :undoc-members:
 
+.. autoclass:: SimPEG.Survey.BaseSrc
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.BaseRx
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.BaseTimeRx
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.Data
+    :members:
+    :undoc-members:

docs/content/api_core/api_InnerProducts.rst

@@ -4,7 +4,10 @@
 Inner Products
 **************
 
-By using the weak formulation of many of the PDEs in geophysical applications, we can rapidly develop discretizations. Much of this work, however, needs a good understanding of how to approximate inner products on our discretized meshes. We will define the inner product as:
+By using the weak formulation of many of the PDEs in geophysical applications,
+we can rapidly develop discretizations. Much of this work, however, needs a
+good understanding of how to approximate inner products on our discretized
+meshes. We will define the inner product as:
 
 .. math::
 
@@ -14,12 +17,15 @@ where a and b are either scalars or vectors.
 
 .. note::
 
-    The InnerProducts class is a base class providing inner product matrices for meshes and cannot run on its own.
+    The InnerProducts class is a base class providing inner product matrices
+    for meshes and cannot run on its own.
 
 
 Example problem for DC resistivity
 ----------------------------------
-We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics.
+
+We will start with the formulation of the Direct Current (DC) resistivity
+problem in geophysics.
 
 
 .. math::
@@ -28,12 +34,13 @@ We will start with the formulation of the Direct Current (DC) resistivity proble
 
         \nabla\cdot \vec{j} = q
 
-In the following discretization, \\\( \\sigma \\\) and \\\( \\phi \\\)
+In the following discretization, :math:`\sigma` and :math:`\phi`
-will be discretized on the cell-centers and the flux, \\\(\\vec{j}\\\),
+will be discretized on the cell-centers and the flux, :math:`\vec{j}`,
 will be on the faces. We will use the weak formulation to discretize
 the DC resistivity equation.
 
-We can define in weak form by integrating with a general face function \\\(\\vec{f}\\\):
+We can define in weak form by integrating with a general face function
+:math:`\vec{f}`:
 
 .. math::
 
@@ -61,9 +68,16 @@ We can then discretize for every cell:
 
 .. note::
 
-    We have discretized the dot product above, but remember that we do not really have a single vector \\\(\\mathbf{J}\\\), but approximations of \\\(\\vec{j}\\\) on each face of our cell. In 2D that means 2 approximations of \\\(\\mathbf{J}_x\\\) and 2 approximations of \\\(\\mathbf{J}_y\\\). In 3D we also have 2 approximations of \\\(\\mathbf{J}_z\\\).
+    We have discretized the dot product above, but remember that we do not
+    really have a single vector :math:`\mathbf{J}`, but approximations of
+    :math:`\vec{j}` on each face of our cell. In 2D that means 2
+    approximations of :math:`\mathbf{J}_x` and 2 approximations of
+    :math:`\mathbf{J}_y`. In 3D we also have 2 approximations of
+    :math:`\mathbf{J}_z`.
 
-Regardless of how we choose to approximate this dot product, we can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma.
+Regardless of how we choose to approximate this dot product, we can represent
+this in vector form (again this is for every cell), and will generalize for
+the case of anisotropic (tensor) sigma.
 
 .. math::
 
@@ -71,14 +85,17 @@ Regardless of how we choose to approximate this dot product, we can represent th
     -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F})
     + \text{BC}
 
-We multiply by  square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh:
+We multiply by  square-root of volume on each side of the tensor conductivity
+to keep symmetry in the system. Here :math:`\mathbf{J}_c` is the Cartesian
+:math:`\mathbf{J}` (on the faces that we choose to use in our approximation)
+and must be calculated differently depending on the mesh:
 
 .. math::
     \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
     \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.
+Here the :math:`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:
+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 :math:`\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:
 
 .. math::
 
@@ -107,10 +124,12 @@ By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in
         \sum_{i=1}^{2^d}
         \mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)}
 
-Where \\\(d\\\) is the dimension of the mesh.
+Where :math:`d` is the dimension of the mesh.
-The \\\( \\mathbf{M}^f \\\) is returned when given the input of \\\( \\Sigma^{-1} \\\).
+The :math:`\mathbf{M}^f` is returned when given the input of :math:`\Sigma^{-1}`.
 
-Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined):
+Here each :math:`\mathbf{P} ~ \in ~ \mathbb{R}^{(d*nC, nF)}` is a combination
+of the projection, volume, and any normalization to Cartesian coordinates
+(where the dot product is well defined):
 
 .. math::
 
@@ -129,7 +148,10 @@ If ``returnP=True`` is requested in any of these methods the projection matrices
     # In 1D
     P = [P0, P1]
 
-The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same.
+The derivation for ``edgeInnerProducts`` is exactly the same, however, when we
+approximate the integral using the fields around each node, the projection
+matrices look a bit different because we have 12 edges in 3D instead of just 6
+faces. The interface to the code is exactly the same.
 
 
 Defining Tensor Properties
@@ -137,7 +159,8 @@ Defining Tensor Properties
 
 **For 3D:**
 
-Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
+Depending on the number of columns (either 1, 3, or 6) of mu, the material
+property is interpreted as follows:
 
 .. math::
 
@@ -188,13 +211,16 @@ Which is nice and easy to invert if necessary, however, in the fully anisotropic
 Taking Derivatives
 ------------------
 
-We will take the derivative of the fully anisotropic tensor for a 3D mesh, the other cases are easier and will not be discussed here. Let us start with one part of the sum which makes up \\\(\\mathbf{M}^f_\\Sigma\\\) and take the derivative when this is multiplied by some vector \\\(\\mathbf{v}\\\):
+We will take the derivative of the fully anisotropic tensor for a 3D mesh, the
+other cases are easier and will not be discussed here. Let us start with one
+part of the sum which makes up :math:`\mathbf{M}^f_\Sigma` and take the
+derivative when this is multiplied by some vector :math:`\mathbf{v}`:
 
 .. math::
 
     \mathbf{P}^\top \boldsymbol{\Sigma} \mathbf{Pv}
 
-Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) will have the form:
+Here we will let :math:`\mathbf{Pv} = \mathbf{y}` and :math:`\mathbf{y}` will have the form:
 
 .. math::
 
@@ -233,7 +259,9 @@ Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) wi
         \end{matrix}
     \right]
 
-Now it is easy to take the derivative with respect to any one of the parameters, for example, \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_1}\\\)
+Now it is easy to take the derivative with respect to any one of the
+parameters, for example,
+:math:`\frac{\partial}{\partial\boldsymbol{\sigma}_1}`
 
 .. math::
     \frac{\partial}{\partial \boldsymbol{\sigma}_1}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -247,7 +275,8 @@ Now it is easy to take the derivative with respect to any one of the parameters,
         \end{matrix}
     \right]
 
-Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example, is:
+Whereas :math:`\frac{\partial}{\partial\boldsymbol{\sigma}_4}`, for
+example, is:
 
 .. math::
     \frac{\partial}{\partial \boldsymbol{\sigma}_4}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -261,11 +290,12 @@ Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example
         \end{matrix}
     \right]
 
-These are computed for each of the 8 projections, horizontally concatenated, and returned.
+These are computed for each of the 8 projections, horizontally concatenated,
+and returned.
 
 The API
 -------
 
-.. automodule:: SimPEG.Mesh.InnerProducts
+.. autoclass:: SimPEG.Mesh.InnerProducts.InnerProducts
     :members:
     :undoc-members:

docs/content/api_core/api_Inversion.rst

@@ -3,7 +3,7 @@
 InvProblem
 **********
 
-.. automodule:: SimPEG.InvProblem
+.. autoclass:: SimPEG.InvProblem.BaseInvProblem
     :show-inheritance:
     :members:
     :undoc-members:
@@ -12,7 +12,7 @@ InvProblem
 Inversion
 *********
 
-.. automodule:: SimPEG.Inversion
+.. autoclass:: SimPEG.Inversion.BaseInversion
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/api_core/api_Maps.rst

@@ -27,7 +27,8 @@ back to conductivity. This is a relatively trivial example (we are just taking
 the exponential!) but by defining maps we can start to combine and manipulate
 exactly what we think about as our model, \\\(m\\\). In code, this looks like
 
-::
+.. code-block:: python
+    :linenos:
 
     M = Mesh.TensorMesh([100]) # Create a mesh
     expMap = Maps.ExpMap(M)    # Create a mapping
@@ -46,14 +47,15 @@ We will use an example where we want a 1D layered earth as
 our model, but we want to map this to a 2D discretization to do our forward
 modeling. We will also assume that we are working in log conductivity still,
 so after the transformation we want to map to conductivity space.
-To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.Vertical1DMap`),
+To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
 which does the first part of what we just described. The second part will be
 done by the :class:`SimPEG.Maps.ExpMap` described above.
 
-::
+.. code-block:: python
+    :linenos:
 
     M = Mesh.TensorMesh([7,5])
-    v1dMap = Maps.Vertical1DMap(M)
+    v1dMap = Maps.SurjectVertical1D(M)
     expMap = Maps.ExpMap(M)
     myMap = expMap * v1dMap
     m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
@@ -61,26 +63,8 @@ done by the :class:`SimPEG.Maps.ExpMap` described above.
 
 .. plot::
 
-    from SimPEG import *
+    from SimPEG import Examples
-    import matplotlib.pyplot as plt
+    Examples.Maps_ComboMaps.run()
-    M = Mesh.TensorMesh([7,5])
-    v1dMap = Maps.Vertical1DMap(M)
-    expMap = Maps.ExpMap(M)
-    myMap = expMap * v1dMap
-    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
-    sig = myMap * m
-    figs, axs = plt.subplots(1,2)
-    axs[0].plot(m, M.vectorCCy, 'b-o')
-    axs[0].set_title('Model')
-    axs[0].set_ylabel('Depth, y')
-    axs[0].set_xlabel('Value, $m_i$')
-    axs[0].set_xlim(0,3)
-    axs[0].set_ylim(0,1)
-    clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
-    axs[1].set_title('Physical Property')
-    axs[1].set_ylabel('Depth, y')
-    clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
-    plt.tight_layout()
 
 If you noticed, it was pretty easy to combine maps. What is even cooler is
 that the derivatives also are made for you (if everything goes right).
@@ -122,6 +106,8 @@ When these are used in the inverse problem, this is extremely important!!
 The API
 =======
 
+The :code:`IdentityMap` is the base class for all mappings, and it does absolutely nothing.
+
 .. autoclass:: SimPEG.Maps.IdentityMap
     :members:
     :undoc-members:
@@ -130,7 +116,6 @@ The API
 Common Maps
 ===========
 
-
 Exponential Map
 ---------------
 
@@ -148,7 +133,7 @@ lives (i.e. it varies logarithmically).
 Vertical 1D Map
 ---------------
 
-.. autoclass:: SimPEG.Maps.Vertical1DMap
+.. autoclass:: SimPEG.Maps.SurjectVertical1D
     :members:
     :undoc-members:
 
@@ -164,31 +149,10 @@ Map 2D Cross-Section to 3D Model
 Mesh to Mesh Map
 ----------------
 
-
 .. plot::
 
-    from SimPEG import *
+    from SimPEG import Examples
-    import matplotlib.pyplot as plt
+    Examples.Maps_Mesh2Mesh.run()
-    M = Mesh.TensorMesh([100,100])
-    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
-    h1 = h1/h1.sum()
-    M2 = Mesh.TensorMesh([h1,h1])
-    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
-    v = Utils.mkvc(V)
-    modh = Maps.Mesh2Mesh([M,M2])
-    modH = Maps.Mesh2Mesh([M2,M])
-    H = modH * v
-    h = modh * H
-    ax = plt.subplot(131)
-    M.plotImage(v, ax=ax)
-    ax.set_title('Fine Mesh (Original)')
-    ax = plt.subplot(132)
-    M2.plotImage(H,clim=[0,1],ax=ax)
-    ax.set_title('Course Mesh')
-    ax = plt.subplot(133)
-    M.plotImage(h,clim=[0,1],ax=ax)
-    ax.set_title('Fine Mesh (Interpolated)')
-    plt.show()
 
 
 .. autoclass:: SimPEG.Maps.Mesh2Mesh
@@ -196,8 +160,8 @@ Mesh to Mesh Map
     :undoc-members:
 
 
-Some Extras
+Under the Hood
-===========
+==============
 
 Combo Map
 ---------

docs/content/api_core/api_Mesh.rst

@@ -188,6 +188,6 @@ other types of meshes in this SimPEG framework.
 The API
 =======
 
-.. automodule:: SimPEG.Mesh.BaseMesh
+.. autoclass:: SimPEG.Mesh.BaseMesh.BaseMesh
     :members:
     :undoc-members:

docs/content/api_core/api_MeshCode.rst

@@ -0,0 +1,68 @@
+.. _api_MeshCode:
+
+Tensor Mesh
+===========
+
+.. autoclass:: SimPEG.Mesh.TensorMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Cylindrical Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.CylMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Tree Mesh
+=========
+
+.. autoclass:: SimPEG.Mesh.TreeMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Curvilinear Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.CurvilinearMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Base Rectangular Mesh
+=====================
+
+.. autoclass:: SimPEG.Mesh.BaseMesh.BaseRectangularMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Base Tensor Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.TensorMesh.BaseTensorMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Mesh IO
+=======
+
+.. automodule:: SimPEG.Mesh.MeshIO
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Mesh Viewing
+============
+
+.. automodule:: SimPEG.Mesh.View
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/content/api_core/api_PropMaps.rst

@@ -0,0 +1,29 @@
+SimPEG PropMaps
+***************
+
+The API
+=======
+
+Property
+--------
+
+.. autoclass:: SimPEG.PropMaps.Property
+    :members:
+    :undoc-members:
+
+
+PropMap
+-------
+
+.. autoclass:: SimPEG.PropMaps.PropMap
+    :members:
+    :undoc-members:
+
+
+PropModel
+---------
+
+.. autoclass:: SimPEG.PropMaps.PropModel
+    :members:
+    :undoc-members:
+

docs/content/api_core/api_Regularization.rst

@@ -91,10 +91,21 @@ The API
     :members:
     :undoc-members:
 
+.. autoclass:: SimPEG.Regularization.Simple
+    :show-inheritance:
+    :members:
 
 .. autoclass:: SimPEG.Regularization.Tikhonov
     :show-inheritance:
     :members:
 
+.. autoclass:: SimPEG.Regularization.Sparse
+    :show-inheritance:
+    :members:
+
+.. autoclass:: SimPEG.Regularization.RegularizationMesh
+    :show-inheritance:
+    :members:
+
 
 

docs/content/api_core/api_Solver.rst

@@ -46,6 +46,8 @@ The API
 =======
 
 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
+    :noindex:
 
 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI
+    :noindex:
 

docs/content/api_core/api_Utilities.rst

@@ -6,5 +6,6 @@ Utilities
 
    api_Solver
    api_Maps
+   api_PropMaps
    api_Utils
    api_Tests

docs/content/api_core/api_Utils.rst

@@ -21,7 +21,7 @@ Solver Utilities
     :undoc-members:
 
 Curv Utilities
-=============
+==============
 
 .. automodule:: SimPEG.Utils.curvutils
     :members:
@@ -51,7 +51,9 @@ Interpolation Utilities
 Counter Utilities
 =================
 
-::
+.. code-block:: python
+    :linenos:
+
     class MyClass(object):
         def __init__(self, url):
             self.counter = Counter()
@@ -69,7 +71,9 @@ Counter Utilities
     for i in range(300): c.MySecondMethod()
     c.counter.summary()
 
-::
+
+.. code-block:: text
+    :linenos:
 
     Counters:
       MyClass.MyMethod                        :      100
@@ -77,6 +81,8 @@ Counter Utilities
     Times:                                        mean      sum
       MyClass.MySecondMethod                  : 1.70e-06, 5.10e-04,  300x
 
+
+
 The API
 -------
 

docs/content/api_core/api_bigPicture.rst

@@ -35,7 +35,7 @@ The Big Picture
 Defining a well-posed inverse problem and solving it is a complex task that requires many components that must interact. It is helpful
 to view this task as a workflow in which various elements are explicitly identified and integrated. The figure below outlines the inversion components that consists of inputs, implementation, and evaluation. The inputs are composed of the geophysical data, the equations which are a mathematical description of the governing physics, and prior knowledge or assumptions about the setting. The implementation consists of two broad categories: the forward simulation and the inversion. The **forward simulation** is the means by which we solve the governing equations given a model and the **inversion components** evaluate and update this model. We are considering a gradient based approach, which updates the model through an optimization routine. The output of this implementation is a model, which, prior to interpretation, must be evaluated. This requires considering, and often re-assessing, the choices and assumptions made in both the input and implementation stages.
 
-.. image:: InversionWorkflow-PreSimPEG.png
+.. image:: ../../images/InversionWorkflow-PreSimPEG.png
    :width: 400 px
    :alt: Components
    :align: center
@@ -46,24 +46,24 @@ A Comprehensive Framework
 
 There are an overwhelming amount of choices to be made as one works through the forward modeling and inversion process (see figure above). As a result, software implementations of this workflow often become complex and highly interdependent, making it difficult to interact with and to ask other scientists to pick up and change. Our approach to handling this complexity is to propose a framework, (see below), that compartmentalizes the implementation of inversions into various units. We present it in this specific modular style, as each unit contains a targeted subset of choices crucial to the inversion process.
 
-.. image:: InversionWorkflow.png
+.. image:: ../../images/InversionWorkflow.png
    :width: 400 px
    :alt: Framework
    :align: center
 
-The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The (:class:`SimPEG.Mesh.BaseMesh`) class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the (:class:`SimPEG.Survey.BaseSurvey`) and the (:class:`SimPEG.Problem.BaseProblem`). The (:class:`SimPEG.Survey.BaseSurvey`) class handles the geometry of a geophysical problem as well as sources. The (:class:`SimPEG.Problem.BaseProblem`) class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The (:class:`SimPEG.Problem.BaseProblem`) creates geophysical fields given a source from the (:class:`SimPEG.Survey.BaseSurvey`). The (:class:`SimPEG.Survey.BaseSurvey`) interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a (:class:`SimPEG.DataMisfit.BaseDataMisfit`) is chosen to capture the goodness of fit of the predicted data and a (:class:`SimPEG.Regularization.BaseRegularization`) is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class (:class:`SimPEG.InvProblem.BaseInvProblem`). (:class:`SimPEG.InvProblem.BaseInvProblem`) is the mathematical statement that will be numerically solved by running an Inversion. The (:class:`SimPEG.Inversion.BaseInversion`) class handles organization and dispatch of directives between all of the various pieces of the framework.
+The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the :class:`SimPEG.Survey.BaseSurvey` and the :class:`SimPEG.Problem.BaseProblem`. The :class:`SimPEG.Survey.BaseSurvey` class handles the geometry of a geophysical problem as well as sources. The :class:`SimPEG.Problem.BaseProblem` class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The :class:`SimPEG.Problem.BaseProblem` creates geophysical fields given a source from the :class:`SimPEG.Survey.BaseSurvey`. The :class:`SimPEG.Survey.BaseSurvey` interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a :class:`SimPEG.DataMisfit.BaseDataMisfit` is chosen to capture the goodness of fit of the predicted data and a :class:`SimPEG.Regularization.BaseRegularization` is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class :class:`SimPEG.InvProblem.BaseInvProblem`. :class:`SimPEG.InvProblem.BaseInvProblem` is the mathematical statement that will be numerically solved by running an Inversion. The :class:`SimPEG.Inversion.BaseInversion` class handles organization and dispatch of directives between all of the various pieces of the framework.
 
-The arrows in the figure above indicate what each class takes as a primary argument. For example, both the (:class:`SimPEG.Problem.BaseProblem`) and (:class:`SimPEG.Regularization.BaseRegularization`) classes take a (:class:`SimPEG.Mesh.BaseMesh`) class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The (:class:`SimPEG.Mesh.BaseMesh`) class, for example, could be a regular Cartesian mesh (:class:`SimPEG.Mesh.TensorMesh`) or a cylindrical coordinate mesh (:class:`SimPEG.Mesh.CylMesh`), which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
+The arrows in the figure above indicate what each class takes as a primary argument. For example, both the :class:`SimPEG.Problem.BaseProblem` and :class:`SimPEG.Regularization.BaseRegularization` classes take a :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class, for example, could be a regular Cartesian mesh :class:`SimPEG.Mesh.TensorMesh` or a cylindrical coordinate mesh :class:`SimPEG.Mesh.CylMesh`, which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
 
 
-.. include:: ../CITATION.rst
+.. include:: ../../../CITATION.rst
 
 Authors
 -------
 
-.. include:: ../AUTHORS.rst
+.. include:: ../../../AUTHORS.rst
 
 License
 -------
 
-.. include:: ../LICENSE
+.. include:: ../../../LICENSE

docs/content/api_core/api_installing.rst

@@ -66,7 +66,7 @@ Numpy and Matlab
 Lessons in Python
 -----------------
 
-* `Software Carpentry <http://software-carpentry.org/v4/python/index.html>`_
+* `Software Carpentry <http://swcarpentry.github.io/python-novice-inflammation/>`_
 * `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_
 
 Editing Python

docs/content/dc/index.rst

@@ -1,5 +1,3 @@
-.. _api_DC:
-
 .. math::
 
     \renewcommand{\div}{\nabla\cdot\,}
@@ -38,10 +36,27 @@
     \renewcommand {\u}  { {\vec u} }
     \newcommand{\I}{\vec{I}}
 
-DC resistivity survey
-*********************
 
-Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
+Direct Current Resistivity
+**************************
+
+`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
+direct current (DC) resistivity and induced polarization (IP) geophysical problems.
+
+
+DC resistivity survey
+=====================
+
+Electrical resistivity of subsurface materials is measured by causing an
+electrical current to flow in the earth between one pair of electrodes while
+the voltage across a second pair of electrodes is measured. The result is an
+"apparent" resistivity which is a value representing the weighted average
+resistivity over a volume of the earth. Variations in this measurement are
+caused by variations in the soil, rock, and pore fluid electrical resistivity.
+Surveys require contact with the ground, so they can be labour intensive.
+Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D
+models are estimated using inversion procedures (`GPG
+<http://gpg.geosci.xyz>`_).
 
 
 Background
@@ -55,7 +70,7 @@ As direct current (DC) implies, in DC resistivity survey, we assume steady-state
 
     \curl \e = 0
 
-Then by taking \\(\\curl\\) for the first equation, we have
+Then by taking \\(\\div\\) of the first equation, we have
 
 .. math::
 
@@ -137,13 +152,14 @@ Comparing to the analytic function:
 
 .. plot::
 
-    import simpegDC as DC
+    from SimPEG import Examples
-    DC.Examples.Verification.run(plotIt=True)
+    Examples.DC_Analytic_Dipole.run(plotIt=True)
 
-API
-===
 
-.. automodule:: simpegDC.BaseDC
+API for DC codes
+================
+
+.. automodule:: SimPEG.DCIP.BaseDC
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/em/api_FDEM.rst

@@ -9,17 +9,28 @@
 Frequency Domain Electromagnetics
 *********************************
 
-Electromagnetic (EM) geophysical methods are used in a variety of applications from resource exploration, including for hydrocarbons and minerals, to environmental applications, such as groundwater monitoring. The primary physical property of interest in EM is electrical conductivity, which describes the ease with which electric current flows through a material.
+Electromagnetic (EM) geophysical methods are used in a variety of applications
+from resource exploration, including for hydrocarbons and minerals, to
+environmental applications, such as groundwater monitoring. The primary
+physical property of interest in EM is electrical conductivity, which
+describes the ease with which electric current flows through a material.
 
 
 Background
 ==========
 
-Electromagnetic phenomena are governed by Maxwell's equations. They describe the behavior of EM fields and fluxes. Electromagnetic theory for geophysical applications by Ward and Hohmann (1988) is a highly recommended resource on this topic.
+Electromagnetic phenomena are governed by Maxwell's equations. They describe
+the behavior of EM fields and fluxes. Electromagnetic theory for geophysical
+applications by Ward and Hohmann (1988) is a highly recommended resource on
+this topic.
 
 Fourier Transform Convention
 ----------------------------
-In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the :math:`e^{i \omega t}` convention, so we define our Fourier Transform pair as
+
+In order to examine Maxwell's equations in the frequency domain, we must first
+define our choice of harmonic time-dependence by choosing a Fourier transform
+convention. We use the :math:`e^{i \omega t}` convention, so we define our
+Fourier Transform pair as
 
 .. math ::
     F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\
@@ -31,6 +42,7 @@ where :math:`\omega` is angular frequency, :math:`t` is time, :math:`F(\omega)`
 
 Maxwell's Equations
 ===================
+
 In the frequency domain, Maxwell's equations are given by
 
 .. math ::
@@ -104,19 +116,20 @@ The H-J formulation is in terms of the current density and the magnetic field:
 
 Discretizing
 ------------
+
 For both formulations, we use a finite volume discretization
 and discretize fields on cell edges, fluxes on cell faces and
 physical properties in cell centers. This is particularly
 important when using symmetry to reduce the dimensionality of a problem
 (for instance on a 2D CylMesh, there are :math:`r`, :math:`z` faces and :math:`\theta` edges)
 
-.. figure:: ../images/finitevolrealestate.png
+.. figure:: ../../images/finitevolrealestate.png
     :align: center
     :scale: 60 %
 
 For the two formulations, the discretization of the physical properties, fields and fluxes are summarized below.
 
-.. figure:: ../images/ebjhdiscretizations.png
+.. figure:: ../../images/ebjhdiscretizations.png
     :align: center
     :scale: 60 %
 
@@ -150,7 +163,7 @@ API
 FDEM Problem
 ------------
 
-.. automodule:: SimPEG.EM.FDEM.FDEM
+.. automodule:: SimPEG.EM.FDEM.ProblemFDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -169,6 +182,11 @@ FDEM Survey
     :members:
     :undoc-members:
 
+.. automodule:: SimPEG.EM.FDEM.RxFDEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
 FDEM Fields
 -----------
 

docs/content/em/api_TDEM.rst

@@ -359,7 +359,7 @@ TDEM - B formulation
 Field Storage
 =============
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
+.. autoclass:: SimPEG.EM.TDEM.BaseTDEM.FieldsTDEM
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/em/api_basic.rst

@@ -0,0 +1,33 @@
+Overview of Electromagnetics in SimPEG
+**************************************
+
+
+The API
+=======
+
+Physical Properties
+-------------------
+
+.. autoclass:: SimPEG.EM.Base.EMPropMap
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+Problem
+-------
+
+.. autoclass:: SimPEG.EM.Base.BaseEMProblem
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+
+Survey
+------
+
+.. autoclass:: SimPEG.EM.Base.BaseEMSurvey
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+

docs/content/em/index.rst

@@ -3,22 +3,23 @@ Electromagnetics
 ================
 
 `SimPEG.EM` uses SimPEG as the framework for the forward and inverse
-electromagnetics geophysical problems. 
+electromagnetics geophysical problems.
 
 To solve for predicted data, we follow the framework shown below. The model is
 what we invert for. This is mapped to a physical property on the simulation
 mesh. A source which is used to excite the system is specified. Having a model
 and a source, we can solve Maxwell's equations for fields. We sample these
-fields with recievers to give us predicted data. 
+fields with recievers to give us predicted data.
 
 
-.. image:: ../images/simpegEM_noMath.png
+.. image:: ../../images/simpegEM_noMath.png
     :scale: 50%
 
 
 .. toctree::
    :maxdepth: 2
 
+   api_basic
    api_FDEM
    api_TDEM
    api_Utils

docs/content/examples/DC_Analytic_Dipole.rst

@@ -16,6 +16,6 @@ DC Analytic Dipole
     from SimPEG import Examples
     Examples.DC_Analytic_Dipole.run()
 
-.. literalinclude:: ../../SimPEG/Examples/DC_Analytic_Dipole.py
+.. literalinclude:: ../../../SimPEG/Examples/DC_Analytic_Dipole.py
     :language: python
     :linenos:

docs/content/examples/DC_Forward_PseudoSection.rst

@@ -31,6 +31,6 @@ Created by @fourndo
     from SimPEG import Examples
     Examples.DC_Forward_PseudoSection.run()
 
-.. literalinclude:: ../../SimPEG/Examples/DC_Forward_PseudoSection.py
+.. literalinclude:: ../../../SimPEG/Examples/DC_Forward_PseudoSection.py
     :language: python
     :linenos:

docs/content/examples/EM_FDEM_1D_Inversion.rst

@@ -21,6 +21,6 @@ Here we will create and run a FDEM 1D inversion.
     from SimPEG import Examples
     Examples.EM_FDEM_1D_Inversion.run()
 
-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
     :language: python
     :linenos:

docs/content/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst

@@ -21,6 +21,6 @@ Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
     from SimPEG import Examples
     Examples.EM_FDEM_Analytic_MagDipoleWholespace.run()
 
-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
     :language: python
     :linenos: