From 150cbc7df3b7c94e4173629b5130776079289aaa Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 12:11:45 -0700 Subject: [PATCH] Test for edge inner products working. --- SimPEG/getEdgeInnerProducts.py | 122 +++++++++++++++--------------- SimPEG/tests/test_massMatrices.py | 47 ++++++------ 2 files changed, 86 insertions(+), 83 deletions(-) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py index 5a95f146..d3b91211 100644 --- a/SimPEG/getEdgeInnerProducts.py +++ b/SimPEG/getEdgeInnerProducts.py @@ -6,7 +6,7 @@ from numpy import * from TensorMesh import * # [A] = getEdgeInnerProduct(X,Y,Z,sigma) -# +# # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) # / / @@ -35,147 +35,147 @@ from TensorMesh import * # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k -def subarray(T,i1,i2,i3): - return take(take(take(T,i1,0),i2,1),i3,2) - +def subarray(T, i1, i2, i3): + return take(take(take(T, i1, 0), i2, 1), i3, 2) -def getEdgeInnerProduct(mesh,sigma): + +def getEdgeInnerProduct(mesh, sigma): h = mesh.h - m = array([size(h[0]),size(h[1]),size(h[2])]) + m = array([size(h[0]), size(h[1]), size(h[2])]) nc = prod(m) - + me1 = m + array([0, 1, 1]); ne1 = prod(me1) me2 = m + array([1, 0, 1]); ne2 = prod(me2) me3 = m + array([1, 1, 0]); ne3 = prod(me3) - + i = int64(linspace(0,m[0]-1,m[0])) j = int64(linspace(0,m[1]-1,m[1])) k = int64(linspace(0,m[2]-1,m[2])) - - ii,jj,kk = ndgrid(i,j,k,vector=False) + + ii,jj,kk = ndgrid(i,j,k,vector=False) ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk) - + ## -------- # no | node | e1 | e2 | e3 # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - ind1 = sub2ind(me1,c_[ii,jj,kk]) + ind1 = sub2ind(me1,c_[ii,jj,kk]) ind2 = sub2ind(me2,c_[ii,jj,kk]) + ne1 ind3 = sub2ind(me3,c_[ii,jj,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P000 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k - ind1 = sub2ind(me1,c_[ii,jj,kk]) + ind1 = sub2ind(me1,c_[ii,jj,kk]) ind2 = sub2ind(me2,c_[ii+1,jj,kk]) + ne1 ind3 = sub2ind(me3,c_[ii+1,jj,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P100 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k - ind1 = sub2ind(me1,c_[ii,jj+1,kk]) + ind1 = sub2ind(me1,c_[ii,jj+1,kk]) ind2 = sub2ind(me2,c_[ii,jj,kk]) + ne1 ind3 = sub2ind(me3,c_[ii,jj+1,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P010 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k - ind1 = sub2ind(me1,c_[ii,jj+1,kk]) + ind1 = sub2ind(me1,c_[ii,jj+1,kk]) ind2 = sub2ind(me2,c_[ii+1,jj,kk]) + ne1 ind3 = sub2ind(me3,c_[ii+1,jj+1,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P110 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ###### - + ## -------- # no | node | e1 | e2 | e3 # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k - ind1 = sub2ind(me1,c_[ii,jj,kk+1]) + ind1 = sub2ind(me1,c_[ii,jj,kk+1]) ind2 = sub2ind(me2,c_[ii,jj,kk+1]) + ne1 ind3 = sub2ind(me3,c_[ii,jj,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P001 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k+1 - ind1 = sub2ind(me1,c_[ii,jj,kk+1]) + ind1 = sub2ind(me1,c_[ii,jj,kk+1]) ind2 = sub2ind(me2,c_[ii+1,jj,kk+1]) + ne1 ind3 = sub2ind(me3,c_[ii+1,jj,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P101 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k+1 - ind1 = sub2ind(me1,c_[ii,jj+1,kk+1]) + ind1 = sub2ind(me1,c_[ii,jj+1,kk+1]) ind2 = sub2ind(me2,c_[ii,jj,kk+1]) + ne1 ind3 = sub2ind(me3,c_[ii,jj+1,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P011 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - + ## -------- # no | node | e1 | e2 | e3 # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k+1 - ind1 = sub2ind(me1,c_[ii,jj+1,kk+1]) + ind1 = sub2ind(me1,c_[ii,jj+1,kk+1]) ind2 = sub2ind(me2,c_[ii+1,jj,kk+1]) + ne1 ind3 = sub2ind(me3,c_[ii+1,jj+1,kk]) + ne1 + ne2 - + IND = vstack((vstack((ind1,ind2)),ind3)) IND = array(IND).flatten() - + P111 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr() - - - + + + # Cell volume - v = sqrt(mesh.vol) - 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]))) + 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)) - - v3 = r_[v,v,v] + + v = sqrt(mesh.vol) + v3 = r_[v, v, v] V = sdiag(v3)*Sigma*sdiag(v3) - - 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 - + return A - - + + if __name__ == '__main__': - h = [array([1,2,3,4]),array([1,2,1,4,2]),array([1,1,4,1])] + h = [array([1, 2, 3, 4]), array([1, 2, 1, 4, 2]), array([1, 1, 4, 1])] mesh = TensorMesh(h) - sigma = ones((mesh.nC,6)) - A = getEdgeInnerProduct(mesh,sigma) \ No newline at end of file + sigma = ones((mesh.nC, 6)) + A = getEdgeInnerProduct(mesh, sigma) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 08fa4caf..68ed8516 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -2,18 +2,21 @@ import numpy as np import unittest import sys sys.path.append('../') -from TensorMesh import TensorMesh from OrderTest import OrderTest -from scipy.sparse.linalg import dsolve -from getEdgeInnerProducts import getEdgeInnerProducts +from getEdgeInnerProducts import * -class TestNodalGrad(OrderTest): - name = "Nodal Gradient" - +class TestEdgeInnerProduct(OrderTest): + """Integrate a function over a unit cube domain.""" + + name = "Edge Inner Product" + meshSizes = [4, 8, 16, 32] def getError(self): + + call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + ex = lambda x, y, z: x**2+y*z ey = lambda x, y, z: (z**2)*x+y*z ez = lambda x, y, z: y**2+x*z @@ -24,25 +27,25 @@ class TestNodalGrad(OrderTest): sigma4 = lambda x, y, z: 0.1*x*y*z sigma5 = lambda x, y, z: 0.2*x*y sigma6 = lambda x, y, z: 0.1*z - - Ex = ex(self.M.gridEx[:, 0],self.M.gridEx[:, 1],self.M.gridEx[:, 2]) - Ey = ey(self.M.gridEy[:, 0],self.M.gridEy[:, 1],self.M.gridEy[:, 2]) - Ez = ez(self.M.gridEz[:, 0],self.M.gridEz[:, 1],self.M.gridEz[:, 2]) - - E = np.r_[Ex,Ey,Ez] + + Ex = call(ex, self.M.gridEx) + Ey = call(ey, self.M.gridEy) + Ez = call(ez, self.M.gridEz) + + E = np.matrix(mkvc(np.r_[Ex, Ey, Ez], 2)) Gc = self.M.gridCC - sigma = np.c_[sigma1(Gc[:,0],Gc[:,1],Gc[:,2]), - sigma2(Gc[:,0],Gc[:,1],Gc[:,2]), - sigma3(Gc[:,0],Gc[:,1],Gc[:,2]), - sigma4(Gc[:,0],Gc[:,1],Gc[:,2]), - sigma5(Gc[:,0],Gc[:,1],Gc[:,2]), - sigma6(Gc[:,0],Gc[:,1],Gc[:,2])] - - A = getEdgeInnerProducts(self.M, sigma) - - err = np.abs(E.T*A*E - 69881./21600) + 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) + numeric = E.T*A*E + analytic = 69881./21600 # Found using matlab symbolic toolbox. + err = np.abs(numeric - analytic) return err def test_order(self): self.orderTest() + + +if __name__ == '__main__': + unittest.main()