diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index f9d1f0aa..a755662f 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -1,6 +1,16 @@ import numpy as np -from scipy import sparse -from sputils import ddx, sdiag, speye, kron3, spzeros, av +from scipy import sparse as sp +from sputils import sdiag, speye, kron3, spzeros + + +def ddx(n): + """Define 1D derivatives""" + return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") + + +def av(n): + """Define 1D averaging operator""" + return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") class DiffOperators(object): @@ -16,20 +26,20 @@ class DiffOperators(object): def fget(self): if(self._faceDiv is None): # The number of cell centers in each direction - n = [hk.size for hk in self.h] + n = self.n # Compute faceDivergence operator on faces dd = [ddx(k) for k in n] if(self.dim == 1): D = dd[0] elif(self.dim == 2): - D1 = sparse.kron(speye(n[1]), dd[0]) - D2 = sparse.kron(dd[1], speye(n[0])) - D = sparse.hstack((D1, D2), format="csr") + D1 = sp.kron(speye(n[1]), dd[0]) + D2 = sp.kron(dd[1], speye(n[0])) + D = sp.hstack((D1, D2), format="csr") elif(self.dim == 3): D1 = kron3(speye(n[2]), speye(n[1]), dd[0]) D2 = kron3(speye(n[2]), dd[1], speye(n[0])) D3 = kron3(dd[2], speye(n[1]), speye(n[0])) - D = sparse.hstack((D1, D2, D3), format="csr") + D = sp.hstack((D1, D2, D3), format="csr") # Compute areas of cell faces S = self.area # Compute cell volumes @@ -62,7 +72,7 @@ class DiffOperators(object): D2 = kron3(speye(n3+1), d2, speye(n1+1)) D3 = kron3(d3, speye(n2+1), speye(n1+1)) - G = sparse.vstack((D1, D2, D3), format="csr") + G = sp.vstack((D1, D2, D3), format="csr") self._nodalGrad = sdiag(1/L)*G return self._nodalGrad return locals() @@ -101,9 +111,9 @@ class DiffOperators(object): O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - C = sparse.vstack((sparse.hstack((O1, -D32, D23)), - sparse.hstack((D31, O2, -D13)), - sparse.hstack((-D21, D12, O3))), format="csr") + C = sp.vstack((sp.hstack((O1, -D32, D23)), + sp.hstack((D31, O2, -D13)), + sp.hstack((-D21, D12, O3))), format="csr") self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) return self._edgeCurl @@ -116,18 +126,16 @@ class DiffOperators(object): def fget(self): if(self._faceAve is None): - # The number of cell centers in each direction - n1 = np.size(self.hx) - n2 = np.size(self.hy) - n3 = np.size(self.hz) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - self._faceAve = sparse.hstack(kron3(speye(n3), speye(n2), av1), - kron3(speye(n3), av2, speye(n3)), - kron3(av3, speye(n2), speye(n3)), format="csr") + n = self.n + if(self.dim == 1): + self._faceAve = av(n[0]) + elif(self.dim == 2): + self._faceAve = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") + elif(self.dim == 3): + self._faceAve = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._faceAve return locals() _faceAve = None @@ -139,17 +147,16 @@ class DiffOperators(object): def fget(self): if(self._edgeAve is None): # The number of cell centers in each direction - n1 = np.size(self.hx) - n2 = np.size(self.hy) - n3 = np.size(self.hz) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - self._edgeAve = sparse.hstack(kron3(av3, av2, speye(n1)), - kron3(av3, speye(n2), av1), - kron3(speye(n3), av2, av1), format="csr") + n = self.n + if(self.dim == 1): + raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + elif(self.dim == 2): + self._edgeAve = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + sp.kron(speye(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._edgeAve = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._edgeAve return locals() _edgeAve = None diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index 023362ae..8395c54d 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -22,15 +22,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): """ def __init__(self, h, x0=None): - super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0) + super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" for i, h_i in enumerate(h): assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) + assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) # Ensure h contains 1D vectors - self._h = [x.ravel() for x in h] + self._h = [mkvc(x) for x in h] def h(): doc = "h is a list containing the cell widths of the tensor mesh in each dimension." @@ -193,17 +194,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): def fget(self): if(self._vol is None): - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # Compute cell volumes if(self.dim == 1): self._vol = mkvc(vh[0]) elif(self.dim == 2): # Cell sizes in each direction - self._vol = mkvc(vh[0]*vh[1].T) + self._vol = mkvc(np.outer(vh[0], vh[1])) elif(self.dim == 3): # Cell sizes in each direction - v12 = vh[0]*vh[1].T - self._vol = mkvc(mkvc(v12, 2)*vh[2].T) + self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2])) return self._vol return locals() _vol = None @@ -215,20 +215,20 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): def fget(self): if(self._area is None): # Ensure that we are working with column vectors - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # The number of cell centers in each direction - n = [x.size for x in self.h] + n = self.n # Compute areas of cell faces if(self.dim == 1): - self._area = np.ones((n[0]+1, 1)) + self._area = np.ones(n[0]+1) elif(self.dim == 2): - area1 = np.ones((n[0]+1, 1))*vh[1].T - area2 = vh[0]*np.ones((n[1]+1, 1)).T + area1 = np.outer(np.ones(n[0]+1), vh[1]) + area2 = np.outer(vh[0], np.ones(n[1]+1)) self._area = np.r_[mkvc(area1), mkvc(area2)] elif(self.dim == 3): - area1 = np.ones((n[0]+1, 1))*mkvc(vh[1]*vh[2].T) - area2 = vh[0]*mkvc(np.ones((n[1]+1, 1))*vh[2].T) - area3 = vh[0]*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1)))) self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] return self._area return locals() @@ -239,22 +239,22 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): doc = "Construct edge legnths of the 3D model as 1d array." def fget(self): - if(self._area is None): + if(self._edge is None): # Ensure that we are working with column vectors - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # The number of cell centers in each direction - n = [x.size for x in self.h] + n = self.n # Compute edge lengths if(self.dim == 1): self._edge = mkvc(vh[0]) elif(self.dim == 2): - l1 = vh[0]*np.ones((n[1]+1, 1)).T - l2 = np.ones((n[0]+1, 1))*vh[1].T + l1 = np.outer(vh[0], np.ones(n[1]+1)) + l2 = np.outer(np.ones(n[0]+1), vh[1]) self._edge = np.r_[mkvc(l1), mkvc(l2)] elif(self.dim == 3): - l1 = vh[0]*mkvc(np.ones((n[1]+1, 1))*np.ones((n[2]+1, 1)).T) - l2 = np.ones((n[0]+1, 1))*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) - l3 = np.ones((n[0]+1, 1))*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2]))) self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] return self._edge return locals() diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index f1e9a677..79bbdf04 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -1,43 +1,30 @@ import numpy as np -from scipy import sparse -from numpy import ones +from scipy import sparse as sp -def ddx(n): - """Define 1D derivatives""" - return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1, format="csr") - def sdiag(h): """Sparse diagonal matrix""" - return sparse.spdiags(h, 0, np.size(h), np.size(h), format="csr") + return sp.spdiags(h, 0, np.size(h), np.size(h), format="csr") + def speye(n): """Sparse identity""" - return sparse.identity(n, format="csr") + return sp.identity(n, format="csr") + def kron3(A, B, C): """Two kron prods""" - return sparse.kron(sparse.kron(A, B), C, format="csr") + return sp.kron(sp.kron(A, B), C, format="csr") + def spzeros(n1, n2): """spzeros""" - return sparse.coo_matrix((n1, n2)).tocsr() + return sp.coo_matrix((n1, n2)).tocsr() -def av(n): - """Define 1D averaging operator""" - return sparse.spdiags((0.5*np.ones((n+1,1))*[1,1]).T, [0,1], n, n+1, format="csr") - -def av(n): - """Define 1D average""" - return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)) - -def spzeros(n1, n2): - """spzeros""" - return sparse.coo_matrix((n1, n2)) def appendBottom(A, B): """append on bottom""" - C = sparse.vstack((A, B)) + C = sp.vstack((A, B)) C = C.tocsr() return C @@ -51,7 +38,7 @@ def appendBottom3(A, B, C): def appendRight(A, B): """append on right""" - C = sparse.hstack((A, B)) + C = sp.hstack((A, B)) C = C.tocsr() return C @@ -65,9 +52,9 @@ def appendRight3(A, B, C): def blkDiag(A, B): """blockdigonal""" - O12 = sparse.coo_matrix((np.shape(A)[0], np.shape(B)[1])) - O21 = sparse.coo_matrix((np.shape(B)[0], np.shape(A)[1])) - C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B)))) + O12 = sp.coo_matrix((np.shape(A)[0], np.shape(B)[1])) + O21 = sp.coo_matrix((np.shape(B)[0], np.shape(A)[1])) + C = sp.vstack((sp.hstack((A, O12)), sp.hstack((O21, B)))) C = C.tocsr() return C diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py index d8ee0018..e90f8cfa 100644 --- a/SimPEG/tests/test_curl.py +++ b/SimPEG/tests/test_curl.py @@ -11,12 +11,11 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] - x0 = np.zeros((3, 1)) - M = TensorMesh(h, x0) + M = TensorMesh(h) #n = M.plotGrid() # Generate DIV matrix @@ -41,7 +40,7 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index dd469a3f..53ca506b 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -11,9 +11,9 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] x0 = np.zeros((3, 1)) M = TensorMesh(h, x0) @@ -24,9 +24,9 @@ for i in range(4): #Test function fun = lambda x: np.sin(x) - Fx = fun(M.gridFx[:,0]) - Fy = fun(M.gridFy[:,1]) - Fz = fun(M.gridFz[:,2]) + Fx = fun(M.gridFx[:, 0]) + Fy = fun(M.gridFy[:, 1]) + Fz = fun(M.gridFz[:, 2]) F = np.concatenate((Fx,Fy,Fz)) divF = DIV*F @@ -39,7 +39,7 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py index 58401dc7..fbc36aeb 100644 --- a/SimPEG/tests/test_grad.py +++ b/SimPEG/tests/test_grad.py @@ -11,12 +11,11 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] - x0 = np.zeros((3, 1)) - M = TensorMesh(h, x0) + M = TensorMesh(h) #n = M.plotGrid() # Generate DIV matrix @@ -38,8 +37,8 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err