mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-15 11:26:09 +08:00
Eldad's initial code for logically orthogonal mesh simulation.
This commit is contained in:
Rowan Cockett
@@ -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
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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)))
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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])
|
||||
@@ -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))
|
||||
|
||||
+222
@@ -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;
|
||||
+124
@@ -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)
|
||||
|
||||
+119
@@ -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)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user