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