mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-27 11:27:23 +08:00
renamed folder 'code' to 'SimPEG'
new folder for ipython notebooks improved 2D plots
This commit is contained in:
Lars Ruthotto
@@ -0,0 +1,205 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
class BaseMesh(object):
|
||||
"""BaseMesh does all the counting you don't want to do.
|
||||
|
||||
x0 origin ndarray (dim, )
|
||||
n number of cells ndarray (dim, )
|
||||
dim dimension of mesh int 1, 2, or 3
|
||||
|
||||
nCx num cells in x dir int
|
||||
nCy num cells in y dir int
|
||||
nCz num cells in z dir int
|
||||
nC total number of cells int
|
||||
|
||||
nNx num nodes in x dir int
|
||||
nNy num nodes in y dir int
|
||||
nNz num nodes in z dir int
|
||||
nN total number of nodes int
|
||||
|
||||
nEx num edges in x dir ndarray [nEx_x, nEx_y, nEx_z]
|
||||
nEy num edges in y dir ndarray [nEy_x, nEy_y, nEy_z]
|
||||
nEz num edges in z dir ndarray [nEz_x, nEz_y, nEz_z]
|
||||
nE total number of edges ndarray (dim, )
|
||||
|
||||
nFx num faces in x dir ndarray [nFx_x, nFx_y, nFx_z]
|
||||
nFy num faces in y dir ndarray [nFy_x, nFy_y, nFy_z]
|
||||
nFz num faces in z dir ndarray [nFz_x, nFz_y, nFz_z]
|
||||
nF total number of faces ndarray (dim, )
|
||||
|
||||
"""
|
||||
def __init__(self, n, x0=None):
|
||||
|
||||
# Check inputs
|
||||
if x0 is None:
|
||||
x0 = np.zeros(len(n))
|
||||
|
||||
if not len(n) == len(x0):
|
||||
raise Exception("Dimension mismatch. x0 != len(n)")
|
||||
|
||||
if len(n) > 3:
|
||||
raise Exception("Dimensions higher than 3 are not supported.")
|
||||
|
||||
# Ensure x0 & n are 1D vectors
|
||||
self._n = np.array(n, dtype=int).ravel()
|
||||
self._x0 = np.array(x0).ravel()
|
||||
self._dim = len(n)
|
||||
|
||||
def x0():
|
||||
doc = "Origin of the mesh"
|
||||
fget = lambda self: self._x0
|
||||
return locals()
|
||||
x0 = property(**x0())
|
||||
|
||||
def n():
|
||||
doc = "Number of Cells in each dimension (array of integers)"
|
||||
fget = lambda self: self._n
|
||||
return locals()
|
||||
n = property(**n())
|
||||
|
||||
def dim():
|
||||
doc = "The dimension of the mesh (1, 2, or 3)."
|
||||
fget = lambda self: self._dim
|
||||
return locals()
|
||||
dim = property(**dim())
|
||||
|
||||
def nCx():
|
||||
doc = "Number oc cells in the x direction"
|
||||
fget = lambda self: self.n[0]
|
||||
return locals()
|
||||
nCx = property(**nCx())
|
||||
|
||||
def nCy():
|
||||
doc = "Number of cells in the y direction"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 1:
|
||||
return self.n[1]
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nCy = property(**nCy())
|
||||
|
||||
def nCz():
|
||||
doc = "Number of cells in the z direction"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 2:
|
||||
return self.n[2]
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nCz = property(**nCz())
|
||||
|
||||
def nC():
|
||||
doc = "Total number of cells"
|
||||
fget = lambda self: np.prod(self.n)
|
||||
return locals()
|
||||
nC = property(**nC())
|
||||
|
||||
def nNx():
|
||||
doc = "Number of nodes in the x-direction"
|
||||
fget = lambda self: self.nCx + 1
|
||||
return locals()
|
||||
nNx = property(**nNx())
|
||||
|
||||
def nNy():
|
||||
doc = "Number of noes in the y-direction"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 1:
|
||||
return self.n[1] + 1
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nNy = property(**nNy())
|
||||
|
||||
def nNz():
|
||||
doc = "Number of nodes in the z-direction"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 2:
|
||||
return self.n[2] + 1
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nNz = property(**nNz())
|
||||
|
||||
def nN():
|
||||
doc = "Total number of nodes"
|
||||
fget = lambda self: np.prod(self.n + 1)
|
||||
return locals()
|
||||
nN = property(**nN())
|
||||
|
||||
def nEx():
|
||||
doc = "Number of x-edges"
|
||||
fget = lambda self: np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
|
||||
return locals()
|
||||
nEx = property(**nEx())
|
||||
|
||||
def nEy():
|
||||
doc = "Number of y-edges"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 1:
|
||||
return np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nEy = property(**nEy())
|
||||
|
||||
def nEz():
|
||||
doc = "Number of z-edges"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 2:
|
||||
return np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nEz = property(**nEz())
|
||||
|
||||
def nE():
|
||||
doc = "Total number of edges"
|
||||
fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
|
||||
return locals()
|
||||
nE = property(**nE())
|
||||
|
||||
def nFx():
|
||||
doc = "Number of x-faces"
|
||||
fget = lambda self: np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
|
||||
return locals()
|
||||
nFx = property(**nFx())
|
||||
|
||||
def nFy():
|
||||
doc = "Number of y-faces"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 1:
|
||||
return np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nFy = property(**nFy())
|
||||
|
||||
def nFz():
|
||||
doc = "Number of z-faces"
|
||||
|
||||
def fget(self):
|
||||
if self.dim > 2:
|
||||
return np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
|
||||
else:
|
||||
return None
|
||||
return locals()
|
||||
nFz = property(**nFz())
|
||||
|
||||
def nF():
|
||||
doc = "Total number of faces in each dimension"
|
||||
fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
|
||||
return locals()
|
||||
nF = property(**nF())
|
||||
|
||||
if __name__ == '__main__':
|
||||
m = BaseMesh([3, 2, 4])
|
||||
print m.n
|
||||
@@ -0,0 +1,141 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from pylab import norm
|
||||
|
||||
def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3,
|
||||
tolG=1e-3, eps=1e-16, xStop=[]):
|
||||
"""
|
||||
GaussNewton Optimization
|
||||
|
||||
Input:
|
||||
------
|
||||
fctn - objective Function (lambda function)
|
||||
x0 - starting guess
|
||||
|
||||
Output:
|
||||
-------
|
||||
xOpt - numerical optimizer
|
||||
"""
|
||||
# initial output
|
||||
print "%s GaussNewton %s" % ('='*22,'='*22)
|
||||
print "iter\tJc\t\tnorm(dJ)\tLS"
|
||||
print "%s" % '-'*57
|
||||
|
||||
# evaluate stopping criteria
|
||||
if xStop==[]:
|
||||
xStop=x0
|
||||
Jstop = fctn(xStop)
|
||||
print "%3d\t%1.2e" % (-1, Jstop[0])
|
||||
|
||||
# initialize
|
||||
xc = x0
|
||||
STOP = np.zeros((5,1),dtype=bool)
|
||||
iterLS=0; iter=0
|
||||
|
||||
Jold = Jstop
|
||||
xOld=xc
|
||||
while 1:
|
||||
# evaluate objective function
|
||||
Jc,dJ,H = fctn(xc)
|
||||
print "%3d\t%1.2e\t%1.2e\t%d" % (iter, Jc[0],norm(dJ),iterLS)
|
||||
|
||||
# check stopping rules
|
||||
STOP[0] = (iter>0) & (abs(Jc[0]-Jold[0]) <= tolJ*(1+abs(Jstop[0])))
|
||||
STOP[1] = (iter>0) & (norm(xc-xOld) <= tolX*(1+norm(x0)))
|
||||
STOP[2] = norm(dJ) <= tolG*(1+abs(Jstop[0]))
|
||||
STOP[3] = norm(dJ) <= 1e3*eps
|
||||
STOP[4] = (iter >= maxIter)
|
||||
if all(STOP[0:3]) | any(STOP[3:]):
|
||||
break
|
||||
|
||||
# get search direction
|
||||
dx = np.linalg.solve(H,-dJ)
|
||||
|
||||
# Armijo linesearch
|
||||
descent = np.dot(dJ.T,dx)
|
||||
LS =0; t = 1; iterLS=1
|
||||
while (iterLS<maxIterLS):
|
||||
xt = xc + t*dx
|
||||
Jt = fctn(xt)
|
||||
LS = Jt[0]<Jc[0]+t*LSreduction*descent
|
||||
if LS:
|
||||
break
|
||||
iterLS = iterLS+1
|
||||
t = .5*t
|
||||
|
||||
# store old values
|
||||
Jold = Jc; xOld = xc
|
||||
# update
|
||||
xc = xt
|
||||
iter = iter +1
|
||||
|
||||
print "%s STOP! %s" % ('-'*25,'-'*25)
|
||||
print "%d : |Jc-Jold| = %1.4e <= tolJ*(1+|Jstop|) = %1.4e" % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0])))
|
||||
print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0)))
|
||||
print "%d : |dJ| = %1.4e <= tolG*(1+|Jstop|) = %1.4e" % (STOP[2],norm(dJ),tolG*(1+abs(Jstop[0])))
|
||||
print "%d : |dJ| = %1.4e <= 1e3*eps = %1.4e" % (STOP[3],norm(dJ),1e3*eps)
|
||||
print "%d : iter = %3d\t <= maxIter\t = %3d" % (STOP[4],iter,maxIter)
|
||||
print "%s DONE! %s\n" % ('='*25,'='*25)
|
||||
|
||||
return xc
|
||||
|
||||
def Rosenbrock(x):
|
||||
"""
|
||||
Rosenbrock function for testing GaussNewton scheme
|
||||
"""
|
||||
J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
|
||||
dJ = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1),200*(x[1]-x[0]**2)])
|
||||
H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]],[ -400*x[0], 200]],dtype=float);
|
||||
|
||||
return J,dJ,H
|
||||
|
||||
def checkDerivative(fctn,x0):
|
||||
"""
|
||||
Basic derivative check
|
||||
|
||||
Compares error decay of 0th and 1st order Taylor approximation at point
|
||||
x0 for a randomized search direction.
|
||||
|
||||
Input:
|
||||
------
|
||||
fctn - function handle
|
||||
x0 - point at which to check derivative
|
||||
"""
|
||||
|
||||
print "%s checkDerivative %s" % ('='*20,'='*20)
|
||||
print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|"
|
||||
|
||||
Jc,dJ,H = fctn(x0)
|
||||
|
||||
dx = np.random.randn(len(x0),1)
|
||||
|
||||
t = np.logspace(-1,-10,10)
|
||||
E0 = np.zeros(t.shape)
|
||||
E1 = np.zeros(t.shape)
|
||||
|
||||
for i in range(0,10):
|
||||
Jt = fctn(x0+t[i]*dx)
|
||||
E0[i] = norm(Jt[0]-Jc[0]) # 0th order Taylor
|
||||
E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx)) # 1st order Taylor
|
||||
|
||||
print "%d\t%1.2e\t%1.3e\t%1.3e" % (i,t[i],E0[i],E1[i])
|
||||
|
||||
print "%s DONE! %s\n" % ('='*25,'='*25)
|
||||
plt.figure()
|
||||
plt.clf()
|
||||
plt.loglog(t,E0,'b')
|
||||
plt.loglog(t,E1,'g--')
|
||||
plt.title('checkDerivative')
|
||||
plt.xlabel('h')
|
||||
plt.ylabel('error of Taylor approximation')
|
||||
plt.legend(['0th order', '1st order'],loc='upper left')
|
||||
plt.show()
|
||||
return
|
||||
|
||||
if __name__ == '__main__':
|
||||
x0 = np.array([[2.6],[3.7]])
|
||||
fctn = lambda x:Rosenbrock(x)
|
||||
checkDerivative(fctn,x0)
|
||||
xOpt = GaussNewton(fctn,x0,maxIter=20)
|
||||
print "xOpt=[%f,%f]" % (xOpt[0],xOpt[1])
|
||||
|
||||
@@ -0,0 +1,36 @@
|
||||
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
|
||||
@@ -0,0 +1,58 @@
|
||||
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
|
||||
@@ -0,0 +1,106 @@
|
||||
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)
|
||||
@@ -0,0 +1,213 @@
|
||||
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)
|
||||
@@ -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,36 @@
|
||||
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,94 @@
|
||||
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)))
|
||||
@@ -0,0 +1,77 @@
|
||||
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))
|
||||
@@ -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;
|
||||
@@ -0,0 +1,87 @@
|
||||
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
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -0,0 +1,213 @@
|
||||
import numpy as np
|
||||
from BaseMesh import BaseMesh
|
||||
from TensorView import TensorView
|
||||
from utils import ndgrid
|
||||
|
||||
|
||||
class TensorMesh(BaseMesh, TensorView):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
|
||||
Any Mesh that has a constant width along the entire axis
|
||||
such that it can defined by a single width vector, called 'h'.
|
||||
|
||||
e.g.
|
||||
|
||||
hx = np.array([1,1,1])
|
||||
hy = np.array([1,2])
|
||||
hz = np.array([1,1,1,1])
|
||||
|
||||
mesh = TensorMesh([hx, hy, hz])
|
||||
|
||||
"""
|
||||
def __init__(self, h, x0=None):
|
||||
super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0)
|
||||
|
||||
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
|
||||
|
||||
for i, h_i in enumerate(h):
|
||||
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
|
||||
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [x.ravel() for x in h]
|
||||
|
||||
def h():
|
||||
doc = "h is a list containing the cell widths of the tensor mesh in each dimension."
|
||||
fget = lambda self: self._h
|
||||
return locals()
|
||||
h = property(**h())
|
||||
|
||||
def hx():
|
||||
doc = "Width of cells in the x direction"
|
||||
fget = lambda self: self._h[0]
|
||||
return locals()
|
||||
hx = property(**hx())
|
||||
|
||||
def hy():
|
||||
doc = "Width of cells in the y direction"
|
||||
fget = lambda self: None if self.dim < 2 else self._h[1]
|
||||
return locals()
|
||||
hy = property(**hy())
|
||||
|
||||
def hz():
|
||||
doc = "Width of cells in the z direction"
|
||||
fget = lambda self: None if self.dim < 3 else self._h[2]
|
||||
return locals()
|
||||
hz = property(**hz())
|
||||
|
||||
def vectorNx():
|
||||
doc = "Nodal grid vector (1D) in the x direction."
|
||||
fget = lambda self: np.r_[0., self.hx.cumsum()] + self.x0[0]
|
||||
return locals()
|
||||
vectorNx = property(**vectorNx())
|
||||
|
||||
def vectorNy():
|
||||
doc = "Nodal grid vector (1D) in the y direction."
|
||||
fget = lambda self: None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1]
|
||||
return locals()
|
||||
vectorNy = property(**vectorNy())
|
||||
|
||||
def vectorNz():
|
||||
doc = "Nodal grid vector (1D) in the z direction."
|
||||
fget = lambda self: None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2]
|
||||
return locals()
|
||||
vectorNz = property(**vectorNz())
|
||||
|
||||
def vectorCCx():
|
||||
doc = "Cell-centered grid vector (1D) in the x direction."
|
||||
fget = lambda self: np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0]
|
||||
return locals()
|
||||
vectorCCx = property(**vectorCCx())
|
||||
|
||||
def vectorCCy():
|
||||
doc = "Cell-centered grid vector (1D) in the y direction."
|
||||
fget = lambda self: None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1]
|
||||
return locals()
|
||||
vectorCCy = property(**vectorCCy())
|
||||
|
||||
def vectorCCz():
|
||||
doc = "Cell-centered grid vector (1D) in the z direction."
|
||||
fget = lambda self: None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2]
|
||||
return locals()
|
||||
vectorCCz = property(**vectorCCz())
|
||||
|
||||
def gridCC():
|
||||
doc = "Cell-centered grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridCC is None:
|
||||
self._gridCC = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorCCz] if not x is None])
|
||||
return self._gridCC
|
||||
return locals()
|
||||
_gridCC = None # Store grid by default
|
||||
gridCC = property(**gridCC())
|
||||
|
||||
def gridN():
|
||||
doc = "Nodal grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridN is None:
|
||||
self._gridN = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorNz] if not x is None])
|
||||
return self._gridN
|
||||
return locals()
|
||||
_gridN = None # Store grid by default
|
||||
gridN = property(**gridN())
|
||||
|
||||
def gridFx():
|
||||
doc = "Face staggered grid in the x direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFx is None:
|
||||
self._gridFx = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorCCz] if not x is None])
|
||||
return self._gridFx
|
||||
return locals()
|
||||
_gridFx = None # Store grid by default
|
||||
gridFx = property(**gridFx())
|
||||
|
||||
def gridFy():
|
||||
doc = "Face staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFy is None:
|
||||
self._gridFy = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorCCz] if not x is None])
|
||||
return self._gridFy
|
||||
return locals()
|
||||
_gridFy = None # Store grid by default
|
||||
gridFy = property(**gridFy())
|
||||
|
||||
def gridFz():
|
||||
doc = "Face staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFz is None:
|
||||
self._gridFz = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorNz] if not x is None])
|
||||
return self._gridFz
|
||||
return locals()
|
||||
_gridFz = None # Store grid by default
|
||||
gridFz = property(**gridFz())
|
||||
|
||||
def gridEx():
|
||||
doc = "Edge staggered grid in the x direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEx is None:
|
||||
self._gridEx = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorNz] if not x is None])
|
||||
return self._gridEx
|
||||
return locals()
|
||||
_gridEx = None # Store grid by default
|
||||
gridEx = property(**gridEx())
|
||||
|
||||
def gridEy():
|
||||
doc = "Edge staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEy is None:
|
||||
self._gridEy = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorNz] if not x is None])
|
||||
return self._gridEy
|
||||
return locals()
|
||||
_gridEy = None # Store grid by default
|
||||
gridEy = property(**gridEy())
|
||||
|
||||
def gridEz():
|
||||
doc = "Edge staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEz is None:
|
||||
self._gridEz = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorCCz] if not x is None])
|
||||
return self._gridEz
|
||||
return locals()
|
||||
_gridEz = None # Store grid by default
|
||||
gridEz = property(**gridEz())
|
||||
|
||||
def getBoundaryIndex(self, gridType):
|
||||
"""Needed for faces edges and cells"""
|
||||
pass
|
||||
|
||||
def getCellNumbering(self):
|
||||
pass
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print('Welcome to tensor mesh!')
|
||||
|
||||
testDim = 1
|
||||
h1 = 0.3*np.ones((1, 7))
|
||||
h1[:, 0] = 0.5
|
||||
h1[:, -1] = 0.6
|
||||
h2 = .5 * np.ones((1, 4))
|
||||
h3 = .4 * np.ones((1, 6))
|
||||
x0 = np.zeros((3, 1))
|
||||
|
||||
if testDim == 1:
|
||||
h = [h1]
|
||||
x0 = x0[0]
|
||||
elif testDim == 2:
|
||||
h = [h1, h2]
|
||||
x0 = x0[0:2]
|
||||
else:
|
||||
h = [h1, h2, h3]
|
||||
|
||||
I = np.linspace(0, 1, 8)
|
||||
M = TensorMesh(h, x0)
|
||||
|
||||
xn = M.plotGrid()
|
||||
@@ -0,0 +1,124 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
|
||||
class TensorView(object):
|
||||
"""
|
||||
Provides viewing functions for TensorMesh
|
||||
|
||||
This class is inherited by TensorMesh
|
||||
"""
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def plotImage(self, I, imageType='CC', figNum=1,ax=None):
|
||||
|
||||
assert type(I) == np.ndarray, "I must be a numpy array"
|
||||
assert imageType in ["CC", "N"], "imageType must be 'CC' or 'N'"
|
||||
|
||||
if imageType == 'CC':
|
||||
assert I.size == self.nC, "Incorrect dimensions for CC."
|
||||
elif imageType == 'N':
|
||||
assert I.size == self.nN, "Incorrect dimensions for N."
|
||||
|
||||
if ax is None:
|
||||
fig = plt.figure(figNum)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
else:
|
||||
assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!"
|
||||
fig = ax.figure
|
||||
|
||||
if self.dim == 1:
|
||||
if imageType == 'CC':
|
||||
ph = ax.plot(self.vectorCCx, I, '-ro')
|
||||
elif imageType == 'N':
|
||||
ph = ax.plot(self.vectorNx, I, '-bs')
|
||||
ax.set_xticks(self.vectorNx)
|
||||
ax.set_xlabel("x")
|
||||
ax.axis('tight')
|
||||
elif self.dim == 2:
|
||||
if imageType == 'CC':
|
||||
C = I[:].reshape(self.n, order='F')
|
||||
elif imageType == 'N':
|
||||
C = I[:].reshape(self.n+1, order='F')
|
||||
C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:])
|
||||
|
||||
ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T)
|
||||
ax.axis('tight')
|
||||
ax.set_xlabel("x")
|
||||
ax.set_ylabel("y")
|
||||
ax.set_xticks(self.vectorNx)
|
||||
ax.set_yticks(self.vectorNy)
|
||||
|
||||
|
||||
fig.show()
|
||||
return ph
|
||||
|
||||
|
||||
def plotGrid(self):
|
||||
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions."""
|
||||
if self.dim == 1:
|
||||
fig = plt.figure(1)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
xn = self.gridN
|
||||
xc = self.gridCC
|
||||
ax.hold(True)
|
||||
ax.plot(xn, np.ones(np.shape(xn)), 'bs')
|
||||
ax.plot(xc, np.ones(np.shape(xc)), 'ro')
|
||||
ax.plot(xn, np.ones(np.shape(xn)), 'k--')
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
fig.show()
|
||||
elif self.dim == 2:
|
||||
fig = plt.figure(2)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
xn = self.gridN
|
||||
xc = self.gridCC
|
||||
xs1 = self.gridFx
|
||||
xs2 = self.gridFy
|
||||
|
||||
ax.hold(True)
|
||||
ax.plot(xn[:, 0], xn[:, 1], 'bs')
|
||||
ax.plot(xc[:, 0], xc[:, 1], 'ro')
|
||||
ax.plot(xs1[:, 0], xs1[:, 1], 'g>')
|
||||
ax.plot(xs2[:, 0], xs2[:, 1], 'g^')
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
fig.show()
|
||||
elif self.dim == 3:
|
||||
fig = plt.figure(3)
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
xn = self.gridN
|
||||
xc = self.gridCC
|
||||
xfs1 = self.gridFx
|
||||
xfs2 = self.gridFy
|
||||
xfs3 = self.gridFz
|
||||
|
||||
xes1 = self.gridEx
|
||||
xes2 = self.gridEy
|
||||
xes3 = self.gridEz
|
||||
|
||||
ax.hold(True)
|
||||
ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2])
|
||||
ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2])
|
||||
ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2])
|
||||
ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2])
|
||||
ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2])
|
||||
ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2])
|
||||
ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2])
|
||||
ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2])
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
ax.set_zlabel('x3')
|
||||
fig.show()
|
||||
@@ -0,0 +1 @@
|
||||
from TensorMesh import TensorMesh
|
||||
@@ -0,0 +1,75 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from utils import mkvc
|
||||
from sputils import ddx, sdiag, speye, kron3
|
||||
|
||||
def getvol(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# Compute cell volumes
|
||||
v12 = h1.T*h2
|
||||
V = mkvc(v12.reshape(-1,1)*h3)
|
||||
|
||||
return V
|
||||
|
||||
def getarea(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
# Compute areas of cell faces
|
||||
area1 = np.ones((n1+1,1))*mkvc(h2.T*h3)
|
||||
area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
|
||||
area3 = h1.T*mkvc(h2.T*np.ones(n3+1))
|
||||
area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
|
||||
return area
|
||||
|
||||
def getDivMatrix(h):
|
||||
"""Consturct the 3D divergence operator on Faces."""
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
|
||||
# Compute areas of cell faces
|
||||
#area1 = np.ones((n1+1,1))*mkvc(h2.T*h3)
|
||||
#area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
|
||||
#area3 = h1.T*mkvc(h2.T*np.ones(n3+1))
|
||||
#area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
area = getarea(h)
|
||||
|
||||
S = sdiag(area)
|
||||
|
||||
# Compute cell volumes
|
||||
#v12 = h1.T*h2
|
||||
#V = mkvc(v12.reshape(-1,1)*h3)
|
||||
V = getvol(h)
|
||||
|
||||
# Compute divergence operator on faces
|
||||
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))
|
||||
|
||||
D = sparse.hstack((sparse.hstack((D1, D2)), D3))
|
||||
return sdiag(1/V)*D*S
|
||||
|
||||
@@ -0,0 +1,18 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
|
||||
def ddx(n):
|
||||
"""Define 1D derivatives"""
|
||||
return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1)
|
||||
|
||||
def sdiag(h):
|
||||
"""Sparse diagonal matrix"""
|
||||
return sparse.spdiags(h, 0, np.size(h), np.size(h))
|
||||
|
||||
def speye(n):
|
||||
"""Sparse identity"""
|
||||
return sparse.identity(n)
|
||||
|
||||
def kron3(A, B, C):
|
||||
"""Two kron prods"""
|
||||
return sparse.kron(sparse.kron(A, B), C)
|
||||
@@ -0,0 +1,11 @@
|
||||
import glob
|
||||
import unittest
|
||||
|
||||
# This code will run all tests in directory named test_*.py
|
||||
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
text_runner = unittest.TextTestRunner().run(testSuite)
|
||||
@@ -0,0 +1,3 @@
|
||||
#!/bin/sh
|
||||
|
||||
python -m unittest discover
|
||||
@@ -0,0 +1,88 @@
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from BaseMesh import BaseMesh
|
||||
import numpy as np
|
||||
|
||||
|
||||
class TestBaseMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh = BaseMesh([6, 2, 3])
|
||||
|
||||
def test_meshDimensions(self):
|
||||
self.assertTrue(self.mesh.dim, 3)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.n == [6, 2, 3]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
x = np.all(self.mesh.nCx == 6)
|
||||
y = np.all(self.mesh.nCy == 2)
|
||||
z = np.all(self.mesh.nCz == 3)
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_nf(self):
|
||||
x = np.all(self.mesh.nFx == [7, 2, 3])
|
||||
y = np.all(self.mesh.nFy == [6, 3, 3])
|
||||
z = np.all(self.mesh.nFz == [6, 2, 4])
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_ne(self):
|
||||
x = np.all(self.mesh.nEx == [6, 3, 4])
|
||||
y = np.all(self.mesh.nEy == [7, 2, 4])
|
||||
z = np.all(self.mesh.nEz == [7, 3, 3])
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
c = self.mesh.nC == 36
|
||||
f = np.all(self.mesh.nF == [42, 54, 48])
|
||||
e = np.all(self.mesh.nE == [72, 56, 63])
|
||||
|
||||
self.assertTrue(np.all([c, f, e]))
|
||||
|
||||
|
||||
class TestMeshNumbers2D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh = BaseMesh([6, 2])
|
||||
|
||||
def test_meshDimensions(self):
|
||||
self.assertTrue(self.mesh.dim, 2)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.n == [6, 2]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
x = np.all(self.mesh.nCx == 6)
|
||||
y = np.all(self.mesh.nCy == 2)
|
||||
z = self.mesh.nCz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_nf(self):
|
||||
x = np.all(self.mesh.nFx == [7, 2])
|
||||
y = np.all(self.mesh.nFy == [6, 3])
|
||||
z = self.mesh.nFz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_ne(self):
|
||||
x = np.all(self.mesh.nEx == [6, 3])
|
||||
y = np.all(self.mesh.nEy == [7, 2])
|
||||
z = self.mesh.nEz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
c = self.mesh.nC == 12
|
||||
f = np.all(self.mesh.nF == [14, 18])
|
||||
e = np.all(self.mesh.nE == [18, 14])
|
||||
|
||||
self.assertTrue(np.all([c, f, e]))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,45 @@
|
||||
import numpy as np
|
||||
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
from getDIV import getDivMatrix, getarea, getvol
|
||||
|
||||
# Define the mesh
|
||||
err=0.
|
||||
for i in range(4):
|
||||
icount=i+1;
|
||||
nc = 2*icount;
|
||||
h1 = np.pi/nc*np.ones((1,nc))
|
||||
h2 = np.pi/nc*np.ones((1,nc))
|
||||
h3 = np.pi/nc*np.ones((1,nc))
|
||||
h = [h1, h2, h3]
|
||||
x0 = -np.pi/2*np.ones((3, 1))
|
||||
M = TensorMesh(h, x0)
|
||||
#n = M.plotGrid()
|
||||
|
||||
# Generate DIV matrix
|
||||
DIV = getDivMatrix(h)
|
||||
|
||||
#Test function
|
||||
fun = lambda x: np.sin(x)
|
||||
Fx = fun(M.gridFx[:,0])
|
||||
Fy = fun(M.gridFy[:,1])
|
||||
Fz = fun(M.gridFz[:,2])
|
||||
|
||||
F = np.concatenate((Fx,Fy,Fz))
|
||||
divF = DIV*F
|
||||
sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2])
|
||||
|
||||
area = getarea(h)
|
||||
vol = getvol(h)
|
||||
err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
|
||||
if icount == 1:
|
||||
err1 = err
|
||||
print 'h | 2 norm | error ratio'
|
||||
print '---------------------------------------'
|
||||
print '%6.4f | %8.2e |'% (h1[0,0], err)
|
||||
else:
|
||||
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err1/err)
|
||||
|
||||
@@ -0,0 +1,34 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
a = np.array([1, 1, 1])
|
||||
b = np.array([1, 2])
|
||||
x0 = np.array([3, 5])
|
||||
self.mesh2 = TensorMesh([a, b], x0)
|
||||
|
||||
def test_vectorN_2D(self):
|
||||
testNx = np.array([3, 4, 5, 6])
|
||||
testNy = np.array([5, 6, 8])
|
||||
|
||||
xtest = np.all(self.mesh2.vectorNx == testNx)
|
||||
ytest = np.all(self.mesh2.vectorNy == testNy)
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
def test_vectorCC_2D(self):
|
||||
testNx = np.array([3.5, 4.5, 5.5])
|
||||
testNy = np.array([5.5, 7])
|
||||
|
||||
xtest = np.all(self.mesh2.vectorCCx == testNx)
|
||||
ytest = np.all(self.mesh2.vectorCCy == testNy)
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,53 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from utils import mkvc, ndgrid
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.a = np.array([1, 2, 3])
|
||||
self.b = np.array([1, 2])
|
||||
self.c = np.array([1, 2, 3, 4])
|
||||
|
||||
def test_mkvc1(self):
|
||||
x = mkvc(self.a)
|
||||
self.assertTrue(x.shape, (3,))
|
||||
|
||||
def test_mkvc2(self):
|
||||
x = mkvc(self.a, 2)
|
||||
self.assertTrue(x.shape, (3, 1))
|
||||
|
||||
def test_mkvc3(self):
|
||||
x = mkvc(self.a, 3)
|
||||
self.assertTrue(x.shape, (3, 1, 1))
|
||||
|
||||
def test_ndgrid_2D(self):
|
||||
XY = ndgrid([self.a, self.b])
|
||||
|
||||
X1_test = np.array([1, 2, 3, 1, 2, 3])
|
||||
X2_test = np.array([1, 1, 1, 2, 2, 2])
|
||||
|
||||
xtest = np.all(XY[:, 0] == X1_test)
|
||||
ytest = np.all(XY[:, 1] == X2_test)
|
||||
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
def test_ndgrid_3D(self):
|
||||
XYZ = ndgrid([self.a, self.b, self.c])
|
||||
|
||||
X1_test = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3])
|
||||
X2_test = np.array([1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
|
||||
X3_test = np.array([1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4])
|
||||
|
||||
xtest = np.all(XYZ[:, 0] == X1_test)
|
||||
ytest = np.all(XYZ[:, 1] == X2_test)
|
||||
ztest = np.all(XYZ[:, 2] == X3_test)
|
||||
|
||||
self.assertTrue(xtest and ytest and ztest)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,99 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def reshapeF(x, size):
|
||||
return np.reshape(x, size, order='F')
|
||||
|
||||
|
||||
def mkvc(x, numDims=1):
|
||||
"""Creates a vector with the number of dimension specified
|
||||
|
||||
e.g.:
|
||||
|
||||
a = np.array([1, 2, 3])
|
||||
|
||||
mkvc(a, 1).shape
|
||||
> (3, )
|
||||
|
||||
mkvc(a, 2).shape
|
||||
> (3, 1)
|
||||
|
||||
mkvc(a, 3).shape
|
||||
> (3, 1, 1)
|
||||
|
||||
"""
|
||||
assert type(x) == np.ndarray, "Vector must be a numpy array"
|
||||
|
||||
if numDims == 1:
|
||||
return x.flatten(order='F')
|
||||
elif numDims == 2:
|
||||
return x.flatten(order='F')[:, np.newaxis]
|
||||
elif numDims == 3:
|
||||
return x.flatten(order='F')[:, np.newaxis, np.newaxis]
|
||||
|
||||
|
||||
def ndgrid(*args, **kwargs):
|
||||
"""
|
||||
Form tensorial grid for 1, 2, or 3 dimensions.
|
||||
|
||||
Returns as column vectors by default.
|
||||
|
||||
To return as matrix input:
|
||||
|
||||
ndgrid(..., vector=False)
|
||||
|
||||
The inputs can be a list or separate arguments.
|
||||
|
||||
e.g.
|
||||
|
||||
a = np.array([1, 2, 3])
|
||||
b = np.array([1, 2])
|
||||
|
||||
XY = ndgrid(a, b)
|
||||
> [[1 1]
|
||||
[2 1]
|
||||
[3 1]
|
||||
[1 2]
|
||||
[2 2]
|
||||
[3 2]]
|
||||
|
||||
X, Y = ndgrid(a, b, vector=False)
|
||||
> X = [[1 1]
|
||||
[2 2]
|
||||
[3 3]]
|
||||
> Y = [[1 2]
|
||||
[1 2]
|
||||
[1 2]]
|
||||
|
||||
"""
|
||||
|
||||
# Read the keyword arguments, and only accept a vector=True/False
|
||||
vector = kwargs.pop('vector', True)
|
||||
assert type(vector) == bool, "'vector' keyword must be a bool"
|
||||
assert len(kwargs) == 0, "Only 'vector' keyword accepted"
|
||||
|
||||
# you can either pass a list [x1, x2, x3] or each seperately
|
||||
if type(args[0]) == list:
|
||||
xin = args[0]
|
||||
else:
|
||||
xin = args
|
||||
|
||||
# Each vector needs to be a numpy array
|
||||
assert np.all([type(x) == np.ndarray for x in xin]), "All vectors must be numpy arrays."
|
||||
|
||||
if len(xin) == 1:
|
||||
return xin[0]
|
||||
elif len(xin) == 2:
|
||||
XY = np.broadcast_arrays(mkvc(xin[1], 1), mkvc(xin[0], 2))
|
||||
if vector:
|
||||
X2, X1 = [mkvc(x) for x in XY]
|
||||
return np.c_[X1, X2]
|
||||
else:
|
||||
return XY[1], XY[0]
|
||||
elif len(xin) == 3:
|
||||
XYZ = np.broadcast_arrays(mkvc(xin[2], 1), mkvc(xin[1], 2), mkvc(xin[0], 3))
|
||||
if vector:
|
||||
X3, X2, X1 = [mkvc(x) for x in XYZ]
|
||||
return np.c_[X1, X2, X3]
|
||||
else:
|
||||
return XYZ[2], XYZ[1], XYZ[0]
|
||||
Reference in New Issue
Block a user