diff --git a/simpegDCIP/Dev/DC2D_Demo_HTML_anim.py b/simpegDCIP/Dev/DC2D_Demo_HTML_anim.py
new file mode 100644
index 00000000..9a00257d
--- /dev/null
+++ b/simpegDCIP/Dev/DC2D_Demo_HTML_anim.py
@@ -0,0 +1,217 @@
+import os
+from SimPEG import np, sp, Utils, Mesh, mkvc
+import simpegDCIP as DC
+import pylab as plt
+#from ipywidgets import interact, IntSlider
+from matplotlib import animation
+from JSAnimation import HTMLWriter
+import time
+import re
+from readUBC_DC2DMesh import readUBC_DC2DMesh
+from readUBC_DC2DModel import readUBC_DC2DModel
+from readUBC_DC2DLoc import readUBC_DC2DLoc
+from convertObs_DC3D_to_2D import convertObs_DC3D_to_2D
+from readUBC_DC3Dobs import readUBC_DC3Dobs
+
+#%%
+home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Two_Sphere'
+msh_file = 'Mesh_2D.msh'
+mod_file = 'Model_2D.con'
+obs_file = 'FWR_data3D.dat'
+
+dsep = '\\'
+
+# Forward solver
+slvr = 'BiCGStab' #'LU'
+
+# Preconditioner
+pcdr = 'Jacobi'#'Gauss-Seidel'#
+
+# Number of padding cells to remove from plotting
+padc = 15
+
+# Load UBC mesh 2D
+mesh = readUBC_DC2DMesh(home_dir + dsep + msh_file)
+
+# Load model
+model = readUBC_DC2DModel(home_dir + dsep + mod_file)
+
+# load obs file
+[Tx,Rx,d,wd] = readUBC_DC3Dobs(home_dir + dsep + obs_file)
+[Tx, Rx] = convertObs_DC3D_to_2D(Tx,Rx)
+#%% Create system
+#Set boundary conditions
+mesh.setCellGradBC('neumann')
+
+Div = mesh.faceDiv
+Grad = mesh.cellGrad
+Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
+
+A = Div*Msig*Grad
+
+# Change one corner to deal with nullspace
+A[0,0] = 1
+A = sp.csc_matrix(A)
+
+start_time = time.time()
+
+if re.match(slvr,'BiCGStab'):
+ # Create Jacobi Preconditioner
+ if re.match(pcdr,'Jacobi'):
+ dA = A.diagonal()
+ P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
+
+ # Create Gauss-Seidel Preconditioner
+ elif re.match(pcdr,'Gauss-Seidel'):
+ LD = sp.tril(A,k=0)
+ #LDinv = sp.linalg.splu(LD)
+
+elif re.match(slvr,'LU'):
+ # Factor A matrix
+ Ainv = sp.linalg.splu(A)
+ print("LU DECOMP--- %s seconds ---" % (time.time() - start_time))
+
+#%% Create SimPEG objects
+
+
+# Create sub-mesh for plotting
+hx = mesh.hx
+hy = mesh.hy
+
+hx_sub = hx[padc:-padc]
+hy_sub = hy[padc:]
+mesh_sub = Mesh.TensorMesh([hx_sub,hy_sub],(mesh.vectorNx[padc], mesh.vectorNy[padc]))
+model_sub = model.reshape(mesh.nCy,mesh.nCx)
+model_sub = mkvc(model_sub[padc:,padc:-padc].T)
+
+xx = mesh_sub.vectorCCx
+yy = mesh_sub.vectorCCy
+
+#%% Solve
+#txii = range(50,1950,100)
+#jx_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
+#jy_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
+
+fig = plt.figure(figsize=(10,5))
+axs = plt.axes(ylim = (yy[0],yy[-1]+mesh.hy[-1]*2), xlim = (xx[0],xx[-1]))#
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.ylim(yy[0],yy[-1]+mesh.hy[-1]*2)
+plt.xlim(xx[0],xx[-1])
+#im1 = axs.pcolormesh([],[],[], alpha=0.75,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',vmin=-1e-2, vmax=1e-2)
+#im2 = axs.pcolormesh([],[],[],alpha=0.2,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',cmap='gray')
+im1 = axs.pcolormesh(mesh_sub.vectorCCx,mesh_sub.vectorCCy,np.zeros((mesh_sub.nCy,mesh_sub.nCx)), alpha=0.75,vmin=-1e-2, vmax=1e-2)
+im2 = axs.pcolormesh(mesh_sub.vectorCCx,mesh_sub.vectorCCy,np.zeros((mesh_sub.nCy,mesh_sub.nCx)), alpha=0.75,vmin=-1e-2, vmax=1e-2)
+im3 = axs.streamplot(xx, yy, np.zeros((mesh_sub.nCy,mesh_sub.nCx)), np.zeros((mesh_sub.nCy,mesh_sub.nCx)),color='k')
+im4 = axs.scatter([],[], c='r', s=200)
+im5 = axs.scatter([],[], c='r', s=200)
+
+
+
+#==============================================================================
+# def init():
+# im1.set_data([[],[],[]])
+# im2.set_data([[],[],[]])
+#
+# return [im1]+[im2]
+#==============================================================================
+
+def animate(ii):
+
+
+
+#for ii in range(len(txii)):
+
+ removeStream()
+ tx = np.asarray(np.c_[Tx[ii],np.ones(Tx[ii].shape[0])*mesh.vectorNy[-1]-1])
+ inds = Utils.closestPoints(mesh, tx )
+ RHS = mesh.getInterpolationMat( tx , 'CC').T*( [-1,1] / mesh.vol[inds] )
+
+ if re.match(slvr,'BiCGStab'):
+
+ if re.match(pcdr,'Jacobi'):
+ dA = A.diagonal()
+ P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
+
+ # Iterative Solve
+ phi = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
+ phi = mkvc(phi[0])
+
+ elif re.match(slvr,'LU'):
+ #Direct Solve
+ phi = Ainv.solve(RHS)
+
+
+ j = -Msig*Grad*phi
+ j_CC = mesh.aveF2CCV*j
+
+ # Compute charge density solving div*grad*phi
+ Q = -mesh.faceDiv*mesh.cellGrad*phi
+
+ jx_CC = j_CC[0:mesh.nC].reshape(mesh.nCy,mesh.nCx)
+ jy_CC = j_CC[mesh.nC:].reshape(mesh.nCy,mesh.nCx)
+
+ #%% Grab only the core for presentation
+ jx_CC_sub = jx_CC[padc:,padc:-padc]
+
+ jy_CC_sub = jy_CC[padc:,padc:-padc]
+
+ Q_sub = Q.reshape(mesh.nCy,mesh.nCx)
+ Q_sub = Q_sub[padc:,padc:-padc]
+
+ J_rho = np.sqrt(jx_CC_sub**2 + jy_CC_sub**2)
+ lw = np.log10(J_rho/J_rho.min())
+
+
+ #axs.imshow(Q_sub,alpha=0.75,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',vmin=-1e-2, vmax=1e-2)
+ #axs.imshow(np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)),alpha=0.2,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',cmap='gray')
+ global im1
+ im1 = axs.pcolormesh(mesh_sub.vectorCCx,mesh_sub.vectorCCy,Q_sub, alpha=0.75,vmin=-1e-2, vmax=1e-2)
+
+ global im2
+ im2 = axs.pcolormesh(mesh_sub.vectorCCx,mesh_sub.vectorCCy,np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)), alpha=0.25)
+
+ global im3
+ im3 = axs.streamplot(xx, yy, jx_CC_sub, jy_CC_sub,color='k',linewidth = lw,density=0.5)
+
+ global im4
+ im4 = axs.scatter(tx[0,0],mesh.vectorNy[-1], c='r', s=75, marker='v' )
+
+ global im5
+ im5 = axs.scatter(tx[1,0],mesh.vectorNy[-1], c='b', s=75, marker='v' )
+
+
+ #plt.show()
+ #im1.set_array(Q_sub)
+ #im2.set_array(np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)))
+ #im2.set_array(mesh_sub.vectorCCx, mesh_sub.vectorCCy,jx_CC_sub.T,jy_CC_sub.T)
+
+ #return [im1] + [im2]
+#%% Create widget
+def removeStream():
+ global im1
+ im1.remove()
+
+ global im2
+ im2.remove()
+
+ global im3
+ im3.lines.remove()
+ axs.patches = []
+
+ global im4
+ im4.remove()
+
+ global im5
+ im5.remove()
+#def viewInv(msh,iteration):
+
+
+
+#, linewidth=lw.T
+#%%
+#interact(viewInv,msh = mesh_sub, iteration = IntSlider(min=0, max=len(txii)-1 ,step=1, value=0))
+# set embed_frames=True to embed base64-encoded frames directly in the HTML
+anim = animation.FuncAnimation(fig, animate,
+ frames=len(Tx), interval=5)
+
+anim.save(home_dir + '\\animation.html', writer=HTMLWriter(embed_frames=True))
diff --git a/simpegDCIP/Dev/DC2D_fwr_obs.py b/simpegDCIP/Dev/DC2D_fwr_obs.py
deleted file mode 100644
index a9e161a6..00000000
--- a/simpegDCIP/Dev/DC2D_fwr_obs.py
+++ /dev/null
@@ -1,139 +0,0 @@
-import os
-
-home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegdc\simpegDCIP\Dev'
-
-os.chdir(home_dir)
-
-#%%
-from SimPEG import np, Utils, Mesh, mkvc, SolverLU
-import simpegDCIP as DC
-import pylab as plt
-#from ipywidgets import interact, IntSlider
-from matplotlib import animation
-from JSAnimation import HTMLWriter
-
-from readUBC_DC2DMesh import readUBC_DC2DMesh
-from readUBC_DC2DModel import readUBC_DC2DModel
-from readUBC_DC2DLoc import readUBC_DC2DLoc
-
-# Number of padding cells to remove from plotting
-padc = 16
-
-# Load UBC mesh 2D
-mesh = readUBC_DC2DMesh('mesh2d_fine.txt')
-
-# Load model
-model = readUBC_DC2DModel('model2d_fine.con')
-
-# load obs file
-[txLoc,rxLoc,d,wd] = readUBC_DC2DLoc('obs2d_East.loc')
-
-# Create SimPEG objects
-rx = DC.RxDipole(rxLoc[:,0], rxLoc[:,1])
-#tx = DC.SrcDipole([rx],txLoc[200,0],txLoc[200,1])
-
-# Create sub-mesh for plotting
-hx = mesh.hx
-hy = mesh.hy
-
-hx_sub = hx[padc:-padc]
-hy_sub = hy[padc:]
-mesh_sub = Mesh.TensorMesh([hx_sub,hy_sub],(hx_sub[0], -sum(hy_sub)))
-model_sub = model.reshape(mesh.nCy,mesh.nCx)
-model_sub = mkvc(model_sub[padc:,padc:-padc].T)
-
-xx = mesh_sub.vectorCCx
-yy = mesh_sub.vectorCCy
-
-#%% Solve
-txii = range(50,1950,100)
-#jx_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
-#jy_CC_sub = np.zeros((len(txii),mesh_sub.nCx,mesh_sub.nCy))
-
-fig = plt.figure(figsize=(10,5))
-axs = plt.axes(ylim=(-800,50), xlim=(25,2000))
-plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
-im1 = axs.imshow([[],[]], alpha=0.75,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',vmin=-1e-2, vmax=1e-2)
-im2 = axs.imshow([[],[]],alpha=0.2,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',cmap='gray')
-im3 = axs.streamplot(mesh_sub.vectorCCx, mesh_sub.vectorCCy, np.zeros((mesh_sub.nCy,mesh_sub.nCx)), np.zeros((mesh_sub.nCy,mesh_sub.nCx)),color='k')
-im4 = axs.scatter([],[], c='r', s=200)
-
-def init():
- im1.set_data([[],[]])
- im2.set_data([[],[]])
-
- return [im1]+[im2]
-
-def animate(ii):
-
-
-
-#for ii in range(len(txii)):
-
- removeStream()
-
- tx = DC.SrcDipole([rx],txii[ii],txii[ii])
-
- survey = DC.SurveyDC([tx])
- problem = DC.ProblemDC_CC(mesh)
- problem.pair(survey)
- problem.Solver = SolverLU
-
- u1 = problem.fields(model)
-
- Msig1 = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
-
- j = -Msig1*mesh.cellGrad*u1[tx, 'phi_sol']
- j_CC = mesh.aveF2CCV*j
-
- # Compute charge density solving div*grad*phi
- Q = -mesh.faceDiv*mesh.cellGrad*u1[tx, 'phi_sol']
-
- jx_CC = j_CC[0:mesh.nC].reshape(mesh.nCy,mesh.nCx).T
- jy_CC = j_CC[mesh.nC:].reshape(mesh.nCy,mesh.nCx).T
-
- #%% Grab only the core for presentation
- jx_CC_sub = jx_CC[padc:-padc,padc:]
- jy_CC_sub = jy_CC[padc:-padc,padc:]
- Q_sub = Q.reshape(mesh.nCy,mesh.nCx)
- Q_sub = Q_sub[padc:,padc:-padc]
-
- J_rho = np.sqrt(jx_CC_sub**2 + jy_CC_sub**2)
- lw = np.log10(J_rho/J_rho.min())
-
- #axs.imshow(Q_sub,alpha=0.75,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',vmin=-1e-2, vmax=1e-2)
- #axs.imshow(np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)),alpha=0.2,extent = (xx[0],xx[-1],yy[-1],yy[0]),interpolation='nearest',cmap='gray')
-
-
- global im3
- im3 = axs.streamplot(mesh_sub.vectorCCx, mesh_sub.vectorCCy, jx_CC_sub.T, jy_CC_sub.T,color='k',linewidth = lw.T,density=1.25)
-
- global im4
- im4 = axs.scatter(txii[ii],10, c='r', s=60, marker='+' )
-
- #plt.show()
- im1.set_array(Q_sub)
- im2.set_array(np.log10(model_sub.reshape(mesh_sub.nCy,mesh_sub.nCx)))
- #im2.set_array(mesh_sub.vectorCCx, mesh_sub.vectorCCy,jx_CC_sub.T,jy_CC_sub.T)
-
- return [im1] + [im2]
-#%% Create widget
-def removeStream():
- global im3
- im3.lines.remove()
- axs.patches = []
-
- global im4
- im4.remove()
-#def viewInv(msh,iteration):
-
-
-
-#, linewidth=lw.T
-#%%
-#interact(viewInv,msh = mesh_sub, iteration = IntSlider(min=0, max=len(txii)-1 ,step=1, value=0))
-# set embed_frames=True to embed base64-encoded frames directly in the HTML
-anim = animation.FuncAnimation(fig, animate, init_func=init,
- frames=len(txii), interval=10)
-
-anim.save('animation.html', writer=HTMLWriter(embed_frames=True))
diff --git a/simpegDCIP/Dev/DC3D_Demo_TwoSpheres.py b/simpegDCIP/Dev/DC3D_Demo_TwoSpheres.py
index b6aea621..5b661469 100644
--- a/simpegDCIP/Dev/DC3D_Demo_TwoSpheres.py
+++ b/simpegDCIP/Dev/DC3D_Demo_TwoSpheres.py
@@ -43,7 +43,7 @@ dsep = '\\'
#from scipy.linalg import solve_banded
# Load UBC mesh 3D
-mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_10m.msh')
+mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_5m.msh')
#mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\MtIsa_20m.msh')
#mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_50m.msh')
@@ -74,6 +74,9 @@ flr = 1e-4
chifact = 100
ref_mod = 1e-2
+# DOI threshold
+cutoff = 0.8
+
#%% Create system
#Set boundary conditions
mesh.setCellGradBC('neumann')
@@ -276,7 +279,7 @@ if not re.match(stype,'gradient'):
# m3D = np.reshape(model, (mesh.nCz, mesh.nCy, mesh.nCx))
# m2D = m3D[:,1,:]
#==============================================================================
-
+ #%%
plt.figure()
axs = plt.subplot(1,1,1)
@@ -284,8 +287,8 @@ if not re.match(stype,'gradient'):
plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]+2*dx])
plt.gca().set_aspect('equal', adjustable='box')
- circle1=plt.Circle((150,1500),50,color='w',fill=False, lw=3)
- circle2=plt.Circle((325,1500),50,color='k',fill=False, lw=3)
+ circle1=plt.Circle((144,1500),50,color='w',fill=False, lw=3)
+ circle2=plt.Circle((344,1500),50,color='k',fill=False, lw=3)
axs.add_artist(circle1)
axs.add_artist(circle2)
plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(m2D))#axes = [mesh2d.vectorNx[0],mesh2d.vectorNx[-1],mesh2d.vectorNy[0],mesh2d.vectorNy[-1]])
@@ -307,66 +310,130 @@ if not re.match(stype,'gradient'):
plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]+2*dx])
plt.gca().set_aspect('equal', adjustable='box')
- circle1=plt.Circle((150,1500),50,color='w',fill=False, lw=3)
- circle2=plt.Circle((325,1500),50,color='k',fill=False, lw=3)
+ circle1=plt.Circle((144,1500),50,color='w',fill=False, lw=3)
+ circle2=plt.Circle((344,1500),50,color='k',fill=False, lw=3)
axs.add_artist(circle1)
axs.add_artist(circle2)
plot_pseudoSection(Tx2d,Rx2d,data,nz[-1],stype)
plt.show()
- #%% Create dcin2d inversion files and run
- inv_dir = home_dir + '\Inv2D'
- if not os.path.exists(inv_dir):
- os.makedirs(inv_dir)
+ #%% Run two inversions with different reference models and compute a DOI
+
+ invmod = []
+ refmod = []
+ plt.figure()
+
+ for jj in range(2):
+
+ # Create dcin2d inversion files and run
+ inv_dir = home_dir + '\Inv2D'
+ if not os.path.exists(inv_dir):
+ os.makedirs(inv_dir)
+
+ mshfile2d = 'Mesh_2D.msh'
+ modfile2d = 'Model_2D.con'
+ obsfile2d = 'FWR_3D_2_2D.dat'
+ inp_file = 'dcinv2d.inp'
- mshfile2d = 'Mesh_2D.msh'
- modfile2d = 'MtIsa_2D.con'
- obsfile2d = 'FWR_3D_2_2D.dat'
- inp_file = 'dcinv2d.inp'
+
+ # Export 2D mesh
+ fid = open(inv_dir + dsep + mshfile2d,'w')
+ fid.write('%i\n'% mesh2d.nCx)
+ fid.write('%f %f 1\n'% (mesh2d.vectorNx[0],mesh2d.vectorNx[1]))
+ np.savetxt(fid, np.c_[mesh2d.vectorNx[2:],np.ones(mesh2d.nCx-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+ fid.write('\n')
+ fid.write('%i\n'% mesh2d.nCy)
+ fid.write('%f %f 1\n'%( 0,mesh2d.hy[-1]))
+ np.savetxt(fid, np.c_[np.cumsum(mesh2d.hy[-2::-1])+mesh2d.hy[-1],np.ones(mesh2d.nCy-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+ fid.close()
+
+ # Export 2D model
+ fid = open(inv_dir + dsep + modfile2d,'w')
+ fid.write('%i %i\n'% (mesh2d.nCx,mesh2d.nCy))
+ np.savetxt(fid, mkvc(m2D[::-1,:].T), fmt='%e',delimiter=' ',newline='\n')
+ fid.close()
+
+ # Export data file
+ writeUBC_DCobs(inv_dir + dsep + obsfile2d,Tx2d,Rx2d,data,unct,'2D')
+
+ # Write input file
+ fid = open(inv_dir + dsep + inp_file,'w')
+ fid.write('OBS LOC_X %s \n'% obsfile2d)
+ fid.write('MESH FILE %s \n'% mshfile2d)
+ fid.write('CHIFACT 1 %f\n'% chifact)
+ fid.write('TOPO DEFAULT %s \n')
+ fid.write('INIT_MOD DEFAULT\n')
+ fid.write('REF_MOD VALUE %e\n'% (ref_mod*(jj+1)))
+ fid.write('ALPHA DEFAULT\n')
+ fid.write('WEIGHT DEFAULT\n')
+ fid.write('STORE_ALL_MODELS FALSE\n')
+ fid.write('INVMODE SVD\n')
+ fid.write('USE_MREF TRUE\n')
+ fid.close()
+
+ os.chdir(inv_dir)
+ os.system('dcinv2d ' + inp_file)
+
+
+ #Load model
+ minv = readUBC_DC2DModel(inv_dir + dsep + 'dcinv2d.con')
+
+ axs = plt.subplot(2,1,jj+1)
+
+ plt.xlim([-dx,nc*dx+dx])
+ plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]+2*dx])
+ plt.gca().set_aspect('equal', adjustable='box')
+
+ minv = np.reshape(minv,(mesh2d.nCy,mesh2d.nCx))
+ #plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(m2D),alpha=0.5, cmap='gray')
+
+ circle1=plt.Circle((144,1500),50,color='w',fill=False, lw=3)
+ circle2=plt.Circle((344,1500),50,color='k',fill=False, lw=3)
+ axs.add_artist(circle1)
+ axs.add_artist(circle2)
+
+
+ axp = plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(minv),alpha=1,vmin = -2.25, vmax = -1.5)
+
+ plt.show()
+
+ if jj == 1:
+
+ plt.ylabel('(b)',rotation=360)
+ plt.xlabel('Distance (m)')
+
+ else:
+ plt.ylabel('(a)',rotation=360)
+
+
+ cbar = plt.colorbar(format = '%.2f',fraction=0.05,orientation='vertical',pad=0.02)
+ cmin,cmax = cbar.get_clim()
+ ticks = np.linspace(cmin,cmax,3)
+ cbar.set_ticks(ticks)
+ #cbar.set_ticklabels('%.2f')
+
+ invmod.append(minv)
+ refmod.append(ref_mod*(jj+1))
+
+ #%% Compute DOI
+ DOI = np.abs(invmod[0] - invmod[1]) / np.abs(refmod[0] - refmod[1])
+ # Normalize between [0 1]
+ DOI = DOI - np.min(DOI)
+ DOI = (1.- DOI/np.max(DOI))
+ DOI[DOI > cutoff] = 1
+ plt.figure()
+ plt.xlim([-dx,nc*dx+dx])
+ plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]+2*dx])
+ plt.gca().set_aspect('equal', adjustable='box')
- # Export 2D mesh
- fid = open(inv_dir + dsep + mshfile2d,'w')
- fid.write('%i\n'% mesh2d.nCx)
- fid.write('%f %f 1\n'% (mesh2d.vectorNx[0],mesh2d.vectorNx[1]))
- np.savetxt(fid, np.c_[mesh2d.vectorNx[2:],np.ones(mesh2d.nCx-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
- fid.write('\n')
- fid.write('%i\n'% mesh2d.nCy)
- fid.write('%f %f 1\n'%( 0,mesh2d.hy[-1]))
- np.savetxt(fid, np.c_[np.cumsum(mesh2d.hy[-2::-1])+mesh2d.hy[-1],np.ones(mesh2d.nCy-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
- fid.close()
+ plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,DOI,alpha=1)
+ cbar = plt.colorbar(format = '%.2f',fraction=0.02)
- # Export 2D model
- fid = open(inv_dir + dsep + modfile2d,'w')
- fid.write('%i %i\n'% (mesh2d.nCx,mesh2d.nCy))
- np.savetxt(fid, mkvc(m2D[::-1,:].T), fmt='%e',delimiter=' ',newline='\n')
- fid.close()
-
- # Export data file
- writeUBC_DCobs(inv_dir + dsep + obsfile2d,Tx2d,Rx2d,data,unct,'2D')
-
- # Write input file
- fid = open(inv_dir + dsep + inp_file,'w')
- fid.write('OBS LOC_X %s \n'% obsfile2d)
- fid.write('MESH FILE %s \n'% mshfile2d)
- fid.write('CHIFACT 1 %f\n'% chifact)
- fid.write('TOPO DEFAULT %s \n')
- fid.write('INIT_MOD DEFAULT\n')
- fid.write('REF_MOD VALUE %e\n'% ref_mod)
- fid.write('ALPHA DEFAULT\n')
- fid.write('WEIGHT DEFAULT\n')
- fid.write('STORE_ALL_MODELS FALSE\n')
- fid.write('INVMODE SVD\n')
- fid.write('USE_MREF TRUE\n')
- fid.close()
-
- os.chdir(inv_dir)
- os.system('dcinv2d ' + inp_file)
-
- #%%
- #Load model
- minv = readUBC_DC2DModel(inv_dir + dsep + 'dcinv2d.con')
+ #%% Replace alpha values from inversion
+ #rgba_plt = axp.get_facecolor()
+ #rgba_plt[:,3] = mkvc(DOI)/2
plt.figure()
axs = plt.subplot(1,1,1)
@@ -374,22 +441,20 @@ if not re.match(stype,'gradient'):
plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]+2*dx])
plt.gca().set_aspect('equal', adjustable='box')
- minv = np.reshape(minv,(mesh2d.nCy,mesh2d.nCx))
- #plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(m2D),alpha=0.5, cmap='gray')
-
- circle1=plt.Circle((150,1500),50,color='w',fill=False, lw=3)
- circle2=plt.Circle((325,1500),50,color='k',fill=False, lw=3)
+ circle1=plt.Circle((144,1500),50,color='w',fill=False, lw=3)
+ circle2=plt.Circle((344,1500),50,color='k',fill=False, lw=3)
axs.add_artist(circle1)
axs.add_artist(circle2)
- axp = plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(minv),alpha=1,vmin = np.min(np.log10(minv)), vmax = np.max(np.log10(minv)))
- #t = [-3, -2, -1]
+ axs = plt.pcolor(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(invmod[0]),edgecolor="none")
+ plt.draw()
cbar = plt.colorbar(format = '%.2f',fraction=0.02)
- cmin,cmax = cbar.get_clim()
- ticks = np.linspace(cmin,cmax,3)
- cbar.set_ticks(ticks)
- #cbar.set_ticklabels('%.2f')
-
+ aa = axs.get_facecolors()
+ aa[:,3] = mkvc(DOI.T)
+ axs.set_facecolor(aa)
+
+ plt.draw()
+
#%% Othrwise it is a gradient array, plot surface of apparent resisitivty
elif re.match(stype,'gradient'):
diff --git a/simpegDCIP/Dev/DC3D_fwr_live_test.py b/simpegDCIP/Dev/DC3D_fwr_live_test.py
index 7641e679..a4a96d56 100644
--- a/simpegDCIP/Dev/DC3D_fwr_live_test.py
+++ b/simpegDCIP/Dev/DC3D_fwr_live_test.py
@@ -37,29 +37,30 @@ from gen_DCIPsurvey import gen_DCIPsurvey
from convertObs_DC3D_to_2D import convertObs_DC3D_to_2D
import os
-home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Two_Sphere'
+#home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Two_Sphere'
+home_dir ='C:\\LC\Private\\dominiquef\\Projects\\4414_Minsim\\Model'
dsep = '\\'
#from scipy.linalg import solve_banded
# Load UBC mesh 3D
-mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_10m.msh')
+#mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_10m.msh')
#mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\MtIsa_20m.msh')
-#mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_50m.msh')
+mesh = Utils.meshutils.readUBCTensorMesh(home_dir + '\Mesh_50m.msh')
# Load model
#model = Utils.meshutils.readUBCTensorModel(home_dir + '\MtIsa_3D.con',mesh)
#model = Utils.meshutils.readUBCTensorModel(home_dir + '\Synthetic.con',mesh)
-#model = Utils.meshutils.readUBCTensorModel(home_dir + '\Lalor_model_50m.con',mesh)
-model = Utils.meshutils.readUBCTensorModel(home_dir + '\TwoSpheres.con',mesh)
+model = Utils.meshutils.readUBCTensorModel(home_dir + '\Lalor_model_50m.con',mesh)
+#model = Utils.meshutils.readUBCTensorModel(home_dir + '\TwoSpheres.con',mesh)
#model = model**0 * 1e-2
# Specify survey type
-stype = 'dpdp'
+stype = 'pdp'
# Survey parameters
-a = 30
-b = 30
-n = 20
+a = 150
+b = 150
+n = 40
# Forward solver
slvr = 'BiCGStab' #'LU'
@@ -71,7 +72,7 @@ pcdr = 'Jacobi'#'Gauss-Seidel'#
pct = 0.01
flr = 1e-4
chifact = 100
-ref_mod = 1e-2
+ref_mod = 1e-3
#%% Create system
#Set boundary conditions
@@ -112,8 +113,8 @@ top = int(mesh.nCz)-1
plt.figure()
ax_prim = plt.subplot(1,1,1)
mesh.plotSlice(model, ind=top, normal='Z', grid=False, pcolorOpts={'alpha':0.5}, ax =ax_prim)
-plt.xlim([423000,424000])
-plt.ylim([546200,547000])
+#plt.xlim([423000,424000])
+#plt.ylim([546200,547000])
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
@@ -129,7 +130,7 @@ plt.sca(ax_prim)
# Takes two points from ginput and create survey
#if re.match(stype,'gradient'):
-# gin = [(423187. , 546311.), (423867. , 546991.)]
+#gin = [(425347, 6079766), (427792, 6081806)]
#else:
gin = plt.ginput(2, timeout = 0)
@@ -280,7 +281,7 @@ if not re.match(stype,'gradient'):
axs = plt.subplot(2,1,1)
plt.xlim([0,nc*dx])
- plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]])
+ plt.ylim([mesh2d.vectorNy[-1]-dl_len,mesh2d.vectorNy[-1]])
plt.gca().set_aspect('equal', adjustable='box')
plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(m2D),alpha=0.5, cmap='gray')#axes = [mesh2d.vectorNx[0],mesh2d.vectorNx[-1],mesh2d.vectorNy[0],mesh2d.vectorNy[-1]])
@@ -348,27 +349,29 @@ if not re.match(stype,'gradient'):
axs = plt.subplot(2,1,2)
plt.xlim([0,nc*dx])
- plt.ylim([mesh2d.vectorNy[-1]-dl_len/2,mesh2d.vectorNy[-1]])
+ plt.ylim([mesh2d.vectorNy[-1]-dl_len,mesh2d.vectorNy[-1]])
plt.gca().set_aspect('equal', adjustable='box')
minv = np.reshape(minv,(mesh2d.nCy,mesh2d.nCx))
plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(m2D),alpha=0.5, cmap='gray')
plt.pcolormesh(mesh2d.vectorNx,mesh2d.vectorNy,np.log10(minv),alpha=0.5, clim=(np.min(np.log10(m2D)),np.max(np.log10(m2D))))
- plt.colorbar
-
+ cbar = plt.colorbar(format = '%.2f',fraction=0.02)
+ cmin,cmax = cbar.get_clim()
+ ticks = np.linspace(cmin,cmax,3)
+ cbar.set_ticks(ticks)
#%% Othrwise it is a gradient array, plot surface of apparent resisitivty
elif re.match(stype,'gradient'):
rC1P1 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,0],Rx[0].shape[0], 1) - Rx[0][:,0:2])**2, axis=1 ))
rC2P1 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,1],Rx[0].shape[0], 1) - Rx[0][:,0:2])**2, axis=1 ))
- rC1P2 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,1],Rx[0].shape[0], 1) - Rx[0][:,3:5])**2, axis=1 ))
- rC2P2 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,0],Rx[0].shape[0], 1) - Rx[0][:,3:5])**2, axis=1 ))
+ rC1P2 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,0],Rx[0].shape[0], 1) - Rx[0][:,3:5])**2, axis=1 ))
+ rC2P2 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,1],Rx[0].shape[0], 1) - Rx[0][:,3:5])**2, axis=1 ))
rC1C2 = np.sqrt( np.sum( (npm.repmat(Tx[0][0:2,0]-Tx[0][0:2,1],Rx[0].shape[0], 1) )**2, axis=1 ))
rP1P2 = np.sqrt( np.sum( (Rx[0][:,0:2] - Rx[0][:,3:5])**2, axis=1 ))
- rho = np.abs(data[0]) * np.pi *((rC1P1)**2 / rP1P2)#/ ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 )
+ rho = np.abs(data[0]) *np.pi *2. / ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 )#*((rC1P1)**2 / rP1P2)#
Pmid = (Rx[0][:,0:2] + Rx[0][:,3:5])/2
@@ -378,7 +381,10 @@ elif re.match(stype,'gradient'):
#plt.subplot(2,1,2)
+ plt.figure()
plt.imshow(grid_rho.T, extent = (np.min(grid_x),np.max(grid_x),np.min(grid_z),np.max(grid_z)) ,origin='lower')
+
var = 'Gradient Array - a-spacing: ' + str(a) + ' m'
plt.title(var)
plt.colorbar()
+ plt.contour(grid_x,grid_z,grid_rho, colors='k')
\ No newline at end of file
diff --git a/simpegDCIP/Dev/animation.html b/simpegDCIP/Dev/animation.html
new file mode 100644
index 00000000..22cc20f8
--- /dev/null
+++ b/simpegDCIP/Dev/animation.html
@@ -0,0 +1,191 @@
+
+
+
+
+
+
+
+
+
–
+
+
+
+
+
+
+
+
+
+
+