From 87331c4c9238db17f8f9ca44591485b8dcbc7f0f Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Fri, 26 Jul 2013 11:22:52 -0700 Subject: [PATCH 1/5] mass matrices for anisotropic sigma --- SimPEG/getEdgeInnerProducts.py | 181 ++++++++++++++++++++++++++++++ SimPEG/massMatrices.py | 94 ++++++++++++++++ SimPEG/subArray.py | 8 ++ SimPEG/tests/test_massMatrices.py | 48 ++++++++ SimPEG/utils.py | 31 ++++- 5 files changed, 361 insertions(+), 1 deletion(-) create mode 100644 SimPEG/getEdgeInnerProducts.py create mode 100644 SimPEG/massMatrices.py create mode 100644 SimPEG/subArray.py create mode 100644 SimPEG/tests/test_massMatrices.py diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py new file mode 100644 index 00000000..5a95f146 --- /dev/null +++ b/SimPEG/getEdgeInnerProducts.py @@ -0,0 +1,181 @@ +from scipy.sparse import linalg +from scipy import sparse +from sputils import * +from utils import * +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) +# / / +# / / | +# 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 + + +def subarray(T,i1,i2,i3): + return take(take(take(T,i1,0),i2,1),i3,2) + + +def getEdgeInnerProduct(mesh,sigma): + + h = mesh.h + 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 = 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]))) + Sigma = sp.vstack((row1, row2, row3)) + + 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 = 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])] + mesh = TensorMesh(h) + sigma = ones((mesh.nC,6)) + A = getEdgeInnerProduct(mesh,sigma) \ No newline at end of file diff --git a/SimPEG/massMatrices.py b/SimPEG/massMatrices.py new file mode 100644 index 00000000..ead21e3a --- /dev/null +++ b/SimPEG/massMatrices.py @@ -0,0 +1,94 @@ +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 diff --git a/SimPEG/subArray.py b/SimPEG/subArray.py new file mode 100644 index 00000000..2811f350 --- /dev/null +++ b/SimPEG/subArray.py @@ -0,0 +1,8 @@ +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/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py new file mode 100644 index 00000000..08fa4caf --- /dev/null +++ b/SimPEG/tests/test_massMatrices.py @@ -0,0 +1,48 @@ +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 + + +class TestNodalGrad(OrderTest): + name = "Nodal Gradient" + + meshSizes = [4, 8, 16, 32] + + def getError(self): + 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 + + sigma1 = lambda x, y, z: x*y+1 + sigma2 = lambda x, y, z: x*z+2 + sigma3 = lambda x, y, z: 3+z*y + 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] + 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) + + return err + + def test_order(self): + self.orderTest() diff --git a/SimPEG/utils.py b/SimPEG/utils.py index 4937062b..b9e17847 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') @@ -97,3 +97,32 @@ def ndgrid(*args, **kwargs): return np.c_[X1, X2, X3] else: return XYZ[2], XYZ[1], XYZ[0] + + +def ind2sub(shape, ind): + """From the given shape, returns the subscrips of the given index""" + revshp = [] + revshp.extend(shape) + mult = [1] + for i in range(0, len(revshp)-1): + mult.extend([mult[i]*revshp[i]]) + mult = 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]) + return sub + + +def sub2ind(shape, subs): + """From the given shape, returns the index of the given subscript""" + revshp = list(shape) + mult = [1] + for i in range(0, len(revshp)-1): + mult.extend([mult[i]*revshp[i]]) + mult = array(mult).reshape(len(mult), 1) + + idx = dot((subs), (mult)) + return idx \ No newline at end of file From 150cbc7df3b7c94e4173629b5130776079289aaa Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 12:11:45 -0700 Subject: [PATCH 2/5] 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() From 90bf6d31390d6d89e8ae795cd46b9eae1f622eeb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 13:56:05 -0700 Subject: [PATCH 3/5] Refactored and cleaned up code. --- SimPEG/getEdgeInnerProducts.py | 198 ++++++++---------------------- SimPEG/tests/test_massMatrices.py | 2 +- 2 files changed, 51 insertions(+), 149 deletions(-) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py index d3b91211..858d5a72 100644 --- a/SimPEG/getEdgeInnerProducts.py +++ b/SimPEG/getEdgeInnerProducts.py @@ -1,160 +1,61 @@ -from scipy.sparse import linalg -from scipy import sparse -from sputils import * -from utils import * -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) -# / / -# / / | -# 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 - - -def subarray(T, i1, i2, i3): - return take(take(take(T, i1, 0), i2, 1), i3, 2) +from scipy import sparse as sp +from sputils import sdiag +from utils import sub2ind, ndgrid +import numpy as np def getEdgeInnerProduct(mesh, sigma): - h = mesh.h - m = array([size(h[0]), size(h[1]), size(h[2])]) - nc = prod(m) + m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) + nc = mesh.nC - 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, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) - 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])) + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - ii,jj,kk = ndgrid(i,j,k,vector=False) - ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk) + 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), (np.linspace(0, 3*nc-1, 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 - 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]) - 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]) - 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]) - 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]) - 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]) - 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]) - 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]) - 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() - - + # 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]]) # Cell volume row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) @@ -162,8 +63,8 @@ def getEdgeInnerProduct(mesh, sigma): row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) Sigma = sp.vstack((row1, row2, row3)) - v = sqrt(mesh.vol) - v3 = r_[v, v, v] + v = np.sqrt(mesh.vol) + v3 = np.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 @@ -175,7 +76,8 @@ def getEdgeInnerProduct(mesh, sigma): if __name__ == '__main__': - h = [array([1, 2, 3, 4]), array([1, 2, 1, 4, 2]), array([1, 1, 4, 1])] + from TensorMesh import * + h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] mesh = TensorMesh(h) - sigma = ones((mesh.nC, 6)) + sigma = np.ones((mesh.nC, 6)) A = getEdgeInnerProduct(mesh, sigma) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 68ed8516..87310a47 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -32,7 +32,7 @@ class TestEdgeInnerProduct(OrderTest): Ey = call(ey, self.M.gridEy) Ez = call(ez, self.M.gridEz) - E = np.matrix(mkvc(np.r_[Ex, Ey, Ez], 2)) + 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)] From c8633881fbdfeaaaa0cec92a60b24c155d1cf4a8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 14:19:44 -0700 Subject: [PATCH 4/5] Support for isotropic sigma --- SimPEG/getEdgeInnerProducts.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py index 858d5a72..7e123ab6 100644 --- a/SimPEG/getEdgeInnerProducts.py +++ b/SimPEG/getEdgeInnerProducts.py @@ -1,6 +1,6 @@ from scipy import sparse as sp from sputils import sdiag -from utils import sub2ind, ndgrid +from utils import sub2ind, ndgrid, mkvc import numpy as np @@ -21,7 +21,7 @@ def getEdgeInnerProduct(mesh, sigma): IND = np.r_[ind1, ind2, ind3].flatten() - return sp.coo_matrix((np.ones(3*nc), (np.linspace(0, 3*nc-1, 3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr() + 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) # / / @@ -57,15 +57,21 @@ def getEdgeInnerProduct(mesh, sigma): P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) - # Cell volume - 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)) + if sigma.size == mesh.nC: # Isotropic! + sigma = mkvc(sigma) + 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) + 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 @@ -73,10 +79,8 @@ def getEdgeInnerProduct(mesh, sigma): return A - if __name__ == '__main__': - - from TensorMesh import * + 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)) From 5c83095781f67e1c315bf2a7889ec09004964526 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 14:30:04 -0700 Subject: [PATCH 5/5] Minor changes and comments. --- SimPEG/getEdgeInnerProducts.py | 2 +- SimPEG/tests/test_massMatrices.py | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py index 7e123ab6..68423a3f 100644 --- a/SimPEG/getEdgeInnerProducts.py +++ b/SimPEG/getEdgeInnerProducts.py @@ -58,7 +58,7 @@ def getEdgeInnerProduct(mesh, sigma): P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) if sigma.size == mesh.nC: # Isotropic! - sigma = mkvc(sigma) + 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]]) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 87310a47..0aae0615 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,18 +1,16 @@ import numpy as np import unittest +from OrderTest import OrderTest import sys sys.path.append('../') -from OrderTest import OrderTest from getEdgeInnerProducts import * class TestEdgeInnerProduct(OrderTest): - """Integrate a function over a unit cube domain.""" + """Integrate an edge function over a unit cube domain using edgeInnerProducts.""" name = "Edge Inner Product" - meshSizes = [4, 8, 16, 32] - def getError(self): call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])