Initial commit of DCIP.

.travis.yml

@@ -18,6 +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/flow
     - TEST_DIR=tests/examples
     - TEST_DIR=tests/em/fdem/inverse/adjoint

SimPEG/DCIP/BaseDC.py

@@ -0,0 +1,853 @@
+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.
+
+    """
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+        # self._rhsDict = {}
+        self._Ps = {}
+
+    def projectFields(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.array m: model
+            :rtype: scipy.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, u=None):
+        """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: 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 u is None:
+            # Run forward simulation if $u$ not provided
+            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
+        else:
+            u = u[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, u=None):
+
+        self.curModel = m
+        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
+        if u is None:
+            # Run forward simulation if $u$ not provided
+            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
+        else:
+            u = u[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
+
+def readUBC_DC2DModel(fileName):
+
+    from SimPEG import np, mkvc
+    """
+        Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
+
+        Input:
+        :param fileName, path to the UBC GIF 2D model file
+
+        Output:
+        :param SimPEG TensorMesh 2D object
+        :return
+        
+        Created on Thu Nov 12 13:14:10 2015
+
+        @author: dominiquef
+
+    """
+    
+    # Open fileand skip header... assume that we know the mesh already
+
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+    
+    dim = np.array(obsfile[0].split(),dtype=float)
+    
+    temp = np.array(obsfile[1].split(),dtype=float)
+    
+    if len(temp) > 1:
+        model = np.zeros(dim)
+        
+        for ii in range(len(obsfile)-1):
+            mm = np.array(obsfile[ii+1].split(),dtype=float)
+            model[:,ii] = mm
+            
+        model = model[:,::-1]
+        
+    else:
+        
+        if len(obsfile[1:])==1:
+            mm = np.array(obsfile[1:].split(),dtype=float)
+            
+        else:
+            mm = np.array(obsfile[1:],dtype=float)
+            
+        # Permute the second dimension to flip the order
+        model = mm.reshape(dim[1],dim[0])
+    
+        model = model[::-1,:]
+        model = np.transpose(model, (1, 0))
+        
+    model = mkvc(model)
+
+
+    return model
+
+def plot_pseudoSection(Tx,Rx,data,z0, stype):
+    
+    from SimPEG import np, mkvc
+    from scipy.interpolate import griddata
+    from matplotlib.colors import LogNorm
+    import pylab as plt
+    import re
+    """
+        Read list of 2D tx-rx location and plot a speudo-section of apparent
+        resistivity.
+        
+        Assumes flat topo for now...
+    
+        Input:
+        :param d2D, z0
+        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
+    
+        Output:
+        :figure scatter plot overlayed on image
+        
+        Created on Mon December 7th, 2015
+    
+        @author: dominiquef
+    
+    """
+    #d2D = np.asarray(d2D)
+    
+    midl = []
+    midz = []
+    rho = []
+    
+    for ii in range(len(Tx)):
+        # Get distances between each poles
+        rC1P1 = np.abs(Tx[ii][0] - Rx[ii][:,0]) 
+        rC2P1 = np.abs(Tx[ii][1] - Rx[ii][:,0])
+        rC1P2 = np.abs(Tx[ii][1] - Rx[ii][:,1])
+        rC2P2 = np.abs(Tx[ii][0] - Rx[ii][:,1])
+        rP1P2 = np.abs(Rx[ii][:,1] - Rx[ii][:,0])    
+    
+        # Compute apparent resistivity
+        if re.match(stype,'pdp'):
+            rho = np.hstack([rho, data[ii] * 2*np.pi  * rC1P1 * ( rC1P1 + rP1P2 ) / rP1P2] )
+            
+        elif re.match(stype,'dpdp'):
+            rho = np.hstack([rho, data[ii] * 2*np.pi / ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 ) ])
+    
+        Cmid = (Tx[ii][0] + Tx[ii][1])/2
+        Pmid = (Rx[ii][:,0] + Rx[ii][:,1])/2
+    
+        midl = np.hstack([midl, ( Cmid + Pmid )/2 ])
+        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
+    
+   
+    # Grid points
+    grid_x, grid_z = np.mgrid[np.min(midl):np.max(midl), np.min(midz):np.max(midz)]
+    grid_rho = griddata(np.c_[midl,midz], np.log10(abs(1/rho.T)), (grid_x, grid_z), method='linear')
+    
+    
+    #plt.subplot(2,1,2)
+    plt.imshow(grid_rho.T, extent = (np.min(midl),np.max(midl),np.min(midz),np.max(midz)), origin='lower', alpha=0.8)
+    cbar = plt.colorbar(format = '%.2f',fraction=0.02)
+    cmin,cmax = cbar.get_clim()
+    ticks = np.linspace(cmin,cmax,3)
+    cbar.set_ticks(ticks)
+    
+    # Plot apparent resistivity
+    plt.scatter(midl,midz,s=50,c=np.log10(abs(1/rho.T)))
+
+def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
+    
+    from SimPEG import np
+    import re
+    """
+        Load in endpoints and survey specifications to generate Tx, Rx location
+        stations.
+        
+        Assumes flat topo for now...
+    
+        Input:
+        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
+        :object mesh -> SimPEG mesh object
+        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
+        : param a, n -> pole seperation, number of rx dipoles per tx
+        
+        Output:
+        :param Tx, Rx -> List objects for each tx location
+            Lines: P1x, P1y, P1z, P2x, P2y, P2z
+        
+        Created on Wed December 9th, 2015
+    
+        @author: dominiquef
+    
+    """
+    def xy_2_r(x1,x2,y1,y2):
+        r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) )
+        return r 
+        
+    ## Evenly distribute electrodes and put on surface
+    # Mesure survey length and direction
+    dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
+    
+    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
+    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
+       
+    nstn = np.floor( dl_len / a )
+    
+    # Compute discrete pole location along line
+    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
+    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
+    
+    # Create line of P1 locations
+    M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+    
+    # Create line of P2 locations
+    N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+    
+    ## Build list of Tx-Rx locations depending on survey type
+    # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
+    # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B]
+    Tx = []
+    Rx = []
+    
+    if not re.match(stype,'gradient'):
+        
+        for ii in range(0, int(nstn)-1): 
+            
+            
+            if re.match(stype,'dpdp'):
+                tx = np.c_[M[ii,:],N[ii,:]]
+            elif re.match(stype,'pdp'):
+                tx = np.c_[M[ii,:],M[ii,:]]
+                
+            #Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
+            
+            # Current elctrode seperation
+            AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
+            
+            # Number of receivers to fit
+            nstn = np.min([np.floor( (AB - b) / a ) , n])
+            
+            # Check if there is enough space, else break the loop
+            if nstn <= 0:
+                continue
+            
+            # Compute discrete pole location along line
+            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
+            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
+            
+            # Create receiver poles
+            # Create line of P1 locations
+            P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            
+            # Create line of P2 locations
+            P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            
+            Rx.append(np.c_[P1,P2])
+            Tx.append(tx)            
+            
+#==============================================================================
+#     elif re.match(stype,'dpdp'):
+#         
+#         for ii in range(0, int(nstn)-2):  
+#             
+#             indx = np.min([ii+n+1,nstn])
+#             Tx.append(np.c_[M[ii,:],N[ii,:]])
+#             Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
+#==============================================================================
+            
+    elif re.match(stype,'gradient'):
+        
+        # Gradient survey only requires Tx at end of line and creates a square
+        # grid of receivers at in the middle at a pre-set minimum distance
+        Tx.append(np.c_[M[0,:],N[-1,:]])
+              
+        # Get the edge limit of survey area
+        min_x = endl[0,0] + dl_x * b
+        min_y = endl[0,1] + dl_y * b
+            
+        max_x = endl[1,0] - dl_x * b
+        max_y = endl[1,1] - dl_y * b
+        
+        box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
+        box_w = box_l/2.
+        
+        nstn = np.floor( box_l / a )
+        
+        # Compute discrete pole location along line
+        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
+        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
+        
+        # Define number of cross lines
+        nlin = int(np.floor( box_w / a ))
+        lind = range(-nlin,nlin+1) 
+        
+        ngrad = nstn * len(lind)
+        
+        rx = np.zeros([ngrad,6])
+        for ii in range( len(lind) ):
+            
+            # Move line in perpendicular direction by dipole spacing
+            lxx = stn_x - lind[ii]*a*dl_y
+            lyy = stn_y + lind[ii]*a*dl_x
+            
+            
+            M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
+            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            
+            rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
+            
+        Rx.append(rx)
+        
+    else:
+        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
+
+        
+   
+    return Tx, Rx  
+
+def writeUBC_DCobs(fileName,Tx,Rx,d,wd, dtype):
+    
+    from SimPEG import np, mkvc
+    import re
+    """
+        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
+    
+        Input:
+        :param fileName, path to the UBC GIF 3D obs file
+    
+        Output:
+        :param rx, tx, d, wd
+        :return
+        
+        Created on Mon December 7th, 2015
+    
+        @author: dominiquef
+    
+    """
+    fid = open(fileName,'w')
+    fid.write('! GENERAL FORMAT\n')  
+    
+    for ii in range(len(Tx)):
+        
+        tx = np.asarray(Tx[ii])
+        rx = np.asarray(Rx[ii])
+        nrx = rx.shape[0]
+        
+        fid.write('\n')
+        
+        if re.match(dtype,'2D'):
+            
+            for jj in range(nrx):
+                
+                fid.writelines("%e " % ii for ii in mkvc(tx))
+                fid.writelines("%e " % ii for ii in mkvc(rx[jj]))
+                fid.write('%e %e\n'% (d[ii][jj],wd[ii][jj]))
+                #np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') 
+        
+        elif re.match(dtype,'3D'):
+                                
+            fid.write('\n')
+            fid.writelines("%e " % ii for ii in mkvc(tx))
+            fid.write('%i\n'% nrx)
+            np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') 
+        
+        
+    fid.close()  
+
+def convertObs_DC3D_to_2D(Tx,Rx):
+    
+    from SimPEG import np
+    import numpy.matlib as npm
+    """
+        Read list of 3D Tx Rx location and change coordinate system to distance
+        along line assuming all data is acquired along line
+        First transmitter pole is assumed to be at the origin
+
+        Assumes flat topo for now...
+    
+        Input:
+        :param Tx, Rx
+        
+        Output:
+        :figure Tx2d, Rx2d
+        
+        Created on Mon December 7th, 2015
+    
+        @author: dominiquef
+    
+    """
+    
+                
+    Tx2d = []
+    Rx2d = []
+
+    for ii in range(len(Tx)):
+        
+        if ii == 0:
+            endp = Tx[0][0:2,0]
+        
+        nrx = Rx[ii].shape[0]
+                  
+        rP1 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,0] )**2 , axis=0))
+        rP2 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,1] )**2 , axis=0))
+        rC1 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,0:2] )**2 , axis=1))
+        rC2 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,3:5] )**2 , axis=1))
+        
+        Tx2d.append( np.r_[rP1, rP2] )
+        Rx2d.append( np.c_[rC1, rC2] )
+            #np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n')
+
+    return Tx2d, Rx2d
+
+def readUBC_DC3Dobs(fileName):
+    
+    from SimPEG import np
+    """
+        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
+    
+        Input:
+        :param fileName, path to the UBC GIF 3D obs file
+    
+        Output:
+        :param rx, tx, d, wd
+        :return
+        
+        Created on Mon December 7th, 2015
+    
+        @author: dominiquef
+    
+    """
+       
+    # Load file
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+    
+    # Pre-allocate
+    Tx = []
+    Rx = []
+    d = []
+    wd = []
+    
+    # Countdown for number of obs/tx
+    count = 0
+    for ii in range(obsfile.shape[0]):
+        
+        if not obsfile[ii]:
+            continue
+        
+        # First line is transmitter with number of receivers
+        if count==0:
+    
+            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+            count = int(temp[-1])
+            temp = np.reshape(temp[0:-1],[2,3]).T
+            
+            Tx.append(temp)
+            rx = []
+            continue
+        
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+        
+            
+        rx.append(temp)          
+        
+        count = count -1        
+        
+        # Reach the end of  
+        if count == 0:
+            temp = np.asarray(rx)
+            Rx.append(temp[:,0:6])
+            
+            # Check for data + uncertainties
+            if temp.shape[1]==8:
+                d.append(temp[:,6])
+                wd.append(temp[:,7])
+                
+            # Check for data only    
+            elif temp.shape[1]==7:
+                d.append(temp[:,6])
+            
+    return Tx, Rx, d, wd
+
+def readUBC_DC2DLoc(fileName):
+
+    from SimPEG import np
+    """
+        Read UBC GIF 2D observation file and generate arrays for tx-rx location
+
+        Input:
+        :param fileName, path to the UBC GIF 2D model file
+
+        Output:
+        :param rx, tx
+        :return
+        
+        Created on Thu Nov 12 13:14:10 2015
+
+        @author: dominiquef
+
+    """
+    
+    # Open fileand skip header... assume that we know the mesh already
+#==============================================================================
+#     fopen = open(fileName,'r')
+#     lines = fopen.readlines()
+#     fopen.close()
+#==============================================================================
+
+    # Load file
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+    
+    # Check first line and figure out if 2D or 3D file format
+    line = np.array(obsfile[0].split(),dtype=float)   
+    
+    tx_A  = []
+    tx_B  = []
+    rx_M  = []
+    rx_N  = []
+    d   = []
+    wd  = []
+    
+    for ii in range(obsfile.shape[0]):
+        
+        # If len==3, then simple format where tx-rx is listed on each line
+        if len(line) == 4:
+        
+            temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+            tx_A = np.hstack((tx_A,temp[0]))
+            tx_B = np.hstack((tx_B,temp[1]))
+            rx_M = np.hstack((rx_M,temp[2]))
+            rx_N = np.hstack((rx_N,temp[3]))
+            
+        
+    rx = np.transpose(np.array((rx_M,rx_N)))    
+    tx = np.transpose(np.array((tx_A,tx_B)))
+    
+    return tx, rx, d, wd
+
+def readUBC_DC2DMesh(fileName):
+
+    from SimPEG import np
+    """
+        Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
+
+        Input:
+        :param fileName, path to the UBC GIF mesh file
+
+        Output:
+        :param SimPEG TensorMesh 2D object
+        :return
+        
+        Created on Thu Nov 12 13:14:10 2015
+
+        @author: dominiquef
+
+    """
+    
+    # Open file 
+    fopen = open(fileName,'r')
+    
+    # Read down the file and unpack dx vector
+    def unpackdx(fid,nrows):
+        for ii in range(nrows):
+            
+            line = fid.readline()
+            var = np.array(line.split(),dtype=float)
+            
+            if ii==0:
+                x0= var[0]
+                xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2])
+                xend = var[1]            
+                
+            else:
+                xvec = np.hstack((xvec,np.ones(int(var[1])) * (var[0] - xend) / int(var[1])))
+                xend = var[0] 
+                
+        return x0, xvec
+        
+    #%% Start with dx block
+    # First line specifies the number of rows for x-cells
+    line = fopen.readline()
+    nl = np.array(line.split(),dtype=float)
+
+    [x0, dx] = unpackdx(fopen,nl) 
+    
+
+    #%% Move down the file until reaching the z-block
+    line = fopen.readline()
+    if not line:
+        line = fopen.readline()
+        
+    #%% End with dz block
+    # First line specifies the number of rows for z-cells
+    line = fopen.readline()
+    nl = np.array(line.split(),dtype=float)
+
+    [z0, dz] = unpackdx(fopen,nl) 
+    
+    # Flip z0 to be the bottom of the mesh for SimPEG
+    z0 = z0 - sum(dz)
+    dz = dz[::-1]
+    #%% Make the mesh using SimPEG
+    
+    from SimPEG import Mesh
+    tensMsh = Mesh.TensorMesh([dx,dz],(x0, z0))
+    return tensMsh
+

SimPEG/DCIP/BaseIP.py

@@ -0,0 +1,182 @@
+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, u=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pu(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.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, u=None):
+        return self.forward(v)
+
+    def Jtvec(self, m, v, u=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
+

SimPEG/DCIP/__init__.py

@@ -0,0 +1,2 @@
+from BaseDC import *
+from BaseIP import *

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -0,0 +1,69 @@
+from SimPEG import *
+import SimPEG.DCIP as DC
+import matplotlib.pyplot as plt
+
+
+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)]
+    hz = [(cs,7, -1.3),(cs,20)]
+    mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
+    sighalf = 1e-2
+    sigma = np.ones(mesh.nC)*sighalf
+    xtemp = np.linspace(-150, 150, 21)
+    ytemp = np.linspace(-150, 150, 21)
+    xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
+    xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
+    xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
+
+    # if plotIt:
+    #     fig, ax = plt.subplots(1,1, figsize = (5,5))
+    #     mesh.plotSlice(sigma, grid=True, ax = ax)
+    #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
+    #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
+
+    rx = DC.RxDipole(xyz_rxP, xyz_rxN)
+    src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
+    survey = DC.SurveyDC([src])
+    problem = DC.ProblemDC_CC(mesh)
+    problem.pair(survey)
+    try:
+        from pymatsolver import MumpsSolver
+        problem.Solver = MumpsSolver
+    except Exception, e:
+        pass
+    data = survey.dpred(sigma)
+
+    def DChalf(srclocP, srclocN, rxloc, sigma, I=1.):
+        rp = (srclocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
+        rn = (srclocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
+        rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
+        rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
+        return I/(sigma*2.*np.pi)*(1/rP-1/rN)
+
+    data_anaP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
+    data_anaN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
+    data_ana = data_anaP-data_anaN
+    Data_ana = data_ana.reshape((21, 21), order = 'F')
+    Data = data.reshape((21, 21), order = 'F')
+    X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
+    Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
+
+    if plotIt:
+        fig, ax = plt.subplots(1,2, figsize = (12, 5))
+        vmin = np.r_[data, data_ana].min()
+        vmax = np.r_[data, data_ana].max()
+        dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
+        dat0 = ax[0].contourf(X, Y, Data_ana, 60, vmin = vmin, vmax = vmax)
+        cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
+        cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
+        ax[1].set_title('Analytic')
+        ax[0].set_title('Computed')
+        plt.show()
+
+    return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
+
+
+if __name__ == '__main__':
+    print run(plotIt=True)

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -0,0 +1,179 @@
+from SimPEG import *
+import SimPEG.DCIP as DC
+import scipy.interpolate as interpolation
+import matplotlib.pyplot as plt
+import time
+import re
+
+def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp',  plotIt=True):
+    """
+    DC Forward Simulation
+
+    Forward model conductive spheres in a half-space and plot a pseudo-section
+
+    Created on Mon Feb 01 19:28:06 2016
+
+    @fourndo
+    """
+
+    # First we need to create a mesh and a model.
+
+    # This is our mesh
+    dx    = 5.
+
+    hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)]
+    hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)]
+    hzind = [(dx,15,-1.3),(dx, 15)]
+
+    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
+
+
+    # Set background conductivity
+    model = np.ones(mesh.nC) * sig[0]
+
+    # First anomaly
+    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
+    model[ind] = sig[1]
+
+    # Second anomaly
+    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC)
+    model[ind] = sig[2]
+
+    # Get index of the center
+    indy = int(mesh.nCy/2)
+
+
+    # Plot the model for reference
+    # Define core mesh extent
+    xlim = 200
+    zlim = 125
+
+    # Specify the survey type: "pdp" | "dpdp"
+
+
+    # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
+    ends = [(-175,0),(175,0)]
+    ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]]
+
+    # Snap the endpoints to the grid. Easier to create 2D section.
+    indx = Utils.closestPoints(mesh, ends )
+    locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
+
+    # 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, stype, param[0], param[1], param[2])
+
+    # Define some global geometry
+    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
+    dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
+    dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
+    azm =  np.arctan(dl_y/dl_x)
+
+    #Set boundary conditions
+    mesh.setCellGradBC('neumann')
+
+    # Define the differential operators needed for the DC problem
+    Div = mesh.faceDiv
+    Grad = mesh.cellGrad
+    Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
+
+    A = Div*Msig*Grad
+
+    # Change one corner to deal with nullspace
+    A[0,0] = 1
+    A = sp.csc_matrix(A)
+
+    # We will solve the system iteratively, so a pre-conditioner is helpful
+    # This is simply a Jacobi preconditioner (inverse of the main diagonal)
+    dA = A.diagonal()
+    P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
+
+    # Now we can solve the system for all the transmitters
+    # We want to store the data
+    data = []
+
+    # There is probably a more elegant way to do this, but we can just for-loop through the transmitters
+    for ii in range(len(Tx)):
+
+        start_time = time.time() # Let's time the calculations
+
+        #print("Transmitter %i / %i\r" % (ii+1,len(Tx)))
+
+        # Select dipole locations for receiver
+        rxloc_M = np.asarray(Rx[ii][:,0:3])
+        rxloc_N = np.asarray(Rx[ii][:,3:])
+
+
+        # For usual cases "dpdp" or "gradient"
+        if not re.match(stype,'pdp'):
+            inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
+            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
+
+        else:
+
+            # Create an "inifinity" pole
+            tx =  np.squeeze(Tx[ii][:,0:1])
+            tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
+            inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
+            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
+
+
+        # Iterative Solve
+        Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
+
+        # We now have the potential everywhere
+        phi = mkvc(Ainvb[0])
+
+        # Solve for phi on pole locations
+        P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
+        P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
+
+        # Compute the potential difference
+        dtemp = (P1*phi - P2*phi)*np.pi
+
+        data.append( dtemp )
+        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
+
+    print 'Transmitter {0} of {1}'.format(ii,len(Tx))
+    print 'Forward completed'
+
+
+    # Let's just convert the 3D format into 2D (distance along line) and plot
+    [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx)
+
+
+    # Here is an example for the first tx-rx array
+    if plotIt:
+        fig = plt.figure()
+        ax = plt.subplot(2,1,1, aspect='equal')
+        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
+        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
+        plt.gca().set_aspect('equal', adjustable='box')
+
+        plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
+        plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
+        plt.xlim([-xlim,xlim])
+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+
+
+        ax = plt.subplot(2,1,2, aspect='equal')
+
+        # Plot the location of the spheres for reference
+        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        ax.add_artist(circle1)
+        ax.add_artist(circle2)
+
+        # Add the speudo section
+        DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype)
+
+        plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
+        plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
+        plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+
+        plt.show()
+
+        return fig, ax
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/DC_Forward_WennerArray.py

@@ -0,0 +1,65 @@
+from SimPEG import *
+import SimPEG.DCIP as DC
+import matplotlib.pyplot as plt
+
+
+def getSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+
+    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
+
+
+
+def run(plotIt=False,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')
+    if plotIt:
+        mesh.plotGrid(showIt=True)
+
+    srcList = getSrcList(nElecs, aSpacing, in2D=True)
+    survey = DC.SurveyDC(srcList)
+    problem = DC.ProblemDC_CC(mesh)
+    problem.pair(survey)
+
+    return mesh, survey, problem
+
+
+if __name__ == '__main__':
+    run(plotIt=True)

docs/api_DC.rst

@@ -0,0 +1,150 @@
+.. _api_DC:
+
+.. math::
+
+    \renewcommand{\div}{\nabla\cdot\,}
+    \newcommand{\grad}{\vec \nabla}
+    \newcommand{\curl}{{\vec \nabla}\times\,}
+    \newcommand{\dcurl}{{\mathbf C}}
+    \newcommand{\dgrad}{{\mathbf G}}
+    \newcommand{\Acf}{{\mathbf A_c^f}}
+    \newcommand{\Ace}{{\mathbf A_c^e}}
+    \renewcommand{\S}{{\mathbf \Sigma}}
+    \renewcommand{\Div}{{\mathbf {Div}}}
+    \renewcommand{\Grad}{{\mathbf {Grad}}}
+    \newcommand{\St}{{\mathbf \Sigma_\tau}}
+    \newcommand{\diag}{\mathbf{diag}}
+    \newcommand{\M}{{\mathbf M}}
+    \newcommand{\Me}{{\M^e}}
+    \newcommand{\Mes}[1]{{\M^e_{#1}}}
+    \newcommand{\be}{\mathbf{e}}
+    \newcommand{\bj}{\mathbf{j}}
+    \newcommand{\bphi}{\mathbf{\phi}}
+    \newcommand{\bq}{\mathbf{q}}
+    \newcommand{\bJ}{\mathbf{J}}
+    \newcommand{\bG}{\mathbf{G}}
+    \newcommand{\bP}{\mathbf{P}}
+    \newcommand{\bA}{\mathbf{A}}
+    \newcommand{\bm}{\mathbf{m}}
+    \newcommand{\B}{\vec{B}}
+    \newcommand{\D}{\vec{D}}
+    \renewcommand{\H}{\vec{H}}
+    \renewcommand {\j}  { {\vec j} }
+    \newcommand {\h}  { {\vec h} }
+    \renewcommand {\b}  { {\vec b} }
+    \newcommand {\e}  { {\vec e} }
+    \newcommand {\c}  { {\vec c} }
+    \renewcommand {\d}  { {\vec d} }
+    \renewcommand {\u}  { {\vec u} }
+    \newcommand{\I}{\vec{I}}
+
+DC resistivity survey
+*********************
+
+Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
+
+
+Background
+==========
+
+As direct current (DC) implies, in DC resistivity survey, we assume steady-state. We consider Maxwell's equations in steady state as
+
+.. math::
+
+    \curl \frac{1}{\mu} \vec{b} - \j = \j_s \\
+
+    \curl \e = 0
+
+Then by taking \\(\\curl\\) for the first equation, we have
+
+.. math::
+
+     - \div\j = q \\
+
+
+where
+
+.. math::
+
+    \div \j_s = q = I(\delta(\vec{r}-\vec{r}_{s+})-\delta(\vec{r}-\vec{r}_{s-}))
+
+Since \\(\\curl \\e = 0\\), we have
+
+.. math::
+
+    \e = \grad \phi
+
+And by Ohm's law, we have
+
+.. math::
+
+    \j = \sigma \grad \phi
+
+Finally, we can compute the solution of the system:
+
+.. math::
+
+    - \div\j = q
+
+    \j = \sigma \grad \phi
+
+    \frac{\partial \phi}{\partial r}\Big|_{\partial \Omega_{BC}} = 0
+
+
+Discretization
+==============
+
+By using finite volume method (FVM), we discretize our system as
+
+.. math::
+
+    -\Div \bj = \bq
+
+    \diag(\Acf^{T}\sigma^{-1}) \bj = \Grad \bphi
+
+Here boundary condtions are embedded in the discrete differential operators. With some linear algebra we have
+
+.. math::
+
+    \bA\bphi  = -\bq
+
+where
+
+.. math::
+
+    \bA = \Div (\diag(\Acf^{T}\sigma^{-1}))^{-1} \Grad
+
+By solving this linear equation, we can compute the solution of \\(\\phi\\). Based on this discretization, we derive sensitivity in discretized space. Sensitivity matrix can be in general can be written as
+
+.. math ::
+
+    \bJ = -\bP\bA^{-1}\bG
+
+where
+
+.. math ::
+
+    \bP: \text{Projection}
+
+    \bJ = \bP\frac{\partial \phi}{\partial \bm}
+
+Here \\(\\bm\\) indicates model parameters in discretized space.
+
+Verification
+============
+
+Comparing to the analytic function:
+
+.. plot::
+
+    import simpegDC as DC
+    DC.Examples.Verification.run(plotIt=True)
+
+API
+===
+
+.. automodule:: simpegDC.BaseDC
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:

tests/dcip/__init__.py

@@ -0,0 +1,12 @@
+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)

tests/dcip/test_forward_DCproblem.py

@@ -0,0 +1,63 @@
+import unittest
+from SimPEG import *
+import simpegDCIP as DC
+
+
+class DCProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        mesh, survey, problem = DC.Examples.WennerArray.example()
+
+
+        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()

tests/dcip/test_forward_DCproblem_analytic.py

@@ -0,0 +1,12 @@
+import unittest
+import simpegDCIP as DC
+
+
+class DCAnalyticTests(unittest.TestCase):
+
+    def test_forwardAnalytic(self):
+        self.assertTrue(DC.Examples.Verification.run() < 0.1)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/dcip/test_forward_IPproblem.py

@@ -0,0 +1,57 @@
+import unittest
+import simpegDCIP as DC
+from SimPEG import *
+from pymatsolver import MumpsSolver
+
+class IPforwardTests(unittest.TestCase):
+
+    def test_IPforward(self):
+
+		cs = 12.5
+		nc = 500/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(-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.Examples.WennerArray.getSrcList(nElecs,aSpacing)
+		survey = DC.SurveyDC(srcList)
+
+		imap   = Maps.IdentityMap(mesh)
+		problem = DC.ProblemDC_CC(mesh, mapping= imap )
+		problem.Solver = MumpsSolver
+		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 = MumpsSolver
+		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()

tests/dcip/test_sens_IPproblem.py

@@ -0,0 +1,80 @@
+import unittest
+from SimPEG import *
+import simpegDCIP as DC
+from pymatsolver import MumpsSolver
+
+
+
+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.Examples.WennerArray.getSrcList(nElecs,aSpacing)
+        survey = DC.SurveyIP(srcList)
+        imap   = Maps.IdentityMap(mesh)
+        problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap)
+        problem.pair(survey)
+        problem.Solver = MumpsSolver
+
+        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()