From 2554318427a0a4cdd0f65f735e7d434df2b10346 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 1 Jun 2013 14:18:31 -0700 Subject: [PATCH] Eldad's initial code for logically orthogonal mesh simulation. --- code/MFDdriver.py | 43 +++++++ code/getCellVolume.py | 58 ++++++++++ code/getDiffOps.py | 103 +++++++++++++++++ code/getEdgeInnerProduct.py | 214 ++++++++++++++++++++++++++++++++++ code/getEdgeTangent.py | 60 ++++++++++ code/getFaceInnerProduct.py | 86 ++++++++++++++ code/getFaceNormals.py | 73 ++++++++++++ code/getVolume.py | 34 ++++++ code/inv3X3BlockDiagonal.py | 37 ++++++ code/meshUtils.py | 101 ++++++++++++++++ code/ndgrid.py | 29 +++++ code/sputils.py | 67 +++++++++++ code/tools.py | 222 ++++++++++++++++++++++++++++++++++++ code/utils.py | 124 ++++++++++++++++++++ code/zevel.py | 119 +++++++++++++++++++ 15 files changed, 1370 insertions(+) create mode 100644 code/MFDdriver.py create mode 100644 code/getCellVolume.py create mode 100644 code/getDiffOps.py create mode 100644 code/getEdgeInnerProduct.py create mode 100644 code/getEdgeTangent.py create mode 100644 code/getFaceInnerProduct.py create mode 100644 code/getFaceNormals.py create mode 100644 code/getVolume.py create mode 100644 code/inv3X3BlockDiagonal.py create mode 100644 code/meshUtils.py create mode 100644 code/ndgrid.py create mode 100644 code/sputils.py create mode 100644 code/tools.py create mode 100644 code/utils.py create mode 100644 code/zevel.py diff --git a/code/MFDdriver.py b/code/MFDdriver.py new file mode 100644 index 00000000..c0b8edb4 --- /dev/null +++ b/code/MFDdriver.py @@ -0,0 +1,43 @@ +# from scipy.sparse import linalg +from numpy import * +#from numpy.linalg import * +from numpy.random import randn +from utils import * +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 * + +n1 = 14 +n2 = 14 +n3 = 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]) +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 = shape(A) +b = 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/code/getCellVolume.py b/code/getCellVolume.py new file mode 100644 index 00000000..95d01427 --- /dev/null +++ b/code/getCellVolume.py @@ -0,0 +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]) + + 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) \ No newline at end of file diff --git a/code/getDiffOps.py b/code/getDiffOps.py new file mode 100644 index 00000000..40f8239e --- /dev/null +++ b/code/getDiffOps.py @@ -0,0 +1,103 @@ +from scipy.sparse import linalg +from scipy import sparse +from sputils import * +from utils import * +from numpy import * +from getEdgeTangent import * +from getCellVolume import getCellVolume +from getFaceNormals import getFaceNormals + + +#============= Face DIV =========================== +def getDivMatrix(X,Y,Z): + + 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 + + #============= Edge CURL =========================== + +def getCurlMatrix(X,Y,Z): + + 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) + +#============= Nodal Gradients =========================== +def getNodalGradient(X,Y,Z): + + 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) \ No newline at end of file diff --git a/code/getEdgeInnerProduct.py b/code/getEdgeInnerProduct.py new file mode 100644 index 00000000..9391b425 --- /dev/null +++ b/code/getEdgeInnerProduct.py @@ -0,0 +1,214 @@ +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 + +# [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(X,Y,Z,sigma): + + 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/code/getEdgeTangent.py b/code/getEdgeTangent.py new file mode 100644 index 00000000..a8426df3 --- /dev/null +++ b/code/getEdgeTangent.py @@ -0,0 +1,60 @@ +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/code/getFaceInnerProduct.py b/code/getFaceInnerProduct.py new file mode 100644 index 00000000..2870eb88 --- /dev/null +++ b/code/getFaceInnerProduct.py @@ -0,0 +1,86 @@ +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/code/getFaceNormals.py b/code/getFaceNormals.py new file mode 100644 index 00000000..2038a301 --- /dev/null +++ b/code/getFaceNormals.py @@ -0,0 +1,73 @@ +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/code/getVolume.py b/code/getVolume.py new file mode 100644 index 00000000..3b4cd250 --- /dev/null +++ b/code/getVolume.py @@ -0,0 +1,34 @@ +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/code/inv3X3BlockDiagonal.py b/code/inv3X3BlockDiagonal.py new file mode 100644 index 00000000..cf523ddd --- /dev/null +++ b/code/inv3X3BlockDiagonal.py @@ -0,0 +1,37 @@ +from scipy.sparse import linalg +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 \ No newline at end of file diff --git a/code/meshUtils.py b/code/meshUtils.py new file mode 100644 index 00000000..fbf5cc30 --- /dev/null +++ b/code/meshUtils.py @@ -0,0 +1,101 @@ +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:]) + + 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:]) + + 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) + +#----- 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) + + + 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 projectFaceVectorField(FV1,FV2,FV3,X,Y,Z): + + 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))) + + + + \ No newline at end of file diff --git a/code/ndgrid.py b/code/ndgrid.py new file mode 100644 index 00000000..d051db8f --- /dev/null +++ b/code/ndgrid.py @@ -0,0 +1,29 @@ +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 new file mode 100644 index 00000000..d5a83145 --- /dev/null +++ b/code/sputils.py @@ -0,0 +1,67 @@ +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 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)) + +#======== Diagonal matrix ============= +def sdiag(h): + return sparse.spdiags(h,0,size(h),size(h)) + +#======== sparse identity ============= +def speye(n): + 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 + +#======== 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)))) + C = C.tocsr() + return C + +#======== blockdigonal 3 ============= +def blkDiag3(A,B,C): + ABC = blkDiag(blkDiag(A,B),C) + ABC = ABC.tocsr() + return ABC + +#======== spzeros ============= +def spzeros(n1,n2): + return sparse.coo_matrix((n1,n2)) + diff --git a/code/tools.py b/code/tools.py new file mode 100644 index 00000000..132962a7 --- /dev/null +++ b/code/tools.py @@ -0,0 +1,222 @@ +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/code/utils.py b/code/utils.py new file mode 100644 index 00000000..ebd028e3 --- /dev/null +++ b/code/utils.py @@ -0,0 +1,124 @@ +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:] + else: + return A[0:,0:,1:] - A[0:,0:,0:end] + #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:] + elif(d1 == 1 and d2 == 3): + return A[1:,0:,1:] - A[0:end,0:,0:end] + else: + return A[0:,1:,1:] - A[0:,0:end,0:end] + + +def diffm(A, d1, d2): + end = -1 + if(d1 == 3 and d2 == 2 ): + return A[:,0:end,1:] - A[:,1:,0:end] + elif(d1 == 1 and d2 == 3): + return A[1:, :, 0:end] - A[0:end,:,1:] + elif(d1 == 2 and d2 == 1): + return A[0:end, 1:, :] - A[1:, 0:end, :] + 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): + + 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) + + +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; + + +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)) + + + + +if __name__ == '__main__': + + X, Y, Z = mgrid[0:4, 0:5, 0:6] + + t = ave(X, 1) + \ No newline at end of file diff --git a/code/zevel.py b/code/zevel.py new file mode 100644 index 00000000..c3e9ae2d --- /dev/null +++ b/code/zevel.py @@ -0,0 +1,119 @@ +#============= 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) + +