Updates to documentation and DC utilities.

Merge remote-tracking branch 'refs/remotes/origin/dev'
Conflicts: SimPEG/DCIP/BaseDC.py SimPEG/DCIP/BaseIP.py
2026-09-14 11:35:32 +08:00 · 2016-06-26 16:07:35 -06:00 · 2016-06-26 14:08:28 -06:00 · 2016-05-31 14:09:57 -07:00 · 2016-05-31 08:42:52 -07:00
84 changed files with 1112 additions and 3305 deletions
@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.12
+current_version = 0.1.11
 files = setup.py SimPEG/__init__.py docs/conf.py
@@ -40,4 +40,3 @@ nosetests.xml
 docs/_build/
 Makefile
 docs/warnings.txt
 .DS_Store
@@ -18,7 +18,7 @@ env:
    - TEST_DIR="tests/mesh tests/base tests/utils"
    - TEST_DIR=tests/em/fdem/inverse/derivs
    - TEST_DIR=tests/em/tdem
-    - TEST_DIR=tests/dcip
+    - TEST_DIR=tests/em/static
    - TEST_DIR=tests/flow
    - TEST_DIR=tests/mt
    - TEST_DIR=tests/examples
@@ -58,7 +58,8 @@ script:
 # Calculate coverage
 after_success:
-  - bash <(curl -s https://codecov.io/bash)
+  - coveralls --config_file .coveragerc
  - 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}");
@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
    :alt: SimPEG Logo
 ======
@@ -21,21 +21,14 @@ 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&amp;utm_medium=referral&amp;utm_content=simpeg/simpeg&amp;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:
@@ -1,292 +0,0 @@
 from SimPEG import *
 class FieldsDC_CC(Problem.Fields):
    knownFields = {'phi_sol':'CC'}
    aliasFields = {
                    'phi' : ['phi_sol','CC','_phi'],
                    'e' : ['phi_sol','F','_e'],
                    'j' : ['phi_sol','F','_j']
                  }
    def __init__(self,mesh,survey,**kwargs):
        super(FieldsDC_CC, self).__init__(mesh, survey, **kwargs)
    def startup(self):
        self._cellGrad = self.survey.prob.mesh.cellGrad
        self._Mfinv = self.survey.prob.mesh.getFaceInnerProduct(invMat=True)
    def _phi(self, phi_sol, srcList):
        phi = phi_sol
        # for i, src in enumerate(srcList):
        #     phi_p = src.phi_p(self.survey.prob)
        #     if phi_p is not None:
        #         phi[:,i] += phi_p
        return phi
    def _e(self, phi_sol, srcList):
        e = -self._cellGrad*phi_sol
        # for i, src in enumerate(srcList):
        #     e_p = src.e_p(self.survey.prob)
        #     if e_p is not None:
        #         e[:,i] += e_p
        return e
    def _j(self, phi_sol, srcList):
        j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol
        # for i, src in enumerate(srcList):
        #     j_p = src.j_p(self.survey.prob)
        #     if j_p is not None:
        #         j[:,i] += j_p
        return j
 class SrcDipole(Survey.BaseSrc):
    """A dipole source, locA and locB are moved to the closest cell-centers"""
    current = 1
    loc  = None
    # _rhsDict = None
    def __init__(self, rxList, locA, locB, **kwargs):
        self.loc = (locA, locB)
        super(SrcDipole, self).__init__(rxList, **kwargs)
    def eval(self, prob):
        # Recompute rhs
        # if getattr(self, '_rhsDict', None) is None:
        #     self._rhsDict = {}
        # if mesh not in self._rhsDict:
        pts = [self.loc[0], self.loc[1]]
        inds = Utils.closestPoints(prob.mesh, pts)
        q = np.zeros(prob.mesh.nC)
        q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
        # self._rhsDict[mesh] = q
        # return self._rhsDict[mesh]
        return q
 class RxDipole(Survey.BaseRx):
    """A dipole source, locA and locB are moved to the closest cell-centers"""
    def __init__(self, locsM, locsN, **kwargs):
        locs = (locsM, locsN)
        assert locsM.shape == locsN.shape, 'locs must be the same shape.'
        super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)
    @property
    def nD(self):
        """Number of data in the receiver."""
        return self.locs[0].shape[0]
    def getP(self, mesh):
        P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc)
        P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc)
        return P0 - P1
 class SurveyDC(Survey.BaseSurvey):
    """
        **SurveyDC**
        Geophysical DC resistivity data.
    """
    uncert = None
    def __init__(self, srcList, **kwargs):
        self.srcList = srcList
        Survey.BaseSurvey.__init__(self, **kwargs)
        # self._rhsDict = {}
        self._Ps = {}
    def eval(self, u):
        """
            Predicted data.
            .. math::
                d_\\text{pred} = Pu(m)
        """
        P = self.getP(self.prob.mesh)
        return P*mkvc(u[self.srcList, 'phi_sol'])
    def getP(self, mesh):
        if mesh in self._Ps:
            return self._Ps[mesh]
        P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList]
        self._Ps[mesh] = sp.block_diag(P_src)
        return self._Ps[mesh]
 class ProblemDC_CC(Problem.BaseProblem):
    """
        **ProblemDC**
        Geophysical DC resistivity problem.
    """
    surveyPair = SurveyDC
    Solver     = Solver
    fieldsPair = FieldsDC_CC
    Ainv = None
    def __init__(self, mesh, **kwargs):
        Problem.BaseProblem.__init__(self, mesh)
        self.mesh.setCellGradBC('neumann')
        Utils.setKwargs(self, **kwargs)
    deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
    @property
    def Msig(self):
        if getattr(self, '_Msig', None) is None:
            sigma = self.curModel.transform
            Av = self.mesh.aveF2CC
            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
        return self._Msig
    @property
    def dMdsig(self):
        if getattr(self, '_dMdsig', None) is None:
            sigma = self.curModel.transform
            Av = self.mesh.aveF2CC
            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
        return self._dMdsig
    @property
    def A(self):
        """
            Makes the matrix A(m) for the DC resistivity problem.
            :param numpy.ndarray m: model
            :rtype: scipy.sparse.csc_matrix
            :return: A(m)
            .. math::
                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
            Where M() is the mass matrix and mT is the model transform.
        """
        if getattr(self, '_A', None) is None:
            D = self.mesh.faceDiv
            G = self.mesh.cellGrad
            self._A = D*self.Msig*G
            # Remove the null space from the matrix.
            self._A[0,0] /= self.mesh.vol[0]
            self._A = self._A.tocsc()
        return self._A
    def getRHS(self):
        # if self.mesh not in self._rhsDict:
        RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
        # self._rhsDict[mesh] = RHS
        # return self._rhsDict[mesh]
        return RHS
    def fields(self, m):
        F = self.fieldsPair(self.mesh, self.survey)
        self.curModel = m
        A    = self.A
        self.Ainv = self.Solver(A, **self.solverOpts)
        RHS    = self.getRHS()
        Phi  = self.Ainv * RHS
        Srcs = self.survey.srcList
        F[Srcs, 'phi_sol'] = Phi
        return F
    def Jvec(self, m, v, f=None):
        """
            :param numpy.array m: model
            :param numpy.array v: vector to multiply
            :param Fields f: fields
            :rtype: numpy.array
            :return: Jv
            .. math::
                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
                \\nabla_u (A(m)u - q) = A(m)
                \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
            Where M() is the mass matrix and mT is the model transform.
            .. math::
                J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
                J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
        """
        # Set current model; clear dependent property $\mathbf{A(m)}$
        self.curModel = m
        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        if f is None:
            # Run forward simulation if $u$ not provided
            f = self.fields(self.curModel)
        u = f[self.survey.srcList, 'phi_sol']
        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        # Derivative of model transform, $\deriv{\sigma}{\m}$
        dsigdm_x_v = self.curModel.transformDeriv * v
        # Take derivative of $C(m,u)$ w.r.t. $m$
        dCdm_x_v = np.empty_like(u)
        # loop over fields for each source
        for i in range(self.survey.nSrc):
            # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
            dAdsig         = D * self.dMdsig( G * u[:,i] )
            dCdm_x_v[:, i] = dAdsig *  dsigdm_x_v
        # Take derivative of $C(m,u)$ w.r.t. $u$
        dA_du = self.A
        # Solve for $\deriv{u}{m}$
        # dCdu_inv = self.Solver(dCdu, **self.solverOpts)
        if self.Ainv is None:
            self.Ainv = self.Solver(dA_du, **self.solverOpts)
        P     = self.survey.getP(self.mesh)
        Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
        return Jv
    def Jtvec(self, m, v, f=None):
        self.curModel = m
        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        if f is None:
            # Run forward simulation if $f$ not provided
            f = self.fields(self.curModel)
        u = f[self.survey.srcList, 'phi_sol']
        shp = u.shape
        P = self.survey.getP(self.mesh)
        PT_x_v = (P.T*v).reshape(shp, order='F')
        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        dA_du = self.A
        mT_dm = self.mapping.deriv(m)
        # We probably always need this due to the linesearch .. (?)
        self.Ainv = self.Solver(dA_du.T, **self.solverOpts)
        # if self.Ainv is None:
        #     self.Ainv = self.Solver(dCdu, **self.solverOpts)
        w = self.Ainv * PT_x_v
        Jtv = 0
        for i, ui in enumerate(u.T):  # loop over each column
            Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
        Jtv = - mT_dm.T * ( Jtv )
        return Jtv
@@ -1,182 +0,0 @@
 from SimPEG import *
 from BaseDC import SurveyDC, FieldsDC_CC
 class SurveyIP(SurveyDC):
    """
        **SurveyDC**
        Geophysical DC resistivity data.
    """
    def __init__(self, srcList, **kwargs):
        self.srcList = srcList
        Survey.BaseSurvey.__init__(self, **kwargs)
        self._Ps = {}
    def dpred(self, m, f=None):
        """
            Predicted data.
            .. math::
                d_\\text{pred} = Pf(m)
        """
        return self.prob.forward(m)
 class ProblemIP(Problem.BaseProblem):
    """
        **ProblemIP**
        Geophysical IP resistivity problem.
    """
    surveyPair = SurveyDC
    Solver     = Solver
    sigma = None
    Ainv = None
    u = None
    def __init__(self, mesh, **kwargs):
        Problem.BaseProblem.__init__(self, mesh)
        self.mesh.setCellGradBC('neumann')
        Utils.setKwargs(self, **kwargs)
    # deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
    @property
    def Msig(self):
        if getattr(self, '_Msig', None) is None:
            # sigma = self.curModel.transform
            sigma = self.sigma
            Av = self.mesh.aveF2CC
            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
        return self._Msig
    @property
    def dMdsig(self):
        if getattr(self, '_dMdsig', None) is None:
            # sigma = self.curModel.transform
            sigma = self.sigma
            Av = self.mesh.aveF2CC
            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
        return self._dMdsig
    @property
    def A(self):
        """
            Makes the matrix A(m) for the DC resistivity problem.
            :param numpy.array m: model
            :rtype: scipy.sparse.csc_matrix
            :return: A(m)
            .. math::
                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
            Where M() is the mass matrix and mT is the model transform.
        """
        if getattr(self, '_A', None) is None:
            D = self.mesh.faceDiv
            G = self.mesh.cellGrad
            self._A = D*self.Msig*G
            # Remove the null space from the matrix.
            self._A[-1,-1] /= self.mesh.vol[-1]
            self._A = self._A.tocsc()
        return self._A
    def getRHS(self):
        # if self.mesh not in self._rhsDict:
        RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
        # self._rhsDict[mesh] = RHS
        # return self._rhsDict[mesh]
        return RHS
    def fields(self, m):
        if self.u is None:
 	        A    = self.A
 	        if self.Ainv == None:
 	        	self.Ainv = self.Solver(A, **self.solverOpts)
 	        Q    = self.getRHS()
 	        self.u  = self.Ainv * Q
        return self.u
    def forward(self, m, u=None):
        # Set current model; clear dependent property $\mathbf{A(m)}$
        self.curModel = m
        # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        sigma = self.sigma
        if self.u is None:
            # Run forward simulation if $u$ not provided
            u = self.fields(sigma)
        shp = (self.mesh.nC, self.survey.nSrc)
        u = self.u.reshape(shp, order='F')
        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        # Derivative of model transform, $\deriv{\sigma}{\m}$
        # dsigdm_x_v = self.curModel.transformDeriv * v
        dsigdm_x_v = Utils.sdiag(sigma) * self.curModel.transformDeriv *  m
        # Take derivative of $C(m,u)$ w.r.t. $m$
        dCdm_x_v = np.empty_like(u)
        # loop over fields for each source
        for i in range(self.survey.nSrc):
            # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
            dAdsig         = D * self.dMdsig( G * u[:,i] )
            dCdm_x_v[:, i] = dAdsig *  dsigdm_x_v
        # Take derivative of $C(m,u)$ w.r.t. $u$
        if self.Ainv == None:
        	self.Ainv = self.Solver(A, **self.solverOpts)
        # dCdu = self.A
        # Solve for $\deriv{u}{m}$
        # dCdu_inv = self.Solver(dCdu, **self.solverOpts)
        P        = self.survey.getP(self.mesh)
        J_x_v    = - P * mkvc( self.Ainv * dCdm_x_v )
        return -J_x_v
    def Jvec(self, m, v, f=None):
        return self.forward(v)
    def Jtvec(self, m, v, f=None):
        self.curModel = m
        # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
        sigma = self.sigma
        if self.u is None:
            u = self.fields(sigma)
        else:
        	u = self.u
        shp = (self.mesh.nC, self.survey.nSrc)
        u = u.reshape(shp, order='F')
        P = self.survey.getP(self.mesh)
        PT_x_v = (P.T*v).reshape(shp, order='F')
        D = self.mesh.faceDiv
        G = self.mesh.cellGrad
        A = self.A
        mT_dm = Utils.sdiag(sigma)*self.mapping.deriv(m)
        # mT_dm = self.mapping.deriv(m)
        # dCdu = A.T
        # Ainv = self.Solver(dCdu, **self.solverOpts)
        # if self.Ainv == None:
    	self.Ainv = self.Solver(A.T, **self.solverOpts)
        w = self.Ainv * PT_x_v
        Jtv = 0
        for i, ui in enumerate(u.T):  # loop over each column
            Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
        Jtv = - mT_dm.T * ( Jtv )
        return -Jtv
@@ -1,38 +0,0 @@
 import numpy as np
 def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
    import SimPEG.DCIP as DC
    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
    elocs -= (nElecs*aSpacing - aSpacing)/2
    space = 1
    WENNER = np.zeros((0,),dtype=int)
    for ii in range(nElecs):
        for jj in range(nElecs):
            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
            if np.any(test >= nElecs):
                break
            WENNER = np.r_[WENNER, test]
        space += 1
    WENNER = WENNER.reshape((-1,4))
    if plotIt:
        for i, s in enumerate('rbkg'):
            plt.plot(elocs[WENNER[:,i]],s+'.')
        plt.show()
    # Create sources and receivers
    i = 0
    if in2D:
        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
    else:
        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
    srcList = []
    for i in range(WENNER.shape[0]):
        rx = DC.RxDipole(getLoc(i,1),getLoc(i,2))
        src = DC.SrcDipole([rx], getLoc(i,0),getLoc(i,3))
        srcList += [src]
    return srcList
@@ -1,4 +0,0 @@
 from BaseDC import *
 from BaseIP import *
 from DCIPUtils import *
 import Utils
@@ -15,7 +15,7 @@ class InversionDirective(object):
    @inversion.setter
    def inversion(self, i):
        if getattr(self,'_inversion',None) is not None:
-            print 'Warning: InversionDirective {0!s} has switched to a new inversion.'.format(self.__name__)
+            print 'Warning: InversionDirective %s has switched to a new inversion.' % 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 {0!s}'.format(d.__name__)
+            assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not %s' % 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: {0!s} has switched to a new inversion.'.format(self.__name__)
+            print 'Warning: %s has switched to a new inversion.' % 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 ["{0!s}"]'.format('", "'.join(directives))
+        assert ruleType in directives, 'Directive type must be in ["%s"]' % '", "'.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: {0:d}'.format(self.opt.iter)
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
            self.invProb.beta /= self.coolingFactor
@@ -181,7 +181,7 @@ class SaveEveryIteration(InversionDirective):
    def fileName(self):
        if getattr(self, '_fileName', None) is None:
            from datetime import datetime
-            self._fileName = '{0!s}-{1!s}'.format(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
+            self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
        return self._fileName
    @fileName.setter
    def fileName(self, value):
@@ -192,31 +192,31 @@ class SaveModelEveryIteration(SaveEveryIteration):
    """SaveModelEveryIteration"""
    def initialize(self):
-        print "SimPEG.SaveModelEveryIteration will save your models as: '###-{0!s}.npy'".format(self.fileName)
+        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
    def endIter(self):
-        np.save('{0:03d}-{1!s}'.format(self.opt.iter, self.fileName), self.opt.xc)
+        np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
 class SaveOutputEveryIteration(SaveEveryIteration):
    """SaveModelEveryIteration"""
    def initialize(self):
-        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-{0!s}.txt'".format(self.fileName)
+        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
        f = open(self.fileName+'.txt', 'w')
        f.write("  #     beta     phi_d     phi_m       f\n")
        f.close()
    def endIter(self):
        f = open(self.fileName+'.txt', 'a')
-        f.write(' {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.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
        f.close()
 class SaveOutputDictEveryIteration(SaveEveryIteration):
    """SaveOutputDictEveryIteration"""
    def initialize(self):
-        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-{0!s}.npz'".format(self.fileName)
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
    def endIter(self):
        # Save the data.
@@ -253,7 +253,8 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
 class Update_IRLS(InversionDirective):
    eps_min = None
-    eps = None
+    eps_p = None
    eps_q = None
    norms = [2.,2.,2.,2.]
    factor = None
    gamma = None
@@ -262,7 +263,6 @@ class Update_IRLS(InversionDirective):
    f_old = None
    f_min_change = 1e-2
    beta_tol = 5e-2
    prctile = 95
    # Solving parameter for IRLS (mode:2)
    IRLSiter   = 0
@@ -297,22 +297,9 @@ class Update_IRLS(InversionDirective):
            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
            self.mode = 2
-
+            print self.eps_p, self.eps_q, self.norms
-            # Either use the supplied epsilon, or fix base on distribution of
+            self.reg.eps_p = self.eps_p
-            # model values
+            self.reg.eps_q = self.eps_q
            if getattr(self, 'reg.eps', None) is None:
                self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
            else:                 
                self.reg.eps_p = self.eps[0]
            if getattr(self, 'reg.eps', None) is None:
                self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
            else:                 
                self.reg.eps_q = self.eps[1]
            print "L[p qx qy qz]-norm : " + str(self.reg.norms)
            print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
            self.reg.norms = self.norms
            self.coolingFactor = 1.
            self.coolingRate = 1
@@ -328,7 +315,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: {0:d}'.format(self.opt.iter)
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
            self.invProb.beta /= self.coolingFactor
@@ -340,11 +327,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: {0:6.3e}".format((self.f_change))
+            print "Regularization decrease: %6.3e" % (self.f_change)
            # Check for maximum number of IRLS cycles
            if self.IRLSiter == self.maxIRLSiter:
-                print "Reach maximum number of IRLS cycles: {0:d}".format(self.maxIRLSiter)
+                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
                self.opt.stopNextIteration = True
                return
@@ -356,14 +343,14 @@ class Update_IRLS(InversionDirective):
            else:
                self.f_old = phim_new
-#            # Cool the threshold parameter if required
+            # Cool the threshold parameter if required
-#            if getattr(self, 'factor', None) is not None:
+            if getattr(self, 'factor', None) is not None:
-#                eps = self.reg.eps / self.factor
+                eps = self.reg.eps / self.factor
-#
+
-#                if getattr(self, 'eps_min', None) is not None:
+                if getattr(self, 'eps_min', None) is not None:
-#                    self.reg.eps = np.max([self.eps_min,eps])
+                    self.reg.eps = np.max([self.eps_min,eps])
-#                else:
+                else:
-#                    self.reg.eps = eps
+                    self.reg.eps = eps
            # Get phi_m at the end of current iteration
            self.phi_m_last = self.invProb.phi_m_last
@@ -114,7 +114,7 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l
    """
    mu = mu_0*(1+kappa)
    epsilon = epsilon_0*epsr
-    sig_hat = sig + 1j*omega(f)*epsilon
+    sig_hat = sig + 1j*omeg*epsilon
    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
    # Check
@@ -160,7 +160,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1.,
            Add description of parameters
    """
-    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
    Jx = sig*Ex
    Jy = sig*Ey
    Jz = sig*Ez
@@ -175,7 +175,7 @@ def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le
            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)
+    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
    Jx_galvanic = sig*Ex_galvanic
    Jy_galvanic = sig*Ey_galvanic
    Jz_galvanic = sig*Ez_galvanic
@@ -190,7 +190,7 @@ def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l
            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)
+    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
    Jx_inductive = sig*Ex_inductive
    Jy_inductive = sig*Ey_inductive
    Jz_inductive = sig*Ez_inductive
@@ -248,7 +248,7 @@ def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1.,
            Add description of parameters
    """
-    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
    Bx = mu*Hx
    By = mu*Hy
    Bz = mu*Hz
@@ -42,7 +42,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
        return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -56,7 +56,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
        return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -70,7 +70,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
        return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -84,7 +84,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
        return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -100,7 +100,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -118,7 +118,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -136,7 +136,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -154,7 +154,7 @@ class FieldsFDEM(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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
@@ -11,8 +11,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
    """
    def __init__(self, locs, orientation=None, component=None):
-        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(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(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not {0!s}".format(component)
+        assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
        self.projComp = orientation
        self.component = component
@@ -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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
@@ -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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %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 {0!s} is not implemented'.format(self.knownFields.keys()[0]))
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
        if adjoint:
            return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
@@ -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, '_{0!s}Deriv'.format(rx.projField), None)
+                df_dmFun = getattr(f, '_%sDeriv'%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, '_{0!s}Deriv'.format(rx.projField), None)
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
                df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
                ATinvdf_duT = self.Ainv * df_duT
@@ -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, '_{0!s}Deriv'.format(rx.projField), None)
+                    df_dmFun = getattr(f, '_%sDeriv'%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, '_{0!s}Deriv'.format(rx.projField), None)
+                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
                    df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
                    ATinvdf_duT = self.Ainv[iky] * df_duT
@@ -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, '_{0!s}Deriv'.format(rx.projField), None)
+                df_dmFun = getattr(f, '_%sDeriv'%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, '_{0!s}Deriv'.format(rx.projField), None)
+                df_duTFun = getattr(f, '_%sDeriv'%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)
@@ -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, '_{0!s}Deriv'.format(rx.projField), None)
+                        df_dmFun = getattr(f, '_%sDeriv'%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, '_{0!s}Deriv'.format(rx.projField), None)
+                        df_dmFun = getattr(f, '_%sDeriv'%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, '_{0!s}Deriv'.format(rx.projField), None)
+                        df_duTFun = getattr(f, '_%sDeriv'%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)
@@ -219,10 +219,12 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
        for ii in range(0, int(nstn)-1):
-            if stype == 'dpdp':
+            if stype == 'dipole-dipole':
                tx = np.c_[M[ii,:],N[ii,:]]
-            elif stype == 'pdp':
+            elif stype == 'pole-dipole':
                tx = np.c_[M[ii,:],M[ii,:]]
            else:
                raise Exception('The stype must be "dipole-dipole" or "pole-dipole"')
            # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
@@ -256,10 +258,10 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
                P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
                rxClass = DC.Rx.Dipole_ky(P1, P2)
-            if stype == 'dpdp':
+            if stype == 'dipole-dipole':
-                    srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
+                srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
-            elif stype == 'pdp':
+            elif stype == 'pole-dipole':
-                    srcClass = DC.Src.Pole([rxClass], M[ii,:])
+                srcClass = DC.Src.Pole([rxClass], M[ii,:])
            SrcList.append(srcClass)
    elif stype == 'gradient':
@@ -310,7 +312,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
            srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
        SrcList.append(srcClass)
    else:
-        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
+        print """stype must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
    return SrcList
@@ -1,3 +1,4 @@
 import DC
 import IP
 import SIP
 import Utils
@@ -47,7 +47,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
        self.waveformType = "GENERAL"
    def fields(self, m):
-        if self.verbose: print '{0!s}\nCalculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*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 '{0!s}\nDone calculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*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 = {0:e})'.format(dt)
+                if self.verbose: print 'Factoring...   (dt = %e)'%dt
                Ainv = self.Solver(A, **self.solverOpts)
                if self.verbose: print 'Done'
            rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving...   (tInd = {0:d})'.format(tInd)
+            if self.verbose: print '    Solving...   (tInd = %d)'%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 = {0:e})'.format(dt)
+                if self.verbose: print 'Factoring (Adjoint)...   (dt = %e)'%dt
                Ainv = self.Solver(A, **self.solverOpts)
                if self.verbose: print 'Done'
            rhs = RHS(tInd, F)
-            if self.verbose: print '    Solving (Adjoint)...   (tInd = {0:d})'.format(tInd)
+            if self.verbose: print '    Solving (Adjoint)...   (tInd = %d)'%tInd
            sol = Ainv * rhs
            if self.verbose: print '    Done...'
            if sol.ndim == 1:
@@ -123,14 +123,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
            * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
        """
-        if self.verbose: print '{0!s}\nCalculating J(v)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*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 '{0!s}\nDone calculating J(v)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
        return - mkvc(Jv)
    def Jtvec(self, m, v, f=None):
@@ -148,7 +148,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
            * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
        """
-        if self.verbose: print '{0!s}\nCalculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*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 '{0!s}\nDone calculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
+        if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
        return - mkvc(w)
@@ -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 {0!s} problem'.format((fdemType))
+        print '  Fetching %s problem' % (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: {0!s}, {1!s} formulations - {2!s}'.format(fdemType1, fdemType2, comp)
+    print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
    logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
    mu = np.ones(mesh.nC)*MU
@@ -1,7 +1,7 @@
 from SimPEG import *
 import SimPEG.EM.Static.DC as DC
-def run(plotIt=True):
+def run(plotIt=False):
    cs = 25.
    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -65,4 +65,4 @@ def run(plotIt=True):
 if __name__ == '__main__':
-    print run()
+    print run(plotIt=True)
@@ -1,5 +1,5 @@
 from SimPEG import Mesh, Utils, np, sp
-import SimPEG.DCIP as DC
+from SimPEG.EM.Static import DC, Utils as StaticUtils
 import time
 def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
@@ -74,7 +74,13 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
    # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
    # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
-    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
+    srcList = StaticUtils.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
    print(srcList)
    raise Exception('The output of the function changed. We need to create the survey and get Rx to continue.')
    # Something like this!
    survey, Tx, Rx
    # Define some global geometry
    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
@@ -110,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: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
+    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
    print 'Number of cells', mesh.nC
    if plotIt is True:
@@ -42,33 +42,55 @@ def run(N=100, plotIt=True):
    survey = Survey.LinearSurvey()
    survey.pair(prob)
    survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
    #survey.makeSyntheticData(mtrue, std=std_noise)
    wd = np.ones(nk) * std_noise
    #print survey.std[0]
    #M = prob.mesh
    # Distance weighting
    wr = np.sum(prob.G**2.,axis=0)**0.5
    wr = ( wr/np.max(wr) )
 #    reg = Regularization.Simple(mesh)
 #    reg.mref = mref
 #    reg.cell_weights = wr
 #
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmis.Wd = 1./wd
-
+#
 #    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
 #    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
 #    invProb.curModel = m0
 #
 #    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
 #    target = Directives.TargetMisfit()
 #
    betaest = Directives.BetaEstimate_ByEig()
 #    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
 #
 #
 #    mrec = inv.run(m0)
 #    ml2 = mrec
 #    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
 #
 #    # Switch regularization to sparse
 #    phim = invProb.phi_m_last
 #    phid =  invProb.phi_d
    reg = Regularization.Sparse(mesh)
    reg.mref = mref
    reg.cell_weights = wr
    reg.mref = np.zeros(mesh.nC)
-    
+    eps_p = 5e-2
    eps_q = 5e-2
    norms   = [0., 0., 2., 2.]
    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
    update_Jacobi = Directives.Update_lin_PreCond()
-    
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
    # Set the IRLS directive, penalize the lowest 25 percentile of model values
    # Start with an l2-l2, then switch to lp-norms
    norms   = [0., 0., 2., 2.]    
    IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
@@ -1,62 +0,0 @@
 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()
@@ -1,41 +0,0 @@
 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()
@@ -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: {0:4.1f}'.format(PHIS[ii][0]))
+            ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
        plt.show()
 if __name__ == '__main__':
@@ -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],'{0:d}'.format(ii), color='r')
+            axes[0].text(loc[0]+0.2,loc[1],'%d'%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],'{0:d}'.format(ii), color='g')
+            axes[0].text(loc[0]+0.2,loc[1],'%d'%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,'{0:d}'.format((ii+M.nFx)), color='m')
+            axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
        axes[1].spy(M.faceDiv)
        axes[1].set_title('Face Divergence')
@@ -2,7 +2,7 @@ from SimPEG import *
 from SimPEG.Utils import surface2ind_topo
-def run(plotIt=True, nx=5, ny=5):
+def run(plotIt=False, nx=5, ny=5):
    """
        Utils: surface2ind_topo
@@ -10,8 +10,6 @@ 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
@@ -24,7 +22,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", "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"]
+__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"]
 ##### AUTOIMPORTS #####
@@ -59,7 +57,7 @@ if __name__ == '__main__':
        if line == "##### AUTOIMPORTS #####\n":
            inimports = not inimports
            if inimports:
-                out += '\n'.join(["import {0!s}".format(_) for _ in exfiles])
+                out += '\n'.join(["import %s"%_ for _ in exfiles])
                out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
                out += '\n##### AUTOIMPORTS #####\n'
    f.close()
@@ -76,11 +74,11 @@ if __name__ == '__main__':
        docstr = runFunction.__doc__
        if docstr is None:
-            doc = '{0!s}\n{1!s}'.format(name.replace('_',' '), '='*len(name))
+            doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
        else:
            doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
-        out = """.. _examples_{0!s}:
+        out = """.. _examples_%s:
 .. --------------------------------- ..
 ..                                   ..
@@ -90,21 +88,21 @@ if __name__ == '__main__':
 ..                                   ..
 .. --------------------------------- ..
-{1!s}
+%s
 .. plot::
    from SimPEG import Examples
-    Examples.{2!s}.run()
+    Examples.%s.run()
-.. literalinclude:: ../../../SimPEG/Examples/{3!s}.py
+.. literalinclude:: ../../../SimPEG/Examples/%s.py
    :language: python
    :linenos:
-""".format(name, doc, name, name)
+"""%(name,doc,name,name)
        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
-        print 'Creating: {0!s}.rst'.format(name)
+        print 'Creating: %s.rst'%name
        f = open(rst, 'w')
        f.write(out)
        f.close()
@@ -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 {0!s} class.".format((pair.__name__))
+        assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__)
@@ -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 ({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)
+            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)
        return u
    @Utils.timeIt
@@ -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 "{0:e} MB".format(sz)
+        return "%e MB"%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 ({0!s}) for setter, you can't set an aliased property".format(name))
+                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
            else:
-                raise KeyError('Invalid field name ({0!s}) for setter'.format(name))
+                raise KeyError('Invalid field name (%s) for setter'%name)
        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
-            raise KeyError('Invalid field name ({0!s}) for getter'.format(name))
+            raise KeyError('Invalid field name (%s) for getter'%name)
        return name
    def _indexAndNameFromKey(self, key, accessType):
@@ -101,7 +101,7 @@ class IdentityMap(object):
            :return: passed the test?
        """
-        print 'Testing {0!s}'.format(str(self))
+        print 'Testing %s' % 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 {0!s} class.".format((pair.__name__))
+            ), "Mapping object must be an instance of a %s class."%(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 {0!s} and {1!s}.'.format(str(self), str(val)))
+                raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
            return ComboMap([self, val])
        elif isinstance(val, np.ndarray):
            if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
-                raise ValueError('Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(str(self), str(val.shape)))
+                raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
            return self._transform(val)
        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
    def __str__(self):
-        return "{0!s}({1!s},{2!s})".format(self.__class__.__name__, self.shape[0], self.shape[1])
+        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
 class ComboMap(IdentityMap):
@@ -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[{0!s}] ({1!s}, {2!s}) and map[{3!s}] ({4!s}, {5!s}).'.format(*errArgs))
+                raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
            if isinstance(m, ComboMap):
                self.maps += m.maps
@@ -173,7 +173,7 @@ class ComboMap(IdentityMap):
        return deriv
    def __str__(self):
-        return 'ComboMap[{0!s}]({1!s},{2!s})'.format(' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
+        return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
 class ExpMap(IdentityMap):
@@ -502,9 +502,7 @@ class InjectActiveCells(IdentityMap):
        if Utils.isScalar(valInactive):
            self.valInactive = np.ones(self.nC)*float(valInactive)
        else:
-            self.valInactive = np.ones(self.nC)
+            self.valInactive = valInactive.copy()
            self.valInactive[self.indInactive] = valInactive.copy()
        self.valInactive[self.indActive] = 0
        inds = np.nonzero(self.indActive)[0]
@@ -522,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[{0:d}] must be a numpy array".format(i)
+                assert isinstance(x, np.ndarray), "x[%i] must be a numpy array" % i
                assert xi.size == x[0].size, "Number of elements in list must not change."
            x_array = np.ones((x.size, len(x)))
@@ -4,28 +4,14 @@ from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 from View import CurvView
 # Some helper functions.
-def length2D(x):
+length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
-    return (x[:, 0]**2 + x[:, 1]**2)**0.5
+length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**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))
-def length3D(x):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
    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.
@@ -45,16 +31,12 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
    _meshType = 'Curv'
    def __init__(self, nodes):
-        assert type(nodes) == list, ("'nodes' variable must be a list of "
+        assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
                                     "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[{0:d}] is not a"
+            assert isinstance(nodes_i, np.ndarray), ("nodes[%i] is not a numpy array." % i)
-                                                     "numpy array.".format(i))
+            assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % 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."
@@ -66,113 +48,121 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
        for i, node_i in enumerate(nodes):
            self._gridN[:, i] = Utils.mkvc(node_i.astype(float))
-    @property
+    def gridCC():
-    def gridCC(self):
+        doc = "Cell-centered grid."
        """
        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
-    @property
+        def fget(self):
-    def gridN(self):
+            if self._gridCC is None:
-        """
+                self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F')
-        Nodal grid.
+            return self._gridCC
-        """
+        return locals()
-        if getattr(self, '_gridN', None) is None:
+    _gridCC = None  # Store grid by default
-            raise Exception("Someone deleted this. I blame you.")
+    gridCC = property(**gridCC())
        return self._gridN
-    @property
+    def gridN():
-    def gridFx(self):
+        doc = "Nodal grid."
        """
        Face staggered grid in the x direction.
        """
-        if getattr(self, '_gridFx', None) is None:
+        def fget(self):
-            N = self.r(self.gridN, 'N', 'N', 'M')
+            if self._gridN is None:
-            if self.dim == 2:
+                raise Exception("Someone deleted this. I blame you.")
-                XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+            return self._gridN
-                self._gridFx = np.c_[XY[0], XY[1]]
+        return locals()
-            elif self.dim == 3:
+    _gridN = None  # Store grid by default
-                XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] +
+    gridN = property(**gridN())
                       n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
                self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
        return self._gridFx
-    @property
+    def gridFx():
-    def gridFy(self):
+        doc = "Face staggered grid in the x direction."
        """
        Face staggered grid in the y direction.
        """
-        if getattr(self, '_gridFy', None) is None:
+        def fget(self):
-            N = self.r(self.gridN, 'N', 'N', 'M')
+            if self._gridFx is None:
-            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]
+                if self.dim == 2:
-                self._gridFy = np.c_[XY[0], XY[1]]
+                    XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
-            elif self.dim == 3:
+                    self._gridFx = np.c_[XY[0], XY[1]]
-                XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] +
+                elif self.dim == 3:
-                       n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
+                    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]]
+                    self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-        return self._gridFy
+            return self._gridFx
        return locals()
    _gridFx = None  # Store grid by default
    gridFx = property(**gridFx())
-    @property
+    def gridFy():
-    def gridFz(self):
+        doc = "Face staggered grid in the y direction."
        """
        Face staggered grid in the y direction.
        """
-        if getattr(self, '_gridFz', None) is None:
+        def fget(self):
-            N = self.r(self.gridN, 'N', 'N', 'M')
+            if self._gridFy is None:
-            XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] +
+                N = self.r(self.gridN, 'N', 'N', 'M')
-                   n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
+                if self.dim == 2:
-            self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+                    XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
-        return self._gridFz
+                    self._gridFy = np.c_[XY[0], XY[1]]
                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())
-    @property
+    def gridFz():
-    def gridEx(self):
+        doc = "Face staggered grid in the z direction."
        """
        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
-    @property
+        def fget(self):
-    def gridEy(self):
+            if self._gridFz is None and self.dim == 3:
-        """
+                N = self.r(self.gridN, 'N', 'N', 'M')
-        Edge staggered grid in the y direction.
+                XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
-        """
+                self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
-        if getattr(self, '_gridEy', None) is None:
+            return self._gridFz
-            N = self.r(self.gridN, 'N', 'N', 'M')
+        return locals()
-            if self.dim == 2:
+    _gridFz = None  # Store grid by default
-                XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+    gridFz = property(**gridFz())
                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
-    @property
+    def gridEx():
-    def gridEz(self):
+        doc = "Edge staggered grid in the x direction."
-        """
+
-        Edge staggered grid in the z direction.
+        def fget(self):
-        """
+            if self._gridEx is None:
-        if getattr(self, '_gridEz', None) is None and self.dim == 3:
+                N = self.r(self.gridN, 'N', 'N', 'M')
-            N = self.r(self.gridN, 'N', 'N', 'M')
+                if self.dim == 2:
-            XYZ = [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._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
+                    self._gridEx = np.c_[XY[0], XY[1]]
-        return self._gridEz
+                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
        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 ---------------------
    #
@@ -204,94 +194,78 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
    #            | /              | /
    #            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."
-    @property
+        def fget(self):
-    def vol(self):
+            if(self._vol is None):
-        """
+                if self.dim == 2:
-        Construct cell volumes of the 3D model as 1d array
+                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
-        """
+                    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)
-        if getattr(self, '_vol', None) is None:
+                    vol1 = (Utils.volTetra(self.gridN, A, B, D, E) +  # cutted edge top
-            if self.dim == 2:
+                            Utils.volTetra(self.gridN, B, E, F, G) +  # cutted edge top
-                A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
+                            Utils.volTetra(self.gridN, B, D, E, G) +  # middle
-                normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros(
+                            Utils.volTetra(self.gridN, B, C, D, G) +  # cutted edge bottom
-                                              (self.nN, 1))], A, B, C, D)
+                            Utils.volTetra(self.gridN, D, E, G, H))   # cutted edge bottom
                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
+                    vol2 = (Utils.volTetra(self.gridN, A, F, B, C) +  # cutted edge top
-                        Utils.volTetra(self.gridN, B, E, F, G) +  # cutted edge top
+                            Utils.volTetra(self.gridN, A, E, F, H) +  # cutted edge top
-                        Utils.volTetra(self.gridN, B, D, E, G) +  # middle
+                            Utils.volTetra(self.gridN, A, H, F, C) +  # middle
-                        Utils.volTetra(self.gridN, B, C, D, G) +  # cutted edge bottom
+                            Utils.volTetra(self.gridN, C, H, D, A) +  # cutted edge bottom
-                        Utils.volTetra(self.gridN, D, E, G, H))   # cutted edge bottom
+                            Utils.volTetra(self.gridN, C, G, H, F))   # cutted edge bottom
-                vol2 = (Utils.volTetra(self.gridN, A, F, B, C) +  # cutted edge top
+                    self._vol = (vol1 + vol2)/2
-                        Utils.volTetra(self.gridN, A, E, F, H) +  # cutted edge top
+            return self._vol
-                        Utils.volTetra(self.gridN, A, H, F, C) +  # middle
+        return locals()
-                        Utils.volTetra(self.gridN, C, H, D, A) +  # cutted edge bottom
+    _vol = None
-                        Utils.volTetra(self.gridN, C, G, H, F))   # cutted edge bottom
+    vol = property(**vol())
-                self._vol = (vol1 + vol2)/2
+    def area():
-        return self._vol
+        doc = "Face areas."
-    @property
+        def fget(self):
-    def area(self):
+            if(self._area is None or self._normals is None):
-        if (getattr(self, '_area', None) is None or
+                # Compute areas of cell faces
-            getattr(self, '_normals', None) is None):
+                if(self.dim == 2):
-            # Compute areas of cell faces
+                    xy = self.gridN
-            if(self.dim == 2):
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
-                xy = self.gridN
+                    edge1 = xy[B, :] - xy[A, :]
-                A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
+                    normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
-                                       self.nCy]))
+                    area1 = length2D(edge1)
-                edge1 = xy[B, :] - xy[A, :]
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
-                normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
+                    # Note that we are doing A-D to make sure the normal points the right way.
-                area1 = length2D(edge1)
+                    # Think about it. Look at the picture. Normal points towards C iff you do this.
-                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
+                    edge2 = xy[A, :] - xy[D, :]
-                                       self.nNy]))
+                    normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
-                # Note that we are doing A-D to make sure the normal points the
+                    area2 = length2D(edge2)
-                # right way.
+                    self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
-                # Think about it. Look at the picture. Normal points towards C
+                    self._normals = [normalize2D(normal1), normalize2D(normal2)]
-                # iff you do this.
+                elif(self.dim == 3):
                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)]
-            elif(self.dim == 3):
+                    A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz]))
                    normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
-                A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array(
+                    A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz]))
-                                             [self.nNx, self.nCy, self.nCz]))
+                    normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
                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(
+                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz]))
-                                             [self.nCx, self.nNy, self.nCz]))
+                    normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
                normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E,
                                                average=False,
                                                normalizeNormals=False)
-                A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array(
+                    self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
-                                             [self.nCx, self.nCy, self.nNz]))
+                    self._normals = [normal1, normal2, normal3]
-                normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D,
+            return self._area
-                                                average=False,
+        return locals()
-                                                normalizeNormals=False)
+    _area = None
    area = property(**area())
-                self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2),
+    def normals():
-                                   Utils.mkvc(area3)]
+        doc = """Face normals: calling this will average
                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
@@ -302,64 +276,58 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
            NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
        """
-        if getattr(self, '_normals', None) is None:
+        def fget(self):
-            self.area  # calling .area will create the face normals
+            if(self._normals is None):
-        if self.dim == 2:
+                self.area  # calling .area will create the face normals
-            return normalize2D(np.r_[self._normals[0], self._normals[1]])
+            if self.dim == 2:
-        elif self.dim == 3:
+                return normalize2D(np.r_[self._normals[0], self._normals[1]])
-            normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
+            elif self.dim == 3:
-            normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
+                normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
-            normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
+                normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
-            return normalize3D(np.r_[normal1, normal2, normal3])
+                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 locals()
    _normals = None
    normals = property(**normals())
-    @property
+    def edge():
-    def edge(self):
+        doc = "Edge legnths."
-        """
+
-        Edge lengths
+        def fget(self):
-        """
+            if(self._edge is None or self._tangents is None):
-        if getattr(self, '_edge', None) is None:
+                if(self.dim == 2):
-            if(self.dim == 2):
+                    xy = self.gridN
-                xy = self.gridN
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
-                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
+                    edge1 = xy[D, :] - xy[A, :]
-                                                                  self.nNy]))
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
-                edge1 = xy[D, :] - xy[A, :]
+                    edge2 = xy[B, :] - xy[A, :]
-                A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
+                    self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
-                                                                   self.nCy]))
+                    self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
-                edge2 = xy[B, :] - xy[A, :]
+                elif(self.dim == 3):
-                self._edge = np.r_[Utils.mkvc(length2D(edge1)),
+                    xyz = self.gridN
-                                   Utils.mkvc(length2D(edge2))]
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
-                self._tangents = np.r_[edge1, edge2]/np.c_[self._edge,
+                    edge1 = xyz[D, :] - xyz[A, :]
-                                                           self._edge]
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz]))
-            elif(self.dim == 3):
+                    edge2 = xyz[B, :] - xyz[A, :]
-                xyz = self.gridN
+                    A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
-                A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
+                    edge3 = xyz[E, :] - xyz[A, :]
-                                                                   self.nNy,
+                    self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
-                                                                   self.nNz]))
+                    self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
                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 self._edge
+        return locals()
    _edge = None
    edge = property(**edge())
-    @property
+    def tangents():
-    def tangents(self):
+        doc = "Edge tangents."
-        """
+
-        Edge tangents
+        def fget(self):
-        """
+            if(self._tangents is None):
-        if getattr(self, '_tangents', None) is None:
+                self.edge  # calling .edge will create the tangents
-            self.edge  # calling .edge will create the tangents
+            return self._tangents
-        return self._tangents
+        return locals()
    _tangents = None
    tangents = property(**tangents())
@@ -18,15 +18,13 @@ 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,"
+        assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
                                                  "'dirichlet' or 'neumann'")
    return bc
 def ddxCellGrad(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::
@@ -54,8 +52,7 @@ def ddxCellGrad(n, bc):
    """
    bc = checkBC(bc)
-    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n,
+    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
                   format="csr")
    # Set the first side
    if(bc[0] == 'dirichlet'):
        D[0, 0] = 2
@@ -68,11 +65,10 @@ 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
+
-        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::
@@ -103,7 +99,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
@@ -116,7 +112,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
@@ -125,166 +121,175 @@ class DiffOperators(object):
        Class creates the differential operators that you need!
    """
    def __init__(self):
-        raise Exception('DiffOperators is a base class providing differential'
+        raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
                        'operators on meshes and cannot run on its own.'
                        'Inherit to your favorite Mesh class.')
-    @property
+    def faceDiv():
-    def faceDiv(self):
+        doc = "Construct divergence operator (face-stg to cell-centres)."
-        """
+
-        Construct divergence operator (face-stg to cell-centres).
+        def fget(self):
-        """
+            if(self._faceDiv is None):
-        if getattr(self, '_faceDiv', None) is None:
+                # The number of cell centers in each direction
-            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]))
-                D = sp.hstack((D1, D2, D3), format="csr")
+                # Compute areas of cell faces & volumes
-            # Compute areas of cell faces & volumes
+                S = self.r(self.area, 'F', 'Fz', 'V')
-            S = self.area
+                V = self.vol
-            V = self.vol
+                self._faceDivz = sdiag(1/V)*D3*sdiag(S)
            self._faceDiv = sdiag(1/V)*D*sdiag(S)
        return self._faceDiv
-    @property
+            return self._faceDivz
-    def faceDivx(self):
+        return locals()
-        """
+    _faceDivz = None
-        Construct divergence operator in the x component (face-stg to
+    faceDivz = property(**faceDivz())
        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)
-        return self._faceDivx
+    def nodalGrad():
        doc = "Construct gradient operator (nodes to edges)."
-    @property
+        def fget(self):
-    def faceDivy(self):
+            if(self._nodalGrad is None):
-        if(self.dim < 2):
+                # The number of cell centers in each direction
-            return None
+                n = self.vnC
-        if getattr(self, '_faceDivy', None) is None:
+                # Compute divergence operator on faces
-            # The number of cell centers in each direction
+                if(self.dim == 1):
-            n = self.vnC
+                    G = ddx(n[0])
-            # Compute faceDivergence operator on faces
+                elif(self.dim == 2):
-            if(self.dim == 2):
+                    D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
-                D2 = sp.kron(ddx(n[1]), speye(n[0]))
+                    D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
-            elif(self.dim == 3):
+                    G = sp.vstack((D1, D2), format="csr")
-                D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
+                elif(self.dim == 3):
-            # Compute areas of cell faces & volumes
+                    D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
-            S = self.r(self.area, 'F', 'Fy', 'V')
+                    D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
-            V = self.vol
+                    D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
-            self._faceDivy = sdiag(1/V)*D2*sdiag(S)
+                    G = sp.vstack((D1, D2, D3), format="csr")
-        return self._faceDivy
+                # Compute lengths of cell edges
                L = self.edge
                self._nodalGrad = sdiag(1/L)*G
            return self._nodalGrad
        return locals()
    _nodalGrad = None
    nodalGrad = property(**nodalGrad())
-    @property
+    def nodalLaplacian():
-    def faceDivz(self):
+        doc = "Construct laplacian operator (nodes to edges)."
        """
        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
-    @property
+        def fget(self):
-    def nodalGrad(self):
+            if(self._nodalLaplacian is None):
-        """
+                print 'Warning: Laplacian has not been tested rigorously.'
-        Construct gradient operator (nodes to edges).
+                # The number of cell centers in each direction
-        """
+                n = self.vnC
-        if getattr(self, '_nodalGrad', None) is None:
+                # Compute divergence operator on faces
-            # The number of cell centers in each direction
+                if(self.dim == 1):
-            n = self.vnC
+                    D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
-            # Compute divergence operator on faces
+                    L  = - D1.T*D1
-            if(self.dim == 1):
+                elif(self.dim == 2):
-                G = ddx(n[0])
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
-            elif(self.dim == 2):
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
-                D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
+                    L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
-                D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
+                    L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
-                G = sp.vstack((D1, D2), format="csr")
+                    L  = L1 + L2
-            elif(self.dim == 3):
+                elif(self.dim == 3):
-                D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
-                D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
-                D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
+                    D3 = sdiag(1./self.hz) * ddx(n[2])
-                G = sp.vstack((D1, D2, D3), format="csr")
+                    L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
-            # Compute lengths of cell edges
+                    L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
-            L = self.edge
+                    L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
-            self._nodalGrad = sdiag(1/L)*G
+                    L  = L1 + L2 + L3
-        return self._nodalGrad
+                self._nodalLaplacian = L
-
+            return self._nodalLaplacian
-    @property
+        return locals()
-    def nodalLaplacian(self):
+    _nodalLaplacian = None
-        """
+    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
+        Function that sets the boundary conditions for cell-centred derivative operators.
        operators.
        Examples::
            # Neumann in all directions
            BC = 'neumann'
-            # 3D, Dirichlet in y Neumann else
+            BC = 'neumann'                                            # Neumann in all directions
-            BC = ['neumann', 'dirichlet', 'neumann']
+            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
            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):
@@ -318,69 +323,47 @@ class DiffOperators(object):
            G = sp.vstack((G1, G2, G3), format="csr")
        return G
-    @property
+    def cellGrad():
-    def cellGrad(self):
+        doc = "The cell centered Gradient, takes you to cell faces."
        """
        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
-    @property
+        def fget(self):
-    def cellGradBC(self):
+            if(self._cellGrad is None):
-        """
+                G = self._cellGradStencil()
-        The cell centered Gradient boundary condition matrix
+                # Compute areas of cell faces & volumes
-        """
+                S = self.area
-        if getattr(self, '_cellGradBC', None) is None:
+                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
-            BC = self.setCellGradBC(self._cellGradBC_list)
+                self._cellGrad = sdiag(S/V)*G
-            n = self.vnC
+            return self._cellGrad
-            if(self.dim == 1):
+        return locals()
-                G = ddxCellGradBC(n[0], BC[0])
+    _cellGrad = None
-            elif(self.dim == 2):
+    cellGrad = property(**cellGrad())
                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']
@@ -393,19 +376,20 @@ class DiffOperators(object):
            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
        return G1
-    @property
+
-    def cellGradx(self):
+    def cellGradx():
-        """
+        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
-        Cell centered Gradient in the x dimension. Has neumann boundary
+
-        conditions.
+        def fget(self):
-        """
+            if getattr(self, '_cellGradx', None) is None:
-        if getattr(self, '_cellGradx', None) is None:
+                G1 = self._cellGradxStencil()
-            G1 = self._cellGradxStencil()
+                # 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','Fx', 'V')
-            L = self.r(self.area/V, 'F','Fx', 'V')
+                self._cellGradx = sdiag(L)*G1
-            self._cellGradx = sdiag(L)*G1
+            return self._cellGradx
-        return self._cellGradx
+        return locals()
    cellGradx = property(**cellGradx())
    def _cellGradyStencil(self):
        if self.dim < 2: return None
@@ -417,17 +401,19 @@ class DiffOperators(object):
            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
        return G2
-    @property
+    def cellGrady():
-    def cellGrady(self):
+        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
-        if self.dim < 2:
+        def fget(self):
-            return None
+            if self.dim < 2: 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
@@ -436,61 +422,66 @@ class DiffOperators(object):
        G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
        return G3
-    @property
+    def cellGradz():
-    def cellGradz(self):
+        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
-        """
+        def fget(self):
-        Cell centered Gradient in the x dimension. Has neumann boundary
+            if self.dim < 3: return None
-        conditions.
+            if getattr(self, '_cellGradz', None) is None:
-        """
+                G3 = self._cellGradzStencil()
-        if self.dim < 3:
+                # Compute areas of cell faces & volumes
-            return None
+                V = self.aveCC2F*self.vol
-        if getattr(self, '_cellGradz', None) is None:
+                L = self.r(self.area/V, 'F','Fz', 'V')
-            G3 = self._cellGradzStencil()
+                self._cellGradz = sdiag(L)*G3
-            # Compute areas of cell faces & volumes
+            return self._cellGradz
-            V = self.aveCC2F*self.vol
+        return locals()
-            L = self.r(self.area/V, 'F', 'Fz', 'V')
+    cellGradz = property(**cellGradz())
            self._cellGradz = sdiag(L)*G3
        return self._cellGradz
-    @property
+    def edgeCurl():
-    def edgeCurl(self):
+        doc = "Construct the 3D curl operator."
        """
        Construct the 3D curl operator.
        """
        if getattr(self, '_edgeCurl', None) is None:
            assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
-            n = self.vnC  # The number of cell centers in each direction
+        def fget(self):
-            L = self.edge  # Compute lengths of cell edges
+            if(self._edgeCurl is None):
-            S = self.area # Compute areas of cell faces
+                assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
                # The number of cell centers in each direction
                n = self.vnC
-            # Compute divergence operator on faces
+                # Compute lengths of cell edges
-            if self.dim == 2:
+                L = self.edge
-                D21 = sp.kron(ddx(n[1]), speye(n[0]))
+                # Compute areas of cell faces
-                D12 = sp.kron(speye(n[1]), ddx(n[0]))
+                S = self.area
                C = sp.hstack((-D21, D12), format="csr")
                self._edgeCurl = C*sdiag(1/S)
-            elif self.dim == 3:
+                # Compute divergence operator on faces
                if self.dim == 2:
-                D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
+                    D21 = sp.kron(ddx(n[1]), speye(n[0]))
-                D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
+                    D12 = sp.kron(speye(n[1]), ddx(n[0]))
-                D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
+                    C = sp.hstack((-D21, D12), format="csr")
-                D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0]))
+                    self._edgeCurl = C*sdiag(1/S)
                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])
+                elif self.dim == 3:
                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
-                C = sp.vstack((sp.hstack((O1, -D32, D23)),
+                    D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
-                               sp.hstack((D31, O2, -D13)),
+                    D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
-                               sp.hstack((-D21, D12, O3))), format="csr")
+                    D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
                    D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0]))
                    D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]))
                    D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0]))
-                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
+                    O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
-        return self._edgeCurl
+                    O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
                    O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
                    C = sp.vstack((sp.hstack((O1, -D32, D23)),
                                   sp.hstack((D31, O2, -D13)),
                                   sp.hstack((-D21, D12, O3))), format="csr")
                    self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
            return self._edgeCurl
        return locals()
    _edgeCurl = None
    edgeCurl = property(**edgeCurl())
    def getBCProjWF(self, BC, discretization='CC'):
        """
@@ -498,19 +489,16 @@ class DiffOperators(object):
        The weak form boundary condition projection matrices.
        Examples::
            # Neumann in all directions
            BC = 'neumann'
-            # 3D, Dirichlet in y Neumann else
+            BC = 'neumann'                                            # Neumann in all directions
-            BC = ['neumann', 'dirichlet', 'neumann']
+            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
            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'
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
                                      'for CC discretization.')
        if(type(BC) is str):
            BC = [BC for _ in self.vnC]  # Repeat the str self.dim times
@@ -522,34 +510,35 @@ 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
@@ -561,7 +550,6 @@ 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]))
@@ -576,14 +564,12 @@ 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 = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
                          indF[5])]
            Pbc = Pbc*sdiag(self.area[indF])
            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
@@ -600,36 +586,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'
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
                                      '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'],
+        BC = [['dirichlet','dirichlet'],['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])
@@ -667,11 +653,9 @@ class DiffOperators(object):
        if(self.dim == 1):
            return self.aveFx2CC
        elif(self.dim == 2):
-            return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC),
+            return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
                                    format="csr")
        elif(self.dim == 3):
-            return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC,
+            return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
                                      self.aveFz2CC), format="csr")
    @property
    def aveF2CCV(self):
@@ -681,16 +665,11 @@ 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),
+            return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
                                 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):
@@ -703,12 +682,8 @@ class DiffOperators(object):
    @property
    def aveFy2CC(self):
-        """
+        "Construct the averaging operator on cell faces in the y direction to cell centers."
-        Construct the averaging operator on cell faces in the y direction to
+        if self.dim < 2: return None
        cell centers.
        """
        if self.dim < 2:
            return None
        if getattr(self, '_aveFy2CC', None) is None:
            n = self.vnC
            if(self.dim == 2):
@@ -719,10 +694,7 @@ 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
@@ -739,18 +711,12 @@ 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]),
+                self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
-                                           avExtrap(n[0])),
+                                           sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
                                           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]),
+                self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
-                                                 avExtrap(n[0])),
+                                           kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
-                                           kron3(speye(n[2]), avExtrap(n[1]),
+                                           kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
                                                 speye(n[0])),
                                           kron3(avExtrap(n[2]), speye(n[1]),
                                                 speye(n[0]))),
                                          format="csr")
        return self._aveCC2F
    @property
@@ -761,8 +727,7 @@ 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,
+            return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
                                     self.aveEz2CC), format="csr")
    @property
    def aveE2CCV(self):
@@ -772,15 +737,11 @@ 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),
+            return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
                                 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
@@ -794,12 +755,8 @@ class DiffOperators(object):
    @property
    def aveEy2CC(self):
-        """
+        "Construct the averaging operator on cell edges in the y direction to cell centers."
-        Construct the averaging operator on cell edges in the y direction to
+        if self.dim < 2: return None
        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
@@ -811,12 +768,8 @@ class DiffOperators(object):
    @property
    def aveEz2CC(self):
-        """
+        "Construct the averaging operator on cell edges in the z direction to cell centers."
-        Construct the averaging operator on cell edges in the z direction to
+        if self.dim < 3: return None
        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
@@ -840,10 +793,7 @@ 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
@@ -852,24 +802,16 @@ 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))),
+                                          sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
                                         format="csr")
            elif(self.dim == 3):
-                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1),
+                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
-                                                av(n[0])),
+                                          kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
-                                          kron3(speye(n[2]+1), av(n[1]),
+                                          kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
                                                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
@@ -877,14 +819,9 @@ 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]))),
+                                          sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
                                         format="csr")
            elif(self.dim == 3):
-                self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]),
+                self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
-                                                speye(n[0]+1)),
+                                          kron3(av(n[2]), speye(n[1]+1), av(n[0])),
-                                          kron3(av(n[2]), speye(n[1]+1),
+                                          kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
                                                av(n[0])),
                                          kron3(speye(n[2]+1), av(n[1]),
                                                av(n[0]))),
                                         format="csr")
        return self._aveN2F
@@ -421,7 +421,7 @@ class InnerProducts(object):
    def _getEdgePx(M):
        """Returns a function for creating projection matrices"""
        def Px(xEdge):
-            assert xEdge == 'eX0', 'xEdge = {0!s}, not eX0'.format(xEdge)
+            assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
            return sp.identity(M.nC)
        return Px
@@ -198,11 +198,11 @@ class TensorMeshIO(object):
        """
        assert mesh.dim == 3
        s = ''
-        s += '{0:d} {1:d} {2:d}\n'.format(*tuple(mesh.vnC))
+        s += '%i %i %i\n' %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 += '{0:.2f} {1:.2f} {2:.2f}\n'.format(*tuple(origin))
+        s += '%.2f %.2f %.2f\n' %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])
@@ -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[{0:d}] is not a numpy array.".format(i))
+            assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i)
-            assert len(h_i.shape) == 1, ("h[{0:d}] must be a 1D numpy array.".format(i))
+            assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % 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[{0:d}] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative.".format(i))
+                    raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i)
        if isinstance(self, BaseRectangularMesh):
            BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
@@ -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 {0!s} interpolation, as this variable does not exist.'.format(locType))
+            raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
        loc = Utils.asArray_N_x_Dim(loc, self.dim)
@@ -177,7 +177,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
        return l
    def __str__(self):
-        outStr = '  ---- {0!s}TreeMesh ----  '.format(('Oc' if self.dim == 3 else 'Quad'))
+        outStr = '  ---- %sTreeMesh ----  '%('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: {0:2.2f}%'.format((self.fill*100))
+        outStr += '\n  Fill: %2.2f%%'%(self.fill*100)
        return outStr
    @property
@@ -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 {0:d}, {1!s} = {2:4.2f}'.format(ind, normal, indLoc))
+        ax.set_title('Slice %d, %s = %4.2f' % (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 ({0:d}) is out of range.".format(key)
+                raise IndexError, "The index (%d) is out of range."%key
            self._numberCells() # no-op if numbered
            index   = self._i2cc[key]
@@ -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,
-                                     '#{0:.0f}'.format(iz),color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
+                                     '#%i'%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 ['{0!s}']".format("','".join(vTypeOpts))
+        assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
        assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.'
-        assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
+        assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
-        assert normal in normalOpts, "normal must be in ['{0!s}']".format("','".join(normalOpts))
+        assert normal in normalOpts, "normal must be in ['%s']" % "','".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 {0:.0f}'.format(ind))
+        ax.set_title('Slice %d' % 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 ['{0!s}'] when view='vec'".format("','".join(vTypeOptsV))
+            assert vType in vTypeOptsV, "vType must be in ['%s'] when view='vec'" % "','".join(vTypeOptsV)
-        assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
+        assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
        viewOpts = ['real','imag','abs','vec']
-        assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
+        assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
        if ax is None:
@@ -121,7 +121,7 @@ class Minimize(object):
    @callback.setter
    def callback(self, value):
        if self.callback is not None:
-            print 'The callback on the {0!s} Optimization was replaced.'.format(self.__name__)
+            print 'The callback on the %s Optimization was replaced.' % self.__name__
        self._callback = value
@@ -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: {0:e}\n'.format(norm(rt))
+                if self.comments: print '\t\tResid: %e\n'%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 ({0:d}). norm: {1:4.4e}'.format(self.maxIter, norm(rt))
+                print 'NewtonRoot stopped by maxIters (%d). norm: %4.4e' % (self.maxIter, norm(rt))
                break
        return x
@@ -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 {0!s} class.".format((self.surveyPair.__name__))
+        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
        if d.ispaired:
            raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
        self._survey = d
@@ -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: ("{0!s}", Maps.ReciprocalMap)'.format(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: ("%s", Maps.ReciprocalMap)'%self.name
        self._propertyLink = value
    def _getMapProperty(self):
        prop = self
        def fget(self):
-            return getattr(self, '_{0!s}Map'.format(prop.name), None)
+            return getattr(self, '_%sMap'%prop.name, None)
        def fset(self, val):
            if prop.propertyLink is not None:
                linkName, linkMap = prop.propertyLink
-                assert getattr(self, '{0!s}Map'.format(linkName), None) is None, 'Cannot set both sides of a linked property.'
+                assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.'
            # TODO: Check if the mapping can be correct
-            setattr(self, '_{0!s}Map'.format(prop.name), val)
+            setattr(self, '_%sMap'%prop.name, val)
        return property(fget=fget, fset=fset, doc=prop.doc)
    def _getIndexProperty(self):
        prop = self
        def fget(self):
-            return getattr(self, '_{0!s}Index'.format(prop.name), slice(None))
+            return getattr(self, '_%sIndex'%prop.name, slice(None))
        def fset(self, val):
-            setattr(self, '_{0!s}Index'.format(prop.name), val)
+            setattr(self, '_%sIndex'%prop.name, val)
        return property(fget=fget, fset=fset, doc=prop.doc)
    def _getProperty(self):
        prop = self
        def fget(self):
-            mapping = getattr(self, '{0!s}Map'.format(prop.name))
+            mapping = getattr(self, '%sMap'%prop.name)
            if mapping is None and prop.propertyLink is None:
                return prop.defaultVal
            if mapping is None and prop.propertyLink is not None:
                linkName, linkMapClass = prop.propertyLink
                linkMap = linkMapClass(None)
-                if getattr(self, '{0!s}Map'.format(linkName), None) is None:
+                if getattr(self, '%sMap'%linkName, None) is None:
                    return prop.defaultVal
-                m = getattr(self, '{0!s}'.format(linkName))
+                m = getattr(self, '%s'%linkName)
                return linkMap * m
-            m = getattr(self, '{0!s}Model'.format(prop.name))
+            m = getattr(self, '%sModel'%prop.name)
            return mapping * m
        return property(fget=fget)
    def _getModelDerivProperty(self):
        prop = self
        def fget(self):
-            mapping = getattr(self, '{0!s}Map'.format(prop.name))
+            mapping = getattr(self, '%sMap'%prop.name)
            if mapping is None and prop.propertyLink is None:
                return None
            if mapping is None and prop.propertyLink is not None:
                linkName, linkMapClass = prop.propertyLink
-                linkedMap = getattr(self, '{0!s}Map'.format(linkName))
+                linkedMap = getattr(self, '%sMap'%linkName)
                if linkedMap is None:
                    return None
                linkMap = linkMapClass(None) * linkedMap
-                m = getattr(self, '{0!s}Model'.format(linkName))
+                m = getattr(self, '%sModel'%linkName)
                return linkMap.deriv( m )
-            m = getattr(self, '{0!s}Model'.format(prop.name))
+            m = getattr(self, '%sModel'%prop.name)
            return mapping.deriv( m )
        return property(fget=fget)
    def _getModelProperty(self):
        prop = self
        def fget(self):
-            mapping = getattr(self, '{0!s}Map'.format(prop.name))
+            mapping = getattr(self, '%sMap'%prop.name)
            if mapping is None:
                return None
-            index = getattr(self.propMap, '{0!s}Index'.format(prop.name))
+            index = getattr(self.propMap, '%sIndex'%prop.name)
            return self.vector[index]
        return property(fget=fget)
    def _getModelProjProperty(self):
        prop = self
        def fget(self):
-            mapping = getattr(self, '{0!s}Map'.format(prop.name))
+            mapping = getattr(self, '%sMap'%prop.name)
            if mapping is None:
                return None
-            inds = getattr(self.propMap, '{0!s}Index'.format(prop.name))
+            inds = getattr(self.propMap, '%sIndex'%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, '_{0!s}Map'.format(prop.name), None)
+            return getattr(self.propMap, '_%sMap'%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, '{0!s}Index'.format(name), None)
+                index = getattr(self.propMap, '%sIndex'%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: '{0!s}' is trying to link to '{1!s}'".format(prop.name, prop.propertyLink[0])
+                assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0])
        if len(defaultInvProps) > 1:
-            raise Exception('You have more than one default inversion property: {0!s}'.format(defaultInvProps))
+            raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
        newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
@@ -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: [{0!s}]".format((', '.join(["('{0!s}', {1!s}Map)".format(p, p) for p in self._properties])))
+                for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(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, '{0!s}Map'.format(name), mapping)
+            setattr(self, '%sMap'%name, mapping)
-            setattr(self, '{0!s}Index'.format(name), slices.get(name, slice(nP, nP + mapping.nP)))
+            setattr(self, '%sIndex'%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, '{0!s}Map'.format(name), None)
+            setattr(self, '%sMap'%name, None)
-            setattr(self, '{0!s}Index'.format(name), None)
+            setattr(self, '%sIndex'%name, None)
    def __call__(self, vec):
        return self.PropModel(self, vec)
    def __contains__(self, val):
-        activeMaps = [name for name in self._properties if getattr(self, '{0!s}Map'.format(name)) is not None]
+        activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None]
        return val in activeMaps
@@ -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 ['{0!s}']".format(("', '".join(known)))
+            assert value in known, "rxType must be in ['%s']" % ("', '".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 {0!s}'.format(self.rxPair.__name__)
+            assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
        assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
        self.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 {0!s}'.format(self.srcPair.__name__)
+        assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
        assert len(set(value)) == len(value), 'The srcList must be unique'
        self._srcList = value
        self._sourceOrder = dict()
@@ -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: {0!s}'.format(str(src)))
+                raise KeyError('Source does not have a uid: %s'%str(src))
        inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
        if None in inds:
-            raise KeyError('Some of the sources specified are not in this survey. {0!s}'.format(str(inds)))
+            raise KeyError('Some of the sources specified are not in this survey. %s'%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 {0!s} class.".format((p.surveyPair.__name__))
+        assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__)
        if p.ispaired:
            raise Exception("The problem object is already paired to a survey. Use prob.unpair()")
        self._prob = p
@@ -4,6 +4,7 @@ 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
@@ -199,10 +200,10 @@ class OrderTest(unittest.TestCase):
                    print '_____________________________________________'
                    print '   h  |    error    | e(i-1)/e(i) |  order'
                    print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
-                    print '{0:4d}  |  {1:8.2e}   |'.format(nc, err)
+                    print '%4i  |  %8.2e   |' % (nc, err)
                else:
                    order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
-                    print '{0:4d}  |  {1:8.2e}   |   {2:6.4f}    |  {3:6.4f}'.format(nc, err, err_old/err, order[-1])
+                    print '%4i  |  %8.2e   |   %6.4f    |  %6.4f' % (nc, err, err_old/err, order[-1])
                err_old = err
                max_h_old = max_h
            print '---------------------------------------------'
@@ -257,8 +258,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
            Tests.checkDerivative(simplePass, np.random.randn(5))
    """
-    print "{0!s} checkDerivative {1!s}".format('='*20, '='*20)
+    print "%s checkDerivative %s" % ('='*20, '='*20)
-    print "iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order\n{0!s}".format(('-'*57))
+    print "iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order\n%s" % ('-'*57)
    f0, J0 = fctn(x0)
@@ -289,7 +290,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 " {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])
+        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])
    # Ensure we are about precision
    order0 = order0[E0[1:] > eps]
@@ -301,10 +302,10 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
    passTest = belowTol or correctOrder
    if passTest:
-        print "{0!s} PASS! {1!s}".format('='*25, '='*25)
+        print "%s PASS! %s" % ('='*25, '='*25)
        print happiness[np.random.randint(len(happiness))]+'\n'
    else:
-        print "{0!s}\n{1!s} FAIL! {2!s}\n{3!s}".format('*'*57, '<'*25, '>'*25, '*'*57)
+        print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57)
        print sadness[np.random.randint(len(sadness))]+'\n'
@@ -313,7 +314,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 - {0!s}'.format(('PASSED :)' if passTest else 'FAILED :(')))
+        ax.set_title('Check Derivative - %s' % ('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',
@@ -8,7 +8,7 @@ def _checkAccuracy(A, b, X, accuracyTol):
    if nrm_b > 0:
        nrm /= nrm_b
    if nrm > accuracyTol:
-        msg = '### SolverWarning ###: Accuracy on solve is above tolerance: {0:e} > {1:e}'.format(nrm, accuracyTol)
+        msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
        print msg
        warnings.warn(msg, RuntimeWarning)
@@ -32,7 +32,7 @@ def memProfileWrapper(towrap, *funNames):
        if hasattr(towrap,f):
            attrs[f] = profile(getattr(towrap,f))
        else:
-            print '{0!s} not found in {1!s} Class'.format(f, towrap.__name__)
+            print '%s not found in %s Class' % (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('{0!s} attr is not recognized'.format(attr))
+            raise Exception('%s attr is not recognized' % 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 += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['title']) + ''
+        titles += ('{:^%i}'%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 += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['format'] % printer['value'](obj))
+        values += ('{:^%i}'%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 + "{0!s}{1!s}{2!s}".format('-'*25, stop, '-'*25)
+    print pad + "%s%s%s" % ('-'*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 + "{0!s}{1!s}{2!s}".format('-'*25, done, '-'*25)
+    print pad + "%s%s%s" % ('-'*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 {0!s}, you can create a method::
+            If you have things that also need to run in the method %s, you can create a method::
-                def _{1!s}*(self, ... ):
+                def _%s*(self, ... ):
                    pass
-            Where the * can be any string. If present, _{2!s}* will be called at the start of the default {3!s} call.
+            Where the * can be any string. If present, _%s* will be called at the start of the default %s call.
            You may also completely overwrite this function.
-        """.format(match, match, match, match)
+        """ % (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, {0:d}) not ({1:d}, {2:d})".format(*((dim,)+pts.shape))
+        assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, %d) not (%d, %d)" % ((dim,)+pts.shape)
        return pts
@@ -207,17 +207,17 @@ def requires(var):
        .. note::
-            To use survey.{0!s}(), SimPEG requires that a problem be bound to the survey.
+            To use survey.%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)
-            """.format(f.__name__)
+            """ % f.__name__
        else:
            extra = """
-                To use *{0!s}* method, SimPEG requires that the {1!s} be specified.
+                To use *%s* method, SimPEG requires that the %s be specified.
-            """.format(f.__name__, var)
+            """ % (f.__name__, var)
        @wraps(f)
        def requiresVarWrapper(self,*args,**kwargs):
            if getattr(self, var, None) is None:
@@ -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: '{0!s}'".format(possibleNodes)
+        assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
    dim = gridSize.size
    if n is None:
        n = gridSize - 1
@@ -278,7 +278,7 @@ class TensorType(object):
        else:
            raise Exception('Unexpected shape of tensor')
    def __str__(self):
-        return 'TensorType[{0:d}]: {1!s}'.format(self._tt, self._tts)
+        return 'TensorType[%i]: %s' % (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
@@ -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 {0!s} mesh'.format(mesh._meshType))
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%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 {0!s} mesh'.format(mesh._meshType))
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%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')
@@ -15,7 +15,7 @@ import Directives
 import Inversion
 import Tests
-__version__   = '0.1.12'
+__version__   = '0.1.11'
 __author__    = 'Rowan Cockett'
 __license__   = 'MIT'
 __copyright__ = 'Copyright 2014 Rowan Cockett'
@@ -51,9 +51,9 @@ copyright = u'2013 - 2016, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.12'
+version = '0.1.11'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.12'
+release = '0.1.11'
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -266,7 +266,7 @@ 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)))
+        self._warnfunc(msg, '%s:%s' % get_source_line(node))
 supress_nonlocal_image_warn()
@@ -63,8 +63,26 @@ done by the :class:`SimPEG.Maps.ExpMap` described above.
 .. plot::
-    from SimPEG import Examples
+    from SimPEG import *
-    Examples.Maps_ComboMaps.run()
+    import matplotlib.pyplot as plt
    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
    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).
@@ -149,10 +167,31 @@ Map 2D Cross-Section to 3D Model
 Mesh to Mesh Map
 ----------------
 .. plot::
-    from SimPEG import Examples
+    from SimPEG import *
-    Examples.Maps_Mesh2Mesh.run()
+    import matplotlib.pyplot as plt
    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
@@ -66,7 +66,7 @@ Numpy and Matlab
 Lessons in Python
 -----------------
-* `Software Carpentry <http://swcarpentry.github.io/python-novice-inflammation/>`_
+* `Software Carpentry <http://software-carpentry.org/v4/python/index.html>`_
 * `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_
 Editing Python
@@ -40,23 +40,14 @@
 Direct Current Resistivity
 **************************
-`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
+`SimPEG.EM.Static.DC` and `SimPEG.EM.Static.IP` 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 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/>`_).
 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
@@ -159,7 +150,7 @@ Comparing to the analytic function:
 API for DC codes
 ================
-.. automodule:: SimPEG.DCIP.BaseDC
+.. automodule:: SimPEG.EM.Static.DC
    :show-inheritance:
    :members:
    :undoc-members:
@@ -1,26 +0,0 @@
 .. _examples_Inversion_IRLS:
 .. --------------------------------- ..
 ..                                   ..
 ..    THIS FILE IS AUTO GENEREATED   ..
 ..                                   ..
 ..    SimPEG/Examples/__init__.py    ..
 ..                                   ..
 .. --------------------------------- ..
 Inversion: Linear Problem
 =========================
 Here we go over the basics of creating a linear problem and inversion.
 .. plot::
    from SimPEG import Examples
    Examples.Inversion_IRLS.run()
 .. literalinclude:: ../../../SimPEG/Examples/Inversion_IRLS.py
    :language: python
    :linenos:
@@ -1,48 +0,0 @@
 .. _examples_Maps_ComboMaps:
 .. --------------------------------- ..
 ..                                   ..
 ..    THIS FILE IS AUTO GENEREATED   ..
 ..                                   ..
 ..    SimPEG/Examples/__init__.py    ..
 ..                                   ..
 .. --------------------------------- ..
 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.
 .. plot::
    from SimPEG import Examples
    Examples.Maps_ComboMaps.run()
 .. literalinclude:: ../../../SimPEG/Examples/Maps_ComboMaps.py
    :language: python
    :linenos:
@@ -1,27 +0,0 @@
 .. _examples_Maps_Mesh2Mesh:
 .. --------------------------------- ..
 ..                                   ..
 ..    THIS FILE IS AUTO GENEREATED   ..
 ..                                   ..
 ..    SimPEG/Examples/__init__.py    ..
 ..                                   ..
 .. --------------------------------- ..
 Maps: Mesh2Mesh
 ===============
 This mapping allows you to go from one mesh to another.
 .. plot::
    from SimPEG import Examples
    Examples.Maps_Mesh2Mesh.run()
 .. literalinclude:: ../../../SimPEG/Examples/Maps_Mesh2Mesh.py
    :language: python
    :linenos:
@@ -7,7 +7,7 @@ Todo: docs for IP!
 API for IP codes
 ================
-.. automodule:: SimPEG.DCIP.BaseIP
+.. automodule:: SimPEG.EM.Static.IP
    :show-inheritance:
    :members:
    :undoc-members:
@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
    :alt: SimPEG Logo
 SimPEG Documentation
@@ -16,25 +16,13 @@ SimPEG Documentation
    :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
    :alt: BSD 3 clause license.
-.. image:: https://api.travis-ci.org/simpeg/simpeg.svg?branch=master
+.. image:: https://img.shields.io/travis/simpeg/simpeg.svg
-    :target: https://travis-ci.org/simpeg/simpeg
+    :target: https://travis-ci.org/simpeg/simpeg?branch=master
    :alt: Travis CI build status
-.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
+.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
-    :alt: gitter chat room at https://gitter.im/simpeg/simpeg
+    :target: https://coveralls.io/r/simpeg/simpeg?branch=master
-    :target: https://gitter.im/simpeg/simpeg
+    :alt: Coverage status
 .. 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&amp;utm_medium=referral&amp;utm_content=simpeg/simpeg&amp;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.
@@ -54,7 +54,7 @@ class Images(webapp2.RequestHandler):
 class Redirect(webapp2.RequestHandler):
    def get(self):
        path = str(self.request.path).split(os.path.sep)[3:]
-        self.redirect(('/{0!s}'.format(os.path.sep.join(path))), permanent=True)
+        self.redirect(('/%s'%os.path.sep.join(path)), permanent=True)
 class MainPage(webapp2.RequestHandler):
@@ -83,7 +83,7 @@ with open("README.rst") as f:
 setup(
    name = "SimPEG",
-    version = "0.1.12",
+    version = "0.1.11",
    packages = find_packages(),
    install_requires = ['numpy>=1.7',
                        'scipy>=0.13',
@@ -1,16 +1,14 @@
 import unittest
-from SimPEG import Mesh, Problem, Fields, Survey, Utils
+from SimPEG import *
 import numpy as np
 class FieldsTest(unittest.TestCase):
    def setUp(self):
-        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
-                               [0, 0, -30])
+        x = np.linspace(5,10,3)
-        x = np.linspace(5, 10, 3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        XYZ = Utils.ndgrid(x, x, np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
        srcLoc = np.r_[0., 0., 0.]
        rxList0 = Survey.BaseRx(XYZ, 'exi')
        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
        rxList1 = Survey.BaseRx(XYZ, 'bxi')
@@ -20,13 +18,10 @@ class FieldsTest(unittest.TestCase):
        rxList3 = Survey.BaseRx(XYZ, 'bxi')
        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
-        srcList = [Src0, Src1, Src2, Src3, Src4]
+        srcList = [Src0,Src1,Src2,Src3,Src4]
        survey = Survey.BaseSurvey(srcList=srcList)
        self.D = Survey.Data(survey)
-        self.F = Problem.Fields(mesh, survey, knownFields={'phi': 'CC',
+        self.F = Problem.Fields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}, dtype={"phi":float,"e":complex,"b":complex})
                                                           'e': 'E', 'b': 'F'},
                                dtype={"phi": float, "e": complex,
                                       "b": complex})
        self.Src0 = Src0
        self.Src1 = Src1
        self.mesh = mesh
@@ -43,94 +38,83 @@ class FieldsTest(unittest.TestCase):
        self.assertTrue('e' in F)
    def test_overlappingFields(self):
-        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh,
+        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
-                          self.F.survey, knownFields={'b': 'F'},
+                            knownFields={'b':'F'},
-                          aliasFields={'b': ['b', (lambda F, b, ind: b)]})
+                            aliasFields={'b':['b',(lambda F, b, ind: b)]})
    def test_SetGet(self):
        F = self.F
        nSrc = F.survey.nSrc
-        e = (np.random.rand(F.mesh.nE, nSrc) +
+        e = np.random.rand(F.mesh.nE, nSrc) + np.random.rand(F.mesh.nE, nSrc)*1j
             np.random.rand(F.mesh.nE, nSrc)*1j)
        F[:, 'e'] = e
-        b = (np.random.rand(F.mesh.nF, nSrc) +
+        b = np.random.rand(F.mesh.nF, nSrc) + np.random.rand(F.mesh.nF, nSrc)*1j
             np.random.rand(F.mesh.nF, nSrc)*1j)
        F[:, 'b'] = b
        self.assertTrue(np.all(F[:, 'e'] == e))
        self.assertTrue(np.all(F[:, 'b'] == b))
-        F[:] = {'b': b, 'e': e}
+        F[:] = {'b':b,'e':e}
        self.assertTrue(np.all(F[:, 'e'] == e))
        self.assertTrue(np.all(F[:, 'b'] == b))
-        for s in [0, 0.0, np.r_[0], long(0)]:
+        for s in [0,0.0,np.r_[0],long(0)]:
            F[:, 'b'] = s
            self.assertTrue(np.all(F[:, 'b'] == b*0))
-        b = np.random.rand(F.mesh.nF, 1)
+        b = np.random.rand(F.mesh.nF,1)
        F[self.Src0, 'b'] = b
        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
-        b = np.random.rand(F.mesh.nF, 1)
+        b = np.random.rand(F.mesh.nF,1)
        F[self.Src0, 'b'] = b
        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
-        phi = np.random.rand(F.mesh.nC, 2)
+        phi = np.random.rand(F.mesh.nC,2)
-        F[[self.Src0, self.Src1], 'phi'] = phi
+        F[[self.Src0,self.Src1], 'phi'] = phi
-        self.assertTrue(np.all(F[[self.Src0, self.Src1], 'phi'] == phi))
+        self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
-        fdict = F[:, :]
+        fdict = F[:,:]
        self.assertTrue(type(fdict) is dict)
-        self.assertTrue(sorted([k for k in fdict]) == ['b', 'e', 'phi'])
+        self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
        b = np.random.rand(F.mesh.nF, 2)
-        F[[self.Src0, self.Src1], 'b'] = b
+        F[[self.Src0, self.Src1],'b'] = b
-        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF, 1))
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1))
-        self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, 1))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1))
-        self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b[:, 0], 2)))
+        self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2)))
-        self.assertTrue(np.all(F[self.Src1, 'b'] == Utils.mkvc(b[:, 1], 2)))
+        self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2)))
    def test_assertions(self):
        freq = [self.Src0, self.Src1]
        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
        def fun(): self.F[freq, 'b'] = bWrongSize
        self.assertRaises(ValueError, fun)
        def fun(): self.F[-999.]
        self.assertRaises(KeyError, fun)
        def fun(): self.F['notRight']
        self.assertRaises(KeyError, fun)
-
+        def fun(): self.F[freq,'notThere']
        def fun(): self.F[freq, 'notThere']
        self.assertRaises(KeyError, fun)
 class FieldsTest_Alias(unittest.TestCase):
    def setUp(self):
-        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
-                               [0, 0, -30])
+        x = np.linspace(5,10,3)
-        x = np.linspace(5, 10, 3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        XYZ = Utils.ndgrid(x, x, np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
        srcLoc = np.r_[0, 0, 0.]
        rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
+        Src0 = Survey.BaseSrc([rxList0],loc=srcLoc)
        rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
+        Src1 = Survey.BaseSrc([rxList1],loc=srcLoc)
        rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
+        Src2 = Survey.BaseSrc([rxList2],loc=srcLoc)
        rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src3 = Survey.BaseSrc([rxList3],loc=srcLoc)
-        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
-        srcList = [Src0, Src1, Src2, Src3, Src4]
+        srcList = [Src0,Src1,Src2,Src3,Src4]
        survey = Survey.BaseSurvey(srcList=srcList)
-        self.F = Problem.Fields(mesh, survey, knownFields={'e': 'E'},
+        self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
                                aliasFields={'b': ['e', 'F',
                                                   (lambda e, ind:
                                                    self.F.mesh.edgeCurl *
                                                    e)]})
        self.Src0 = Src0
        self.Src1 = Src1
        self.mesh = mesh
@@ -151,35 +135,33 @@ class FieldsTest_Alias(unittest.TestCase):
        nSrc = F.survey.nSrc
        e = np.random.rand(F.mesh.nE, nSrc)
        F[:, 'e'] = e
-        self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e))
+        self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e ))
-        e = np.random.rand(F.mesh.nE, 1)
+        e = np.random.rand(F.mesh.nE,1)
        F[self.Src0, 'e'] = e
        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e))
        def f():
            F[self.Src0, 'b'] = F[self.Src0, 'b']
-        self.assertRaises(KeyError, f)  # can't set a alias attr.
+        self.assertRaises(KeyError, f) # can't set a alias attr.
    def test_aliasFunction(self):
        def alias(e, ind):
            self.assertTrue(ind[0] is self.Src0)
            return self.F.mesh.edgeCurl * e
-        F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
+        F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
-                           aliasFields={'b': ['e', 'F', alias]})
+        e = np.random.rand(F.mesh.nE,1)
        e = np.random.rand(F.mesh.nE, 1)
        F[self.Src0, 'e'] = e
        F[self.Src0, 'b']
        def alias(e, ind):
            self.assertTrue(type(ind) is list)
            self.assertTrue(ind[0] is self.Src0)
            self.assertTrue(ind[1] is self.Src1)
            return self.F.mesh.edgeCurl * e
-
+        F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
-        F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
+        e = np.random.rand(F.mesh.nE,2)
                           aliasFields={'b': ['e', 'F', alias]})
        e = np.random.rand(F.mesh.nE, 2)
        F[[self.Src0, self.Src1], 'e'] = e
        F[[self.Src0, self.Src1], 'b']
@@ -187,11 +169,10 @@ class FieldsTest_Alias(unittest.TestCase):
 class FieldsTest_Time(unittest.TestCase):
    def setUp(self):
-        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
-                               [0, 0, -30])
+        x = np.linspace(5,10,3)
-        x = np.linspace(5, 10, 3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        XYZ = Utils.ndgrid(x, x, np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
        srcLoc = np.r_[0, 0, 0.]
        rxList0 = Survey.BaseRx(XYZ, 'exi')
        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
        rxList1 = Survey.BaseRx(XYZ, 'bxi')
@@ -201,13 +182,11 @@ class FieldsTest_Time(unittest.TestCase):
        rxList3 = Survey.BaseRx(XYZ, 'bxi')
        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
-        srcList = [Src0, Src1, Src2, Src3, Src4]
+        srcList = [Src0,Src1,Src2,Src3,Src4]
        survey = Survey.BaseSurvey(srcList=srcList)
-        prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
+        prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
        survey.pair(prob)
-        self.F = Problem.TimeFields(mesh, survey, knownFields={'phi': 'CC',
+        self.F = Problem.TimeFields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'})
                                                               'e': 'E',
                                                               'b': 'F'})
        self.Src0 = Src0
        self.Src1 = Src1
        self.mesh = mesh
@@ -237,92 +216,81 @@ class FieldsTest_Time(unittest.TestCase):
        self.assertTrue(np.all(F[:, 'e'] == e))
        self.assertTrue(np.all(F[:, 'b'] == b))
-        F[:] = {'b': b, 'e': e}
+        F[:] = {'b':b,'e':e}
        self.assertTrue(np.all(F[:, 'e'] == e))
        self.assertTrue(np.all(F[:, 'b'] == b))
-        for s in [0, 0.0, np.r_[0], long(0)]:
+        for s in [0,0.0,np.r_[0],long(0)]:
            F[:, 'b'] = s
            self.assertTrue(np.all(F[:, 'b'] == b*0))
-        b = np.random.rand(F.mesh.nF, 1, nT)
+        b = np.random.rand(F.mesh.nF,1,nT)
        F[self.Src0, 'b'] = b
-        self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :]))
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b[:,0,:]))
-        b = np.random.rand(F.mesh.nF, 1, nT)
+        b = np.random.rand(F.mesh.nF,1,nT)
-        F[self.Src0, 'b', 0] = b[:, :, 0]
+        F[self.Src0, 'b', 0] = b[:,:,0]
-        self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:, 0, 0],
+        self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2)))
                                                                  2)))
-        phi = np.random.rand(F.mesh.nC, 2, nT)
+        phi = np.random.rand(F.mesh.nC,2,nT)
-        F[[self.Src0, self.Src1], 'phi'] = phi
+        F[[self.Src0,self.Src1], 'phi'] = phi
-        self.assertTrue(np.all(F[[self.Src0, self.Src1], 'phi'] == phi))
+        self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
        fdict = F[:]
        self.assertTrue(type(fdict) is dict)
-        self.assertTrue(sorted([k for k in fdict]) == ['b', 'e', 'phi'])
+        self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
        b = np.random.rand(F.mesh.nF, 2, nT)
-        F[[self.Src0, self.Src1], 'b'] = b
+        F[[self.Src0, self.Src1],'b'] = b
-        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF, nT))
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,nT))
-        self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, nT))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
-        self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :]))
+        self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
-        self.assertTrue(np.all(F[self.Src1, 'b'] == b[:, 1, :]))
+        self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
-        self.assertTrue(np.all(F[self.Src0, 'b', 1] ==
+        self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2)))
-                        Utils.mkvc(b[:, 0, 1], 2)))
+        self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2)))
-        self.assertTrue(np.all(F[self.Src1, 'b', 1] ==
+        self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2)))
-                        Utils.mkvc(b[:, 1, 1], 2)))
+        self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2)))
-        self.assertTrue(np.all(F[self.Src0, 'b', 4] ==
+
                        Utils.mkvc(b[:, 0, 4], 2)))
        self.assertTrue(np.all(F[self.Src1, 'b', 4] ==
                        Utils.mkvc(b[:, 1, 4], 2)))
        b = np.random.rand(F.mesh.nF, 2, nT)
-        F[[self.Src0, self.Src1], 'b', 0] = b[:, :, 0]
+        F[[self.Src0, self.Src1],'b', 0] = b[:,:,0]
    def test_assertions(self):
        freq = [self.Src0, self.Src1]
        bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
        def fun(): self.F[freq, 'b'] = bWrongSize
        self.assertRaises(ValueError, fun)
        def fun(): self.F[-999.]
        self.assertRaises(KeyError, fun)
        def fun(): self.F['notRight']
        self.assertRaises(KeyError, fun)
-
+        def fun(): self.F[freq,'notThere']
        def fun(): self.F[freq, 'notThere']
        self.assertRaises(KeyError, fun)
 class FieldsTest_Time_Aliased(unittest.TestCase):
    def setUp(self):
-        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
-                               [0, 0, -30])
+        x = np.linspace(5,10,3)
-        x = np.linspace(5, 10, 3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
-        XYZ = Utils.ndgrid(x, x, np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
        srcLoc = np.r_[0, 0, 0.]
        rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
+        Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc)
        rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
+        Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc)
        rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
+        Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc)
        rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src3 = Survey.BaseSrc( [rxList3],loc=srcLoc)
-        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc( [rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
-        srcList = [Src0, Src1, Src2, Src3, Src4]
+        srcList = [Src0,Src1,Src2,Src3,Src4]
        survey = Survey.BaseSurvey(srcList=srcList)
-        prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
+        prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
        survey.pair(prob)
        def alias(b, srcInd, timeInd):
            return self.F.mesh.edgeCurl.T * b + timeInd
-        self.F = Problem.TimeFields(mesh, survey, knownFields={'b': 'F'},
+        self.F = Problem.TimeFields(mesh, survey, knownFields={'b':'F'}, aliasFields={'e':['b','E',alias]})
                                    aliasFields={'e': ['b', 'E', alias]})
        self.Src0 = Src0
        self.Src1 = Src1
        self.mesh = mesh
@@ -339,58 +307,55 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
        self.assertTrue('e' in F)
        self.assertTrue('b' in F)
    def test_simpleAlias(self):
        F = self.F
        nSrc = F.survey.nSrc
        nT = F.survey.prob.nT + 1
        b = np.random.rand(F.mesh.nF, nSrc, nT)
        F[:, 'b', :] = b
-        self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:, :, 0]))
+        self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:,:,0] ))
        e = range(nT)
        for i in range(nT):
-            e[i] = F.mesh.edgeCurl.T*b[:, :, i] + i
+            e[i] = F.mesh.edgeCurl.T*b[:,:,i] + i
-            e[i] = e[i][:, :, np.newaxis]
+            e[i] = e[i][:,:,np.newaxis]
        e = np.concatenate(e, axis=2)
-        self.assertTrue(np.all(F[:, 'e', :] == e))
+        self.assertTrue(np.all(F[:, 'e', :] == e ))
-        self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:, 0, :]))
+        self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
-        self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:, 1, :]))
+        self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
        for t in range(nT):
-            self.assertTrue(np.all(F[self.Src1, 'e', t] ==
+            self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) ))
                            Utils.mkvc(e[:, 1, t], 2)))
-        b = np.random.rand(F.mesh.nF, nT)
+        b = np.random.rand(F.mesh.nF,nT)
-        F[self.Src0, 'b', :] = b
+        F[self.Src0, 'b',:] = b
        Cb = F.mesh.edgeCurl.T * b
        for i in range(Cb.shape[1]):
-            Cb[:, i] += i
+            Cb[:,i] += i
-        self.assertTrue(np.all(F[self.Src0, 'e', :] == Cb))
+        self.assertTrue(np.all(F[self.Src0, 'e',:] == Cb))
        def f():
            F[self.Src0, 'e'] = F[self.Src0, 'e']
-        self.assertRaises(KeyError, f)  # can't set a alias attr.
+        self.assertRaises(KeyError, f) # can't set a alias attr.
    def test_aliasFunction(self):
        nT = self.F.survey.prob.nT + 1
        count = [0]
        def alias(e, srcInd, timeInd):
            count[0] += 1
            self.assertTrue(srcInd[0] is self.Src0)
            return self.F.mesh.edgeCurl * e
-        F = Problem.TimeFields(self.F.mesh, self.F.survey,
+        F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
-                               knownFields={'e': 'E'},
+        e = np.random.rand(F.mesh.nE,1,nT)
                               aliasFields={'b': ['e', 'F', alias]})
        e = np.random.rand(F.mesh.nE, 1, nT)
        F[self.Src0, 'e', :] = e
        F[self.Src0, 'b', :]
-        # ensure that this is called for every time separately.
+        self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
-        self.assertTrue(count[0] == nT)
+        e = np.random.rand(F.mesh.nE,1,1)
        e = np.random.rand(F.mesh.nE, 1, 1)
        F[self.Src0, 'e', 1] = e
        count[0] = 0
        F[self.Src0, 'b', 1]
-        self.assertTrue(count[0] == 1)  # ensure that this is called only once.
+        self.assertTrue(count[0] == 1) # ensure that this is called only once.
        def alias(e, srcInd, timeInd):
            count[0] += 1
@@ -398,21 +363,17 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
            self.assertTrue(srcInd[0] is self.Src0)
            self.assertTrue(srcInd[1] is self.Src1)
            return self.F.mesh.edgeCurl * e
-        F = Problem.TimeFields(self.F.mesh, self.F.survey,
+        F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
-                               knownFields={'e': 'E'},
+        e = np.random.rand(F.mesh.nE,2, nT)
                               aliasFields={'b': ['e', 'F', alias]})
        e = np.random.rand(F.mesh.nE, 2, nT)
        F[[self.Src0, self.Src1], 'e', :] = e
        count[0] = 0
        F[[self.Src0, self.Src1], 'b', :]
-
+        self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
-        # ensure that this is called for every time separately.
+        e = np.random.rand(F.mesh.nE,2, 1)
        self.assertTrue(count[0] == nT)
        e = np.random.rand(F.mesh.nE, 2, 1)
        F[[self.Src0, self.Src1], 'e', 1] = e
        count[0] = 0
        F[[self.Src0, self.Src1], 'b', 1]
-        self.assertTrue(count[0] == 1)  # ensure that this is called only once.
+        self.assertTrue(count[0] == 1) # ensure that this is called only once.
 if __name__ == '__main__':
@@ -28,14 +28,14 @@ class RegularizationTests(unittest.TestCase):
                for i, mesh in enumerate(self.meshlist):
-                    print 'Testing {0:d}D'.format(mesh.dim)
+                    print 'Testing %iD'%mesh.dim
                    mapping = r.mapPair(mesh)
                    reg = r(mesh, mapping=mapping)
                    m = np.random.rand(mapping.nP)
                    reg.mref = np.ones_like(m)*np.mean(m)
-                    print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
                    passed = reg.eval(reg.mref) < TOL
                    self.assertTrue(passed)
@@ -56,7 +56,7 @@ class RegularizationTests(unittest.TestCase):
                for i, mesh in enumerate(self.meshlist):
-                    print 'Testing Active Cells {0:d}D'.format((mesh.dim))
+                    print 'Testing Active Cells %iD'%(mesh.dim)
                    if mesh.dim == 1:
                        indActive = Utils.mkvc(mesh.gridCC <= 0.8)
@@ -70,7 +70,7 @@ class RegularizationTests(unittest.TestCase):
                        m = np.random.rand(mesh.nC)[indAct]
                        reg.mref = np.ones_like(m)*np.mean(m)
-                    print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
                    passed = reg.eval(reg.mref) < TOL
                    self.assertTrue(passed)
@@ -87,7 +87,7 @@ class RegularizationTests(unittest.TestCase):
            for i, mesh in enumerate(self.meshlist):
-                print 'Testing {0:d}D'.format(mesh.dim)
+                print 'Testing %iD'%mesh.dim
                # mapping = r.mapPair(mesh)
                # reg = r(mesh, mapping=mapping)
@@ -1,12 +0,0 @@
 import os
 import glob
 import unittest
 if __name__ == '__main__':
    test_file_strings = glob.glob('test_*.py')
    module_strings = [str[0:len(str)-3] for str in test_file_strings]
    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
              in module_strings]
    testSuite = unittest.TestSuite(suites)
    unittest.TextTestRunner(verbosity=2).run(testSuite)
@@ -1,77 +0,0 @@
 import unittest
 from SimPEG import *
 import SimPEG.DCIP as DC
 class DCProblemTests(unittest.TestCase):
    def setUp(self):
        aSpacing=2.5
        nElecs=10
        surveySize = nElecs*aSpacing - aSpacing
        cs = surveySize/nElecs/4
        mesh = Mesh.TensorMesh([
                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
                [(cs,3, -1.3),(cs,3,1.3)],
        #         [(cs,5, -1.3),(cs,10)]
            ],'CN')
        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
        survey = DC.SurveyDC(srcList)
        problem = DC.ProblemDC_CC(mesh)
        problem.pair(survey)
        mSynth = np.ones(mesh.nC)
        survey.makeSyntheticData(mSynth)
        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)
        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis
    def test_misfit(self):
        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
        self.assertTrue(passed)
    def test_adjoint(self):
        # Adjoint Test
        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
        v = np.random.rand(self.mesh.nC)
        w = np.random.rand(self.survey.dobs.shape[0])
        wtJv = w.dot(self.p.Jvec(self.m0, v))
        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
        passed = np.abs(wtJv - vtJtw) < 1e-10
        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
        self.assertTrue(passed)
    def test_dataObj(self):
        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
        self.assertTrue(passed)
    def test_massMatrices(self):
        Gu = np.random.rand(self.mesh.nF)
        def derChk(m):
            self.p.curModel = m
            return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
        self.assertTrue(passed)
 if __name__ == '__main__':
    unittest.main()
@@ -1,65 +0,0 @@
 import unittest
 import SimPEG.DCIP as DC
 from SimPEG import *
 class IPforwardTests(unittest.TestCase):
    def test_IPforward(self):
        cs = 12.5
        nc = 200/cs+1
        hx = [(cs,7, -1.3),(cs,nc),(cs,7, 1.3)]
        hy = [(cs,7, -1.3),(cs,int(nc/2+1)),(cs,7, 1.3)]
        hz = [(cs,7, -1.3),(cs,int(nc/2+1))]
        mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
        sighalf = 1e-2
        sigma = np.ones(mesh.nC)*sighalf
        p0 = np.r_[-50., 50., -50.]
        p1 = np.r_[ 50.,-50., -150.]
        blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
        sigma[blk_ind] = 1e-3
        eta = np.zeros_like(sigma)
        eta[blk_ind] = 0.1
        sigmaInf = sigma.copy()
        sigma0 = sigma*(1-eta)
        nElecs = 11
        x_temp = np.linspace(-100, 100, nElecs)
        aSpacing = x_temp[1]-x_temp[0]
        y_temp = 0.
        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
        survey = DC.SurveyDC(srcList)
        imap   = Maps.IdentityMap(mesh)
        problem = DC.ProblemDC_CC(mesh, mapping=imap)
        try:
            from pymatsolver import MumpsSolver
            solver = MumpsSolver
        except ImportError, e:
            solver = SolverLU
        problem.Solver = solver
        problem.pair(survey)
        phi0 = survey.dpred(sigma0)
        phiInf = survey.dpred(sigmaInf)
        phiIP_true = phi0-phiInf
        surveyIP = DC.SurveyIP(srcList)
        problemIP = DC.ProblemIP(mesh, sigma=sigma)
        problemIP.pair(surveyIP)
        problemIP.Solver = solver
        phiIP_approx = surveyIP.dpred(eta)
        err =  np.linalg.norm(phiIP_true-phiIP_approx) / np.linalg.norm(phiIP_true)
        self.assertTrue(err < 0.02)
 if __name__ == '__main__':
    unittest.main()
@@ -1,82 +0,0 @@
 import unittest
 from SimPEG import *
 import SimPEG.DCIP as DC
 class IPProblemTests(unittest.TestCase):
    def setUp(self):
        cs = 12.5
        nc = 500/cs+1
        hx = [(cs,0, -1.3),(cs,nc),(cs,0, 1.3)]
        hy = [(cs,0, -1.3),(cs,int(nc/2+1)),(cs,0, 1.3)]
        hz = [(cs,0, -1.3),(cs,int(nc/2+1))]
        mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
        sighalf = 1e-2
        sigma = np.ones(mesh.nC)*sighalf
        p0 = np.r_[-50., 50., -50.]
        p1 = np.r_[ 50.,-50., -150.]
        blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
        sigma[blk_ind] = 1e-3
        eta = np.zeros_like(sigma)
        eta[blk_ind] = 0.1
        nElecs = 5
        x_temp = np.linspace(-250, 250, nElecs)
        aSpacing = x_temp[1]-x_temp[0]
        y_temp = 0.
        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
        survey = DC.SurveyIP(srcList)
        imap   = Maps.IdentityMap(mesh)
        problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap)
        problem.pair(survey)
        try:
            from pymatsolver import MumpsSolver
            problem.Solver = MumpsSolver
        except ImportError, e:
            problem.Solver = SolverLU
        mSynth = eta
        survey.makeSyntheticData(mSynth)
        # Now set up the problem to do some minimization
        dmis = DataMisfit.l2_DataMisfit(survey)
        reg = Regularization.Tikhonov(mesh)
        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
        inv = Inversion.BaseInversion(invProb)
        self.inv = inv
        self.reg = reg
        self.p = problem
        self.mesh = mesh
        self.m0 = mSynth
        self.survey = survey
        self.dmis = dmis
    def test_misfit(self):
        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
        passed = Tests.checkDerivative(derChk, self.m0*0, plotIt=False)
        self.assertTrue(passed)
    def test_adjoint(self):
        # Adjoint Test
        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
        v = np.random.rand(self.mesh.nC)
        w = np.random.rand(self.survey.dobs.shape[0])
        wtJv = w.dot(self.p.Jvec(self.m0, v))
        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
        passed = np.abs(wtJv - vtJtw) < 1e-10
        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
        self.assertTrue(passed)
    def test_dataObj(self):
        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
        self.assertTrue(passed)
 if __name__ == '__main__':
    unittest.main()
@@ -14,9 +14,9 @@ class Doc_Test(unittest.TestCase):
        html_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['html'])
        check = subprocess.call(["sphinx-build", "-nW", "-b", "html", "-d",
-            "{0!s}".format((doctrees_path)) ,
+            "%s"%(doctrees_path) ,
-            "{0!s}".format((self.path_to_docs)),
+            "%s"%(self.path_to_docs),
-            "{0!s}".format((html_path))])
+            "%s"%(html_path)])
        assert check == 0
    # def test_latex(self):
@@ -29,15 +29,15 @@ class Doc_Test(unittest.TestCase):
    #         "%s"%(latex_path)])
    #     assert check == 0
-    def test_linkcheck(self):
+    # def test_linkcheck(self):
-        doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees'])
+    #     doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees'])
-        link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build'])
+    #     link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build'])
-        check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d",
+    #     check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d",
-            "%s"%(doctrees_path),
+    #         "%s"%(doctrees_path),
-            "%s"%(self.path_to_docs),
+    #         "%s"%(self.path_to_docs),
-            "%s"%(link_path)])
+    #         "%s"%(link_path)])
-        assert check == 0
+    #     assert check == 0
 if __name__ == '__main__':
    unittest.main()
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
 def adjointTest(fdemType, comp):
    prb = getFDEMProblem(fdemType, comp, SrcList, freq)
-    print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
+    print 'Adjoint %s formulation - %s' % (fdemType, comp)
    m  = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
    mu = np.ones(prb.mesh.nC)*MU
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
 def adjointTest(fdemType, comp):
    prb = getFDEMProblem(fdemType, comp, SrcList, freq)
-    print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
+    print 'Adjoint %s formulation - %s' % (fdemType, comp)
    m  = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
    mu = np.ones(prb.mesh.nC)*MU
@@ -26,7 +26,7 @@ SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
 def derivTest(fdemType, comp):
    prb = getFDEMProblem(fdemType, comp, SrcType, freq)
-    print '{0!s} formulation - {1!s}'.format(fdemType, comp)
+    print '%s formulation - %s' % (fdemType, comp)
    x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
    mu = np.log(np.ones(prb.mesh.nC)*MU)
@@ -55,7 +55,7 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
    if showIt == True:
        plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
        plt.loglog(rx.times, abs(bz_ana), 'b*')
-        plt.title('sig_half = {0:e}'.format(sig_half))
+        plt.title('sig_half = %e'%sig_half)
        plt.show()
    return log10diff
@@ -25,7 +25,7 @@ class compareInitFiles(unittest.TestCase):
 def get(test):
    def test_func(self):
-        print '\nTesting {0!s}.run(plotIt=False)\n'.format(test)
+        print '\nTesting %s.run(plotIt=False)\n'%test
        getattr(Examples, test).run(plotIt=False)
        self.assertTrue(True)
    return test_func
@@ -121,7 +121,7 @@ class RichardsTests1D(unittest.TestCase):
        tol = TOL*(10**int(np.log10(np.abs(zJv))))
        passed = np.abs(vJz - zJv) < tol
        print 'Richards Adjoint Test - PressureHead'
-        print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
+        print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
        self.assertTrue(passed,True)
    def test_Sensitivity(self):
@@ -193,7 +193,7 @@ class RichardsTests2D(unittest.TestCase):
        tol = TOL*(10**int(np.log10(np.abs(zJv))))
        passed = np.abs(vJz - zJv) < tol
        print '2D: Richards Adjoint Test - PressureHead'
-        print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
+        print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
        self.assertTrue(passed,True)
    def test_Sensitivity(self):
@@ -265,7 +265,7 @@ class RichardsTests3D(unittest.TestCase):
        tol = TOL*(10**int(np.log10(np.abs(zJv))))
        passed = np.abs(vJz - zJv) < tol
        print '3D: Richards Adjoint Test - PressureHead'
-        print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
+        print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
        self.assertTrue(passed,True)
    def test_Sensitivity(self):
@@ -10,9 +10,7 @@ class BasicCurvTests(unittest.TestCase):
        a = np.array([1, 1, 1])
        b = np.array([1, 2])
        c = np.array([1, 4])
-
+        gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
        def gridIt(h): return [np.cumsum(np.r_[0, x]) for x in h]
        X, Y = ndgrid(gridIt([a, b]), vector=False)
        self.TM2 = TensorMesh([a, b])
        self.Curv2 = CurvilinearMesh([X, Y])
@@ -21,10 +19,7 @@ class BasicCurvTests(unittest.TestCase):
        self.Curv3 = CurvilinearMesh([X, Y, Z])
    def test_area_3D(self):
-        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8,
+        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
                              1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4,
                              4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1,
                              1, 2, 2, 2])
        self.assertTrue(np.all(self.Curv3.area == test_area))
    def test_vol_3D(self):
@@ -38,85 +33,54 @@ class BasicCurvTests(unittest.TestCase):
        self.assertTrue(t1)
    def test_edge_3D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
                              1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
                              2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2,
                              2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4,
                              4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
        t1 = np.all(self.Curv3.edge == test_edge)
        self.assertTrue(t1)
    def test_edge_2D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
                              2])
        t1 = np.all(self.Curv2.edge == test_edge)
        self.assertTrue(t1)
    def test_tangents(self):
        T = self.Curv2.tangents
-        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx)))
-                        np.ones(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx)))
-        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy)))
-                        np.zeros(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy)))
        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] ==
                        np.zeros(self.Curv2.nEy)))
        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] ==
                        np.ones(self.Curv2.nEy)))
        T = self.Curv3.tangents
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx)))
-                        np.ones(self.Curv3.nEx)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx)))
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx)))
                        np.zeros(self.Curv3.nEx)))
        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] ==
                        np.zeros(self.Curv3.nEx)))
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy)))
-                        np.zeros(self.Curv3.nEy)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy)))
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy)))
                        np.ones(self.Curv3.nEy)))
        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] ==
                        np.zeros(self.Curv3.nEy)))
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz)))
-                        np.zeros(self.Curv3.nEz)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz)))
-        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz)))
                        np.zeros(self.Curv3.nEz)))
        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] ==
                        np.ones(self.Curv3.nEz)))
    def test_normals(self):
        N = self.Curv2.normals
-        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx)))
-                        np.ones(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx)))
-        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy)))
-                        np.zeros(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy)))
        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] ==
                        np.zeros(self.Curv2.nFy)))
        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] ==
                        np.ones(self.Curv2.nFy)))
        N = self.Curv3.normals
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx)))
-                        np.ones(self.Curv3.nFx)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx)))
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx)))
                        np.zeros(self.Curv3.nFx)))
        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] ==
                        np.zeros(self.Curv3.nFx)))
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy)))
-                        np.zeros(self.Curv3.nFy)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy)))
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy)))
                        np.ones(self.Curv3.nFy)))
        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] ==
                        np.zeros(self.Curv3.nFy)))
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz)))
-                        np.zeros(self.Curv3.nFz)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz)))
-        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] ==
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz)))
                        np.zeros(self.Curv3.nFz)))
        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] ==
                        np.ones(self.Curv3.nFz)))
    def test_grid(self):
        self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
@@ -2,12 +2,12 @@ import numpy as np
 import scipy.sparse as sp
 import unittest
 import matplotlib.pyplot as plt
-from SimPEG import Mesh, Tests, Utils, Solver
+from SimPEG import *
 MESHTYPES = ['uniformTensorMesh']
 def getxBCyBC_CC(mesh, alpha, beta, gamma):
 # def getxBCyBC(mesh, alpha, beta, gamma):
    """
    This is a subfunction generating mixed-boundary condition:
@@ -17,19 +17,17 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        \rho \vec{j} = -\nabla \phi \phi
-        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r
+        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
        = \partial \Omega
        xBC = f_1(\alpha, \beta, \gamma)
        yBC = f(\alpha, \beta, \gamma)
    Computes xBC and yBC for cell-centered discretizations
    """
-
+    if mesh.dim == 1: #1D
    if mesh.dim == 1:  # 1D
        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
            raise Exception("Lenght of list, alpha should be 2")
-        fCCxm, fCCxp = mesh.cellBoundaryInd
+        fCCxm,fCCxp = mesh.cellBoundaryInd
        nBC = fCCxm.sum()+fCCxp.sum()
        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
@@ -52,11 +50,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        xBC = np.r_[xBC_xm, xBC_xp]
        yBC = np.r_[yBC_xm, yBC_xp]
-    elif mesh.dim == 2:  # 2D
+    elif mesh.dim == 2: #2D
        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
            raise Exception("Lenght of list, alpha should be 4")
-        fxm, fxp, fym, fyp = mesh.faceBoundaryInd
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
@@ -67,10 +65,8 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
-        h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
        h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
        h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
@@ -91,10 +87,8 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        yBC_ym = 0.5*(1.-b_ym)
        yBC_yp = 0.5*(1.-1./b_yp)
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-                                np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
                                np.arange(mesh.nFy)[fyp]])
        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
@@ -104,11 +98,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        xBC = np.r_[xBC_x, xBC_y]
        yBC = np.r_[yBC_x, yBC_y]
-    elif mesh.dim == 3:  # 3D
+    elif mesh.dim == 3: #3D
        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
            raise Exception("Lenght of list, alpha should be 6")
        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
-        fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
@@ -122,12 +116,9 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
-        h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
-        h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
+        h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
        h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
        h_zm = mesh.hz[0]*np.ones_like(alpha_zm)
        h_zp = mesh.hz[-1]*np.ones_like(alpha_zp)
        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
@@ -157,12 +148,9 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
        yBC_zm = 0.5*(1.-b_zm)
        yBC_zp = 0.5*(1.-1./b_zp)
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-                                np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
-        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
                                np.arange(mesh.nFy)[fyp]])
        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm],
                                np.arange(mesh.nFz)[fzp]])
        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
@@ -177,7 +165,6 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
    return xBC, yBC
 class Test1D_InhomogeneousMixed(Tests.OrderTest):
    name = "1D - Mixed"
    meshTypes = MESHTYPES
@@ -186,14 +173,11 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
    meshSizes = [4, 8, 16, 32]
    def getError(self):
-        # Test function
+        #Test function
-        def phi_fun(x): return np.cos(np.pi*x)
+        phi_fun = lambda x: np.cos(np.pi*x)
-
+        j_fun = lambda x: np.pi*np.sin(np.pi*x)
-        def j_fun(x): return np.pi*np.sin(np.pi*x)
+        phi_deriv = lambda x: -j_fun(x)
-
+        q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
        def phi_deriv(x): return -j_fun(x)
        def q_fun(x): return (np.pi**2)*np.cos(np.pi*x)
        xc_ana = phi_fun(self.M.gridCC)
        q_ana = q_fun(self.M.gridCC)
@@ -205,16 +189,17 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
        # Setup Mixed B.C (alpha, beta, gamma)
        alpha_xm, alpha_xp = 1., 1.
-        beta_xm, beta_xp = 1., 1.
+        beta_xm, beta_xp =  1., 1.
        alpha = np.r_[alpha_xm, alpha_xp]
        beta = np.r_[beta_xm, beta_xp]
        vecN = self.M.vectorNx
        vecC = self.M.vectorCCx
-        phi_bc = phi_fun(vecN[[0, -1]])
+        phi_bc = phi_fun(vecN[[0,-1]])
-        phi_deriv_bc = phi_deriv(vecN[[0, -1]])
+        phi_deriv_bc = phi_deriv(vecN[[0,-1]])
        gamma = alpha*phi_bc + beta*phi_deriv_bc
        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
        sigma = np.ones(self.M.nC)
        Mfrho = self.M.getFaceInnerProduct(1./sigma)
        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -229,7 +214,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
        A = Div*MfrhoI*G
        if self.myTest == 'xc':
-            # TODO: fix the null space
+            #TODO: fix the null space
            Ainv = Solver(A)
            xc = Ainv*rhs
            err = np.linalg.norm((xc-xc_ana), np.inf)
@@ -237,13 +222,13 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
            NotImplementedError
        return err
    def test_order(self):
        print "==== Testing Mixed boudary conduction for CC-problem ===="
        self.name = "1D"
        self.myTest = 'xc'
        self.orderTest()
 class Test2D_InhomogeneousMixed(Tests.OrderTest):
    name = "2D - Mixed"
    meshTypes = MESHTYPES
@@ -252,59 +237,40 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
    meshSizes = [4, 8, 16, 32]
    def getError(self):
-        # Test function
+        #Test function
-        def phi_fun(x):
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
-            return np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
-
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
-        def j_funX(x):
+        phideriv_funX = lambda x: -j_funX(x)
-            return +np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1])
+        phideriv_funY = lambda x: -j_funY(x)
-
+        q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
        def j_funY(x):
            return +np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1])
        def phideriv_funX(x):
            return -j_funX(x)
        def phideriv_funY(x):
            return -j_funY(x)
        def q_fun(x):
            return +2*(np.pi**2)*phi_fun(x)
        xc_ana = phi_fun(self.M.gridCC)
        q_ana = q_fun(self.M.gridCC)
        jX_ana = j_funX(self.M.gridFx)
        jY_ana = j_funY(self.M.gridFy)
-        j_ana = np.r_[jX_ana, jY_ana]
+        j_ana = np.r_[jX_ana,jY_ana]
        # Get boundary locations
-        fxm, fxp, fym, fyp = self.M.faceBoundaryInd
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
-        gBFxm = self.M.gridFx[fxm, :]
+        gBFxm = self.M.gridFx[fxm,:]
-        gBFxp = self.M.gridFx[fxp, :]
+        gBFxp = self.M.gridFx[fxp,:]
-        gBFym = self.M.gridFy[fym, :]
+        gBFym = self.M.gridFy[fym,:]
-        gBFyp = self.M.gridFy[fyp, :]
+        gBFyp = self.M.gridFy[fyp,:]
        # Setup Mixed B.C (alpha, beta, gamma)
-        alpha_xm = np.ones_like(gBFxm[:, 0])
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        alpha_xp = np.ones_like(gBFxp[:, 0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        beta_xm = np.ones_like(gBFxm[:, 0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        beta_xp = np.ones_like(gBFxp[:, 0])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
        alpha_ym = np.ones_like(gBFym[:, 1])
        alpha_yp = np.ones_like(gBFyp[:, 1])
        beta_ym = np.ones_like(gBFym[:, 1])
        beta_yp = np.ones_like(gBFyp[:, 1])
        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
-        phiderivX_bc_xm = phideriv_funX(gBFxm)
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
-        phiderivX_bc_xp = phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
        phiderivY_bc_ym = phideriv_funY(gBFym)
        phiderivY_bc_yp = phideriv_funY(gBFyp)
        def gamma_fun(alpha, beta, phi, phi_deriv):
            return alpha*phi + beta*phi_deriv
        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
@@ -316,6 +282,7 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
        sigma = np.ones(self.M.nC)
        Mfrho = self.M.getFaceInnerProduct(1./sigma)
        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -336,13 +303,13 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
            NotImplementedError
        return err
    def test_order(self):
        print "==== Testing Mixed boudary conduction for CC-problem ===="
        self.name = "2D"
        self.myTest = 'xc'
        self.orderTest()
 class Test3D_InhomogeneousMixed(Tests.OrderTest):
    name = "3D - Mixed"
    meshTypes = MESHTYPES
@@ -351,74 +318,51 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
    meshSizes = [4, 8, 16]
    def getError(self):
-        # Test function
+        #Test function
-        def phi_fun(x):
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
-            return (np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
-                    np.cos(np.pi*x[:, 2]))
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
        j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
-        def j_funX(x):
+        phideriv_funX = lambda x: -j_funX(x)
-            return (np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
+        phideriv_funY = lambda x: -j_funY(x)
-                    np.cos(np.pi*x[:, 2]))
+        phideriv_funZ = lambda x: -j_funZ(x)
-        def j_funY(x):
+        q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
            return (np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1]) *
                    np.cos(np.pi*x[:, 2]))
        def j_funZ(x):
            return (np.pi*np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
                    np.sin(np.pi*x[:, 2]))
        def phideriv_funX(x): return -j_funX(x)
        def phideriv_funY(x): return -j_funY(x)
        def phideriv_funZ(x): return -j_funZ(x)
        def q_fun(x): return 3*(np.pi**2)*phi_fun(x)
        xc_ana = phi_fun(self.M.gridCC)
        q_ana = q_fun(self.M.gridCC)
        jX_ana = j_funX(self.M.gridFx)
        jY_ana = j_funY(self.M.gridFy)
-        j_ana = np.r_[jX_ana, jY_ana, jY_ana]
+        j_ana = np.r_[jX_ana,jY_ana,jY_ana]
        # Get boundary locations
-        fxm, fxp, fym, fyp, fzm, fzp = self.M.faceBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
-        gBFxm = self.M.gridFx[fxm, :]
+        gBFxm = self.M.gridFx[fxm,:]
-        gBFxp = self.M.gridFx[fxp, :]
+        gBFxp = self.M.gridFx[fxp,:]
-        gBFym = self.M.gridFy[fym, :]
+        gBFym = self.M.gridFy[fym,:]
-        gBFyp = self.M.gridFy[fyp, :]
+        gBFyp = self.M.gridFy[fyp,:]
-        gBFzm = self.M.gridFz[fzm, :]
+        gBFzm = self.M.gridFz[fzm,:]
-        gBFzp = self.M.gridFz[fzp, :]
+        gBFzp = self.M.gridFz[fzp,:]
        # Setup Mixed B.C (alpha, beta, gamma)
-        alpha_xm = np.ones_like(gBFxm[:, 0])
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        alpha_xp = np.ones_like(gBFxp[:, 0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        beta_xm = np.ones_like(gBFxm[:, 0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        beta_xp = np.ones_like(gBFxp[:, 0])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        alpha_ym = np.ones_like(gBFym[:, 1])
+        alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-        alpha_yp = np.ones_like(gBFyp[:, 1])
+        beta_zm, beta_zp =  np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-        beta_ym = np.ones_like(gBFym[:, 1])
+
        beta_yp = np.ones_like(gBFyp[:, 1])
        alpha_zm = np.ones_like(gBFzm[:, 2])
        alpha_zp = np.ones_like(gBFzp[:, 2])
        beta_zm = np.ones_like(gBFzm[:, 2])
        beta_zp = np.ones_like(gBFzp[:, 2])
        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
        phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
-        phiderivX_bc_xm = phideriv_funX(gBFxm)
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
-        phiderivX_bc_xp = phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
-        phiderivY_bc_ym = phideriv_funY(gBFym)
+        phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp)
        phiderivY_bc_yp = phideriv_funY(gBFyp)
        phiderivY_bc_zm = phideriv_funZ(gBFzm)
        phiderivY_bc_zp = phideriv_funZ(gBFzp)
        def gamma_fun(alpha, beta, phi, phi_deriv):
            return alpha*phi + beta*phi_deriv
        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
@@ -432,6 +376,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
        sigma = np.ones(self.M.nC)
        Mfrho = self.M.getFaceInnerProduct(1./sigma)
        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -445,7 +390,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
        A = Div*MfrhoI*G
        if self.myTest == 'xc':
-            # TODO: fix the null space
+            #TODO: fix the null space
            Ainv = Solver(A)
            xc = Ainv*rhs
            err = np.linalg.norm((xc-xc_ana), np.inf)
@@ -453,11 +398,14 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
            NotImplementedError
        return err
    def test_order(self):
        print "==== Testing Mixed boudary conduction for CC-problem ===="
        self.name = "3D"
        self.myTest = 'xc'
        self.orderTest()
 if __name__ == '__main__':
    unittest.main()
@@ -2,6 +2,8 @@ import numpy as np
 import unittest
 from SimPEG.Utils import mkvc
 from SimPEG import Mesh, Tests
 import unittest
 MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
 TOLERANCES = [0.9, 0.5, 0.5]
@@ -52,7 +52,7 @@ class Tests(unittest.TestCase):
        assert o <= -1
        assert not (o > -1)
        assert o >= -1
-        assert -1.*(-o)*o == -o
+        assert -(-o)*o == -o
        o = Identity()
        assert +o == o
        assert -o == -o
Author	SHA1	Message	Date
Rowan Cockett	51c112a529	Updates to documentation and DC utilities.	2016-06-26 16:07:35 -06:00
Rowan Cockett	57f8e35db7	Merge remote-tracking branch 'refs/remotes/origin/dev' Conflicts: SimPEG/DCIP/BaseDC.py SimPEG/DCIP/BaseIP.py	2016-06-26 14:08:28 -06:00
Rowan Cockett	16a016e9ac	switch coverage and delete old tests.	2016-05-31 14:09:57 -07:00
Rowan Cockett	c9017bd180	Delete DCIP folder.	2016-05-31 08:42:52 -07:00