diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index d4b2e9ad..1ff5cea4 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -169,14 +169,14 @@ def getEdgeInnerProduct(mesh, sigma): # Cell volume v = np.sqrt(mesh.vol) v3 = np.r_[v, v, v] - #V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry + V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry - #A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111 + A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111 - #A = 0.125*A + A = 0.125*A P = sp.vstack((sdiag(v3)*P000,sdiag(v3)*P001,sdiag(v3)*P010,sdiag(v3)*P011, sdiag(v3)*P100,sdiag(v3)*P101,sdiag(v3)*P110,sdiag(v3)*P111)) - A = 0.125* (P.T*Sigma*P) + return A, P @@ -186,4 +186,6 @@ if __name__ == '__main__': h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] mesh = TensorMesh(h) mu = np.ones((mesh.nC, 6)) - A = mesh.getFaceInnerProduct(mu) + A = getFaceInnerProduct(mesh,mu) + B = getEdgeInnerProduct(mesh,mu) + diff --git a/SimPEG/interpmat.py b/SimPEG/interpmat.py index 68b6c4c1..bfb80291 100644 --- a/SimPEG/interpmat.py +++ b/SimPEG/interpmat.py @@ -1,6 +1,7 @@ from scipy import sparse as sp import numpy as np + def interpmat(x,y,z,xr,yr,zr): # # This function does local linear interpolation @@ -23,7 +24,8 @@ def interpmat(x,y,z,xr,yr,zr): 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)) + im = np.argmin(abs(xr[i]-x)) + print i,im if xr[i] - x[im] >= 0: # Point on the left ind_x[0] = im; ind_x[1] = im+1 else: # Point on the right @@ -33,7 +35,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 = np.amin(abs(yr[i] - y)) + im = np.argmin(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 @@ -43,7 +45,7 @@ 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 = np.amin(abs(zr[i] - z)); + im = np.argmin(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 @@ -80,9 +82,9 @@ def interpmat(x,y,z,xr,yr,zr): 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]) + x = np.array([1.1, 2.1, 3.6, 4.9]) + y = np.array([1.2, 2.2, 3.3, 4.9, 5.6]) + z = np.array([0.8, 1.7, 4.9, 6.5]) xr = np.array([2.5,3.2]) yr = np.array([2.4,3.6])