From 7e5ea835e321c664db45ca22aa436c41a0e54f59 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 6 Aug 2013 17:56:21 -0700 Subject: [PATCH] Removed eldadCode directory as it is now integrated. Added GaussNewton to the main directory. --- SimPEG/EldadsCode/MFDdriver.py | 36 ---- SimPEG/EldadsCode/getCellVolume.py | 58 ------ SimPEG/EldadsCode/getDiffOps.py | 106 ----------- SimPEG/EldadsCode/getEdgeInnerProduct.py | 213 ---------------------- SimPEG/EldadsCode/getEdgeTangent.py | 60 ------ SimPEG/EldadsCode/getFaceInnerProduct.py | 86 --------- SimPEG/EldadsCode/getFaceNormals.py | 73 -------- SimPEG/EldadsCode/getVolume.py | 34 ---- SimPEG/EldadsCode/inv3X3BlockDiagonal.py | 36 ---- SimPEG/EldadsCode/meshUtils.py | 94 ---------- SimPEG/EldadsCode/sputils.py | 77 -------- SimPEG/EldadsCode/tools.py | 222 ----------------------- SimPEG/EldadsCode/utils.py | 87 --------- SimPEG/EldadsCode/zevel.py | 119 ------------ SimPEG/{EldadsCode => }/GaussNewton.py | 0 15 files changed, 1301 deletions(-) delete mode 100644 SimPEG/EldadsCode/MFDdriver.py delete mode 100644 SimPEG/EldadsCode/getCellVolume.py delete mode 100644 SimPEG/EldadsCode/getDiffOps.py delete mode 100644 SimPEG/EldadsCode/getEdgeInnerProduct.py delete mode 100644 SimPEG/EldadsCode/getEdgeTangent.py delete mode 100644 SimPEG/EldadsCode/getFaceInnerProduct.py delete mode 100644 SimPEG/EldadsCode/getFaceNormals.py delete mode 100644 SimPEG/EldadsCode/getVolume.py delete mode 100644 SimPEG/EldadsCode/inv3X3BlockDiagonal.py delete mode 100644 SimPEG/EldadsCode/meshUtils.py delete mode 100644 SimPEG/EldadsCode/sputils.py delete mode 100644 SimPEG/EldadsCode/tools.py delete mode 100644 SimPEG/EldadsCode/utils.py delete mode 100644 SimPEG/EldadsCode/zevel.py rename SimPEG/{EldadsCode => }/GaussNewton.py (100%) diff --git a/SimPEG/EldadsCode/MFDdriver.py b/SimPEG/EldadsCode/MFDdriver.py deleted file mode 100644 index 16d18fe3..00000000 --- a/SimPEG/EldadsCode/MFDdriver.py +++ /dev/null @@ -1,36 +0,0 @@ -import numpy as np -from numpy.random import randn -from utils import ndgrid -from getDiffOps import getCurlMatrix, getNodalGradient -from getFaceInnerProduct import getFaceInnerProduct -from getEdgeInnerProduct import getEdgeInnerProduct -from scipy.sparse.linalg import dsolve -from pylab import norm - -n = np.array([14, 14, 15]) - -X, Y, Z = ndgrid(*[np.linspace(0, 1, x) for x in n]) -sigma = 1e-2*np.ones(n-1) -sigma[:, :, (n[2]-1)/2:] = 1e-6 -mu = 4*np.pi*1e-7*np.ones(n-1) -w = 10 - -CURL = getCurlMatrix(X, Y, Z) -GRAD = getNodalGradient(X, Y, Z) -Mf = getFaceInnerProduct(X, Y, Z, 1/mu) -Me = getEdgeInnerProduct(X, Y, Z, sigma) - -A = CURL.T * Mf * CURL + 1j * w * Me - -ne = np.shape(A) -b = np.matrix(randn(ne[0])).T -# clean b -DIVb = GRAD.T*b -p = dsolve.spsolve(GRAD.T*GRAD, DIVb, use_umfpack=True).T -b = b - GRAD*p - -#x = spsolve(A, b) -x = dsolve.spsolve(A, b, use_umfpack=True).T - -t = norm(A*x-b)/norm(b) -print t diff --git a/SimPEG/EldadsCode/getCellVolume.py b/SimPEG/EldadsCode/getCellVolume.py deleted file mode 100644 index bb536842..00000000 --- a/SimPEG/EldadsCode/getCellVolume.py +++ /dev/null @@ -1,58 +0,0 @@ -from sputils import * -from utils import * -from sputils import * -from numpy import * -from getEdgeTangent import * - - -def volTetra(y, m, I, A, B, C, D): - - a11 = array(y[A, 0]-y[B, 0]); a12 = array(y[A, 0]-y[C, 0]); a13 = array(y[A, 0]-y[D, 0]) - a21 = array(y[A, 1]-y[B, 1]); a22 = array(y[A, 1]-y[C, 1]); a23 = array(y[A, 1]-y[D, 1]) - a31 = array(y[A, 2]-y[B, 2]); a32 = array(y[A, 2]-y[C, 2]); a33 = array(y[A, 2]-y[D, 2]) - - return abs(a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a31*a22*a13 - a32*a23*a11 - a33*a21*a12) - - -def getCellVolume(X, Y, Z): - - m = array(shape(X))-1 - y = hstack3(mkvc(X), mkvc(Y), mkvc(Z)) - - 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) - ii = mkvc(ii) - jj = mkvc(jj) - kk = mkvc(kk) - - I = int64(sub2ind(m, hstack3(ii, jj, kk))) - A = int64(sub2ind(m+1, hstack3(ii, jj, kk))) - B = int64(sub2ind(m+1, hstack3(ii, jj+1, kk))) - C = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk))) - D = int64(sub2ind(m+1, hstack3(ii+1, jj, kk))) - E = int64(sub2ind(m+1, hstack3(ii, jj, kk+1))) - F = int64(sub2ind(m+1, hstack3(ii, jj+1, kk+1))) - G = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk+1))) - H = int64(sub2ind(m+1, hstack3(ii+1, jj, kk+1))) - - v1 = volTetra(y, m, I, A, B, D, E) - v2 = volTetra(y, m, I, B, E, F, G) - v3 = volTetra(y, m, I, B, D, E, G) - v4 = volTetra(y, m, I, B, C, D, G) - v5 = volTetra(y, m, I, D, E, G, H) - - v = 1.0/6.0 * (v1 + v2 + v3 + v4 + v5) - return v.flatten() - - -if __name__ == '__main__': - - X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3)) - Z[2, 2, 2] = 2.5; Z[0, 0, 0] = -0.5 - X[2, 2, 2] = 2.5; X[0, 0, 0] = -0.5 - - v = getCellVolume(X, Y, Z) - print v diff --git a/SimPEG/EldadsCode/getDiffOps.py b/SimPEG/EldadsCode/getDiffOps.py deleted file mode 100644 index dbb4e5b0..00000000 --- a/SimPEG/EldadsCode/getDiffOps.py +++ /dev/null @@ -1,106 +0,0 @@ -from sputils import * -from utils import * -from numpy import * -from getEdgeTangent import * -from getCellVolume import getCellVolume -from getFaceNormals import getFaceNormals - - -def getDivMatrix(X, Y, Z): - """Face DIV""" - - n = array(shape(X))-1 - n1 = n[0] - n2 = n[1] - n3 = n[2] - - n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z) - - area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - S = sdiag(area) - V = getCellVolume(X, Y, Z) - - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - # divergence on faces - D = appendRight3(D1, D2, D3) - - return sdiag(1/V)*D*S - - -def getCurlMatrix(X, Y, Z): - """Edge CURL """ - - n = array(shape(X))-1 - n1 = n[0]; n2 = n[1]; n3 = n[2] - - d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3) - # derivatives on x-edge variables - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(shape(D32)[0], shape(D31)[1]) - O2 = spzeros(shape(D31)[0], shape(D32)[1]) - O3 = spzeros(shape(D21)[0], shape(D13)[1]) - - CURL = appendBottom3( - appendRight3(O1, -D32, D23), - appendRight3(D31, O2, -D13), - appendRight3(-D21, D12, O3)) - - # scale for non-uniform mesh - e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z) - n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z) - - area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - S = sdiag(1/area) - lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3))) - L = sdiag(lngth) - - return S*(CURL*L) - - -def getNodalGradient(X, Y, Z): - """Nodal Gradients""" - - n = array(shape(X))-1 - n1 = n[0]; n2 = n[1]; n3 = n[2] - - D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1)) - D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1)) - D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1)) - - # topological gradient - GRAD = appendBottom3(D1, D2, D3) - - # scale for non-uniform mesh - e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z) - lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3))) - L = sdiag(1/lngth) - - return L*GRAD - - -if __name__ == '__main__': - - X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3)) - Z[2, 2, 2] = 2.5 - Z[0, 0, 0] = -0.5 - X[2, 2, 2] = 2.5 - X[0, 0, 0] = -0.5 - sig = ones([2, 2, 2]) - C = getCurlMatrix(X, Y, Z) - - G = getNodalGradient(X, Y, Z) - - tt = C*G - print(tt) diff --git a/SimPEG/EldadsCode/getEdgeInnerProduct.py b/SimPEG/EldadsCode/getEdgeInnerProduct.py deleted file mode 100644 index 776f31f7..00000000 --- a/SimPEG/EldadsCode/getEdgeInnerProduct.py +++ /dev/null @@ -1,213 +0,0 @@ -from scipy.sparse import linalg -from scipy import sparse -from sputils import * -from utils import * -from sputils import * -from numpy import * -from getEdgeTangent import * -from inv3X3BlockDiagonal import * -from getCellVolume import getCellVolume - - -def subarray(T, i1, i2, i3): - return take(take(take(T, i1, 0), i2, 1), i3, 2) - - -def getEdgeInnerProduct(X, Y, Z, sigma): - """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 - """ - m = array(shape(X))-1 - 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) - - e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,norme1,norme2,norme3 = getEdgeTangent(X,Y,Z) - - 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) - 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,hstack3(ii,jj,kk)) - ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT000 = inv3X3BlockDiagonal(subarray(e1x,i,j,k) , subarray(e1y,i,j,k), subarray(e1z,i,j,k), - subarray(e2x,i,j,k) , subarray(e2y,i,j,k), subarray(e2z,i,j,k) , - subarray(e3x,i,j,k) , subarray(e3y,i,j,k), subarray(e3z,i,j,k) ) - - ## -------- - # no | node | e1 | e2 | e3 - # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k - ind1 = sub2ind(me1,hstack3(ii,jj,kk)) - ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT100 = inv3X3BlockDiagonal(subarray(e1x,i,j,k), subarray(e1y,i,j,k), subarray(e1z,i,j,k), - subarray(e2x,i+1,j,k), subarray(e2y,i+1,j,k), subarray(e2z,i+1,j,k), - subarray(e3x,i+1,j,k) , subarray(e3y,i+1,j,k), subarray(e3z,i+1,j,k)) - - ## -------- - # no | node | e1 | e2 | e3 - # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k - ind1 = sub2ind(me1,hstack3(ii,jj+1,kk)) - ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT010 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) , subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) , - subarray(e2x,i,j,k) , subarray(e2y,i,j,k) , subarray(e2z,i,j,k) , - subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) ,subarray( e3z,i,j+1,k) ) - - ## -------- - # 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,hstack3(ii,jj+1,kk)) - ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT110 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) ,subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) , - subarray(e2x,i+1,j,k) ,subarray(e2y,i+1,j,k) , subarray(e2z,i+1,j,k), - subarray(e3x,i+1,j+1,k) ,subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) ) - - ###### - - ## -------- - # no | node | e1 | e2 | e3 - # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k - ind1 = sub2ind(me1,hstack3(ii,jj,kk+1)) - ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT001 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1) ,subarray(e1y,i,j,k+1) , subarray(e1z,i,j,k+1) , - subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) , - subarray(e3x,i,j,k) , subarray(e3y,i,j,k) , subarray(e3z,i,j,k) ) - - ## -------- - # 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,hstack3(ii,jj,kk+1)) - ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT101 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1), subarray(e1y,i,j,k+1), subarray(e1z,i,j,k+1) , - subarray(e2x,i+1,j,k+1), subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) , - subarray(e3x,i+1,j,k), subarray(e3y,i+1,j,k) , subarray(e3z,i+1,j,k) ) - - ## -------- - # 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,hstack3(ii,jj+1,kk+1)) - ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT011 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) , - subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) , - subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) , subarray(e3z,i,j+1,k) ) - - ## -------- - # 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,hstack3(ii,jj+1,kk+1)) - ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1 - ind3 = sub2ind(me3,hstack3(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() - - invT111 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) , - subarray(e2x,i+1,j,k+1) , subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) , - subarray(e3x,i+1,j+1,k) , subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) ) - - # Cell volume - v = mkvc(getCellVolume(X,Y,Z)) #mkvc(getVolume(X,Y,Z)) - vsig = v*mkvc(sigma) - v3 = vstack((vstack((vsig,vsig)),vsig)) - v3 = v3.flatten() - - V = sdiag(v3) - - A = P000.T*invT000.T*V*invT000*P000 + P001.T*invT001.T*V*invT001*P001 + P010.T*invT010.T*V*invT010*P010 + P011.T*invT011.T*V*invT011*P011 + P100.T*invT100.T*V*invT100*P100 + P101.T*invT101.T*V*invT101*P101 + P110.T*invT110.T*V*invT110*P110 + P111.T*invT111.T*V*invT111*P111 - - A = 0.125*A - - return A - - -if __name__ == '__main__': - - X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3)) - Z[2,2,2] = 2.5; Z[0,0,0] = -0.5 - X[2,2,2] = 2.5; X[0,0,0] = -0.5 - sig = ones([2,2,2]) - A = getEdgeInnerProduct(X,Y,Z,sig) \ No newline at end of file diff --git a/SimPEG/EldadsCode/getEdgeTangent.py b/SimPEG/EldadsCode/getEdgeTangent.py deleted file mode 100644 index a8426df3..00000000 --- a/SimPEG/EldadsCode/getEdgeTangent.py +++ /dev/null @@ -1,60 +0,0 @@ -from numpy import * -from utils import diff - -#function[t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,normt1,normt2,normt3] = getEdgeTangent(X,Y,Z) -#%[t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,normt1,normt2,normt3] = getEdgeTangent(X,Y,Z) -#% -#% node(i,j,k+1) ------ edgt2(i,j,k+1) ----- node(i,j+1,k+1) -#% / / -#% / / | -#% edgt3(i,j,k) fact1(i,j,k) edgt3(i,j+1,k) -#% / / | -#% / / | -#% node(i,j,k) ------ edgt2(i,j,k) ----- node(i,j+1,k) -#% | | | -#% | | node(i+1,j+1,k+1) -#% | | / -#% edgt1(i,j,k) fact3(i,j,k) edgt1(i,j+1.k) -#% | | / -#% | | / -#% | |/ -#% node(i+1,j,k) ------ edgt2(i+1,j,k) ----- node(i+1,j+1,k) - - -def getEdgeTangent(X, Y, Z): - - t1x = diff(X, 1) - t1y = diff(Y, 1) - t1z = diff(Z, 1) - - normt1 = sqrt(t1x**2+t1y**2+t1z**2) - t1x = t1x/normt1 - t1y = t1y/normt1 - t1z = t1z/normt1 - - t2x = diff(X, 2) - t2y = diff(Y, 2) - t2z = diff(Z, 2) - normt2 = sqrt(t2x**2 + t2y**2 + t2z**2) - t2x = t2x/normt2 - t2y = t2y/normt2 - t2z = t2z/normt2 - - t3x = diff(X, 3) - t3y = diff(Y, 3) - t3z = diff(Z, 3) - normt3 = sqrt(t3x**2+t3y**2+t3z**2) - t3x = t3x/normt3 - t3y = t3y/normt3 - t3z = t3z/normt3 - - # print t3x - - return (t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z, normt1, normt2, normt3) - - -if __name__ == '__main__': - - X, Y, Z = mgrid[0:4, 0:5, 0:6] - - t = getEdgeTangent(X, Y, Z) diff --git a/SimPEG/EldadsCode/getFaceInnerProduct.py b/SimPEG/EldadsCode/getFaceInnerProduct.py deleted file mode 100644 index 2870eb88..00000000 --- a/SimPEG/EldadsCode/getFaceInnerProduct.py +++ /dev/null @@ -1,86 +0,0 @@ -from scipy.sparse import linalg -from scipy import sparse -from sputils import * -from utils import * -from numpy import * -from getEdgeTangent import * -from inv3X3BlockDiagonal import * -from getCellVolume import getCellVolume -from getFaceNormals import getFaceNormals - - -#----------------------- -def subarray(T,i1,i2,i3): - return take(take(take(T,i1,0),i2,1),i3,2) - -#----------------------- - -def getFaceInnerProduct(X,Y,Z,sigma): - - m = array(shape(X))-1 - nc = prod(m) - mf1 = m+[1, 0, 0] - mf2 = m+[0, 1, 0] - mf3 = m+[0, 0, 1] - - nf1 = prod(m+[1, 0, 0]) - nf2 = prod(m+[0, 1, 0]) - nf3 = prod(m+[0, 0, 1]) - - # compute the normals - n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3 = getFaceNormals(X,Y,Z) - - 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) - ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk) - - ind1 = sub2ind(mf1,hstack3(ii,jj,kk)) - ind2 = sub2ind(mf2,hstack3(ii,jj,kk)) + nf1 - ind3 = sub2ind(mf3,hstack3(ii,jj,kk)) + nf1 + nf2 - - IND = vstack((vstack((ind1,ind2)),ind3)) - IND = array(IND).flatten() - - P1 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,nf1+nf2+nf3)).tocsr() - - ind1 = sub2ind(mf1,hstack3(ii+1,jj,kk)) - ind2 = sub2ind(mf2,hstack3(ii,jj+1,kk)) + nf1 - ind3 = sub2ind(mf3,hstack3(ii,jj,kk+1)) + nf1 + nf2 - - IND = vstack((vstack((ind1,ind2)),ind3)) - IND = array(IND).flatten() - - P2 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,nf1+nf2+nf3)).tocsr() - - - invN1 = inv3X3BlockDiagonal(subarray(n1x,i,j,k) , subarray(n1y,i,j,k), subarray(n1z,i,j,k), - subarray(n2x,i,j,k) , subarray(n2y,i,j,k), subarray(n2z,i,j,k), - subarray(n3x,i,j,k) , subarray(n3y,i,j,k), subarray(n3z,i,j,k) ) - - - invN2 = inv3X3BlockDiagonal(subarray(n1x,i+1,j,k) , subarray(n1y,i+1,j,k), subarray(n1z,i+1,j,k), - subarray(n2x,i,j+1,k) , subarray(n2y,i,j+1,k), subarray(n2z,i,j+1,k), - subarray(n3x,i,j,k+1) , subarray(n3y,i,j,k+1), subarray(n3z,i,j,k+1) ) - - # Cell volume - v = mkvc(getCellVolume(X,Y,Z)) #mkvc(getVolume(X,Y,Z)) - vsig = v*mkvc(sigma) - v3 = vstack((vstack((vsig,vsig)),vsig)) - v3 = v3.flatten() - - V = sdiag(v3) - - return (P1.T*invN1.T*V*invN1*P1 + P2.T*invN2.T*V*invN2*P2)/2.0 - - -if __name__ == '__main__': - - X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3)) - Z[2,2,2] = 2.5; Z[0,0,0] = -0.5 - X[2,2,2] = 2.5; X[0,0,0] = -0.5 - sigma = ones([2,2,2]) - A = getFaceInnerProduct(X,Y,Z,sigma) - print(A) \ No newline at end of file diff --git a/SimPEG/EldadsCode/getFaceNormals.py b/SimPEG/EldadsCode/getFaceNormals.py deleted file mode 100644 index 2038a301..00000000 --- a/SimPEG/EldadsCode/getFaceNormals.py +++ /dev/null @@ -1,73 +0,0 @@ -from numpy import * -from utils import * - - -def getFaceNormals(X, Y, Z): -# compute the x normals - d1xp = diffp(X,2,3) - d1yp = diffp(Y,2,3) - d1zp = diffp(Z,2,3) - - d1xm = diffm(X,3,2) - d1ym = diffm(Y,3,2) - d1zm = diffm(Z,3,2) - - # normals - n1x = d1yp*d1zm - d1zp*d1ym - n1y = d1zp*d1xm - d1xp*d1zm - n1z = d1xp*d1ym - d1yp*d1xm - normn1 = sqrt(n1x**2 + n1y**2 + n1z**2) - n1x = n1x / normn1; - n1y = n1y / normn1; - n1z = n1z / normn1; - - area1 = normn1/2 - - - # compute the y normals - d2xp = diffp(X,1,3) - d2yp = diffp(Y,1,3) - d2zp = diffp(Z,1,3) - - d2xm = diffm(X,1,3) - d2ym = diffm(Y,1,3) - d2zm = diffm(Z,1,3) - - # normals - n2x = d2yp*d2zm - d2zp*d2ym - n2y = d2zp*d2xm - d2xp*d2zm - n2z = d2xp*d2ym - d2yp*d2xm - normn2 = sqrt(n2x**2 + n2y**2 + n2z**2) - n2x = n2x / normn2 - n2y = n2y / normn2 - n2z = n2z / normn2 - - area2 = normn2/2 - - # compute the z normals - d3xp = diffp(X,1,2) - d3yp = diffp(Y,1,2) - d3zp = diffp(Z,1,2) - - d3xm = diffm(X,2,1) - d3ym = diffm(Y,2,1) - d3zm = diffm(Z,2,1) - - # normals - n3x = d3yp*d3zm - d3zp*d3ym - n3y = d3zp*d3xm - d3xp*d3zm - n3z = d3xp*d3ym - d3yp*d3xm; - normn3 = sqrt(n3x**2 + n3y**2 + n3z**2); - n3x = n3x / normn3; - n3y = n3y / normn3; - n3z = n3z / normn3; - - area3 = normn3/2; - - return (n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3) - -if __name__ == '__main__': - - X, Y, Z = mgrid[0:4, 0:5, 0:6] - - t = getFaceNormals(X, Y, Z) diff --git a/SimPEG/EldadsCode/getVolume.py b/SimPEG/EldadsCode/getVolume.py deleted file mode 100644 index 3b4cd250..00000000 --- a/SimPEG/EldadsCode/getVolume.py +++ /dev/null @@ -1,34 +0,0 @@ -from numpy import * -from utils import diff, ave - - -def getVolume(X,Y,Z): - - # compute edge vectors - t1x = ave(ave(diff(X, 1),2),3) - t1y = ave(ave(diff(Y, 1),2),3) - t1z = ave(ave(diff(Z, 1),2),3) - - t2x = ave(ave(diff(X, 2),1),3) - t2y = ave(ave(diff(Y, 2),1),3) - t2z = ave(ave(diff(Z, 2),1),3) - - t3x = ave(ave(diff(X, 3),1),2) - t3y = ave(ave(diff(Y, 3),1),2) - t3z = ave(ave(diff(Z, 3),1),2) - - # v = [t1x t1y t1z][i j k] - # [t2x t2y t2z] - # [t3x t3y t3z] - - v = t1x*(t2y*t3z - t2z*t3y) - t1y*(t2x*t3z - t2z*t3x) + t1z*(t2x*t3y-t2y*t3x) - - return v - - -if __name__ == '__main__': - - X, Y, Z = mgrid[0:4, 0:5, 0:6] - X = (1.0*X)/2 - v = getVolume(X, Y, Z) - print v \ No newline at end of file diff --git a/SimPEG/EldadsCode/inv3X3BlockDiagonal.py b/SimPEG/EldadsCode/inv3X3BlockDiagonal.py deleted file mode 100644 index 2bfd86c4..00000000 --- a/SimPEG/EldadsCode/inv3X3BlockDiagonal.py +++ /dev/null @@ -1,36 +0,0 @@ -from utils import * -from sputils import * - - -def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33): - - a11 = mkvc(a11) - a12 = mkvc(a12) - a13 = mkvc(a13) - a21 = mkvc(a21) - a22 = mkvc(a22) - a23 = mkvc(a23) - a31 = mkvc(a31) - a32 = mkvc(a32) - a33 = mkvc(a33) - - detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32 - - b11 = +(a22*a33 - a23*a32)/detA - b12 = -(a12*a33 - a13*a32)/detA - b13 = +(a12*a23 - a13*a22)/detA - - b21 = +(a31*a23 - a21*a33)/detA - b22 = -(a31*a13 - a11*a33)/detA - b23 = +(a21*a13 - a11*a23)/detA - - b31 = -(a31*a22 - a21*a32)/detA - b32 = +(a31*a12 - a11*a32)/detA - b33 = -(a21*a12 - a11*a22)/detA - - B = appendBottom3( - appendRight3(sdiag(b11), sdiag(b12), sdiag(b13)), - appendRight3(sdiag(b21), sdiag(b22), sdiag(b23)), - appendRight3(sdiag(b31), sdiag(b32), sdiag(b33))) - - return B diff --git a/SimPEG/EldadsCode/meshUtils.py b/SimPEG/EldadsCode/meshUtils.py deleted file mode 100644 index 5d3f071d..00000000 --- a/SimPEG/EldadsCode/meshUtils.py +++ /dev/null @@ -1,94 +0,0 @@ -from sputils import * -from utils import * -from numpy import * - - -def getCellCenterFromNodal(X, Y, Z): - """Cell Centers from Nodal locations""" - XC = 1.0/8.0 * (X[:-1, :-1, :-1] + X[1:, :-1, :-1] + X[:-1, 1:, :-1] + X[1:, 1:, :-1] + - X[:-1, :-1, 1:] + X[1:, :-1, 1:] + X[:-1, 1:, 1:] + X[1:, 1:, 1:]) - - YC = 1.0/8.0 * (Y[:-1, :-1, :-1] + Y[1:, :-1, :-1] + Y[:-1, 1:, :-1] + Y[1:, 1:, :-1] + - Y[:-1, :-1, 1:] + Y[1:, :-1, 1:] + Y[:-1, 1:, 1:] + Y[1:, 1:, 1:]) - - ZC = 1.0/8.0 * (Z[:-1, :-1, :-1] + Z[1:, :-1, :-1] + Z[:-1, 1:, :-1] + Z[1:, 1:, :-1] + - Z[:-1, :-1, 1:] + Z[1:, :-1, 1:] + Z[:-1, 1:, 1:] + Z[1:, 1:, 1:]) - - return (XC, YC, ZC) - - -def getEdgesFromNodal(X, Y, Z): - """Edges from Nodal locations - - 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) - """ - - XE1 = (X[1:, :, :]+X[:-1, :, :])/2.0 - YE1 = (Y[1:, :, :]+Y[:-1, :, :])/2.0 - ZE1 = (Z[1:, :, :]+Z[:-1, :, :])/2.0 - - XE2 = (X[:, 1:, :]+X[:, :-1, :])/2.0 - YE2 = (Y[:, 1:, :]+Y[:, :-1, :])/2.0 - ZE2 = (Z[:, 1:, :]+Z[:, :-1, :])/2.0 - - XE3 = (X[:, :, 1:]+X[:, :, :-1])/2.0 - YE3 = (Y[:, :, 1:]+Y[:, :, :-1])/2.0 - ZE3 = (Z[:, :, 1:]+Z[:, :, :-1])/2.0 - - return (XE1, YE1, ZE1, XE2, YE2, ZE2, XE3, YE3, ZE3) - - -def getFacesFromNodal(X, Y, Z): - """Get faces from nodal --""" - - XF1 = 1.0/4.0*(X[:, :-1, :-1]+X[:, 1:, :-1]+X[:, :-1, 1:]+X[:, 1:, 1:]) - YF1 = 1.0/4.0*(Y[:, :-1, :-1]+Y[:, 1:, :-1]+Y[:, :-1, 1:]+Y[:, 1:, 1:]) - ZF1 = 1.0/4.0*(Z[:, :-1, :-1]+Z[:, 1:, :-1]+Z[:, :-1, 1:]+Z[:, 1:, 1:]) - - XF2 = 1.0/4.0*(X[:-1, :, :-1]+X[1:, :, :-1]+X[:-1, :, 1:]+X[1:, :, 1:]) - YF2 = 1.0/4.0*(Y[:-1, :, :-1]+Y[1:, :, :-1]+Y[:-1, :, 1:]+Y[1:, :, 1:]) - ZF2 = 1.0/4.0*(Z[:-1, :, :-1]+Z[1:, :, :-1]+Z[:-1, :, 1:]+Z[1:, :, 1:]) - - XF3 = 1.0/4.0*(X[:-1, :-1, :]+X[1:, :-1, :]+X[:-1, 1:, :]+X[1:, 1:, :]) - YF3 = 1.0/4.0*(Y[:-1, :-1, :]+Y[1:, :-1, :]+Y[:-1, 1:, :]+Y[1:, 1:, :]) - ZF3 = 1.0/4.0*(Z[:-1, :-1, :]+Z[1:, :-1, :]+Z[:-1, 1:, :]+Z[1:, 1:, :]) - - return (XF1, YF1, ZF1, XF2, YF2, ZF2, XF3, YF3, ZF3) - - -def projectEdgeVectorField(EV1, EV2, EV3, X, Y, Z): - """Project Edge vector field""" - - t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z, nrm1, nrm2, nrm3 = getEdgeTangent(X, Y, Z) - - E1 = EV1[:, 0]*mkvc(t1x) + EV1[:, 1]*mkvc(t1y) + EV1[:, 2]*mkvc(t1z) - E2 = EV2[:, 0]*mkvc(t2x) + EV2[:, 1]*mkvc(t2y) + EV2[:, 2]*mkvc(t2z) - E3 = EV3[:, 0]*mkvc(t3x) + EV3[:, 1]*mkvc(t3y) + EV3[:, 2]*mkvc(t3z) - - return hstack((hstack((mkvc(E1), mkvc(E2))), mkvc(E3))) - - -def projectFaceVectorField(FV1, FV2, FV3, X, Y, Z): - """Prolect Face vector field""" - - n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, ar1, ar2, ar3 = getFaceNormals(X, Y, Z) - - F1 = FV1[:, 0]*mkvc(n1x) + FV1[:, 1]*mkvc(n1y) + FV1[:, 2]*mkvc(n1z) - F2 = FV2[:, 0]*mkvc(n2x) + FV2[:, 1]*mkvc(n2y) + FV2[:, 2]*mkvc(n2z) - F3 = FV3[:, 0]*mkvc(n3x) + FV3[:, 1]*mkvc(n3y) + FV3[:, 2]*mkvc(n3z) - - return hstack((hstack((mkvc(F1), mkvc(F2))), mkvc(F3))) diff --git a/SimPEG/EldadsCode/sputils.py b/SimPEG/EldadsCode/sputils.py deleted file mode 100644 index bf54b15a..00000000 --- a/SimPEG/EldadsCode/sputils.py +++ /dev/null @@ -1,77 +0,0 @@ -from scipy import sparse -from numpy import * - - -def ddx(n): - """Define 1D derivatives""" - # ddx = lambda n: sparse.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format='csr') - return sparse.spdiags(-ones(n), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1) - - -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 sdiag(h): - """Diagonal matrix""" - return sparse.spdiags(h, 0, size(h), size(h)) - - -def speye(n): - """sparse identity""" - return sparse.spdiags(ones(n), 0, n, n) - - -def kron3(A, B, C): - """two kron prods""" - return sparse.kron(sparse.kron(A, B), C) - - -def appendBottom(A, B): - """append on bottom""" - C = sparse.vstack((A, B)) - C = C.tocsr() - return C - - -def appendBottom3(A, B, C): - """append on bottom""" - C = appendBottom(appendBottom(A, B), C) - C = C.tocsr() - return C - - -def appendRight(A, B): - """append on right""" - C = sparse.hstack((A, B)) - C = C.tocsr() - return C - - -def appendRight3(A, B, C): - """append on right""" - C = appendRight(appendRight(A, B), C) - C = C.tocsr() - return C - - -def blkDiag(A, B): - """blockdigonal""" - O12 = sparse.coo_matrix((shape(A)[0], shape(B)[1])) - O21 = sparse.coo_matrix((shape(B)[0], shape(A)[1])) - C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B)))) - C = C.tocsr() - return C - - -def blkDiag3(A, B, C): - """blockdigonal 3""" - ABC = blkDiag(blkDiag(A, B), C) - ABC = ABC.tocsr() - return ABC - - -def spzeros(n1, n2): - """spzeros""" - return sparse.coo_matrix((n1, n2)) diff --git a/SimPEG/EldadsCode/tools.py b/SimPEG/EldadsCode/tools.py deleted file mode 100644 index 132962a7..00000000 --- a/SimPEG/EldadsCode/tools.py +++ /dev/null @@ -1,222 +0,0 @@ -import numpy; -import cmath; -import math; - -def prod(arg): - """ returns the product of elements in arg. - arg can be list, tuple, set, and array with numerical values. """ - ret = 1; - for i in range(0,len(arg)): - ret = ret * arg[i]; - return ret; - - -def allIndices(dim): - """ From the given shape of dimenions (e.g. (2,3,4)), - generate a numpy.array of all, sorted indices.""" - - length = len(dim); - - sub = numpy.arange(dim[length-1]).reshape(dim[length-1],1); - - for d in range(length-2, -1, -1): - for i in range(0, dim[d]): - temp = numpy.ndarray([len(sub), 1]); - temp.fill(i); - temp = numpy.concatenate((temp,sub), axis=1); - if(i == 0): - newsub = temp; - else: - newsub = numpy.concatenate((newsub, temp), axis = 0); - - sub = newsub; - - return sub; - -def find(nda, obj): - """returns the index of the obj in the given nda(ndarray, list, or tuple)""" - for i in range(0, len(nda)): - if(nda[i] == obj): - return i; - return -1; - - -def notin(n, vector): - """returns a numpy.array object that contains - elements in [0,1, ... n-1] but not in vector.""" - ret = numpy.arange(n).tolist(); - for i in vector: - if (0 <= i and i < n): - ret.remove(i); - return numpy.array(ret); - - - -def getelts(nda, indices): - """From the given nda(ndarray, list, or tuple), returns the list located at the given indices""" - ret = []; - for i in indices: - ret.extend([nda[i]]); - return numpy.array(ret); - -def sub2ind(shape, subs): - """ From the given shape, returns the index of the given subscript""" - revshp = list(shape); - revshp.reverse(); - mult = [1]; - for i in range(0, len(revshp)-1): - mult.extend([mult[i]*revshp[i]]); - mult.reverse(); - mult = numpy.array(mult).reshape(len(mult),1); - - idx = numpy.dot((subs) , (mult)); - return idx; - -def ind2sub(shape, ind): - """ From the given shape, returns the subscrips of the given index""" - revshp = []; - revshp.extend(shape); - revshp.reverse(); - mult = [1]; - for i in range(0, len(revshp)-1): - mult.extend([mult[i]*revshp[i]]); - mult.reverse(); - mult = numpy.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 tt_dimscehck(dims, N, M = None, exceptdims = False): - """Checks whether the specified dimensions are valid in a tensor of N-dimension. - If M is given, then it will also retuns an index for M multiplicands. - If exceptdims == True, then it will compute for the dimensions not specified.""" - - # if exceptdims is true - if(exceptdims): - dims = listdiff(range(0,N), dims); - - #check vals in between 0 and N-1 - for i in range(0, len(dims)): - if(dims[i] < 0 or dims[i] >= N): - raise ValueError("invalid dimensions specified"); - - # number of dimensions in dims - p = len(dims); - - sdims = []; - sdims.extend(dims); - sdims.sort(); - - #indices of the elements in the sorted array - sidx = []; - #table that denotes whether the index is used - table = numpy.ndarray([len(sdims)]); - table.fill(0); - - for i in range(0, len(sdims)): - for j in range(0, len(dims)): - if(sdims[i] == dims[j] and table[j] == 0): - sidx.extend([j]); - table[j] = 1; - break; - - if (M == None): - return sdims - - if(M > N): - raise ValueError("Cannot have more multiplicands than dimensions"); - - if(M != N and M != p): - raise ValueError("invalid number of multiplicands"); - - if(M == p): - vidx = sidx; - else: - vidx = sdims; - - return (sdims, vidx); - -def listtimes(list, c): - """multiplies the elements in the list by the given scalar value c""" - ret = [] - for i in range(0, len(list)): - ret.extend([list[i]]*c); - return ret; - -def listdiff(list1, list2): - """returns the list of elements that are in list 1 but not in list2""" - if(list1.__class__ == numpy.ndarray): - list1 = list1.tolist(); - if(list2.__class__ == numpy.ndarray): - list2 = list2.tolist(); - ret = [] - for i in range(0,len(list1)): - ok = true - for j in range(0, len(list2)): - if(list[i] == list[j]): - ok = false; - break; - if(ok): - ret.extend([list[i]]); - return ret; - - - -def tt_subscheck(subs): - """Check whether the given list of subscripts are valid. Used for sptensor""" - isOk = True; - if(subs.size == 0): - isOk = True; - - elif(subs.ndim != 2): - isOk = False; - - else: - for i in range(0, (subs.size / subs[0].size)): - for j in range(0, (subs[0].size)): - val = subs[i][j]; - if( cmath.isnan(val) or cmath.isinf(val) or val < 0 or val != round(val) ): - isOk = False; - - if(not isOk): - raise ValueError("Subscripts must be a matrix of non-negative integers"); - - return isOk; - - -def tt_valscheck(vals): - """Check whether the given list of values are valid. Used for sptensor""" - isOk = True; - - if(vals.size == 0): - isOk = True; - - elif(vals.ndim != 2 or vals[0].size != 1): - isOk = False; - - if(not isOk): - raise ValueError("values must be a column array"); - - return isOk; - -def tt_sizecheck(size): - """Check whether the given size is valid. Used for sptensor""" - size = numpy.array(size); - isOk = True; - - if(size.ndim != 1): - isOk = False; - else: - for i in range(0, len(size)): - val = size[i]; - if(cmath.isnan(val) or cmath.isinf(val) - or val <= 0 or val != round(val)): - isOk = False; - - if(not isOk): - raise ValueError("size must be a row vector of real positive integers"); - return isOk; diff --git a/SimPEG/EldadsCode/utils.py b/SimPEG/EldadsCode/utils.py deleted file mode 100644 index e912c5b8..00000000 --- a/SimPEG/EldadsCode/utils.py +++ /dev/null @@ -1,87 +0,0 @@ -from numpy import * -import numpy as np - - -def diff(A, d): - if(d == 1): - return A[1:, :, :] - A[:-1, :, :] - elif(d == 2): - return A[:, 1:, :] - A[:, :-1, :] - else: - return A[:, :, 1:] - A[:, :, :-1] - #else: - # print('d must be 1,2 or 3') - - -def diffp(A, d1, d2): - if(d1 == 1 and d2 == 2): - return A[1:, 1:, :] - A[:-1, :-1, :] - elif(d1 == 1 and d2 == 3): - return A[1:, :, 1:] - A[:-1, :, :-1] - else: - return A[:, 1:, 1:] - A[:, :-1, :-1] - - -def diffm(A, d1, d2): - if(d1 == 3 and d2 == 2): - return A[:, :-1, 1:] - A[:, 1:, :-1] - elif(d1 == 1 and d2 == 3): - return A[1:, :, :-1] - A[:-1, :, 1:] - elif(d1 == 2 and d2 == 1): - return A[:-1, 1:, :] - A[1:, :-1, :] - else: - print('d must be 1, 2 or 3') - - -def ave(A, d): - if(d == 1): - return 0.5*(A[1:, :, :] + A[:-1, :, :]) - elif(d == 2): - return 0.5*(A[:, 1:, :] + A[:, :-1, :]) - elif(d == 3): - return 0.5*(A[:, :, 1:] + A[:, :, :-1]) - else: - print('d must be 1,2 or 3') - - -def mkmat(x): - return reshape(matrix(x), (size(x), 1), 'F') - - -def hstack3(a, b, c): - a = mkvc(a) - b = mkvc(b) - c = mkvc(c) - a = mkmat(a) - b = mkmat(b) - c = mkmat(c) - return hstack((hstack((a, b)), c)) - - -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 diff --git a/SimPEG/EldadsCode/zevel.py b/SimPEG/EldadsCode/zevel.py deleted file mode 100644 index c3e9ae2d..00000000 --- a/SimPEG/EldadsCode/zevel.py +++ /dev/null @@ -1,119 +0,0 @@ -#============= Nodal Gradients =========================== -def getNodalGradient(h1,h2,h3): - - n1 = size(h1) - n2 = size(h2) - n3 = size(h3) - D1 = kron3(speye(n3+1),speye(n2+1),ddx(n1)) - D2 = kron3(speye(n3+1),ddx(n2),speye(n1+1)) - D3 = kron3(ddx(n3),speye(n2+1),speye(n1+1)) - - # topological gradient - GRAD = appendBottom3(D1,D2,D3) - - # scale for non-uniform mesh - L = blkDiag3(kron3(speye(n3+1),speye(n2+1),sdiag(1/h1)), - kron3(speye(n3+1),sdiag(1/h2),speye(n1+1)), - kron3(sdiag(1/h3),speye(n2+1),speye(n1+1))) - - return L*GRAD - -#============= Edge CURL =========================== -def getCurlMatrix(h1,h2,h3): - - n1 = size(h1) - n2 = size(h2) - n3 = size(h3) - - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - # derivatives on x-edge variables - D32 = kron3(d3,speye(n2),speye(n1+1)) - D23 = kron3(speye(n3),d2,speye(n1+1)) - D31 = kron3(d3,speye(n2+1),speye(n1)) - D13 = kron3(speye(n3),speye(n2+1),d1) - D21 = kron3(speye(n3+1),d2,speye(n1)) - D12 = kron3(speye(n3+1),speye(n2),d1) - - O1 = spzeros(shape(D32)[0],shape(D31)[1]) - O2 = spzeros(shape(D31)[0],shape(D32)[1]) - O3 = spzeros(shape(D21)[0],shape(D13)[1]) - - CURL = appendBottom3( - appendRight3(O1, -D32, D23), - appendRight3(D31, O2, -D13), - appendRight3(-D21, D12, O3)) - - # scale for non-uniform mesh - F = blkDiag3(kron3(sdiag(1/h3),sdiag(1/h2),speye(n1+1)), - kron3(sdiag(1/h3),speye(n2+1),sdiag(1/h1)), - kron3(speye(n3+1),sdiag(1/h2),sdiag(1/h1))) - - L = blkDiag3(kron3(speye(n3+1),speye(n2+1),sdiag(h1)), - kron3(speye(n3+1),sdiag(h2),speye(n1+1)), - kron3(sdiag(h3),speye(n2+1),speye(n1+1))) - - - return F*(CURL*L) - -#============= Face DIV =========================== -def getDivMatrix(h1,h2,h3): - - n1 = size(h1) - n2 = size(h2) - n3 = size(h3) - - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3),speye(n2),d1) - D2 = kron3(speye(n3),d2,speye(n1)) - D3 = kron3(d3,speye(n2),speye(n1)) - - # divergence on faces - D = appendRight3(D1, D2, D3) - - # scale for non-uniform mesh - F = blkDiag3(kron3(sdiag(h3),sdiag(h2),speye(n1+1)), - kron3(sdiag(h3),speye(n2+1),sdiag(h1)), - kron3(speye(n3+1),sdiag(h2),sdiag(h1))) - - V = kron3(sdiag(1/h3),sdiag(1/h2),sdiag(1/h1)) - - return V*(D*F) - -#====== Face Averageing ================= -def getFaceAverage(n1,n2,n3): - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - Af = appendRight3(kron3(speye(n3),speye(n2),av1), - kron3(speye(n3),av2,speye(n1)), - kron3(av3,speye(n2),speye(n1))) - return Af - -#====== Edge Averageing ================= -def getEdgeAverage(n1,n2,n3): - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - Ae = appendRight3(kron3(av3,av2,speye(n1)), - kron3(av3,speye(n2),av1), - kron3(speye(n3),av2,av1)) - return Ae - -#====== Node Averageing ================= -def getNodeAverage(n1,n2,n3): - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - return kron3(av3,av2,av1) - - diff --git a/SimPEG/EldadsCode/GaussNewton.py b/SimPEG/GaussNewton.py similarity index 100% rename from SimPEG/EldadsCode/GaussNewton.py rename to SimPEG/GaussNewton.py