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