From 05ddf93d647ad6f54949b4c716e5a791711dbf00 Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Sat, 27 Jul 2013 16:11:24 -0700 Subject: [PATCH 1/9] Added getFaceInnerProduct and getInnerProduct that puts the edge and face mass matrices together --- SimPEG/getFaceInnerProducts.py | 87 +++++++++++++++++ SimPEG/getInnerProducts.py | 165 +++++++++++++++++++++++++++++++++ 2 files changed, 252 insertions(+) create mode 100644 SimPEG/getFaceInnerProducts.py create mode 100644 SimPEG/getInnerProducts.py diff --git a/SimPEG/getFaceInnerProducts.py b/SimPEG/getFaceInnerProducts.py new file mode 100644 index 00000000..961013c8 --- /dev/null +++ b/SimPEG/getFaceInnerProducts.py @@ -0,0 +1,87 @@ +from scipy import sparse as sp +from sputils import sdiag +from utils import sub2ind, ndgrid, mkvc +import numpy as np + + +def getFaceInnerProduct(mesh, mu): + + m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) + nc = mesh.nC + + i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + + def Pxxx(pos): + ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] + ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] + + IND = np.r_[ind1, ind2, ind3].flatten() + + return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr() + + # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) + # / / + # / / | + # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) + # / / | + # / / | + # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) + # | | | + # | | node(i+1,j+1,k+1) + # | | / + # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) + # | | / + # | | / + # | |/ + # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + + # no | node | e1 | e2 | e3 + # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k + # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k + # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k + # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k + # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 + # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 + # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 + # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 + P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) + P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) + P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) + P110 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]]) + P001 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]]) + P101 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]]) + P011 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]]) + P111 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + + if mu.size == mesh.nC: # Isotropic! + mu = mkvc(mu) # ensure it is a vector. + mu = sdiag(np.r_[mu, mu, mu]) + elif mu.shape[1] == 3: # Diagonal tensor + mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]]) + elif mu.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4]))) + row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5]))) + row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2]))) + mu = sp.vstack((row1, row2, row3)) + + # Cell volume + v = np.sqrt(mesh.vol) + v3 = np.r_[v, v, v] + V = sdiag(v3)*mu*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 = 0.125*A + + return A + +if __name__ == '__main__': + from TensorMesh import TensorMesh + 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 = getFaceInnerProduct(mesh, mu) diff --git a/SimPEG/getInnerProducts.py b/SimPEG/getInnerProducts.py new file mode 100644 index 00000000..94c02826 --- /dev/null +++ b/SimPEG/getInnerProducts.py @@ -0,0 +1,165 @@ +from scipy import sparse as sp +from sputils import sdiag +from utils import sub2ind, ndgrid, mkvc +import numpy as np + + +def getFaceInnerProduct(mesh, mu): + + m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) + nc = mesh.nC + + i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + + def Pxxx(pos): + ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] + ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] + + IND = np.r_[ind1, ind2, ind3].flatten() + + return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr() + + # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) + # / / + # / / | + # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) + # / / | + # / / | + # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) + # | | | + # | | node(i+1,j+1,k+1) + # | | / + # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) + # | | / + # | | / + # | |/ + # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + + # no | node | e1 | e2 | e3 + # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k + # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k + # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k + # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k + # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 + # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 + # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 + # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 + P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) + P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) + P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) + P110 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]]) + P001 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]]) + P101 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]]) + P011 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]]) + P111 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + + if mu.size == mesh.nC: # Isotropic! + mu = mkvc(mu) # ensure it is a vector. + mu = sdiag(np.r_[mu, mu, mu]) + elif mu.shape[1] == 3: # Diagonal tensor + mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]]) + elif mu.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4]))) + row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5]))) + row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2]))) + mu = sp.vstack((row1, row2, row3)) + + # Cell volume + v = np.sqrt(mesh.vol) + v3 = np.r_[v, v, v] + V = sdiag(v3)*mu*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 = 0.125*A + + return A + + +def getEdgeInnerProduct(mesh, sigma): + + m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) + nc = mesh.nC + + i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + + def Pxxx(pos): + ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0] + ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] + + IND = np.r_[ind1, ind2, ind3].flatten() + + return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr() + + # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) + # / / + # / / | + # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) + # / / | + # / / | + # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) + # | | | + # | | node(i+1,j+1,k+1) + # | | / + # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) + # | | / + # | | / + # | |/ + # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + + # no | node | e1 | e2 | e3 + # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k + # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k + # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k + # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k + # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k + # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k + # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k + # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k + P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) + P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]]) + P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]]) + P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]]) + P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]]) + P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]]) + P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) + P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) + + if sigma.size == mesh.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma, sigma]) + elif sigma.shape[1] == 3: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) + elif sigma.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) + row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) + row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) + Sigma = sp.vstack((row1, row2, row3)) + + # Cell volume + v = np.sqrt(mesh.vol) + v3 = np.r_[v, v, v] + 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 = 0.125*A + + return A + + + +if __name__ == '__main__': + from TensorMesh import TensorMesh + 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 = getFaceInnerProduct(mesh, mu) From f7713656ef47bfc8a0363358cfa6fee7e4a8ea15 Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Sat, 27 Jul 2013 17:36:00 -0700 Subject: [PATCH 2/9] Added interpolation matrix. Not sure its working since my canopy is playing tricks --- SimPEG/interpmat.py | 86 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 SimPEG/interpmat.py diff --git a/SimPEG/interpmat.py b/SimPEG/interpmat.py new file mode 100644 index 00000000..35e998ff --- /dev/null +++ b/SimPEG/interpmat.py @@ -0,0 +1,86 @@ +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): +# +# This function does local linear interpolation +# computed for each receiver point in turn +# +# [Q] = linint(x,y,z,xr,yr,zr) +# Interpolation matrix +# + + nx = size(x) + ny = size(y) + nz = size(z) + + 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)) + if xr[i] - x[im] >= 0: # Point on the left + ind_x[0] = im; ind_x[1] = im+1 + else: # Point on the right + ind_x[0] = im-1; ind_x[1] = im + + + dx[0] = xr[i] - x[ind_x[0]] + dx[1] = x[ind_x[1]] - xr[i] + + im = 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 + ind_y[0] = im-1; ind_y[1] = im + + + 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 + ind_z[0] = im; ind_z[1] = im+1 + else: # Point on the right + ind_z[0] = im-1; ind_z[1] = im; + + dz[0] = zr[i] - z[ind_z[0]]; dz[1] = z[ind_z[1]] - zr[i] + + 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 + + # Get the row in the matrix + v = 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); + + + Q[i,:] = v.flatten('F') + + return Q + + +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]) + + xr = np.array([2.5,3.2]) + yr = np.array([2.4,3.6]) + zr = np.array([2.5,3.9]) + + A = interpmat(x,y,z,xr,yr,zr) + \ No newline at end of file From 1fca1b894e6e165e8e6309916bc9daabdad3c1a0 Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Mon, 29 Jul 2013 12:47:41 -0700 Subject: [PATCH 3/9] fixed interpmat need more testing --- SimPEG/interpmat.py | 52 +++++++++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/SimPEG/interpmat.py b/SimPEG/interpmat.py index 35e998ff..68b6c4c1 100644 --- a/SimPEG/interpmat.py +++ b/SimPEG/interpmat.py @@ -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]) From d8f646c7d99216fe860daf861794c6ddb39afb15 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 21:49:23 -0700 Subject: [PATCH 4/9] Updates to OrderTest so that it runs random meshes. --- SimPEG/tests/OrderTest.py | 42 +++++++++++++++++++++++++++------------ 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 7f9f97ca..59dce6d6 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -20,7 +20,7 @@ class OrderTest(unittest.TestCase): Note that you can provide any norm. - Test is passed when estimated rate order of convergence is at least 90% of the + Test is passed when estimated rate order of convergence is at least within the specified tolerance of the estimated rate supplied by the user. Minimal example for a curl operator: @@ -63,17 +63,32 @@ class OrderTest(unittest.TestCase): name = "Order Test" expectedOrder = 2 - meshSizes = [4, 8, 16, 32] + tolerance = 0.85 + meshSizes = [4, 8, 16, 32, 64] + meshType = 'uniformTensorMesh' + meshDimension = 3 def setupMesh(self, nc): """ For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc. """ - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] - self.M = TensorMesh(h) + if 'TensorMesh' in self.meshType: + if 'uniform' in self.meshType: + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc + h = [h1, h2, h3] + elif 'random' in self.meshType: + h1 = np.random.rand(nc) + h2 = np.random.rand(nc) + h3 = np.random.rand(nc) + h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize + else: + raise Exception('Unexpected meshType') + + self.M = TensorMesh(h[:self.meshDimension]) + max_h = max([np.max(hi) for hi in self.M.h]) + return max_h def getError(self): """For given h, generate A[h], f and A(f) and return norm of error.""" @@ -89,9 +104,9 @@ class OrderTest(unittest.TestCase): """ order = [] err_old = 0. - nc_old = 0. + max_h_old = 0. for ii, nc in enumerate(self.meshSizes): - self.setupMesh(nc) + max_h = self.setupMesh(nc) err = self.getError() if ii == 0: print '' @@ -101,13 +116,14 @@ class OrderTest(unittest.TestCase): print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' print '%4i | %8.2e |' % (nc, err) else: - order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc))) + order.append(np.log(err/err_old)/np.log(max_h/max_h_old)) print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) err_old = err - nc_old = nc + max_h_old = max_h print '---------------------------------------------' - self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order))) - + passTest = np.mean(np.array(order)) > self.tolerance*self.expectedOrder + # passTest = len(np.where(np.array(order) > self.tolerance*self.expectedOrder)[0]) > np.floor(0.75*len(order)) + self.assertTrue(passTest) if __name__ == '__main__': unittest.main() From c259a6651a7760ff2cf93aca20642516f48673bd Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:12:53 -0700 Subject: [PATCH 5/9] Integrated getEdge/FaceInnerProduct into the tensor mesh class. --- SimPEG/DiffOperators.py | 2 +- .../{getInnerProducts.py => InnerProducts.py} | 42 ++++++--- SimPEG/TensorMesh.py | 5 +- SimPEG/getEdgeInnerProducts.py | 87 ------------------- SimPEG/getFaceInnerProducts.py | 87 ------------------- SimPEG/tests/OrderTest.py | 2 +- SimPEG/tests/test_massMatrices.py | 5 +- 7 files changed, 39 insertions(+), 191 deletions(-) rename SimPEG/{getInnerProducts.py => InnerProducts.py} (82%) delete mode 100644 SimPEG/getEdgeInnerProducts.py delete mode 100644 SimPEG/getFaceInnerProducts.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index c5daa96d..6688f254 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -50,7 +50,7 @@ class DiffOperators(object): Class creates the differential operators that you need! """ def __init__(self): - raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') + raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') def faceDiv(): doc = "Construct divergence operator (face-stg to cell-centres)." diff --git a/SimPEG/getInnerProducts.py b/SimPEG/InnerProducts.py similarity index 82% rename from SimPEG/getInnerProducts.py rename to SimPEG/InnerProducts.py index 94c02826..bc482de4 100644 --- a/SimPEG/getInnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -4,6 +4,28 @@ from utils import sub2ind, ndgrid, mkvc import numpy as np +class InnerProducts(object): + """ + Class creates the inner product matrices that you need! + """ + def __init__(self): + raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') + + def getFaceInnerProduct(self, mu): + if self._meshType == 'TENSOR': + pass + elif self._meshType == 'LOM': + pass # todo: we should be doing something slightly different here! + return getFaceInnerProduct(self, mu) + + def getEdgeInnerProduct(self, sigma): + if self._meshType == 'TENSOR': + pass + elif self._meshType == 'LOM': + pass # todo: we should be doing something slightly different here! + return getEdgeInnerProduct(self, sigma) + + def getFaceInnerProduct(mesh, mu): m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) @@ -39,15 +61,15 @@ def getFaceInnerProduct(mesh, mu): # | |/ # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k - # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k - # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k - # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 - # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 - # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 - # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 + # no | node | f1 | f2 | f3 + # 000 | i ,j ,k | i , j, k | i, j , k | i, j, k + # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k + # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k + # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k + # 001 | i ,j ,k | i , j, k | i, j , k | i, j, k+1 + # 101 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k+1 + # 011 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k+1 + # 111 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k+1 P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) @@ -162,4 +184,4 @@ 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 = getFaceInnerProduct(mesh, mu) + A = mesh.getFaceInnerProduct(mu) diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index b247138b..acfdfbac 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -2,10 +2,11 @@ import numpy as np from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators +from InnerProducts import InnerProducts from utils import ndgrid, mkvc -class TensorMesh(BaseMesh, TensorView, DiffOperators): +class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -21,6 +22,8 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): mesh = TensorMesh([hx, hy, hz]) """ + _meshType = 'TENSOR' + def __init__(self, h, x0=None): super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py deleted file mode 100644 index 68423a3f..00000000 --- a/SimPEG/getEdgeInnerProducts.py +++ /dev/null @@ -1,87 +0,0 @@ -from scipy import sparse as sp -from sputils import sdiag -from utils import sub2ind, ndgrid, mkvc -import numpy as np - - -def getEdgeInnerProduct(mesh, sigma): - - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC - - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) - - iijjkk = ndgrid(i, j, k) - ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - - def Pxxx(pos): - ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0] - ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] - - IND = np.r_[ind1, ind2, ind3].flatten() - - return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr() - - # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) - # / / - # / / | - # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) - # / / | - # / / | - # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) - # | | | - # | | node(i+1,j+1,k+1) - # | | / - # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) - # | | / - # | | / - # | |/ - # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k - # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k - # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k - # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k - # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k - # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k - # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k - P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]]) - P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]]) - P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]]) - P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]]) - P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]]) - P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) - P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) - - if sigma.size == mesh.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.shape[1] == 3: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) - row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) - row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) - Sigma = sp.vstack((row1, row2, row3)) - - # Cell volume - v = np.sqrt(mesh.vol) - v3 = np.r_[v, v, v] - 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 = 0.125*A - - return A - -if __name__ == '__main__': - from TensorMesh import TensorMesh - h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] - mesh = TensorMesh(h) - sigma = np.ones((mesh.nC, 6)) - A = getEdgeInnerProduct(mesh, sigma) diff --git a/SimPEG/getFaceInnerProducts.py b/SimPEG/getFaceInnerProducts.py deleted file mode 100644 index 961013c8..00000000 --- a/SimPEG/getFaceInnerProducts.py +++ /dev/null @@ -1,87 +0,0 @@ -from scipy import sparse as sp -from sputils import sdiag -from utils import sub2ind, ndgrid, mkvc -import numpy as np - - -def getFaceInnerProduct(mesh, mu): - - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC - - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) - - iijjkk = ndgrid(i, j, k) - ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - - def Pxxx(pos): - ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] - ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] - - IND = np.r_[ind1, ind2, ind3].flatten() - - return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr() - - # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) - # / / - # / / | - # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) - # / / | - # / / | - # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) - # | | | - # | | node(i+1,j+1,k+1) - # | | / - # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) - # | | / - # | | / - # | |/ - # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k - # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k - # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k - # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 - # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 - # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 - # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 - P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) - P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) - P110 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]]) - P001 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]]) - P101 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]]) - P011 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]]) - P111 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) - - if mu.size == mesh.nC: # Isotropic! - mu = mkvc(mu) # ensure it is a vector. - mu = sdiag(np.r_[mu, mu, mu]) - elif mu.shape[1] == 3: # Diagonal tensor - mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]]) - elif mu.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4]))) - row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5]))) - row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2]))) - mu = sp.vstack((row1, row2, row3)) - - # Cell volume - v = np.sqrt(mesh.vol) - v3 = np.r_[v, v, v] - V = sdiag(v3)*mu*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 = 0.125*A - - return A - -if __name__ == '__main__': - from TensorMesh import TensorMesh - 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 = getFaceInnerProduct(mesh, mu) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 59dce6d6..55ac7d34 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -64,7 +64,7 @@ class OrderTest(unittest.TestCase): name = "Order Test" expectedOrder = 2 tolerance = 0.85 - meshSizes = [4, 8, 16, 32, 64] + meshSizes = [4, 8, 16, 32] meshType = 'uniformTensorMesh' meshDimension = 3 diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 0aae0615..4840abed 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,9 +1,6 @@ import numpy as np import unittest from OrderTest import OrderTest -import sys -sys.path.append('../') -from getEdgeInnerProducts import * class TestEdgeInnerProduct(OrderTest): @@ -35,7 +32,7 @@ class TestEdgeInnerProduct(OrderTest): sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)] - A = getEdgeInnerProduct(self.M, sigma) + A = self.M.getEdgeInnerProduct(sigma) numeric = E.T*A*E analytic = 69881./21600 # Found using matlab symbolic toolbox. err = np.abs(numeric - analytic) From 9573c6c230571b21f5cc8998dff4bb89efac9c29 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:40:10 -0700 Subject: [PATCH 6/9] Tested face and edge inner products for anisotropic, isotropic and tensor sigma. --- SimPEG/tests/test_massMatrices.py | 94 ++++++++++++++++++++++++++----- 1 file changed, 81 insertions(+), 13 deletions(-) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 4840abed..f8fc99b0 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -3,8 +3,34 @@ import unittest from OrderTest import OrderTest -class TestEdgeInnerProduct(OrderTest): - """Integrate an edge function over a unit cube domain using edgeInnerProducts.""" +# MATLAB code: + +# syms x y z + +# ex = x.^2+y.*z; +# ey = (z.^2).*x+y.*z; +# ez = y.^2+x.*z; + +# e = [ex;ey;ez]; + +# sigma1 = x.*y+1; +# sigma2 = x.*z+2; +# sigma3 = 3+z.*y; +# sigma4 = 0.1.*x.*y.*z; +# sigma5 = 0.2.*x.*y; +# sigma6 = 0.1.*z; + +# S1 = [sigma1,0,0;0,sigma1,0;0,0,sigma1]; +# S2 = [sigma1,0,0;0,sigma2,0;0,0,sigma3]; +# S3 = [sigma1,sigma4,sigma5;sigma4,sigma2,sigma6;sigma5,sigma6,sigma3]; + +# i1 = int(int(int(e.'*S1*e,x,0,1),y,0,1),z,0,1); +# i2 = int(int(int(e.'*S2*e,x,0,1),y,0,1),z,0,1); +# i3 = int(int(int(e.'*S3*e,x,0,1),y,0,1),z,0,1); + + +class TestInnerProducts(OrderTest): + """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" name = "Edge Inner Product" @@ -23,22 +49,64 @@ class TestEdgeInnerProduct(OrderTest): sigma5 = lambda x, y, z: 0.2*x*y sigma6 = lambda x, y, z: 0.1*z - Ex = call(ex, self.M.gridEx) - Ey = call(ey, self.M.gridEy) - Ez = call(ez, self.M.gridEz) - - E = np.matrix(np.r_[Ex, Ey, Ez]).T Gc = self.M.gridCC - sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), - call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)] + if self.sigmaTest == 1: + sigma = np.c_[call(sigma1, Gc)] + analytic = 647./360 # Found using matlab symbolic toolbox. + elif self.sigmaTest == 3: + sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + analytic = 37./12 # Found using matlab symbolic toolbox. + elif self.sigmaTest == 6: + sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), + call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)] + analytic = 69881./21600 # Found using matlab symbolic toolbox. + + if self.location == 'edges': + Ex = call(ex, self.M.gridEx) + Ey = call(ey, self.M.gridEy) + Ez = call(ez, self.M.gridEz) + E = np.matrix(np.r_[Ex, Ey, Ez]).T + A = self.M.getEdgeInnerProduct(sigma) + numeric = E.T*A*E + elif self.location == 'faces': + Fx = call(ex, self.M.gridFx) + Fy = call(ey, self.M.gridFy) + Fz = call(ez, self.M.gridFz) + F = np.matrix(np.r_[Fx, Fy, Fz]).T + A = self.M.getFaceInnerProduct(sigma) + numeric = F.T*A*F - A = self.M.getEdgeInnerProduct(sigma) - numeric = E.T*A*E - analytic = 69881./21600 # Found using matlab symbolic toolbox. err = np.abs(numeric - analytic) return err - def test_order(self): + def test_order1_edges(self): + self.location = 'edges' + self.sigmaTest = 1 + self.orderTest() + + def test_order3_edges(self): + self.location = 'edges' + self.sigmaTest = 3 + self.orderTest() + + def test_order6_edges(self): + self.location = 'edges' + self.sigmaTest = 6 + self.orderTest() + + def test_order1_faces(self): + self.location = 'faces' + self.sigmaTest = 1 + self.orderTest() + + def test_order3_faces(self): + self.location = 'faces' + self.sigmaTest = 3 + self.orderTest() + + def test_order6_faces(self): + self.location = 'faces' + self.sigmaTest = 6 self.orderTest() From 4660b177a1f0ceca9336c8263063fa21cdcea1fe Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:48:27 -0700 Subject: [PATCH 7/9] Names in tests. --- SimPEG/tests/test_massMatrices.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index f8fc99b0..c304dda5 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -32,8 +32,6 @@ from OrderTest import OrderTest class TestInnerProducts(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - name = "Edge Inner Product" - def getError(self): call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) @@ -80,31 +78,37 @@ class TestInnerProducts(OrderTest): return err def test_order1_edges(self): + self.name = "Edge Inner Product - Isotropic" self.location = 'edges' self.sigmaTest = 1 self.orderTest() def test_order3_edges(self): + self.name = "Edge Inner Product - Anisotropic" self.location = 'edges' self.sigmaTest = 3 self.orderTest() def test_order6_edges(self): + self.name = "Edge Inner Product - Full Tensor" self.location = 'edges' self.sigmaTest = 6 self.orderTest() def test_order1_faces(self): + self.name = "Face Inner Product - Isotropic" self.location = 'faces' self.sigmaTest = 1 self.orderTest() def test_order3_faces(self): + self.name = "Face Inner Product - Anisotropic" self.location = 'faces' self.sigmaTest = 3 self.orderTest() def test_order6_faces(self): + self.name = "Face Inner Product - Full Tensor" self.location = 'faces' self.sigmaTest = 6 self.orderTest() From 68e44cca5a253c9277ebd640c04ce294e772c565 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:52:37 -0700 Subject: [PATCH 8/9] Deleted massMatrices.py This code is not working and is repeated in InnerProducts.py --- SimPEG/massMatrices.py | 94 ------------------------------------------ 1 file changed, 94 deletions(-) delete mode 100644 SimPEG/massMatrices.py diff --git a/SimPEG/massMatrices.py b/SimPEG/massMatrices.py deleted file mode 100644 index ead21e3a..00000000 --- a/SimPEG/massMatrices.py +++ /dev/null @@ -1,94 +0,0 @@ -import numpy as np -from scipy import sparse as sp -from sputils import sdiag, speye, kron3, spzeros -from utils import mkvc - - - -def getEdgeMassMatrix(sigma,mesh): - """Get anisotropic mass matrix""" - - n = array([size(mesh.h[0]),size(mesh.h[1]),size(mesh.h[2])]) - nx = prod(n + [1, 0, 0]) - ex = reshape(arange(0,nx),[n[0]+1,n[1],n[2]]) - ny = prod(n + [0, 1, 0]) - ey = reshape(arange(0,ny),[n[0],n[1]+1,n[2]]) - nz = prod(n + [0, 0, 1]); - ez = reshape(arange(0,nz),[n[0],n[1],n[2]+1]) - - - i = arange(0,n[0]-1); j = arange(0,n[1]-1); k = arange(0,n[2]-1) - - # corner i,j,k - Px1 = take(ex,[i,j,k]); Py1 = take(ey,[i,j,k]); Pz1 = take(ez,[i,j,k]) - # corner i+1,j,k - Px2 = take(ex,[i,j,k]); Py2 = take(ey,[i+1,j,k]); Pz2 = take(ez,[i+1,j,k]) - # corner i,j+1,k - Px3 = take(ex,[i,j+1,k]); Py3 = take(ey,[i,j,k]); Pz3 = take(ez,[i,j+1,k]) - # corner i+1,j+1,k - Px4 = take(ex,[i,j+1,k]); Py4 = take(ey,[i+1,j,k]); Pz4 = take(ez,[i+1,j+1,k]); - - # corner i,j,k+1 - Px5 = take(ex,[i,j,k+1]); Py5 = take(ey,[i,j,k+1]); Pz5 = take(ez,[i,j,k]) - # corner i+1,j,k+1 - Px6 = take(ex,[i,j,k+1]); Py6 = take(ey,[i+1,j,k+1]); Pz6 = take(ez,[i+1,j,k]) - # corner i,j+1,k+1 - Px7 = take(ex,[i,j+1,k+1]); Py7 = take(ey,[i,j,k+1]); Pz7 = take(ez,[i,j+1,k]) - # corner i+1,j+1,k+1 - Px8 = take(ex,[i,j+1,k+1]); Py8 = take(ey,[i+1,j,k+1]); Pz8 = take(ez,[i+1,j+1,k]) - - - nx1 = size(Px1); ny1 = size(Py1); nz1 = size(Pz1) - #sparse.coo_matrix((V,(I,J)),shape=(4,4)) - P1 = block_diag(( sparse.coo_matrix(arange(0,nx1),Px1(:), e(nx1), nx1,nx), - sparse.coo_matrix(arange(0,ny1),Py1(:),e(ny1), ny1,ny), - sparse.coo_matrix(arange(0,nz1),Pz1(:),e(nz1), nz1,nz))) - - nx2 = numel(Px2); ny2 = numel(Py2); nz2 = numel(Pz2); - P2 = blkdiag( sparse(1:nx2,Px2(:), e(nx2), nx2,nx) , ... - sparse(1:ny2,Py2(:),e(ny2), ny2,ny), ... - sparse(1:nz2,Pz2(:),e(nz2), nz2,nz)); - - nx3 = numel(Px3); ny3 = numel(Py3); nz3 = numel(Pz3); - P3 = blkdiag( sparse(1:nx3,Px3(:), e(nx3), nx3,nx) , ... - sparse(1:ny3,Py3(:),e(ny3), ny3,ny), ... - sparse(1:nz3,Pz3(:),e(nz3), nz3,nz)); - - nx4 = numel(Px4); ny4 = numel(Py4); nz4 = numel(Pz4); - P4 = blkdiag( sparse(1:nx4,Px4(:), e(nx4), nx4,nx) , ... - sparse(1:ny4,Py4(:), e(ny4), ny4,ny), ... - sparse(1:nz4,Pz4(:), e(nz4), nz4,nz)); - - nx5 = numel(Px5); ny5 = numel(Py5); nz5 = numel(Pz5); - P5 = blkdiag( sparse(1:nx5,Px5(:), e(nx5), nx5,nx) , ... - sparse(1:ny5,Py5(:), e(ny5), ny5,ny), ... - sparse(1:nz5,Pz5(:), e(nz5), nz5,nz)); - - nx6 = numel(Px6); ny6 = numel(Py6); nz6 = numel(Pz6); - P6 = blkdiag( sparse(1:nx6,Px6(:), e(nx6), nx6,nx) , ... - sparse(1:ny6,Py6(:), e(ny6), ny6,ny), ... - sparse(1:nz6,Pz6(:), e(nz6), nz6,nz)); - - nx7 = numel(Px7); ny7 = numel(Py7); nz7 = numel(Pz7); - P7 = blkdiag( sparse(1:nx7,Px7(:), e(nx7), nx7,nx) , ... - sparse(1:ny7,Py7(:), e(ny7), ny7,ny), ... - sparse(1:nz7,Pz7(:), e(nz7), nz7,nz)); - - nx8 = numel(Px8); ny8 = numel(Py8); nz8 = numel(Pz8); - P8 = blkdiag( sparse(1:nx8,Px8(:), e(nx8), nx8,nx) , ... - sparse(1:ny8,Py8(:), e(ny8), ny8,ny), ... - sparse(1:nz8,Pz8(:), e(nz8), nz8,nz)); - - V = sdiag(sqrt([v(:); v(:); v(:)])); - - # generate the conductivity - S = [sdiag(sig(:,1)) , sdiag(sig(:,4)) , sdiag(sig(:,5)); ... - sdiag(sig(:,4)) , sdiag(sig(:,2)) , sdiag(sig(:,6)); ... - sdiag(sig(:,5)) , sdiag(sig(:,6)) , sdiag(sig(:,3))]; - - # scale by the volume - S = V*S*V; - - M = 1/8*(P1'*S*P1 + P2'*S*P2 + P3'*S*P3 + P4'*S*P4 + ... - P5'*S*P5 + P6'*S*P6 + P7'*S*P7 + P8'*S*P8); - \ No newline at end of file From 7e169bbb7988e28c578a8c3b7eaeb6345cf8dc6b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:56:36 -0700 Subject: [PATCH 9/9] Cleaned utils code, merged in subArray. --- SimPEG/subArray.py | 8 -------- SimPEG/utils.py | 19 ++++++++++++------- 2 files changed, 12 insertions(+), 15 deletions(-) delete mode 100644 SimPEG/subArray.py diff --git a/SimPEG/subArray.py b/SimPEG/subArray.py deleted file mode 100644 index 2811f350..00000000 --- a/SimPEG/subArray.py +++ /dev/null @@ -1,8 +0,0 @@ -import numpy as np - -def getSubArray(A,ind): - """subArray""" - i = ind[0]; j = ind[1]; k = ind[2] - - return A[i,:,:][:,j,:][:,:,k] - \ No newline at end of file diff --git a/SimPEG/utils.py b/SimPEG/utils.py index b9e17847..ec89739b 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -1,5 +1,5 @@ import numpy as np -from numpy import * + def reshapeF(x, size): return np.reshape(x, size, order='F') @@ -106,13 +106,13 @@ def ind2sub(shape, ind): mult = [1] for i in range(0, len(revshp)-1): mult.extend([mult[i]*revshp[i]]) - mult = array(mult).reshape(len(mult)) + mult = np.array(mult).reshape(len(mult)) sub = [] for i in range(0, len(shape)): - sub.extend([math.floor(ind / mult[i])]) - ind = ind - (math.floor(ind/mult[i]) * mult[i]) + sub.extend([np.math.floor(ind / mult[i])]) + ind = ind - (np.math.floor(ind/mult[i]) * mult[i]) return sub @@ -122,7 +122,12 @@ def sub2ind(shape, subs): mult = [1] for i in range(0, len(revshp)-1): mult.extend([mult[i]*revshp[i]]) - mult = array(mult).reshape(len(mult), 1) + mult = np.array(mult).reshape(len(mult), 1) - idx = dot((subs), (mult)) - return idx \ No newline at end of file + idx = np.dot((subs), (mult)) + return idx + + +def getSubArray(A, ind): + """subArray""" + return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]