Test script for multiple *.ts surfaces to block model (Lalor model)

simpegPF/BaseMag.py

@@ -401,6 +401,10 @@ def read_GOCAD_ts(tsfile):
     
     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 = []
     
@@ -486,13 +490,6 @@ def gocad2vtk(gcFile,mesh,bcflag,inflag):
     insideGrid = extractImpDistRectGridFilt.GetOutput()
     insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))
     
-    # Get index surface intersect
-    extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
-    extractImpDistRectGridFilt.ExtractInsideOff()
-    extractImpDistRectGridFilt.Update()
-    
-    bcGrid = extractImpDistRectGridFilt.GetOutput()
-    bcGrid = npsup.vtk_to_numpy(bcGrid.GetCellData().GetArray('Index'))
     
     # Return the indexes inside
-    return insideGrid, bcGrid
+    return insideGrid

simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py

@@ -222,12 +222,15 @@ for ii in range(d_iter):
     
 
 #%% Plot the forward solution from integral
-plt.figure()
+plt.figure(figsize=[8,4])
 ax = plt.subplot()
 plt.imshow(d2d, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
-plt.colorbar(fraction=0.04)
+plt.colorbar(fraction=0.02)
+plt.contour(X,Y, d2d,10)
+plt.scatter(X,Y, s=5)
 
-plt.figure()
+
+plt.figure(figsize=[8,4])
 plt.contour(X,Y, np.reshape(d,X.shape),10)
 plt.scatter(X,Y, c=np.reshape(d,X.shape), s=10)
 plt.scatter(xx,yy, c='k', s=20, marker='o')
@@ -236,26 +239,29 @@ ax.set_title('Numerical')
 
 #%%
 
-plt.figure()
+plt.figure(figsize=[7,9])
 ax = plt.subplot()
 plt.plot(xx,d_line,c='r', linewidth=3)
-plt.plot(xx,d_i2d_lin,c='b')
-plt.plot(xx,d_i2d_cub,c='g')
-plt.plot(xx,d_i2d_qui,c='m')
-plt.plot(xx,d_i2d_nnb,c='k')
-plt.plot(xx,d_i2d_CTI,c='c')
+plt.plot(xx,d_i2d_lin)
+plt.plot(xx,d_i2d_cub)
+plt.plot(xx,d_i2d_qui)
+plt.plot(xx,d_i2d_nnb)
+plt.plot(xx,d_i2d_CTI)
+plt.xlim(xx.min(),xx.max())
 
+plt.legend(['True','linear','Cubic','Quintic','NearestN','CloughTorcher'],bbox_to_anchor=(0.75, 0.25))
 # Plot interpolation from true value on line
 F = interpolation.interp1d(xx,d_line)
 dtrue = F(xr[::indx])
 plt.plot(xr[::indx],dtrue,c='r',linewidth=0.,marker='o')
-ax.set_title('Analytical')
+ax.set_title('Interpolated data profile')
 
 #%% Write result to file
-with file('l2_residual.dat','w') as fid:
-    fid.write('NearestN \t Linear \t Cubic \t Quintic \t FFT\n')
-    np.savetxt(fid, l2_r, fmt='%e',delimiter=' ',newline='\n')
-
-with file('l1_residual.dat','w') as fid:
-    fid.write('NearestN \t Linear \t Cubic \t Quintic \t FFT\n')
-    np.savetxt(fid, l1_r, fmt='%e',delimiter=' ',newline='\n')
+with file('Interp_residual.dat','w') as fid:
+    fid.write('NearestN\tLinear\tCubic\tQuintic\tCloughTocher\n')
+    fid.write('\nL2-norm\n')
+    np.savetxt(fid, l2_r, fmt='%5.3e',delimiter='\t',newline='\n')
+    fid.write('\nL1-norm\n')
+    np.savetxt(fid, l1_r, fmt='%5.3e',delimiter='\t',newline='\n')
+    fid.write('\nLinf-norm\n')
+    np.savetxt(fid, linf_r, fmt='%5.3e',delimiter='\t',newline='\n')

simpegPF/Dev/Intgrl_MAG_Fwr_VTK_3DModel.py

@@ -0,0 +1,105 @@
+import os
+
+home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG\\Lalor'
+
+inpfile = 'PYMAG3C_fwr.inp'
+
+dsep = '\\'
+
+os.chdir(home_dir)
+
+#%%
+from SimPEG import np, sp, Utils, Mesh 
+import scipy.interpolate as interpolation
+import simpegPF as PF
+import pylab as plt
+import time as tm
+
+## New scripts to be added to basecode
+#from fwr_MAG_data import fwr_MAG_data
+#from read_MAGfwr_inp import read_MAGfwr_inp
+
+#%%
+[mshfile, obsfile, modfile, magfile, topofile] = PF.BaseMag.read_MAGfwr_inp(inpfile)
+
+# Load mesh file
+mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
+
+# Load in observation file
+#[B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
+
+# Read in topo surface
+topsurf = "Topography_Local.ts"
+geosurf = [['_Top_Bollock_Rhyodacite_SUB.ts',False,True],
+['_Top_Bollock_Rhyodacite_SUB_B.ts',False,True],
+['_Top_Mafic_Intrusion.ts',False,True],
+['_Top_Mafic_Volcaniclastic.ts',False,True]]
+
+vals = np.asarray([0.005,0.005,0.01,0.0025])
+
+# Background density
+bkgr = 0.0001
+
+# Offset data above topo 
+zoffset = 2
+
+
+#%% Script starts here       
+# # Create a grid of observations and offset the z from topo
+#xr = np.linspace(mesh.vectorCCx[0], mesh.vectorCCx[-1], 99)
+#yr = np.linspace(mesh.vectorCCy[0], mesh.vectorCCy[-1], 74)
+
+xr = mesh.vectorCCx[::2]
+yr = mesh.vectorCCy[::2]
+X, Y = np.meshgrid(xr, yr)
+
+topo = np.genfromtxt(topofile,skip_header=1)
+F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
+Z = F(X,Y) + zoffset
+rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
+    
+ndata = rxLoc.shape[0]
+
+B = np.array(([90.,0.,50000.]))
+  
+M = np.array(([90.,0.,315.]))     
+
+model= np.ones(mesh.nC) * bkgr
+# Load GOCAD surf
+#[vrtx, trgl] = PF.BaseMag.read_GOCAD_ts(tsfile)
+# Find active cells from surface
+
+for ii in range(len(geosurf)):
+    tin = tm.time()
+    print "Computing indices with VTK: " + geosurf[ii][0]
+    indx = PF.BaseMag.gocad2vtk(geosurf[ii][0],mesh, bcflag = geosurf[ii][1], inflag = geosurf[ii][2])
+    print "VTK operation completed in " + str(tm.time() - tin)
+    
+    model[indx] = vals[ii]
+    
+indx = PF.BaseMag.gocad2vtk(topsurf,mesh, bcflag = False, inflag = True) 
+actv = np.zeros(mesh.nC)
+actv[indx] = 1
+
+model[actv==0] = -100
+
+Utils.meshutils.writeUBCTensorModel('VTKout.dat',mesh,model)
+
+Utils.meshutils.writeUBCTensorMesh('Mesh_temp.msh',mesh)
+
+start_time = tm.time()
+    
+d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
+
+timer = (tm.time() - start_time)
+
+#%% Plot data
+plt.figure()
+ax = plt.subplot()
+plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.clim(0,25)
+plt.colorbar(fraction=0.02)
+plt.contour(X,Y, np.reshape(d,X.shape),10)
+plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
+
+ax.set_title('Forward data')

simpegPF/Dev/Intgrl_MAG_Fwr_VTK_topo_test.py

@@ -60,9 +60,10 @@ dx = 10
 
 nc = int(sclx/dx)
 
-hxind = [(dx, 2*nc)]
-hyind = [(dx, nc)]
-hzind = [(dx, nc)]
+
+hxind = np.ones(2*nc)*dx
+hyind = np.ones(nc)*dx
+hzind = np.ones(nc)*dx
 
 mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
     
@@ -91,49 +92,88 @@ d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
 timer = (tm.time() - start_time)
 
 #%% Plot data
-plt.figure(1)
+plt.figure()
 ax = plt.subplot()
 plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+plt.clim(0,25)
 plt.colorbar(fraction=0.02)
 plt.contour(X,Y, np.reshape(d,X.shape),10)
 plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
+
 ax.set_title('Forward data')
 
-# Load model file
-#model = Utils.meshutils.readUBCTensorModel(modfile,mesh)
-  
-# Load in topofile or create flat surface
-#==============================================================================
-# if topofile == 'null':
-#  
-#     actv = np.ones(mesh.nC)   
-#     
-# else: 
-#     topo = np.genfromtxt(topofile,skip_header=1)
-#     actv = PF.Magnetics.getActiveTopo(mesh,topo,'N')
-# 
-# 
-# Utils.writeUBCTensorModel('nullcell.dat',mesh,actv)
-#          
-# # Load in observation file
-# [B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
-# 
-# rxLoc = dobs[:,0:3]
-# #rxLoc[:,2] += 5 # Temporary change for test
-# ndata = rxLoc.shape[0]
-#==============================================================================
-
-
-#%% Run forward modeling
-# Compute forward model using integral equation
-#==============================================================================
-# d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
-# 
-# # Form data object with coordinates and write to file
-# wd =  np.zeros((ndata,1))
-# 
-# # Save forward data to file
-# PF.Magnetics.writeUBCobs(home_dir + dsep + 'FWR_data.dat',B,M,rxLoc,d,wd)
-#==============================================================================
+#%% First test, just brake the intersecting cells
 
+ii = 1
+ddata = 9999.
 
+while ddata > 2.:
+    
+    ii += 1
+    
+    indx = np.unravel_index(bc,(mesh.nCx,mesh.nCy,mesh.nCz), order = 'F')
+    
+    xbreak = np.unique(indx[0])
+    ybreak = np.unique(indx[1])
+    zbreak = np.unique(indx[2])
+    
+    
+    # Compute the new distance
+    dl = np.sum(hxind[xbreak])
+    dx = float(hxind[xbreak].min()/2)
+    nx = dl/dx
+    
+    hxind = np.r_[hxind[0:xbreak.min()],np.ones(nx)*dx,hxind[xbreak.max():]]
+    
+    dl = np.sum(hyind[ybreak])
+    dy = float(hyind[ybreak].min()/2)
+    ny = dl/dy
+    
+    hyind = np.r_[hyind[0:ybreak.min()],np.ones(ny)*dy,hyind[ybreak.max():]]
+    
+    dl = np.sum(hzind[zbreak])
+    dz = float(hzind[zbreak].min()/2)
+    nz = dl/dz
+    
+    hzind = np.r_[hzind[0:zbreak.min()],np.ones(nz)*dz,hzind[zbreak.max():]]
+    
+    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
+        
+    # Load GOCAD surf
+    #[vrtx, trgl] = PF.BaseMag.read_GOCAD_ts(tsfile)
+    # Find active cells from surface
+    tin = tm.time()
+    print "Computing indices with VTK: "
+    [indx, bc] = PF.BaseMag.gocad2vtk(tsfile,mesh, bcflag = False, inflag = True)
+    print "VTK operation completed in " + str(tm.time() - tin)
+    
+    actv = np.zeros(mesh.nC)
+    actv[indx] = 1
+    
+    model= np.zeros(mesh.nC)
+    model[indx]= chibkg
+    
+    Utils.meshutils.writeUBCTensorModel('VTKout' + str(ii) + '.dat',mesh,model)
+    
+    Utils.meshutils.writeUBCTensorMesh('Mesh_temp' + str(ii) + '.msh',mesh)
+    
+    start_time = tm.time()
+        
+    d_temp = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
+    
+    timer = (tm.time() - start_time)
+    
+    ddata = np.max(abs(d-d_temp))
+    
+    d = d_temp
+    
+    #%% Plot data
+    plt.figure()
+    ax = plt.subplot()
+    plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
+    plt.clim(0,25)    
+    plt.colorbar(fraction=0.02)
+    plt.contour(X,Y, np.reshape(d,X.shape),10)
+    plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
+    
+    ax.set_title('Forward data')

simpegPF/Dev/Intgrl_MAG_Inv_Driver.py

@@ -65,15 +65,29 @@ if magfile=='DEFAULT':
 else:
     M_xyz = np.genfromtxt(magfile,delimiter=' \n',dtype=np.str,comments='!')
 
+# Get index of the center
+midx = int(mesh.nCx/2)
+midy = int(mesh.nCy/2)
 
 # Create forward operator
 F = PF.Magnetics.Intrgl_Fwr_Op(mesh,B,M_xyz,rxLoc,actv,'tmi')
 
 # Get distance weighting function
 wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,3.,np.min(mesh.hx)/4)
+wrMap = PF.BaseMag.WeightMap(mesh, wr)
+
 Utils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr)
 
 # Write out the predicted
 pred = F.dot(mstart)
 PF.Magnetics.writeUBCobs(home_dir + dsep + 'Pred.dat',B,M,rxLoc,pred,wd)
 
+#%%
+plt.figure()
+ax = plt.subplot()
+mesh.plotSlice(wr, ax = ax, normal = 'Y', ind=midx)
+plt.title('Distance weighting')
+plt.xlabel('x');plt.ylabel('z')
+plt.gca().set_aspect('equal', adjustable='box')
+
+

simpegPF/Dev/Test_tria_2_grid/Mesh_temp2.msh

@@ -0,0 +1,5 @@
+41 21 15
+-105.00 -55.00 0.00
+5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00 
+5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00 
+10.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00 10.00 10.00 10.00 10.00 10.00 

simpegPF/Dev/Test_tria_2_grid/Mesh_temp3.msh

@@ -0,0 +1,5 @@
+85 45 26
+-110.00 -60.00 0.00
+2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 10.00 
+2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 10.00 
+10.00 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 5.00 10.00 10.00 10.00 10.00 10.00 10.00 

simpegPF/Dev/Test_tria_2_grid/Mesh_temp4.msh

@@ -0,0 +1,5 @@
+20 10 45
+-100.00 -50.00 0.00
+10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 
+10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 
+10.00 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 2.50 2.50 2.50 2.50 5.00 5.00 10.00 10.00 10.00 10.00 10.00 10.00