wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-02 08:30:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

fixed interpmat need more testing

Browse Source
This commit is contained in:
ehaber99
2013-07-29 12:47:41 -07:00
parent f7713656ef
commit 1fca1b894e
1 changed files with 29 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/interpmat.py
+29 -23
View File
@@ -1,6 +1,4 @@
from scipy import sparse as sp
from sputils import sdiag
from utils import sub2ind, ndgrid, mkvc
import numpy as np
def interpmat(x,y,z,xr,yr,zr):
@@ -12,17 +10,20 @@ def interpmat(x,y,z,xr,yr,zr):
# Interpolation matrix
#
nx = size(x)
ny = size(y)
nz = size(z)
nx = np.size(x)
ny = np.size(y)
nz = np.size(z)
np = size(xr)
nps = np.size(xr)
#Q = spalloc(np,nx*ny*nz,8*np);
Q = sparse.coo_matrix((0.0,(0,0)),shape=(nx*ny*nz,8*np))
for i in range(0, np):
im = amin(abs(xr[i]-x))
Q = sp.lil_matrix((nps,nx*ny*nz))
ind_x = np.array([0,0])
ind_y = np.array([0,0])
ind_z = np.array([0,0])
dx, dy, dz = np.zeros(2), np.zeros(2), np.zeros(2)
for i in range(0, nps):
im = np.amin(abs(xr[i]-x))
if xr[i] - x[im] >= 0: # Point on the left
ind_x[0] = im; ind_x[1] = im+1
else: # Point on the right
@@ -32,7 +33,7 @@ def interpmat(x,y,z,xr,yr,zr):
dx[0] = xr[i] - x[ind_x[0]]
dx[1] = x[ind_x[1]] - xr[i]
im = amin(abs(yr[i] - y))
im = np.amin(abs(yr[i] - y))
if yr[i] - y[im] >= 0: # Point on the left
ind_y[0] = im; ind_y[1] = im+1
else: # Point on the right
@@ -42,8 +43,8 @@ def interpmat(x,y,z,xr,yr,zr):
dy[0] = yr[i] - y[ind_y[0]]
dy[1] = y[ind_y[1]] - yr[i];
im = amin(abs(zr[i] - z));
if zr(i) -z(im) >= 0: # Point on the left
im = np.amin(abs(zr[i] - z));
if zr[i] -z[im] >= 0: # Point on the left
ind_z[0] = im; ind_z[1] = im+1
else: # Point on the right
ind_z[0] = im-1; ind_z[1] = im;
@@ -53,27 +54,32 @@ def interpmat(x,y,z,xr,yr,zr):
Dx = x[ind_x[1]] - x[ind_x[0]]
Dy = y[ind_y[1]] - y[ind_y[0]]
Dz = z[ind_z[1]] - z[ind_z[0]]
dv = Dx*Dy*Dz
#dv = Dx*Dy*Dz
# Get the row in the matrix
v = zeros([nx, ny,nz]);
v = np.zeros([nx, ny,nz])
v[ ind_x[0], ind_y[0], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
v[ ind_x[0], ind_y[1], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz);
v[ ind_x[1], ind_y[0], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz);
v[ ind_x[1], ind_y[1], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz);
v[ ind_x[0], ind_y[0], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz);
v[ ind_x[0], ind_y[1], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz);
v[ ind_x[1], ind_y[0], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz);
v[ ind_x[1], ind_y[1], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz);
v[ ind_x[0], ind_y[1], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
v[ ind_x[1], ind_y[0], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
v[ ind_x[1], ind_y[1], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
v[ ind_x[0], ind_y[0], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
v[ ind_x[0], ind_y[1], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
v[ ind_x[1], ind_y[0], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
v[ ind_x[1], ind_y[1], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
print(np.shape(v.flatten('F')))
print(np.shape(Q))
Q[i,:] = v.flatten('F')
return Q
return Q.tocsr()
if __name__ == '__main__':
x = np.array([1, 2, 3, 4])
y = np.array([1, 2, 3, 4, 5])
z = np.array([0, 1, 4, 6])