don't test 3D mesh2mesh (too big at the moment)

# Conflicts:
#	SimPEG/EM/Utils/testingUtils.py

README.rst

@@ -25,10 +25,6 @@ SimPEG
     :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
-    
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
 
 The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

SimPEG/DCIP/DCIPUtils.py

@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):
 
     return model
 
-def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
+def plot_pseudoSection(DCsurvey, axs, stype):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -179,7 +179,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
         Input:
         :param d2D, z0
         :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
-        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
+
         Output:
         :figure scatter plot overlayed on image
 
@@ -192,6 +192,9 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
     from scipy.interpolate import griddata
     import pylab as plt
 
+    # Set depth to 0 for now
+    z0 = 0.
+
     # Pre-allocate
     midx = []
     midz = []
@@ -218,97 +221,47 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
         Cmid = (Tx[0][0] + Tx[1][0])/2
         Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
 
-        # Change output for dtype
-        if dtype == 'volt':
+        # Compute pant leg of apparent rho
+        if stype == 'pdp':
+            leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
 
-            rho = np.hstack([rho,data])
+            leg = np.log10(abs(1/leg))
 
-        else:
+        elif stype == 'dpdp':
+            leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
 
-            # Compute pant leg of apparent rho
-            if stype == 'pdp':
-
-                leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
-
-            elif stype == 'dpdp':
-
-                leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
-
-            else:
-                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
-                break
-
-
-            if dtype == 'appc':
-
-                leg = np.log10(abs(1./leg))
-                rho = np.hstack([rho,leg])
-
-            elif dtype == 'appr':
-
-                leg = np.log10(abs(leg))
-                rho = np.hstack([rho,leg])
-
-            else:
-                print """dtype must be 'appr' | 'appc' | 'volt' """
-                break
 
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
-        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
+        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
+        rho = np.hstack([rho,leg])
+
+
+    ax = axs
 
     # Grid points
     grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
     grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
-    
-    # Scale the color scheme
-    if clim == None:
-        vmin, vmax = rho.min(), rho.max()
-    else:
-        vmin, vmax = clim[0], clim[1]
-
-    # Plot data
-    grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
-    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax)
-    plt.gca().tick_params(axis='both', which='major', labelsize=8)
-    
-    if contour is not None:
-        plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax)
-        
-    # Add scatter points
-    axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax)
-    
-    if colorbar:
-        
-        if dtype == 'volt':
-            cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
-
-        else:               
-            cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
-        
-        cmin,cmax = cbar.get_clim()
-        ticks = np.linspace(cmin,cmax,3)
-        cbar.set_ticks(ticks)
-        cbar.ax.tick_params(labelsize=10)
-        
-        if cblabel:
-            if dtype == 'appc':
-                cbar.set_label("App.Cond",size=12)
-            elif dtype == 'appr':
-                cbar.set_label("App.Res.",size=12)
-            elif dtype == 'volt':
-                cbar.set_label("Potential (V)",size=12)
 
 
+    plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin = np.min(rho), vmax = np.max(rho))
+    cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal")
 
-    if not axlabel:
-        axs.set_xticklabels([])
-        axs.set_yticklabels([])
+    cmin,cmax = cbar.get_clim()
+    ticks = np.linspace(cmin,cmax,3)
+    cbar.set_ticks(ticks)
 
+    # Plot apparent resistivity
+    plt.scatter(midx,midz,s=50,c=rho.T)
+
+    ax.set_xticklabels([])
+
+    ax.set_ylabel('Z')
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position('right')
     plt.gca().set_aspect('equal', adjustable='box')
 
 
-
-    return ph
+    return ax
 
 def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     """
@@ -408,6 +361,16 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
                 srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
             SrcList.append(srcClass)
 
+#==============================================================================
+#     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 stype == 'gradient':
 
         # Gradient survey only requires Tx at end of line and creates a square
@@ -460,15 +423,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
+def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
     """
         Write UBC GIF DCIP 2D or 3D observation file
 
         Input:
-        :string     fileName -> including path where the file is written out
-        :DCsurvey   DC survey class object
-        :string     dtype ->  either '2D' | '3D'
-        :string     stype ->  either 'SURFACE' | 'GENERAL'
+        :string fileName -> including path where the file is written out
+        :DCsurvey -> DC survey class object
+        :string dtype ->  either '2D' | '3D'
+        :string  stype ->  either 'SURFACE' | 'GENERAL'
 
         Output:
         :param UBC2D-Data file
@@ -483,16 +446,10 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
     assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
     assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-    
+
     fid = open(fileName,'w')
-    
-    
-    if iptype!=0:
-        fid.write('IPTYPE=%i\n'%iptype)
-        
-    else:
-        fid.write('! ' + stype + ' FORMAT\n')
-        
+    fid.write('! ' + stype + ' FORMAT\n')
+
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -516,7 +473,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
                 B = np.repeat(tx[0,1],M.shape[0],axis=0)
                 M = M[:,0]
                 N = N[:,0]
-                
+
                 np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
 
@@ -524,25 +481,18 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
                 if stype == 'SURFACE':
 
-                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
                     M = M[:,0]
                     N = N[:,0]
 
                 if stype == 'GENERAL':
 
-                    # Flip sign for z-elevation to depth
-                    tx[2::2,:] = -tx[2::2,:]
-                    
                     fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
                     M = M[:,0::2]
                     N = N[:,0::2]
 
-                    # Flip sign for z-elevation to depth
-                    M[:,1::2] = -M[:,1::2]
-                    N[:,1::2] = -N[:,1::2]
-                    
                 fid.write('%i\n'% nD)
-                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
         if dtype=='3D':
 
@@ -554,32 +504,31 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
             if stype == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
+                fid.writelines("%e " % ii for ii in mkvc(tx))
 
             fid.write('%i\n'% nD)
             np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
-            fid.write('\n')
-            
+
         count += nD
 
     fid.close()
 
-def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
+def convertObs_DC3D_to_2D(DCsurvey,lineID):
     """
-        Read DC survey and projects the coordinate system
-        according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
-        In the 'local' system, station coordinates are referenced
-        to distance from the first srcLoc[0].loc[0]
+        Read DC survey and data 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
 
-        The Z value is preserved, but Y coordinates zeroed.
+        Assumes flat topo for now...
 
         Input:
-        :param survey3D
+        :param Tx, Rx
 
         Output:
-        :figure survey2D
+        :figure Tx2d, Rx2d
 
-        Edited April 6th, 2016
+        Edited Feb 17th, 2016
 
         @author: dominiquef
 
@@ -621,39 +570,25 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
             Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs
             nrx = Rx[0].shape[0]
 
-            if flag == 'local':
-                # Find A electrode along line
-                vec, r = r_unit(x0,Tx[ii][0,0:2])
-                A = stn_id(vecTx,vec,r)
+            # Find A electrode along line
+            vec, r = r_unit(x0,Tx[ii][0,0:2])
+            A = stn_id(vecTx,vec,r)
 
-                # Find B electrode along line
-                vec, r = r_unit(x0,Tx[ii][1,0:2])
-                B = stn_id(vecTx,vec,r)
+            # Find B electrode along line
+            vec, r = r_unit(x0,Tx[ii][1,0:2])
+            B = stn_id(vecTx,vec,r)
 
-                M = np.zeros(nrx)
-                N = np.zeros(nrx)
-                for kk in range(nrx):
+            M = np.zeros(nrx)
+            N = np.zeros(nrx)
+            for kk in range(nrx):
 
-                    # Find all M electrodes along line
-                    vec, r = r_unit(x0,Rx[0][kk,0:2])
-                    M[kk] = stn_id(vecTx,vec,r)
+                # Find all M electrodes along line
+                vec, r = r_unit(x0,Rx[0][kk,0:2])
+                M[kk] = stn_id(vecTx,vec,r)
 
-                    # Find all N electrodes along line
-                    vec, r = r_unit(x0,Rx[1][kk,0:2])
-                    N[kk] = stn_id(vecTx,vec,r)
-            elif flag == 'Yloc':
-                """ Flip the XY axis locs"""
-                A = Tx[ii][0,1]
-                B = Tx[ii][1,1]
-                M = Rx[0][:,1]
-                N = Rx[1][:,1]
-
-            elif flag == 'Xloc':
-                """ Copy the rx-tx locs"""
-                A = Tx[ii][0,0]
-                B = Tx[ii][1,0]
-                M = Rx[0][:,0]
-                N = Rx[1][:,0]
+                # Find all N electrodes along line
+                vec, r = r_unit(x0,Rx[1][kk,0:2])
+                N[kk] = stn_id(vecTx,vec,r)
 
             Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]])
 
@@ -666,59 +601,51 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     DCsurvey2D.std = np.asarray(DCsurvey.std)
 
     return DCsurvey2D
-    
-def readUBC_DC3Dobs(fileName, dtype = 'DC'):
+
+def readUBC_DC3Dobs(fileName):
     """
-        Read UBC GIF IP 3D observation file and generate survey
+        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 IPsurvey
+        :param rx, tx, d, wd
         :return
 
+        Created on Mon December 7th, 2015
+
         @author: dominiquef
 
     """
-    zflag = True # Flag for z value provided
-    
+
     # Load file
-    if dtype == 'IP':
-        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
-        
-    elif dtype == 'DC':
-        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-        
-    else:
-        print "dtype must be 'DC'(default) | 'IP'"
-        
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+
     # Pre-allocate
     srcLists = []
     Rx = []
     d = []
     wd = []
-    
+    zflag = True # Flag for z value provided
 
     # Countdown for number of obs/tx
     count = 0
     for ii in range(obsfile.shape[0]):
 
-        # Skip if blank line
         if not obsfile[ii]:
             continue
 
-        # First line or end of a transmitter block, read transmitter info
+        # First line is transmitter with number of receivers
         if count==0:
-            # Read the line
-            temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
+
+            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
             count = int(temp[-1])
 
             # Check if z value is provided, if False -> nan
             if len(temp)==5:
-                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
-                
-                zflag = False # Pass on the flag to the receiver loc
+                tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
+                zflag = False
 
             else:
                 tx = temp[:-1]
@@ -726,16 +653,8 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
             rx = []
             continue
 
-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
 
-        # Filter out negative IP
-#        if temp[-2] < 0:       
-#            count = count -1
-#            print "Negative!"
-#            
-#        else:
-            
-        # If the Z-location is provided, otherwise put nan
         if zflag:
 
             rx.append(temp[:-2])
@@ -745,7 +664,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
                 wd.append(temp[-1])
 
         else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
             # Check if there is data with the location
             if len(temp)==6:
                 d.append(temp[-2])
@@ -753,7 +672,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
 
         count = count -1
 
-        # Reach the end of transmitter block, append the src, rx and continue
+        # Reach the end of transmitter block
         if count == 0:
             rx = np.asarray(rx)
             Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
@@ -769,7 +688,6 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
 
 def readUBC_DC2Dobs(fileName):
     """
-        ------- NEEDS TO BE UPDATED ------
         Read UBC GIF 2D observation file and generate arrays for tx-rx location
 
         Input:
@@ -817,73 +735,6 @@ def readUBC_DC2Dobs(fileName):
 
     return tx, rx, d, wd
 
-def readUBC_DC2Dpre(fileName):
-    """
-        Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location
-
-        Input:
-        :param fileName, path to the UBC GIF 3D obs file
-
-        Output:
-        DCsurvey
-        :return
-
-        Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>
-
-        @author: dominiquef
-
-    """
-
-    # Load file
-    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-
-    # Pre-allocate
-    srcLists = []
-    Rx = []
-    d = []
-    zflag = True # Flag for z value provided
-
-    for ii in range(obsfile.shape[0]):
-
-        if not obsfile[ii]:
-            continue
-
-        # First line is transmitter with number of receivers
-
-
-        temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
-
-
-        # Check if z value is provided, if False -> nan
-        if len(temp)==5:
-            tx = np.r_[temp[0],np.nan,np.nan,temp[1],np.nan,np.nan]
-            zflag = False
-
-        else:
-            tx = np.r_[temp[0],np.nan,temp[1],temp[2],np.nan,temp[3]]
-
-
-        if zflag:
-            rx = np.c_[temp[4],np.nan,temp[5],temp[6],np.nan,temp[7]]
-
-
-        else:
-            rx = np.c_[temp[2],np.nan,np.nan,temp[3],np.nan,np.nan]
-            # Check if there is data with the location
-
-        d.append(temp[-1])
-
-
-        Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
-        srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) )
-
-    # Create survey class
-    survey = DC.SurveyDC(srcLists)
-
-    survey.dobs = np.asarray(d)
-
-    return {'DCsurvey':survey}
-
 def readUBC_DC2DMesh(fileName):
     """
         Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
@@ -1077,6 +928,7 @@ def getSrc_locs(DCsurvey):
 
     srcMat = np.zeros((DCsurvey.nSrc,2,3))
     for ii in range(DCsurvey.nSrc):
+        print np.asarray(DCsurvey.srcList[ii].loc).shape
         srcMat[ii,:,:] =  np.asarray(DCsurvey.srcList[ii].loc)
 
     return srcMat

SimPEG/Directives.py

@@ -216,13 +216,13 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the data.
         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
         phi_ms = 0.5*ms.dot(ms)
-        if self.reg.mrefInSmooth == True:
+        if self.reg.smoothModel == True:
             mref = self.reg.mref
         else:
             mref = 0
         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
         phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim >= 2:
+        if self.prob.mesh.dim==2:
             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
             phi_my = 0.5 * my.dot(my)
         else:
@@ -237,6 +237,40 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration
+    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
+    """
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+
 
 # class UpdateReferenceModel(Parameter):
 
@@ -249,112 +283,3 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 #             mref = self.mref0
 #         self.m_prev = self.invProb.m_current
 #         return mref
-
-class Update_IRLS(InversionDirective):
-
-    eps_min = None
-    factor = None
-    gamma = None
-    phi_m_last = None
-    phi_d_last = None
-
-
-    def initialize(self):
-
-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
-
-    def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor
-
-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
-            else:
-                self.reg.eps = eps
-
-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
-
-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
-
-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
-
-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
-
-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
-
-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
-
-class Update_lin_PreCond(InversionDirective):
-
-
-    def endIter(self):
-        # Cool the threshold parameter
-
-        if getattr(self.opt, 'approxHinv', None) is not None:
-            # Update the pre-conditioner
-            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
-            self.opt.approxHinv = PC
-
-
-class Update_Wj(InversionDirective):
-    """
-        Create approx-sensitivity base weighting using the probing method
-    """
-    k = None # Number of probing cycles
-    itr = None # Iteration number to update Wj, or always update if None
-
-    def endIter(self):
-
-        if self.itr is None or self.itr == self.opt.iter:
-
-            m = self.invProb.curModel
-            if self.k is None:
-                self.k = int(self.survey.nD/10)
-
-            def JtJv(v):
-
-                Jv = self.prob.Jvec(m, v)
-
-                return self.prob.Jtvec(m,Jv)
-
-            JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k)
-            JtJdiag = JtJdiag / max(JtJdiag)
-
-            self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

SimPEG/EM/Base.py

@@ -7,7 +7,7 @@ class EMPropMap(Maps.PropMap):
     """
 
     sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
-    mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
+    mu = Maps.Property("Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
 
     rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
     mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))

SimPEG/EM/FDEM/FDEM.py

@@ -8,29 +8,29 @@ from SimPEG.EM.Utils import omega
 
 class BaseFDEMProblem(BaseEMProblem):
     """
-        We start by looking at Maxwell's equations in the electric
-        field \\\(\\\mathbf{e}\\\) and the magnetic flux
-        density \\\(\\\mathbf{b}\\\)
+    We start by looking at Maxwell's equations in the electric
+    field \\\(\\\mathbf{e}\\\) and the magnetic flux
+    density \\\(\\\mathbf{b}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
-            {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
+        \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
+        {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
+    if using the E-B formulation (:code:`Problem_e`
+    or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
-        If we write Maxwell's equations in terms of
-        \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
+    If we write Maxwell's equations in terms of
+    \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
-            \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
+        \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
+        \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
+    if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
-        The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+    The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
 
     surveyPair = SurveyFDEM
@@ -204,6 +204,17 @@ class Problem_e(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
+
     def getA(self, freq):
         """
         System matrix
@@ -283,7 +294,6 @@ class Problem_e(BaseFDEMProblem):
         if adjoint:
             dRHS = MfMui * (C * v)
             return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v)
-
         else:
             return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
 
@@ -315,6 +325,16 @@ class Problem_b(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
     def getA(self, freq):
         """
         System matrix
@@ -463,6 +483,16 @@ class Problem_j(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
     def getA(self, freq):
         """
         System matrix
@@ -601,6 +631,17 @@ class Problem_h(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
+
     def getA(self, freq):
         """
         System matrix

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -193,16 +193,6 @@ class Fields_e(Fields):
         self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
         self._MfMui = self.survey.prob.MfMui
 
-    def _GLoc(self, fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
 
     def _ePrimary(self, eSolution, srcList):
         """
@@ -465,17 +455,6 @@ class Fields_b(Fields):
         self._nC = self.survey.prob.mesh.nC
 
 
-
-    def _GLoc(self,fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _bPrimary(self, bSolution, srcList):
         """
         Primary magnetic flux density from source
@@ -729,16 +708,6 @@ class Fields_j(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _jPrimary(self, jSolution, srcList):
         """
         Primary current density from source
@@ -1024,16 +993,6 @@ class Fields_h(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _hPrimary(self, hSolution, srcList):
         """
         Primary magnetic field from source

SimPEG/EM/FDEM/SrcFDEM.py

@@ -9,8 +9,10 @@ class BaseSrc(Survey.BaseSrc):
     """
 
     freq = None
-    # rxPair = RxFDEM
-    integrate = True
+    integrate = False
+
+    def __init__(self, rxList, **kwargs):
+        Survey.BaseSrc.__init__(self, rxList, **kwargs)
 
     def eval(self, prob):
         """
@@ -103,7 +105,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def s_mDeriv(self, prob, v, adjoint = False):
+    def s_mDeriv(self, prob, v, adjoint=False):
         """
         Derivative of magnetic source term with respect to the inversion model
 
@@ -116,7 +118,7 @@ class BaseSrc(Survey.BaseSrc):
 
         return Zero()
 
-    def s_eDeriv(self, prob, v, adjoint = False):
+    def s_eDeriv(self, prob, v, adjoint=False):
         """
         Derivative of electric source term with respect to the inversion model
 
@@ -136,13 +138,12 @@ class RawVec_e(BaseSrc):
     :param list rxList: receiver list
     :param float freq: frequency
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+    def __init__(self, rxList, freq, s_e):
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
@@ -166,13 +167,12 @@ class RawVec_m(BaseSrc):
     :param float freq: frequency
     :param rxList: receiver list
     :param numpy.array s_m: magnetic source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
     def __init__(self, rxList, freq, s_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
         self._s_m = np.array(s_m, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
@@ -197,14 +197,13 @@ class RawVec(BaseSrc):
     :param float freq: frequency
     :param numpy.array s_m: magnetic source term
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
-    def __init__(self, rxList, freq, s_m, s_e, integrate=True):
+    def __init__(self, rxList, freq, s_m, s_e, **kwargs):
         self._s_m = np.array(s_m, dtype=complex)
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_m(self, prob):
         """
@@ -278,14 +277,13 @@ class MagDipole(BaseSrc):
     :param float mu: background magnetic permeability
     """
 
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
         self.moment = moment
         self.mu = mu
-        self.integrate = False
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
@@ -543,10 +541,10 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -605,5 +603,231 @@ class CircularLoop(BaseSrc):
             return -C.T * (MMui_s * self.bPrimary(prob))
 
 
+class PrimSec(BaseSrc):
+    """
+    Primary-Secondary source in the physical properties. A primary problem is
+    first solved, and the fields from this problem are used to construct a
+    source term for the secondary problem. Either a mesh and
+    fields need to be provided or a prob and a survey.
+
+    For the EB formulation, we start the derivation from Maxwell's equations:
+
+    .. math::
+        \\nabla \\times \\vec{E} + i \omega \\vec{B} = \\vec{s_m} \\\\
+        \\nabla \\times \\mu^{-1} \\vec{B} - \sigma \\vec{E} = \\vec{s_e}
+
+    we consider the physical properties, fields, and fluxes to be composed of
+    two parts, a primary and a secondary:
+
+        - :math:`\sigma   = \sigma_p + \sigma_s`
+        - :math:`\mu^{-1} = \mu^{-1}_p + \mu^{-1}_s`
+        - :math:`\\vec{E} = \\vec{E_p} + \\vec{E_s}`
+        - :math:`\\vec{B} = \\vec{B_p} + \\vec{B_s}`
+
+    and choose our primary such that
+
+    .. math::
+        \\nabla \\times \\vec{E}_p + i \omega \\vec{B}_p = \\vec{s_m} \\\\
+        \\nabla \\times \\mu^{-1}_p \\vec{B}_p - \sigma_p \\vec{E}_p = \\vec{s_e}_p
+
+    so the secondary problem is then
+
+    .. math::
+        \\nabla \\times \\vec{E}_s + i \omega \\vec{B}_s = 0 \\\\
+        \\nabla \\times \\mu^{-1} \\vec{B}_s - \sigma \\vec{E}_s = - \\nabla \\times \\mu^{-1}_s \\vec{B}_p + \sigma_s \\vec{E}_p
+
+
+    If instead, HJ formulation is considered, then we start off with
+
+    .. math::
+        \\nabla \\times \\rho \\vec{J} + i \omega \\mu \\vec{H} = \\vec{s_m} \\\\
+        \\nabla \\times \\vec{H} - \\vec{J} = \\vec{s_e}
+
+    and we define the primary secondary problem in terms of
+
+        - :math:`\\rho  = \\rho_p + \\rho_s`
+        - :math:`\mu = \mu_p + \mu_s`
+        - :math:`\\vec{J} = \\vec{J_p} + \\vec{J_s}`
+        - :math:`\\vec{H} = \\vec{H_p} + \\vec{H_s}`
+
+    with the primary being defined by
+
+    .. math::
+        \\nabla \\times \\rho_p \\vec{J}_p + i \omega \\mu_p \\vec{H}_p = \\vec{s_m} \\\\
+        \\nabla \\times \\vec{H}_p - \\vec{J}_p = \\vec{s_e}
+
+    so the secondary problem is given by
+
+    .. math::
+        \\nabla \\times \\rho \\vec{J}_s + i \omega \\mu \\vec{H} = - \\nabla \\times \\rho_s \\vec{J}_p - i \omega \\mu_s \\vec{H}_p \\
+        \\nabla \\times \\vec{H}_p - \\vec{J}_p = 0
+
+    Note: if different meshes are employed for the primary and secondary
+    problems, then we need to interpolate the fields from the primary mesh to
+    the secondary mesh. We do this by always interpolating the field and
+    computing a flux if need be in order to ensure that fluxes remain
+    numerically divergence free.
+
+    :param list rxList: Receiver list
+    :param float freq: frequency
+    :param numpy.array m: primary model
+    :param Problem prob: primary problem
+    :param Survey survey: primary survey
+    """
+
+
+    def __init__(self, rxList, freq, m, prob, survey):
+        self.freq = float(freq)
+        self.m = m
+        self.prob = prob
+        self.survey = survey
+        self.fields = None
+
+        if self.survey.ispaired:
+            if self.survey.prob is not self.prob:
+                raise Exception('The survey object is already paired to a problem. Use survey.unpair()')
+        else:
+            self.prob.pair(self.survey)
+
+        self.mesh = self.prob.mesh
+        self.prob.curModel = self.m
+
+        BaseSrc.__init__(self, rxList)
+
+    def MeSigma(self, prob):
+        if getattr(self, '_MeSigma', None) is None:
+            sigmaprimary = self.prob.curModel.sigma
+            if self.mesh != prob.mesh:
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                sigmaprimary = P * sigmaprimary
+            self._MeSigma = prob.mesh.getEdgeInnerProduct(sigmaprimary)
+        return self._MeSigma
+
+    def MfMui(self, prob):
+        if getattr(self, '_MfMui', None) is None:
+            muiprimary = self.prob.curModel.mui
+            if self.mesh != prob.mesh and not isinstance(muiprimary,float): # if different meshes and mu is a vector --> need to interpolate
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                muiprimary = P * muiprimary
+            self._MfMui = prob.mesh.getFaceInnerProduct(muiprimary)
+        return self._MfMui
+
+    def MfRho(self, prob):
+        if getattr(self, '_MfRho', None) is None:
+            rhoprimary = self.prob.curModel.rho
+            if self.mesh != prob.mesh:
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                rhoprimary = P * rhoprimary
+            self._MfRho = prob.mesh.getFaceInnerProduct(rhoprimary)
+        return self._MfRho
+
+    def MeMu(self, prob):
+        if getattr(self, '_MeMu', None) is None:
+            muprimary = self.prob.curModel.mu
+            if self.mesh != prob.mesh and not isinstance(muiprimary,float): # if different meshes and mu is a vector --> need to interpolate
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                muprimary = P * muprimary
+            self._MeMu = prob.mesh.getEdgeInnerProduct(muprimary)
+        return self._MeMu
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def ePrimary(self,prob):
+        if getattr(self, '_ePrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            ePrimary = self.fields[:,'e']
+
+            if self.mesh != prob.mesh:
+                if self.prob._formulation == 'HJ':
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('e'), locTypeFrom='CCV')
+                else:
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('e'))
+                ePrimary = Utils.mkvc(P * ePrimary)
+            self._ePrimary = Utils.mkvc(ePrimary)
+
+        return self._ePrimary
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def bPrimary(self, prob):
+        if getattr(self, '_bPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            if self.mesh == prob.mesh:
+                bPrimary = self.fields[:,'b']
+            else:
+                bPrimary = prob.mesh.edgeCurl * self.ePrimary(prob)
+
+            self._bPrimary = Utils.mkvc(bPrimary)
+
+        return self._bPrimary
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def hPrimary(self, prob):
+        if getattr(self, '_hPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            hPrimary = self.fields[:,'h']
+
+            if self.mesh != prob.mesh:
+                if self.prob._formulation == 'EB':
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('h'), locTypeFrom='CCV')
+                else:
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('h'))
+                print P.shape, hPrimary.shape, prob._GLoc('h')
+                hPrimary = Utils.mkvc(P * hPrimary)
+            self._hPrimary = Utils.mkvc(hPrimary)
+
+        return self._hPrimary
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def jPrimary(self, prob):
+        if getattr(self, '_jPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            if self.mesh == prob.mesh:
+                jPrimary = self.fields[:,'j']
+            else:
+                jPrimary = prob.mesh.edgeCurl * self.hPrimary(prob)
+
+            self._jPrimary = Utils.mkvc(jPrimary)
+
+        return self._jPrimary
+
+    def s_e(self,prob):
+        if prob._formulation == 'EB':
+            # - \\nabla \\times \\mu^{-1}_s \\vec{B}_p + \sigma_s \\vec{E}_p
+            s_e =  -prob.mesh.edgeCurl.T * ((prob.MfMui - self.MfMui(prob)) * self.bPrimary(prob)) +  (prob.MeSigma - self.MeSigma(prob)) * self.ePrimary(prob)
+            return Utils.mkvc(s_e)
+        else:
+            return Zero()
+
+    def s_eDeriv(self, prob, v, adjoint=False):
+        if prob._formulation == 'EB':
+            if adjoint is True:
+                return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
+            return prob.MeSigmaDeriv(self.ePrimary(prob)) * v
+        else:
+            return Zero()
+
+    def s_m(self,prob):
+        if prob._formulation == 'HJ':
+            # - \\nabla \\times \\rho_s \\vec{J}_p - i \omega \\mu_s \\vec{H}_p
+            s_m = - prob.mesh.edgeCurl.T * (prob.MfRho - self.MfRho(prob)) * self.jPrimary(prob) - 1j * omega(self.freq) * ((prob.MeMu - self.MeMu(prob)) * self.hPrimary(prob))
+            return s_m
+        else:
+            return Zero()
+
+    def s_mDeriv(self, prob, v, adjoint=False):
+        if prob._formulation == 'HJ':
+            if adjoint is True:
+                return - prob.MfRhoDeriv(self.jPrimary(prob)).T * (prob.mesh.edgeCurl * v)
+            return - prob.mesh.edgeCurl.T * (prob.MfRhoDeriv(self.jPrimary(prob)) * v)
+        else:
+            return Zero()
+
 
 

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -65,7 +65,7 @@ class Rx(SimPEG.Survey.BaseRx):
 
     def projGLoc(self, u):
         """Grid Location projection (e.g. Ex Fy ...)"""
-        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+        return u.prob._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
     def eval(self, src, mesh, f):
         """
@@ -77,8 +77,6 @@ class Rx(SimPEG.Survey.BaseRx):
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
-        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
-        # projGLoc += self.knownRxTypes[self.rxType][1]
 
         P = self.getP(mesh, self.projGLoc(f))
         f_part_complex = f[src, self.projField]

SimPEG/EM/Utils/testingUtils.py

@@ -20,7 +20,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
 
     if useMu is True:
-        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
     else:
         mapping = Maps.ExpMap(mesh)
 
@@ -37,20 +37,39 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
             Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'CircularLoop':
             Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+
         elif SrcType is 'RawVec':
             if fdemType is 'e' or fdemType is 'b':
                 S_m = np.zeros(mesh.nF)
                 S_e = np.zeros(mesh.nE)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e, integrate=True))
 
             elif fdemType is 'h' or fdemType is 'j':
                 S_m = np.zeros(mesh.nE)
                 S_e = np.zeros(mesh.nF)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e, integrate=True))
+
+        elif SrcType is 'PrimSec':
+                primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.])
+                primarySurvey = EM.FDEM.Survey([primSrc])
+                primaryProblem = EM.FDEM.Problem_e(mesh,mapping=mapping)
+                mPrimary = np.ones(mapping.nP)*np.log(CONDUCTIVITY)
+                Src.append(EM.FDEM.Src.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
+
+        elif SrcType is 'PrimSecCyl':
+                hx = [(cs,ncx + 2), (cs,npad + 2,1.3)]
+                hz = [(cs,npad + 2 ,-1.3), (cs,ncz+2), (cs,npad+2,1.3)]
+                primmesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+                primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.])
+                primarySurvey = EM.FDEM.Survey([primSrc])
+                primaryProblem = EM.FDEM.Problem_e(primmesh)
+                mPrimary = np.ones(primmesh.nC)*CONDUCTIVITY
+                Src.append(EM.FDEM.Src.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
@@ -90,7 +109,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
     print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
-    
+
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU
 

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -2,27 +2,19 @@ from SimPEG import Mesh, Utils, np, sp
 import SimPEG.DCIP as DC
 import time
 
-def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
+def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
     """
         DC Forward Simulation
         =====================
 
-        Forward model two conductive spheres in a half-space and plot a
-        pseudo-section. Assumes an infinite line source and measures along the
-        center of the spheres.
+        Forward model conductive spheres in a half-space and plot a pseudo-section
 
-        INPUT:
-        loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
-        radi    = Radius of spheres [r1,r2]
-        param   = Conductivity of background and two spheres [m0,m1,m2]
-        stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
-        dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
-        Created by @fourndo
+        Created by @fourndo on Mon Feb 01 19:28:06 2016
 
     """
 
     assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
-    assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)"
+
 
     if loc is None:
         loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
@@ -35,6 +27,7 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
 
 
     # First we need to create a mesh and a model.
+
     # This is our mesh
     dx    = 5.
 
@@ -59,10 +52,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
     # Get index of the center
     indy = int(mesh.nCy/2)
 
+
     # Plot the model for reference
     # Define core mesh extent
     xlim = 200
-    zlim = 100
+    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)]
@@ -80,13 +77,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
     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)
+    azm =  np.arctan(dl_y/dl_x)
 
     #Set boundary conditions
     mesh.setCellGradBC('neumann')
 
-    # Define the linear system needed for the DC problem. We assume an infitite
-    # line source for simplicity.
+    # Define the differential operators needed for the DC problem
     Div = mesh.faceDiv
     Grad = mesh.cellGrad
     Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
@@ -149,23 +145,16 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
     print 'Forward completed'
 
     # Let's just convert the 3D format into 2D (distance along line) and plot
-    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
+    # [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
+    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
     survey2D.dobs =np.hstack(data)
-
+    # Here is an example for the first tx-rx array
     if plotIt:
         import matplotlib.pyplot as plt
-        fig = plt.figure(figsize=(7,7))
+        fig = plt.figure()
         ax = plt.subplot(2,1,1, aspect='equal')
-        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax.add_artist(circle1)
-        ax.add_artist(circle2)
-
-        dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
-                             ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))
-
-        ax.set_title('3-D model')
+        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')
@@ -174,34 +163,22 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
         plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
 
 
-        pos =  ax.get_position()
-        ax.set_position([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height])
-        pos =  ax.get_position()
-        cbarax = fig.add_axes([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height * 0.04])  ## the parameters are the specified position you set
-        cb = fig.colorbar(dat[0],cax=cbarax, orientation="horizontal",
-                          ax = ax, ticks=np.linspace(np.log10(sig.min()),
-                         np.log10(sig.max()), 3), format="$10^{%.1f}$")
-        cb.set_label("Conductivity (S/m)",size=12)
-        cb.ax.tick_params(labelsize=12)
-
-        # Second plot for the predicted apparent resistivity data
-        ax2 = plt.subplot(2,1,2, aspect='equal')
+        ax = plt.subplot(2,1,2, aspect='equal')
 
         # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax2.add_artist(circle1)
-        ax2.add_artist(circle2)
+        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
-        dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)
+        DC.plot_pseudoSection(survey2D,ax,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')
-        ax2.set_title('Apparent Conductivity data')
-
         plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+
         plt.show()
 
         return fig, ax

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -48,7 +48,8 @@ def run(plotIt=True):
     freqs = np.logspace(1,3,10)
     srcLoc = np.array([0., 0., 10.])
 
-    srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]
+    srcList = []
+    [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
 
     survey = EM.FDEM.Survey(srcList)
     prb = EM.FDEM.Problem_b(mesh, mapping=mapping)

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -1,275 +0,0 @@
-from SimPEG import *
-from SimPEG.EM import FDEM, Analytics, mu_0
-import time
-
-try:
-    from pymatsolver import MumpsSolver
-    solver = MumpsSolver
-except Exception:
-    solver = SolverLU
-    pass
-
-def run(plotIt=True):
-    """
-        EM: Schenkel and Morrison Casing Model
-        ======================================
-
-        Here we create and run a FDEM forward simulation to calculate the vertical
-        current inside a steel-cased. The model is based on the Schenkel and
-        Morrison Casing Model, and the results are used in a 2016 SEG abstract by
-        Yang et al.
-
-        - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
-
-
-        The model consists of:
-        - Air: Conductivity 1e-8 S/m, above z = 0
-        - Background: conductivity 1e-2 S/m, below z = 0
-        - Casing: conductivity 1e6 S/m
-            - 300m long
-            - radius of 0.1m
-            - thickness of 6e-3m
-
-        Inside the casing, we take the same conductivity as the background.
-
-        We are using an EM code to simulate DC, so we use frequency low enough
-        that the skin depth inside the casing is longer than the casing length (f
-        = 1e-6 Hz). The plot produced is of the current inside the casing.
-
-        These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
-        resistivity modeling of steel casing for reservoir monitoring using
-        equivalent resistor network. The solver used to produce these results and
-        achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
-
-        .. _pymatsolver: https://github.com/rowanc1/pymatsolver
-
-        This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
-
-        If you would use this example for a code comparison, or build upon it, a
-        citation would be much appreciated!
-
-    """
-
-    if plotIt:
-        import matplotlib.pylab as plt
-
-    # ------------------ MODEL ------------------
-    sigmaair = 1e-8 # air
-    sigmaback = 1e-2 # background
-    sigmacasing = 1e6 # casing
-    sigmainside = sigmaback # inside the casing
-
-
-    casing_t = 0.006  # 1cm thickness
-    casing_l = 300   # length of the casing
-
-    casing_r = 0.1
-    casing_a = casing_r - casing_t/2. # inner radius
-    casing_b = casing_r + casing_t/2. # outer radius
-    casing_z = np.r_[-casing_l,0.]
-
-
-    # ------------------ SURVEY PARAMETERS ------------------
-    freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies
-    dsz = -300 # down-hole z source location
-    src_loc = np.r_[0.,0.,dsz]
-    inf_loc = np.r_[0.,0.,1e4]
-
-    print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
-
-
-    # ------------------ MESH ------------------
-    # fine cells near well bore
-    csx1, csx2 = 2e-3, 60.
-    pfx1, pfx2 = 1.3, 1.3
-    ncx1 = np.ceil(casing_b/csx1+2)
-
-    # pad nicely to second cell size
-    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
-    hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
-    dx1 = sum(hx1a)+sum(hx1b)
-    dx1 = np.floor(dx1/csx2)
-    hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
-
-    # second chunk of mesh
-    dx2 = 300. # uniform mesh out to here
-    ncx2 = np.ceil((dx2 - dx1)/csx2)
-    npadx2 = 45
-    hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
-    hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
-
-    # z-direction
-    csz = 0.05
-    nza = 10
-    ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
-    hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
-
-    # Mesh
-    mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
-
-    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
-    print 'Number of cells', mesh.nC
-
-    if plotIt is True:
-        fig, ax = plt.subplots(1, 1, figsize=(6, 4))
-        ax.set_title('Simulation Mesh')
-        mesh.plotGrid(ax=ax)
-        plt.show()
-
-    # Put the model on the mesh
-    sigWholespace = sigmaback*np.ones((mesh.nC))
-
-    sigBack = sigWholespace.copy()
-    sigBack[mesh.gridCC[:,2] > 0.] = sigmaair
-
-    sigCasing = sigBack.copy()
-    iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0])
-    iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b)
-    iCasing = iCasingX & iCasingZ
-    sigCasing[iCasing] = sigmacasing
-
-
-    if plotIt is True:
-
-        # plotting parameters
-        xlim = np.r_[0., 0.2]
-        zlim = np.r_[-350., 10.]
-        clim_sig = np.r_[-8,6]
-
-        # plot models
-        fig, ax = plt.subplots(1,1,figsize=(4,4))
-
-        f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax)
-        ax.grid(which='both')
-        ax.set_title('Log_10 (Sigma)')
-        ax.set_xlim(xlim)
-        ax.set_ylim(zlim)
-        f.set_clim(clim_sig)
-
-        plt.show()
-
-
-    # -------------- Sources --------------------
-    # Define Custom Current Sources
-
-    # surface source
-    sg_x = np.zeros(mesh.vnF[0],dtype=complex)
-    sg_y = np.zeros(mesh.vnF[1],dtype=complex)
-    sg_z = np.zeros(mesh.vnF[2],dtype=complex)
-
-    nza = 2 # put the wire two cells above the surface
-    ncin = 2
-
-    # vertically directed wire
-    sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1)  # hook it up to casing at the surface
-    sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
-    sgv_ind = sgv_indx & sgv_indz
-    sg_z[sgv_ind] = -1.
-
-    # horizontally directed wire
-    sgh_indx = (mesh.gridFx[:,0] > casing_a)    & (mesh.gridFx[:,0] <= inf_loc[2])
-    sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5))
-    sgh_ind = sgh_indx & sgh_indz
-    sg_x[sgh_ind] = -1.
-
-    sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2)  # hook it up to casing at the surface
-    sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
-    sgv2_ind = sgv2_indx & sgv2_indz
-    sg_z[sgv2_ind] = 1.
-
-    # assemble the source
-    sg = np.hstack([sg_x,sg_y,sg_z])
-    sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
-
-    # downhole source
-    dg_x = np.zeros(mesh.vnF[0],dtype=complex)
-    dg_y = np.zeros(mesh.vnF[1],dtype=complex)
-    dg_z = np.zeros(mesh.vnF[2],dtype=complex)
-
-    # vertically directed wire
-    dgv_indx = (mesh.gridFz[:,0] < csx1)  # go through the center of the well
-    dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
-    dgv_ind = dgv_indx & dgv_indz
-    dg_z[dgv_ind] = -1.
-
-    # couple to the casing downhole
-    dgh_indx = mesh.gridFx[:,0] < casing_a + csx1
-    dgh_indz = (mesh.gridFx[:,2] < dsz + csz)  & (mesh.gridFx[:,2] >= dsz)
-    dgh_ind = dgh_indx & dgh_indz
-    dg_x[dgh_ind] = 1.
-
-    # horizontal part at surface
-    dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2
-    dgh2_indz = sgh_indz.copy()
-    dgh2_ind = dgh2_indx & dgh2_indz
-    dg_x[dgh2_ind] = -1.
-
-    # vertical part at surface
-    dgv2_ind = sgv2_ind.copy()
-    dg_z[dgv2_ind] = 1.
-
-    # assemble the source
-    dg = np.hstack([dg_x,dg_y,dg_z])
-    dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
-
-    # ------------ Problem and Survey ---------------
-    survey = FDEM.Survey(sg_p + dg_p)
-    mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem_h(mesh, mapping=mapping)
-    problem.pair(survey)
-
-    # ------------- Solve ---------------------------
-    t0 = time.time()
-    fieldsCasing = problem.fields(sigCasing)
-    print 'Time to solve 2 sources', time.time() - t0
-
-    # Plot current
-
-    # current density
-    jn0 = fieldsCasing[dg_p,'j']
-    jn1 = fieldsCasing[sg_p,'j']
-
-    # current
-    in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))]
-    in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))]
-
-    in0 = np.vstack(in0).T
-    in1 = np.vstack(in1).T
-
-    # integrate to get z-current inside casing
-    inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b)
-    inds_inz = (mesh.gridFz[:,2] >= dsz )     & (mesh.gridFz[:,2] <= 0)
-    inds_fz = inds_inx & inds_inz
-
-    indsx = [False]*mesh.nFx
-    inds = list(indsx) + list(inds_fz)
-
-    in0_in = in0[np.r_[inds]]
-    in1_in = in1[np.r_[inds]]
-    z_in = mesh.gridFz[inds_fz,2]
-
-    in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
-    in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
-    z_in = z_in.reshape([z_in.shape[0]/3,3])
-
-    I0 = in0_in.sum(1).real
-    I1 = in1_in.sum(1).real
-    z_in = z_in[:,0]
-
-    if plotIt is True:
-        fig, ax = plt.subplots(1,2,figsize=(12,4))
-
-        ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1))
-        ax[0].legend(['top casing', 'bottom casing'],loc='best')
-        ax[0].set_title('Magnitude of Vertical Current in Casing')
-
-        ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1))
-        ax[1].legend(['top casing', 'bottom casing'],loc='best')
-        ax[1].set_title('Magnitude of Vertical Current in Casing')
-        ax[1].set_ylim([1e-2, 1.])
-
-        plt.show()
-
-if __name__ == '__main__':
-    run()
-

SimPEG/Examples/Inversion_IRLS.py

@@ -1,132 +0,0 @@
-from SimPEG import *
-
-
-def run(N=200, plotIt=True):
-    """
-        Inversion: Linear Problem
-        =========================
-
-        Here we go over the basics of creating a linear problem and inversion.
-
-    """
-
-
-    np.random.seed(1)
-
-    std_noise = 1e-2
-
-    mesh = Mesh.TensorMesh([N])
-
-    m0 = np.ones(mesh.nC) * 1e-4
-    nk = 10
-    jk = np.linspace(1.,nk,nk)
-    p = -2.
-    q = 1.
-
-    g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
-
-    G = np.empty((nk, mesh.nC))
-
-    for i in range(nk):
-        G[i,:] = g(i)
-
-    mtrue = np.zeros(mesh.nC)
-    mtrue[mesh.vectorCCx > 0.3] = 1.
-    mtrue[mesh.vectorCCx > 0.45] = -0.5
-    mtrue[mesh.vectorCCx > 0.6] = 0
-
-
-    prob = Problem.LinearProblem(mesh, G)
-    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.wght = wr
-
-    dmis = DataMisfit.l2_DataMisfit(survey)
-    dmis.Wd = 1./wd
-
-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, 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)
-
-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
-    reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
-
-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
-
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
-
-    # Run inversion
-    mrec = inv.run(m0)
-
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-
-    if plotIt:
-        import matplotlib.pyplot as plt
-
-        fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-        for i in range(prob.G.shape[0]):
-            axes[0].plot(prob.G[i,:])
-        axes[0].set_title('Columns of matrix G')
-
-        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
-        #axes[1].legend(('True Model', 'Recovered Model'))
-        axes[1].set_ylim(-1.0,1.25)
-
-        axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2)
-        axes[1].legend(('True Model', 'Smooth l2-l2',
-        'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12)
-        plt.show()
-
-    return prob, survey, mesh, mrec
-
-if __name__ == '__main__':
-    run()

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -100,7 +100,7 @@ def run(plotIt=True):
     # Regularization - with a regularization mesh
     regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
-    reg.mrefInSmooth = True
+    reg.smoothModel = True
     reg.alpha_s = 1e-7
     reg.alpha_x = 1.
     # Inversion problem

SimPEG/Examples/__init__.py

@@ -5,11 +5,9 @@ import DC_Analytic_Dipole
 import DC_Forward_PseudoSection
 import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
-import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
-import Inversion_IRLS
 import Inversion_Linear
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
@@ -21,7 +19,7 @@ import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 
-__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", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "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"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "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"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Inversion.py

@@ -33,9 +33,7 @@ class BaseInversion(object):
         self._directiveList = value
         self._directiveList.inversion = self
 
-    def __init__(self, invProb, directiveList=None, **kwargs):
-        if directiveList is None:
-            directiveList = []
+    def __init__(self, invProb, directiveList=[], **kwargs):
         self.directiveList = directiveList
         Utils.setKwargs(self, **kwargs)
 

SimPEG/MT/Utils/ediFilesUtils.py

@@ -7,16 +7,17 @@ from SimPEG.MT.Utils.dataUtils import rec2ndarr
 # Import modules
 import numpy as np
 import os, sys, re
-
+try:
+    import osr
+except ImportError as e:
+    print 'Could not import osr, missing the gdal package'
+    pass
 
 class EDIimporter:
     """
     A class to import EDIfiles.
 
     """
-
-
-    # Define data converters
     _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
     _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
 
@@ -25,8 +26,8 @@ class EDIimporter:
     comps = None
 
     # Hidden properties
-    _outEPSG = None # Project info
-    _2out = None    # The projection operator
+    _outEPSG = None
+    _2out = None
 
 
     def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -112,12 +113,6 @@ class EDIimporter:
     # nOutData=length(obj.data);
     # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
     def _transfromPoints(self,longD,latD):
-        # Import the coordinate projections
-        try:
-            import osr
-        except ImportError as e:
-            print 'Could not import osr, missing the gdal package\nCan not project coordinates'
-            raise e
         # Coordinates convertor
         if self._2out is None:
             src = osr.SpatialReference()

SimPEG/Maps.py

@@ -759,29 +759,15 @@ class PolyMap(IdentityMap):
 
             m = [\sigma_1, \sigma_2, c]
 
-        Can take in an actInd vector to account for topography.
-
     """
-    def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
+    def __init__(self, mesh, order, logSigma=True, normal='X'):
         IdentityMap.__init__(self, mesh)
         self.logSigma = logSigma
         self.order = order
         self.normal = normal
-        self.actInd = actInd
-
-        if getattr(self, 'actInd', None) is None:
-            self.actInd = range(self.mesh.nC)
-            self.nC = self.mesh.nC
-
-        else:
-            self.nC = len(self.actInd)
 
     slope = 1e4
 
-    @property
-    def shape(self):
-        return (self.nC, self.nP)
-
     @property
     def nP(self):
         if np.isscalar(self.order):
@@ -799,8 +785,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
             elif self.normal =='Y':
@@ -809,9 +795,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
             elif self.normal =='Y':
@@ -820,7 +806,6 @@ class PolyMap(IdentityMap):
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
-
         else:
             raise(Exception("Only supports 2D"))
 
@@ -834,8 +819,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
 
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
@@ -847,9 +832,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
 
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X

SimPEG/Mesh/BaseMesh.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.sparse as sp
 from SimPEG import Utils
 
 
@@ -594,3 +595,63 @@ class BaseRectangularMesh(BaseMesh):
             return out
         else:
             return switchKernal(x)
+
+
+    def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC', locTypeFrom=None):
+        """
+        Interpolates variables from the current mesh to a new mesh (mesh2)
+
+        :param Mesh mesh2: SimPEG mesh which we interpolate values to
+        :param string locType: location of variables 'CC', 'E', 'F', 'N'
+        :rtype: scipy.sparse.csr_matrix
+        :return P: interpolation matrix
+        """
+
+        # import warnings
+        # warnings.warn(
+        #      "`getInterpolationMatMesh2Mesh` will be slow. If you want to interpolate a vector from one mesh to another, use `InterpolateVecMesh2Mesh`",
+        #      RuntimeWarning)
+
+        if locTypeFrom is None:
+            locTypeFrom = locType # assume that we are interpolating to and from the same place
+
+        # Error Checking
+        if self._meshType == 'CYL':
+            assert self.isSymmetric, "Currently, we do not support non-symmetric cyl meshes"
+        if mesh2._meshType == 'CYL':
+            assert self._meshType == 'CYL', "Interpolation from 3D mesh to Cyl mesh is not supported"
+
+        # if Cyl to cart call
+        if self._meshType == 'CYL' and mesh2._meshType != 'CYL':
+            return self.getInterpolationMatCartMesh(mesh2, locType)
+
+        # Scalars
+        if locType in ['CC', 'CCVx', 'CCVy', 'CCVz', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
+            grid = getattr(mesh2, 'grid%s'%locTypeFrom)
+            return self.getInterpolationMat(grid, locType)
+
+        # Vectors
+        else:
+            if self._meshType == 'CYL':
+                if locType == 'F':
+                    X = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fx', locTypeFrom=locTypeFrom+'x')
+                    Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fz', locTypeFrom=locTypeFrom+'z')
+                    return sp.block_diag([X, Z])
+                elif locType == 'E':
+                    return self.getInterpolationMatMesh2Mesh(mesh2, locType='Ey', locTypeFrom=locTypeFrom+'y')
+
+            if self.dim == 1:
+                return self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+            elif self.dim == 2:
+                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType, locTypeFrom=locTypeFrom+'y')
+                return sp.block_diag([X, Y])
+            elif self.dim == 3:
+                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType, locTypeFrom=locTypeFrom+'y')
+                Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sz'%locType, locTypeFrom=locTypeFrom+'z')
+                return sp.block_diag([X, Y, Z])
+
+
+
+

SimPEG/Mesh/DiffOperators.py

@@ -307,28 +307,24 @@ class DiffOperators(object):
         return BC
     _cellGradBC_list = 'neumann'
 
-    def _cellGradStencil(self):
-        BC = self.setCellGradBC(self._cellGradBC_list)
-        n = self.vnC
-        if(self.dim == 1):
-            G = ddxCellGrad(n[0], BC[0])
-        elif(self.dim == 2):
-            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
-            G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
-            G = sp.vstack((G1, G2), format="csr")
-        elif(self.dim == 3):
-            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
-            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
-            G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
-            G = sp.vstack((G1, G2, G3), format="csr")
-        return G
-
     def cellGrad():
         doc = "The cell centered Gradient, takes you to cell faces."
 
         def fget(self):
             if(self._cellGrad is None):
-                G = self._cellGradStencil()
+                BC = self.setCellGradBC(self._cellGradBC_list)
+                n = self.vnC
+                if(self.dim == 1):
+                    G = ddxCellGrad(n[0], BC[0])
+                elif(self.dim == 2):
+                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
+                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
+                    G = sp.vstack((G1, G2), format="csr")
+                elif(self.dim == 3):
+                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
+                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
+                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
+                    G = sp.vstack((G1, G2, G3), format="csr")
                 # Compute areas of cell faces & volumes
                 S = self.area
                 V = self.aveCC2F*self.vol  # Average volume between adjacent cells
@@ -365,24 +361,19 @@ class DiffOperators(object):
     _cellGradBC = None
     cellGradBC = property(**cellGradBC())
 
-    def _cellGradxStencil(self):
-        BC = ['neumann', 'neumann']
-        n = self.vnC
-        if(self.dim == 1):
-            G1 = ddxCellGrad(n[0], BC)
-        elif(self.dim == 2):
-            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
-        elif(self.dim == 3):
-            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
-        return G1
-
-
     def cellGradx():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
 
         def fget(self):
             if getattr(self, '_cellGradx', None) is None:
-                G1 = self._cellGradxStencil()
+                BC = ['neumann', 'neumann']
+                n = self.vnC
+                if(self.dim == 1):
+                    G1 = ddxCellGrad(n[0], BC)
+                elif(self.dim == 2):
+                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
+                elif(self.dim == 3):
+                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fx', 'V')
@@ -391,22 +382,17 @@ class DiffOperators(object):
         return locals()
     cellGradx = property(**cellGradx())
 
-    def _cellGradyStencil(self):
-        if self.dim < 2: return None
-        BC = ['neumann', 'neumann']
-        n = self.vnC
-        if(self.dim == 2):
-            G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
-        elif(self.dim == 3):
-            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
-        return G2
-
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 2: return None
             if getattr(self, '_cellGrady', None) is None:
-                G2 = self._cellGradyStencil()
+                BC = ['neumann', 'neumann']
+                n = self.vnC
+                if(self.dim == 2):
+                    G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
+                elif(self.dim == 3):
+                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fy', 'V')
@@ -415,19 +401,14 @@ class DiffOperators(object):
         return locals()
     cellGrady = property(**cellGrady())
 
-    def _cellGradzStencil(self):
-        if self.dim < 3: return None
-        BC = ['neumann', 'neumann']
-        n = self.vnC
-        G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
-        return G3
-
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 3: return None
             if getattr(self, '_cellGradz', None) is None:
-                G3 = self._cellGradzStencil()
+                BC = ['neumann', 'neumann']
+                n = self.vnC
+                G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fz', 'V')

SimPEG/Mesh/MeshIO.py

@@ -21,9 +21,10 @@ class TensorMeshIO(object):
                 if '*' in seg:
                     st = seg
                     sp = seg.split('*')
-                    re = int(sp[0])*(' ' + sp[1])
+                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
                     line = line.replace(st,re.strip())
             return np.array(line.split(),dtype=float)
+
         # Read the file as line strings, remove lines with comment = !
         msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
 

SimPEG/Mesh/TreeMesh.py

@@ -2131,16 +2131,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
     def plotSlice(self, v, vType='CC',
         normal='Z', ind=None, grid=True, view='real',
         ax=None, clim=None, showIt=False,
-        pcolorOpts=None,
-        streamOpts=None,
-        gridOpts=None):
+        pcolorOpts={},
+        streamOpts={'color':'k'},
+        gridOpts={'color':'k', 'alpha':0.5}):
 
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k', 'alpha':0.5}
         assert vType in ['CC','F','E']
         assert self.dim == 3
 

SimPEG/Mesh/View.py

@@ -42,9 +42,9 @@ class TensorView(object):
 
     def plotImage(self, v, vType='CC', grid=False, view='real',
               ax=None, clim=None, showIt=False,
-              pcolorOpts=None,
-              streamOpts=None,
-              gridOpts=None,
+              pcolorOpts={},
+              streamOpts={'color':'k'},
+              gridOpts={'color':'k'},
               numbering=True, annotationColor='w'
               ):
         """
@@ -84,12 +84,6 @@ class TensorView(object):
             M.plotImage(v, annotationColor='k', showIt=True)
 
         """
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k'}
 
         if ax is None:
             fig = plt.figure()
@@ -180,9 +174,9 @@ class TensorView(object):
     def plotSlice(self, v, vType='CC',
                   normal='Z', ind=None, grid=False, view='real',
                   ax=None, clim=None, showIt=False,
-                  pcolorOpts=None,
-                  streamOpts=None,
-                  gridOpts=None
+                  pcolorOpts={},
+                  streamOpts={'color':'k'},
+                  gridOpts={'color':'k', 'alpha':0.5}
                   ):
 
         """
@@ -203,12 +197,6 @@ class TensorView(object):
             M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
 
         """
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k', 'alpha':0.5}
         if type(vType) in [list, tuple]:
             assert ax is None, "cannot specify an axis to plot on with this function."
             fig, axs = plt.subplots(1,len(vType))
@@ -301,17 +289,11 @@ class TensorView(object):
 
     def _plotImage2D(self, v, vType='CC', grid=False, view='real',
               ax=None, clim=None, showIt=False,
-              pcolorOpts=None,
-              streamOpts=None,
-              gridOpts=None
+              pcolorOpts={},
+              streamOpts={'color':'k'},
+              gridOpts={'color':'k'}
               ):
 
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k'}
         vTypeOptsCC = ['N','CC','Fx','Fy','Ex','Ey']
         vTypeOptsV = ['CCv','F','E']
         vTypeOpts = vTypeOptsCC + vTypeOptsV

SimPEG/Optimization.py

@@ -888,8 +888,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
     maxIterCG = 5
     tolCG = 1e-1
 
-    stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds
-
     lower = -np.inf
     upper = np.inf
 
@@ -992,20 +990,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                     cgFlag = 1
                 # End CG Iterations
 
-            # Take a gradient step on the active cells if exist
-            if temp != self.xc.size:
-
-                rhs_a = (Active) * -self.g
-
-                dm_i = max( abs( delx ) )
-                dm_a = max( abs(rhs_a) )
-
-                # perturb inactive set off of bounds so that they are included in the step
-                delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
-
-
-            # Only keep gradients going in the right direction on the active set
-            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
-            delx[indx] = 0.
-
             return delx

SimPEG/Regularization.py

@@ -1,289 +1,5 @@
 import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
 
-class RegularizationMesh(object):
-    """
-    **Regularization Mesh**
-
-    This contains the operators used in the regularization. Note that these
-    are not necessarily true differential operators, but are constructed from
-    a SimPEG Mesh.
-
-    :param Mesh mesh: problem mesh
-    :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
-    """
-
-    def __init__(self, mesh, indActive=None):
-        self.mesh = mesh
-        assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool'
-        self.indActive = indActive
-
-    @property
-    def vol(self):
-        """
-        reduced volume vector
-        :rtype: numpy.array
-        :return: reduced cell volume
-        """
-        if getattr(self, '_vol', None) is None:
-            self._vol = self._Pac.T * self.mesh.vol
-        return self._vol
-
-    @property
-    def nC(self):
-        """
-        reduced number of cells
-        :rtype: int
-        :return: number of cells being regularized
-        """
-        if getattr(self, '_nC', None) is None:
-            if self.indActive is None:
-                self._nC = self.mesh.nC
-            else:
-                self._nC = sum(self.indActive)
-        return self._nC
-
-    @property
-    def dim(self):
-        """
-        dimension of regularization mesh (1D, 2D, 3D)
-        :rtype: int
-        :return: dimension
-        """
-        if getattr(self, '_dim', None) is None:
-            self._dim = self.mesh.dim
-        return self._dim
-
-
-    @property
-    def _Pac(self):
-        """
-        projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)
-        :rtype: scipy.sparse.csr_matrix
-        :return: active cell projection matrix
-        """
-        if getattr(self, '__Pac', None) is None:
-            if self.indActive is None:
-                self.__Pac = Utils.speye(self.mesh.nC)
-            else:
-                self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-        return self.__Pac
-
-    @property
-    def _Pafx(self):
-        """
-        projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )
-        :rtype: scipy.sparse.csr_matrix
-        :return: active face-x projection matrix
-        """
-        if getattr(self, '__Pafx', None) is None:
-            if self.indActive is None:
-                self.__Pafx = Utils.speye(self.mesh.nFx)
-            else:
-                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
-                self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
-        return self.__Pafx
-
-    @property
-    def _Pafy(self):
-        """
-        projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )
-        :rtype: scipy.sparse.csr_matrix
-        :return: active face-y projection matrix
-        """
-        if getattr(self, '__Pafy', None) is None:
-            if self.indActive is None:
-                self.__Pafy = Utils.speye(self.mesh.nFy)
-            else:
-                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
-                self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
-        return self.__Pafy
-
-    @property
-    def _Pafz(self):
-        """
-        projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )
-        :rtype: scipy.sparse.csr_matrix
-        :return: active face-z projection matrix
-        """
-        if getattr(self, '__Pafz', None) is None:
-            if self.indActive is None:
-                self.__Pafz = Utils.speye(self.mesh.nFz)
-            else:
-                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
-                self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
-        return self.__Pafz
-
-    @property
-    def aveFx2CC(self):
-        """
-        averaging from active cell centers to active x-faces
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging from active cell centers to active x-faces
-        """
-        if getattr(self, '_aveFx2CC', None) is None:
-            self._aveFx2CC =  self._Pac.T * self.mesh.aveFx2CC * self._Pafx
-        return self._aveFx2CC
-
-    @property
-    def aveCC2Fx(self):
-        """
-        averaging from active x-faces to active cell centers
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging matrix from active x-faces to active cell centers
-        """
-        if getattr(self, '_aveCC2Fx', None) is None:
-            self._aveCC2Fx =  Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
-        return self._aveCC2Fx
-
-    @property
-    def aveFy2CC(self):
-        """
-        averaging from active cell centers to active y-faces
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging from active cell centers to active y-faces
-        """
-        if getattr(self, '_aveFy2CC', None) is None:
-            self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
-        return self._aveFy2CC
-
-    @property
-    def aveCC2Fy(self):
-        """
-        averaging from active y-faces to active cell centers
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging matrix from active y-faces to active cell centers
-        """
-        if getattr(self, '_aveCC2Fy', None) is None:
-            self._aveCC2Fy =  Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
-        return self._aveCC2Fy
-
-    @property
-    def aveFz2CC(self):
-        """
-        averaging from active cell centers to active z-faces
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging from active cell centers to active z-faces
-        """
-        if getattr(self, '_aveFz2CC', None) is None:
-            self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
-        return self._aveFz2CC
-
-    @property
-    def aveCC2Fz(self):
-        """
-        averaging from active z-faces to active cell centers
-        :rtype: scipy.sparse.csr_matrix
-        :return: averaging matrix from active z-faces to active cell centers
-        """
-        if getattr(self, '_aveCC2Fz', None) is None:
-            self._aveCC2Fz =  Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
-        return self._aveCC2Fz
-
-    @property
-    def cellDiffx(self):
-        """
-        cell centered difference in the x-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the x-direction
-        """
-        if getattr(self, '_cellDiffx', None) is None:
-            self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac
-        return self._cellDiffx
-
-    @property
-    def cellDiffy(self):
-        """
-        cell centered difference in the y-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the y-direction
-        """
-        if getattr(self, '_cellDiffy', None) is None:
-            self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac
-        return self._cellDiffy
-
-    @property
-    def cellDiffz(self):
-        """
-        cell centered difference in the z-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the z-direction
-        """
-        if getattr(self, '_cellDiffz', None) is None:
-            self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac
-        return self._cellDiffz
-
-    @property
-    def faceDiffx(self):
-        """
-        x-face differences
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active faces in the x-direction
-        """
-        if getattr(self, '_faceDiffx', None) is None:
-            self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx
-        return self._faceDiffx
-
-    @property
-    def faceDiffy(self):
-        """
-        y-face differences
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active faces in the y-direction
-        """
-        if getattr(self, '_faceDiffy', None) is None:
-            self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy
-        return self._faceDiffy
-
-    @property
-    def faceDiffz(self):
-        """
-        z-face differences
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active faces in the z-direction
-        """
-        if getattr(self, '_faceDiffz', None) is None:
-            self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz
-        return self._faceDiffz
-
-    @property
-    def cellDiffxStencil(self):
-        """
-        cell centered difference stencil (no cell lengths include) in the x-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the x-direction
-        """
-        if getattr(self, '_cellDiffxStencil', None) is None:
-
-            self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
-        return self._cellDiffxStencil
-
-    @property
-    def cellDiffyStencil(self):
-        """
-        cell centered difference stencil (no cell lengths include) in the y-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the y-direction
-        """
-        if self.dim < 2: return None
-        if getattr(self, '_cellDiffyStencil', None) is None:
-
-            self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
-        return self._cellDiffyStencil
-
-    @property
-    def cellDiffzStencil(self):
-        """
-        cell centered difference stencil (no cell lengths include) in the y-direction
-        :rtype: scipy.sparse.csr_matrix
-        :return: differencing matrix for active cells in the y-direction
-        """
-        if self.dim < 3: return None
-        if getattr(self, '_cellDiffzStencil', None) is None:
-
-            self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
-        return self._cellDiffzStencil
-
-
 class BaseRegularization(object):
     """
     **Base Regularization Class**
@@ -302,16 +18,12 @@ class BaseRegularization(object):
 
     mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
-    mref    = None    #: Reference model.
+    mref = None       #: Reference model.
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
+        self.mesh = mesh
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        if indActive is not None and indActive.dtype != 'bool':
-            tmp = indActive
-            indActive = np.zeros(mesh.nC, dtype=bool)
-            indActive[tmp] = True
-        self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
@@ -343,7 +55,8 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.regmesh.nC)
+        return sp.identity(self.mapping.nP)
+
 
     @Utils.timeIt
     def eval(self, m):
@@ -374,12 +87,11 @@ class BaseRegularization(object):
     @Utils.timeIt
     def eval2Deriv(self, m, v=None):
         """
-        Second derivative
 
-        :param numpy.array m: geophysical model
-        :param numpy.array v: vector to multiply
-        :rtype: scipy.sparse.csr_matrix or numpy.ndarray
-        :return: WtW or WtW*v
+            :param numpy.array m: geophysical model
+            :param numpy.array v: vector to multiply
+            :rtype: scipy.sparse.csr_matrix or numpy.ndarray
+            :return: WtW or WtW*v
 
         The regularization is:
 
@@ -400,94 +112,112 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
+
 class Tikhonov(BaseRegularization):
     """
-    L2 Tikhonov regularization with both smallness and smoothness (first order
-    derivative) contributions.
-
-    .. math::
-        \phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-
-    Note if the key word argument `mrefInSmooth` is False, then mref is not
-    included in the smoothness contribution.
-
-    :param Mesh mesh: SimPEG mesh
-    :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
-    :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
-    :param bool mrefInSmooth: (default = False) put mref in the smoothness component?
-    :param float alpha_s: (default 1e-6) smallness weight
-    :param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction
-    :param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction
-    :param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction
-    :param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction
-    :param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction
-    :param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction
     """
-    mrefInSmooth = False  # put mref in the smoothness contribution
-    alpha_s      = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x      = Utils.dependentProperty('_alpha_x',  1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
-    alpha_y      = Utils.dependentProperty('_alpha_y',  1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
-    alpha_z      = Utils.dependentProperty('_alpha_z',  1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
-    alpha_xx     = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'],    "Weight for the second derivative in the x direction")
-    alpha_yy     = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'],    "Weight for the second derivative in the y direction")
-    alpha_zz     = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'],    "Weight for the second derivative in the z direction")
+    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
+    alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+    alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+    alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
+    alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
+    alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
     def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+        self.indActive = indActive
 
     @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
-        return self._Wsmall
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Ws = Pac.T * self._Ws * Pac
+        return self._Ws
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
-            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx
+            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
+
+            if self.indActive is not None:
+                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
+                self._Wx = Pafx.T*self._Wx*Pac
+
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
-            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy
+            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
+
+            if self.indActive is not None:
+                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
+                self._Wy = Pafy.T*self._Wy*Pac
+
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
-            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz
+            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
+
+            if self.indActive is not None:
+                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
+                self._Wz = Pafz.T*self._Wz*Pac
+
         return self._Wz
 
     @property
     def Wxx(self):
         """Regularization matrix Wxx"""
         if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx
+            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wxx = Pac.T*self._Wxx*Pac
+
         return self._Wxx
 
     @property
     def Wyy(self):
         """Regularization matrix Wyy"""
         if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy
+            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wyy = Pac.T*self._Wyy*Pac
+
         return self._Wyy
 
     @property
     def Wzz(self):
         """Regularization matrix Wzz"""
         if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
+            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
+
+            if self.indActive is not None:
+                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+                self._Wzz = Pac.T*self._Wzz*Pac
+
         return self._Wzz
 
     @property
@@ -495,9 +225,9 @@ class Tikhonov(BaseRegularization):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
             wlist = (self.Wx, self.Wxx)
-            if self.regmesh.dim > 1:
+            if self.mesh.dim > 1:
                 wlist += (self.Wy, self.Wyy)
-            if self.regmesh.dim > 2:
+            if self.mesh.dim > 2:
                 wlist += (self.Wz, self.Wzz)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
@@ -506,44 +236,25 @@ class Tikhonov(BaseRegularization):
     def W(self):
         """Full regularization matrix W"""
         if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
+            wlist = (self.Ws, self.Wsmooth)
             self._W = sp.vstack(wlist)
         return self._W
 
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
-        return 0.5 * r.dot(r)
-
     @Utils.timeIt
     def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
+        if self.smoothModel == True:
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.smoothModel == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
 
-    @Utils.timeIt
-    def _evalSmallDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self,m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
 
     @Utils.timeIt
     def evalDeriv(self, m):
         """
+
         The regularization is:
 
         .. math::
@@ -557,185 +268,257 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        if self.smoothModel == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.smoothModel == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
+
+# <<<<<<< HEAD
+
+# class Simple(BaseRegularization):
+#     """
+#         Only for tensor mesh
+#     """
+
+#     smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
+#     alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
+#     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+#     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+#     alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+#     alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
+#     alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
+#     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
+
+#     def __init__(self, mesh, mapping=None, **kwargs):
+#         BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
 
 
-class Simple(Tikhonov):
-    """
-    Simple regularization that does not include length scales in the derivatives.
-    """
 
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    wght        = 1.
+#     @property
+#     def Ws(self):
+#         """Regularization matrix Ws"""
+#         if getattr(self,'_Ws', None) is None:
+#                 self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+#         return self._Ws
 
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+#     @property
+#     def Wx(self):
+#         """Regularization matrix Wx"""
+#         if getattr(self, '_Wx', None) is None:
+#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
+#         return self._Wx
 
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+#     @property
+#     def Wy(self):
+#         """Regularization matrix Wy"""
+#         if getattr(self, '_Wy', None) is None:
+#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
+#         return self._Wy
 
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
-        return self._Wsmall
+#     @property
+#     def Wz(self):
+#         """Regularization matrix Wz"""
+#         if getattr(self, '_Wz', None) is None:
+#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
+#         return self._Wz
 
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
-        return self._Wx
+#     @property
+#     def Wxx(self):
+#         """Regularization matrix Wxx"""
+#         if getattr(self, '_Wxx', None) is None:
+#             self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
+#         return self._Wxx
 
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
-        return self._Wy
+#     @property
+#     def Wyy(self):
+#         """Regularization matrix Wyy"""
+#         if getattr(self, '_Wyy', None) is None:
+#             self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
+#         return self._Wyy
 
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
-        return self._Wz
+#     @property
+#     def Wzz(self):
+#         """Regularization matrix Wzz"""
+#         if getattr(self, '_Wzz', None) is None:
+#             self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
+#         return self._Wzz
 
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx,)
-            if self.regmesh.dim > 1:
-                wlist += (self.Wy,)
-            if self.regmesh.dim > 2:
-                wlist += (self.Wz,)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
+#     @property
+#     def Wsmooth(self):
+#         """Full smoothness regularization matrix W"""
+#         if getattr(self, '_Wsmooth', None) is None:
+#             wlist = (self.Wx, self.Wxx)
+#             if self.mesh.dim > 1:
+#                 wlist += (self.Wy, self.Wyy)
+#             if self.mesh.dim > 2:
+#                 wlist += (self.Wz, self.Wzz)
+#             self._Wsmooth = sp.vstack(wlist)
+#         return self._Wsmooth
 
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
+#     @property
+#     def W(self):
+#         """Full regularization matrix W"""
+#         if getattr(self, '_W', None) is None:
+#             wlist = (self.Ws, self.Wsmooth)
+#             self._W = sp.vstack(wlist)
+#         return self._W
 
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
+#     @Utils.timeIt
+#     def eval(self, m):
+#         if self.smoothModel == True:
+#             r1 = self.Wsmooth * ( self.mapping * (m) )
+#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
+#             return 0.5*(r1.dot(r1)+r2.dot(r2))
+#         elif self.smoothModel == False:
+#             r = self.W * ( self.mapping * (m - self.mref) )
+#             return 0.5*r.dot(r)
 
 
-class Sparse(Simple):
+#     @Utils.timeIt
+#     def evalDeriv(self, m):
+#         """
 
-    # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
-    curModel = None # use a model to compute the weights
-    gamma    = 1.
-    norms    = [0., 2., 2., 2.]
-    wght     = 1.
+#         The regularization is:
 
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+#         .. math::
 
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+#             R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
 
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self, 'curModel', None) is None:
-            self.Rs = Utils.speye(self.regmesh.nC)
+#         So the derivative is straight forward:
 
-        else:
-            f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
+#         .. math::
 
-        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
+#             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+#         """
+#         if self.smoothModel == True:
+#             mD1 = self.mapping.deriv(m)
+#             mD2 = self.mapping.deriv(m - self.mref)
+#             r1 = self.Wsmooth * ( self.mapping * (m))
+#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
+#             out1 = mD1.T * ( self.Wsmooth.T * r1 )
+#             out2 = mD2.T * ( self.Ws.T * r2 )
+#             out = out1+out2
+#         elif self.smoothModel == False:
+#             mD = self.mapping.deriv(m - self.mref)
+#             r = self.W * ( self.mapping * (m - self.mref) )
+#             out = mD.T * ( self.W.T * r )
+#         return out
+
+# class SparseRegularization(Simple):
+
+#     eps = 1e-1
+
+#     m = None
+#     gamma = 1.
+#     p = 0.
+#     qx = 2.
+#     qy = 2.
+#     qz = 2.
+
+#     def __init__(self, mesh, mapping=None, **kwargs):
+#         Simple.__init__(self, mesh, mapping=mapping, **kwargs)
 
 
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
+#     @property
+#     def Wsmooth(self):
+#         """Full smoothness regularization matrix W"""
+#         if getattr(self, '_Wsmooth', None) is None:
+#             wlist = (self.Wx, self.Wxx)
+#             if self.mesh.dim > 1:
+#                 wlist += (self.Wy, self.Wyy)
+#             if self.mesh.dim > 2:
+#                 wlist += (self.Wz, self.Wzz)
+#             self._Wsmooth = sp.vstack(wlist)
+#         return self._Wsmooth
 
-        if getattr(self, 'curModel', None) is None:
-            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+#     @property
+#     def W(self):
+#         """Full regularization matrix W"""
+#         if getattr(self, '_W', None) is None:
+#             wlist = (self.Ws, self.Wsmooth)
+#             self._W = sp.vstack(wlist)
+#         return self._W
 
-        else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.eps_q, self.norms[1])
-            self.Rx = Utils.sdiag( self.rx )
+#     @property
+#     def Ws(self):
+#         """Regularization matrix Ws"""
+#         if getattr(self, 'm', None) is None:
+#             self.Rs = Utils.speye(self.mesh.nC)
 
-        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+#         else:
+#             f_m = self.m
+#             self.rs = self.R(f_m , self.p, self.eps)
+#             #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
+#             self.Rs = Utils.sdiag( self.rs )
 
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
+#         self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
 
-        if getattr(self, 'curModel', None) is None:
-            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+#         return self._Ws
 
-        else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
-            self.Ry = Utils.sdiag( self.ry )
+#     @property
+#     def Wx(self):
+#         """Regularization matrix Wx"""
 
-        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+#         if getattr(self, 'm', None) is None:
+#             self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
 
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
+#         else:
+#             f_m = self.mesh.unitCellGradx * self.m
+#             self.rx = self.R( f_m , self.qx, self.eps)
+#             self.Rx = Utils.sdiag( self.rx )
 
-        if getattr(self, 'curModel', None) is None:
-            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+#         if getattr(self, '_Wx', None) is None:
+#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
+#         return self._Wx
 
-        else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
-            self.Rz = Utils.sdiag( self.rz )
+#     @property
+#     def Wy(self):
+#         """Regularization matrix Wy"""
 
-        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+#         if getattr(self, 'm', None) is None:
+#             self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
 
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        #if getattr(self, '_Wsmooth', None) is None:
-        wlist = (self.Wx,)
-        if self.regmesh.dim > 1:
-            wlist += (self.Wy,)
-        if self.regmesh.dim > 2:
-            wlist += (self.Wz,)
-        #self._Wsmooth = sp.vstack(wlist)
-        return sp.vstack(wlist)
+#         else:
+#             f_m = self.mesh.unitCellGrady * self.m
+#             self.ry = self.R( f_m , self.qy, self.eps)
+#             self.Ry = Utils.sdiag( self.ry )
 
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
+#         if getattr(self, '_Wy', None) is None:
+#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
+#         return self._Wy
 
-    def R(self, f_m , eps, exponent):
+#     @property
+#     def Wz(self):
+#         """Regularization matrix Wz"""
 
-        eta = (eps**(1.-exponent/2.))**0.5
-        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
+#         if getattr(self, 'm', None) is None:
+#             self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
 
-        return r
+#         else:
+#             f_m = self.mesh.unitCellGradz * self.m
+#             self.rz = self.R( f_m , self.qz, self.eps)
+#             self.Rz = Utils.sdiag( self.rz )
+
+#         if getattr(self, '_Wz', None) is None:
+#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
+#         return self._Wz
+
+
+#     def R(self, f_m , p, dec):
+
+#         eta = (self.eps**(1-p/2.))**0.5
+#         r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
+
+#         return r
+# =======
+# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001

SimPEG/Tests.py

@@ -82,14 +82,14 @@ class OrderTest(unittest.TestCase):
     _meshType = meshTypes[0]
     meshDimension = 3
 
-    def setupMesh(self, nc):
+    def makeMesh(self, nc, meshType=_meshType, meshDimension=meshDimension):
         """
         For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
         """
-        if 'TensorMesh' in self._meshType:
-            if 'uniform' in self._meshType:
+        if 'TensorMesh' in meshType:
+            if 'uniform' in meshType:
                 h = [nc, nc, nc]
-            elif 'random' in self._meshType:
+            elif 'random' in meshType:
                 h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
@@ -97,46 +97,46 @@ class OrderTest(unittest.TestCase):
             else:
                 raise Exception('Unexpected meshType')
 
-            self.M = TensorMesh(h[:self.meshDimension])
-            max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            M = TensorMesh(h[:meshDimension])
+            max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
 
-        elif 'CylMesh' in self._meshType:
-            if 'uniform' in self._meshType:
+        elif 'CylMesh' in meshType:
+            if 'uniform' in meshType:
                 h = [nc, nc, nc]
             else:
                 raise Exception('Unexpected meshType')
 
-            if self.meshDimension == 2:
-                self.M = CylMesh([h[0], 1, h[2]])
-                max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]])
-            elif self.meshDimension == 3:
-                self.M = CylMesh(h)
-                max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            if meshDimension == 2:
+                M = CylMesh([h[0], 1, h[2]])
+                max_h = max([np.max(hi) for hi in [M.hx, M.hz]])
+            elif meshDimension == 3:
+                M = CylMesh(h)
+                max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
 
-        elif 'Curv' in self._meshType:
-            if 'uniform' in self._meshType:
+        elif 'Curv' in meshType:
+            if 'uniform' in meshType:
                 kwrd = 'rect'
-            elif 'rotate' in self._meshType:
+            elif 'rotate' in meshType:
                 kwrd = 'rotate'
             else:
                 raise Exception('Unexpected meshType')
-            if self.meshDimension == 1:
+            if meshDimension == 1:
                 raise Exception('Lom not supported for 1D')
-            elif self.meshDimension == 2:
+            elif meshDimension == 2:
                 X, Y = Utils.exampleLrmGrid([nc, nc], kwrd)
-                self.M = CurvilinearMesh([X, Y])
-            elif self.meshDimension == 3:
+                M = CurvilinearMesh([X, Y])
+            elif meshDimension == 3:
                 X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
-                self.M = CurvilinearMesh([X, Y, Z])
-            return 1./nc
+                M = CurvilinearMesh([X, Y, Z])
+            return M, 1./nc
 
-        elif 'Tree' in self._meshType:
+        elif 'Tree' in meshType:
             nc *= 2
-            if 'uniform' in self._meshType or 'notatree' in self._meshType:
+            if 'uniform' in meshType or 'notatree' in meshType:
                 h = [nc, nc, nc]
-            elif 'random' in self._meshType:
+            elif 'random' in meshType:
                 h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
@@ -145,20 +145,29 @@ class OrderTest(unittest.TestCase):
                 raise Exception('Unexpected meshType')
 
             levels = int(np.log(nc)/np.log(2))
-            self.M = Tree(h[:self.meshDimension], levels=levels)
+            M = Tree(h[:meshDimension], levels=levels)
             def function(cell):
-                if 'notatree' in self._meshType:
+                if 'notatree' in meshType:
                     return levels - 1
                 r = cell.center - np.array([0.5]*len(cell.center))
                 dist = np.sqrt(r.dot(r))
                 if dist < 0.2:
                     return levels
                 return levels - 1
-            self.M.refine(function,balance=False)
-            self.M.number(balance=False)
-            # self.M.plotGrid(showIt=True)
-            max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            M.refine(function,balance=False)
+            M.number(balance=False)
+            # M.plotGrid(showIt=True)
+            max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
+
+
+    def setupMesh(self, nc):
+        """
+        For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
+        """
+        M, h = self.makeMesh(nc, meshType=self._meshType, meshDimension=self.meshDimension)
+        self.M = M
+        return h
 
     def getError(self):
         """For given h, generate A[h], f and A(f) and return norm of error."""

SimPEG/Utils/ModelBuilder.py

@@ -88,14 +88,12 @@ def getIndicesBlock(p0,p1,ccMesh):
     # Return a tuple
     return ind
 
-def defineBlock(ccMesh,p0,p1,vals=None):
+def defineBlock(ccMesh,p0,p1,vals=[0,1]):
     """
         Build a block with the conductivity specified by condVal.  Returns an array.
         vals[0]  conductivity of the block
         vals[1]  conductivity of the ground
     """
-    if vals is None:
-        vals = [0,1]
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
     ind   = getIndicesBlock(p0,p1,ccMesh)
 
@@ -103,11 +101,7 @@ def defineBlock(ccMesh,p0,p1,vals=None):
 
     return mkvc(sigma)
 
-def defineElipse(ccMesh, center=None, anisotropy=None, slope=10., theta=0.):
-    if center is None:
-        center = [0,0,0]
-    if anisotropy is None:
-        anisotropy = [1,1,1]
+def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.):
     G = ccMesh.copy()
     dim = ccMesh.shape[1]
     for i in range(dim):
@@ -162,7 +156,7 @@ def getIndicesSphere(center,radius,ccMesh):
     # Return a tuple
     return ind
 
-def defineTwoLayers(ccMesh,depth,vals=None):
+def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     """
     Define a two layered model.  Depth of the first layer must be specified.
     CondVals vector with the conductivity values of the layers.  Eg:
@@ -173,8 +167,6 @@ def defineTwoLayers(ccMesh,depth,vals=None):
         0                          depth                                 zf
              1st layer                       2nd layer
     """
-    if vals is None:
-        vals = [0,1]
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
 
     dim = np.size(ccMesh[0,:])
@@ -260,7 +252,7 @@ def layeredModel(ccMesh, layerTops, layerValues):
 
 
 
-def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
+def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
         Create a random model by convolving a kernel with a
         uniformly distributed model.
@@ -284,8 +276,6 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
 
 
     """
-    if bounds is None:
-        bounds = [0,1]
 
     if seed is None:
         seed = np.random.randint(1e3)

SimPEG/Utils/codeutils.py

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

SimPEG/Utils/io_utils.py

@@ -1,137 +0,0 @@
-from SimPEG import np, Mesh
-import time as tm
-import vtk, vtk.util.numpy_support as npsup
-import re
-
-def read_GOCAD_ts(tsfile):
-    """
-
-    Read GOCAD triangulated surface (*.ts) file
-    INPUT:
-    tsfile: Triangulated surface
-
-    OUTPUT:
-    vrts : Array of vertices in XYZ coordinates [n x 3]
-    trgl : Array of index for triangles [m x 3]. The order of the vertices
-            is important and describes the normal
-            n = cross( (P2 - P1 ) , (P3 - P1) )
-
-    Author: @fourndo
-
-
-    .. note::
-
-        Remove all attributes from the GoCAD surface before exporting it!
-
-    """
-
-
-    fid = open(tsfile,'r')
-    line = fid.readline()
-
-    # Skip all the lines until the vertices
-    while re.match('TFACE',line)==None:
-        line = fid.readline()
-
-    line = fid.readline()
-    vrtx = []
-
-    # Run down all the vertices and save in array
-    while re.match('VRTX',line):
-        l_input  = re.split('[\s*]',line)
-        temp = np.array(l_input[2:5])
-        vrtx.append(temp.astype(np.float))
-
-        # Read next line
-        line = fid.readline()
-
-    vrtx = np.asarray(vrtx)
-
-    # Skip lines to the triangles
-    while re.match('TRGL',line)==None:
-        line = fid.readline()
-
-    # Run down the list of triangles
-    trgl = []
-
-    # Run down all the vertices and save in array
-    while re.match('TRGL',line):
-        l_input  = re.split('[\s*]',line)
-        temp = np.array(l_input[1:4])
-        trgl.append(temp.astype(np.int))
-
-        # Read next line
-        line = fid.readline()
-
-    trgl = np.asarray(trgl)
-
-    return vrtx, trgl
-
-def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True):
-    """"
-    Function to read gocad polystructure file and output indexes of mesh with in the structure.
-
-    """
-    # Adjust the index
-    trgl = trgl - 1
-
-    # Make vtk pts
-    ptsvtk = vtk.vtkPoints()
-    ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1))
-
-    # Make the polygon connection
-    polys = vtk.vtkCellArray()
-    for face in trgl:
-        poly = vtk.vtkPolygon()
-        poly.GetPointIds().SetNumberOfIds(len(face))
-        for nrv, vert in enumerate(face):
-            poly.GetPointIds().SetId(nrv,vert)
-        polys.InsertNextCell(poly)
-
-    # Make the polydata, structure of connections and vrtx
-    polyData = vtk.vtkPolyData()
-    polyData.SetPoints(ptsvtk)
-    polyData.SetPolys(polys)
-
-    # Make implicit func
-    ImpDistFunc = vtk.vtkImplicitPolyDataDistance()
-    ImpDistFunc.SetInput(polyData)
-
-    # Convert the mesh
-    vtkMesh = vtk.vtkRectilinearGrid()
-    vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz)
-    vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1))
-    vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1))
-    vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1))
-    # Add indexes
-    vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1)
-    vtkInd.SetName('Index')
-    vtkMesh.GetCellData().AddArray(vtkInd)
-
-    extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor
-    extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) #
-    extractImpDistRectGridFilt.SetInputData(vtkMesh)
-
-    if boundaries is True:
-        extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
-
-    else:
-        extractImpDistRectGridFilt.ExtractBoundaryCellsOff()
-
-    if internal is True:
-        extractImpDistRectGridFilt.ExtractInsideOn()
-
-    else:
-        extractImpDistRectGridFilt.ExtractInsideOff()
-
-    print "Extracting indices from grid..."
-    # Executing the pipe
-    extractImpDistRectGridFilt.Update()
-
-    # Get index inside
-    insideGrid = extractImpDistRectGridFilt.GetOutput()
-    insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))
-
-
-    # Return the indexes inside
-    return insideGrid

SimPEG/Utils/matutils.py

@@ -27,7 +27,7 @@ def mkvc(x, numDims=1):
 
     if isinstance(x, Zero):
         return x
-    
+
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -422,9 +422,9 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
-    @property 
+    @property
     def transpose(self): return Zero()
-    
+
     @property
     def T(self): return Zero()
 

docs/examples/DC_Forward_PseudoSection.rst

@@ -12,16 +12,8 @@
 DC Forward Simulation
 =====================
 
-Forward model two conductive spheres in a half-space and plot a
-pseudo-section. Assumes an infinite line source and measures along the
-center of the spheres.
+Forward model conductive spheres in a half-space and plot a pseudo-section
 
-INPUT:
-loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
-radi    = Radius of spheres [r1,r2]
-param   = Conductivity of background and two spheres [m0,m1,m2]
-stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
-dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
 Created by @fourndo on Mon Feb 01 19:28:06 2016
 
 

docs/examples/EM_Schenkel_Morrison_Casing.rst

@@ -1,58 +0,0 @@
-.. _examples_EM_Schenkel_Morrison_Casing:
-
-.. --------------------------------- ..
-..                                   ..
-..    THIS FILE IS AUTO GENEREATED   ..
-..                                   ..
-..    SimPEG/Examples/__init__.py    ..
-..                                   ..
-.. --------------------------------- ..
-
-
-EM: Schenkel and Morrison Casing Model
-======================================
-
-Here we create and run a FDEM forward simulation to calculate the vertical
-current inside a steel-cased. The model is based on the Schenkel and
-Morrison Casing Model, and the results are used in a 2016 SEG abstract by
-Yang et al.
-
-- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
-
-
-The model consists of:
-- Air: Conductivity 1e-8 S/m, above z = 0
-- Background: conductivity 1e-2 S/m, below z = 0
-- Casing: conductivity 1e6 S/m
-    - 300m long
-    - radius of 0.1m
-    - thickness of 6e-3m
-
-Inside the casing, we take the same conductivity as the background.
-
-We are using an EM code to simulate DC, so we use frequency low enough
-that the skin depth inside the casing is longer than the casing length (f
-= 1e-6 Hz). The plot produced is of the current inside the casing.
-
-These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
-resistivity modeling of steel casing for reservoir monitoring using
-equivalent resistor network. The solver used to produce these results and
-achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
-
-.. _pymatsolver: https://github.com/rowanc1/pymatsolver
-
-This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
-
-If you would use this example for a code comparison, or build upon it, a
-citation would be much appreciated!
-
-
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.EM_Schenkel_Morrison_Casing.run()
-
-.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
-    :language: python
-    :linenos:

docs/examples/Inversion_IRLS.rst

@@ -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:

tests/base/test_regularization.py

@@ -5,8 +5,6 @@ from scipy.sparse.linalg import dsolve
 import inspect
 
 TOL = 1e-20
-testReg = True
-testRegMesh = True
 
 class RegularizationTests(unittest.TestCase):
 
@@ -18,82 +16,44 @@ class RegularizationTests(unittest.TestCase):
         mesh3 = Mesh.TensorMesh([hx, hy, hz])
         self.meshlist = [mesh1,mesh2, mesh3]
 
-    if testReg:
-        def test_regularization(self):
-            for R in dir(Regularization):
-                r = getattr(Regularization, R)
-                if not inspect.isclass(r): continue
-                if not issubclass(r, Regularization.BaseRegularization):
-                    continue
-
-                for i, mesh in enumerate(self.meshlist):
-
-                    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) = %f' %reg.eval(reg.mref)
-                    passed = reg.eval(reg.mref) < TOL
-                    self.assertTrue(passed)
-
-                    print 'Check:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                    self.assertTrue(passed)
-
-                    print 'Check 2 Deriv:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                    self.assertTrue(passed)
-
-        def test_regularization_ActiveCells(self):
-            for R in dir(Regularization):
-                r = getattr(Regularization, R)
-                if not inspect.isclass(r): continue
-                if not issubclass(r, Regularization.BaseRegularization):
-                    continue
-
-                for i, mesh in enumerate(self.meshlist):
-
-                    print 'Testing Active Cells %iD'%(mesh.dim)
-
-                    if mesh.dim == 1:
-                        indActive = Utils.mkvc(mesh.gridCC <= 0.8)
-                    elif mesh.dim == 2:
-                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
-                    elif mesh.dim == 3:
-                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
-
-                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
-
-                    for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
-                        reg = r(mesh, mapping=mapping, indActive=indAct)
-                        m = np.random.rand(mesh.nC)[indAct]
-                        reg.mref = np.ones_like(m)*np.mean(m)
-
-                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                    passed = reg.eval(reg.mref) < TOL
-                    self.assertTrue(passed)
-
-                    print 'Check:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                    self.assertTrue(passed)
-
-                    print 'Check 2 Deriv:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                    self.assertTrue(passed)
-
-    if testRegMesh:
-        def test_regularizationMesh(self):
+    def test_regularization(self):
+        for R in dir(Regularization):
+            r = getattr(Regularization, R)
+            if not inspect.isclass(r): continue
+            if not issubclass(r, Regularization.BaseRegularization):
+                continue
 
             for i, mesh in enumerate(self.meshlist):
 
                 print 'Testing %iD'%mesh.dim
 
-                # mapping = r.mapPair(mesh)
-                # reg = r(mesh, mapping=mapping)
-                # m = np.random.rand(mapping.nP)
+                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) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+    def test_regularization_ActiveCells(self):
+        for R in dir(Regularization):
+            r = getattr(Regularization, R)
+            if not inspect.isclass(r): continue
+            if not issubclass(r, Regularization.BaseRegularization):
+                continue
+
+            for i, mesh in enumerate(self.meshlist):
+
+                print 'Testing Active Cells %iD'%(mesh.dim)
 
                 if mesh.dim == 1:
                     indAct = Utils.mkvc(mesh.gridCC <= 0.8)
@@ -102,9 +62,23 @@ class RegularizationTests(unittest.TestCase):
                 elif mesh.dim == 3:
                     indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct)
+                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
 
-                assert (regmesh.vol == mesh.vol[indAct]).all()
+                reg = r(mesh, mapping=mapping, indActive=indAct)
+                m = np.random.rand(mesh.nC)[indAct]
+                reg.mref = np.ones_like(m)*np.mean(m)
+
+                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
 
 
 if __name__ == '__main__':

tests/em/fdem/inverse/derivs/test_FDEM_derivs.py

@@ -8,8 +8,8 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem
 
 testE = True
 testB = True
-testH = True
-testJ = True
+testH = False
+testJ = False
 
 verbose = False
 
@@ -20,8 +20,8 @@ MU = mu_0
 freq = 1e-1
 addrandoms = True
 
-SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
-
+# SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+SrcType = ['PrimSecCyl']
 
 def derivTest(fdemType, comp):
 

tests/em/tdem/test_TDEM_forward_Analytic.py

@@ -10,9 +10,7 @@ except ImportError, e:
     MumpsSolver = SolverLU
 
 
-def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None, showIt=False):
-    if bounds is None:
-        bounds = [1e-5,1e-3]
+def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
     if meshType == 'CYL':
         cs, ncx, ncz, npad = 5., 30, 10, 15
         hx = [(cs,ncx), (cs,npad,1.3)]

tests/mesh/test_InterpolationMesh2Mesh.py

@@ -0,0 +1,353 @@
+import numpy as np
+import unittest
+from SimPEG.Utils import mkvc
+from SimPEG import Mesh, Tests
+import unittest
+
+test1D = True
+test2D = True
+test3D = False
+
+call1 = lambda fun, xyz: fun(xyz)
+call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, -1])
+call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
+cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
+cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
+cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
+cartF2Cyl = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFz, fx, fy)))
+cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
+cartE2Cyl = lambda M, ex, ey: cart_row2(M.gridEy, ex, ey)
+cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
+cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
+
+TOL = 1e-7
+
+if test1D:
+    class TestInterpolationMesh2Mesh_Tensor1D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor1D'
+        meshSizes = [8, 16, 32]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 1
+
+        def getError(self):
+            funX = lambda x: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+            ana = call1(funX, getattr(mesh2, 'grid%s'%self.type))
+
+            v = call1(funX, getattr(self.M, 'grid%s'%self.type))
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_1D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor1D: CC'
+            self.orderTest()
+
+        def test_orderN_1D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor1D: N'
+            self.orderTest()
+
+        def test_orderEx_1D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor1D: Ex'
+            self.orderTest()
+
+        def test_orderFx_1D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor1D: Fx'
+            self.orderTest()
+
+if test2D:
+    class TestInterpolationMesh2Mesh_Tensor2D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor2D'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 2
+
+        def getError(self):
+            funX = lambda x, y: np.cos(2*np.pi*y)
+            funY = lambda x, y: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF2(mesh2, funX, funY)
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE2(mesh2, funX, funY)
+                ana = mesh2.projectEdgeVector(ana)
+            else:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF2(self.M, funX, funY)
+                if 'x' in self.type or 'y' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE2(self.M, funX, funY)
+                if 'x' in self.type or 'y' in self.type:
+                    v = self.M.projectEdgeVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'CC' == self.type:
+                v = call2(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call2(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            # print P.shape, v.shape
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_2D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor2D: CC'
+            self.orderTest()
+
+        def test_orderN_2D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor2D: N'
+            self.orderTest()
+
+        def test_orderE_2D(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor2D: E'
+            self.orderTest()
+
+        def test_orderEx_2D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor2D: Ex'
+            self.orderTest()
+
+        def test_orderEy_2D(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor2D: Ey'
+            self.orderTest()
+
+        def test_orderF_2D(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor2D: F'
+            self.orderTest()
+
+        def test_orderFx_2D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor2D: Fx'
+            self.orderTest()
+
+        def test_orderFy_2D(self):
+            self.type = 'Fy'
+            self.name = 'Mesh2Mesh Tensor2D: Fy'
+            self.orderTest()
+
+    class TestInterpolationMesh2Mesh_Cyl(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Cyl'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformCylMesh']
+        meshDimension = 2
+
+        def getError(self):
+            funX = lambda x, y: np.cos(2*np.pi*y)
+            funY = lambda x, y: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'z' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF2Cyl(mesh2, funX, funY)
+                ana = np.c_[ana[:,0], np.zeros_like(ana[:,0]), ana[:,1]]
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE2Cyl(mesh2, funX, funY)
+                ana = np.c_[np.zeros_like(ana[:,1]),ana[:,1],np.zeros_like(ana[:,1])]
+                ana = mesh2.projectEdgeVector(ana)
+            else:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF2Cyl(self.M, funX, funY)
+                v = np.c_[v[:,0], np.zeros_like(v[:,0]),v[:,1]]
+                if 'x' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = np.c_[v[:,0], v[:,2]]
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE2Cyl(self.M, funX, funY)
+                v = np.c_[np.zeros_like(v[:,1]), v[:,1],np.zeros_like(v[:,1])]
+                v = self.M.projectEdgeVector(v)
+
+            elif 'CC' == self.type:
+                v = call2(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call2(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_Cyl(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor2D: CC'
+            self.orderTest()
+
+        def test_orderN_Cyl(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor2D: N'
+            self.orderTest()
+
+        def test_orderE_Cyl(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor2D: E'
+            self.orderTest()
+
+        def test_orderEy_Cyl(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor2D: Ey'
+            self.orderTest()
+
+        def test_orderF_Cyl(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor2D: F'
+            self.orderTest()
+
+        def test_orderFx_Cyl(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor2D: Fx'
+            self.orderTest()
+
+        def test_orderFz_Cyl(self):
+            self.type = 'Fz'
+            self.name = 'Mesh2Mesh Tensor2D: Fz'
+            self.orderTest()
+
+if test3D:
+    class TestInterpolationMesh2Mesh_Tensor3D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor3D'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 3
+
+        def getError(self):
+            funX = lambda x, y, z: np.cos(2*np.pi*y)
+            funY = lambda x, y, z: np.cos(2*np.pi*z)
+            funZ = lambda x, y, z: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call3(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call3(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'z' in self.type:
+                ana = call3(funZ, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF3(mesh2, funX, funY, funZ)
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE3(mesh2, funX, funY, funZ)
+                ana = mesh2.projectFaceVector(ana)
+            else:
+                ana = call3(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF3(self.M, funX, funY, funZ)
+                if 'x' in self.type or 'y' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE3(self.M, funX, funY, funZ)
+                if 'x' in self.type or 'y' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'CC' == self.type:
+                v = call3(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call3(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            # print P.shape, v.shape
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_3D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor3D: CC'
+            self.orderTest()
+
+        def test_orderN_3D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor3D: N'
+            self.orderTest()
+
+        def test_orderE_3D(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor3D: E'
+            self.orderTest()
+
+        def test_orderEx_3D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor3D: Ex'
+            self.orderTest()
+
+        def test_orderEy_3D(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor3D: Ey'
+            self.orderTest()
+
+        def test_orderEz_3D(self):
+            self.type = 'Ez'
+            self.name = 'Mesh2Mesh Tensor3D: Ez'
+            self.orderTest()
+
+        def test_orderF_3D(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor3D: F'
+            self.orderTest()
+
+        def test_orderFx_3D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor3D: Fx'
+            self.orderTest()
+
+        def test_orderFy_3D(self):
+            self.type = 'Fy'
+            self.name = 'Mesh2Mesh Tensor3D: Fy'
+            self.orderTest()
+
+        def test_orderFz_3D(self):
+            self.type = 'Fz'
+            self.name = 'Mesh2Mesh Tensor3D: Fz'
+            self.orderTest()
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/test_Problem3D_againstAnalytic.py

@@ -242,6 +242,9 @@ class TestAnalytics(unittest.TestCase):
     def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
     def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
 
+    # Do a derivative test
+    def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
+
     # Do a derivative test of Jvec
     # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
     # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))