From e073eaeb8b35f5fdb205b13edc92aae191dc8d8f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 5 Aug 2013 11:51:03 -0700 Subject: [PATCH] Innerproduct work. Testing, visualization, and 2D - NOTE there are bugs in anything but a uniform LOM Not very helpful yet! --- SimPEG/InnerProducts.py | 94 ++++++++++++++++++- SimPEG/LomView.py | 32 ++++++- SimPEG/sputils.py | 45 +++++++++- SimPEG/tests/OrderTest.py | 61 +++++++------ SimPEG/tests/test_massMatrices.py | 90 ++++++++++++++++++- SimPEG/tests/test_operators.py | 144 ++++++++++++++++++++++++++++++ SimPEG/tests/test_tensorMesh.py | 80 +---------------- SimPEG/tests/test_utils.py | 25 ++++++ SimPEG/utils.py | 10 ++- 9 files changed, 465 insertions(+), 116 deletions(-) create mode 100644 SimPEG/tests/test_operators.py diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index d1d0d623..5cc391c4 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from sputils import sdiag, inv3X3BlockDiagonal +from sputils import sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal from utils import sub2ind, ndgrid, mkvc, getSubArray import numpy as np @@ -16,7 +16,10 @@ class InnerProducts(object): pass elif self._meshType == 'LOM': pass # todo: we should be doing something slightly different here! - return getFaceInnerProduct(self, mu, returnP) + if self.dim == 2: + return getFaceInnerProduct2D(self, mu, returnP) + elif self.dim == 3: + return getFaceInnerProduct(self, mu, returnP) def getEdgeInnerProduct(self, sigma=None, returnP=False): if self._meshType == 'TENSOR': @@ -59,6 +62,11 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + if mesh._meshType == 'LOM': + fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M') + fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M') + fN3 = mesh.r(mesh.normals, 'F', 'Fz', 'M') + 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] @@ -66,7 +74,15 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): 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() + PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr() + + if mesh._meshType == 'LOM': + I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), + getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), + getSubArray(fN3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]])) + PXXX = I3x3 * PXXX + + return PXXX # no | node | f1 | f2 | f3 # 000 | i ,j ,k | i , j, k | i, j , k | i, j, k @@ -110,6 +126,78 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): return A +def getFaceInnerProduct2D(mesh, mu=None, returnP=False): + + if mu is None: # default is ones + mu = np.ones((mesh.nC, 1)) + + m = np.array([mesh.nCx, mesh.nCy]) + nc = mesh.nC + + i, j = np.int64(range(m[0])), np.int64(range(m[1])) + + iijj = ndgrid(i, j) + ii, jj = iijj[:, 0], iijj[:, 1] + + if mesh._meshType == 'LOM': + fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M') + fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M') + + def Pxx(pos): + ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]]) + ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0] + + IND = np.r_[ind1, ind2].flatten() + + PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nF))).tocsr() + + if mesh._meshType == 'LOM': + # print fN1[0].shape + # print fN2[0].shape + # print np.c_[i+pos[0][0],j+pos[0][1],i+pos[1][0],j+pos[1][1]] + # print fN1[1].shape + I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1]]), + getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1]])) + PXX = I2x2 * PXX + + # import matplotlib.pyplot as plt + # plt.spy(PXX) + # plt.show() + return PXX + + # no | node | f1 | f2 + # 00 | i ,j | i , j | i, j + # 10 | i+1,j | i+1, j | i, j + # 01 | i ,j+1 | i , j | i, j+1 + # 11 | i+1,j+1 | i+1, j | i, j+1 + + # Square root of cell volume multiplied by 1/4 + v = np.sqrt(0.25*mesh.vol) + V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry + + P00 = V2*Pxx([[0, 0], [0, 0]]) + P10 = V2*Pxx([[1, 0], [0, 0]]) + P01 = V2*Pxx([[0, 0], [0, 1]]) + P11 = V2*Pxx([[1, 0], [0, 1]]) + + if mu.size == mesh.nC: # Isotropic! + mu = mkvc(mu) # ensure it is a vector. + Mu = sdiag(np.r_[mu, mu]) + elif mu.shape[1] == 2: # Diagonal tensor + Mu = sdiag(np.r_[mu[:, 0], mu[:, 1]]) + elif mu.shape[1] == 3: # Fully anisotropic + row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 2]))) + row2 = sp.hstack((sdiag(mu[:, 2]), sdiag(mu[:, 1]))) + Mu = sp.vstack((row1, row2)) + + A = P00.T*Mu*P00 + P10.T*Mu*P10 + P01.T*Mu*P01 + P11.T*Mu*P11 + P = [P00, P10, P01, P11] + if returnP: + return A, P + else: + return A + + def getEdgeInnerProduct(mesh, sigma=None, returnP=False): if sigma is None: # default is ones diff --git a/SimPEG/LomView.py b/SimPEG/LomView.py index 1596dfd9..4c9b1dab 100644 --- a/SimPEG/LomView.py +++ b/SimPEG/LomView.py @@ -14,7 +14,7 @@ class LomView(object): def __init__(self): pass - def plotGrid(self): + def plotGrid(self, length=0.05): """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.""" NN = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: @@ -31,6 +31,36 @@ class LomView(object): Y = np.r_[Y1, Y2] plt.plot(X, Y) + + plt.hold(True) + Nx = self.r(self.normals, 'F', 'Fx', 'V') + Ny = self.r(self.normals, 'F', 'Fy', 'V') + Tx = self.r(self.tangents, 'E', 'Ex', 'V') + Ty = self.r(self.tangents, 'E', 'Ey', 'V') + + plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') + + nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() + nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() + plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') + plt.plot(nX, nY, 'r-') + + nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() + nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() + #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') + plt.plot(nX, nY, 'g-') + + tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() + tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() + plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') + plt.plot(tX, tY, 'r-') + + nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() + nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() + #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') + plt.plot(nX, nY, 'g-') + plt.axis('equal') + elif self.dim == 3: fig = plt.figure(3) fig.clf() diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index ae41a210..79a1abf9 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -4,7 +4,7 @@ from utils import mkvc def sdiag(h): """Sparse diagonal matrix""" - return sp.spdiags(h, 0, h.size, h.size, format="csr") + return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr") def speye(n): @@ -23,6 +23,16 @@ def spzeros(n1, n2): def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33): + """ B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33) + + inverts a stack of 3x3 matrices + + Input: + A - a11, a12, a13, a21, a22, a23, a31, a32, a33 + + Output: + B - inverse + """ a11 = mkvc(a11) a12 = mkvc(a12) @@ -53,3 +63,36 @@ def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33): sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33))))) return B + + +def inv2X2BlockDiagonal(a11, a12, a21, a22): + """ B = inv2X2BlockDiagonal(a11, a12, a21, a22) + + Inverts a stack of 2x2 matrices by using the inversion formula + + inv(A) = (1/det(A)) * cof(A)^T + + Input: + A - a11, a12, a13, a21, a22, a23, a31, a32, a33 + + Output: + B - inverse + """ + + a11 = mkvc(a11) + a12 = mkvc(a12) + a21 = mkvc(a21) + a22 = mkvc(a22) + + # compute inverse of the determinant. + detAinv = 1./(a11*a22 - a21*a12) + + b11 = +detAinv*a22 + b12 = -detAinv*a12 + b21 = -detAinv*a21 + b22 = +detAinv*a11 + + B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))), + sp.hstack((sdiag(b21), sdiag(b22))))) + + return B diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 64d37a32..107a8d9f 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -65,20 +65,21 @@ class OrderTest(unittest.TestCase): expectedOrder = 2 tolerance = 0.85 meshSizes = [4, 8, 16, 32] - meshType = 'uniformTensorMesh' + meshTypes = ['uniformTensorMesh'] + _meshType = meshTypes[0] 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. """ - if 'TensorMesh' in self.meshType: - if 'uniform' in self.meshType: + 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: + elif 'random' in self._meshType: h1 = np.random.rand(nc) h2 = np.random.rand(nc) h3 = np.random.rand(nc) @@ -90,10 +91,10 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h - elif 'LOM' in self.meshType: - if 'uniform' in self.meshType: + elif 'LOM' in self._meshType: + if 'uniform' in self._meshType: kwrd = 'rect' - elif 'rotate' in self.meshType: + elif 'rotate' in self._meshType: kwrd = 'rotate' else: raise Exception('Unexpected meshType') @@ -117,28 +118,30 @@ class OrderTest(unittest.TestCase): """ - order = [] - err_old = 0. - max_h_old = 0. - for ii, nc in enumerate(self.meshSizes): - max_h = self.setupMesh(nc) - err = self.getError() - if ii == 0: - print '' - print 'Testing convergence on ' + self.M._meshType + ' of: ' + self.name - print '_____________________________________________' - print ' h | error | e(i-1)/e(i) | order' - print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' - print '%4i | %8.2e |' % (nc, err) - else: - 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 - max_h_old = max_h - print '---------------------------------------------' - 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) + for meshType in self.meshTypes: + self._meshType = meshType + order = [] + err_old = 0. + max_h_old = 0. + for ii, nc in enumerate(self.meshSizes): + max_h = self.setupMesh(nc) + err = self.getError() + if ii == 0: + print '' + print 'Testing convergence on ' + self.M._meshType + ' of: ' + self.name + print '_____________________________________________' + print ' h | error | e(i-1)/e(i) | order' + print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' + print '%4i | %8.2e |' % (nc, err) + else: + 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 + max_h_old = max_h + print '---------------------------------------------' + 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 021a307b..83d8bf3e 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -32,9 +32,9 @@ from OrderTest import OrderTest class TestInnerProducts(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshType = 'rotateLOM' - meshDimension = 2 - meshSizes = [16, 32, 64] + meshTypes = ['uniformTensorMesh', 'uniformLOM'] + meshDimension = 3 + meshSizes = [16, 32] def getError(self): @@ -118,5 +118,89 @@ class TestInnerProducts(OrderTest): self.orderTest() +class TestInnerProducts2D(OrderTest): + """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" + + meshTypes = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] + meshDimension = 2 + meshSizes = [4, 8, 16, 32, 64, 128] + + def getError(self): + + z = 5 # Because 5 is just such a great number. + + call = lambda fun, xy: fun(xy[:, 0], xy[:, 1]) + + ex = lambda x, y: x**2+y*z + ey = lambda x, y: (z**2)*x+y*z + + sigma1 = lambda x, y: x*y+1 + sigma2 = lambda x, y: x*z+2 + sigma3 = lambda x, y: 3+z*y + + Gc = self.M.gridCC + if self.sigmaTest == 1: + sigma = np.c_[call(sigma1, Gc)] + analytic = 144877./360 # Found using matlab symbolic toolbox. z=5 + elif self.sigmaTest == 2: + sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)] + analytic = 189959./120 # Found using matlab symbolic toolbox. z=5 + elif self.sigmaTest == 3: + sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + analytic = 781427./360 # Found using matlab symbolic toolbox. z=5 + + if self.location == 'edges': + Ex = call(ex, self.M.gridEx) + Ey = call(ey, self.M.gridEy) + E = np.matrix(np.r_[Ex, Ey]).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) + F = np.matrix(np.r_[Fx, Fy]).T + A = self.M.getFaceInnerProduct(sigma) + numeric = F.T*A*F + + err = np.abs(numeric - analytic) + return err + + # def test_order1_edges(self): + # self.name = "2D Edge Inner Product - Isotropic" + # self.location = 'edges' + # self.sigmaTest = 1 + # self.orderTest() + + # def test_order3_edges(self): + # self.name = "2D Edge Inner Product - Anisotropic" + # self.location = 'edges' + # self.sigmaTest = 2 + # self.orderTest() + + # def test_order6_edges(self): + # self.name = "2D Edge Inner Product - Full Tensor" + # self.location = 'edges' + # self.sigmaTest = 3 + # self.orderTest() + + def test_order1_faces(self): + self.name = "2D Face Inner Product - Isotropic" + self.location = 'faces' + self.sigmaTest = 1 + self.orderTest() + + def test_order3_faces(self): + self.name = "2D Face Inner Product - Anisotropic" + self.location = 'faces' + self.sigmaTest = 2 + self.orderTest() + + def test_order6_faces(self): + self.name = "2D Face Inner Product - Full Tensor" + self.location = 'faces' + self.sigmaTest = 3 + self.orderTest() + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py new file mode 100644 index 00000000..03ea36d1 --- /dev/null +++ b/SimPEG/tests/test_operators.py @@ -0,0 +1,144 @@ +import numpy as np +import unittest +import sys +sys.path.append('../') +from OrderTest import OrderTest + +MESHTYPES = ['uniformTensorMesh', 'uniformLOM'] # , 'rotateLOM' + + +class TestCurl(OrderTest): + name = "Curl" + meshTypes = MESHTYPES + + def getError(self): + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + E = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + curlE_anal = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + CURL = self.M.edgeCurl + + curlE = CURL*E + err = np.linalg.norm((curlE-curlE_anal), np.inf) + return err + + def test_order(self): + self.orderTest() + + +class TestFaceDiv(OrderTest): + name = "Face Divergence" + meshTypes = MESHTYPES + + def getError(self): + DIV = self.M.faceDiv + + #Test function + fun = lambda x: np.sin(x) + Fx = fun(self.M.gridFx[:, 0]) + Fy = fun(self.M.gridFy[:, 1]) + Fz = fun(self.M.gridFz[:, 2]) + + F = np.concatenate((Fx, Fy, Fz)) + divF = DIV*F + sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) + + err = np.linalg.norm((divF-divF_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestFaceDiv2D(OrderTest): + name = "Face Divergence 2D" + meshTypes = MESHTYPES + meshDimension = 2 + + def getError(self): + DIV = self.M.faceDiv + + #Test function + fun = lambda x: np.sin(x) + Fx = fun(self.M.gridFx[:, 0]) + Fy = fun(self.M.gridFy[:, 1]) + + F = np.concatenate((Fx, Fy)) + divF = DIV*F + sol = lambda x, y: (np.cos(x)+np.cos(y)) + divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1]) + + err = np.linalg.norm((divF-divF_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad(OrderTest): + name = "Nodal Gradient" + meshTypes = MESHTYPES + + def getError(self): + GRAD = self.M.nodalGrad + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) + gradE = GRAD*phi + + Ex = sol(self.M.gridEx[:, 0]) + Ey = sol(self.M.gridEy[:, 1]) + Ez = sol(self.M.gridEz[:, 2]) + + gradE_anal = np.concatenate((Ex, Ey, Ez)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad2D(OrderTest): + name = "Nodal Gradient 2D" + meshTypes = MESHTYPES + meshDimension = 2 + + def getError(self): + GRAD = self.M.nodalGrad + #Test function + fun = lambda x, y: (np.cos(x)+np.cos(y)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1]) + gradE = GRAD*phi + + Ex = sol(self.M.gridEx[:, 0]) + Ey = sol(self.M.gridEy[:, 1]) + + gradE_anal = np.concatenate((Ex, Ey)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index bc034e0b..32d47353 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -58,84 +58,6 @@ class BasicTensorMeshTests(unittest.TestCase): self.assertTrue(t1) -class TestCurl(OrderTest): - name = "Curl" - - def getError(self): - fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) - sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - - Ex = fun(self.M.gridEx[:, 1]) - Ey = fun(self.M.gridEy[:, 2]) - Ez = fun(self.M.gridEz[:, 0]) - E = np.concatenate((Ex, Ey, Ez)) - - Fx = sol(self.M.gridFx[:, 2]) - Fy = sol(self.M.gridFy[:, 0]) - Fz = sol(self.M.gridFz[:, 1]) - curlE_anal = np.concatenate((Fx, Fy, Fz)) - - # Generate DIV matrix - CURL = self.M.edgeCurl - - curlE = CURL*E - err = np.linalg.norm((curlE-curlE_anal), np.inf) - return err - - def test_order(self): - self.orderTest() - - -class TestFaceDiv(OrderTest): - name = "Face Divergence" - - def getError(self): - DIV = self.M.faceDiv - - #Test function - fun = lambda x: np.sin(x) - Fx = fun(self.M.gridFx[:, 0]) - Fy = fun(self.M.gridFy[:, 1]) - Fz = fun(self.M.gridFz[:, 2]) - - F = np.concatenate((Fx, Fy, Fz)) - divF = DIV*F - sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) - - err = np.linalg.norm((divF-divF_anal), np.inf) - - return err - - def test_order(self): - self.orderTest() - - -class TestNodalGrad(OrderTest): - name = "Nodal Gradient" - - def getError(self): - GRAD = self.M.nodalGrad - #Test function - fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) - - phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) - gradE = GRAD*phi - - Ex = sol(self.M.gridEx[:, 0]) - Ey = sol(self.M.gridEy[:, 1]) - Ez = sol(self.M.gridEz[:, 2]) - - gradE_anal = np.concatenate((Ex, Ey, Ez)) - err = np.linalg.norm((gradE-gradE_anal), np.inf) - - return err - - def test_order(self): - self.orderTest() - - class TestPoissonEqn(OrderTest): name = "Poisson Equation" meshSizes = [16, 20, 24] @@ -168,5 +90,7 @@ class TestPoissonEqn(OrderTest): self.name = "Poisson Equation - Backward" self.forward = False self.orderTest() + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index 59b70147..393cd78d 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -3,6 +3,7 @@ import unittest import sys sys.path.append('../') from utils import mkvc, ndgrid, indexCube +from sputils import sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, sp class TestSequenceFunctions(unittest.TestCase): @@ -62,5 +63,29 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.all(indexCube('G', nN) == np.array([13, 14, 16, 17, 22, 23, 25, 26]))) self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23]))) + def test_invXXXBlockDiagonal(self): + + a = [np.random.rand(5, 1) for i in range(4)] + + B = inv2X2BlockDiagonal(*a) + + A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]))), + sp.hstack((sdiag(a[2]), sdiag(a[3]))))) + + Z2 = B*A - sp.eye(10, 10) + self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-12) + + a = [np.random.rand(5, 1) for i in range(9)] + B = inv3X3BlockDiagonal(*a) + + A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]), sdiag(a[2]))), + sp.hstack((sdiag(a[3]), sdiag(a[4]), sdiag(a[5]))), + sp.hstack((sdiag(a[6]), sdiag(a[7]), sdiag(a[8]))))) + + Z3 = B*A - sp.eye(15, 15) + + self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/utils.py b/SimPEG/utils.py index 0741b7bb..1b5be6bd 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -319,4 +319,12 @@ def sub2ind(shape, subs): def getSubArray(A, ind): """subArray""" - return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]] + assert type(ind) == list, "ind must be a list of vectors" + assert len(A.shape) == len(ind), "ind must have the same length as the dimension of A" + + if len(A.shape) == 2: + return A[ind[0], :][:, ind[1]] + elif len(A.shape) == 3: + return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]] + else: + raise Exception("getSubArray does not support dimension asked.")