remove import of ipywidgets

# Conflicts:
#	SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
#	SimPEG/Examples/__init__.py
#	SimPEG/Examples/sphereElectrostatic_example.py
#	SimPEG/Optimization.py
#	docs/examples/DC_PseudoSection_Simulation.rst
#	docs/examples/Inversion_IRLS.rst
#	docs/examples/MT_1D_analytic_nlayer_Earth.rst

SimPEG/Directives.py

@@ -305,15 +305,28 @@ class Update_IRLS(InversionDirective):
         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.)
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
 
 

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -0,0 +1,427 @@
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+import numpy as np
+from SimPEG.EM.Utils import k, omega
+
+"""
+MT1D: n layered earth problem
+*****************************
+
+Author: Thibaut Astic
+Contact: thast@eos.ubc.ca
+
+
+This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
+
+We start by looking at Maxwell's equations in the electric
+field \\\(\\\mathbf{E}\\) and the magnetic flux
+\\\(\\\mathbf{H}\\) to write the wave equations
+\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
+\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
+
+Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
+\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
+\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
+
+With U and D the Up and Down components of the E-field.
+
+The iteration from one layer to another is ensure by:
+
+\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
+    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
+
+And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
+and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
+
+The layer 0 is assumed to be the air layer.
+
+"""
+
+#Define a frquency range for a survey
+frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
+
+#Functions to create random physical Properties for a n-layered earth
+thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
+                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
+                                                            ,decimals =1))
+
+sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
+                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
+                                                      ,decimals=3))
+
+mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
+                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
+                                                    ,decimals=1))
+
+eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
+                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
+                                                                ,decimals=1))
+
+#Evaluate Impedance Z of a layer
+ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
+
+#Complex Cole-Cole Conductivity - EM utils
+PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
+
+#Converted thickness array into top of layer array
+top = lambda thick: np.cumsum(thick)
+
+#Propagation Matrix and theirs inverses
+
+#matrix T for transition of Up and Down components accross a layer
+T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
+
+Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
+
+#transition of Up and Down components accross a layer
+UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
+
+
+#matrix P relating Up and Down components with E and H fields
+P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
+
+Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
+
+
+#Time Variation of E and H
+E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
+H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
+
+#Plot the configuration of the problem
+def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
+
+    topn = top(thick)
+    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
+
+    ax.set_ylim([z.min(),z.max()])
+    ax.set_xlim([-widthg,widthg])
+
+    ax.set_ylabel("Depth (m)", fontsize=16.)
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+    #define filling for the different layers
+    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
+
+    #Write the physical properties of air
+    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
+            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
+    for i in range(1,len(topn)-1,1):
+        if topn[i] == topn[i+1]:
+            pass
+        else:
+            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
+                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
+                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
+                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
+            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
+
+    #Write the physical properties of the n-th layer and fill it with pattern
+    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
+    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
+
+    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
+            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
+            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
+            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    #plot Trees!
+    ax.annotate("",
+            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+
+    ax.invert_yaxis()
+
+    return ax
+
+#Propagate Up and Down component for a certain frequency & evaluate E and H field
+
+def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
+
+    sigcm = np.zeros_like(sig,dtype='complex_')
+
+    for j in range(1,len(sig)):
+        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
+
+    K = k(f, sigcm, mu, eps)
+    Z = ImpZ(f,mu,K)
+
+    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+
+    UD[1,-1] = 1.
+
+    for i in range(-2,-(n+2),-1):
+
+        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
+        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
+
+    for j in range(0,n+1):
+        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
+
+    return UD, EH, Z ,K
+
+
+#Evaluate the apparent resistivity and phase for a frequency range
+def appres(F,H,sig,chg,taux,c,mu,eps,n):
+
+    Res = np.zeros_like(F)
+    Phase = np.zeros_like(F)
+    App_ImpZ= np.zeros_like(F,dtype='complex_')
+
+    for i in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
+
+        App_ImpZ[i] = EH[0,1]/EH[1,1]
+
+        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
+        Phase[i] = np.angle(App_ImpZ[i], deg = True)
+
+    return Res,Phase
+
+#Evaluate Up, Down components, E and H field, for a frequency range,
+#a discretized depth range and a time range (use to calculate envelope)
+def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
+
+    topc = top(H)
+
+    layer = np.zeros(len(zsample),dtype=np.int)-1
+
+    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
+
+    for i in range(0,n+1,1):
+        layer = layer+(zsample>=topc[i])*1
+
+    for j in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
+
+        for p in range(0,len(zsample)):
+
+            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
+
+            for q in range(0,len(tsample)):
+
+                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
+                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
+                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+
+    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
+
+
+#Function to Plot Apparent Resistivity and Phase
+def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
+
+    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
+
+    fig,ax = plt.subplots(1,2,figsize=(16,10))
+
+    ax[0].scatter(Res,F,color='black')
+    ax[0].set_xscale('Log')
+    ax[0].set_yscale('Log')
+    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
+    ax[0].set_ylim([F.min(),F.max()])
+    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
+    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
+    ax[0].grid(which='major')
+
+    ax0 = ax[0].twiny()
+
+    ax0.set_xlim([0.,90.])
+    ax0.set_ylim([F.min(),F.max()])
+    ax0.scatter(Phase,F,color='purple')
+    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
+
+    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
+
+    ax[0].tick_params(labelsize=16)
+    ax[1].tick_params(labelsize=16)
+    ax0.tick_params(labelsize=16)
+
+    if PlotEnvelope:
+
+        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
+        fenvelope1n=np.ones(100)*fenvelope
+        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
+
+        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
+        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
+
+        ax1=ax[1].twiny()
+
+        ax[1].tick_params(labelsize=16)
+        ax1.tick_params(labelsize=16)
+
+        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
+
+        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
+
+        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
+                      color='blue', alpha=0.1)
+
+        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
+                      color='red', alpha=0.1)
+
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+    else:
+        print 'No envelop (if True, might be slow)'
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
+        ax[1].get_xaxis().set_ticks([])
+
+    plt.show()
+
+#Interactive MT for Notebook
+def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
+
+    frangn=frange(-5,5,100.)
+    sig3= np.array([0.,0.001,0.1, 0.001])
+    thick3 = np.array([120000.,50.,50.])
+    eps3=np.array([1.,1.,1.,1])
+    mu3=np.array([1.,1.,1.,1])
+    chg3=np.array([0.,0.1,0.,0.2])
+    chg3_0=np.array([0.,0.1,0.,0.])
+    taux3=np.array([0.,0.1,0.,0.1])
+    c3=np.array([1.,1.,1.,1.])
+
+    sig3[1]=sigl1
+    sig3[1]=10.**sig3[1]
+    sig3[2]=sigl2
+    sig3[2]=10.**sig3[2]
+    sig3[3]=sigl3
+    sig3[3]=10.**sig3[3]
+    mu3[1]=mul1
+    mu3[2]=mul2
+    mu3[3]=mul3
+    eps3[1]=epsl1
+    eps3[2]=epsl2
+    eps3[3]=epsl3
+    thick3[1]=h1
+    thick3[2]=h2
+
+    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
+
+
+def run(plotIt=True, n=3):
+    # something to make a plot
+
+    F = frange(-5.,5.,20)
+    H = thick(50.,100.,n)
+    sign = sig(-5.,0.,n)
+    mun = mu(1.,2.,n)
+    epsn = eps(1.,9.,n)
+    chg = np.zeros_like(sign)
+    taux = np.zeros_like(sign)
+    c = np.zeros_like(sign)
+
+    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
+
+    if plotIt:
+
+        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
+
+    return Res, Phase
+
+if __name__ == '__main__':
+    run(plotIt=True)
+
+
+
+
+

SimPEG/Examples/__init__.py

@@ -18,10 +18,12 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
+import MT_1D_analytic_nlayer_Earth
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
+import sphereElectrostatic_example
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "sphereElectrostatic_example"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Examples/sphereElectrostatic_example.py

@@ -0,0 +1,785 @@
+from scipy.constants import epsilon_0
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import numpy as np
+from SimPEG.Utils import ndgrid, mkvc
+
+'''
+Authors: Thibaut Astic, Lindsey Heagy, Sanna Tyrvainen, Ronghua Peng
+
+
+This code defines function to resolve analytically the electrostatic sphere problem.
+We first define a problem configuration, with a conductive or resistive sphere in a
+wholespace background.
+We then calculate the potential, then the electric field, then the current density and
+finally the charges accumulation.
+
+Several plotting functions are defined for data visualisation.
+
+
+'''
+
+# Plot options
+ftsize_title = 18      #font size for titles
+ftsize_axis  = 14      #font size for axis ticks
+ftsize_label = 14      #font size for axis labels
+
+# Radius function, useful sigma ratio, and log scale converter
+r  = lambda x,y,z: np.sqrt(x**2.+y**2.+z**2.)
+sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0)
+
+#tools to convert log conductivity in conductivity
+def conductivity_log_wrapper(log_sig0,log_sig1):
+    sig0 = 10.**log_sig0
+    sig1 = 10.**log_sig1
+
+    return sig0,sig1
+
+# Examples
+#Plot the configuration. Label=False is used to generate a general case figure
+def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
+    '''
+    XYZ: ndgrid
+    sig0: conductivity of the background
+    sig1: conductivity of the sphere
+    R: radius of the sphere
+    E0: Amplitude of the uniform electrostatic field
+    ax: ax where to plot the configuration
+    label: True: plot real values, False: plot general case
+    colorsphere: color of the sphere, format [x,x,x]
+    '''
+
+    xplt = np.linspace(-R, R, num=100)
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    dx = xr[1]-xr[0]
+    top = np.sqrt(R**2-xplt**2)
+    bot = -np.sqrt(R**2-xplt**2)
+
+    if R != 0:
+        ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5)
+        ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 )
+        ax.arrow(0.,0.,np.sqrt(2.)*R/2.,np.sqrt(2.)*R/2.,head_width=0.,head_length=0.)
+
+        if label:
+            ax.annotate(("$\sigma_1$=%3.3f mS/m")%(sig1*10.**(3.)),
+            xy=(0.,-R/2.), xycoords='data',
+            xytext=(0.,-R/2.), textcoords='data',
+            fontsize=14.)
+            ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
+            xy=(0.,-1.5*R), xycoords='data',
+            xytext=(0.,-1.5*R), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$R$ = %1i m')%(R),
+            xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
+            xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
+            fontsize=14.)
+            ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+            ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+            ax.tick_params(labelsize=ftsize_axis)
+
+        else:
+            ax.set_xticklabels([])
+            ax.set_yticklabels([])
+            ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14)
+            ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14)
+            ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$R$'),
+            xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
+            xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
+            fontsize=14.)
+            ax.set_xlabel('x',fontsize=12)
+            ax.set_ylabel('y',fontsize=12)
+
+    else:
+        if label:
+            ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
+            xy=(0.,-1.5*R), xycoords='data',
+            xytext=(0.,-1.5*R), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$\mathbf{E_0} = %1i  \mathbf{\hat{x}}$ V/m')%(E0),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+            ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+            ax.tick_params(labelsize=ftsize_axis)
+
+        else:
+            ax.set_xticklabels([])
+            ax.set_yticklabels([])
+            ax.text(-0.05,-10,'$\sigma_0$',fontsize=14)
+            ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14)
+            ax.set_xlabel('x',fontsize=12)
+            ax.set_ylabel('y',fontsize=12)
+
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    [ax.arrow(xr.min(),_,np.abs(xr.max()-xr.min())/20.,0.,head_width=5.,head_length=2.,color='k') for _ in np.linspace(yr.min(),yr.max(),num=10)]
+    ax.patch.set_facecolor([0.4,0.7,0.4])
+    ax.patch.set_alpha(0.2)
+
+    ax.set_aspect('equal')
+
+
+
+    return ax
+
+def get_Conductivity(XYZ,sig0,sig1,R):
+    '''
+    Define the conductivity for each point of the space
+    '''
+    x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
+    r_view=r(x,y,z)
+
+    ind0= (r_view>R)
+    ind1= (r_view<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Sigma = np.zeros_like(x)
+
+    Sigma[ind0] = sig0
+    Sigma[ind1] = sig1
+
+    return Sigma
+
+
+def get_Potential(XYZ,sig0,sig1,R,E0):
+
+    '''
+    Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin
+    :input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
+    '''
+
+    x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
+
+    sig_cur = sigf(sig0,sig1)
+
+    r_cur = r(x,y,z)  # current radius
+
+    ind0 = (r_cur > R)
+    ind1 = (r_cur <= R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Vt = np.zeros_like(x)
+    Vp = np.zeros_like(x)
+    Vs = np.zeros_like(x)
+
+    Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere
+    Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0)             # inside the sphere
+
+
+    Vp = - E0*x  # primary potential
+
+    Vs = Vt - Vp # secondary potential
+
+    return Vt,Vp,Vs
+
+#plot the primary potential on ax
+def Plot_Primary_Potential(XYZ,Vp,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Primary Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+#plot the total potential on ax
+def Plot_Total_Potential(XYZ,Vt,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Total Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+#plot the secondary potential on ax
+def Plot_Secondary_Potential(XYZ,Vs,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Secondary Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+
+def get_ElectricField(XYZ,sig0,sig1,R,E0):
+    '''
+    Function that returns the total, the primary and the secondary electric fields,
+    input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    r_cur=r(x,y,z)  # current radius
+
+    ind0= (r_cur>R)
+    ind1= (r_cur<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Ep = np.zeros(shape=(len(x),3))
+    Ep[:,0] = E0
+
+    Et = np.zeros(shape=(len(x),3))
+
+    Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.);
+    Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1);
+    Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1);
+
+    Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0;
+    Et[ind1,1] = 0.;
+    Et[ind1,2] = 0.;
+
+    Es = Et - Ep
+
+    return Et, Ep, Es
+
+#plot the total electric field on ax
+def Plot_Total_ElectricField(XYZ,Et,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    EtXr = Et[:,0].reshape(xr.size, yr.size)
+    EtYr = Et[:,1].reshape(xr.size, yr.size)
+    EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Eplot = ax.pcolor(xr,yr,EtAmp)
+    cb = plt.colorbar(Eplot,ax=ax)
+    cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05)
+    ax.set_title('Total Field',fontsize=ftsize_title)
+
+
+    return ax
+
+#plot the secondary electric field on ax
+def Plot_Secondary_ElectricField(XYZ,Es,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    EsXr = Es[:,0].reshape(xr.size, yr.size)
+    EsYr = Es[:,1].reshape(xr.size, yr.size)
+    EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Eplot = ax.pcolor(xr,yr,EsAmp)
+    cb = plt.colorbar(Eplot,ax=ax)
+    cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Secondary Field',fontsize=ftsize_title)
+
+    return ax
+
+
+def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es):
+    '''
+    Function that returns the total, the primary and the secondary current densities,
+    :input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields,
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    r_cur=r(x,y,z)
+
+    ind0= (r_cur>R)
+    ind1= (r_cur<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Jt = np.zeros(shape=(len(x),3))
+    J0 = np.zeros(shape=(len(x),3))
+    Js = np.zeros(shape=(len(x),3))
+
+
+    Jp = sig0*Ep
+
+    Jt[ind0,:] = sig0*Et[ind0,:]
+    Jt[ind1,:] = sig1*Et[ind1,:]
+
+    Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:])
+    Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:]
+
+    return Jt,Jp,Js
+
+#plot the total currents density on ax
+def Plot_Total_Currents(XYZ,Jt,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    JtXr = Jt[:,0].reshape(xr.size, yr.size)
+    JtYr = Jt[:,1].reshape(xr.size, yr.size)
+    JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size))
+    cb = plt.colorbar(Jplot,ax=ax)
+    cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
+    ax.set_title('Total Current Density',fontsize=ftsize_title)
+
+    return ax
+
+
+#plot the secondary currents density on ax
+def Plot_Secondary_Currents(XYZ,Js,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    JsXr = Js[:,0].reshape(xr.size, yr.size)
+    JsYr = Js[:,1].reshape(xr.size, yr.size)
+    JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Js[:,2]**2).reshape(xr.size,yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size))
+    cb = plt.colorbar(Jplot,ax=ax)
+    cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
+    ax.set_title('Secondary Current Density',fontsize=ftsize_title)
+
+    return ax
+
+
+def get_ChargesDensity(XYZ,sig0,sig1,R,Ep):
+    '''
+    Function that returns the charges accumulation at the background/sphere interface,
+    :input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields,
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    dx = x[1]-x[0]
+
+    r_cur=r(x,y,z)
+
+    ind0 = (r_cur > R)
+    ind1 = (r_cur < R)
+    ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) )
+
+    assert (ind0 + ind1 + ind2).all(), 'Some indicies not included'
+
+    rho = np.zeros_like(x)
+
+    rho[ind0] = 0
+    rho[ind1] = 0
+    rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.))
+
+    return rho
+
+#Plot charges density on ax
+def Plot_ChargesDensity(XYZ,rho,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_aspect('equal')
+    Cplot = ax.pcolor(xr,yr,rho.reshape(xr.size, yr.size))
+    cb1 = plt.colorbar(Cplot,ax=ax)
+    cb1.set_label(label= 'Charge Density ($C/m^2$)',size=ftsize_label) #weight='bold')
+    cb1.ax.tick_params(labelsize=ftsize_axis)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_title('Charges Density', fontsize=ftsize_title)
+
+    return ax
+
+def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
+
+    '''
+    Function that return array of midpoints electrodes, electrodes positions,
+    potentials differences for total and secondary potentials fields, unormalized and
+    normalized to electrodes distances.
+    sig0: background conductivity
+    sig1: sphere conductivity
+    R: Sphere's radius
+    E0: uniform E field value
+    start: start point for the profile start.shape = (2,)
+    end: end point for the profile end.shape = (2,)
+    nbmp: number of dipoles
+    mn: Space between the M and N electrodes
+    '''
+
+    #D: total distance from start to end
+    D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.)
+
+    #MP: dipoles'midpoint positions (x,y)
+    MP = np.zeros(shape=(nbmp,2))
+    MP[:,0] = np.linspace(start[0],end[0],nbmp)
+    MP[:,1] = np.linspace(start[1],end[1],nbmp)
+
+    #Dipoles'Electrodes positions around each midpoints
+    EL = np.zeros(shape=(2*nbmp,2))
+    for n in range(0,len(EL),2):
+        EL[n,0]   = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2.
+        EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2.
+        EL[n,1]   = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2.
+        EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2.
+
+    VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL)
+    VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL)
+    dVtMP = np.zeros(nbmp)  #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP)
+    dVtMPn = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+    dVsMP = np.zeros(nbmp)  #Diffence (d-) of Secondaty Potential (Vt-) at each dipole (-MP)
+    dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+    dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP)
+    dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+
+    #Computing VtEL
+    for m in range(0,2*nbmp):
+        if (r(EL[m,0],EL[m,1],0) > R):
+            VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.)
+        else:
+            VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0)
+
+    #Computing VsEL
+    VsEL = VtEL + E0*EL[:,0]
+
+    #Computing dVtMP, dVsMP
+    for p in range(0,nbmp):
+        dVtMP[p] = VtEL[2*p]-VtEL[2*p+1]
+        dVtMPn[p] = dVtMP[p]/mn
+        dVsMP[p] = VsEL[2*p]-VsEL[2*p+1]
+        dVsMPn[p] = dVsMP[p]/mn
+
+    return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn
+
+#Compare the DC response of two configurations
+def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor):
+
+    #Define the mesh
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    #Defining the Profile
+    start = np.array([xstart,ystart])
+    end = np.array([xend,yend])
+
+    #Calculating the data from the defined survey line for Configuration 0 and 1
+    MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing)
+    MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing)
+
+
+    # Initializing the figure
+    fig = plt.figure(figsize=(20,20))
+    ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6)
+    ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6)
+    ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4)
+    ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6)
+    ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6)
+
+    #Plotting the Configuration 0
+    ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1])
+
+    #Plotting the Configuration 1
+    ax1   = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6])
+
+    #Plotting the Data (Legends)
+    ax2.set_title('Potential Differences',fontsize=ftsize_title)
+    ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
+    ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
+    ax2.tick_params(labelsize=ftsize_axis)
+    ax2.grid()
+
+    #Calculating the potential
+    Vt0,Vp0,Vs0 = get_Potential(XYZ,sig0,sig1,R0,E0)
+    Vt1,Vp1,Vs1 = get_Potential(XYZ,sig0,sig2,R1,E0)
+
+    if PlotOpt == 'Total':
+        ax3= Plot_Total_Potential(XYZ,Vt0,R0,ax3)
+        ax4= Plot_Total_Potential(XYZ,Vt1,R1,ax4)
+
+        #Plot the Data (from Configuration 0)
+        gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0
+                         ,marker='o',color='blue',linewidth=3.,label ='Left Model Response' )
+
+        #Plot the Data (from Configuration 1)
+        gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1
+                ,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
+        ax2.legend(('Left Model Response','Right Model Response'),loc=4)
+
+    elif PlotOpt == 'Secondary':
+        #plot the secondary potentials
+        ax3= Plot_Secondary_Potential(XYZ,Vt0,R0,ax3)
+        ax4= Plot_Secondary_Potential(XYZ,Vt1,R1,ax3)
+
+        #Plot the data(from configuration 0)
+        gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue'
+                ,marker='o',linewidth=3.,label ='Left Model Response' )
+
+
+        #Plot the Data (from Configuration 1)
+        gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1
+                 ,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
+        ax2.legend(('Left Model Response','Right Model Response'),loc=4 )
+
+    else:
+        print('What dont you get? Total or Secondary?')
+
+    #Legends
+    ax3.plot(MP0[:,0],MP0[:,1],color='gray')
+    Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black')
+    Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red')
+    ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
+
+    ax4.plot(MP1[:,0],MP1[:,1],color='gray')
+    Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black')
+    Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red')
+    ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
+
+    return fig
+
+#Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field
+def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b):
+
+    sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
+    E0   = 1.           # inducing field strength in V/m
+    n = 100             #level of discretisation
+    xr = np.linspace(-200., 200., n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+
+    Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
+
+
+    fig, ax = plt.subplots(1,2,figsize=(18,6))
+
+    #Setup figure 1 with options Configuration, Total or Secondary,
+    #then Potential, ElectricField, Current Density or Charges Density
+    if Figure1a == 'Configuration':
+        ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6])
+
+    elif Figure1a == 'Total':
+
+        if Figure1b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Total_Potential(XYZ,Vt,R,ax[0])
+
+        elif Figure1b == 'ElectricField':
+            ax[0] = Plot_Total_ElectricField(XYZ,Et,R,ax[0])
+
+        elif Figure1b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Total_Currents(XYZ,Jt,R,ax[0])
+
+        elif Figure1b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[0])
+
+    elif Figure1a == 'Secondary':
+
+        if Figure1b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Secondary_Potential(XYZ,Vs,R,ax[0])
+
+        elif Figure1b == 'ElectricField':
+            ax[0] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[0])
+
+        elif Figure1b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Secondary_Currents(XYZ,Js,R,ax[0])
+
+        elif Figure1b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[0])
+
+
+    if Figure1a== 'Configuration':
+        Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+        ax[1] = Plot_Primary_Potential(XYZ,Vp,R,ax[1])
+        print 'While figure1 is plotting Configuration, figure2 plots the primary field'
+
+    elif Figure2a == 'Total':
+    
+        if Figure2b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Total_Potential(XYZ,Vt,R,ax[1])
+
+        elif Figure2b == 'ElectricField':
+            ax[0] = Plot_Total_ElectricField(XYZ,Et,R,ax[1])
+
+        elif Figure2b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Total_Currents(XYZ,Jt,R,ax[1])
+
+        elif Figure2b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[1])
+
+
+    elif Figure2a == 'Secondary':
+        
+        if Figure2b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Secondary_Potential(XYZ,Vs,R,ax[1])
+
+        elif Figure2b == 'ElectricField':
+            ax[0] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[1])
+
+        elif Figure2b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Secondary_Currents(XYZ,Js,R,ax[1])
+
+        elif Figure2b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[1])
+
+    plt.tight_layout(True)
+    plt.show()
+
+#Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey
+def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example):
+
+    sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
+    sig2 = 10.**log_sig2
+    E0   = 1.           # inducing field strength in V/m
+    n = 100             #level of discretisation
+    xr = np.linspace(-200., 200., n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+    PlotOpt = 'Total'
+
+    if matching_spheres_example:
+        sig0 = 10.**(-3)
+        sig1 = 10.**(-2)
+        sig2 = 1.310344828 * 10**(-3)
+        R0   = 20.
+        R1   = 40.
+
+        two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
+
+    else:
+        two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
+
+    plt.tight_layout(True)
+    plt.show()
+
+def run(plotIt=True):
+    sig0 = -3.          # conductivity of the wholespace
+    sig1 = -1.         # conductivity of the sphere
+    sig0, sig1 = conductivity_log_wrapper(sig0,sig1)
+    R    = 50.          # radius of the sphere
+    E0   = 1.           # inducing field strength
+    n = 100             #level of discretisation
+    xr = np.linspace(-2.*R, 2.*R, n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+            
+    Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+    Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
+    Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+    rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+
+    if plotIt:
+        fig, ax = plt.subplots(2,5,figsize=(50,10))
+        ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1])
+        ax[1,0] = Plot_Primary_Potential(XYZ,Vp,R,ax[1,0])
+        ax[0,1] = Plot_Total_Potential(XYZ,Vt,R,ax[0,1])
+        ax[1,1] = Plot_Secondary_Potential(XYZ,Vs,R,ax[1,1])
+        ax[0,2] = Plot_Total_ElectricField(XYZ,Et,R,ax[0,2])
+        ax[1,2] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[1,2])
+        ax[0,3] = Plot_Total_Currents(XYZ,Jt,R,ax[0,3])
+        ax[1,3] = Plot_Secondary_Currents(XYZ,Js,R,ax[1,3])
+        ax[0,4] = Plot_Primary_Potential(XYZ,Vp,R,ax[0,4])
+        ax[1,4] = Plot_ChargesDensity(XYZ,rho,R,ax[1,4])
+
+
+        plt.show()
+
+if __name__ == '__main__':
+    run()
+
+

SimPEG/Mesh/CylMesh.py

@@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                 raise NotImplementedError('wrapping in the averaging is not yet implemented')
         return self._aveF2CCV
 
-    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
         """
             Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
         """
@@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
 
 
+        if locTypeTo is None:
+            locTypeTo = locType
+
         if locType == 'F':
             # do this three times for each component
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
-            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
             return sp.vstack((X,Y,Z))
         if locType == 'E':
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
-            Z = spzeros(Mrect.nEz, self.nE)
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
+            Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
             return sp.vstack((X,Y,Z))
 
-        grid = getattr(Mrect, 'grid' + locType)
+        grid = getattr(Mrect, 'grid' + locTypeTo)
         # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
         theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
         theta[theta < 0] += np.pi*2.0
@@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                         'Ex': Mrect.tangents[:Mrect.nEx,:],
                         'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
                         'Ez': Mrect.tangents[-Mrect.nEz:,:],
-                    }[locType]
+                    }[locTypeTo]
             if 'F' in locType:
                 normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
                 proj = ( normals * dotMe ).sum(axis=1)

SimPEG/Optimization.py

@@ -1003,7 +1003,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 # perturb inactive set off of bounds so that they are included in the step
                 delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
 
-
             # Only keep gradients going in the right direction on the active set
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.

SimPEG/Regularization.py

@@ -311,6 +311,9 @@ class BaseRegularization(object):
             tmp = indActive
             indActive = np.zeros(mesh.nC, dtype=bool)
             indActive[tmp] = True
+        if indActive is not None and mapping is None:
+            mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
         self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)

SimPEG/Utils/__init__.py

@@ -7,3 +7,4 @@ from CounterUtils import *
 import ModelBuilder
 import SolverUtils
 from coordutils import *
+from plottingUtils import *

SimPEG/Utils/plottingUtils.py

@@ -0,0 +1,3 @@
+# Plot Tree!
+# Plot SphereSetup
+# Plot LayerEarth

docs/examples/DC_Forward_PseudoSection.rst

@@ -22,7 +22,7 @@ radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
 stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
 dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
-Created by @fourndo on Mon Feb 01 19:28:06 2016
+Created by @fourndo
 
 
 

docs/examples/MT_1D_analytic_nlayer_Earth.rst

@@ -0,0 +1,21 @@
+.. _examples_MT_1D_analytic_nlayer_Earth:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+MT 1D analytic nlayer Earth
+===========================
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.MT_1D_analytic_nlayer_Earth.run()
+
+.. literalinclude:: ../../SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+    :language: python
+    :linenos:

docs/examples/sphereElectrostatic_example.rst

@@ -0,0 +1,21 @@
+.. _examples_sphereElectrostatic_example:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+sphereElectrostatic example
+===========================
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.sphereElectrostatic_example.run()
+
+.. literalinclude:: ../../SimPEG/Examples/sphereElectrostatic_example.py
+    :language: python
+    :linenos:

tests/base/test_regularization.py

@@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase):
                     elif mesh.dim == 3:
                         indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
-
                     for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
-                        reg = r(mesh, mapping=mapping, indActive=indAct)
+                        reg = r(mesh, indActive=indAct)
                         m = np.random.rand(mesh.nC)[indAct]
                         reg.mref = np.ones_like(m)*np.mean(m)
 

tests/mesh/test_cylMesh.py

@@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase):
 
         assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
 
+    def test_getInterpMatCartMesh_Cells2Nodes(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
     def test_getInterpMatCartMesh_Faces(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase):
         assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
+    def test_getInterpMatCartMesh_Faces2Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E')
+        mf = np.ones(Mc.nF)
+
+        ecart = Pf2e * mf
+
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 1) < TOL
+        assert np.abs(float(excc[indY]) - 0) < TOL
+        assert np.abs(float(eycc[indX]) - 0) < TOL
+        assert np.abs(float(eycc[indY]) - 1) < TOL
+        assert np.abs((ezcc - 1).sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
     def test_getInterpMatCartMesh_Edges(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase):
         Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
         me = np.ones(Mc.nE)
 
-        erect = Pe * me
+        ecart = Pe * me
 
-        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
-        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
-        ezcc = Mr.r(erect, 'E', 'Ez')
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges2Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F')
+        me = np.ones(Mc.nE)
+
+        frect = Pe2f * me
+
+        excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        ezcc = Mr.r(frect, 'F', 'Fz')
 
         indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
         indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])