From 125144002120f41a9dcf5a2be088d65a0f2ee861 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 11 Jun 2013 16:18:15 +0200 Subject: [PATCH] Initial Cleanup of Code. No major changes. (duplicate code, and python conventions) --- code/MFDdriver.py | 29 +++--- code/getCellVolume.py | 78 ++++++++-------- code/getDiffOps.py | 137 +++++++++++++++-------------- code/getEdgeInnerProduct.py | 171 ++++++++++++++++++------------------ code/inv3X3BlockDiagonal.py | 19 ++-- code/meshUtils.py | 167 +++++++++++++++++------------------ code/ndgrid.py | 29 ------ code/sputils.py | 108 ++++++++++++----------- code/utils.py | 163 +++++++++++++++++----------------- 9 files changed, 436 insertions(+), 465 deletions(-) delete mode 100644 code/ndgrid.py diff --git a/code/MFDdriver.py b/code/MFDdriver.py index c0b8edb4..16d18fe3 100644 --- a/code/MFDdriver.py +++ b/code/MFDdriver.py @@ -1,25 +1,18 @@ -# from scipy.sparse import linalg -from numpy import * -#from numpy.linalg import * +import numpy as np from numpy.random import randn -from utils import * +from utils import ndgrid from getDiffOps import getCurlMatrix, getNodalGradient -from sputils import * -from meshUtils import * from getFaceInnerProduct import getFaceInnerProduct from getEdgeInnerProduct import getEdgeInnerProduct -#from scipy.sparse.linalg import spsolve -from scipy.sparse.linalg import * -from pylab import * +from scipy.sparse.linalg import dsolve +from pylab import norm -n1 = 14 -n2 = 14 -n3 = 15 +n = np.array([14, 14, 15]) -X, Y, Z = ndgrid(linspace(0, 1, n1), linspace(0, 1, n2), linspace(0, 1, n3)) -sigma = 1e-2*ones([n1-1, n2-1, n3-1]) -sigma[:, :, (n3-1)/2:] = 1e-6 -mu = 4*pi*1e-7*ones([n1-1, n2-1, n3-1]) +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) @@ -29,8 +22,8 @@ Me = getEdgeInnerProduct(X, Y, Z, sigma) A = CURL.T * Mf * CURL + 1j * w * Me -ne = shape(A) -b = matrix(randn(ne[0])).T +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 diff --git a/code/getCellVolume.py b/code/getCellVolume.py index 95d01427..bb536842 100644 --- a/code/getCellVolume.py +++ b/code/getCellVolume.py @@ -1,58 +1,58 @@ -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 * -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]) +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): +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 ) + 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 + 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) \ No newline at end of file + v = getCellVolume(X, Y, Z) + print v diff --git a/code/getDiffOps.py b/code/getDiffOps.py index 40f8239e..dbb4e5b0 100644 --- a/code/getDiffOps.py +++ b/code/getDiffOps.py @@ -1,5 +1,3 @@ -from scipy.sparse import linalg -from scipy import sparse from sputils import * from utils import * from numpy import * @@ -8,96 +6,101 @@ from getCellVolume import getCellVolume from getFaceNormals import getFaceNormals -#============= Face DIV =========================== -def getDivMatrix(X,Y,Z): +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) - + 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)) - + 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 - - #============= Edge CURL =========================== - -def getCurlMatrix(X,Y,Z): - + + +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) + 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]) + 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) -#============= Nodal Gradients =========================== -def getNodalGradient(X,Y,Z): - + # 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)) - + 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) - + 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) + 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) + 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) - G = getNodalGradient(X,Y,Z) - tt = C*G - print(tt) \ No newline at end of file + print(tt) diff --git a/code/getEdgeInnerProduct.py b/code/getEdgeInnerProduct.py index 9391b425..776f31f7 100644 --- a/code/getEdgeInnerProduct.py +++ b/code/getEdgeInnerProduct.py @@ -8,203 +8,202 @@ from getEdgeTangent import * from inv3X3BlockDiagonal import * from getCellVolume import getCellVolume -# [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 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) -def getEdgeInnerProduct(X,Y,Z,sigma): - m = array(shape(X))-1 + 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,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)) + 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)) + 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)) + 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)) + 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)) + 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)) + 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)) + 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)) + 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 = 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)) diff --git a/code/inv3X3BlockDiagonal.py b/code/inv3X3BlockDiagonal.py index cf523ddd..2bfd86c4 100644 --- a/code/inv3X3BlockDiagonal.py +++ b/code/inv3X3BlockDiagonal.py @@ -1,12 +1,11 @@ -from scipy.sparse import linalg from utils import * from sputils import * -def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,a33): +def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33): - a11 = mkvc(a11) - a12 = mkvc(a12) + a11 = mkvc(a11) + a12 = mkvc(a12) a13 = mkvc(a13) a21 = mkvc(a21) a22 = mkvc(a22) @@ -17,16 +16,16 @@ def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,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 + b11 = +(a22*a33 - a23*a32)/detA b12 = -(a12*a33 - a13*a32)/detA - b13 = (a12*a23 - a13*a22)/detA + b13 = +(a12*a23 - a13*a22)/detA - b21 = (a31*a23 - a21*a33)/detA + b21 = +(a31*a23 - a21*a33)/detA b22 = -(a31*a13 - a11*a33)/detA - b23 = (a21*a13 - a11*a23)/detA + b23 = +(a21*a13 - a11*a23)/detA b31 = -(a31*a22 - a21*a32)/detA - b32 = (a31*a12 - a11*a32)/detA + b32 = +(a31*a12 - a11*a32)/detA b33 = -(a21*a12 - a11*a22)/detA B = appendBottom3( @@ -34,4 +33,4 @@ def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,a33): appendRight3(sdiag(b21), sdiag(b22), sdiag(b23)), appendRight3(sdiag(b31), sdiag(b32), sdiag(b33))) - return B \ No newline at end of file + return B diff --git a/code/meshUtils.py b/code/meshUtils.py index fbf5cc30..5d3f071d 100644 --- a/code/meshUtils.py +++ b/code/meshUtils.py @@ -1,101 +1,94 @@ -from scipy.sparse import linalg -from scipy import sparse from sputils import * from utils import * from numpy import * -#----- Cell Centers from Nodal locations ----- -def getCellCenterFromNodal(X,Y,Z): - XC = 1.0/8.0 * (X[0:-1,0:-1,0:-1] + X[1:,0:-1,0:-1] + X[0:-1,1:,0:-1] + X[1:,1:,0:-1] + - X[0:-1,0:-1,1:] + X[1:,0:-1,1:] + X[0:-1,1:,1:] + X[1:,1:,1:]) +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[0:-1,0:-1,0:-1] + Y[1:,0:-1,0:-1] + Y[0:-1,1:,0:-1] + Y[1:,1:,0:-1] + - Y[0:-1,0:-1,1:] + Y[1:,0:-1,1:] + Y[0:-1,1:,1:] + Y[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[0:-1,0:-1,0:-1] + Z[1:,0:-1,0:-1] + Z[0:-1,1:,0:-1] + Z[1:,1:,0:-1] + - Z[0:-1,0:-1,1:] + Z[1:,0:-1,1:] + Z[0:-1,1:,1:] + Z[1:,1:,1:]) - - return (XC,YC,ZC) + 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:]) -#----- Edges from Nodal locations ----- - -def getEdgesFromNodal(X,Y,Z): -# -# 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[0:-1,:,:])/2.0 - YE1 = (Y[1:,:,:]+Y[0:-1,:,:])/2.0 - ZE1 = (Z[1:,:,:]+Z[0:-1,:,:])/2.0 - - XE2 = (X[:,1:,:]+X[:,0:-1,:])/2.0 - YE2 = (Y[:,1:,:]+Y[:,0:-1,:])/2.0 - ZE2 = (Z[:,1:,:]+Z[:,0:-1,:])/2.0 - - XE3 = (X[:,:,1:]+X[:,:,0:-1])/2.0 - YE3 = (Y[:,:,1:]+Y[:,:,0:-1])/2.0 - ZE3 = (Z[:,:,1:]+Z[:,:,0:-1])/2.0 - - return (XE1,YE1,ZE1,XE2,YE2,ZE2,XE3,YE3,ZE3) - - #-- Get faces from nodal -- - -def getFacesFromNodal(X,Y,Z): - - XF1 = 1.0/4.0*(X[:,0:-1,0:-1]+X[:,1:,0:-1]+X[:,0:-1,1:]+X[:,1:,1:]) - YF1 = 1.0/4.0*(Y[:,0:-1,0:-1]+Y[:,1:,0:-1]+Y[:,0:-1,1:]+Y[:,1:,1:]) - ZF1 = 1.0/4.0*(Z[:,0:-1,0:-1]+Z[:,1:,0:-1]+Z[:,0:-1,1:]+Z[:,1:,1:]) - - XF2 = 1.0/4.0*(X[0:-1,:,0:-1]+X[1:,:,0:-1]+X[0:-1,:,1:]+X[1:,:,1:]) - YF2 = 1.0/4.0*(Y[0:-1,:,0:-1]+Y[1:,:,0:-1]+Y[0:-1,:,1:]+Y[1:,:,1:]) - ZF2 = 1.0/4.0*(Z[0:-1,:,0:-1]+Z[1:,:,0:-1]+Z[0:-1,:,1:]+Z[1:,:,1:]) - - XF3 = 1.0/4.0*(X[0:-1,0:-1,:]+X[1:,0:-1,:]+X[0:-1,1:,:]+X[1:,1:,:]) - YF3 = 1.0/4.0*(Y[0:-1,0:-1,:]+Y[1:,0:-1,:]+Y[0:-1,1:,:]+Y[1:,1:,:]) - ZF3 = 1.0/4.0*(Z[0:-1,0:-1,:]+Z[1:,0:-1,:]+Z[0:-1,1:,:]+Z[1:,1:,:]) - - return (XF1,YF1,ZF1,XF2,YF2,ZF2,XF3,YF3,ZF3) - - #-- Project Edge vector field - -def projectEdgeVectorField(EV1,EV2,EV3,X,Y,Z): - - t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,nrm1,nrm2,nrm3 = getEdgeTangent(X,Y,Z) + return (XC, YC, ZC) - 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))) - - -#-- Prolect Face vector field +def getEdgesFromNodal(X, Y, Z): + """Edges from Nodal locations -def projectFaceVectorField(FV1,FV2,FV3,X,Y,Z): + 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) + """ - n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,ar1,ar2,ar3 = getFaceNormals(X,Y,Z) + XE1 = (X[1:, :, :]+X[:-1, :, :])/2.0 + YE1 = (Y[1:, :, :]+Y[:-1, :, :])/2.0 + ZE1 = (Z[1:, :, :]+Z[:-1, :, :])/2.0 - 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))) - + XE2 = (X[:, 1:, :]+X[:, :-1, :])/2.0 + YE2 = (Y[:, 1:, :]+Y[:, :-1, :])/2.0 + ZE2 = (Z[:, 1:, :]+Z[:, :-1, :])/2.0 - - \ No newline at end of file + 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/code/ndgrid.py b/code/ndgrid.py deleted file mode 100644 index d051db8f..00000000 --- a/code/ndgrid.py +++ /dev/null @@ -1,29 +0,0 @@ -from numpy import * - -def ndgrid(x,y,z): - - n1 = size(x) - n2 = size(y) - n3 = size(z) - X = zeros([n1,n2,n3]) - Y = zeros([n1,n2,n3]) - Z = zeros([n1,n2,n3]) - for i in range(0, n2): - for j in range(0,n3): - X[:,i,j] = x - - for i in range(0, n1): - for j in range(0,n3): - Y[i,:,j] = y - - for i in range(0, n1): - for j in range(0,n2): - Z[i,j,:] = z - - - return (X,Y,Z) - - -if __name__ == '__main__': - - X = ndgrid([1,2,3],[2,4,5,6],[4,6,7,8]) diff --git a/code/sputils.py b/code/sputils.py index d5a83145..bf54b15a 100644 --- a/code/sputils.py +++ b/code/sputils.py @@ -1,67 +1,77 @@ -from scipy.sparse import linalg from scipy import sparse from numpy import * -#======== Define 1D derivatives ============= def ddx(n): - return sparse.spdiags(-ones(n),0,n,n+1) + sparse.spdiags(ones(n+1),1,n,n+1) + """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) + -#======== Define 1D average ============= def av(n): - return 0.5*(sparse.spdiags(ones(n+1),0,n,n+1) + sparse.spdiags(ones(n+1),1,n,n+1)) + """Define 1D average""" + return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)) + -#======== Diagonal matrix ============= def sdiag(h): - return sparse.spdiags(h,0,size(h),size(h)) + """Diagonal matrix""" + return sparse.spdiags(h, 0, size(h), size(h)) + -#======== sparse identity ============= def speye(n): - return sparse.spdiags(ones(n),0,n,n) + """sparse identity""" + return sparse.spdiags(ones(n), 0, n, n) -#======== two kron prods ============= -def kron3(A,B,C): - return sparse.kron(sparse.kron(A,B),C) - -#======== append on bottom ============= -def appendBottom(A,B): - C = sparse.vstack((A,B)) - C = C.tocsr() - return C - -#======== append on bottom ============= -def appendBottom3(A,B,C): - C = appendBottom(appendBottom(A,B),C) - C = C.tocsr() - return C -#======== append on right ============= -def appendRight(A,B): - C = sparse.hstack((A,B)) - C = C.tocsr() - return C +def kron3(A, B, C): + """two kron prods""" + return sparse.kron(sparse.kron(A, B), C) -#======== append on right ============= -def appendRight3(A,B,C): - C = appendRight(appendRight(A,B),C) - C = C.tocsr() - return C -#======== blockdigonal ============= -def blkDiag(A,B): - 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)))) +def appendBottom(A, B): + """append on bottom""" + C = sparse.vstack((A, B)) C = C.tocsr() - return C - -#======== blockdigonal 3 ============= -def blkDiag3(A,B,C): - ABC = blkDiag(blkDiag(A,B),C) + 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 - -#======== spzeros ============= -def spzeros(n1,n2): - return sparse.coo_matrix((n1,n2)) - + + +def spzeros(n1, n2): + """spzeros""" + return sparse.coo_matrix((n1, n2)) diff --git a/code/utils.py b/code/utils.py index ebd028e3..7ea54f17 100644 --- a/code/utils.py +++ b/code/utils.py @@ -1,124 +1,127 @@ from numpy import * -def diff(A,d): - end = -1 - if(d==1): - return A[1:,0:,0:] - A[0:end,0:,0:] - elif(d==2): - return A[0:,1:,0:] - A[0:,0:end,0:] +def diff(A, d): + if(d == 1): + return A[1:, :, :] - A[:-1, :, :] + elif(d == 2): + return A[:, 1:, :] - A[:, :-1, :] else: - return A[0:,0:,1:] - A[0:,0:,0:end] + return A[:, :, 1:] - A[:, :, :-1] #else: # print('d must be 1,2 or 3') def diffp(A, d1, d2): - end = -1 - if(d1 == 1 and d2 == 2 ): - return A[1:,1:, 0:] - A[0:end,0:end,0:] + if(d1 == 1 and d2 == 2): + return A[1:, 1:, :] - A[:-1, :-1, :] elif(d1 == 1 and d2 == 3): - return A[1:,0:,1:] - A[0:end,0:,0:end] + return A[1:, :, 1:] - A[:-1, :, :-1] else: - return A[0:,1:,1:] - A[0:,0:end,0:end] + return A[:, 1:, 1:] - A[:, :-1, :-1] def diffm(A, d1, d2): - end = -1 - if(d1 == 3 and d2 == 2 ): - return A[:,0:end,1:] - A[:,1:,0:end] + if(d1 == 3 and d2 == 2): + return A[:, :-1, 1:] - A[:, 1:, :-1] elif(d1 == 1 and d2 == 3): - return A[1:, :, 0:end] - A[0:end,:,1:] + return A[1:, :, :-1] - A[:-1, :, 1:] elif(d1 == 2 and d2 == 1): - return A[0:end, 1:, :] - A[1:, 0:end, :] + 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 ave(A,d): - end = 0 - if(d==1): - return 0.5*(A[1:,:,:] + A[0:end-1,0:,:]) - elif(d==2): - return 0.5*(A[:,1:,:] + A[0:,0:end-1,:]) - elif(d==3): - return 0.5*(A[:,:,1:] + A[0:,0:,0:end-1]) - else: - print('d must be 1,2 or 3') - -def reshapeF(sp,d): - return reshape(sp,d,'F') -def mkvc(A): - return reshape(A,[size(A),1],'F').flatten() - - -def ndgrid(x,y,z): - +def reshapeF(sp, d): + return reshape(sp, d, 'F') + + +def mkvc(A): + return reshape(A, [size(A), 1], 'F').flatten() + + +def ndgrid(x, y, z): + n1 = size(x) n2 = size(y) n3 = size(z) - X = zeros([n1,n2,n3]) - Y = zeros([n1,n2,n3]) - Z = zeros([n1,n2,n3]) + X = zeros([n1, n2, n3]) + Y = zeros([n1, n2, n3]) + Z = zeros([n1, n2, n3]) for i in range(0, n2): - for j in range(0,n3): - X[:,i,j] = x - + for j in range(0, n3): + X[:, i, j] = x + for i in range(0, n1): - for j in range(0,n3): - Y[i,:,j] = y - + for j in range(0, n3): + Y[i, :, j] = y + for i in range(0, n1): - for j in range(0,n2): - Z[i,j,:] = z - - - return (X,Y,Z) - + for j in range(0, n2): + Z[i, j, :] = z + + return (X, Y, Z) + def ind2sub(shape, ind): # From the given shape, returns the subscrips of the given index - revshp = []; - revshp.extend(shape); - mult = [1]; + 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; + 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]; + 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; + mult.extend([mult[i]*revshp[i]]) + mult = array(mult).reshape(len(mult), 1) + + idx = dot((subs), (mult)) + return idx 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)) - - + 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)) if __name__ == '__main__': X, Y, Z = mgrid[0:4, 0:5, 0:6] + print Z + t = ave(X, 1) - \ No newline at end of file + print t