Minor changes...test runs fine from here

Fixed the boolean error in Directives.

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):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
     from scipy.interpolate import griddata
     import pylab as plt
 
-    # Set depth to 0 for now
-    z0 = 0.
-
     # Pre-allocate
     midx = []
     midz = []
@@ -259,38 +256,53 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
         midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
 
-    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, clim=(vmin, vmax))
-    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
+    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")
 
-    cmin,cmax = cbar.get_clim()
-    ticks = np.linspace(cmin,cmax,3)
-    cbar.set_ticks(ticks)
-    cbar.ax.tick_params(labelsize=10)
+        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)
 
-    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)
 
-    # Plot apparent resistivity
-    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
 
-    #ax.set_xticklabels([])
-    #ax.set_yticklabels([])
+    if not axlabel:
+        axs.set_xticklabels([])
+        axs.set_yticklabels([])
 
     plt.gca().set_aspect('equal', adjustable='box')
 
@@ -448,15 +460,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
     """
         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
@@ -471,10 +483,16 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
     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')
-    fid.write('! ' + stype + ' FORMAT\n')
-
+    
+    
+    if iptype!=0:
+        fid.write('IPTYPE=%i\n'%iptype)
+        
+    else:
+        fid.write('! ' + stype + ' FORMAT\n')
+        
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -498,7 +516,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
                 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')
 
 
@@ -506,18 +524,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
                 if stype == 'SURFACE':
 
-                    fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("%f " % 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='%e',delimiter=' ',newline='\n')
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
 
         if dtype=='3D':
 
@@ -529,11 +554,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
             if stype == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx))
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
 
             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()
@@ -640,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     DCsurvey2D.std = np.asarray(DCsurvey.std)
 
     return DCsurvey2D
-
-def readUBC_DC3Dobs(fileName):
+    
+def readUBC_DC3Dobs(fileName, dtype = 'DC'):
     """
-        Read UBC GIF DCIP 3D observation file and generate survey
+        Read UBC GIF IP 3D observation file and generate survey
 
         Input:
         :param fileName, path to the UBC GIF 3D obs file
 
         Output:
-        :param DCIPsurvey
+        :param IPsurvey
         :return
 
-        Created on Mon April 6th, 2015
-
         @author: dominiquef
 
     """
-
+    zflag = True # Flag for z value provided
+    
     # Load file
-    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-
+    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'"
+        
     # 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 is transmitter with number of receivers
+        # First line or end of a transmitter block, read transmitter info
         if count==0:
-
-            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+            # Read the line
+            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[0:2],np.nan]
-                zflag = False
+                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
+                
+                zflag = False # Pass on the flag to the receiver loc
 
             else:
                 tx = temp[:-1]
@@ -692,8 +726,16 @@ def readUBC_DC3Dobs(fileName):
             rx = []
             continue
 
-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
 
+        # 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])
@@ -703,7 +745,7 @@ def readUBC_DC3Dobs(fileName):
                 wd.append(temp[-1])
 
         else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
             # Check if there is data with the location
             if len(temp)==6:
                 d.append(temp[-2])
@@ -711,7 +753,7 @@ def readUBC_DC3Dobs(fileName):
 
         count = count -1
 
-        # Reach the end of transmitter block
+        # Reach the end of transmitter block, append the src, rx and continue
         if count == 0:
             rx = np.asarray(rx)
             Rx = DC.RxDipole(rx[:,:3],rx[:,3:])

SimPEG/Directives.py

@@ -222,7 +222,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
             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,41 +237,6 @@ 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.mrefInSmooth == 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):
 
@@ -293,6 +258,7 @@ class Update_IRLS(InversionDirective):
     phi_m_last = None
     phi_d_last = None
 
+
     def initialize(self):
 
         # Scale the regularization for changes in norm
@@ -305,12 +271,13 @@ class Update_IRLS(InversionDirective):
 
         self.reg.curModel = self.invProb.curModel
         self.reg.gamma = self.gamma
-
+        
+        self.reg._W = None
         if getattr(self, 'phi_d_last', None) is None:
             self.phi_d_last = self.invProb.phi_d
 
     def endIter(self):
-        # Cool the threshold parameter
+        # Cool the threshold parameter if required
         if getattr(self, 'factor', None) is not None:
             eps = self.reg.eps / self.factor
 
@@ -321,32 +288,49 @@ class Update_IRLS(InversionDirective):
 
         # Get phi_m at the end of current iteration
         self.phi_m_last = self.invProb.phi_m_last
-
+        
+        self.reg._W = None
          # Update the model used for the IRLS weights
         self.reg.curModel = self.invProb.curModel
 
-        # Temporarely set gamma to 1.
+        # Temporarely set gamma to 1. to get raw phi_m
         self.reg.gamma = 1.
 
-        # Compute change in model objective function and update scaling
+        # 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
-
-        self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+        
+        # 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):
-
-
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
+    
+    def initialize(self):
+    
+        if getattr(self.opt, 'approxHinv', None) is 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((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+            
     def endIter(self):
         # Cool the threshold parameter
-
+        if self.onlyOnStart==True:
+            return
+            
         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.)
+            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((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
-            print 'Updated pre-cond'
+
 
 class Update_Wj(InversionDirective):
     """
@@ -373,3 +357,19 @@ class Update_Wj(InversionDirective):
             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/Analytics/DC.py

@@ -1,119 +0,0 @@
-import numpy as np
-from scipy.constants import mu_0, pi
-from scipy import special
-
-def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
-    """
-        Analytic solution for electric potential from a postive pole
-
-        Input variables:
-
-            txloc =  a xyz location of A (+) electrode (np.r_[xa, ya, za])
-
-            rxlocs = [M, N]
-                M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
-                N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
-
-            sigma = conductivity (either float or complex)
-            flag = "wholsespace" or "halfspace"
-
-    """
-    M = rxlocs[0]
-    N = rxlocs[1]
-
-    rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
-    rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
-
-    phiM = 1./(4*np.pi*rM*sigma)
-    phiN = 1./(4*np.pi*rN*sigma)
-    phi = phiM - phiN
-
-    if flag == "halfspace":
-        phi *= 2
-
-    return phi
-
-deg2rad  = lambda deg: deg/180.*np.pi
-rad2deg  = lambda rad: rad*180./np.pi
-
-def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
-                 flag = "sec", order=12, halfspace=False):
-# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
-#                  flag = "sec", order=12):
-    """
-
-        Parameters:
-
-            txloc (array) : current electrode location (x,y,z)
-            xc (float)    : x center of depressed sphere
-            rxloc (array) : electrode locations
-                              (Nx3 array, # of electrodes)
-            radius (float): radius of the sphere (m)
-            rho (float)   : resistivity of the background (ohm-m)
-            rho1 (float)  : resistivity of the sphere
-            flag (string) : "sec", "total", "prim"
-                              (default="sec")
-                              "sec": secondary potential only due to sphere
-                              "prim": primary potential from the point source
-                              "total": "sec"+"prim"
-            order (float) : maximum order of Legendre polynomial
-                              (default=12)
-
-        Written by Seogi Kang (skang@eos.ubc.ca)
-        Ph.D. Candidate of University of British Columbia, Canada
-
-    """
-
-    Pleg = []
-    # Compute Legendre Polynomial
-    for i in range(order):
-        Pleg.append(special.legendre(i, monic=0))
-
-
-    rho = 1./sigma
-    rho1 = 1./sigma1
-
-    # Center of the sphere should be aligned in txloc in y-direction
-    yc = txloc[1]
-    xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]]
-    r = np.sqrt( (xyz**2).sum(axis=1) )
-
-    x0 = abs(txloc[0]-xc)
-
-    costheta = xyz[:,0]/r * (txloc[0]-xc)/x0
-    phi = np.zeros_like(r)
-    R = (r**2+x0**2.-2.*r*x0*costheta)**0.5
-    # primary potential in a whole space
-    prim = rho*1./(4*np.pi*R)
-
-    if flag =="prim":
-        return prim
-
-    sphind = r < radius
-    out = np.zeros_like(r)
-    for n in range(order):
-        An, Bn = AnBnfun(n, radius, x0, rho, rho1)
-        dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind])
-        out[~sphind] += dumout
-        dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind])
-        out[sphind] += dumin
-
-    out[~sphind] += prim[~sphind]
-
-    if halfspace:
-        scale = 2
-    else:
-        scale = 1
-
-    if flag == "sec":
-        return scale*(out-prim)
-    elif flag == "total":
-        return scale*out
-
-def AnBnfun(n, radius, x0, rho, rho1, I=1.):
-    const = I*rho/(4*np.pi)
-    bunmo = n*rho + (n+1)*rho1
-    An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \
-        (rho1-rho) / bunmo
-    Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo
-    return An, Bn

SimPEG/EM/Analytics/__init__.py

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

SimPEG/EM/Base.py

@@ -1,7 +1,6 @@
 from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 
-
 class EMPropMap(Maps.PropMap):
     """
         Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
@@ -71,12 +70,6 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Mf = self.mesh.getFaceInnerProduct()
         return self._Mf
 
-    @property
-    def Vol(self):
-        if getattr(self, '_Vol', None) is None:
-            self._Vol = Utils.sdiag(self.mesh.vol)
-        return self._Vol
-
 
     # ----- Magnetic Permeability ----- #
     @property
@@ -134,6 +127,7 @@ class BaseEMProblem(Problem.BaseProblem):
         """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
 
+
     @property
     def MeSigmaI(self):
         """
@@ -156,6 +150,7 @@ class BaseEMProblem(Problem.BaseProblem):
         return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
         # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
 
+
     @property
     def MfRho(self):
         """
@@ -188,12 +183,7 @@ class BaseEMProblem(Problem.BaseProblem):
         """
             Derivative of :code:`MfRhoI` with respect to the model.
         """
-
-        dMfRhoI_dI = -self.MfRhoI**2
-        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
-        return dMfRhoI_dI * ( dMf_drho * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) )
-
-        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
 
 class BaseEMSurvey(Survey.BaseSurvey):
 
@@ -202,7 +192,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
         self.srcList = srcList
         Survey.BaseSurvey.__init__(self, **kwargs)
 
-    def eval(self, f):
+    def eval(self, u):
         """
         Project fields to receiver locations
         :param Fields u: fields object
@@ -212,8 +202,8 @@ class BaseEMSurvey(Survey.BaseSurvey):
         data = Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, f)
+                data[src, rx] = rx.eval(src, self.mesh, u)
         return data
 
-    def evalDeriv(self, f):
+    def evalDeriv(self, u):
         raise Exception('Use Receivers to project fields deriv.')

SimPEG/EM/FDEM/FDEM.py

@@ -167,7 +167,6 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
-            #Why are you adding?
             s_m[:,i] = s_m[:,i] + smi
             s_e[:,i] = s_e[:,i] + sei
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -181,7 +181,7 @@ class Fields_e(Fields):
                   }
 
     def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
+        Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob

SimPEG/EM/FDEM/SrcFDEM.py

@@ -14,6 +14,7 @@ class BaseSrc(Survey.BaseSrc):
 
     def eval(self, prob):
         """
+        Evaluate the source terms.
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -63,9 +63,9 @@ class Rx(SimPEG.Survey.BaseRx):
         """Component projection (real/imag)"""
         return self.knownRxTypes[self.rxType][2]
 
-    def projGLoc(self, f):
+    def projGLoc(self, u):
         """Grid Location projection (e.g. Ex Fy ...)"""
-        return f._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
     def eval(self, src, mesh, f):
         """

SimPEG/EM/Static/DC/BoundaryUtils.py

@@ -1,160 +0,0 @@
-import numpy as np
-
-def getxBCyBC_CC(mesh, alpha, beta, gamma):
-# def getxBCyBC(mesh, alpha, beta, gamma):
-    """
-    This is a subfunction generating mixed-boundary condition:
-
-    .. math::
-
-        \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
-
-        \rho \vec{j} = -\nabla \phi \phi
-
-        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
-
-        xBC = f_1(\alpha, \beta, \gamma)
-        yBC = f(\alpha, \beta, \gamma)
-
-    Computes xBC and yBC for cell-centered discretizations
-    """
-    if mesh.dim == 1: #1D
-        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
-            raise Exception("Lenght of list, alpha should be 2")
-        fCCxm,fCCxp = mesh.cellBoundaryInd
-        nBC = fCCxm.sum()+fCCxp.sum()
-        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
-        h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-
-        xBC = np.r_[xBC_xm, xBC_xp]
-        yBC = np.r_[yBC_xm, yBC_xp]
-
-    elif mesh.dim == 2: #2D
-        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
-            raise Exception("Lenght of list, alpha should be 4")
-
-        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
-        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
-        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
-        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
-
-        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
-        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
-        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
-        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-        xBC_ym = 0.5*a_ym
-        xBC_yp = 0.5*a_yp/b_yp
-        yBC_ym = 0.5*(1.-b_ym)
-        yBC_yp = 0.5*(1.-1./b_yp)
-
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
-
-        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
-        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
-        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
-        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
-
-        xBC = np.r_[xBC_x, xBC_y]
-        yBC = np.r_[yBC_x, yBC_y]
-
-    elif mesh.dim == 3: #3D
-        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
-            raise Exception("Lenght of list, alpha should be 6")
-        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
-        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
-        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
-        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
-        alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4]
-        alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
-        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
-        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
-
-        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
-        h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
-        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
-        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
-        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
-
-        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
-        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
-        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
-        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-        xBC_ym = 0.5*a_ym
-        xBC_yp = 0.5*a_yp/b_yp
-        yBC_ym = 0.5*(1.-b_ym)
-        yBC_yp = 0.5*(1.-1./b_yp)
-        xBC_zm = 0.5*a_zm
-        xBC_zp = 0.5*a_zp/b_zp
-        yBC_zm = 0.5*(1.-b_zm)
-        yBC_zp = 0.5*(1.-1./b_zp)
-
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
-        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
-
-        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
-        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
-        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
-
-        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
-        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
-        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
-
-        xBC = np.r_[xBC_x, xBC_y, xBC_z]
-        yBC = np.r_[yBC_x, yBC_y, yBC_z]
-
-    return xBC, yBC

SimPEG/EM/Static/DC/FieldsDC.py

@@ -1,148 +0,0 @@
-import SimPEG
-from SimPEG.Utils import Identity, Zero
-import numpy as np
-from scipy.constants import epsilon_0
-
-class Fields(SimPEG.Problem.Fields):
-    knownFields = {}
-    dtype = float
-
-    def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
-
-        return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float)
-
-    def _eDeriv(self, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
-        return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = float)
-
-    def _jDeriv(self, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
-        return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = float)
-
-
-class Fields_CC(Fields):
-    knownFields = {'phiSolution':'CC'}
-    aliasFields = {
-                    'phi': ['phiSolution','CC','_phi'],
-                    'j' : ['phiSolution','F','_j'],
-                    'e' : ['phiSolution','F','_e'],
-                    'charge' : ['phiSolution','CC','_charge'],
-                  }
-                  # primary - secondary
-                  # CC variables
-
-    def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
-        mesh.setCellGradBC("neumann")
-        cellGrad = mesh.cellGrad
-    def startup(self):
-        self.prob = self.survey.prob
-
-    def _GLoc(self, fieldType):
-        if fieldType == 'phi':
-            return 'CC'
-        elif fieldType == 'e' or fieldType == 'j':
-            return 'F'
-        else:
-            raise Exception('Field type must be phi, e, j')
-
-    def _phi(self, phiSolution, srcList):
-        return phiSolution
-
-    def _phiDeriv_u(self, src, v, adjoint = False):
-        return Identity()*v
-
-    def _phiDeriv_m(self, src, v, adjoint = False):
-        return Zero()
-
-    def _j(self, phiSolution, srcList):
-        """
-            .. math::
-                \mathbf{j} = \mathbf{M}^{f \ -1}_{\rho} \mathbf{G} \phi
-        """
-        return self.prob.MfRhoI*self.prob.Grad*phiSolution
-
-    def _e(self, phiSolution, srcList):
-        """
-            In HJ formulation e is not well-defined!!
-            .. math::
-                \vec{e} = -\nabla \phi
-        """
-        return -self.mesh.cellGrad*phiSolution
-
-    def _charge(self, phiSolution, srcList):
-        """
-            .. math::
-                \int \nabla \codt \vec{e} =  \int \frac{\rho_v }{\epsillon_0}
-        """
-        return epsilon_0*self.prob.Vol*(self.mesh.faceDiv*self._e(phiSolution, srcList))
-
-class Fields_N(Fields):
-    knownFields = {'phiSolution':'N'}
-    aliasFields = {
-                    'phi': ['phiSolution','N','_phi'],
-                    'j' : ['phiSolution','E','_j'],
-                    'e' : ['phiSolution','E','_e'],
-                    'charge' : ['phiSolution','N','_charge'],
-                  }
-                  # primary - secondary
-                  # N variables
-
-    def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
-
-    def startup(self):
-        self.prob = self.survey.prob
-
-    def _GLoc(self, fieldType):
-        if fieldType == 'phi':
-            return 'N'
-        elif fieldType == 'e' or fieldType == 'j':
-            return 'E'
-        else:
-            raise Exception('Field type must be phi, e, j')
-
-    def _phi(self, phiSolution, srcList):
-        return phiSolution
-
-    def _phiDeriv_u(self, src, v, adjoint = False):
-        return Identity()*v
-
-    def _phiDeriv_m(self, src, v, adjoint = False):
-        return Zero()
-
-    def _j(self, phiSolution, srcList):
-        """
-            In EB formulation j is not well-defined!!
-            .. math::
-                \mathbf{j} = - \mathbf{M}^{e}_{\sigma} \mathbf{G} \phi
-        """
-        return self.prob.MeSigma * self._e(phiSolution, srcList)
-
-    def _e(self, phiSolution, srcList):
-        """
-            In HJ formulation e is not well-defined!!
-            .. math::
-                \vec{e} = -\nabla \phi
-        """
-        return -self.mesh.nodalGrad * phiSolution
-
-    def _charge(self, phiSolution, srcList):
-        """
-            .. math::
-                \int \nabla \codt \vec{e} =  \int \frac{\rho_v }{\epsillon_0}
-        """
-        return - epsilon_0*(self.mesh.nodalGrad.T*self.mesh.getEdgeInnerProduct()*self._e(phiSolution, srcList))

SimPEG/EM/Static/DC/FieldsDC_2D.py

@@ -1,146 +0,0 @@
-import SimPEG
-from SimPEG.Utils import Identity, Zero
-import numpy as np
-
-class Fields_ky(SimPEG.Problem.TimeFields):
-
-    """
-
-    Fancy Field Storage for a 2.5D code.
-
-    u[:,'phi', kyInd] = phi
-    print u[src0,'phi']
-
-    Only one field type is stored for
-    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
-    .. code-block:: python
-        f = problem.fields(m)
-        e = f[srcList,'e']
-        j = f[srcList,'j']
-
-    If accessing all sources for a given field, use the :code:`:`
-    .. code-block:: python
-        f = problem.fields(m)
-        phi = f[:,'phi']
-        e = f[:,'e']
-        b = f[:,'b']
-    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
-    """
-
-    knownFields = {}
-    dtype = float
-
-    def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
-            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
-
-        return np.array(self._phiDeriv_u(kyInd, src, du_dm_v, adjoint) + self._phiDeriv_m(kyInd, src, v, adjoint), dtype = float)
-
-    def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
-            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
-        return np.array(self._eDeriv_u(kyInd, src, du_dm_v, adjoint) + self._eDeriv_m(kyInd, src, v, adjoint), dtype = float)
-
-    def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
-        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
-            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
-
-        if adjoint:
-            return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
-        return np.array(self._jDeriv_u(kyInd, src, du_dm_v, adjoint) + self._jDeriv_m(kyInd, src, v, adjoint), dtype = float)
-
-
-    # def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-    #     if adjoint is True:
-    #         return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
-    #     return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
-
-    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-    #     if adjoint is True:
-    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
-    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
-
-
-class Fields_ky_CC(Fields_ky):
-    knownFields = {'phiSolution':'CC'}
-    aliasFields = {
-                    'phi': ['phiSolution','CC','_phi'],
-                    'j' : ['phiSolution','F','_j'],
-                    'e' : ['phiSolution','F','_e'],
-                  }
-                  # primary - secondary
-                  # CC variables
-
-    def __init__(self, mesh, survey, **kwargs):
-        Fields_ky.__init__(self, mesh, survey, **kwargs)
-
-    def startup(self):
-        self.prob = self.survey.prob
-
-    def _GLoc(self, fieldType):
-        if fieldType == 'phi':
-            return 'CC'
-        elif fieldType == 'e' or fieldType == 'j':
-            return 'F'
-        else:
-            raise Exception('Field type must be phi, e, j')
-
-    def _phi(self, phiSolution, src, kyInd):
-        return phiSolution
-
-    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
-        return Identity()*v
-
-    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
-        return Zero()
-
-    def _j(self, phiSolution, srcList):
-        raise NotImplementedError
-
-    def _e(self, phiSolution, srcList):
-        raise NotImplementedError
-
-class Fields_ky_N(Fields_ky):
-    knownFields = {'phiSolution':'N'}
-    aliasFields = {
-                    'phi': ['phiSolution','N','_phi'],
-                    'j' : ['phiSolution','E','_j'],
-                    'e' : ['phiSolution','E','_e'],
-                  }
-                  # primary - secondary
-                  # CC variables
-
-    def __init__(self, mesh, survey, **kwargs):
-        Fields_ky.__init__(self, mesh, survey, **kwargs)
-
-    def startup(self):
-        self.prob = self.survey.prob
-
-    def _GLoc(self, fieldType):
-        if fieldType == 'phi':
-            return 'N'
-        elif fieldType == 'e' or fieldType == 'j':
-            return 'E'
-        else:
-            raise Exception('Field type must be phi, e, j')
-
-    def _phi(self, phiSolution, src, kyInd):
-        return phiSolution
-
-    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
-        return Identity()*v
-
-    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
-        return Zero()
-
-    def _j(self, phiSolution, srcList):
-        raise NotImplementedError
-
-    def _e(self, phiSolution, srcList):
-        raise NotImplementedError

SimPEG/EM/Static/DC/ProblemDC.py

@@ -1,307 +0,0 @@
-from SimPEG import Problem, Utils
-from SimPEG.EM.Base import BaseEMProblem
-from SurveyDC import Survey
-from FieldsDC import Fields, Fields_CC, Fields_N
-from SimPEG.Utils import sdiag
-import numpy as np
-from SimPEG.Utils import Zero
-from BoundaryUtils import getxBCyBC_CC
-
-class BaseDCProblem(BaseEMProblem):
-
-    surveyPair = Survey
-    fieldsPair = Fields
-    Ainv = None
-
-    def fields(self, m):
-        self.curModel = m
-
-        if not self.Ainv == None:
-            self.Ainv.clean()
-
-        f = self.fieldsPair(self.mesh, self.survey)
-        A = self.getA()
-        self.Ainv = self.Solver(A, **self.solverOpts)
-        RHS = self.getRHS()
-        u = self.Ainv * RHS
-        Srcs = self.survey.srcList
-        f[Srcs, self._solutionType] = u
-        return f
-
-    def Jvec(self, m, v, f=None):
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        Jv = self.dataPair(self.survey) #same size as the data
-
-        A = self.getA()
-
-        for src in self.survey.srcList:
-            u_src = f[src, self._solutionType] # solution vector
-            dA_dm_v = self.getADeriv(u_src, v)
-            dRHS_dm_v = self.getRHSDeriv(src, v)
-            du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
-
-            for rx in src.rxList:
-                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
-        return Utils.mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-        AT = self.getA()
-
-
-        for src in self.survey.srcList:
-            u_src = f[src, self._solutionType]
-            for rx in src.rxList:
-                PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
-
-                ATinvdf_duT = self.Ainv * df_duT
-
-                dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
-                dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
-                du_dmT = -dA_dmT + dRHS_dmT
-                Jtv += (df_dmT + du_dmT).astype(float)
-
-        return Utils.mkvc(Jtv)
-
-    def getSourceTerm(self):
-        """
-        takes concept of source and turns it into a matrix
-        """
-        """
-        Evaluates the sources, and puts them in matrix form
-
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: q (nC or nN, nSrc)
-        """
-
-        Srcs = self.survey.srcList
-
-        if self._formulation is 'EB':
-            n = self.mesh.nN
-            # return NotImplementedError
-
-        elif self._formulation is 'HJ':
-            n = self.mesh.nC
-
-        q = np.zeros((n, len(Srcs)))
-
-        for i, src in enumerate(Srcs):
-            q[:,i] = src.eval(self)
-        return q
-
-class Problem3D_CC(BaseDCProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'HJ' # CC potentials means J is on faces
-    fieldsPair    = Fields_CC
-
-    def __init__(self, mesh, **kwargs):
-        BaseDCProblem.__init__(self, mesh, **kwargs)
-        self.setBC()
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        D = self.Div
-        G = self.Grad
-        # TODO: this won't work for full anisotropy
-        MfRhoI = self.MfRhoI
-        A = D * MfRhoI * G
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * A
-        return A
-
-    def getADeriv(self, u, v, adjoint= False):
-
-        D = self.Div
-        G = self.Grad
-        MfRhoIDeriv = self.MfRhoIDeriv
-
-        if adjoint:
-            # if self._makeASymmetric is True:
-            #     v = V * v
-            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
-        return D * (MfRhoIDeriv( G * u ) * v)
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return self.Vol.T * RHS
-
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-    def setBC(self):
-        if self.mesh.dim==3:
-            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-            gBFzm = self.mesh.gridFz[fzm,:]
-            gBFzp = self.mesh.gridFz[fzp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-            beta_zm, beta_zp = temp_zm, temp_zp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
-
-        elif self.mesh.dim==2:
-
-            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
-        V = self.Vol
-        self.Div = V * self.mesh.faceDiv
-        P_BC, B = self.mesh.getBCProjWF_simple()
-        M = B*self.mesh.aveCC2F
-        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
-
-
-class Problem3D_N(BaseDCProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'EB' # N potentials means B is on faces
-    fieldsPair    = Fields_N
-
-    def __init__(self, mesh, **kwargs):
-        BaseDCProblem.__init__(self, mesh, **kwargs)
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        # TODO: this won't work for full anisotropy
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        A = Grad.T * MeSigma * Grad
-
-        # Handling Null space of A
-        A[0,0] = A[0,0] + 1.
-
-        return A
-
-    def getADeriv(self, u, v, adjoint=False):
-        """
-
-        Product of the derivative of our system matrix with respect to the model and a vector
-
-        """
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        if not adjoint:
-            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
-        elif adjoint:
-            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
-
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-
-
-

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -1,351 +0,0 @@
-from SimPEG import Problem, Utils
-from SimPEG.EM.Base import BaseEMProblem
-from SurveyDC import Survey, Survey_ky
-from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
-from SimPEG.Utils import sdiag
-import numpy as np
-from SimPEG.Utils import Zero
-from BoundaryUtils import getxBCyBC_CC
-
-class BaseDCProblem_2D(BaseEMProblem):
-
-    surveyPair = Survey_ky
-    fieldsPair = Fields_ky
-    nky = 15
-    kys = np.logspace(-4, 1, nky)
-    Ainv = [None for i in range(nky)]
-    nT = nky # Only for using TimeFields
-
-    def fields(self, m):
-        self.curModel = m
-
-        if not self.Ainv[0] == None:
-            for i in range(self.nky):
-                self.Ainv[i].clean()
-
-        f = self.fieldsPair(self.mesh, self.survey)
-        Srcs = self.survey.srcList
-        for iky in range(self.nky):
-            ky = self.kys[iky]
-            A = self.getA(ky)
-            self.Ainv[iky] = self.Solver(A, **self.solverOpts)
-            RHS = self.getRHS(ky)
-            u = self.Ainv[iky] * RHS
-            f[Srcs, self._solutionType, iky] = u
-        return f
-
-    def Jvec(self, m, v, f=None):
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        Jv = self.dataPair(self.survey) #same size as the data
-        Jv0 = self.dataPair(self.survey)
-
-        # Assume y=0.
-        # This needs some thoughts to implement in general when src is dipole
-        dky = np.diff(self.kys)
-        dky = np.r_[dky[0], dky]
-        y = 0.
-
-        #TODO: this loop is pretty slow .. (Parellize)
-        for iky in range(self.nky):
-            ky = self.kys[iky]
-            A = self.getA(ky)
-            for src in self.survey.srcList:
-                u_src = f[src, self._solutionType, iky] # solution vector
-                dA_dm_v = self.getADeriv(ky, u_src, v)
-                dRHS_dm_v = self.getRHSDeriv(ky, src, v)
-                du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
-                for rx in src.rxList:
-                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                    df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
-                    # Trapezoidal intergration
-                    Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
-                    if iky==0:
-                        #First assigment
-                        Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y)
-                    else:
-                        Jv[src, rx] += Jv1_temp*dky[iky]    /2.*np.cos(ky*y)
-                        Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y)
-                    Jv0[src, rx] = Jv1_temp.copy()
-        return Utils.mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size, dtype=float)
-
-        # Assume y=0.
-        # This needs some thoughts to implement in general when src is dipole
-        dky = np.diff(self.kys)
-        dky = np.r_[dky[0], dky]
-        y = 0.
-
-        for src in self.survey.srcList:
-            for rx in src.rxList:
-                Jtv_temp1 = np.zeros(m.size, dtype=float)
-                Jtv_temp0 = np.zeros(m.size, dtype=float)
-                #TODO: this loop is pretty slow .. (Parellize)
-                for iky in range(self.nky):
-                    u_src = f[src, self._solutionType, iky]
-                    ky = self.kys[iky]
-                    AT = self.getA(ky)
-                    PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                    df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
-
-                    ATinvdf_duT = self.Ainv[iky] * df_duT
-
-                    dA_dmT = self.getADeriv(ky, u_src, ATinvdf_duT, adjoint=True)
-                    dRHS_dmT = self.getRHSDeriv(ky, src, ATinvdf_duT, adjoint=True)
-                    du_dmT = -dA_dmT + dRHS_dmT
-                    Jtv_temp1 = 1./np.pi*(df_dmT + du_dmT).astype(float)
-                    # Trapezoidal intergration
-                    if iky==0:
-                        #First assigment
-                        Jtv += Jtv_temp1*dky[iky]*np.cos(ky*y)
-                    else:
-                        Jtv += Jtv_temp1*dky[iky]/2.*np.cos(ky*y)
-                        Jtv += Jtv_temp0*dky[iky]/2.*np.cos(ky*y)
-                    Jtv_temp0 = Jtv_temp1.copy()
-        return Utils.mkvc(Jtv)
-
-    def getSourceTerm(self, ky):
-        """
-        takes concept of source and turns it into a matrix
-        """
-        """
-        Evaluates the sources, and puts them in matrix form
-
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: q (nC or nN, nSrc)
-        """
-
-        Srcs = self.survey.srcList
-
-        if self._formulation is 'EB':
-            n = self.mesh.nN
-            # return NotImplementedError
-
-        elif self._formulation is 'HJ':
-            n = self.mesh.nC
-
-        q = np.zeros((n, len(Srcs)))
-
-        for i, src in enumerate(Srcs):
-            q[:,i] = src.eval(self)
-        return q
-
-class Problem2D_CC(BaseDCProblem_2D):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'HJ' # CC potentials means J is on faces
-    fieldsPair    = Fields_ky_CC
-
-    def __init__(self, mesh, **kwargs):
-        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
-        self.setBC()
-
-    def getA(self, ky):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        D = self.Div
-        G = self.Grad
-        vol = self.mesh.vol
-        # TODO: this won't work for full anisotropy
-        MfRhoI = self.MfRhoI
-        # Get resistivity rho
-        rho = self.curModel.rho
-        A = D * MfRhoI * G + Utils.sdiag(ky**2*vol/rho)
-        return A
-
-    def getADeriv(self, ky, u, v, adjoint= False):
-
-        D = self.Div
-        G = self.Grad
-        vol = self.mesh.vol
-        MfRhoIDeriv = self.MfRhoIDeriv
-        rho = self.curModel.rho
-        if adjoint:
-            return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
-        return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
-
-    def getRHS(self, ky):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm(ky)
-        return RHS
-
-    def getRHSDeriv(self, ky, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-    def setBC(self):
-        if self.mesh.dim==3:
-            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-            gBFzm = self.mesh.gridFz[fzm,:]
-            gBFzp = self.mesh.gridFz[fzp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-            beta_zm, beta_zp = temp_zm, temp_zp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
-
-        elif self.mesh.dim==2:
-
-            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
-        V = self.Vol
-        self.Div = V * self.mesh.faceDiv
-        P_BC, B = self.mesh.getBCProjWF_simple()
-        M = B*self.mesh.aveCC2F
-        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
-
-class Problem2D_N(BaseDCProblem_2D):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'EB' # CC potentials means J is on faces
-    fieldsPair    = Fields_ky_N
-
-    def __init__(self, mesh, **kwargs):
-        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
-        # self.setBC()
-
-    @property
-    def MnSigma(self):
-        """
-            Node inner product matrix for \\(\\sigma\\). Used in the E-B formulation
-        """
-        # TODO: only works isotropic sigma
-        sigma = self.curModel.sigma
-        vol = self.mesh.vol
-        MnSigma = Utils.sdiag(self.mesh.aveN2CC.T*(Utils.sdiag(vol)*sigma))
-
-        return MnSigma
-
-    def MnSigmaDeriv(self, u):
-        """
-            Derivative of MnSigma with respect to the model
-        """
-        sigma = self.curModel.sigma
-        sigmaderiv = self.curModel.sigmaDeriv
-        vol = self.mesh.vol
-        return Utils.sdiag(u)*self.mesh.aveN2CC.T*Utils.sdiag(vol) * self.curModel.sigmaDeriv
-
-    def getA(self, ky):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        # TODO: this won't work for full anisotropy
-        MeSigma = self.MeSigma
-        MnSigma = self.MnSigma
-        Grad = self.mesh.nodalGrad
-        # Get conductivity sigma
-        sigma = self.curModel.sigma
-        A = Grad.T * MeSigma * Grad + ky**2*MnSigma
-
-        # Handling Null space of A
-        A[0,0] = A[0,0] + 1.
-        return A
-
-    def getADeriv(self, ky, u, v, adjoint= False):
-
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        sigma = self.curModel.sigma
-        vol = self.mesh.vol
-
-        if adjoint:
-            return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + ky**2*self.MnSigmaDeriv(u).T*v
-        return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + ky**2*self.MnSigmaDeriv(u)*v
-
-    def getRHS(self, ky):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm(ky)
-        return RHS
-
-    def getRHSDeriv(self, ky, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
-        # return qDeriv
-        return Zero()

SimPEG/EM/Static/DC/RxDC.py

@@ -1,129 +0,0 @@
-import SimPEG
-import numpy as np
-from SimPEG.Utils import Zero, closestPoints
-
-class BaseRx(SimPEG.Survey.BaseRx):
-    locs = None
-    rxType = None
-
-    knownRxTypes = {
-                    'phi':['phi',None],
-                    'ex':['e','x'],
-                    'ey':['e','y'],
-                    'ez':['e','z'],
-                    'jx':['j','x'],
-                    'jy':['j','y'],
-                    'jz':['j','z'],
-                    }
-
-    def __init__(self, locs, rxType, **kwargs):
-        SimPEG.Survey.BaseRx.__init__(self, locs, rxType, **kwargs)
-
-
-    @property
-    def projField(self):
-        """Field Type projection (e.g. e b ...)"""
-        return self.knownRxTypes[self.rxType][0]
-
-    def projGLoc(self, f):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        comp = self.knownRxTypes[self.rxType][1]
-        if comp is not None:
-            return f._GLoc(self.rxType) + comp
-        return f._GLoc(self.rxType)
-
-    def eval(self, src, mesh, f):
-        P = self.getP(mesh, self.projGLoc(f))
-        return P*f[src, self.projField]
-
-    def evalDeriv(self, src, mesh, f, v, adjoint=False):
-        P = self.getP(mesh, self.projGLoc(f))
-        if not adjoint:
-            return P*v
-        elif adjoint:
-            return P.T*v
-
-# DC.Rx.Dipole(locs)
-class Dipole(BaseRx):
-
-    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
-        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
-        locs = [locsM, locsN]
-        # We may not need this ...
-        BaseRx.__init__(self, locs, rxType)
-
-    @property
-    def nD(self):
-        """Number of data in the receiver."""
-        return self.locs[0].shape[0]
-
-        # Not sure why ...
-        # return int(self.locs[0].size / 2)
-
-
-    def getP(self, mesh, Gloc):
-        if mesh in self._Ps:
-            return self._Ps[mesh]
-
-        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
-        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
-        P = P0 - P1
-
-        if self.storeProjections:
-            self._Ps[mesh] = P
-
-        return P
-
-
-class Dipole_ky(BaseRx):
-
-    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
-        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
-        locs = [locsM, locsN]
-        # We may not need this ...
-        BaseRx.__init__(self, locs, rxType)
-
-    @property
-    def nD(self):
-        """Number of data in the receiver."""
-        return self.locs[0].shape[0]
-
-        # Not sure why ...
-        # return int(self.locs[0].size / 2)
-
-    def getP(self, mesh, Gloc):
-        if mesh in self._Ps:
-            return self._Ps[mesh]
-
-        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
-        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
-        P = P0 - P1
-        if self.storeProjections:
-            self._Ps[mesh] = P
-        return P
-
-    def eval(self, kys, src, mesh, f):
-        P = self.getP(mesh, self.projGLoc(f))
-        Pf = P*f[src, self.projField,:]
-        return self.IntTrapezoidal(kys, Pf, y=0.)
-
-    def evalDeriv(self, ky, src, mesh, f, v, adjoint=False):
-        P = self.getP(mesh, self.projGLoc(f))
-        if not adjoint:
-            return P*v
-        elif adjoint:
-            return P.T*v
-
-    def IntTrapezoidal(self, kys, Pf, y=0.):
-        phi = np.zeros(Pf.shape[0])
-        nky = kys.size
-        dky = np.diff(kys)
-        dky = np.r_[dky[0], dky]
-        phi0 = 1./np.pi*Pf[:,0]
-        for iky in range(nky):
-            phi1 = 1./np.pi*Pf[:,iky]
-            phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y)
-            phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y)
-            phi0 = phi1.copy()
-        return phi
-

SimPEG/EM/Static/DC/SrcDC.py

@@ -1,86 +0,0 @@
-import SimPEG
-# from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG.Utils import Zero, closestPoints, mkvc
-import numpy as np
-
-class BaseSrc(SimPEG.Survey.BaseSrc):
-
-    current = 1.0
-    loc = None
-
-    def __init__(self, rxList, **kwargs):
-        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
-
-    def eval(self, prob):
-        raise NotImplementedError
-
-    def evalDeriv(self, prob):
-        return Zero()
-
-
-class Dipole(BaseSrc):
-
-    def __init__(self, rxList, locA, locB, **kwargs):
-        assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
-        self.loc = [locA, locB]
-        BaseSrc.__init__(self, rxList, **kwargs)
-
-    def eval(self, prob):
-        if prob._formulation == 'HJ':
-            inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
-            q = np.zeros(prob.mesh.nC)
-            q[inds] = self.current * np.r_[1., -1.]
-        elif prob._formulation == 'EB':
-            qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
-            qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
-            q = self.current * mkvc(qa+qb)
-        return q
-
-class Pole(BaseSrc):
-
-    def __init__(self, rxList, loc, **kwargs):
-        BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
-
-    def eval(self, prob):
-        if prob._formulation == 'HJ':
-            inds = closestPoints(prob.mesh, self.loc)
-            q = np.zeros(prob.mesh.nC)
-            q[inds] = self.current * np.r_[1.]
-        elif prob._formulation == 'EB':
-            q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
-            q = self.current * mkvc(q)
-        return q
-
-
-# class Dipole_ky(BaseSrc):
-
-#     def __init__(self, rxList, locA, locB, **kwargs):
-#         assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
-#         self.loc = [locA[[0,2]], locB[[0,2]]]
-#         BaseSrc.__init__(self, rxList, **kwargs)
-
-#     def eval(self, prob):
-#         if prob._formulation == 'HJ':
-#             inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
-#             q = np.zeros(prob.mesh.nC)
-#             q[inds] = self.current * np.r_[1., -1.]
-#         elif prob._formulation == 'EB':
-#             qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
-#             qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
-#             q = self.current * mkvc(qa+qb)
-#         return q
-
-# class Pole_ky(BaseSrc):
-
-#     def __init__(self, rxList, loc, **kwargs):
-#         BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
-
-#     def eval(self, prob):
-#         if prob._formulation == 'HJ':
-#             inds = closestPoints(prob.mesh, self.loc[[0,2]])
-#             q = np.zeros(prob.mesh.nC)
-#             q[inds] = self.current * np.r_[1.]
-#         elif prob._formulation == 'EB':
-#             q = prob.mesh.getInterpolationMat(self.loc[[0,2]], locType='N').todense()
-#             q = self.current * mkvc(q)
-#         return q

SimPEG/EM/Static/DC/SurveyDC.py

@@ -1,38 +0,0 @@
-import SimPEG
-from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG import sp, Survey
-from SimPEG.Utils import Zero, Identity
-from RxDC import BaseRx
-from SrcDC import BaseSrc
-
-class Survey(BaseEMSurvey):
-    rxPair  = BaseRx
-    srcPair = BaseSrc
-
-    def __init__(self, srcList, **kwargs):
-        self.srcList = srcList
-        BaseEMSurvey.__init__(self, srcList, **kwargs)
-
-class Survey_ky(BaseEMSurvey):
-    rxPair  = BaseRx
-    srcPair = BaseSrc
-
-    def __init__(self, srcList, **kwargs):
-        self.srcList = srcList
-        BaseEMSurvey.__init__(self, srcList, **kwargs)
-
-    def eval(self, f):
-        """
-        Project fields to receiver locations
-        :param Fields u: fields object
-        :rtype: numpy.ndarray
-        :return: data
-        """
-        data = SimPEG.Survey.Data(self)
-        kys = self.prob.kys
-        for src in self.srcList:
-            for rx in src.rxList:
-                data[src, rx] = rx.eval(kys, src, self.mesh, f)
-        return data
-
-

SimPEG/EM/Static/DC/Utils.py

@@ -1,38 +0,0 @@
-import numpy as np
-
-def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
-
-    import SimPEG.EM.Static.DC as DC
-
-    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
-    elocs -= (nElecs*aSpacing - aSpacing)/2
-    space = 1
-    WENNER = np.zeros((0,),dtype=int)
-    for ii in range(nElecs):
-        for jj in range(nElecs):
-            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
-            if np.any(test >= nElecs):
-                break
-            WENNER = np.r_[WENNER, test]
-        space += 1
-    WENNER = WENNER.reshape((-1,4))
-
-
-    if plotIt:
-        for i, s in enumerate('rbkg'):
-            plt.plot(elocs[WENNER[:,i]],s+'.')
-        plt.show()
-
-    # Create sources and receivers
-    i = 0
-    if in2D:
-        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
-    else:
-        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
-    srcList = []
-    for i in range(WENNER.shape[0]):
-        rx = DC.Rx.Dipole(getLoc(i,1).reshape([1,-1]),getLoc(i,2).reshape([1,-1]))
-        src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3))
-        srcList += [src]
-
-    return srcList

SimPEG/EM/Static/DC/__init__.py

@@ -1,8 +0,0 @@
-from ProblemDC import Problem3D_CC, Problem3D_N
-from ProblemDC_2D import Problem2D_CC, Problem2D_N
-from SurveyDC import Survey, Survey_ky
-import SrcDC as Src #Pole
-import RxDC as Rx
-from FieldsDC import Fields_CC
-from BoundaryUtils import getxBCyBC_CC
-import Utils

SimPEG/EM/Static/IP/ProblemIP.py

@@ -1,374 +0,0 @@
-from SimPEG import Problem, Utils, Maps, Mesh
-from SimPEG.EM.Base import BaseEMProblem
-from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
-from SimPEG.Utils import sdiag
-import numpy as np
-from SimPEG.Utils import Zero
-from SimPEG.EM.Static.DC import getxBCyBC_CC
-from SurveyIP import Survey
-
-class IPPropMap(Maps.PropMap):
-    """
-        Property Map for IP Problems. The electrical chargeability,
-        (\\(\\eta\\)) is the default inversion property
-    """
-    eta = Maps.Property("Electrical Chargeability", defaultInvProp = True)
-
-class BaseIPProblem(BaseEMProblem):
-
-    surveyPair = Survey
-    fieldsPair = Fields
-    PropMap = IPPropMap
-    Ainv = None
-    sigma = None
-    rho = None
-    f = None
-    Ainv = None
-
-    def fields(self, m):
-        self.curModel = m
-        if self.f is None:
-            self.f = self.fieldsPair(self.mesh, self.survey)
-            if self.Ainv == None:
-                A = self.getA()
-                self.Ainv = self.Solver(A, **self.solverOpts)
-            RHS = self.getRHS()
-            u = self.Ainv * RHS
-            Srcs = self.survey.srcList
-            self.f[Srcs, self._solutionType] = u
-        return self.f
-
-    def Jvec(self, m, v, f=None):
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        Jv = self.dataPair(self.survey) #same size as the data
-
-        A = self.getA()
-
-        for src in self.survey.srcList:
-            u_src = f[src, self._solutionType] # solution vector
-            dA_dm_v = self.getADeriv(u_src, v)
-            dRHS_dm_v = self.getRHSDeriv(src, v)
-            du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
-
-            for rx in src.rxList:
-                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
-        # Conductivity (d u / d log sigma)
-        if self._formulation is 'EB':
-            return -Utils.mkvc(Jv)
-        # Conductivity (d u / d log rho)
-        if self._formulation is 'HJ':
-            return Utils.mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-        AT = self.getA()
-
-        for src in self.survey.srcList:
-            u_src = f[src, self._solutionType]
-            for rx in src.rxList:
-                PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
-                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
-                ATinvdf_duT = self.Ainv * df_duT
-                dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
-                dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
-                du_dmT = -dA_dmT + dRHS_dmT
-                Jtv += (df_dmT + du_dmT).astype(float)
-        # Conductivity ((d u / d log sigma).T)
-        if self._formulation is 'EB':
-            return -Utils.mkvc(Jtv)
-        # Conductivity ((d u / d log rho).T)
-        if self._formulation is 'HJ':
-            return Utils.mkvc(Jtv)
-
-    def getSourceTerm(self):
-        """
-        takes concept of source and turns it into a matrix
-        """
-        """
-        Evaluates the sources, and puts them in matrix form
-
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: q (nC or nN, nSrc)
-        """
-
-        Srcs = self.survey.srcList
-
-        if self._formulation is 'EB':
-            n = self.mesh.nN
-            # return NotImplementedError
-
-        elif self._formulation is 'HJ':
-            n = self.mesh.nC
-
-        q = np.zeros((n, len(Srcs)))
-
-        for i, src in enumerate(Srcs):
-            q[:,i] = src.eval(self)
-        return q
-
-    @property
-    def deleteTheseOnModelUpdate(self):
-        toDelete = []
-        return toDelete
-
-    # assume log rho or log cond
-    @property
-    def MeSigma(self):
-        """
-            Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
-        """
-        if getattr(self, '_MeSigma', None) is None:
-            self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
-        return self._MeSigma
-
-    @property
-    def MfRhoI(self):
-        """
-            Inverse of :code:`MfRho`
-        """
-        if getattr(self, '_MfRhoI', None) is None:
-            self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
-        return self._MfRhoI
-
-    def MfRhoIDeriv(self,u):
-        """
-            Derivative of :code:`MfRhoI` with respect to the model.
-        """
-
-        dMfRhoI_dI = -self.MfRhoI**2
-        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
-        drho_dlogrho = Utils.sdiag(self.rho)*self.curModel.etaDeriv
-        return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
-
-    # TODO: This should take a vector
-    def MeSigmaDeriv(self, u):
-        """
-            Derivative of MeSigma with respect to the model
-        """
-        dsigma_dlogsigma = Utils.sdiag(self.sigma)*self.curModel.etaDeriv
-        return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
-
-class Problem3D_CC(BaseIPProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'HJ' # CC potentials means J is on faces
-    fieldsPair    = Fields_CC
-
-    def __init__(self, mesh, **kwargs):
-        BaseIPProblem.__init__(self, mesh, **kwargs)
-        self.setBC()
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        D = self.Div
-        G = self.Grad
-        # TODO: this won't work for full anisotropy
-        MfRhoI = self.MfRhoI
-        A = D * MfRhoI * G
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * A
-        return A
-
-    def getADeriv(self, u, v, adjoint= False):
-
-        D = self.Div
-        G = self.Grad
-        MfRhoIDeriv = self.MfRhoIDeriv
-
-        if adjoint:
-            # if self._makeASymmetric is True:
-            #     v = V * v
-            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
-        return D * (MfRhoIDeriv( G * u ) * v)
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return self.Vol.T * RHS
-
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-    def setBC(self):
-        if self.mesh.dim==3:
-            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-            gBFzm = self.mesh.gridFz[fzm,:]
-            gBFzp = self.mesh.gridFz[fzp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-            beta_zm, beta_zp = temp_zm, temp_zp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
-
-        elif self.mesh.dim==2:
-
-            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
-        V = self.Vol
-        self.Div = V * self.mesh.faceDiv
-        P_BC, B = self.mesh.getBCProjWF_simple()
-        M = B*self.mesh.aveCC2F
-        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
-
-
-class Problem3D_N(BaseIPProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'EB' # N potentials means B is on faces
-    fieldsPair    = Fields_N
-
-    def __init__(self, mesh, **kwargs):
-        BaseIPProblem.__init__(self, mesh, **kwargs)
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        # TODO: this won't work for full anisotropy
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        A = Grad.T * MeSigma * Grad
-
-        # Handling Null space of A
-        A[0,0] = A[0,0] + 1.
-
-        return A
-
-    def getADeriv(self, u, v, adjoint=False):
-        """
-
-        Product of the derivative of our system matrix with respect to the model and a vector
-
-        """
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        if not adjoint:
-            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
-        elif adjoint:
-            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
-
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-if __name__ == '__main__':
-
-
-    cs = 12.5
-    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-    hz = [(cs,7, -1.3),(cs,20)]
-    mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-    sigma = np.ones(mesh.nC)
-    prob = BaseIPProblem(mesh, sigma=sigma)
-
-

SimPEG/EM/Static/IP/SurveyIP.py

@@ -1,23 +0,0 @@
-import SimPEG
-from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG import sp, Survey
-from SimPEG.Utils import Zero, Identity
-from SimPEG.EM.Static.DC.SrcDC import BaseSrc
-from SimPEG.EM.Static.DC.RxDC import BaseRx
-
-class Survey(BaseEMSurvey):
-    rxPair  = BaseRx
-    srcPair = BaseSrc
-
-    def __init__(self, srcList, **kwargs):
-        self.srcList = srcList
-        BaseEMSurvey.__init__(self, srcList, **kwargs)
-
-    def dpred(self, m, f=None):
-        """
-            Predicted data.
-
-            .. math::
-                d_\\text{pred} = Pf(m)
-        """
-        return self.prob.Jvec(m, m, f=f)

SimPEG/EM/Static/IP/__init__.py

@@ -1,2 +0,0 @@
-from ProblemIP import Problem3D_CC, Problem3D_N
-from SurveyIP import Survey

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -1,445 +0,0 @@
-from SimPEG import Problem, Utils, Maps, Mesh
-from SimPEG.EM.Base import BaseEMProblem
-from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
-from SimPEG.Utils import sdiag
-import numpy as np
-from SimPEG.Utils import Zero
-from SimPEG.EM.Static.DC import getxBCyBC_CC
-from SurveySIP import Survey, Data
-
-class ColeColePropMap(Maps.PropMap):
-    """
-        Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
-    """
-
-    eta = Maps.Property("Electrical Conductivity", defaultInvProp=True)
-    tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap))
-    taui   = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap))
-    c   = Maps.Property("Electrical Conductivity", defaultVal=1.)
-
-
-class BaseSIPProblem(BaseEMProblem):
-
-    surveyPair = Survey
-    fieldsPair = Fields
-    dataPair = Data
-    PropMap = ColeColePropMap
-    Ainv = None
-    sigma = None
-    rho = None
-    f = None
-    Ainv = None
-
-    def DebyeTime(self, t):
-        peta = self.curModel.eta*np.exp(-self.curModel.taui*t)
-        return peta
-
-    def EtaDeriv(self, t, v, adjoint=False):
-        v = np.array(v, dtype=float)
-        if adjoint:
-            return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v)
-        else:
-            return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v)
-
-
-    def TauiDeriv(self, t, v, adjoint=False):
-        v = np.array(v, dtype=float)
-        if adjoint:
-            return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v)
-        else:
-            return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v)
-
-    def fields(self, m):
-        self.curModel = m
-        if self.f is None:
-            self.f = self.fieldsPair(self.mesh, self.survey)
-            if self.Ainv == None:
-                A = self.getA()
-                self.Ainv = self.Solver(A, **self.solverOpts)
-            RHS = self.getRHS()
-            u = self.Ainv * RHS
-            Srcs = self.survey.srcList
-            self.f[Srcs, self._solutionType] = u
-        return self.f
-
-    def forward(self, m, f=None):
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-        Jv = self.dataPair(self.survey) #same size as the data
-        # A = self.getA()
-        JvAll = []
-        for tind in range(len(self.survey.times)):
-            #Pseudo-chareability
-            t = self.survey.times[tind]
-            v = self.DebyeTime(t)
-            for src in self.survey.srcList:
-                u_src = f[src, self._solutionType] # solution vector
-                dA_dm_v = self.getADeriv(u_src, v)
-                dRHS_dm_v = self.getRHSDeriv(src, v)
-                du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
-                for rx in src.rxList:
-                    timeindex = rx.getTimeP(self.survey.times)
-                    if timeindex[tind]:
-                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                        df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                        Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
-
-        # Conductivity (d u / d log sigma)
-        if self._formulation is 'EB':
-            return -Utils.mkvc(Jv)
-        # Resistivity (d u / d log rho)
-        if self._formulation is 'HJ':
-            return Utils.mkvc(Jv)
-
-    def Jvec(self, m, v, f=None):
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-        Jv = self.dataPair(self.survey) #same size as the data
-        # A = self.getA()
-        JvAll = []
-        #Assume only eta and tau (eta first then tau)
-        # v = [2*Mx1]
-        v = v.reshape((int(v.size/2), 2), order='F')
-
-        for tind in range(len(self.survey.times)):
-            t = self.survey.times[tind]
-            v0 = self.EtaDeriv(t, v[:,0])
-            v1 = self.TauiDeriv(t, v[:,1])
-            for src in self.survey.srcList:
-                u_src = f[src, self._solutionType] # solution vector
-                dA_dm_v0 = self.getADeriv(u_src, v0)
-                dRHS_dm_v0 = self.getRHSDeriv(src, v0)
-                du_dm_v0 = self.Ainv * ( - dA_dm_v0 + dRHS_dm_v0 )
-                dA_dm_v1 = self.getADeriv(u_src, v1)
-                dRHS_dm_v1 = self.getRHSDeriv(src, v1)
-                du_dm_v1 = self.Ainv * ( - dA_dm_v1 + dRHS_dm_v1 )
-                for rx in src.rxList:
-                    timeindex = rx.getTimeP(self.survey.times)
-                    if timeindex[tind]:
-                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                        df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
-                        df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
-                        Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
-                        Jv[src, rx, t] += rx.evalDeriv(src, self.mesh, f, df_dm_v1)
-        # Conductivity (d u / d log sigma)
-        if self._formulation is 'EB':
-            return -Jv.tovec()
-        # Resistivity (d u / d log rho)
-        if self._formulation is 'HJ':
-            return Jv.tovec()
-
-    def Jtvec(self, m, v, f=None):
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv= np.zeros(m.size)
-        for tind in range(len(self.survey.times)):
-            t = self.survey.times[tind]
-            for src in self.survey.srcList:
-                u_src = f[src, self._solutionType]
-                for rx in src.rxList:
-                    timeindex = rx.getTimeP(self.survey.times)
-                    if timeindex[tind]:
-                        PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
-                        df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
-                        df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
-                        ATinvdf_duT = self.Ainv * df_duT
-                        dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
-                        dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
-                        du_dmT = -dA_dmT + dRHS_dmT
-                        Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)]
-
-        # Conductivity ((d u / d log sigma).T)
-        if self._formulation is 'EB':
-            return -Jtv
-        # Conductivity ((d u / d log rho).T)
-        if self._formulation is 'HJ':
-            return Jtv
-
-    def getSourceTerm(self):
-        """
-        takes concept of source and turns it into a matrix
-        """
-        """
-        Evaluates the sources, and puts them in matrix form
-
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: q (nC or nN, nSrc)
-        """
-
-        Srcs = self.survey.srcList
-
-        if self._formulation is 'EB':
-            n = self.mesh.nN
-            # return NotImplementedError
-
-        elif self._formulation is 'HJ':
-            n = self.mesh.nC
-
-        q = np.zeros((n, len(Srcs)))
-
-        for i, src in enumerate(Srcs):
-            q[:,i] = src.eval(self)
-        return q
-
-    @property
-    def deleteTheseOnModelUpdate(self):
-        toDelete = []
-        return toDelete
-
-    # assume log rho or log cond
-    @property
-    def MeSigma(self):
-        """
-            Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
-        """
-        if getattr(self, '_MeSigma', None) is None:
-            self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
-        return self._MeSigma
-
-    @property
-    def MfRhoI(self):
-        """
-            Inverse of :code:`MfRho`
-        """
-        if getattr(self, '_MfRhoI', None) is None:
-            self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
-        return self._MfRhoI
-
-    def MfRhoIDeriv(self,u):
-        """
-            Derivative of :code:`MfRhoI` with respect to the model.
-        """
-
-        dMfRhoI_dI = -self.MfRhoI**2
-        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
-        drho_dlogrho = Utils.sdiag(self.rho)
-        return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
-
-    # TODO: This should take a vector
-    def MeSigmaDeriv(self, u):
-        """
-            Derivative of MeSigma with respect to the model
-        """
-        dsigma_dlogsigma = Utils.sdiag(self.sigma)
-        return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
-
-class Problem3D_CC(BaseSIPProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'HJ' # CC potentials means J is on faces
-    fieldsPair    = Fields_CC
-
-    def __init__(self, mesh, **kwargs):
-        BaseSIPProblem.__init__(self, mesh, **kwargs)
-        self.setBC()
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        D = self.Div
-        G = self.Grad
-        # TODO: this won't work for full anisotropy
-        MfRhoI = self.MfRhoI
-        A = D * MfRhoI * G
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * A
-        return A
-
-    def getADeriv(self, u, v, adjoint= False):
-
-        D = self.Div
-        G = self.Grad
-        MfRhoIDeriv = self.MfRhoIDeriv
-
-        if adjoint:
-            # if self._makeASymmetric is True:
-            #     v = V * v
-            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
-        return D * (MfRhoIDeriv( G * u ) * v)
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return self.Vol.T * RHS
-
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-    def setBC(self):
-        if self.mesh.dim==3:
-            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-            gBFzm = self.mesh.gridFz[fzm,:]
-            gBFzp = self.mesh.gridFz[fzp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-            beta_zm, beta_zp = temp_zm, temp_zp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
-
-        elif self.mesh.dim==2:
-
-            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
-            gBFxm = self.mesh.gridFx[fxm,:]
-            gBFxp = self.mesh.gridFx[fxp,:]
-            gBFym = self.mesh.gridFy[fym,:]
-            gBFyp = self.mesh.gridFy[fyp,:]
-
-            # Setup Mixed B.C (alpha, beta, gamma)
-            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-
-            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
-            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
-
-            beta_xm, beta_xp = temp_xm, temp_xp
-            beta_ym, beta_yp = temp_ym, temp_yp
-
-            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
-            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
-
-            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
-            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
-            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
-        V = self.Vol
-        self.Div = V * self.mesh.faceDiv
-        P_BC, B = self.mesh.getBCProjWF_simple()
-        M = B*self.mesh.aveCC2F
-        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
-
-
-class Problem3D_N(BaseSIPProblem):
-
-    _solutionType = 'phiSolution'
-    _formulation  = 'EB' # N potentials means B is on faces
-    fieldsPair    = Fields_N
-
-    def __init__(self, mesh, **kwargs):
-        BaseSIPProblem.__init__(self, mesh, **kwargs)
-
-    def getA(self):
-        """
-
-        Make the A matrix for the cell centered DC resistivity problem
-
-        A = D MfRhoI D^\\top V
-
-        """
-
-        # TODO: this won't work for full anisotropy
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        A = Grad.T * MeSigma * Grad
-
-        # Handling Null space of A
-        A[0,0] = A[0,0] + 1.
-
-        return A
-
-    def getADeriv(self, u, v, adjoint=False):
-        """
-
-        Product of the derivative of our system matrix with respect to the model and a vector
-
-        """
-        MeSigma = self.MeSigma
-        Grad = self.mesh.nodalGrad
-        if not adjoint:
-            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
-        elif adjoint:
-            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
-
-
-    def getRHS(self):
-        """
-        RHS for the DC problem
-
-        q
-        """
-
-        RHS = self.getSourceTerm()
-        return RHS
-
-    def getRHSDeriv(self, src, v, adjoint=False):
-        """
-        Derivative of the right hand side with respect to the model
-        """
-        # TODO: add qDeriv for RHS depending on m
-        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        # return qDeriv
-        return Zero()
-
-if __name__ == '__main__':
-
-
-    cs = 12.5
-    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-    hz = [(cs,7, -1.3),(cs,20)]
-    mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-    sigma = np.ones(mesh.nC)
-    prob = BaseSIPProblem(mesh, sigma=sigma)
-
-

SimPEG/EM/Static/SIP/Regularization.py

@@ -1,204 +0,0 @@
-from SimPEG import Utils, Maps, Mesh, sp, np
-from SimPEG.Regularization import BaseRegularization, Simple
-
-class MultiRegularization(Simple):
-    """
-    **MultiRegularization Class**
-
-    This is used to regularize the model space
-    having multiple models [m1, m2, m3, ...] ::
-
-        reg = Regularization(mesh)
-
-    """
-    nModels = None # Number of models
-    ratios = None  # Ratio for different models
-    crossgrad = False # Use cross gradient or not
-    betacross = 1.
-    wx = []
-    wy = []
-    wz = []
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-        if self.nModels == None:
-            raise Exception("Put nModels as a initial input!")
-        if self.ratios == None:
-            self.ratios = [1. for imodel in range(self.nModels)]
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self,'_Wsmall', None) is None:
-            vecs = []
-            for imodel in range(self.nModels):
-                vecs.append((self.regmesh.vol*self.alpha_s*self.wght*self.ratios[imodel])**0.5)
-            self._Wsmall = Utils.sdiag(np.hstack(vecs))
-        return self._Wsmall
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            mats = []
-            for imodel in range(self.nModels):
-                self.wx.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5))
-                mats.append(self.wx[imodel]*self.regmesh.cellDiffxStencil)
-            self._Wx = sp.block_diag(mats)
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            mats = []
-            for imodel in range(self.nModels):
-                self.wy.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5))
-                mats.append(self.wy[imodel]*self.regmesh.cellDiffyStencil)
-            self._Wy = sp.block_diag(mats)
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            mats = []
-            for imodel in range(self.nModels):
-                self.wz.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5))
-                mats.append(self.wz[imodel]*self.regmesh.cellDiffzStencil)
-            self._Wz = sp.block_diag(mats)
-        return self._Wz
-
-    @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 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
-
-
-    @Utils.timeIt
-    def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @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)
-
-    def cross(a,b):
-        ax, ay, az = a[0], a[1], a[2]
-        bx, by, bz = b[0], b[1], b[2]
-        cx = ay*bz - az*by
-        cy = az*bx - ax*bz
-        cz = ax*by - ay*bx
-        return [cx, cy, cz]
-
-     # TODO: Implement Cross Gradients..
-    @Utils.timeIt
-    def _evalCross(self, m):
-        if self.crossgrad == False:
-            return 0.
-        elif self.crossgrad == True:
-            M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F")
-
-            ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0]
-            ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0]
-            az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0]
-            bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
-            by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
-            bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
-            #ab
-            out_ab = cross([ax, ay, az], [bx, by, bz])
-            r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross)
-
-            if self.nModels == 3:
-                cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
-                cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
-                cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
-                #ac
-                out_ac = cross([ax, ay, az], [cx, cy, cz])
-                #bc
-                out_bc = cross([bx, by, bz], [cx, cy, cz])
-                r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)]
-
-            return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def evalDeriv(self, m):
-        """
-        The regularization is:
-
-        .. math::
-
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
-
-        So the derivative is straight forward:
-
-        .. math::
-
-            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
-        """
-        deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
-        if self.crossgrad==True:
-            deriv += self._evalCrossDeriv(m)
-        return deriv
-
-    @Utils.timeIt
-    def _evalCrossDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @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
-
-        The regularization is:
-
-        .. math::
-
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
-
-        So the second derivative is straight forward:
-
-        .. math::
-
-            R(m) = \mathbf{W^\\top W}
-
-        """
-        mD = self.mapping.deriv(m - self.mref)
-        if v is None:
-            return mD.T * self.W.T * self.W * mD
-
-        return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
-
-
-

SimPEG/EM/Static/SIP/RxSIP.py

@@ -1,88 +0,0 @@
-import SimPEG
-import numpy as np
-from SimPEG.Utils import Zero, closestPoints
-
-class BaseRx(SimPEG.Survey.BaseTimeRx):
-    locs = None
-    rxType = None
-
-    knownRxTypes = {
-                    'phi':['phi',None],
-                    'ex':['e','x'],
-                    'ey':['e','y'],
-                    'ez':['e','z'],
-                    'jx':['j','x'],
-                    'jy':['j','y'],
-                    'jz':['j','z'],
-                    }
-
-    def __init__(self, locs, times, rxType, **kwargs):
-        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs)
-
-    @property
-    def projField(self):
-        """Field Type projection (e.g. e b ...)"""
-        return self.knownRxTypes[self.rxType][0]
-
-    def projGLoc(self, f):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        comp = self.knownRxTypes[self.rxType][1]
-        if comp is not None:
-            return f._GLoc(self.rxType) + comp
-        return f._GLoc(self.rxType)
-
-    def getTimeP(self, timesall):
-        """
-            Returns the time projection matrix.
-
-            .. note::
-
-                This is not stored in memory, but is created on demand.
-        """
-        time_inds = np.in1d(timesall, self.times)
-        return time_inds
-
-    def evalDeriv(self, src, mesh, f, v, adjoint=False):
-        P = self.getP(mesh, self.projGLoc(f))
-        if not adjoint:
-            return P*v
-        elif adjoint:
-            return P.T*v
-
-
-# DC.Rx.Dipole(locs)
-class Dipole(BaseRx):
-
-    def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs):
-        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
-        locs = [locsM, locsN]
-        # We may not need this ...
-        BaseRx.__init__(self, locs, times, rxType)
-
-    @property
-    def nD(self):
-        """Number of data in the receiver."""
-        # return self.locs[0].shape[0] * len(self.times)
-        return self.locs[0].shape[0]
-
-    @property
-    def nRx(self):
-        """Number of data in the receiver."""
-        return self.locs[0].shape[0]
-
-        # Not sure why ...
-        # return int(self.locs[0].size / 2)
-
-
-    def getP(self, mesh, Gloc):
-        if mesh in self._Ps:
-            return self._Ps[mesh]
-
-        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
-        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
-        P = P0 - P1
-
-        if self.storeProjections:
-            self._Ps[mesh] = P
-
-        return P

SimPEG/EM/Static/SIP/SrcSIP.py

@@ -1,64 +0,0 @@
-import SimPEG
-# from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG.Utils import Zero, closestPoints, mkvc
-import numpy as np
-
-class BaseSrc(SimPEG.Survey.BaseSrc):
-
-    current = 1.0
-    loc = None
-
-    def __init__(self, rxList, **kwargs):
-        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
-
-    def eval(self, prob):
-        raise NotImplementedError
-
-    def evalDeriv(self, prob):
-        return Zero()
-
-    @property
-    def nD(self):
-        """Number of data"""
-        return self.vnD.sum()
-
-    @property
-    def vnD(self):
-        """Vector number of data"""
-        return np.array([rx.nD*len(rx.times) for rx in self.rxList])
-
-
-
-class Dipole(BaseSrc):
-
-    def __init__(self, rxList, locA, locB, **kwargs):
-        assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
-        self.loc = [locA, locB]
-        BaseSrc.__init__(self, rxList, **kwargs)
-
-    def eval(self, prob):
-        if prob._formulation == 'HJ':
-            inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
-            q = np.zeros(prob.mesh.nC)
-            q[inds] = self.current * np.r_[1., -1.]
-        elif prob._formulation == 'EB':
-            qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
-            qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
-            q = self.current * mkvc(qa+qb)
-        return q
-
-class Pole(BaseSrc):
-
-    def __init__(self, rxList, loc, **kwargs):
-        BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
-
-    def eval(self, prob):
-        if prob._formulation == 'HJ':
-            inds = closestPoints(prob.mesh, self.loc)
-            q = np.zeros(prob.mesh.nC)
-            q[inds] = self.current * np.r_[1.]
-        elif prob._formulation == 'EB':
-            q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
-            q = self.current * mkvc(q)
-        return q
-

SimPEG/EM/Static/SIP/SurveySIP.py

@@ -1,102 +0,0 @@
-import SimPEG
-from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG import np, sp, Survey, Utils
-from SimPEG.Utils import Zero, Identity
-from SimPEG.EM.Static.SIP.SrcSIP import BaseSrc
-from SimPEG.EM.Static.SIP.RxSIP import BaseRx
-import uuid
-
-
-class Survey(BaseEMSurvey):
-    rxPair  = BaseRx
-    srcPair = BaseSrc
-    times = None
-
-    def __init__(self, srcList, **kwargs):
-        self.srcList = srcList
-        BaseEMSurvey.__init__(self, srcList, **kwargs)
-        self.getUniqueTimes()
-
-    def getUniqueTimes(self):
-        time_rx = []
-        for src in self.srcList:
-            for rx in src.rxList:
-                time_rx.append(rx.times)
-        self.times = np.unique(np.hstack(time_rx))
-
-    def dpred(self, m, f=None):
-        """
-            Predicted data.
-
-            .. math::
-                d_\\text{pred} = Pf(m)
-        """
-        return self.prob.forward(m, f=f)
-
-
-class Data(SimPEG.Survey.Data):
-    """Fancy data storage by Src and Rx"""
-
-    def __init__(self, survey, v=None):
-        self.uid = str(uuid.uuid4())
-        self.survey = survey
-        self._dataDict = {}
-        for src in self.survey.srcList:
-            self._dataDict[src] = {}
-            for rx in src.rxList:
-                self._dataDict[src][rx] = {}
-
-        if v is not None:
-            self.fromvec(v)
-
-    def _ensureCorrectKey(self, key):
-        if type(key) is tuple:
-            if len(key) is not 3:
-                raise KeyError('Key must be [Src, Rx, tInd]')
-            if key[0] not in self.survey.srcList:
-                raise KeyError('Src Key must be a source in the survey.')
-            if key[1] not in key[0].rxList:
-                raise KeyError('Rx Key must be a receiver for the source.')
-            return key
-        elif isinstance(key, self.survey.srcPair):
-            if key not in self.survey.srcList:
-                raise KeyError('Key must be a source in the survey.')
-            return key, None, None
-        else:
-            raise KeyError('Key must be [Src] or [Src,Rx] or [Src, Rx, tInd]')
-
-    def __setitem__(self, key, value):
-        src, rx, t = self._ensureCorrectKey(key)
-        assert rx is not None, 'set data using [Src, Rx]'
-        assert isinstance(value, np.ndarray), 'value must by ndarray'
-        assert value.size == rx.nD, "value must have the same number of data as the source."
-        self._dataDict[src][rx][t] = Utils.mkvc(value)
-
-    def __getitem__(self, key):
-        src, rx, t = self._ensureCorrectKey(key)
-        if rx is not None:
-            if rx not in self._dataDict[src]:
-                raise Exception('Data for receiver has not yet been set.')
-            return self._dataDict[src][rx][t]
-
-        return np.concatenate([self[src,rx, t] for rx in src.rxList])
-
-    def tovec(self):
-        val = []
-        for src in self.survey.srcList:
-            for rx in src.rxList:
-                for t in rx.times:
-                    val.append(self[src, rx, t])
-        return np.concatenate(val)
-
-
-    def fromvec(self, v):
-        v = Utils.mkvc(v)
-        assert v.size == self.survey.nD, 'v must have the correct number of data.'
-        indBot, indTop = 0, 0
-        for src in self.survey.srcList:
-            for rx in src.rxList:
-                for t in rx.times:
-                    indTop += rx.nRx
-                    self[src, rx, t] = v[indBot:indTop]
-                    indBot += rx.nRx

SimPEG/EM/Static/SIP/__init__.py

@@ -1,5 +0,0 @@
-from ProblemSIP import Problem3D_CC, Problem3D_N
-from SurveySIP import Survey, Data
-import SrcSIP as Src #Pole
-import RxSIP as Rx
-from Regularization import MultiRegularization

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -1,317 +0,0 @@
-from SimPEG import np
-from SimPEG.EM.Static import DC, IP
-
-def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
-    """
-        Read list of 2D tx-rx location and plot a speudo-section of apparent
-        resistivity.
-
-        Assumes flat topo for now...
-
-        Input:
-        :param d2D, z0
-        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
-        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
-        Output:
-        :figure scatter plot overlayed on image
-
-        Edited Feb 17th, 2016
-
-        @author: dominiquef
-
-    """
-    from SimPEG import np
-    from scipy.interpolate import griddata
-    import pylab as plt
-
-    # Set depth to 0 for now
-    z0 = 0.
-
-    # Pre-allocate
-    midx = []
-    midz = []
-    rho = []
-    LEG = []
-    count = 0 # Counter for data
-    for ii in range(DCsurvey.nSrc):
-
-        Tx = DCsurvey.srcList[ii].loc
-        Rx = DCsurvey.srcList[ii].rxList[0].locs
-
-        nD = DCsurvey.srcList[ii].rxList[0].nD
-
-        data = DCsurvey.dobs[count:count+nD]
-        count += nD
-
-        # Get distances between each poles A-B-M-N
-        if stype == 'pdp':
-            MA = np.abs(Tx[0] - Rx[0][:,0])
-            NA = np.abs(Tx[0] - Rx[1][:,0])
-            MN = np.abs(Rx[1][:,0] - Rx[0][:,0])
-
-            # Create mid-point location
-            Cmid = Tx[0]
-            Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
-            if DCsurvey.mesh.dim == 2:
-                zsrc = Tx[1]
-            elif DCsurvey.mesh.dim ==3:
-                zsrc = Tx[2]
-
-        elif stype == 'dpdp':
-            MA = np.abs(Tx[0][0] - Rx[0][:,0])
-            MB = np.abs(Tx[1][0] - Rx[0][:,0])
-            NA = np.abs(Tx[0][0] - Rx[1][:,0])
-            NB = np.abs(Tx[1][0] - Rx[1][:,0])
-
-            # Create mid-point location
-            Cmid = (Tx[0][0] + Tx[1][0])/2
-            Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
-            if DCsurvey.mesh.dim == 2:
-                zsrc = (Tx[0][1] + Tx[1][1])/2
-            elif DCsurvey.mesh.dim ==3:
-                zsrc = (Tx[0][2] + Tx[1][2])/2
-
-        # Change output for dtype
-        if dtype == 'volt':
-
-            rho = np.hstack([rho,data])
-
-        else:
-
-            # 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 )
-                LEG.append(1./(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 ])
-        if DCsurvey.mesh.dim==3:
-            midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
-        elif DCsurvey.mesh.dim==2:
-            midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
-    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')
-
-    if clim == None:
-        vmin, vmax = rho.min(), rho.max()
-    else:
-        vmin, vmax = clim[0], clim[1]
-
-    grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
-    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax), vmin=vmin, vmax=vmax)
-    cbar = plt.colorbar(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 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)
-
-    # Plot apparent resistivity
-    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
-
-    #ax.set_xticklabels([])
-    #ax.set_yticklabels([])
-
-    plt.gca().set_aspect('equal', adjustable='box')
-
-
-
-    return ph, LEG
-
-def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
-    """
-        Load in endpoints and survey specifications to generate Tx, Rx location
-        stations.
-
-        Assumes flat topo for now...
-
-        Input:
-        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
-        :object mesh -> SimPEG mesh object
-        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
-        : param a, n -> pole seperation, number of rx dipoles per tx
-
-        Output:
-        :param Tx, Rx -> List objects for each tx location
-            Lines: P1x, P1y, P1z, P2x, P2y, P2z
-
-        Created on Wed December 9th, 2015
-
-        @author: dominiquef
-        !! Require clean up to deal with DCsurvey
-    """
-
-    from SimPEG import np
-
-    def xy_2_r(x1,x2,y1,y2):
-        r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) )
-        return r
-
-    ## Evenly distribute electrodes and put on surface
-    # Mesure survey length and direction
-    dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
-
-    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
-    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
-
-    nstn = np.floor( dl_len / a )
-
-    # Compute discrete pole location along line
-    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
-    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
-
-    if mesh.dim==2:
-        ztop = mesh.vectorNy[-1]
-        # Create line of P1 locations
-        M = np.c_[stn_x, np.ones(nstn).T*ztop]
-        # Create line of P2 locations
-        N = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
-
-    elif mesh.dim==3:
-        ztop = mesh.vectorNz[-1]
-        # Create line of P1 locations
-        M = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop]
-        # Create line of P2 locations
-        N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop]
-
-
-    ## Build list of Tx-Rx locations depending on survey type
-    # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
-    # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B]
-    SrcList = []
-
-
-    if stype != 'gradient':
-
-        for ii in range(0, int(nstn)-1):
-
-
-            if stype == 'dpdp':
-                tx = np.c_[M[ii,:],N[ii,:]]
-            elif stype == 'pdp':
-                tx = np.c_[M[ii,:],M[ii,:]]
-
-            # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
-
-            # Current elctrode seperation
-            AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
-
-            # Number of receivers to fit
-            nstn = np.min([np.floor( (AB - b) / a ) , n])
-
-            # Check if there is enough space, else break the loop
-            if nstn <= 0:
-                continue
-
-            # Compute discrete pole location along line
-            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
-            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
-
-            # Create receiver poles
-
-            if mesh.dim==3:
-                # Create line of P1 locations
-                P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop]
-                # Create line of P2 locations
-                P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop]
-                rxClass = DC.Rx.Dipole(P1, P2)
-
-            elif mesh.dim==2:
-                # Create line of P1 locations
-                P1 = np.c_[stn_x, np.ones(nstn).T*ztop]
-                # Create line of P2 locations
-                P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
-                rxClass = DC.Rx.Dipole_ky(P1, P2)
-
-            if stype == 'dpdp':
-                    srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
-            elif stype == 'pdp':
-                    srcClass = DC.Src.Pole([rxClass], M[ii,:])
-            SrcList.append(srcClass)
-
-    elif stype == 'gradient':
-
-        # Gradient survey only requires Tx at end of line and creates a square
-        # grid of receivers at in the middle at a pre-set minimum distance
-
-        # Get the edge limit of survey area
-        min_x = endl[0,0] + dl_x * b
-        min_y = endl[0,1] + dl_y * b
-
-        max_x = endl[1,0] - dl_x * b
-        max_y = endl[1,1] - dl_y * b
-
-        box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
-        box_w = box_l/2.
-
-        nstn = np.floor( box_l / a )
-
-        # Compute discrete pole location along line
-        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
-        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
-
-        # Define number of cross lines
-        nlin = int(np.floor( box_w / a ))
-        lind = range(-nlin,nlin+1)
-
-        ngrad = nstn * len(lind)
-
-        rx = np.zeros([ngrad,6])
-        for ii in range( len(lind) ):
-
-            # Move line in perpendicular direction by dipole spacing
-            lxx = stn_x - lind[ii]*a*dl_y
-            lyy = stn_y + lind[ii]*a*dl_x
-
-
-            M = np.c_[ lxx, lyy , np.ones(nstn).T*ztop]
-            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*ztop]
-            rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
-
-            if mesh.dim==3:
-                rxClass = DC.Rx.Dipole(rx[:,:3], rx[:,3:])
-            elif mesh.dim==2:
-                M = M[:,[0,2]]
-                N = N[:,[0,2]]
-                rxClass = DC.Rx.Dipole_ky(rx[:,[0,2]], rx[:,[3,5]])
-            srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
-        SrcList.append(srcClass)
-    else:
-        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
-
-
-    return SrcList
-

SimPEG/EM/Static/Utils/__init__.py

@@ -1 +0,0 @@
-from StaticUtils import *

SimPEG/EM/Static/__init__.py

@@ -1,3 +0,0 @@
-import DC
-import IP
-import SIP

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -87,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
         if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
+            print ">> Step waveform: Non-zero initial condition"        
             if mesh._meshType is 'CYL':
                 if mesh.isSymmetric:
                     MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -96,8 +96,8 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
             elif mesh._meshType is 'TENSOR':
                 MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
             else:
-                raise Exception('Unknown mesh for VMD')
-            return {"b": mesh.edgeCurl*MVP}
+                raise Exception('Unknown mesh for VMD') 
+            return {"b": mesh.edgeCurl*MVP}               
         elif self.waveformType == "GENERAL":
             print ">> General waveform: Zero initial condition"
             return {"b": np.zeros(mesh.nF)}
@@ -113,7 +113,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
-            raise Exception('Unknown mesh for VMD')
+            raise Exception('Unknown mesh for VMD')            
         return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 
@@ -122,7 +122,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         self.loc = loc
         self.radius = radius
         self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
+        SrcTDEM.__init__(self,rxList)        
 
     def getInitialFields(self, mesh):
         """Circular Loop, magnetic vector potential"""
@@ -153,7 +153,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
         else:
-            raise Exception('Unknown mesh for CircularLoop')
+            raise Exception('Unknown mesh for CircularLoop')            
         return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 

SimPEG/EM/__init__.py

@@ -1,6 +1,5 @@
 import TDEM
 import FDEM
-import Static
 import Base
 import Analytics
 import Utils

SimPEG/Examples/Inversion_IRLS.py

@@ -86,12 +86,12 @@ def run(N=200, plotIt=True):
     #reg.recModel = mrec
     reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 2e-3
-    reg.eps_q = 2e-3
+    reg.eps_p = 5e-2
+    reg.eps_q = 1e-2
     reg.norms   = [0., 0., 2., 2.]
     reg.wght = wr
 
-    opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3)
+    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()

SimPEG/MT/Utils/ediFilesUtils.py

@@ -7,17 +7,16 @@ 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)
 
@@ -26,8 +25,8 @@ class EDIimporter:
     comps = None
 
     # Hidden properties
-    _outEPSG = None
-    _2out = None
+    _outEPSG = None # Project info
+    _2out = None    # The projection operator
 
 
     def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -113,6 +112,12 @@ 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/Mesh/DiffOperators.py

@@ -584,67 +584,7 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
-    def getBCProjWF_simple(self, discretization='CC'):
-        """
 
-        The weak form boundary condition projection matrices
-        when mixed boundary condition is used
-
-
-        """
-
-        if discretization is not 'CC':
-            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
-
-        def projBC(n):
-            ij = ([0,n], [0,1])
-            vals = [0,0]
-            vals[0] = 1
-            vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(n+1,2))
-
-        def projDirichlet(n, bc):
-            bc = checkBC(bc)
-            ij = ([0,n], [0,1])
-            vals = [0,0]
-            if(bc[0] == 'dirichlet'):
-                vals[0] = -1
-            if(bc[1] == 'dirichlet'):
-                vals[1] = 1
-            return sp.csr_matrix((vals, ij), shape=(n+1,2))
-
-        BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
-        n = self.vnC
-        indF = self.faceBoundaryInd
-        if(self.dim == 1):
-            Pbc = projDirichlet(n[0], BC[0])
-            B = projBC(n[0])
-            indF = indF[0] | indF[1]
-            Pbc = Pbc*sdiag(self.area[indF])
-
-        elif(self.dim == 2):
-            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
-            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
-            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
-            B1 = sp.kron(speye(n[1]), projBC(n[0]))
-            B2 = sp.kron(projBC(n[1]), speye(n[0]))
-            B = sp.block_diag((B1, B2), format="csr")
-            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
-            Pbc = Pbc*sdiag(self.area[indF])
-
-        elif(self.dim == 3):
-            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
-            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
-            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
-            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
-            B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0]))
-            B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0]))
-            B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0]))
-            B = sp.block_diag((B1, B2, B3), format="csr")
-            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
-            Pbc = Pbc*sdiag(self.area[indF])
-
-        return Pbc, B.T
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Mesh/MeshIO.py

@@ -24,7 +24,6 @@ class TensorMeshIO(object):
                     re = int(sp[0])*(' ' + 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/View.py

@@ -218,7 +218,7 @@ class TensorView(object):
             return out
         viewOpts = ['real','imag','abs','vec']
         normalOpts = ['X', 'Y', 'Z']
-        vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
+        vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
 
         # Some user error checking
         assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)

SimPEG/Regularization.py

@@ -663,9 +663,8 @@ class Sparse(Simple):
             self.Rs = Utils.speye(self.regmesh.nC)
 
         else:
-            f_m = self.curModel - self.reg.mref
+            f_m = self.mapping * (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 )
 
         return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
@@ -679,7 +678,7 @@ class Sparse(Simple):
             self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
             self.rx = self.R( f_m , self.eps_q, self.norms[1])
             self.Rx = Utils.sdiag( self.rx )
 
@@ -693,7 +692,7 @@ class Sparse(Simple):
             self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
             self.ry = self.R( f_m , self.eps_q, self.norms[2])
             self.Ry = Utils.sdiag( self.ry )
 
@@ -707,7 +706,7 @@ class Sparse(Simple):
             self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
             self.rz = self.R( f_m , self.eps_q, self.norms[3])
             self.Rz = Utils.sdiag( self.rz )
 
@@ -728,14 +727,15 @@ class Sparse(Simple):
     @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 sp.vstack(wlist)
+        if getattr(self, '_W', None) is None:
+    
+            wlist = (self.Wsmall, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
 
     def R(self, f_m , eps, exponent):
 
-        eta = (eps**(1-exponent/2.))**0.5
-        r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
+        eta = (eps**(1.-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
 
         return r

SimPEG/Utils/io_utils.py

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

tests/em/static/__init__.py

@@ -1,12 +0,0 @@
-import os
-import glob
-import unittest
-
-if __name__ == '__main__':
-    test_file_strings = glob.glob('test_*.py')
-    module_strings = [str[0:len(str)-3] for str in test_file_strings]
-    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
-              in module_strings]
-    testSuite = unittest.TestSuite(suites)
-
-    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/static/test_DC_2D_analytic.py

@@ -1,69 +0,0 @@
-import unittest
-from SimPEG import Mesh, Utils, EM, Maps, np
-import SimPEG.EM.Static.DC as DC
-
-class DCProblemAnalyticTests(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 12.5
-        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
-        hy = [(cs,7, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
-        sighalf = 1e-2
-        sigma = np.ones(mesh.nC)*sighalf
-        x = np.linspace(-135, 250., 20)
-        M = Utils.ndgrid(x-12.5, np.r_[0.])
-        N = Utils.ndgrid(x+12.5, np.r_[0.])
-        A0loc = np.r_[-150, 0.]
-        A1loc = np.r_[-130, 0.]
-        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
-        data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, flag="halfspace")
-
-        rx = DC.Rx.Dipole(M, N)
-        src0 = DC.Src.Pole([rx], A0loc)
-        survey = DC.Survey_ky([src0])
-
-        self.survey = survey
-        self.mesh = mesh
-        self.sigma = sigma
-        self.data_anal = data_anal
-
-        try:
-            from pymatsolver import MumpsSolver
-            self.Solver = MumpsSolver
-        except ImportError, e:
-            self.Solver = SolverLU
-
-    def test_Problem3D_N(self):
-
-        problem = DC.Problem2D_N(self.mesh)
-        problem.Solver = self.Solver
-        problem.pair(self.survey)
-        data = self.survey.dpred(self.sigma)
-        err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
-        if err < 0.05:
-            passed = True
-            print ">> DC analytic test for Problem3D_N is passed"
-        else:
-            passed = False
-            print ">> DC analytic test for Problem3D_N is failed"
-        self.assertTrue(passed)
-
-    def test_Problem3D_CC(self):
-        problem = DC.Problem2D_CC(self.mesh)
-        problem.Solver = self.Solver
-        problem.pair(self.survey)
-        data = self.survey.dpred(self.sigma)
-        err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
-        if err < 0.05:
-            passed = True
-            print ">> DC analytic test for Problem3D_CC is passed"
-        else:
-            passed = False
-            print ">> DC analytic test for Problem3D_CC is failed"
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()
-

tests/em/static/test_DC_2D_jvecjtvecadj.py

@@ -1,127 +0,0 @@
-import unittest
-from SimPEG import *
-import SimPEG.EM.Static.DC as DC
-
-
-class DCProblem_2DTestsCC(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 12.5
-        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
-        hy = [(cs,7, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
-        x = np.linspace(-135, 250., 20)
-        M = Utils.ndgrid(x-12.5, np.r_[0.])
-        N = Utils.ndgrid(x+12.5, np.r_[0.])
-        A0loc = np.r_[-150, 0.]
-        A1loc = np.r_[-130, 0.]
-        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
-        rx = DC.Rx.Dipole_ky(M, N)
-        src0 = DC.Src.Pole([rx], A0loc)
-        src1 = DC.Src.Pole([rx], A1loc)
-        survey = DC.Survey_ky([src0, src1])
-        problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
-        problem.pair(survey)
-
-        mSynth = np.ones(mesh.nC)*1.
-        survey.makeSyntheticData(mSynth)
-
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv    = inv
-        self.reg    = reg
-        self.p      = problem
-        self.mesh   = mesh
-        self.m0     = mSynth
-        self.survey = survey
-        self.dmis   = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-10
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-class DCProblemTestsN(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 12.5
-        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
-        hy = [(cs,7, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
-        x = np.linspace(-135, 250., 20)
-        M = Utils.ndgrid(x-12.5, np.r_[0.])
-        N = Utils.ndgrid(x+12.5, np.r_[0.])
-        A0loc = np.r_[-150, 0.]
-        A1loc = np.r_[-130, 0.]
-        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
-        rx = DC.Rx.Dipole_ky(M, N)
-        src0 = DC.Src.Pole([rx], A0loc)
-        src1 = DC.Src.Pole([rx], A1loc)
-        survey = DC.Survey_ky([src0, src1])
-        problem = DC.Problem2D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
-        problem.pair(survey)
-
-        mSynth = np.ones(mesh.nC)*1.
-        survey.makeSyntheticData(mSynth)
-
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-8
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()

tests/em/static/test_DC_analytic.py

@@ -1,71 +0,0 @@
-import unittest
-from SimPEG import Mesh, Utils, EM, Maps, np
-import SimPEG.EM.Static.DC as DC
-
-class DCProblemAnalyticTests(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 25.
-        hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-        hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
-        hz = [(cs,7, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-        sigma = np.ones(mesh.nC)*1e-2
-
-        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
-        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
-        Aloc = np.r_[-200., 0., 0.]
-        Bloc = np.r_[200., 0., 0.]
-        M = Utils.ndgrid(x-25.,y, np.r_[0.])
-        N = Utils.ndgrid(x+25.,y, np.r_[0.])
-        phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace")
-        phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace")
-        data_anal = phiA-phiB
-
-        rx = DC.Rx.Dipole(M, N)
-        src = DC.Src.Dipole([rx], Aloc, Bloc)
-        survey = DC.Survey([src])
-
-        self.survey = survey
-        self.mesh = mesh
-        self.sigma = sigma
-        self.data_anal = data_anal
-
-        try:
-            from pymatsolver import MumpsSolver
-            self.Solver = MumpsSolver
-        except ImportError, e:
-            self.Solver = SolverLU
-
-    def test_Problem3D_N(self):
-        problem = DC.Problem3D_N(self.mesh)
-        problem.Solver = self.Solver
-        problem.pair(self.survey)
-        data = self.survey.dpred(self.sigma)
-        err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
-        if err < 0.2:
-            passed = True
-            print ">> DC analytic test for Problem3D_N is passed"
-        else:
-            passed = False
-            print ">> DC analytic test for Problem3D_N is failed"
-        self.assertTrue(passed)
-
-    def test_Problem3D_CC(self):
-        problem = DC.Problem3D_CC(self.mesh)
-        problem.Solver = self.Solver
-        problem.pair(self.survey)
-        data = self.survey.dpred(self.sigma)
-        err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
-        if err < 0.2:
-            passed = True
-            print ">> DC analytic test for Problem3D_CC is passed"
-        else:
-            passed = False
-            print ">> DC analytic test for Problem3D_CC is failed"
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()
-

tests/em/static/test_DC_jvecjtvecadj.py

@@ -1,127 +0,0 @@
-import unittest
-from SimPEG import *
-import SimPEG.EM.Static.DC as DC
-
-
-class DCProblemTestsCC(unittest.TestCase):
-
-    def setUp(self):
-
-        aSpacing=2.5
-        nElecs=5
-
-        surveySize = nElecs*aSpacing - aSpacing
-        cs = surveySize/nElecs/4
-
-        mesh = Mesh.TensorMesh([
-                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
-                [(cs,3, -1.3),(cs,3,1.3)],
-                # [(cs,5, -1.3),(cs,10)]
-            ],'CN')
-
-        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
-        survey = DC.Survey(srcList)
-        problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
-        problem.pair(survey)
-
-        mSynth = np.ones(mesh.nC)
-        survey.makeSyntheticData(mSynth)
-
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-10
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-class DCProblemTestsN(unittest.TestCase):
-
-    def setUp(self):
-
-        aSpacing=2.5
-        nElecs=10
-
-        surveySize = nElecs*aSpacing - aSpacing
-        cs = surveySize/nElecs/4
-
-        mesh = Mesh.TensorMesh([
-                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
-                [(cs,3, -1.3),(cs,3,1.3)],
-        #         [(cs,5, -1.3),(cs,10)]
-            ],'CN')
-
-        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
-        survey = DC.Survey(srcList)
-        problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
-        problem.pair(survey)
-
-        mSynth = np.ones(mesh.nC)
-        survey.makeSyntheticData(mSynth)
-
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-8
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()

tests/em/static/test_IP_fwd.py

@@ -1,96 +0,0 @@
-import unittest
-from SimPEG import Mesh, Utils, EM, Maps, np
-import SimPEG.EM.Static.DC as DC
-import SimPEG.EM.Static.IP as IP
-
-class IPProblemAnalyticTests(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 12.5
-        npad=2
-        hx = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
-        hy = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
-        hz = [(cs,npad, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-
-        x = mesh.vectorCCx[(mesh.vectorCCx>-80.)&(mesh.vectorCCx<80.)]
-        y = mesh.vectorCCx[(mesh.vectorCCy>-80.)&(mesh.vectorCCy<80.)]
-        Aloc = np.r_[-100., 0., 0.]
-        Bloc = np.r_[100., 0., 0.]
-        M = Utils.ndgrid(x-12.5,y, np.r_[0.])
-        N = Utils.ndgrid(x+12.5,y, np.r_[0.])
-        radius = 50.
-        xc = np.r_[0., 0., -100]
-        blkind = Utils.ModelBuilder.getIndicesSphere(xc, radius, mesh.gridCC)
-        sigmaInf = np.ones(mesh.nC)*1e-2
-        eta = np.zeros(mesh.nC)
-        eta[blkind] = 0.1
-        sigma0 = sigmaInf*(1.-eta)
-
-        rx = DC.Rx.Dipole(M, N)
-        src = DC.Src.Dipole([rx], Aloc, Bloc)
-        surveyDC = DC.Survey([src])
-
-        self.surveyDC = surveyDC
-        self.mesh = mesh
-        self.sigmaInf = sigmaInf
-        self.sigma0 = sigma0
-        self.src = src
-        self.eta = eta
-
-        try:
-            from pymatsolver import MumpsSolver
-            self.Solver = MumpsSolver
-        except ImportError, e:
-            self.Solver = SolverLU
-
-    def test_Problem3D_N(self):
-
-        problemDC = DC.Problem3D_N(self.mesh)
-        problemDC.Solver = self.Solver
-        problemDC.pair(self.surveyDC)
-        data0 = self.surveyDC.dpred(self.sigma0)
-        finf = problemDC.fields(self.sigmaInf)
-        datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
-        problemIP = IP.Problem3D_N(self.mesh, sigma=self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
-        problemIP.Solver = self.Solver
-        surveyIP = IP.Survey([self.src])
-        problemIP.pair(surveyIP)
-        data_full = data0 - datainf
-        data = surveyIP.dpred(self.eta)
-        err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
-        if err < 0.05:
-            passed = True
-            print ">> IP forward test for Problem3D_N is passed"
-        else:
-            passed = False
-            print ">> IP forward test for Problem3D_N is failed"
-        self.assertTrue(passed)
-
-    def test_Problem3D_CC(self):
-
-        problemDC = DC.Problem3D_CC(self.mesh)
-        problemDC.Solver = self.Solver
-        problemDC.pair(self.surveyDC)
-        data0 = self.surveyDC.dpred(self.sigma0)
-        finf = problemDC.fields(self.sigmaInf)
-        datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
-        problemIP = IP.Problem3D_CC(self.mesh, rho=1./self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
-        problemIP.Solver = self.Solver
-        surveyIP = IP.Survey([self.src])
-        problemIP.pair(surveyIP)
-        data_full = data0 - datainf
-        data = surveyIP.dpred(self.eta)
-        err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
-        if err < 0.05:
-            passed = True
-            print ">> IP forward test for Problem3D_CC is passed"
-        else:
-            passed = False
-            print ">> IP forward test for Problem3D_CC is failed"
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()
-

tests/em/static/test_IP_jvecjtvecadj.py

@@ -1,126 +0,0 @@
-import unittest
-from SimPEG import *
-import SimPEG.EM.Static.DC as DC
-import SimPEG.EM.Static.IP as IP
-
-
-class IPProblemTestsCC(unittest.TestCase):
-
-    def setUp(self):
-
-        aSpacing=2.5
-        nElecs=5
-
-        surveySize = nElecs*aSpacing - aSpacing
-        cs = surveySize/nElecs/4
-
-        mesh = Mesh.TensorMesh([
-                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
-                [(cs,3, -1.3),(cs,3,1.3)],
-                # [(cs,5, -1.3),(cs,10)]
-            ],'CN')
-
-        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
-        survey = IP.Survey(srcList)
-        sigma = np.ones(mesh.nC)
-        problem = IP.Problem3D_CC(mesh, rho=1./sigma)
-        problem.pair(survey)
-        mSynth = np.ones(mesh.nC)*0.1
-        survey.makeSyntheticData(mSynth)
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-10
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-class IPProblemTestsN(unittest.TestCase):
-
-    def setUp(self):
-
-        aSpacing=2.5
-        nElecs=5
-
-        surveySize = nElecs*aSpacing - aSpacing
-        cs = surveySize/nElecs/4
-
-        mesh = Mesh.TensorMesh([
-                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
-                [(cs,3, -1.3),(cs,3,1.3)],
-                # [(cs,5, -1.3),(cs,10)]
-            ],'CN')
-
-        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
-        survey = IP.Survey(srcList)
-        sigma = np.ones(mesh.nC)
-        problem = IP.Problem3D_N(mesh, sigma=sigma)
-        problem.pair(survey)
-        mSynth = np.ones(mesh.nC)*0.1
-        survey.makeSyntheticData(mSynth)
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-8
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()

tests/em/static/test_SIP_jvecjtvecadj.py

@@ -1,232 +0,0 @@
-import unittest
-from SimPEG import *
-import SimPEG
-from SimPEG import Mesh, Utils, EM, Maps, np, Survey
-from SimPEG.EM.Static import SIP, DC, IP
-from pymatsolver import MumpsSolver
-
-
-class IPProblemTestsCC(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 25.
-        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hz = [(cs,0, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
-        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
-        sigma = np.ones(mesh.nC)*1e-2
-        eta = np.zeros(mesh.nC)
-        tau = np.ones_like(sigma)*1.
-        eta[blkind0] = 0.1
-        eta[blkind1] = 0.1
-        tau[blkind0] = 0.1
-        tau[blkind1] = 0.01
-
-        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
-        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
-        Aloc = np.r_[-200., 0., 0.]
-        Bloc = np.r_[200., 0., 0.]
-        M = Utils.ndgrid(x-25.,y, np.r_[0.])
-        N = Utils.ndgrid(x+25.,y, np.r_[0.])
-
-        times = np.arange(10)*1e-3 + 1e-3
-        rx = SIP.Rx.Dipole(M, N, times)
-        src = SIP.Src.Dipole([rx], Aloc, Bloc)
-        survey = SIP.Survey([src])
-        colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
-        problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap)
-        problem.Solver = MumpsSolver
-        problem.pair(survey)
-        mSynth = np.r_[eta, 1./tau]
-        survey.makeSyntheticData(mSynth)
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC*2)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-10
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-class IPProblemTestsN(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 25.
-        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hz = [(cs,0, -1.3),(cs,20)]
-        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
-        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
-        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
-        sigma = np.ones(mesh.nC)*1e-2
-        eta = np.zeros(mesh.nC)
-        tau = np.ones_like(sigma)*1.
-        eta[blkind0] = 0.1
-        eta[blkind1] = 0.1
-        tau[blkind0] = 0.1
-        tau[blkind1] = 0.01
-
-        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
-        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
-        Aloc = np.r_[-200., 0., 0.]
-        Bloc = np.r_[200., 0., 0.]
-        M = Utils.ndgrid(x-25.,y, np.r_[0.])
-        N = Utils.ndgrid(x+25.,y, np.r_[0.])
-
-        times = np.arange(10)*1e-3 + 1e-3
-        rx = SIP.Rx.Dipole(M, N, times)
-        src = SIP.Src.Dipole([rx], Aloc, Bloc)
-        survey = SIP.Survey([src])
-        colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
-        problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
-        problem.Solver = MumpsSolver
-        problem.pair(survey)
-        mSynth = np.r_[eta, 1./tau]
-        survey.makeSyntheticData(mSynth)
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        reg = Regularization.Tikhonov(mesh)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC*2)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-8
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-class IPProblemTestsN_air(unittest.TestCase):
-
-    def setUp(self):
-
-        cs = 25.
-        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
-        hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)]
-        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC")
-        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
-        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
-        sigma = np.ones(mesh.nC)*1e-2
-        airind = mesh.gridCC[:,2]>0.
-        sigma[airind] = 1e-8
-        eta = np.zeros(mesh.nC)
-        tau = np.ones_like(sigma)*1.
-        eta[blkind0] = 0.1
-        eta[blkind1] = 0.1
-        tau[blkind0] = 0.1
-        tau[blkind1] = 0.01
-
-        actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.)
-        actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.)
-
-        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
-        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
-        Aloc = np.r_[-200., 0., 0.]
-        Bloc = np.r_[200., 0., 0.]
-        M = Utils.ndgrid(x-25.,y, np.r_[0.])
-        N = Utils.ndgrid(x+25.,y, np.r_[0.])
-
-        times = np.arange(10)*1e-3 + 1e-3
-        rx = SIP.Rx.Dipole(M, N, times)
-        src = SIP.Src.Dipole([rx], Aloc, Bloc)
-        survey = SIP.Survey([src])
-        colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)]
-        problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
-        problem.Solver = MumpsSolver
-        problem.pair(survey)
-        mSynth = np.r_[eta[~airind], 1./tau[~airind]]
-        survey.makeSyntheticData(mSynth)
-        # Now set up the problem to do some minimization
-        dmis = DataMisfit.l2_DataMisfit(survey)
-        regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2))
-        reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind)
-        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
-        inv = Inversion.BaseInversion(invProb)
-
-        self.inv = inv
-        self.reg = reg
-        self.p =     problem
-        self.mesh = mesh
-        self.m0 = mSynth
-        self.survey = survey
-        self.dmis = dmis
-
-    def test_misfit(self):
-        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-    def test_adjoint(self):
-        # Adjoint Test
-        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
-        v = np.random.rand(self.mesh.nC)
-        w = np.random.rand(self.survey.dobs.shape[0])
-        wtJv = w.dot(self.p.Jvec(self.m0, v))
-        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
-        passed = np.abs(wtJv - vtJtw) < 1e-8
-        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
-        self.assertTrue(passed)
-
-    def test_dataObj(self):
-        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
-        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
-        self.assertTrue(passed)
-
-if __name__ == '__main__':
-    unittest.main()

tests/mesh/test_Mixed_boundaryPoisson.py

@@ -1,411 +0,0 @@
-import numpy as np
-import scipy.sparse as sp
-import unittest
-import matplotlib.pyplot as plt
-from SimPEG import *
-
-MESHTYPES = ['uniformTensorMesh']
-
-def getxBCyBC_CC(mesh, alpha, beta, gamma):
-# def getxBCyBC(mesh, alpha, beta, gamma):
-    """
-    This is a subfunction generating mixed-boundary condition:
-
-    .. math::
-
-        \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
-
-        \rho \vec{j} = -\nabla \phi \phi
-
-        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
-
-        xBC = f_1(\alpha, \beta, \gamma)
-        yBC = f(\alpha, \beta, \gamma)
-
-    Computes xBC and yBC for cell-centered discretizations
-    """
-    if mesh.dim == 1: #1D
-        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
-            raise Exception("Lenght of list, alpha should be 2")
-        fCCxm,fCCxp = mesh.cellBoundaryInd
-        nBC = fCCxm.sum()+fCCxp.sum()
-        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
-        h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-
-        xBC = np.r_[xBC_xm, xBC_xp]
-        yBC = np.r_[yBC_xm, yBC_xp]
-
-    elif mesh.dim == 2: #2D
-        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
-            raise Exception("Lenght of list, alpha should be 4")
-
-        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
-        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
-        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
-        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
-
-        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
-        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
-        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
-        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-        xBC_ym = 0.5*a_ym
-        xBC_yp = 0.5*a_yp/b_yp
-        yBC_ym = 0.5*(1.-b_ym)
-        yBC_yp = 0.5*(1.-1./b_yp)
-
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
-
-        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
-        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
-        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
-        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
-
-        xBC = np.r_[xBC_x, xBC_y]
-        yBC = np.r_[yBC_x, yBC_y]
-
-    elif mesh.dim == 3: #3D
-        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
-            raise Exception("Lenght of list, alpha should be 6")
-        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
-        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
-        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
-
-        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
-        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
-        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
-        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
-        alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4]
-        alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]
-
-        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
-        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
-        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
-
-        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
-        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
-        h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
-
-        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
-        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
-        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
-        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
-
-        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
-        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
-        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
-        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
-
-        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
-        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
-        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
-        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
-
-        xBC_xm = 0.5*a_xm
-        xBC_xp = 0.5*a_xp/b_xp
-        yBC_xm = 0.5*(1.-b_xm)
-        yBC_xp = 0.5*(1.-1./b_xp)
-        xBC_ym = 0.5*a_ym
-        xBC_yp = 0.5*a_yp/b_yp
-        yBC_ym = 0.5*(1.-b_ym)
-        yBC_yp = 0.5*(1.-1./b_yp)
-        xBC_zm = 0.5*a_zm
-        xBC_zp = 0.5*a_zp/b_zp
-        yBC_zm = 0.5*(1.-b_zm)
-        yBC_zp = 0.5*(1.-1./b_zp)
-
-        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
-        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
-        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
-
-        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
-        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
-        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
-
-        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
-        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
-        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
-
-        xBC = np.r_[xBC_x, xBC_y, xBC_z]
-        yBC = np.r_[yBC_x, yBC_y, yBC_z]
-
-    return xBC, yBC
-
-class Test1D_InhomogeneousMixed(Tests.OrderTest):
-    name = "1D - Mixed"
-    meshTypes = MESHTYPES
-    meshDimension = 1
-    expectedOrders = 2
-    meshSizes = [4, 8, 16, 32]
-
-    def getError(self):
-        #Test function
-        phi_fun = lambda x: np.cos(np.pi*x)
-        j_fun = lambda x: np.pi*np.sin(np.pi*x)
-        phi_deriv = lambda x: -j_fun(x)
-        q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
-
-        xc_ana = phi_fun(self.M.gridCC)
-        q_ana = q_fun(self.M.gridCC)
-        j_ana = j_fun(self.M.gridFx)
-
-        # Get boundary locations
-        vecN = self.M.vectorNx
-        vecC = self.M.vectorCCx
-
-        # Setup Mixed B.C (alpha, beta, gamma)
-        alpha_xm, alpha_xp = 1., 1.
-        beta_xm, beta_xp =  1., 1.
-        alpha = np.r_[alpha_xm, alpha_xp]
-        beta = np.r_[beta_xm, beta_xp]
-        vecN = self.M.vectorNx
-        vecC = self.M.vectorCCx
-        phi_bc = phi_fun(vecN[[0,-1]])
-        phi_deriv_bc = phi_deriv(vecN[[0,-1]])
-        gamma = alpha*phi_bc + beta*phi_deriv_bc
-        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
-
-
-        sigma = np.ones(self.M.nC)
-        Mfrho = self.M.getFaceInnerProduct(1./sigma)
-        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
-        V = Utils.sdiag(self.M.vol)
-        Div = V*self.M.faceDiv
-        P_BC, B = self.M.getBCProjWF_simple()
-        q = q_fun(self.M.gridCC)
-        M = B*self.M.aveCC2F
-        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
-        # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC
-        rhs = V*q + Div*MfrhoI*P_BC*x_BC
-        A = Div*MfrhoI*G
-
-        if self.myTest == 'xc':
-            #TODO: fix the null space
-            Ainv = Solver(A)
-            xc = Ainv*rhs
-            err = np.linalg.norm((xc-xc_ana), np.inf)
-        else:
-            NotImplementedError
-        return err
-
-
-    def test_order(self):
-        print "==== Testing Mixed boudary conduction for CC-problem ===="
-        self.name = "1D"
-        self.myTest = 'xc'
-        self.orderTest()
-
-class Test2D_InhomogeneousMixed(Tests.OrderTest):
-    name = "2D - Mixed"
-    meshTypes = MESHTYPES
-    meshDimension = 2
-    expectedOrders = 2
-    meshSizes = [4, 8, 16, 32]
-
-    def getError(self):
-        #Test function
-        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
-        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
-        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
-        phideriv_funX = lambda x: -j_funX(x)
-        phideriv_funY = lambda x: -j_funY(x)
-        q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
-
-        xc_ana = phi_fun(self.M.gridCC)
-        q_ana = q_fun(self.M.gridCC)
-        jX_ana = j_funX(self.M.gridFx)
-        jY_ana = j_funY(self.M.gridFy)
-        j_ana = np.r_[jX_ana,jY_ana]
-
-        # Get boundary locations
-        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
-        gBFxm = self.M.gridFx[fxm,:]
-        gBFxp = self.M.gridFx[fxp,:]
-        gBFym = self.M.gridFy[fym,:]
-        gBFyp = self.M.gridFy[fyp,:]
-
-        # Setup Mixed B.C (alpha, beta, gamma)
-        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-
-        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
-        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
-
-        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
-        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
-
-        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
-        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
-        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
-        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
-        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
-
-        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
-        beta = [beta_xm, beta_xp, beta_ym, beta_yp]
-        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
-
-
-        sigma = np.ones(self.M.nC)
-        Mfrho = self.M.getFaceInnerProduct(1./sigma)
-        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
-        V = Utils.sdiag(self.M.vol)
-        Div = V*self.M.faceDiv
-        P_BC, B = self.M.getBCProjWF_simple()
-        q = q_fun(self.M.gridCC)
-        M = B*self.M.aveCC2F
-        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
-        rhs = V*q + Div*MfrhoI*P_BC*x_BC
-        A = Div*MfrhoI*G
-
-        if self.myTest == 'xc':
-            Ainv = Solver(A)
-            xc = Ainv*rhs
-            err = np.linalg.norm((xc-xc_ana), np.inf)
-        else:
-            NotImplementedError
-        return err
-
-
-    def test_order(self):
-        print "==== Testing Mixed boudary conduction for CC-problem ===="
-        self.name = "2D"
-        self.myTest = 'xc'
-        self.orderTest()
-
-class Test3D_InhomogeneousMixed(Tests.OrderTest):
-    name = "3D - Mixed"
-    meshTypes = MESHTYPES
-    meshDimension = 3
-    expectedOrders = 2
-    meshSizes = [4, 8, 16]
-
-    def getError(self):
-        #Test function
-        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
-        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
-        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
-        j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
-
-        phideriv_funX = lambda x: -j_funX(x)
-        phideriv_funY = lambda x: -j_funY(x)
-        phideriv_funZ = lambda x: -j_funZ(x)
-
-        q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
-
-        xc_ana = phi_fun(self.M.gridCC)
-        q_ana = q_fun(self.M.gridCC)
-        jX_ana = j_funX(self.M.gridFx)
-        jY_ana = j_funY(self.M.gridFy)
-        j_ana = np.r_[jX_ana,jY_ana,jY_ana]
-
-        # Get boundary locations
-        fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
-        gBFxm = self.M.gridFx[fxm,:]
-        gBFxp = self.M.gridFx[fxp,:]
-        gBFym = self.M.gridFy[fym,:]
-        gBFyp = self.M.gridFy[fyp,:]
-        gBFzm = self.M.gridFz[fzm,:]
-        gBFzp = self.M.gridFz[fzp,:]
-
-        # Setup Mixed B.C (alpha, beta, gamma)
-        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
-        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
-        alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-        beta_zm, beta_zp =  np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
-
-
-        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
-        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
-        phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
-
-        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
-        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
-        phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp)
-
-        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
-        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
-        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
-        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
-        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
-        gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm)
-        gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp)
-
-        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
-        beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
-        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
-
-        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
-
-
-        sigma = np.ones(self.M.nC)
-        Mfrho = self.M.getFaceInnerProduct(1./sigma)
-        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
-        V = Utils.sdiag(self.M.vol)
-        Div = V*self.M.faceDiv
-        P_BC, B = self.M.getBCProjWF_simple()
-        q = q_fun(self.M.gridCC)
-        M = B*self.M.aveCC2F
-        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
-        rhs = V*q + Div*MfrhoI*P_BC*x_BC
-        A = Div*MfrhoI*G
-
-        if self.myTest == 'xc':
-            #TODO: fix the null space
-            Ainv = Solver(A)
-            xc = Ainv*rhs
-            err = np.linalg.norm((xc-xc_ana), np.inf)
-        else:
-            NotImplementedError
-        return err
-
-
-    def test_order(self):
-        print "==== Testing Mixed boudary conduction for CC-problem ===="
-        self.name = "3D"
-        self.myTest = 'xc'
-        self.orderTest()
-
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/mt/test_Problem3D_againstAnalytic.py

@@ -242,9 +242,6 @@ 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))