diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 5d8356bd..df321e5f 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/InnerProducts.py b/SimPEG/InnerProducts.py new file mode 100644 index 00000000..bc482de4 --- /dev/null +++ b/SimPEG/InnerProducts.py @@ -0,0 +1,187 @@ +from scipy import sparse as sp +from sputils import sdiag +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]) + 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 | 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]]) + 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 = 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/interpmat.py b/SimPEG/interpmat.py new file mode 100644 index 00000000..68b6c4c1 --- /dev/null +++ b/SimPEG/interpmat.py @@ -0,0 +1,92 @@ +from scipy import sparse as sp +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 = np.size(x) + ny = np.size(y) + nz = np.size(z) + + nps = np.size(xr) + + #Q = spalloc(np,nx*ny*nz,8*np); + 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 + 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 = 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 + 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 = 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; + + 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 = 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) + + + print(np.shape(v.flatten('F'))) + print(np.shape(Q)) + + Q[i,:] = v.flatten('F') + + + 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]) + + 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 diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 7f9f97ca..55ac7d34 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 + tolerance = 0.85 meshSizes = [4, 8, 16, 32] + 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() diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 0aae0615..c304dda5 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,15 +1,36 @@ import numpy as np import unittest from OrderTest import OrderTest -import sys -sys.path.append('../') -from getEdgeInnerProducts import * -class TestEdgeInnerProduct(OrderTest): - """Integrate an edge function over a unit cube domain using edgeInnerProducts.""" +# MATLAB code: - name = "Edge Inner Product" +# 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.""" def getError(self): @@ -26,22 +47,70 @@ 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 = getEdgeInnerProduct(self.M, 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.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() diff --git a/SimPEG/utils.py b/SimPEG/utils.py index 9847751f..c4280457 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -1,5 +1,4 @@ import numpy as np -from numpy import * def mkvc(x, numDims=1): @@ -283,11 +282,6 @@ def faceInfo(xyz, A, B, C, D, average=True): return N, area, edgeLengths -def getSubArray(A, ind): - """subArray""" - return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]] - - def ind2sub(shape, ind): """From the given shape, returns the subscrips of the given index""" revshp = [] @@ -295,13 +289,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 @@ -311,7 +305,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)) + idx = np.dot((subs), (mult)) return idx + + +def getSubArray(A, ind): + """subArray""" + return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]