mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 20:23:01 +08:00
Initial Cleanup of Code.
No major changes. (duplicate code, and python conventions)
This commit is contained in:
Rowan Cockett
+11
-18
@@ -1,25 +1,18 @@
|
||||
# from scipy.sparse import linalg
|
||||
from numpy import *
|
||||
#from numpy.linalg import *
|
||||
import numpy as np
|
||||
from numpy.random import randn
|
||||
from utils import *
|
||||
from utils import ndgrid
|
||||
from getDiffOps import getCurlMatrix, getNodalGradient
|
||||
from sputils import *
|
||||
from meshUtils import *
|
||||
from getFaceInnerProduct import getFaceInnerProduct
|
||||
from getEdgeInnerProduct import getEdgeInnerProduct
|
||||
#from scipy.sparse.linalg import spsolve
|
||||
from scipy.sparse.linalg import *
|
||||
from pylab import *
|
||||
from scipy.sparse.linalg import dsolve
|
||||
from pylab import norm
|
||||
|
||||
n1 = 14
|
||||
n2 = 14
|
||||
n3 = 15
|
||||
n = np.array([14, 14, 15])
|
||||
|
||||
X, Y, Z = ndgrid(linspace(0, 1, n1), linspace(0, 1, n2), linspace(0, 1, n3))
|
||||
sigma = 1e-2*ones([n1-1, n2-1, n3-1])
|
||||
sigma[:, :, (n3-1)/2:] = 1e-6
|
||||
mu = 4*pi*1e-7*ones([n1-1, n2-1, n3-1])
|
||||
X, Y, Z = ndgrid(*[np.linspace(0, 1, x) for x in n])
|
||||
sigma = 1e-2*np.ones(n-1)
|
||||
sigma[:, :, (n[2]-1)/2:] = 1e-6
|
||||
mu = 4*np.pi*1e-7*np.ones(n-1)
|
||||
w = 10
|
||||
|
||||
CURL = getCurlMatrix(X, Y, Z)
|
||||
@@ -29,8 +22,8 @@ Me = getEdgeInnerProduct(X, Y, Z, sigma)
|
||||
|
||||
A = CURL.T * Mf * CURL + 1j * w * Me
|
||||
|
||||
ne = shape(A)
|
||||
b = matrix(randn(ne[0])).T
|
||||
ne = np.shape(A)
|
||||
b = np.matrix(randn(ne[0])).T
|
||||
# clean b
|
||||
DIVb = GRAD.T*b
|
||||
p = dsolve.spsolve(GRAD.T*GRAD, DIVb, use_umfpack=True).T
|
||||
|
||||
+39
-39
@@ -1,58 +1,58 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from sputils import *
|
||||
from numpy import *
|
||||
from getEdgeTangent import *
|
||||
from inv3X3BlockDiagonal import *
|
||||
|
||||
def volTetra(y,m,I,A,B,C,D):
|
||||
|
||||
a11 = array(y[A,0]-y[B,0]); a12 = array(y[A,0]-y[C,0]); a13 = array(y[A,0]-y[D,0])
|
||||
a21 = array(y[A,1]-y[B,1]); a22 = array(y[A,1]-y[C,1]); a23 = array(y[A,1]-y[D,1])
|
||||
a31 = array(y[A,2]-y[B,2]); a32 = array(y[A,2]-y[C,2]); a33 = array(y[A,2]-y[D,2])
|
||||
def volTetra(y, m, I, A, B, C, D):
|
||||
|
||||
a11 = array(y[A, 0]-y[B, 0]); a12 = array(y[A, 0]-y[C, 0]); a13 = array(y[A, 0]-y[D, 0])
|
||||
a21 = array(y[A, 1]-y[B, 1]); a22 = array(y[A, 1]-y[C, 1]); a23 = array(y[A, 1]-y[D, 1])
|
||||
a31 = array(y[A, 2]-y[B, 2]); a32 = array(y[A, 2]-y[C, 2]); a33 = array(y[A, 2]-y[D, 2])
|
||||
|
||||
return abs(a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a31*a22*a13 - a32*a23*a11 - a33*a21*a12)
|
||||
|
||||
|
||||
def getCellVolume(X,Y,Z):
|
||||
def getCellVolume(X, Y, Z):
|
||||
|
||||
m = array(shape(X))-1
|
||||
y = hstack3(mkvc(X),mkvc(Y),mkvc(Z))
|
||||
|
||||
i = int64(linspace(0,m[0]-1,m[0]))
|
||||
j = int64(linspace(0,m[1]-1,m[1]))
|
||||
k = int64(linspace(0,m[2]-1,m[2]))
|
||||
|
||||
ii,jj,kk = ndgrid(i,j,k)
|
||||
ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk)
|
||||
|
||||
|
||||
I = int64(sub2ind(m,hstack3(ii,jj,kk)))
|
||||
A = int64(sub2ind(m+1,hstack3(ii,jj,kk)))
|
||||
B = int64(sub2ind(m+1,hstack3(ii,jj+1,kk)))
|
||||
C = int64(sub2ind(m+1,hstack3(ii+1,jj+1,kk)))
|
||||
D = int64(sub2ind(m+1,hstack3(ii+1,jj,kk)))
|
||||
E = int64(sub2ind(m+1,hstack3(ii,jj,kk+1)))
|
||||
F = int64(sub2ind(m+1,hstack3(ii,jj+1,kk+1)))
|
||||
G = int64(sub2ind(m+1,hstack3(ii+1,jj+1,kk+1)))
|
||||
H = int64(sub2ind(m+1,hstack3(ii+1,jj,kk+1)))
|
||||
|
||||
v1 = volTetra(y,m,I,A,B,D,E)
|
||||
v2 = volTetra(y,m,I,B,E,F,G)
|
||||
v3 = volTetra(y,m,I,B,D,E,G)
|
||||
v4 = volTetra(y,m,I,B,C,D,G)
|
||||
v5 = volTetra(y,m,I,D,E,G,H)
|
||||
|
||||
v = 1.0/6.0 * ( v1 + v2 + v3 + v4 + v5 )
|
||||
y = hstack3(mkvc(X), mkvc(Y), mkvc(Z))
|
||||
|
||||
i = int64(linspace(0, m[0]-1, m[0]))
|
||||
j = int64(linspace(0, m[1]-1, m[1]))
|
||||
k = int64(linspace(0, m[2]-1, m[2]))
|
||||
|
||||
ii, jj, kk = ndgrid(i, j, k)
|
||||
ii = mkvc(ii)
|
||||
jj = mkvc(jj)
|
||||
kk = mkvc(kk)
|
||||
|
||||
I = int64(sub2ind(m, hstack3(ii, jj, kk)))
|
||||
A = int64(sub2ind(m+1, hstack3(ii, jj, kk)))
|
||||
B = int64(sub2ind(m+1, hstack3(ii, jj+1, kk)))
|
||||
C = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk)))
|
||||
D = int64(sub2ind(m+1, hstack3(ii+1, jj, kk)))
|
||||
E = int64(sub2ind(m+1, hstack3(ii, jj, kk+1)))
|
||||
F = int64(sub2ind(m+1, hstack3(ii, jj+1, kk+1)))
|
||||
G = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk+1)))
|
||||
H = int64(sub2ind(m+1, hstack3(ii+1, jj, kk+1)))
|
||||
|
||||
v1 = volTetra(y, m, I, A, B, D, E)
|
||||
v2 = volTetra(y, m, I, B, E, F, G)
|
||||
v3 = volTetra(y, m, I, B, D, E, G)
|
||||
v4 = volTetra(y, m, I, B, C, D, G)
|
||||
v5 = volTetra(y, m, I, D, E, G, H)
|
||||
|
||||
v = 1.0/6.0 * (v1 + v2 + v3 + v4 + v5)
|
||||
return v.flatten()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3))
|
||||
Z[2,2,2] = 2.5; Z[0,0,0] = -0.5
|
||||
X[2,2,2] = 2.5; X[0,0,0] = -0.5
|
||||
X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3))
|
||||
Z[2, 2, 2] = 2.5; Z[0, 0, 0] = -0.5
|
||||
X[2, 2, 2] = 2.5; X[0, 0, 0] = -0.5
|
||||
|
||||
v = getCellVolume(X,Y,Z)
|
||||
v = getCellVolume(X, Y, Z)
|
||||
print v
|
||||
|
||||
+70
-67
@@ -1,5 +1,3 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from numpy import *
|
||||
@@ -8,96 +6,101 @@ from getCellVolume import getCellVolume
|
||||
from getFaceNormals import getFaceNormals
|
||||
|
||||
|
||||
#============= Face DIV ===========================
|
||||
def getDivMatrix(X,Y,Z):
|
||||
def getDivMatrix(X, Y, Z):
|
||||
"""Face DIV"""
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]; n2 = n[1]; n3 = n[2]
|
||||
|
||||
n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3 = getFaceNormals(X,Y,Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1),mkvc(area2))),mkvc(area3)))
|
||||
S = sdiag(area)
|
||||
V = getCellVolume(X,Y,Z)
|
||||
|
||||
n1 = n[0]
|
||||
n2 = n[1]
|
||||
n3 = n[2]
|
||||
|
||||
n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
S = sdiag(area)
|
||||
V = getCellVolume(X, Y, Z)
|
||||
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
D1 = kron3(speye(n3),speye(n2),d1)
|
||||
D2 = kron3(speye(n3),d2,speye(n1))
|
||||
D3 = kron3(d3,speye(n2),speye(n1))
|
||||
|
||||
D1 = kron3(speye(n3), speye(n2), d1)
|
||||
D2 = kron3(speye(n3), d2, speye(n1))
|
||||
D3 = kron3(d3, speye(n2), speye(n1))
|
||||
|
||||
# divergence on faces
|
||||
D = appendRight3(D1, D2, D3)
|
||||
|
||||
|
||||
return sdiag(1/V)*D*S
|
||||
|
||||
#============= Edge CURL ===========================
|
||||
|
||||
def getCurlMatrix(X,Y,Z):
|
||||
|
||||
|
||||
|
||||
def getCurlMatrix(X, Y, Z):
|
||||
"""Edge CURL """
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]; n2 = n[1]; n3 = n[2]
|
||||
|
||||
d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3)
|
||||
# derivatives on x-edge variables
|
||||
D32 = kron3(d3,speye(n2),speye(n1+1))
|
||||
D23 = kron3(speye(n3),d2,speye(n1+1))
|
||||
D31 = kron3(d3,speye(n2+1),speye(n1))
|
||||
D13 = kron3(speye(n3),speye(n2+1),d1)
|
||||
D21 = kron3(speye(n3+1),d2,speye(n1))
|
||||
D12 = kron3(speye(n3+1),speye(n2),d1)
|
||||
D32 = kron3(d3, speye(n2), speye(n1+1))
|
||||
D23 = kron3(speye(n3), d2, speye(n1+1))
|
||||
D31 = kron3(d3, speye(n2+1), speye(n1))
|
||||
D13 = kron3(speye(n3), speye(n2+1), d1)
|
||||
D21 = kron3(speye(n3+1), d2, speye(n1))
|
||||
D12 = kron3(speye(n3+1), speye(n2), d1)
|
||||
|
||||
O1 = spzeros(shape(D32)[0],shape(D31)[1])
|
||||
O2 = spzeros(shape(D31)[0],shape(D32)[1])
|
||||
O3 = spzeros(shape(D21)[0],shape(D13)[1])
|
||||
O1 = spzeros(shape(D32)[0], shape(D31)[1])
|
||||
O2 = spzeros(shape(D31)[0], shape(D32)[1])
|
||||
O3 = spzeros(shape(D21)[0], shape(D13)[1])
|
||||
|
||||
CURL = appendBottom3(
|
||||
appendRight3(O1, -D32, D23),
|
||||
appendRight3(D31, O2, -D13),
|
||||
appendRight3(-D21, D12, O3))
|
||||
|
||||
# scale for non-uniform mesh
|
||||
e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,norme1,norme2,norme3 = getEdgeTangent(X,Y,Z)
|
||||
n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3 = getFaceNormals(X,Y,Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1),mkvc(area2))),mkvc(area3)))
|
||||
S = sdiag(1/area)
|
||||
lngth = hstack((hstack((mkvc(norme1),mkvc(norme2))),mkvc(norme3)))
|
||||
L = sdiag(lngth)
|
||||
|
||||
return S*(CURL*L)
|
||||
|
||||
#============= Nodal Gradients ===========================
|
||||
def getNodalGradient(X,Y,Z):
|
||||
|
||||
# scale for non-uniform mesh
|
||||
e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z)
|
||||
n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
S = sdiag(1/area)
|
||||
lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3)))
|
||||
L = sdiag(lngth)
|
||||
|
||||
return S*(CURL*L)
|
||||
|
||||
|
||||
def getNodalGradient(X, Y, Z):
|
||||
"""Nodal Gradients"""
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]; n2 = n[1]; n3 = n[2]
|
||||
|
||||
D1 = kron3(speye(n3+1),speye(n2+1),ddx(n1))
|
||||
D2 = kron3(speye(n3+1),ddx(n2),speye(n1+1))
|
||||
D3 = kron3(ddx(n3),speye(n2+1),speye(n1+1))
|
||||
|
||||
D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1))
|
||||
D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1))
|
||||
D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1))
|
||||
|
||||
# topological gradient
|
||||
GRAD = appendBottom3(D1,D2,D3)
|
||||
|
||||
GRAD = appendBottom3(D1, D2, D3)
|
||||
|
||||
# scale for non-uniform mesh
|
||||
e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,norme1,norme2,norme3 = getEdgeTangent(X,Y,Z)
|
||||
lngth = hstack((hstack((mkvc(norme1),mkvc(norme2))),mkvc(norme3)))
|
||||
L = sdiag(1/lngth)
|
||||
e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z)
|
||||
lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3)))
|
||||
L = sdiag(1/lngth)
|
||||
|
||||
return L*GRAD
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3))
|
||||
Z[2,2,2] = 2.5; Z[0,0,0] = -0.5
|
||||
X[2,2,2] = 2.5; X[0,0,0] = -0.5
|
||||
sig = ones([2,2,2])
|
||||
C = getCurlMatrix(X,Y,Z)
|
||||
return L*GRAD
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3))
|
||||
Z[2, 2, 2] = 2.5
|
||||
Z[0, 0, 0] = -0.5
|
||||
X[2, 2, 2] = 2.5
|
||||
X[0, 0, 0] = -0.5
|
||||
sig = ones([2, 2, 2])
|
||||
C = getCurlMatrix(X, Y, Z)
|
||||
|
||||
G = getNodalGradient(X, Y, Z)
|
||||
|
||||
G = getNodalGradient(X,Y,Z)
|
||||
|
||||
tt = C*G
|
||||
print(tt)
|
||||
print(tt)
|
||||
|
||||
+85
-86
@@ -8,203 +8,202 @@ from getEdgeTangent import *
|
||||
from inv3X3BlockDiagonal import *
|
||||
from getCellVolume import getCellVolume
|
||||
|
||||
# [A] = getEdgeInnerProduct(X,Y,Z,sigma)
|
||||
#
|
||||
|
||||
# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
# / /
|
||||
# / / |
|
||||
# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
# / / |
|
||||
# / / |
|
||||
# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
# | | |
|
||||
# | | node(i+1,j+1,k+1)
|
||||
# | | /
|
||||
# edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k)
|
||||
# | | /
|
||||
# | | /
|
||||
# | |/
|
||||
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
# no | node | e1 | e2 | e3
|
||||
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
|
||||
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
|
||||
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
|
||||
def subarray(T, i1, i2, i3):
|
||||
return take(take(take(T, i1, 0), i2, 1), i3, 2)
|
||||
|
||||
|
||||
def subarray(T,i1,i2,i3):
|
||||
return take(take(take(T,i1,0),i2,1),i3,2)
|
||||
|
||||
def getEdgeInnerProduct(X, Y, Z, sigma):
|
||||
"""A = getEdgeInnerProduct(X, Y, Z, sigma)
|
||||
|
||||
def getEdgeInnerProduct(X,Y,Z,sigma):
|
||||
|
||||
m = array(shape(X))-1
|
||||
node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
/ /
|
||||
/ / |
|
||||
edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
/ / |
|
||||
/ / |
|
||||
node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
| | |
|
||||
| | node(i+1,j+1,k+1)
|
||||
| | /
|
||||
edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k)
|
||||
| | /
|
||||
| | /
|
||||
| |/
|
||||
node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
no | node | e1 | e2 | e3
|
||||
000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
|
||||
011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
|
||||
111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
|
||||
"""
|
||||
m = array(shape(X))-1
|
||||
nc = prod(m)
|
||||
|
||||
|
||||
me1 = m + array([0, 1, 1]); ne1 = prod(me1)
|
||||
me2 = m + array([1, 0, 1]); ne2 = prod(me2)
|
||||
me3 = m + array([1, 1, 0]); ne3 = prod(me3)
|
||||
|
||||
|
||||
e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,norme1,norme2,norme3 = getEdgeTangent(X,Y,Z)
|
||||
|
||||
|
||||
i = int64(linspace(0,m[0]-1,m[0]))
|
||||
j = int64(linspace(0,m[1]-1,m[1]))
|
||||
k = int64(linspace(0,m[2]-1,m[2]))
|
||||
|
||||
ii,jj,kk = ndgrid(i,j,k)
|
||||
|
||||
ii,jj,kk = ndgrid(i,j,k)
|
||||
ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk)
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P000 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT000 = inv3X3BlockDiagonal(subarray(e1x,i,j,k) , subarray(e1y,i,j,k), subarray(e1z,i,j,k),
|
||||
subarray(e2x,i,j,k) , subarray(e2y,i,j,k), subarray(e2z,i,j,k) ,
|
||||
subarray(e3x,i,j,k) , subarray(e3y,i,j,k), subarray(e3z,i,j,k) )
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P100 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT100 = inv3X3BlockDiagonal(subarray(e1x,i,j,k), subarray(e1y,i,j,k), subarray(e1z,i,j,k),
|
||||
subarray(e2x,i+1,j,k), subarray(e2y,i+1,j,k), subarray(e2z,i+1,j,k),
|
||||
subarray(e3x,i+1,j,k) , subarray(e3y,i+1,j,k), subarray(e3z,i+1,j,k))
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P010 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT010 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) , subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) ,
|
||||
subarray(e2x,i,j,k) , subarray(e2y,i,j,k) , subarray(e2z,i,j,k) ,
|
||||
subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) ,subarray( e3z,i,j+1,k) )
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P110 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT110 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) ,subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) ,
|
||||
subarray(e2x,i+1,j,k) ,subarray(e2y,i+1,j,k) , subarray(e2z,i+1,j,k),
|
||||
subarray(e3x,i+1,j+1,k) ,subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) )
|
||||
|
||||
|
||||
######
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P001 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT001 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1) ,subarray(e1y,i,j,k+1) , subarray(e1z,i,j,k+1) ,
|
||||
subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) ,
|
||||
subarray(e3x,i,j,k) , subarray(e3y,i,j,k) , subarray(e3z,i,j,k) )
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P101 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT101 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1), subarray(e1y,i,j,k+1), subarray(e1z,i,j,k+1) ,
|
||||
subarray(e2x,i+1,j,k+1), subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) ,
|
||||
subarray(e3x,i+1,j,k), subarray(e3y,i+1,j,k) , subarray(e3z,i+1,j,k) )
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P011 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT011 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) ,
|
||||
subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) ,
|
||||
subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) , subarray(e3z,i,j+1,k) )
|
||||
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
|
||||
P111 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
|
||||
invT111 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) ,
|
||||
subarray(e2x,i+1,j,k+1) , subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) ,
|
||||
subarray(e3x,i+1,j+1,k) , subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) )
|
||||
|
||||
|
||||
# Cell volume
|
||||
v = mkvc(getCellVolume(X,Y,Z)) #mkvc(getVolume(X,Y,Z))
|
||||
vsig = v*mkvc(sigma)
|
||||
v3 = vstack((vstack((vsig,vsig)),vsig))
|
||||
v3 = v3.flatten()
|
||||
|
||||
|
||||
V = sdiag(v3)
|
||||
|
||||
A = P000.T*invT000.T*V*invT000*P000 + P001.T*invT001.T*V*invT001*P001 + P010.T*invT010.T*V*invT010*P010 + P011.T*invT011.T*V*invT011*P011 + P100.T*invT100.T*V*invT100*P100 + P101.T*invT101.T*V*invT101*P101 + P110.T*invT110.T*V*invT110*P110 + P111.T*invT111.T*V*invT111*P111
|
||||
|
||||
|
||||
A = P000.T*invT000.T*V*invT000*P000 + P001.T*invT001.T*V*invT001*P001 + P010.T*invT010.T*V*invT010*P010 + P011.T*invT011.T*V*invT011*P011 + P100.T*invT100.T*V*invT100*P100 + P101.T*invT101.T*V*invT101*P101 + P110.T*invT110.T*V*invT110*P110 + P111.T*invT111.T*V*invT111*P111
|
||||
|
||||
A = 0.125*A
|
||||
|
||||
|
||||
return A
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3))
|
||||
|
||||
@@ -1,12 +1,11 @@
|
||||
from scipy.sparse import linalg
|
||||
from utils import *
|
||||
from sputils import *
|
||||
|
||||
|
||||
def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,a33):
|
||||
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
|
||||
|
||||
a11 = mkvc(a11)
|
||||
a12 = mkvc(a12)
|
||||
a11 = mkvc(a11)
|
||||
a12 = mkvc(a12)
|
||||
a13 = mkvc(a13)
|
||||
a21 = mkvc(a21)
|
||||
a22 = mkvc(a22)
|
||||
@@ -17,16 +16,16 @@ def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,a33):
|
||||
|
||||
detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
|
||||
|
||||
b11 = (a22*a33 - a23*a32)/detA
|
||||
b11 = +(a22*a33 - a23*a32)/detA
|
||||
b12 = -(a12*a33 - a13*a32)/detA
|
||||
b13 = (a12*a23 - a13*a22)/detA
|
||||
b13 = +(a12*a23 - a13*a22)/detA
|
||||
|
||||
b21 = (a31*a23 - a21*a33)/detA
|
||||
b21 = +(a31*a23 - a21*a33)/detA
|
||||
b22 = -(a31*a13 - a11*a33)/detA
|
||||
b23 = (a21*a13 - a11*a23)/detA
|
||||
b23 = +(a21*a13 - a11*a23)/detA
|
||||
|
||||
b31 = -(a31*a22 - a21*a32)/detA
|
||||
b32 = (a31*a12 - a11*a32)/detA
|
||||
b32 = +(a31*a12 - a11*a32)/detA
|
||||
b33 = -(a21*a12 - a11*a22)/detA
|
||||
|
||||
B = appendBottom3(
|
||||
@@ -34,4 +33,4 @@ def inv3X3BlockDiagonal(a11,a12,a13,a21,a22,a23,a31,a32,a33):
|
||||
appendRight3(sdiag(b21), sdiag(b22), sdiag(b23)),
|
||||
appendRight3(sdiag(b31), sdiag(b32), sdiag(b33)))
|
||||
|
||||
return B
|
||||
return B
|
||||
|
||||
+80
-87
@@ -1,101 +1,94 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from numpy import *
|
||||
|
||||
#----- Cell Centers from Nodal locations -----
|
||||
def getCellCenterFromNodal(X,Y,Z):
|
||||
|
||||
XC = 1.0/8.0 * (X[0:-1,0:-1,0:-1] + X[1:,0:-1,0:-1] + X[0:-1,1:,0:-1] + X[1:,1:,0:-1] +
|
||||
X[0:-1,0:-1,1:] + X[1:,0:-1,1:] + X[0:-1,1:,1:] + X[1:,1:,1:])
|
||||
def getCellCenterFromNodal(X, Y, Z):
|
||||
"""Cell Centers from Nodal locations"""
|
||||
XC = 1.0/8.0 * (X[:-1, :-1, :-1] + X[1:, :-1, :-1] + X[:-1, 1:, :-1] + X[1:, 1:, :-1] +
|
||||
X[:-1, :-1, 1:] + X[1:, :-1, 1:] + X[:-1, 1:, 1:] + X[1:, 1:, 1:])
|
||||
|
||||
YC = 1.0/8.0 * (Y[0:-1,0:-1,0:-1] + Y[1:,0:-1,0:-1] + Y[0:-1,1:,0:-1] + Y[1:,1:,0:-1] +
|
||||
Y[0:-1,0:-1,1:] + Y[1:,0:-1,1:] + Y[0:-1,1:,1:] + Y[1:,1:,1:])
|
||||
YC = 1.0/8.0 * (Y[:-1, :-1, :-1] + Y[1:, :-1, :-1] + Y[:-1, 1:, :-1] + Y[1:, 1:, :-1] +
|
||||
Y[:-1, :-1, 1:] + Y[1:, :-1, 1:] + Y[:-1, 1:, 1:] + Y[1:, 1:, 1:])
|
||||
|
||||
ZC = 1.0/8.0 * (Z[0:-1,0:-1,0:-1] + Z[1:,0:-1,0:-1] + Z[0:-1,1:,0:-1] + Z[1:,1:,0:-1] +
|
||||
Z[0:-1,0:-1,1:] + Z[1:,0:-1,1:] + Z[0:-1,1:,1:] + Z[1:,1:,1:])
|
||||
|
||||
return (XC,YC,ZC)
|
||||
ZC = 1.0/8.0 * (Z[:-1, :-1, :-1] + Z[1:, :-1, :-1] + Z[:-1, 1:, :-1] + Z[1:, 1:, :-1] +
|
||||
Z[:-1, :-1, 1:] + Z[1:, :-1, 1:] + Z[:-1, 1:, 1:] + Z[1:, 1:, 1:])
|
||||
|
||||
#----- Edges from Nodal locations -----
|
||||
|
||||
def getEdgesFromNodal(X,Y,Z):
|
||||
#
|
||||
# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
# / /
|
||||
# / / |
|
||||
# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
# / / |
|
||||
# / / |
|
||||
# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
# | | |
|
||||
# | | node(i+1,j+1,k+1)
|
||||
# | | /
|
||||
# edge1(i,j,k) face3(i,j,k) edge1(i,j+1.k)
|
||||
# | | /
|
||||
# | | /
|
||||
# | |/
|
||||
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
XE1 = (X[1:,:,:]+X[0:-1,:,:])/2.0
|
||||
YE1 = (Y[1:,:,:]+Y[0:-1,:,:])/2.0
|
||||
ZE1 = (Z[1:,:,:]+Z[0:-1,:,:])/2.0
|
||||
|
||||
XE2 = (X[:,1:,:]+X[:,0:-1,:])/2.0
|
||||
YE2 = (Y[:,1:,:]+Y[:,0:-1,:])/2.0
|
||||
ZE2 = (Z[:,1:,:]+Z[:,0:-1,:])/2.0
|
||||
|
||||
XE3 = (X[:,:,1:]+X[:,:,0:-1])/2.0
|
||||
YE3 = (Y[:,:,1:]+Y[:,:,0:-1])/2.0
|
||||
ZE3 = (Z[:,:,1:]+Z[:,:,0:-1])/2.0
|
||||
|
||||
return (XE1,YE1,ZE1,XE2,YE2,ZE2,XE3,YE3,ZE3)
|
||||
|
||||
#-- Get faces from nodal --
|
||||
|
||||
def getFacesFromNodal(X,Y,Z):
|
||||
|
||||
XF1 = 1.0/4.0*(X[:,0:-1,0:-1]+X[:,1:,0:-1]+X[:,0:-1,1:]+X[:,1:,1:])
|
||||
YF1 = 1.0/4.0*(Y[:,0:-1,0:-1]+Y[:,1:,0:-1]+Y[:,0:-1,1:]+Y[:,1:,1:])
|
||||
ZF1 = 1.0/4.0*(Z[:,0:-1,0:-1]+Z[:,1:,0:-1]+Z[:,0:-1,1:]+Z[:,1:,1:])
|
||||
|
||||
XF2 = 1.0/4.0*(X[0:-1,:,0:-1]+X[1:,:,0:-1]+X[0:-1,:,1:]+X[1:,:,1:])
|
||||
YF2 = 1.0/4.0*(Y[0:-1,:,0:-1]+Y[1:,:,0:-1]+Y[0:-1,:,1:]+Y[1:,:,1:])
|
||||
ZF2 = 1.0/4.0*(Z[0:-1,:,0:-1]+Z[1:,:,0:-1]+Z[0:-1,:,1:]+Z[1:,:,1:])
|
||||
|
||||
XF3 = 1.0/4.0*(X[0:-1,0:-1,:]+X[1:,0:-1,:]+X[0:-1,1:,:]+X[1:,1:,:])
|
||||
YF3 = 1.0/4.0*(Y[0:-1,0:-1,:]+Y[1:,0:-1,:]+Y[0:-1,1:,:]+Y[1:,1:,:])
|
||||
ZF3 = 1.0/4.0*(Z[0:-1,0:-1,:]+Z[1:,0:-1,:]+Z[0:-1,1:,:]+Z[1:,1:,:])
|
||||
|
||||
return (XF1,YF1,ZF1,XF2,YF2,ZF2,XF3,YF3,ZF3)
|
||||
|
||||
#-- Project Edge vector field
|
||||
|
||||
def projectEdgeVectorField(EV1,EV2,EV3,X,Y,Z):
|
||||
|
||||
t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,nrm1,nrm2,nrm3 = getEdgeTangent(X,Y,Z)
|
||||
return (XC, YC, ZC)
|
||||
|
||||
|
||||
E1 = EV1[:,0]*mkvc(t1x) + EV1[:,1]*mkvc(t1y) + EV1[:,2]*mkvc(t1z)
|
||||
E2 = EV2[:,0]*mkvc(t2x) + EV2[:,1]*mkvc(t2y) + EV2[:,2]*mkvc(t2z)
|
||||
E3 = EV3[:,0]*mkvc(t3x) + EV3[:,1]*mkvc(t3y) + EV3[:,2]*mkvc(t3z)
|
||||
|
||||
return hstack((hstack((mkvc(E1),mkvc(E2))),mkvc(E3)))
|
||||
|
||||
|
||||
#-- Prolect Face vector field
|
||||
def getEdgesFromNodal(X, Y, Z):
|
||||
"""Edges from Nodal locations
|
||||
|
||||
def projectFaceVectorField(FV1,FV2,FV3,X,Y,Z):
|
||||
node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
/ /
|
||||
/ / |
|
||||
edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
/ / |
|
||||
/ / |
|
||||
node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
| | |
|
||||
| | node(i+1,j+1,k+1)
|
||||
| | /
|
||||
edge1(i,j,k) face3(i,j,k) edge1(i,j+1.k)
|
||||
| | /
|
||||
| | /
|
||||
| |/
|
||||
node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
"""
|
||||
|
||||
n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,ar1,ar2,ar3 = getFaceNormals(X,Y,Z)
|
||||
XE1 = (X[1:, :, :]+X[:-1, :, :])/2.0
|
||||
YE1 = (Y[1:, :, :]+Y[:-1, :, :])/2.0
|
||||
ZE1 = (Z[1:, :, :]+Z[:-1, :, :])/2.0
|
||||
|
||||
F1 = FV1[:,0]*mkvc(n1x) + FV1[:,1]*mkvc(n1y) + FV1[:,2]*mkvc(n1z)
|
||||
F2 = FV2[:,0]*mkvc(n2x) + FV2[:,1]*mkvc(n2y) + FV2[:,2]*mkvc(n2z)
|
||||
F3 = FV3[:,0]*mkvc(n3x) + FV3[:,1]*mkvc(n3y) + FV3[:,2]*mkvc(n3z)
|
||||
|
||||
return hstack((hstack((mkvc(F1),mkvc(F2))),mkvc(F3)))
|
||||
|
||||
XE2 = (X[:, 1:, :]+X[:, :-1, :])/2.0
|
||||
YE2 = (Y[:, 1:, :]+Y[:, :-1, :])/2.0
|
||||
ZE2 = (Z[:, 1:, :]+Z[:, :-1, :])/2.0
|
||||
|
||||
|
||||
|
||||
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,29 +0,0 @@
|
||||
from numpy import *
|
||||
|
||||
def ndgrid(x,y,z):
|
||||
|
||||
n1 = size(x)
|
||||
n2 = size(y)
|
||||
n3 = size(z)
|
||||
X = zeros([n1,n2,n3])
|
||||
Y = zeros([n1,n2,n3])
|
||||
Z = zeros([n1,n2,n3])
|
||||
for i in range(0, n2):
|
||||
for j in range(0,n3):
|
||||
X[:,i,j] = x
|
||||
|
||||
for i in range(0, n1):
|
||||
for j in range(0,n3):
|
||||
Y[i,:,j] = y
|
||||
|
||||
for i in range(0, n1):
|
||||
for j in range(0,n2):
|
||||
Z[i,j,:] = z
|
||||
|
||||
|
||||
return (X,Y,Z)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X = ndgrid([1,2,3],[2,4,5,6],[4,6,7,8])
|
||||
+59
-49
@@ -1,67 +1,77 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from numpy import *
|
||||
|
||||
|
||||
#======== Define 1D derivatives =============
|
||||
def ddx(n):
|
||||
return sparse.spdiags(-ones(n),0,n,n+1) + sparse.spdiags(ones(n+1),1,n,n+1)
|
||||
"""Define 1D derivatives"""
|
||||
# ddx = lambda n: sparse.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format='csr')
|
||||
return sparse.spdiags(-ones(n), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)
|
||||
|
||||
|
||||
#======== Define 1D average =============
|
||||
def av(n):
|
||||
return 0.5*(sparse.spdiags(ones(n+1),0,n,n+1) + sparse.spdiags(ones(n+1),1,n,n+1))
|
||||
"""Define 1D average"""
|
||||
return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1))
|
||||
|
||||
|
||||
#======== Diagonal matrix =============
|
||||
def sdiag(h):
|
||||
return sparse.spdiags(h,0,size(h),size(h))
|
||||
"""Diagonal matrix"""
|
||||
return sparse.spdiags(h, 0, size(h), size(h))
|
||||
|
||||
|
||||
#======== sparse identity =============
|
||||
def speye(n):
|
||||
return sparse.spdiags(ones(n),0,n,n)
|
||||
"""sparse identity"""
|
||||
return sparse.spdiags(ones(n), 0, n, n)
|
||||
|
||||
#======== two kron prods =============
|
||||
def kron3(A,B,C):
|
||||
return sparse.kron(sparse.kron(A,B),C)
|
||||
|
||||
#======== append on bottom =============
|
||||
def appendBottom(A,B):
|
||||
C = sparse.vstack((A,B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
#======== append on bottom =============
|
||||
def appendBottom3(A,B,C):
|
||||
C = appendBottom(appendBottom(A,B),C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
#======== append on right =============
|
||||
def appendRight(A,B):
|
||||
C = sparse.hstack((A,B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
def kron3(A, B, C):
|
||||
"""two kron prods"""
|
||||
return sparse.kron(sparse.kron(A, B), C)
|
||||
|
||||
#======== append on right =============
|
||||
def appendRight3(A,B,C):
|
||||
C = appendRight(appendRight(A,B),C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
#======== blockdigonal =============
|
||||
def blkDiag(A,B):
|
||||
O12 = sparse.coo_matrix((shape(A)[0],shape(B)[1]))
|
||||
O21 = sparse.coo_matrix((shape(B)[0],shape(A)[1]))
|
||||
C = sparse.vstack((sparse.hstack((A,O12)),sparse.hstack((O21,B))))
|
||||
def appendBottom(A, B):
|
||||
"""append on bottom"""
|
||||
C = sparse.vstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
#======== blockdigonal 3 =============
|
||||
def blkDiag3(A,B,C):
|
||||
ABC = blkDiag(blkDiag(A,B),C)
|
||||
return C
|
||||
|
||||
|
||||
def appendBottom3(A, B, C):
|
||||
"""append on bottom"""
|
||||
C = appendBottom(appendBottom(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def appendRight(A, B):
|
||||
"""append on right"""
|
||||
C = sparse.hstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def appendRight3(A, B, C):
|
||||
"""append on right"""
|
||||
C = appendRight(appendRight(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def blkDiag(A, B):
|
||||
"""blockdigonal"""
|
||||
O12 = sparse.coo_matrix((shape(A)[0], shape(B)[1]))
|
||||
O21 = sparse.coo_matrix((shape(B)[0], shape(A)[1]))
|
||||
C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B))))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def blkDiag3(A, B, C):
|
||||
"""blockdigonal 3"""
|
||||
ABC = blkDiag(blkDiag(A, B), C)
|
||||
ABC = ABC.tocsr()
|
||||
return ABC
|
||||
|
||||
#======== spzeros =============
|
||||
def spzeros(n1,n2):
|
||||
return sparse.coo_matrix((n1,n2))
|
||||
|
||||
|
||||
|
||||
def spzeros(n1, n2):
|
||||
"""spzeros"""
|
||||
return sparse.coo_matrix((n1, n2))
|
||||
|
||||
+83
-80
@@ -1,124 +1,127 @@
|
||||
from numpy import *
|
||||
|
||||
def diff(A,d):
|
||||
|
||||
end = -1
|
||||
if(d==1):
|
||||
return A[1:,0:,0:] - A[0:end,0:,0:]
|
||||
elif(d==2):
|
||||
return A[0:,1:,0:] - A[0:,0:end,0:]
|
||||
def diff(A, d):
|
||||
if(d == 1):
|
||||
return A[1:, :, :] - A[:-1, :, :]
|
||||
elif(d == 2):
|
||||
return A[:, 1:, :] - A[:, :-1, :]
|
||||
else:
|
||||
return A[0:,0:,1:] - A[0:,0:,0:end]
|
||||
return A[:, :, 1:] - A[:, :, :-1]
|
||||
#else:
|
||||
# print('d must be 1,2 or 3')
|
||||
|
||||
|
||||
def diffp(A, d1, d2):
|
||||
end = -1
|
||||
if(d1 == 1 and d2 == 2 ):
|
||||
return A[1:,1:, 0:] - A[0:end,0:end,0:]
|
||||
if(d1 == 1 and d2 == 2):
|
||||
return A[1:, 1:, :] - A[:-1, :-1, :]
|
||||
elif(d1 == 1 and d2 == 3):
|
||||
return A[1:,0:,1:] - A[0:end,0:,0:end]
|
||||
return A[1:, :, 1:] - A[:-1, :, :-1]
|
||||
else:
|
||||
return A[0:,1:,1:] - A[0:,0:end,0:end]
|
||||
return A[:, 1:, 1:] - A[:, :-1, :-1]
|
||||
|
||||
|
||||
def diffm(A, d1, d2):
|
||||
end = -1
|
||||
if(d1 == 3 and d2 == 2 ):
|
||||
return A[:,0:end,1:] - A[:,1:,0:end]
|
||||
if(d1 == 3 and d2 == 2):
|
||||
return A[:, :-1, 1:] - A[:, 1:, :-1]
|
||||
elif(d1 == 1 and d2 == 3):
|
||||
return A[1:, :, 0:end] - A[0:end,:,1:]
|
||||
return A[1:, :, :-1] - A[:-1, :, 1:]
|
||||
elif(d1 == 2 and d2 == 1):
|
||||
return A[0:end, 1:, :] - A[1:, 0:end, :]
|
||||
return A[:-1, 1:, :] - A[1:, :-1, :]
|
||||
else:
|
||||
print('d must be 1, 2 or 3')
|
||||
|
||||
|
||||
def ave(A, d):
|
||||
if(d == 1):
|
||||
return 0.5*(A[1:, :, :] + A[:-1, :, :])
|
||||
elif(d == 2):
|
||||
return 0.5*(A[:, 1:, :] + A[:, :-1, :])
|
||||
elif(d == 3):
|
||||
return 0.5*(A[:, :, 1:] + A[:, :, :-1])
|
||||
else:
|
||||
print('d must be 1,2 or 3')
|
||||
|
||||
def ave(A,d):
|
||||
|
||||
end = 0
|
||||
if(d==1):
|
||||
return 0.5*(A[1:,:,:] + A[0:end-1,0:,:])
|
||||
elif(d==2):
|
||||
return 0.5*(A[:,1:,:] + A[0:,0:end-1,:])
|
||||
elif(d==3):
|
||||
return 0.5*(A[:,:,1:] + A[0:,0:,0:end-1])
|
||||
else:
|
||||
print('d must be 1,2 or 3')
|
||||
|
||||
def reshapeF(sp,d):
|
||||
return reshape(sp,d,'F')
|
||||
|
||||
def mkvc(A):
|
||||
return reshape(A,[size(A),1],'F').flatten()
|
||||
|
||||
|
||||
def ndgrid(x,y,z):
|
||||
|
||||
def reshapeF(sp, d):
|
||||
return reshape(sp, d, 'F')
|
||||
|
||||
|
||||
def mkvc(A):
|
||||
return reshape(A, [size(A), 1], 'F').flatten()
|
||||
|
||||
|
||||
def ndgrid(x, y, z):
|
||||
|
||||
n1 = size(x)
|
||||
n2 = size(y)
|
||||
n3 = size(z)
|
||||
X = zeros([n1,n2,n3])
|
||||
Y = zeros([n1,n2,n3])
|
||||
Z = zeros([n1,n2,n3])
|
||||
X = zeros([n1, n2, n3])
|
||||
Y = zeros([n1, n2, n3])
|
||||
Z = zeros([n1, n2, n3])
|
||||
for i in range(0, n2):
|
||||
for j in range(0,n3):
|
||||
X[:,i,j] = x
|
||||
|
||||
for j in range(0, n3):
|
||||
X[:, i, j] = x
|
||||
|
||||
for i in range(0, n1):
|
||||
for j in range(0,n3):
|
||||
Y[i,:,j] = y
|
||||
|
||||
for j in range(0, n3):
|
||||
Y[i, :, j] = y
|
||||
|
||||
for i in range(0, n1):
|
||||
for j in range(0,n2):
|
||||
Z[i,j,:] = z
|
||||
|
||||
|
||||
return (X,Y,Z)
|
||||
|
||||
for j in range(0, n2):
|
||||
Z[i, j, :] = z
|
||||
|
||||
return (X, Y, Z)
|
||||
|
||||
|
||||
def ind2sub(shape, ind):
|
||||
# From the given shape, returns the subscrips of the given index
|
||||
revshp = [];
|
||||
revshp.extend(shape);
|
||||
mult = [1];
|
||||
revshp = []
|
||||
revshp.extend(shape)
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]]);
|
||||
mult = array(mult).reshape(len(mult));
|
||||
|
||||
sub = [];
|
||||
|
||||
for i in range(0,len(shape)):
|
||||
sub.extend([math.floor(ind / mult[i])]);
|
||||
ind = ind - (math.floor(ind/mult[i]) * mult[i]);
|
||||
return sub;
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult))
|
||||
|
||||
sub = []
|
||||
|
||||
for i in range(0, len(shape)):
|
||||
sub.extend([math.floor(ind / mult[i])])
|
||||
ind = ind - (math.floor(ind/mult[i]) * mult[i])
|
||||
return sub
|
||||
|
||||
|
||||
def sub2ind(shape, subs):
|
||||
# From the given shape, returns the index of the given subscript
|
||||
revshp = list(shape);
|
||||
mult = [1];
|
||||
revshp = list(shape)
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]]);
|
||||
mult = array(mult).reshape(len(mult),1);
|
||||
|
||||
idx = dot((subs) , (mult));
|
||||
return idx;
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult), 1)
|
||||
|
||||
idx = dot((subs), (mult))
|
||||
return idx
|
||||
|
||||
|
||||
def mkmat(x):
|
||||
return reshape(matrix(x),(size(x),1),'F')
|
||||
|
||||
def hstack3(a,b,c):
|
||||
a = mkvc(a); b = mkvc(b); c = mkvc(c)
|
||||
a = mkmat(a); b = mkmat(b); c = mkmat(c)
|
||||
return hstack((hstack((a,b)),c))
|
||||
|
||||
|
||||
return reshape(matrix(x), (size(x), 1), 'F')
|
||||
|
||||
|
||||
def hstack3(a, b, c):
|
||||
a = mkvc(a)
|
||||
b = mkvc(b)
|
||||
c = mkvc(c)
|
||||
a = mkmat(a)
|
||||
b = mkmat(b)
|
||||
c = mkmat(c)
|
||||
return hstack((hstack((a, b)), c))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = mgrid[0:4, 0:5, 0:6]
|
||||
|
||||
print Z
|
||||
|
||||
t = ave(X, 1)
|
||||
|
||||
print t
|
||||
|
||||
Reference in New Issue
Block a user