mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 06:33:31 +08:00
rowanc1
@@ -308,6 +308,51 @@ class BaseMesh(object):
|
||||
return locals()
|
||||
nF = property(**nF())
|
||||
|
||||
def normals():
|
||||
doc = "Face Normals"
|
||||
|
||||
def fget(self):
|
||||
if self.dim == 2:
|
||||
nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0])]
|
||||
nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1])]
|
||||
return np.r_[nX, nY]
|
||||
elif self.dim == 3:
|
||||
nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0]), np.zeros(self.nF[0])]
|
||||
nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1]), np.zeros(self.nF[1])]
|
||||
nZ = np.c_[np.zeros(self.nF[2]), np.zeros(self.nF[2]), np.ones(self.nF[2])]
|
||||
return np.r_[nX, nY, nZ]
|
||||
return locals()
|
||||
normals = property(**normals())
|
||||
|
||||
def tangents():
|
||||
doc = "Edge Tangents"
|
||||
|
||||
def fget(self):
|
||||
if self.dim == 2:
|
||||
tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0])]
|
||||
tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1])]
|
||||
return np.r_[tX, tY]
|
||||
elif self.dim == 3:
|
||||
tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0]), np.zeros(self.nE[0])]
|
||||
tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1]), np.zeros(self.nE[1])]
|
||||
tZ = np.c_[np.zeros(self.nE[2]), np.zeros(self.nE[2]), np.ones(self.nE[2])]
|
||||
return np.r_[tX, tY, tZ]
|
||||
return locals()
|
||||
tangents = property(**tangents())
|
||||
|
||||
def projectFaceVector(self, fV):
|
||||
"""Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals"""
|
||||
assert type(fV) == np.ndarray, 'fV must be an ndarray'
|
||||
assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
|
||||
return np.sum(fV*self.normals, 1)
|
||||
|
||||
def projectEdgeVector(self, eV):
|
||||
"""Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents"""
|
||||
assert type(eV) == np.ndarray, 'eV must be an ndarray'
|
||||
assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
|
||||
return np.sum(eV*self.tangents, 1)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
m = BaseMesh([3, 2, 4])
|
||||
print m.n
|
||||
|
||||
+46
-4
@@ -50,7 +50,7 @@ class DiffOperators(object):
|
||||
Class creates the differential operators that you need!
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
|
||||
def faceDiv():
|
||||
doc = "Construct divergence operator (face-stg to cell-centres)."
|
||||
@@ -168,9 +168,9 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._edgeCurl is None):
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(self.hx)
|
||||
n2 = np.size(self.hy)
|
||||
n3 = np.size(self.hz)
|
||||
n1 = self.nCx
|
||||
n2 = self.nCy
|
||||
n3 = self.nCy
|
||||
|
||||
# Compute lengths of cell edges
|
||||
L = self.edge
|
||||
@@ -245,6 +245,48 @@ class DiffOperators(object):
|
||||
_edgeAve = None
|
||||
edgeAve = property(**edgeAve())
|
||||
|
||||
def nodalAve():
|
||||
doc = "Construct the averaging operator on cell nodes to cell centers."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalAve is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
self._nodalAve = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._nodalAve = sp.hstack((sp.kron(av(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._nodalAve = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])),
|
||||
kron3(av(n[2]), av(n[1]), av(n[0])),
|
||||
kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._nodalAve
|
||||
return locals()
|
||||
_nodalAve = None
|
||||
nodalAve = property(**nodalAve())
|
||||
|
||||
def nodalVectorAve():
|
||||
doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension seperate."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalVectorAve is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
self._nodalVectorAve = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._nodalVectorAve = sp.block_diag((sp.kron(av(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._nodalVectorAve = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])),
|
||||
kron3(av(n[2]), av(n[1]), av(n[0])),
|
||||
kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._nodalVectorAve
|
||||
return locals()
|
||||
_nodalVectorAve = None
|
||||
nodalVectorAve = property(**nodalVectorAve())
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matix for products of edge functions w'*M(materialProp)*e"""
|
||||
if(materialProp is None):
|
||||
|
||||
@@ -1,36 +0,0 @@
|
||||
import numpy as np
|
||||
from numpy.random import randn
|
||||
from utils import ndgrid
|
||||
from getDiffOps import getCurlMatrix, getNodalGradient
|
||||
from getFaceInnerProduct import getFaceInnerProduct
|
||||
from getEdgeInnerProduct import getEdgeInnerProduct
|
||||
from scipy.sparse.linalg import dsolve
|
||||
from pylab import norm
|
||||
|
||||
n = np.array([14, 14, 15])
|
||||
|
||||
X, Y, Z = ndgrid(*[np.linspace(0, 1, x) for x in n])
|
||||
sigma = 1e-2*np.ones(n-1)
|
||||
sigma[:, :, (n[2]-1)/2:] = 1e-6
|
||||
mu = 4*np.pi*1e-7*np.ones(n-1)
|
||||
w = 10
|
||||
|
||||
CURL = getCurlMatrix(X, Y, Z)
|
||||
GRAD = getNodalGradient(X, Y, Z)
|
||||
Mf = getFaceInnerProduct(X, Y, Z, 1/mu)
|
||||
Me = getEdgeInnerProduct(X, Y, Z, sigma)
|
||||
|
||||
A = CURL.T * Mf * CURL + 1j * w * Me
|
||||
|
||||
ne = np.shape(A)
|
||||
b = np.matrix(randn(ne[0])).T
|
||||
# clean b
|
||||
DIVb = GRAD.T*b
|
||||
p = dsolve.spsolve(GRAD.T*GRAD, DIVb, use_umfpack=True).T
|
||||
b = b - GRAD*p
|
||||
|
||||
#x = spsolve(A, b)
|
||||
x = dsolve.spsolve(A, b, use_umfpack=True).T
|
||||
|
||||
t = norm(A*x-b)/norm(b)
|
||||
print t
|
||||
@@ -1,58 +0,0 @@
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from sputils import *
|
||||
from numpy import *
|
||||
from getEdgeTangent import *
|
||||
|
||||
|
||||
def volTetra(y, m, I, A, B, C, D):
|
||||
|
||||
a11 = array(y[A, 0]-y[B, 0]); a12 = array(y[A, 0]-y[C, 0]); a13 = array(y[A, 0]-y[D, 0])
|
||||
a21 = array(y[A, 1]-y[B, 1]); a22 = array(y[A, 1]-y[C, 1]); a23 = array(y[A, 1]-y[D, 1])
|
||||
a31 = array(y[A, 2]-y[B, 2]); a32 = array(y[A, 2]-y[C, 2]); a33 = array(y[A, 2]-y[D, 2])
|
||||
|
||||
return abs(a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a31*a22*a13 - a32*a23*a11 - a33*a21*a12)
|
||||
|
||||
|
||||
def getCellVolume(X, Y, Z):
|
||||
|
||||
m = array(shape(X))-1
|
||||
y = hstack3(mkvc(X), mkvc(Y), mkvc(Z))
|
||||
|
||||
i = int64(linspace(0, m[0]-1, m[0]))
|
||||
j = int64(linspace(0, m[1]-1, m[1]))
|
||||
k = int64(linspace(0, m[2]-1, m[2]))
|
||||
|
||||
ii, jj, kk = ndgrid(i, j, k)
|
||||
ii = mkvc(ii)
|
||||
jj = mkvc(jj)
|
||||
kk = mkvc(kk)
|
||||
|
||||
I = int64(sub2ind(m, hstack3(ii, jj, kk)))
|
||||
A = int64(sub2ind(m+1, hstack3(ii, jj, kk)))
|
||||
B = int64(sub2ind(m+1, hstack3(ii, jj+1, kk)))
|
||||
C = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk)))
|
||||
D = int64(sub2ind(m+1, hstack3(ii+1, jj, kk)))
|
||||
E = int64(sub2ind(m+1, hstack3(ii, jj, kk+1)))
|
||||
F = int64(sub2ind(m+1, hstack3(ii, jj+1, kk+1)))
|
||||
G = int64(sub2ind(m+1, hstack3(ii+1, jj+1, kk+1)))
|
||||
H = int64(sub2ind(m+1, hstack3(ii+1, jj, kk+1)))
|
||||
|
||||
v1 = volTetra(y, m, I, A, B, D, E)
|
||||
v2 = volTetra(y, m, I, B, E, F, G)
|
||||
v3 = volTetra(y, m, I, B, D, E, G)
|
||||
v4 = volTetra(y, m, I, B, C, D, G)
|
||||
v5 = volTetra(y, m, I, D, E, G, H)
|
||||
|
||||
v = 1.0/6.0 * (v1 + v2 + v3 + v4 + v5)
|
||||
return v.flatten()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3))
|
||||
Z[2, 2, 2] = 2.5; Z[0, 0, 0] = -0.5
|
||||
X[2, 2, 2] = 2.5; X[0, 0, 0] = -0.5
|
||||
|
||||
v = getCellVolume(X, Y, Z)
|
||||
print v
|
||||
@@ -1,106 +0,0 @@
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from numpy import *
|
||||
from getEdgeTangent import *
|
||||
from getCellVolume import getCellVolume
|
||||
from getFaceNormals import getFaceNormals
|
||||
|
||||
|
||||
def getDivMatrix(X, Y, Z):
|
||||
"""Face DIV"""
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]
|
||||
n2 = n[1]
|
||||
n3 = n[2]
|
||||
|
||||
n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
S = sdiag(area)
|
||||
V = getCellVolume(X, Y, Z)
|
||||
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
D1 = kron3(speye(n3), speye(n2), d1)
|
||||
D2 = kron3(speye(n3), d2, speye(n1))
|
||||
D3 = kron3(d3, speye(n2), speye(n1))
|
||||
|
||||
# divergence on faces
|
||||
D = appendRight3(D1, D2, D3)
|
||||
|
||||
return sdiag(1/V)*D*S
|
||||
|
||||
|
||||
def getCurlMatrix(X, Y, Z):
|
||||
"""Edge CURL """
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]; n2 = n[1]; n3 = n[2]
|
||||
|
||||
d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3)
|
||||
# derivatives on x-edge variables
|
||||
D32 = kron3(d3, speye(n2), speye(n1+1))
|
||||
D23 = kron3(speye(n3), d2, speye(n1+1))
|
||||
D31 = kron3(d3, speye(n2+1), speye(n1))
|
||||
D13 = kron3(speye(n3), speye(n2+1), d1)
|
||||
D21 = kron3(speye(n3+1), d2, speye(n1))
|
||||
D12 = kron3(speye(n3+1), speye(n2), d1)
|
||||
|
||||
O1 = spzeros(shape(D32)[0], shape(D31)[1])
|
||||
O2 = spzeros(shape(D31)[0], shape(D32)[1])
|
||||
O3 = spzeros(shape(D21)[0], shape(D13)[1])
|
||||
|
||||
CURL = appendBottom3(
|
||||
appendRight3(O1, -D32, D23),
|
||||
appendRight3(D31, O2, -D13),
|
||||
appendRight3(-D21, D12, O3))
|
||||
|
||||
# scale for non-uniform mesh
|
||||
e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z)
|
||||
n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, area1, area2, area3 = getFaceNormals(X, Y, Z)
|
||||
|
||||
area = hstack((hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
S = sdiag(1/area)
|
||||
lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3)))
|
||||
L = sdiag(lngth)
|
||||
|
||||
return S*(CURL*L)
|
||||
|
||||
|
||||
def getNodalGradient(X, Y, Z):
|
||||
"""Nodal Gradients"""
|
||||
|
||||
n = array(shape(X))-1
|
||||
n1 = n[0]; n2 = n[1]; n3 = n[2]
|
||||
|
||||
D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1))
|
||||
D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1))
|
||||
D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1))
|
||||
|
||||
# topological gradient
|
||||
GRAD = appendBottom3(D1, D2, D3)
|
||||
|
||||
# scale for non-uniform mesh
|
||||
e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, norme1, norme2, norme3 = getEdgeTangent(X, Y, Z)
|
||||
lngth = hstack((hstack((mkvc(norme1), mkvc(norme2))), mkvc(norme3)))
|
||||
L = sdiag(1/lngth)
|
||||
|
||||
return L*GRAD
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = ndgrid(linspace(0, 2, 3), linspace(0, 2, 3), linspace(0, 2, 3))
|
||||
Z[2, 2, 2] = 2.5
|
||||
Z[0, 0, 0] = -0.5
|
||||
X[2, 2, 2] = 2.5
|
||||
X[0, 0, 0] = -0.5
|
||||
sig = ones([2, 2, 2])
|
||||
C = getCurlMatrix(X, Y, Z)
|
||||
|
||||
G = getNodalGradient(X, Y, Z)
|
||||
|
||||
tt = C*G
|
||||
print(tt)
|
||||
@@ -1,213 +0,0 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from sputils import *
|
||||
from numpy import *
|
||||
from getEdgeTangent import *
|
||||
from inv3X3BlockDiagonal import *
|
||||
from getCellVolume import getCellVolume
|
||||
|
||||
|
||||
def subarray(T, i1, i2, i3):
|
||||
return take(take(take(T, i1, 0), i2, 1), i3, 2)
|
||||
|
||||
|
||||
def getEdgeInnerProduct(X, Y, Z, sigma):
|
||||
"""A = getEdgeInnerProduct(X, Y, Z, sigma)
|
||||
|
||||
|
||||
node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
/ /
|
||||
/ / |
|
||||
edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
/ / |
|
||||
/ / |
|
||||
node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
| | |
|
||||
| | node(i+1,j+1,k+1)
|
||||
| | /
|
||||
edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k)
|
||||
| | /
|
||||
| | /
|
||||
| |/
|
||||
node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
no | node | e1 | e2 | e3
|
||||
000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
|
||||
011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
|
||||
111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
|
||||
"""
|
||||
m = array(shape(X))-1
|
||||
nc = prod(m)
|
||||
|
||||
me1 = m + array([0, 1, 1]); ne1 = prod(me1)
|
||||
me2 = m + array([1, 0, 1]); ne2 = prod(me2)
|
||||
me3 = m + array([1, 1, 0]); ne3 = prod(me3)
|
||||
|
||||
e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,norme1,norme2,norme3 = getEdgeTangent(X,Y,Z)
|
||||
|
||||
i = int64(linspace(0,m[0]-1,m[0]))
|
||||
j = int64(linspace(0,m[1]-1,m[1]))
|
||||
k = int64(linspace(0,m[2]-1,m[2]))
|
||||
|
||||
ii,jj,kk = ndgrid(i,j,k)
|
||||
ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk)
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P000 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT000 = inv3X3BlockDiagonal(subarray(e1x,i,j,k) , subarray(e1y,i,j,k), subarray(e1z,i,j,k),
|
||||
subarray(e2x,i,j,k) , subarray(e2y,i,j,k), subarray(e2z,i,j,k) ,
|
||||
subarray(e3x,i,j,k) , subarray(e3y,i,j,k), subarray(e3z,i,j,k) )
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P100 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT100 = inv3X3BlockDiagonal(subarray(e1x,i,j,k), subarray(e1y,i,j,k), subarray(e1z,i,j,k),
|
||||
subarray(e2x,i+1,j,k), subarray(e2y,i+1,j,k), subarray(e2z,i+1,j,k),
|
||||
subarray(e3x,i+1,j,k) , subarray(e3y,i+1,j,k), subarray(e3z,i+1,j,k))
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P010 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT010 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) , subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) ,
|
||||
subarray(e2x,i,j,k) , subarray(e2y,i,j,k) , subarray(e2z,i,j,k) ,
|
||||
subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) ,subarray( e3z,i,j+1,k) )
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P110 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT110 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k) ,subarray(e1y,i,j+1,k) , subarray(e1z,i,j+1,k) ,
|
||||
subarray(e2x,i+1,j,k) ,subarray(e2y,i+1,j,k) , subarray(e2z,i+1,j,k),
|
||||
subarray(e3x,i+1,j+1,k) ,subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) )
|
||||
|
||||
######
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P001 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT001 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1) ,subarray(e1y,i,j,k+1) , subarray(e1z,i,j,k+1) ,
|
||||
subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) ,
|
||||
subarray(e3x,i,j,k) , subarray(e3y,i,j,k) , subarray(e3z,i,j,k) )
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P101 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT101 = inv3X3BlockDiagonal(subarray(e1x,i,j,k+1), subarray(e1y,i,j,k+1), subarray(e1z,i,j,k+1) ,
|
||||
subarray(e2x,i+1,j,k+1), subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) ,
|
||||
subarray(e3x,i+1,j,k), subarray(e3y,i+1,j,k) , subarray(e3z,i+1,j,k) )
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P011 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT011 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) ,
|
||||
subarray(e2x,i,j,k+1) , subarray(e2y,i,j,k+1) , subarray(e2z,i,j,k+1) ,
|
||||
subarray(e3x,i,j+1,k) , subarray(e3y,i,j+1,k) , subarray(e3z,i,j+1,k) )
|
||||
|
||||
## --------
|
||||
# no | node | e1 | e2 | e3
|
||||
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k+1
|
||||
ind1 = sub2ind(me1,hstack3(ii,jj+1,kk+1))
|
||||
ind2 = sub2ind(me2,hstack3(ii+1,jj,kk+1)) + ne1
|
||||
ind3 = sub2ind(me3,hstack3(ii+1,jj+1,kk)) + ne1 + ne2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P111 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
|
||||
|
||||
invT111 = inv3X3BlockDiagonal(subarray(e1x,i,j+1,k+1) , subarray(e1y,i,j+1,k+1) , subarray(e1z,i,j+1,k+1) ,
|
||||
subarray(e2x,i+1,j,k+1) , subarray(e2y,i+1,j,k+1) , subarray(e2z,i+1,j,k+1) ,
|
||||
subarray(e3x,i+1,j+1,k) , subarray(e3y,i+1,j+1,k) , subarray(e3z,i+1,j+1,k) )
|
||||
|
||||
# Cell volume
|
||||
v = mkvc(getCellVolume(X,Y,Z)) #mkvc(getVolume(X,Y,Z))
|
||||
vsig = v*mkvc(sigma)
|
||||
v3 = vstack((vstack((vsig,vsig)),vsig))
|
||||
v3 = v3.flatten()
|
||||
|
||||
V = sdiag(v3)
|
||||
|
||||
A = P000.T*invT000.T*V*invT000*P000 + P001.T*invT001.T*V*invT001*P001 + P010.T*invT010.T*V*invT010*P010 + P011.T*invT011.T*V*invT011*P011 + P100.T*invT100.T*V*invT100*P100 + P101.T*invT101.T*V*invT101*P101 + P110.T*invT110.T*V*invT110*P110 + P111.T*invT111.T*V*invT111*P111
|
||||
|
||||
A = 0.125*A
|
||||
|
||||
return A
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3))
|
||||
Z[2,2,2] = 2.5; Z[0,0,0] = -0.5
|
||||
X[2,2,2] = 2.5; X[0,0,0] = -0.5
|
||||
sig = ones([2,2,2])
|
||||
A = getEdgeInnerProduct(X,Y,Z,sig)
|
||||
@@ -1,60 +0,0 @@
|
||||
from numpy import *
|
||||
from utils import diff
|
||||
|
||||
#function[t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,normt1,normt2,normt3] = getEdgeTangent(X,Y,Z)
|
||||
#%[t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z,normt1,normt2,normt3] = getEdgeTangent(X,Y,Z)
|
||||
#%
|
||||
#% node(i,j,k+1) ------ edgt2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
#% / /
|
||||
#% / / |
|
||||
#% edgt3(i,j,k) fact1(i,j,k) edgt3(i,j+1,k)
|
||||
#% / / |
|
||||
#% / / |
|
||||
#% node(i,j,k) ------ edgt2(i,j,k) ----- node(i,j+1,k)
|
||||
#% | | |
|
||||
#% | | node(i+1,j+1,k+1)
|
||||
#% | | /
|
||||
#% edgt1(i,j,k) fact3(i,j,k) edgt1(i,j+1.k)
|
||||
#% | | /
|
||||
#% | | /
|
||||
#% | |/
|
||||
#% node(i+1,j,k) ------ edgt2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
|
||||
def getEdgeTangent(X, Y, Z):
|
||||
|
||||
t1x = diff(X, 1)
|
||||
t1y = diff(Y, 1)
|
||||
t1z = diff(Z, 1)
|
||||
|
||||
normt1 = sqrt(t1x**2+t1y**2+t1z**2)
|
||||
t1x = t1x/normt1
|
||||
t1y = t1y/normt1
|
||||
t1z = t1z/normt1
|
||||
|
||||
t2x = diff(X, 2)
|
||||
t2y = diff(Y, 2)
|
||||
t2z = diff(Z, 2)
|
||||
normt2 = sqrt(t2x**2 + t2y**2 + t2z**2)
|
||||
t2x = t2x/normt2
|
||||
t2y = t2y/normt2
|
||||
t2z = t2z/normt2
|
||||
|
||||
t3x = diff(X, 3)
|
||||
t3y = diff(Y, 3)
|
||||
t3z = diff(Z, 3)
|
||||
normt3 = sqrt(t3x**2+t3y**2+t3z**2)
|
||||
t3x = t3x/normt3
|
||||
t3y = t3y/normt3
|
||||
t3z = t3z/normt3
|
||||
|
||||
# print t3x
|
||||
|
||||
return (t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z, normt1, normt2, normt3)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = mgrid[0:4, 0:5, 0:6]
|
||||
|
||||
t = getEdgeTangent(X, Y, Z)
|
||||
@@ -1,86 +0,0 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from numpy import *
|
||||
from getEdgeTangent import *
|
||||
from inv3X3BlockDiagonal import *
|
||||
from getCellVolume import getCellVolume
|
||||
from getFaceNormals import getFaceNormals
|
||||
|
||||
|
||||
#-----------------------
|
||||
def subarray(T,i1,i2,i3):
|
||||
return take(take(take(T,i1,0),i2,1),i3,2)
|
||||
|
||||
#-----------------------
|
||||
|
||||
def getFaceInnerProduct(X,Y,Z,sigma):
|
||||
|
||||
m = array(shape(X))-1
|
||||
nc = prod(m)
|
||||
mf1 = m+[1, 0, 0]
|
||||
mf2 = m+[0, 1, 0]
|
||||
mf3 = m+[0, 0, 1]
|
||||
|
||||
nf1 = prod(m+[1, 0, 0])
|
||||
nf2 = prod(m+[0, 1, 0])
|
||||
nf3 = prod(m+[0, 0, 1])
|
||||
|
||||
# compute the normals
|
||||
n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3 = getFaceNormals(X,Y,Z)
|
||||
|
||||
i = int64(linspace(0,m[0]-1,m[0]))
|
||||
j = int64(linspace(0,m[1]-1,m[1]))
|
||||
k = int64(linspace(0,m[2]-1,m[2]))
|
||||
|
||||
ii,jj,kk = ndgrid(i,j,k)
|
||||
ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk)
|
||||
|
||||
ind1 = sub2ind(mf1,hstack3(ii,jj,kk))
|
||||
ind2 = sub2ind(mf2,hstack3(ii,jj,kk)) + nf1
|
||||
ind3 = sub2ind(mf3,hstack3(ii,jj,kk)) + nf1 + nf2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P1 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,nf1+nf2+nf3)).tocsr()
|
||||
|
||||
ind1 = sub2ind(mf1,hstack3(ii+1,jj,kk))
|
||||
ind2 = sub2ind(mf2,hstack3(ii,jj+1,kk)) + nf1
|
||||
ind3 = sub2ind(mf3,hstack3(ii,jj,kk+1)) + nf1 + nf2
|
||||
|
||||
IND = vstack((vstack((ind1,ind2)),ind3))
|
||||
IND = array(IND).flatten()
|
||||
|
||||
P2 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,nf1+nf2+nf3)).tocsr()
|
||||
|
||||
|
||||
invN1 = inv3X3BlockDiagonal(subarray(n1x,i,j,k) , subarray(n1y,i,j,k), subarray(n1z,i,j,k),
|
||||
subarray(n2x,i,j,k) , subarray(n2y,i,j,k), subarray(n2z,i,j,k),
|
||||
subarray(n3x,i,j,k) , subarray(n3y,i,j,k), subarray(n3z,i,j,k) )
|
||||
|
||||
|
||||
invN2 = inv3X3BlockDiagonal(subarray(n1x,i+1,j,k) , subarray(n1y,i+1,j,k), subarray(n1z,i+1,j,k),
|
||||
subarray(n2x,i,j+1,k) , subarray(n2y,i,j+1,k), subarray(n2z,i,j+1,k),
|
||||
subarray(n3x,i,j,k+1) , subarray(n3y,i,j,k+1), subarray(n3z,i,j,k+1) )
|
||||
|
||||
# Cell volume
|
||||
v = mkvc(getCellVolume(X,Y,Z)) #mkvc(getVolume(X,Y,Z))
|
||||
vsig = v*mkvc(sigma)
|
||||
v3 = vstack((vstack((vsig,vsig)),vsig))
|
||||
v3 = v3.flatten()
|
||||
|
||||
V = sdiag(v3)
|
||||
|
||||
return (P1.T*invN1.T*V*invN1*P1 + P2.T*invN2.T*V*invN2*P2)/2.0
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X,Y,Z = ndgrid(linspace(0,2,3),linspace(0,2,3),linspace(0,2,3))
|
||||
Z[2,2,2] = 2.5; Z[0,0,0] = -0.5
|
||||
X[2,2,2] = 2.5; X[0,0,0] = -0.5
|
||||
sigma = ones([2,2,2])
|
||||
A = getFaceInnerProduct(X,Y,Z,sigma)
|
||||
print(A)
|
||||
@@ -1,73 +0,0 @@
|
||||
from numpy import *
|
||||
from utils import *
|
||||
|
||||
|
||||
def getFaceNormals(X, Y, Z):
|
||||
# compute the x normals
|
||||
d1xp = diffp(X,2,3)
|
||||
d1yp = diffp(Y,2,3)
|
||||
d1zp = diffp(Z,2,3)
|
||||
|
||||
d1xm = diffm(X,3,2)
|
||||
d1ym = diffm(Y,3,2)
|
||||
d1zm = diffm(Z,3,2)
|
||||
|
||||
# normals
|
||||
n1x = d1yp*d1zm - d1zp*d1ym
|
||||
n1y = d1zp*d1xm - d1xp*d1zm
|
||||
n1z = d1xp*d1ym - d1yp*d1xm
|
||||
normn1 = sqrt(n1x**2 + n1y**2 + n1z**2)
|
||||
n1x = n1x / normn1;
|
||||
n1y = n1y / normn1;
|
||||
n1z = n1z / normn1;
|
||||
|
||||
area1 = normn1/2
|
||||
|
||||
|
||||
# compute the y normals
|
||||
d2xp = diffp(X,1,3)
|
||||
d2yp = diffp(Y,1,3)
|
||||
d2zp = diffp(Z,1,3)
|
||||
|
||||
d2xm = diffm(X,1,3)
|
||||
d2ym = diffm(Y,1,3)
|
||||
d2zm = diffm(Z,1,3)
|
||||
|
||||
# normals
|
||||
n2x = d2yp*d2zm - d2zp*d2ym
|
||||
n2y = d2zp*d2xm - d2xp*d2zm
|
||||
n2z = d2xp*d2ym - d2yp*d2xm
|
||||
normn2 = sqrt(n2x**2 + n2y**2 + n2z**2)
|
||||
n2x = n2x / normn2
|
||||
n2y = n2y / normn2
|
||||
n2z = n2z / normn2
|
||||
|
||||
area2 = normn2/2
|
||||
|
||||
# compute the z normals
|
||||
d3xp = diffp(X,1,2)
|
||||
d3yp = diffp(Y,1,2)
|
||||
d3zp = diffp(Z,1,2)
|
||||
|
||||
d3xm = diffm(X,2,1)
|
||||
d3ym = diffm(Y,2,1)
|
||||
d3zm = diffm(Z,2,1)
|
||||
|
||||
# normals
|
||||
n3x = d3yp*d3zm - d3zp*d3ym
|
||||
n3y = d3zp*d3xm - d3xp*d3zm
|
||||
n3z = d3xp*d3ym - d3yp*d3xm;
|
||||
normn3 = sqrt(n3x**2 + n3y**2 + n3z**2);
|
||||
n3x = n3x / normn3;
|
||||
n3y = n3y / normn3;
|
||||
n3z = n3z / normn3;
|
||||
|
||||
area3 = normn3/2;
|
||||
|
||||
return (n1x,n1y,n1z,n2x,n2y,n2z,n3x,n3y,n3z,area1,area2,area3)
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = mgrid[0:4, 0:5, 0:6]
|
||||
|
||||
t = getFaceNormals(X, Y, Z)
|
||||
@@ -1,34 +0,0 @@
|
||||
from numpy import *
|
||||
from utils import diff, ave
|
||||
|
||||
|
||||
def getVolume(X,Y,Z):
|
||||
|
||||
# compute edge vectors
|
||||
t1x = ave(ave(diff(X, 1),2),3)
|
||||
t1y = ave(ave(diff(Y, 1),2),3)
|
||||
t1z = ave(ave(diff(Z, 1),2),3)
|
||||
|
||||
t2x = ave(ave(diff(X, 2),1),3)
|
||||
t2y = ave(ave(diff(Y, 2),1),3)
|
||||
t2z = ave(ave(diff(Z, 2),1),3)
|
||||
|
||||
t3x = ave(ave(diff(X, 3),1),2)
|
||||
t3y = ave(ave(diff(Y, 3),1),2)
|
||||
t3z = ave(ave(diff(Z, 3),1),2)
|
||||
|
||||
# v = [t1x t1y t1z][i j k]
|
||||
# [t2x t2y t2z]
|
||||
# [t3x t3y t3z]
|
||||
|
||||
v = t1x*(t2y*t3z - t2z*t3y) - t1y*(t2x*t3z - t2z*t3x) + t1z*(t2x*t3y-t2y*t3x)
|
||||
|
||||
return v
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
X, Y, Z = mgrid[0:4, 0:5, 0:6]
|
||||
X = (1.0*X)/2
|
||||
v = getVolume(X, Y, Z)
|
||||
print v
|
||||
@@ -1,36 +0,0 @@
|
||||
from utils import *
|
||||
from sputils import *
|
||||
|
||||
|
||||
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
|
||||
|
||||
a11 = mkvc(a11)
|
||||
a12 = mkvc(a12)
|
||||
a13 = mkvc(a13)
|
||||
a21 = mkvc(a21)
|
||||
a22 = mkvc(a22)
|
||||
a23 = mkvc(a23)
|
||||
a31 = mkvc(a31)
|
||||
a32 = mkvc(a32)
|
||||
a33 = mkvc(a33)
|
||||
|
||||
detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
|
||||
|
||||
b11 = +(a22*a33 - a23*a32)/detA
|
||||
b12 = -(a12*a33 - a13*a32)/detA
|
||||
b13 = +(a12*a23 - a13*a22)/detA
|
||||
|
||||
b21 = +(a31*a23 - a21*a33)/detA
|
||||
b22 = -(a31*a13 - a11*a33)/detA
|
||||
b23 = +(a21*a13 - a11*a23)/detA
|
||||
|
||||
b31 = -(a31*a22 - a21*a32)/detA
|
||||
b32 = +(a31*a12 - a11*a32)/detA
|
||||
b33 = -(a21*a12 - a11*a22)/detA
|
||||
|
||||
B = appendBottom3(
|
||||
appendRight3(sdiag(b11), sdiag(b12), sdiag(b13)),
|
||||
appendRight3(sdiag(b21), sdiag(b22), sdiag(b23)),
|
||||
appendRight3(sdiag(b31), sdiag(b32), sdiag(b33)))
|
||||
|
||||
return B
|
||||
@@ -1,94 +0,0 @@
|
||||
from sputils import *
|
||||
from utils import *
|
||||
from numpy import *
|
||||
|
||||
|
||||
def getCellCenterFromNodal(X, Y, Z):
|
||||
"""Cell Centers from Nodal locations"""
|
||||
XC = 1.0/8.0 * (X[:-1, :-1, :-1] + X[1:, :-1, :-1] + X[:-1, 1:, :-1] + X[1:, 1:, :-1] +
|
||||
X[:-1, :-1, 1:] + X[1:, :-1, 1:] + X[:-1, 1:, 1:] + X[1:, 1:, 1:])
|
||||
|
||||
YC = 1.0/8.0 * (Y[:-1, :-1, :-1] + Y[1:, :-1, :-1] + Y[:-1, 1:, :-1] + Y[1:, 1:, :-1] +
|
||||
Y[:-1, :-1, 1:] + Y[1:, :-1, 1:] + Y[:-1, 1:, 1:] + Y[1:, 1:, 1:])
|
||||
|
||||
ZC = 1.0/8.0 * (Z[:-1, :-1, :-1] + Z[1:, :-1, :-1] + Z[:-1, 1:, :-1] + Z[1:, 1:, :-1] +
|
||||
Z[:-1, :-1, 1:] + Z[1:, :-1, 1:] + Z[:-1, 1:, 1:] + Z[1:, 1:, 1:])
|
||||
|
||||
return (XC, YC, ZC)
|
||||
|
||||
|
||||
def getEdgesFromNodal(X, Y, Z):
|
||||
"""Edges from Nodal locations
|
||||
|
||||
node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
/ /
|
||||
/ / |
|
||||
edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
/ / |
|
||||
/ / |
|
||||
node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
| | |
|
||||
| | node(i+1,j+1,k+1)
|
||||
| | /
|
||||
edge1(i,j,k) face3(i,j,k) edge1(i,j+1.k)
|
||||
| | /
|
||||
| | /
|
||||
| |/
|
||||
node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
"""
|
||||
|
||||
XE1 = (X[1:, :, :]+X[:-1, :, :])/2.0
|
||||
YE1 = (Y[1:, :, :]+Y[:-1, :, :])/2.0
|
||||
ZE1 = (Z[1:, :, :]+Z[:-1, :, :])/2.0
|
||||
|
||||
XE2 = (X[:, 1:, :]+X[:, :-1, :])/2.0
|
||||
YE2 = (Y[:, 1:, :]+Y[:, :-1, :])/2.0
|
||||
ZE2 = (Z[:, 1:, :]+Z[:, :-1, :])/2.0
|
||||
|
||||
XE3 = (X[:, :, 1:]+X[:, :, :-1])/2.0
|
||||
YE3 = (Y[:, :, 1:]+Y[:, :, :-1])/2.0
|
||||
ZE3 = (Z[:, :, 1:]+Z[:, :, :-1])/2.0
|
||||
|
||||
return (XE1, YE1, ZE1, XE2, YE2, ZE2, XE3, YE3, ZE3)
|
||||
|
||||
|
||||
def getFacesFromNodal(X, Y, Z):
|
||||
"""Get faces from nodal --"""
|
||||
|
||||
XF1 = 1.0/4.0*(X[:, :-1, :-1]+X[:, 1:, :-1]+X[:, :-1, 1:]+X[:, 1:, 1:])
|
||||
YF1 = 1.0/4.0*(Y[:, :-1, :-1]+Y[:, 1:, :-1]+Y[:, :-1, 1:]+Y[:, 1:, 1:])
|
||||
ZF1 = 1.0/4.0*(Z[:, :-1, :-1]+Z[:, 1:, :-1]+Z[:, :-1, 1:]+Z[:, 1:, 1:])
|
||||
|
||||
XF2 = 1.0/4.0*(X[:-1, :, :-1]+X[1:, :, :-1]+X[:-1, :, 1:]+X[1:, :, 1:])
|
||||
YF2 = 1.0/4.0*(Y[:-1, :, :-1]+Y[1:, :, :-1]+Y[:-1, :, 1:]+Y[1:, :, 1:])
|
||||
ZF2 = 1.0/4.0*(Z[:-1, :, :-1]+Z[1:, :, :-1]+Z[:-1, :, 1:]+Z[1:, :, 1:])
|
||||
|
||||
XF3 = 1.0/4.0*(X[:-1, :-1, :]+X[1:, :-1, :]+X[:-1, 1:, :]+X[1:, 1:, :])
|
||||
YF3 = 1.0/4.0*(Y[:-1, :-1, :]+Y[1:, :-1, :]+Y[:-1, 1:, :]+Y[1:, 1:, :])
|
||||
ZF3 = 1.0/4.0*(Z[:-1, :-1, :]+Z[1:, :-1, :]+Z[:-1, 1:, :]+Z[1:, 1:, :])
|
||||
|
||||
return (XF1, YF1, ZF1, XF2, YF2, ZF2, XF3, YF3, ZF3)
|
||||
|
||||
|
||||
def projectEdgeVectorField(EV1, EV2, EV3, X, Y, Z):
|
||||
"""Project Edge vector field"""
|
||||
|
||||
t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z, nrm1, nrm2, nrm3 = getEdgeTangent(X, Y, Z)
|
||||
|
||||
E1 = EV1[:, 0]*mkvc(t1x) + EV1[:, 1]*mkvc(t1y) + EV1[:, 2]*mkvc(t1z)
|
||||
E2 = EV2[:, 0]*mkvc(t2x) + EV2[:, 1]*mkvc(t2y) + EV2[:, 2]*mkvc(t2z)
|
||||
E3 = EV3[:, 0]*mkvc(t3x) + EV3[:, 1]*mkvc(t3y) + EV3[:, 2]*mkvc(t3z)
|
||||
|
||||
return hstack((hstack((mkvc(E1), mkvc(E2))), mkvc(E3)))
|
||||
|
||||
|
||||
def projectFaceVectorField(FV1, FV2, FV3, X, Y, Z):
|
||||
"""Prolect Face vector field"""
|
||||
|
||||
n1x, n1y, n1z, n2x, n2y, n2z, n3x, n3y, n3z, ar1, ar2, ar3 = getFaceNormals(X, Y, Z)
|
||||
|
||||
F1 = FV1[:, 0]*mkvc(n1x) + FV1[:, 1]*mkvc(n1y) + FV1[:, 2]*mkvc(n1z)
|
||||
F2 = FV2[:, 0]*mkvc(n2x) + FV2[:, 1]*mkvc(n2y) + FV2[:, 2]*mkvc(n2z)
|
||||
F3 = FV3[:, 0]*mkvc(n3x) + FV3[:, 1]*mkvc(n3y) + FV3[:, 2]*mkvc(n3z)
|
||||
|
||||
return hstack((hstack((mkvc(F1), mkvc(F2))), mkvc(F3)))
|
||||
@@ -1,77 +0,0 @@
|
||||
from scipy import sparse
|
||||
from numpy import *
|
||||
|
||||
|
||||
def ddx(n):
|
||||
"""Define 1D derivatives"""
|
||||
# ddx = lambda n: sparse.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format='csr')
|
||||
return sparse.spdiags(-ones(n), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)
|
||||
|
||||
|
||||
def av(n):
|
||||
"""Define 1D average"""
|
||||
return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1))
|
||||
|
||||
|
||||
def sdiag(h):
|
||||
"""Diagonal matrix"""
|
||||
return sparse.spdiags(h, 0, size(h), size(h))
|
||||
|
||||
|
||||
def speye(n):
|
||||
"""sparse identity"""
|
||||
return sparse.spdiags(ones(n), 0, n, n)
|
||||
|
||||
|
||||
def kron3(A, B, C):
|
||||
"""two kron prods"""
|
||||
return sparse.kron(sparse.kron(A, B), C)
|
||||
|
||||
|
||||
def appendBottom(A, B):
|
||||
"""append on bottom"""
|
||||
C = sparse.vstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def appendBottom3(A, B, C):
|
||||
"""append on bottom"""
|
||||
C = appendBottom(appendBottom(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def appendRight(A, B):
|
||||
"""append on right"""
|
||||
C = sparse.hstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def appendRight3(A, B, C):
|
||||
"""append on right"""
|
||||
C = appendRight(appendRight(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def blkDiag(A, B):
|
||||
"""blockdigonal"""
|
||||
O12 = sparse.coo_matrix((shape(A)[0], shape(B)[1]))
|
||||
O21 = sparse.coo_matrix((shape(B)[0], shape(A)[1]))
|
||||
C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B))))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
|
||||
def blkDiag3(A, B, C):
|
||||
"""blockdigonal 3"""
|
||||
ABC = blkDiag(blkDiag(A, B), C)
|
||||
ABC = ABC.tocsr()
|
||||
return ABC
|
||||
|
||||
|
||||
def spzeros(n1, n2):
|
||||
"""spzeros"""
|
||||
return sparse.coo_matrix((n1, n2))
|
||||
@@ -1,222 +0,0 @@
|
||||
import numpy;
|
||||
import cmath;
|
||||
import math;
|
||||
|
||||
def prod(arg):
|
||||
""" returns the product of elements in arg.
|
||||
arg can be list, tuple, set, and array with numerical values. """
|
||||
ret = 1;
|
||||
for i in range(0,len(arg)):
|
||||
ret = ret * arg[i];
|
||||
return ret;
|
||||
|
||||
|
||||
def allIndices(dim):
|
||||
""" From the given shape of dimenions (e.g. (2,3,4)),
|
||||
generate a numpy.array of all, sorted indices."""
|
||||
|
||||
length = len(dim);
|
||||
|
||||
sub = numpy.arange(dim[length-1]).reshape(dim[length-1],1);
|
||||
|
||||
for d in range(length-2, -1, -1):
|
||||
for i in range(0, dim[d]):
|
||||
temp = numpy.ndarray([len(sub), 1]);
|
||||
temp.fill(i);
|
||||
temp = numpy.concatenate((temp,sub), axis=1);
|
||||
if(i == 0):
|
||||
newsub = temp;
|
||||
else:
|
||||
newsub = numpy.concatenate((newsub, temp), axis = 0);
|
||||
|
||||
sub = newsub;
|
||||
|
||||
return sub;
|
||||
|
||||
def find(nda, obj):
|
||||
"""returns the index of the obj in the given nda(ndarray, list, or tuple)"""
|
||||
for i in range(0, len(nda)):
|
||||
if(nda[i] == obj):
|
||||
return i;
|
||||
return -1;
|
||||
|
||||
|
||||
def notin(n, vector):
|
||||
"""returns a numpy.array object that contains
|
||||
elements in [0,1, ... n-1] but not in vector."""
|
||||
ret = numpy.arange(n).tolist();
|
||||
for i in vector:
|
||||
if (0 <= i and i < n):
|
||||
ret.remove(i);
|
||||
return numpy.array(ret);
|
||||
|
||||
|
||||
|
||||
def getelts(nda, indices):
|
||||
"""From the given nda(ndarray, list, or tuple), returns the list located at the given indices"""
|
||||
ret = [];
|
||||
for i in indices:
|
||||
ret.extend([nda[i]]);
|
||||
return numpy.array(ret);
|
||||
|
||||
def sub2ind(shape, subs):
|
||||
""" From the given shape, returns the index of the given subscript"""
|
||||
revshp = list(shape);
|
||||
revshp.reverse();
|
||||
mult = [1];
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]]);
|
||||
mult.reverse();
|
||||
mult = numpy.array(mult).reshape(len(mult),1);
|
||||
|
||||
idx = numpy.dot((subs) , (mult));
|
||||
return idx;
|
||||
|
||||
def ind2sub(shape, ind):
|
||||
""" From the given shape, returns the subscrips of the given index"""
|
||||
revshp = [];
|
||||
revshp.extend(shape);
|
||||
revshp.reverse();
|
||||
mult = [1];
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]]);
|
||||
mult.reverse();
|
||||
mult = numpy.array(mult).reshape(len(mult));
|
||||
|
||||
sub = [];
|
||||
|
||||
for i in range(0,len(shape)):
|
||||
sub.extend([math.floor(ind / mult[i])]);
|
||||
ind = ind - (math.floor(ind/mult[i]) * mult[i]);
|
||||
return sub;
|
||||
|
||||
def tt_dimscehck(dims, N, M = None, exceptdims = False):
|
||||
"""Checks whether the specified dimensions are valid in a tensor of N-dimension.
|
||||
If M is given, then it will also retuns an index for M multiplicands.
|
||||
If exceptdims == True, then it will compute for the dimensions not specified."""
|
||||
|
||||
# if exceptdims is true
|
||||
if(exceptdims):
|
||||
dims = listdiff(range(0,N), dims);
|
||||
|
||||
#check vals in between 0 and N-1
|
||||
for i in range(0, len(dims)):
|
||||
if(dims[i] < 0 or dims[i] >= N):
|
||||
raise ValueError("invalid dimensions specified");
|
||||
|
||||
# number of dimensions in dims
|
||||
p = len(dims);
|
||||
|
||||
sdims = [];
|
||||
sdims.extend(dims);
|
||||
sdims.sort();
|
||||
|
||||
#indices of the elements in the sorted array
|
||||
sidx = [];
|
||||
#table that denotes whether the index is used
|
||||
table = numpy.ndarray([len(sdims)]);
|
||||
table.fill(0);
|
||||
|
||||
for i in range(0, len(sdims)):
|
||||
for j in range(0, len(dims)):
|
||||
if(sdims[i] == dims[j] and table[j] == 0):
|
||||
sidx.extend([j]);
|
||||
table[j] = 1;
|
||||
break;
|
||||
|
||||
if (M == None):
|
||||
return sdims
|
||||
|
||||
if(M > N):
|
||||
raise ValueError("Cannot have more multiplicands than dimensions");
|
||||
|
||||
if(M != N and M != p):
|
||||
raise ValueError("invalid number of multiplicands");
|
||||
|
||||
if(M == p):
|
||||
vidx = sidx;
|
||||
else:
|
||||
vidx = sdims;
|
||||
|
||||
return (sdims, vidx);
|
||||
|
||||
def listtimes(list, c):
|
||||
"""multiplies the elements in the list by the given scalar value c"""
|
||||
ret = []
|
||||
for i in range(0, len(list)):
|
||||
ret.extend([list[i]]*c);
|
||||
return ret;
|
||||
|
||||
def listdiff(list1, list2):
|
||||
"""returns the list of elements that are in list 1 but not in list2"""
|
||||
if(list1.__class__ == numpy.ndarray):
|
||||
list1 = list1.tolist();
|
||||
if(list2.__class__ == numpy.ndarray):
|
||||
list2 = list2.tolist();
|
||||
ret = []
|
||||
for i in range(0,len(list1)):
|
||||
ok = true
|
||||
for j in range(0, len(list2)):
|
||||
if(list[i] == list[j]):
|
||||
ok = false;
|
||||
break;
|
||||
if(ok):
|
||||
ret.extend([list[i]]);
|
||||
return ret;
|
||||
|
||||
|
||||
|
||||
def tt_subscheck(subs):
|
||||
"""Check whether the given list of subscripts are valid. Used for sptensor"""
|
||||
isOk = True;
|
||||
if(subs.size == 0):
|
||||
isOk = True;
|
||||
|
||||
elif(subs.ndim != 2):
|
||||
isOk = False;
|
||||
|
||||
else:
|
||||
for i in range(0, (subs.size / subs[0].size)):
|
||||
for j in range(0, (subs[0].size)):
|
||||
val = subs[i][j];
|
||||
if( cmath.isnan(val) or cmath.isinf(val) or val < 0 or val != round(val) ):
|
||||
isOk = False;
|
||||
|
||||
if(not isOk):
|
||||
raise ValueError("Subscripts must be a matrix of non-negative integers");
|
||||
|
||||
return isOk;
|
||||
|
||||
|
||||
def tt_valscheck(vals):
|
||||
"""Check whether the given list of values are valid. Used for sptensor"""
|
||||
isOk = True;
|
||||
|
||||
if(vals.size == 0):
|
||||
isOk = True;
|
||||
|
||||
elif(vals.ndim != 2 or vals[0].size != 1):
|
||||
isOk = False;
|
||||
|
||||
if(not isOk):
|
||||
raise ValueError("values must be a column array");
|
||||
|
||||
return isOk;
|
||||
|
||||
def tt_sizecheck(size):
|
||||
"""Check whether the given size is valid. Used for sptensor"""
|
||||
size = numpy.array(size);
|
||||
isOk = True;
|
||||
|
||||
if(size.ndim != 1):
|
||||
isOk = False;
|
||||
else:
|
||||
for i in range(0, len(size)):
|
||||
val = size[i];
|
||||
if(cmath.isnan(val) or cmath.isinf(val)
|
||||
or val <= 0 or val != round(val)):
|
||||
isOk = False;
|
||||
|
||||
if(not isOk):
|
||||
raise ValueError("size must be a row vector of real positive integers");
|
||||
return isOk;
|
||||
@@ -1,87 +0,0 @@
|
||||
from numpy import *
|
||||
import numpy as np
|
||||
|
||||
|
||||
def diff(A, d):
|
||||
if(d == 1):
|
||||
return A[1:, :, :] - A[:-1, :, :]
|
||||
elif(d == 2):
|
||||
return A[:, 1:, :] - A[:, :-1, :]
|
||||
else:
|
||||
return A[:, :, 1:] - A[:, :, :-1]
|
||||
#else:
|
||||
# print('d must be 1,2 or 3')
|
||||
|
||||
|
||||
def diffp(A, d1, d2):
|
||||
if(d1 == 1 and d2 == 2):
|
||||
return A[1:, 1:, :] - A[:-1, :-1, :]
|
||||
elif(d1 == 1 and d2 == 3):
|
||||
return A[1:, :, 1:] - A[:-1, :, :-1]
|
||||
else:
|
||||
return A[:, 1:, 1:] - A[:, :-1, :-1]
|
||||
|
||||
|
||||
def diffm(A, d1, d2):
|
||||
if(d1 == 3 and d2 == 2):
|
||||
return A[:, :-1, 1:] - A[:, 1:, :-1]
|
||||
elif(d1 == 1 and d2 == 3):
|
||||
return A[1:, :, :-1] - A[:-1, :, 1:]
|
||||
elif(d1 == 2 and d2 == 1):
|
||||
return A[:-1, 1:, :] - A[1:, :-1, :]
|
||||
else:
|
||||
print('d must be 1, 2 or 3')
|
||||
|
||||
|
||||
def ave(A, d):
|
||||
if(d == 1):
|
||||
return 0.5*(A[1:, :, :] + A[:-1, :, :])
|
||||
elif(d == 2):
|
||||
return 0.5*(A[:, 1:, :] + A[:, :-1, :])
|
||||
elif(d == 3):
|
||||
return 0.5*(A[:, :, 1:] + A[:, :, :-1])
|
||||
else:
|
||||
print('d must be 1,2 or 3')
|
||||
|
||||
|
||||
def mkmat(x):
|
||||
return reshape(matrix(x), (size(x), 1), 'F')
|
||||
|
||||
|
||||
def hstack3(a, b, c):
|
||||
a = mkvc(a)
|
||||
b = mkvc(b)
|
||||
c = mkvc(c)
|
||||
a = mkmat(a)
|
||||
b = mkmat(b)
|
||||
c = mkmat(c)
|
||||
return hstack((hstack((a, b)), c))
|
||||
|
||||
|
||||
def ind2sub(shape, ind):
|
||||
"""From the given shape, returns the subscrips of the given index"""
|
||||
revshp = []
|
||||
revshp.extend(shape)
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult))
|
||||
|
||||
sub = []
|
||||
|
||||
for i in range(0, len(shape)):
|
||||
sub.extend([math.floor(ind / mult[i])])
|
||||
ind = ind - (math.floor(ind/mult[i]) * mult[i])
|
||||
return sub
|
||||
|
||||
|
||||
def sub2ind(shape, subs):
|
||||
"""From the given shape, returns the index of the given subscript"""
|
||||
revshp = list(shape)
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult), 1)
|
||||
|
||||
idx = dot((subs), (mult))
|
||||
return idx
|
||||
@@ -1,119 +0,0 @@
|
||||
#============= Nodal Gradients ===========================
|
||||
def getNodalGradient(h1,h2,h3):
|
||||
|
||||
n1 = size(h1)
|
||||
n2 = size(h2)
|
||||
n3 = size(h3)
|
||||
D1 = kron3(speye(n3+1),speye(n2+1),ddx(n1))
|
||||
D2 = kron3(speye(n3+1),ddx(n2),speye(n1+1))
|
||||
D3 = kron3(ddx(n3),speye(n2+1),speye(n1+1))
|
||||
|
||||
# topological gradient
|
||||
GRAD = appendBottom3(D1,D2,D3)
|
||||
|
||||
# scale for non-uniform mesh
|
||||
L = blkDiag3(kron3(speye(n3+1),speye(n2+1),sdiag(1/h1)),
|
||||
kron3(speye(n3+1),sdiag(1/h2),speye(n1+1)),
|
||||
kron3(sdiag(1/h3),speye(n2+1),speye(n1+1)))
|
||||
|
||||
return L*GRAD
|
||||
|
||||
#============= Edge CURL ===========================
|
||||
def getCurlMatrix(h1,h2,h3):
|
||||
|
||||
n1 = size(h1)
|
||||
n2 = size(h2)
|
||||
n3 = size(h3)
|
||||
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
# derivatives on x-edge variables
|
||||
D32 = kron3(d3,speye(n2),speye(n1+1))
|
||||
D23 = kron3(speye(n3),d2,speye(n1+1))
|
||||
D31 = kron3(d3,speye(n2+1),speye(n1))
|
||||
D13 = kron3(speye(n3),speye(n2+1),d1)
|
||||
D21 = kron3(speye(n3+1),d2,speye(n1))
|
||||
D12 = kron3(speye(n3+1),speye(n2),d1)
|
||||
|
||||
O1 = spzeros(shape(D32)[0],shape(D31)[1])
|
||||
O2 = spzeros(shape(D31)[0],shape(D32)[1])
|
||||
O3 = spzeros(shape(D21)[0],shape(D13)[1])
|
||||
|
||||
CURL = appendBottom3(
|
||||
appendRight3(O1, -D32, D23),
|
||||
appendRight3(D31, O2, -D13),
|
||||
appendRight3(-D21, D12, O3))
|
||||
|
||||
# scale for non-uniform mesh
|
||||
F = blkDiag3(kron3(sdiag(1/h3),sdiag(1/h2),speye(n1+1)),
|
||||
kron3(sdiag(1/h3),speye(n2+1),sdiag(1/h1)),
|
||||
kron3(speye(n3+1),sdiag(1/h2),sdiag(1/h1)))
|
||||
|
||||
L = blkDiag3(kron3(speye(n3+1),speye(n2+1),sdiag(h1)),
|
||||
kron3(speye(n3+1),sdiag(h2),speye(n1+1)),
|
||||
kron3(sdiag(h3),speye(n2+1),speye(n1+1)))
|
||||
|
||||
|
||||
return F*(CURL*L)
|
||||
|
||||
#============= Face DIV ===========================
|
||||
def getDivMatrix(h1,h2,h3):
|
||||
|
||||
n1 = size(h1)
|
||||
n2 = size(h2)
|
||||
n3 = size(h3)
|
||||
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
D1 = kron3(speye(n3),speye(n2),d1)
|
||||
D2 = kron3(speye(n3),d2,speye(n1))
|
||||
D3 = kron3(d3,speye(n2),speye(n1))
|
||||
|
||||
# divergence on faces
|
||||
D = appendRight3(D1, D2, D3)
|
||||
|
||||
# scale for non-uniform mesh
|
||||
F = blkDiag3(kron3(sdiag(h3),sdiag(h2),speye(n1+1)),
|
||||
kron3(sdiag(h3),speye(n2+1),sdiag(h1)),
|
||||
kron3(speye(n3+1),sdiag(h2),sdiag(h1)))
|
||||
|
||||
V = kron3(sdiag(1/h3),sdiag(1/h2),sdiag(1/h1))
|
||||
|
||||
return V*(D*F)
|
||||
|
||||
#====== Face Averageing =================
|
||||
def getFaceAverage(n1,n2,n3):
|
||||
|
||||
av1 = av(n1)
|
||||
av2 = av(n2)
|
||||
av3 = av(n3)
|
||||
|
||||
Af = appendRight3(kron3(speye(n3),speye(n2),av1),
|
||||
kron3(speye(n3),av2,speye(n1)),
|
||||
kron3(av3,speye(n2),speye(n1)))
|
||||
return Af
|
||||
|
||||
#====== Edge Averageing =================
|
||||
def getEdgeAverage(n1,n2,n3):
|
||||
|
||||
av1 = av(n1)
|
||||
av2 = av(n2)
|
||||
av3 = av(n3)
|
||||
|
||||
Ae = appendRight3(kron3(av3,av2,speye(n1)),
|
||||
kron3(av3,speye(n2),av1),
|
||||
kron3(speye(n3),av2,av1))
|
||||
return Ae
|
||||
|
||||
#====== Node Averageing =================
|
||||
def getNodeAverage(n1,n2,n3):
|
||||
|
||||
av1 = av(n1)
|
||||
av2 = av(n2)
|
||||
av3 = av(n3)
|
||||
|
||||
return kron3(av3,av2,av1)
|
||||
|
||||
|
||||
@@ -0,0 +1,336 @@
|
||||
from scipy import sparse as sp
|
||||
from sputils import sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
|
||||
from utils import sub2ind, ndgrid, mkvc, getSubArray
|
||||
import numpy as np
|
||||
|
||||
|
||||
class InnerProducts(object):
|
||||
"""
|
||||
Class creates the inner product matrices that you need!
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
|
||||
def getFaceInnerProduct(self, mu=None, returnP=False):
|
||||
if self.dim == 2:
|
||||
return getFaceInnerProduct2D(self, mu, returnP)
|
||||
elif self.dim == 3:
|
||||
return getFaceInnerProduct(self, mu, returnP)
|
||||
|
||||
def getEdgeInnerProduct(self, sigma=None, returnP=False):
|
||||
if self.dim == 2:
|
||||
return getEdgeInnerProduct2D(self, sigma, returnP)
|
||||
elif self.dim == 3:
|
||||
return getEdgeInnerProduct(self, sigma, returnP)
|
||||
|
||||
# ------------------------ Geometries ------------------------------
|
||||
#
|
||||
#
|
||||
# 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)
|
||||
|
||||
|
||||
def getFaceInnerProduct(mesh, mu=None, returnP=False):
|
||||
|
||||
if mu is None: # default is ones
|
||||
mu = np.ones((mesh.nC, 1))
|
||||
|
||||
m = np.array([mesh.nCx, mesh.nCy, mesh.nCz])
|
||||
nc = mesh.nC
|
||||
|
||||
i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2]))
|
||||
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M')
|
||||
fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M')
|
||||
fN3 = mesh.r(mesh.normals, 'F', 'Fz', 'M')
|
||||
|
||||
def Pxxx(pos):
|
||||
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
|
||||
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0]
|
||||
ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1]
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr()
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]),
|
||||
getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]),
|
||||
getSubArray(fN3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]))
|
||||
PXXX = I3x3 * PXXX
|
||||
|
||||
return PXXX
|
||||
|
||||
# no | node | f1 | f2 | f3
|
||||
# 000 | i ,j ,k | i , j, k | i, j , k | i, j, k
|
||||
# 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k
|
||||
# 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k
|
||||
# 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k
|
||||
# 001 | i ,j ,k | i , j, k | i, j , k | i, j, k+1
|
||||
# 101 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k+1
|
||||
# 011 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k+1
|
||||
# 111 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k+1
|
||||
|
||||
# Square root of cell volume multiplied by 1/8
|
||||
v = np.sqrt(0.125*mesh.vol)
|
||||
V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry
|
||||
|
||||
P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P100 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P010 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
|
||||
P110 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
|
||||
P001 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]])
|
||||
P101 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]])
|
||||
P011 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]])
|
||||
P111 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
|
||||
|
||||
if mu.size == mesh.nC: # Isotropic!
|
||||
mu = mkvc(mu) # ensure it is a vector.
|
||||
Mu = sdiag(np.r_[mu, mu, mu])
|
||||
elif mu.shape[1] == 3: # Diagonal tensor
|
||||
Mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]])
|
||||
elif mu.shape[1] == 6: # Fully anisotropic
|
||||
row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4])))
|
||||
row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5])))
|
||||
row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2])))
|
||||
Mu = sp.vstack((row1, row2, row3))
|
||||
|
||||
A = P000.T*Mu*P000 + P001.T*Mu*P001 + P010.T*Mu*P010 + P011.T*Mu*P011 + P100.T*Mu*P100 + P101.T*Mu*P101 + P110.T*Mu*P110 + P111.T*Mu*P111
|
||||
P = [P000, P001, P010, P011, P100, P101, P110, P111]
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
|
||||
|
||||
def getFaceInnerProduct2D(mesh, mu=None, returnP=False):
|
||||
|
||||
if mu is None: # default is ones
|
||||
mu = np.ones((mesh.nC, 1))
|
||||
|
||||
m = np.array([mesh.nCx, mesh.nCy])
|
||||
nc = mesh.nC
|
||||
|
||||
i, j = np.int64(range(m[0])), np.int64(range(m[1]))
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M')
|
||||
fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M')
|
||||
|
||||
def Pxx(pos):
|
||||
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]])
|
||||
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0]
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nF))).tocsr()
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1]]),
|
||||
getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1]]))
|
||||
PXX = I2x2 * PXX
|
||||
|
||||
return PXX
|
||||
|
||||
# no | node | f1 | f2
|
||||
# 00 | i ,j | i , j | i, j
|
||||
# 10 | i+1,j | i+1, j | i, j
|
||||
# 01 | i ,j+1 | i , j | i, j+1
|
||||
# 11 | i+1,j+1 | i+1, j | i, j+1
|
||||
|
||||
# Square root of cell volume multiplied by 1/4
|
||||
v = np.sqrt(0.25*mesh.vol)
|
||||
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
|
||||
|
||||
P00 = V2*Pxx([[0, 0], [0, 0]])
|
||||
P10 = V2*Pxx([[1, 0], [0, 0]])
|
||||
P01 = V2*Pxx([[0, 0], [0, 1]])
|
||||
P11 = V2*Pxx([[1, 0], [0, 1]])
|
||||
|
||||
if mu.size == mesh.nC: # Isotropic!
|
||||
mu = mkvc(mu) # ensure it is a vector.
|
||||
Mu = sdiag(np.r_[mu, mu])
|
||||
elif mu.shape[1] == 2: # Diagonal tensor
|
||||
Mu = sdiag(np.r_[mu[:, 0], mu[:, 1]])
|
||||
elif mu.shape[1] == 3: # Fully anisotropic
|
||||
row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 2])))
|
||||
row2 = sp.hstack((sdiag(mu[:, 2]), sdiag(mu[:, 1])))
|
||||
Mu = sp.vstack((row1, row2))
|
||||
|
||||
A = P00.T*Mu*P00 + P10.T*Mu*P10 + P01.T*Mu*P01 + P11.T*Mu*P11
|
||||
P = [P00, P10, P01, P11]
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
|
||||
|
||||
def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
|
||||
|
||||
if sigma is None: # default is ones
|
||||
sigma = np.ones((mesh.nC, 1))
|
||||
|
||||
m = np.array([mesh.nCx, mesh.nCy, mesh.nCz])
|
||||
nc = mesh.nC
|
||||
|
||||
i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2]))
|
||||
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
|
||||
eT3 = mesh.r(mesh.tangents, 'E', 'Ez', 'M')
|
||||
|
||||
def Pxxx(pos):
|
||||
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
|
||||
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0]
|
||||
ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1]
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]),
|
||||
getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]),
|
||||
getSubArray(eT3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]))
|
||||
PXXX = I3x3 * PXXX
|
||||
|
||||
return PXXX
|
||||
|
||||
# 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
|
||||
|
||||
# Square root of cell volume multiplied by 1/8
|
||||
v = np.sqrt(0.125*mesh.vol)
|
||||
V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry
|
||||
|
||||
P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P100 = V3*Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]])
|
||||
P010 = V3*Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]])
|
||||
P110 = V3*Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]])
|
||||
P001 = V3*Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]])
|
||||
P101 = V3*Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]])
|
||||
P011 = V3*Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]])
|
||||
P111 = V3*Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
|
||||
|
||||
if sigma.size == mesh.nC: # Isotropic!
|
||||
sigma = mkvc(sigma) # ensure it is a vector.
|
||||
Sigma = sdiag(np.r_[sigma, sigma, sigma])
|
||||
elif sigma.shape[1] == 3: # Diagonal tensor
|
||||
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
|
||||
elif sigma.shape[1] == 6: # Fully anisotropic
|
||||
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
|
||||
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
|
||||
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
|
||||
Sigma = sp.vstack((row1, row2, row3))
|
||||
|
||||
A = P000.T*Sigma*P000 + P001.T*Sigma*P001 + P010.T*Sigma*P010 + P011.T*Sigma*P011 + P100.T*Sigma*P100 + P101.T*Sigma*P101 + P110.T*Sigma*P110 + P111.T*Sigma*P111
|
||||
P = [P000, P001, P010, P011, P100, P101, P110, P111]
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
|
||||
|
||||
def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False):
|
||||
|
||||
if sigma is None: # default is ones
|
||||
sigma = np.ones((mesh.nC, 1))
|
||||
|
||||
m = np.array([mesh.nCx, mesh.nCy])
|
||||
nc = mesh.nC
|
||||
|
||||
i, j = np.int64(range(m[0])), np.int64(range(m[1]))
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
|
||||
|
||||
def Pxx(pos):
|
||||
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]])
|
||||
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0]
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nE))).tocsr()
|
||||
|
||||
if mesh._meshType == 'LOM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1]]),
|
||||
getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1]]))
|
||||
PXX = I2x2 * PXX
|
||||
|
||||
return PXX
|
||||
|
||||
# no | node | e1 | e2
|
||||
# 00 | i ,j | i ,j | i ,j
|
||||
# 10 | i+1,j | i ,j | i+1,j
|
||||
# 01 | i ,j+1 | i ,j+1 | i ,j
|
||||
# 11 | i+1,j+1 | i ,j+1 | i+1,j
|
||||
|
||||
# Square root of cell volume multiplied by 1/4
|
||||
v = np.sqrt(0.25*mesh.vol)
|
||||
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
|
||||
|
||||
P00 = V2*Pxx([[0, 0], [0, 0]])
|
||||
P10 = V2*Pxx([[0, 0], [1, 0]])
|
||||
P01 = V2*Pxx([[0, 1], [0, 0]])
|
||||
P11 = V2*Pxx([[0, 1], [1, 0]])
|
||||
|
||||
if sigma.size == mesh.nC: # Isotropic!
|
||||
sigma = mkvc(sigma) # ensure it is a vector.
|
||||
Sigma = sdiag(np.r_[sigma, sigma])
|
||||
elif sigma.shape[1] == 2: # Diagonal tensor
|
||||
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
|
||||
elif sigma.shape[1] == 3: # Fully anisotropic
|
||||
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
|
||||
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
|
||||
Sigma = sp.vstack((row1, row2))
|
||||
|
||||
A = P00.T*Sigma*P00 + P10.T*Sigma*P10 + P01.T*Sigma*P01 + P11.T*Sigma*P11
|
||||
P = [P00, P10, P01, P11]
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from TensorMesh import TensorMesh
|
||||
h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])]
|
||||
mesh = TensorMesh(h)
|
||||
mu = np.ones((mesh.nC, 6))
|
||||
A, P = mesh.getFaceInnerProduct(mu, returnP=True)
|
||||
B, P = mesh.getEdgeInnerProduct(mu, returnP=True)
|
||||
@@ -0,0 +1,336 @@
|
||||
import numpy as np
|
||||
from BaseMesh import BaseMesh
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
from LomView import LomView
|
||||
from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
|
||||
|
||||
# Some helper functions.
|
||||
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
|
||||
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
|
||||
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), mkvc(length2D(x), 2))
|
||||
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), mkvc(length3D(x), 2))
|
||||
|
||||
|
||||
class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
|
||||
"""
|
||||
LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes.
|
||||
|
||||
"""
|
||||
_meshType = 'LOM'
|
||||
|
||||
def __init__(self, nodes):
|
||||
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
|
||||
|
||||
for i, nodes_i in enumerate(nodes):
|
||||
assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i)
|
||||
assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
|
||||
|
||||
assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
|
||||
assert len(nodes[0].shape) > 1, "Not worth using LOM for a 1D mesh."
|
||||
|
||||
super(LogicallyOrthogonalMesh, self).__init__(np.array(nodes[0].shape)-1, None)
|
||||
|
||||
# Save nodes to private variable _gridN as vectors
|
||||
self._gridN = np.ones((nodes[0].size, self.dim))
|
||||
for i, node_i in enumerate(nodes):
|
||||
self._gridN[:, i] = mkvc(node_i)
|
||||
|
||||
def gridCC():
|
||||
doc = "Cell-centered grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridCC is None:
|
||||
ccV = (self.nodalVectorAve*mkvc(self.gridN))
|
||||
self._gridCC = ccV.reshape((-1, self.dim), order='F')
|
||||
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:
|
||||
raise Exception("Someone deleted this. I blame you.")
|
||||
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:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
|
||||
self._gridFx = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
|
||||
self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
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:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
|
||||
self._gridFy = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
|
||||
self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
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 and self.dim == 3:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
XYZ = [mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
|
||||
self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
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:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
|
||||
self._gridEx = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
|
||||
self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
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:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
|
||||
self._gridEy = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
|
||||
self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
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 and self.dim == 3:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
XYZ = [mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
|
||||
self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridEz
|
||||
return locals()
|
||||
_gridEz = None # Store grid by default
|
||||
gridEz = property(**gridEz())
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
#
|
||||
#
|
||||
# ------------------- 2D -------------------------
|
||||
#
|
||||
# node(i,j) node(i,j+1)
|
||||
# A -------------- B
|
||||
# | |
|
||||
# | cell(i,j) |
|
||||
# | I |
|
||||
# | |
|
||||
# D -------------- C
|
||||
# node(i+1,j) node(i+1,j+1)
|
||||
#
|
||||
# ------------------- 3D -------------------------
|
||||
#
|
||||
#
|
||||
# node(i,j,k+1) node(i,j+1,k+1)
|
||||
# E --------------- F
|
||||
# /| / |
|
||||
# / | / |
|
||||
# / | / |
|
||||
# node(i,j,k) node(i,j+1,k)
|
||||
# A -------------- B |
|
||||
# | H ----------|---- G
|
||||
# | /cell(i,j) | /
|
||||
# | / I | /
|
||||
# | / | /
|
||||
# D -------------- C
|
||||
# node(i+1,j,k) node(i+1,j+1,k)
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._vol is None):
|
||||
if self.dim == 2:
|
||||
A, B, C, D = indexCube('ABCD', self.n+1)
|
||||
normal, area = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
|
||||
self._vol = area
|
||||
elif self.dim == 3:
|
||||
# Each polyhedron can be decomposed into 5 tetrahedrons
|
||||
# However, this presents a choice so we may as well divide in two ways and average.
|
||||
A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', self.n+1)
|
||||
|
||||
vol1 = (volTetra(self.gridN, A, B, D, E) + # cutted edge top
|
||||
volTetra(self.gridN, B, E, F, G) + # cutted edge top
|
||||
volTetra(self.gridN, B, D, E, G) + # middle
|
||||
volTetra(self.gridN, B, C, D, G) + # cutted edge bottom
|
||||
volTetra(self.gridN, D, E, G, H)) # cutted edge bottom
|
||||
|
||||
vol2 = (volTetra(self.gridN, A, F, B, C) + # cutted edge top
|
||||
volTetra(self.gridN, A, E, F, H) + # cutted edge top
|
||||
volTetra(self.gridN, A, H, F, C) + # middle
|
||||
volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
|
||||
volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
|
||||
|
||||
self._vol = (vol1 + vol2)/2
|
||||
return self._vol
|
||||
return locals()
|
||||
_vol = None
|
||||
vol = property(**vol())
|
||||
|
||||
def area():
|
||||
doc = "Face areas."
|
||||
|
||||
def fget(self):
|
||||
if(self._area is None or self._normals is None):
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
|
||||
edge1 = xy[B, :] - xy[A, :]
|
||||
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
|
||||
area1 = length2D(edge1)
|
||||
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
|
||||
# Note that we are doing A-D to make sure the normal points the right way.
|
||||
# Think about it. Look at the picture. Normal points towards C iff you do this.
|
||||
edge2 = xy[A, :] - xy[D, :]
|
||||
normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
|
||||
area2 = length2D(edge2)
|
||||
self._area = np.r_[mkvc(area1), mkvc(area2)]
|
||||
self._normals = [normalize2D(normal1), normalize2D(normal2)]
|
||||
elif(self.dim == 3):
|
||||
|
||||
A, E, F, B = indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
|
||||
normal1, area1 = faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
|
||||
|
||||
A, D, H, E = indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
|
||||
normal2, area2 = faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
|
||||
|
||||
A, B, C, D = indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
|
||||
normal3, area3 = faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
|
||||
|
||||
self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)]
|
||||
self._normals = [normal1, normal2, normal3]
|
||||
return self._area
|
||||
return locals()
|
||||
_area = None
|
||||
area = property(**area())
|
||||
|
||||
def normals():
|
||||
doc = """
|
||||
Face normals: calling this will average
|
||||
the computed normals so that there is one
|
||||
per face. This is especially relevant in
|
||||
3D, as there are up to 4 different normals
|
||||
for each face that will be different.
|
||||
|
||||
To reshape the normals into a matrix and get the y component:
|
||||
|
||||
NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
"""
|
||||
|
||||
def fget(self):
|
||||
if(self._normals is None):
|
||||
self.area # calling .area will create the face normals
|
||||
if self.dim == 2:
|
||||
return normalize2D(np.r_[self._normals[0], self._normals[1]])
|
||||
elif self.dim == 3:
|
||||
normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4
|
||||
normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
|
||||
normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
|
||||
return normalize3D(np.r_[normal1, normal2, normal3])
|
||||
return locals()
|
||||
_normals = None
|
||||
normals = property(**normals())
|
||||
|
||||
def edge():
|
||||
doc = "Edge legnths."
|
||||
|
||||
def fget(self):
|
||||
if(self._edge is None or self._tangents is None):
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
|
||||
edge1 = xy[D, :] - xy[A, :]
|
||||
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
|
||||
edge2 = xy[B, :] - xy[A, :]
|
||||
self._edge = np.r_[mkvc(length2D(edge1)), mkvc(length2D(edge2))]
|
||||
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
|
||||
elif(self.dim == 3):
|
||||
xyz = self.gridN
|
||||
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
|
||||
edge1 = xyz[D, :] - xyz[A, :]
|
||||
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
|
||||
edge2 = xyz[B, :] - xyz[A, :]
|
||||
A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
|
||||
edge3 = xyz[E, :] - xyz[A, :]
|
||||
self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), mkvc(length3D(edge3))]
|
||||
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
|
||||
return self._edge
|
||||
return locals()
|
||||
_edge = None
|
||||
edge = property(**edge())
|
||||
|
||||
def tangents():
|
||||
doc = "Edge tangents."
|
||||
|
||||
def fget(self):
|
||||
if(self._tangents is None):
|
||||
self.edge # calling .edge will create the tangents
|
||||
return self._tangents
|
||||
return locals()
|
||||
_tangents = None
|
||||
tangents = property(**tangents())
|
||||
|
||||
if __name__ == '__main__':
|
||||
nc = 5
|
||||
h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
|
||||
nc = 7
|
||||
h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
|
||||
h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
|
||||
dee3 = True
|
||||
if dee3:
|
||||
X, Y, Z = ndgrid(h1, h2, h3, vector=False)
|
||||
M = LogicallyOrthogonalMesh([X, Y, Z])
|
||||
else:
|
||||
X, Y = ndgrid(h1, h2, vector=False)
|
||||
M = LogicallyOrthogonalMesh([X, Y])
|
||||
|
||||
print M.r(M.normals, 'F', 'Fx', 'V')
|
||||
@@ -0,0 +1,91 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
from utils import mkvc
|
||||
|
||||
|
||||
class LomView(object):
|
||||
"""
|
||||
Provides viewing functions for TensorMesh
|
||||
|
||||
This class is inherited by TensorMesh
|
||||
"""
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def plotGrid(self, length=0.05):
|
||||
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions."""
|
||||
NN = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
fig = plt.figure(2)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2]
|
||||
Y = np.r_[Y1, Y2]
|
||||
|
||||
plt.plot(X, Y)
|
||||
|
||||
plt.hold(True)
|
||||
Nx = self.r(self.normals, 'F', 'Fx', 'V')
|
||||
Ny = self.r(self.normals, 'F', 'Fy', 'V')
|
||||
Tx = self.r(self.tangents, 'E', 'Ex', 'V')
|
||||
Ty = self.r(self.tangents, 'E', 'Ey', 'V')
|
||||
|
||||
plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
|
||||
|
||||
nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
plt.plot(nX, nY, 'r-')
|
||||
|
||||
nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
|
||||
tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
|
||||
tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
|
||||
plt.plot(tX, tY, 'r-')
|
||||
|
||||
nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
plt.axis('equal')
|
||||
|
||||
elif self.dim == 3:
|
||||
fig = plt.figure(3)
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
|
||||
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
|
||||
Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
|
||||
|
||||
X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
|
||||
Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
|
||||
Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2, X3]
|
||||
Y = np.r_[Y1, Y2, Y3]
|
||||
Z = np.r_[Z1, Z2, Z3]
|
||||
|
||||
plt.plot(X, Y, 'b', zs=Z)
|
||||
ax.set_zlabel('x3')
|
||||
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
fig.show()
|
||||
+8
-12
@@ -2,10 +2,11 @@ import numpy as np
|
||||
from BaseMesh import BaseMesh
|
||||
from TensorView import TensorView
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
from utils import ndgrid, mkvc
|
||||
|
||||
|
||||
class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
|
||||
@@ -21,6 +22,8 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
mesh = TensorMesh([hx, hy, hz])
|
||||
|
||||
"""
|
||||
_meshType = 'TENSOR'
|
||||
|
||||
def __init__(self, h, x0=None):
|
||||
super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0)
|
||||
|
||||
@@ -181,7 +184,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
doc = "Face staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFy is None:
|
||||
if self._gridFy is None and self.dim > 1:
|
||||
self._gridFy = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorCCz] if not x is None])
|
||||
return self._gridFy
|
||||
return locals()
|
||||
@@ -192,7 +195,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
doc = "Face staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFz is None:
|
||||
if self._gridFz is None and self.dim > 2:
|
||||
self._gridFz = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorNz] if not x is None])
|
||||
return self._gridFz
|
||||
return locals()
|
||||
@@ -214,7 +217,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
doc = "Edge staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEy is None:
|
||||
if self._gridEy is None and self.dim > 1:
|
||||
self._gridEy = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorNz] if not x is None])
|
||||
return self._gridEy
|
||||
return locals()
|
||||
@@ -225,20 +228,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
doc = "Edge staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEz is None:
|
||||
if self._gridEz is None and self.dim > 2:
|
||||
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
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
|
||||
@@ -1,2 +1,4 @@
|
||||
from TensorMesh import TensorMesh
|
||||
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
|
||||
import utils
|
||||
from exampleGrid import exampleLomGird
|
||||
|
||||
@@ -0,0 +1,25 @@
|
||||
import numpy as np
|
||||
from utils import ndgrid
|
||||
|
||||
|
||||
def exampleLomGird(nC, exType):
|
||||
assert type(nC) == list, "nC must be a list containing the number of nodes"
|
||||
assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions"
|
||||
exType = exType.lower()
|
||||
|
||||
possibleTypes = ['rect', 'rotate']
|
||||
assert exType in possibleTypes, "Not a possible example type."
|
||||
|
||||
if exType == 'rect':
|
||||
return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
|
||||
elif exType == 'rotate':
|
||||
if len(nC) == 2:
|
||||
X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
|
||||
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
|
||||
amt[amt < 0] = 0
|
||||
return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
|
||||
elif len(nC) == 3:
|
||||
X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
|
||||
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
|
||||
amt[amt < 0] = 0
|
||||
return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt
|
||||
+66
-34
@@ -1,10 +1,10 @@
|
||||
import numpy as np
|
||||
from scipy import sparse as sp
|
||||
from utils import mkvc
|
||||
|
||||
|
||||
def sdiag(h):
|
||||
"""Sparse diagonal matrix"""
|
||||
return sp.spdiags(h, 0, np.size(h), np.size(h), format="csr")
|
||||
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
|
||||
|
||||
|
||||
def speye(n):
|
||||
@@ -22,45 +22,77 @@ def spzeros(n1, n2):
|
||||
return sp.coo_matrix((n1, n2)).tocsr()
|
||||
|
||||
|
||||
def appendBottom(A, B):
|
||||
"""append on bottom"""
|
||||
C = sp.vstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
|
||||
""" B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)
|
||||
|
||||
inverts a stack of 3x3 matrices
|
||||
|
||||
Input:
|
||||
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
|
||||
|
||||
Output:
|
||||
B - inverse
|
||||
"""
|
||||
|
||||
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 = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))),
|
||||
sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))),
|
||||
sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33)))))
|
||||
|
||||
return B
|
||||
|
||||
|
||||
def appendBottom3(A, B, C):
|
||||
"""append on bottom"""
|
||||
C = appendBottom(appendBottom(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
def inv2X2BlockDiagonal(a11, a12, a21, a22):
|
||||
""" B = inv2X2BlockDiagonal(a11, a12, a21, a22)
|
||||
|
||||
Inverts a stack of 2x2 matrices by using the inversion formula
|
||||
|
||||
def appendRight(A, B):
|
||||
"""append on right"""
|
||||
C = sp.hstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
inv(A) = (1/det(A)) * cof(A)^T
|
||||
|
||||
Input:
|
||||
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
|
||||
|
||||
def appendRight3(A, B, C):
|
||||
"""append on right"""
|
||||
C = appendRight(appendRight(A, B), C)
|
||||
C = C.tocsr()
|
||||
return C
|
||||
Output:
|
||||
B - inverse
|
||||
"""
|
||||
|
||||
a11 = mkvc(a11)
|
||||
a12 = mkvc(a12)
|
||||
a21 = mkvc(a21)
|
||||
a22 = mkvc(a22)
|
||||
|
||||
def blkDiag(A, B):
|
||||
"""blockdigonal"""
|
||||
O12 = sp.coo_matrix((np.shape(A)[0], np.shape(B)[1]))
|
||||
O21 = sp.coo_matrix((np.shape(B)[0], np.shape(A)[1]))
|
||||
C = sp.vstack((sp.hstack((A, O12)), sp.hstack((O21, B))))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
# compute inverse of the determinant.
|
||||
detAinv = 1./(a11*a22 - a21*a12)
|
||||
|
||||
b11 = +detAinv*a22
|
||||
b12 = -detAinv*a12
|
||||
b21 = -detAinv*a21
|
||||
b22 = +detAinv*a11
|
||||
|
||||
def blkDiag3(A, B, C):
|
||||
"""blockdigonal 3"""
|
||||
ABC = blkDiag(blkDiag(A, B), C)
|
||||
ABC = ABC.tocsr()
|
||||
return ABC
|
||||
B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
|
||||
sp.hstack((sdiag(b21), sdiag(b22)))))
|
||||
|
||||
return B
|
||||
|
||||
+78
-28
@@ -1,6 +1,6 @@
|
||||
import sys
|
||||
sys.path.append('../../')
|
||||
from SimPEG import TensorMesh
|
||||
from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh, exampleLomGird
|
||||
import numpy as np
|
||||
import unittest
|
||||
|
||||
@@ -20,7 +20,7 @@ class OrderTest(unittest.TestCase):
|
||||
|
||||
Note that you can provide any norm.
|
||||
|
||||
Test is passed when estimated rate order of convergence is at least 90% of the
|
||||
Test is passed when estimated rate order of convergence is at least within the specified tolerance of the
|
||||
estimated rate supplied by the user.
|
||||
|
||||
Minimal example for a curl operator:
|
||||
@@ -62,18 +62,50 @@ class OrderTest(unittest.TestCase):
|
||||
"""
|
||||
|
||||
name = "Order Test"
|
||||
expectedOrder = 2
|
||||
expectedOrders = 2. # This can be a list of orders, must be the same length as meshTypes
|
||||
_expectedOrder = 2.
|
||||
tolerance = 0.85
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
_meshType = meshTypes[0]
|
||||
meshDimension = 3
|
||||
|
||||
def setupMesh(self, nc):
|
||||
"""
|
||||
For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
|
||||
"""
|
||||
h1 = np.ones(nc)/nc
|
||||
h2 = np.ones(nc)/nc
|
||||
h3 = np.ones(nc)/nc
|
||||
h = [h1, h2, h3]
|
||||
self.M = TensorMesh(h)
|
||||
if 'TensorMesh' in self._meshType:
|
||||
if 'uniform' in self._meshType:
|
||||
h1 = np.ones(nc)/nc
|
||||
h2 = np.ones(nc)/nc
|
||||
h3 = np.ones(nc)/nc
|
||||
h = [h1, h2, h3]
|
||||
elif 'random' in self._meshType:
|
||||
h1 = np.random.rand(nc)
|
||||
h2 = np.random.rand(nc)
|
||||
h3 = np.random.rand(nc)
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
|
||||
self.M = TensorMesh(h[:self.meshDimension])
|
||||
max_h = max([np.max(hi) for hi in self.M.h])
|
||||
return max_h
|
||||
|
||||
elif 'LOM' in self._meshType:
|
||||
if 'uniform' in self._meshType:
|
||||
kwrd = 'rect'
|
||||
elif 'rotate' in self._meshType:
|
||||
kwrd = 'rotate'
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
if self.meshDimension == 2:
|
||||
X, Y = exampleLomGird([nc, nc], kwrd)
|
||||
self.M = LogicallyOrthogonalMesh([X, Y])
|
||||
if self.meshDimension == 3:
|
||||
X, Y, Z = exampleLomGird([nc, nc, nc], kwrd)
|
||||
self.M = LogicallyOrthogonalMesh([X, Y, Z])
|
||||
return 1./nc
|
||||
|
||||
def getError(self):
|
||||
"""For given h, generate A[h], f and A(f) and return norm of error."""
|
||||
@@ -87,27 +119,45 @@ class OrderTest(unittest.TestCase):
|
||||
|
||||
|
||||
"""
|
||||
order = []
|
||||
err_old = 0.
|
||||
nc_old = 0.
|
||||
for ii, nc in enumerate(self.meshSizes):
|
||||
self.setupMesh(nc)
|
||||
err = self.getError()
|
||||
if ii == 0:
|
||||
print ''
|
||||
print 'Testing order of: ' + self.name
|
||||
print '_____________________________________________'
|
||||
print ' h | error | e(i-1)/e(i) | order'
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
|
||||
print '%4i | %8.2e |' % (nc, err)
|
||||
else:
|
||||
order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc)))
|
||||
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
|
||||
err_old = err
|
||||
nc_old = nc
|
||||
print '---------------------------------------------'
|
||||
self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order)))
|
||||
assert type(self.meshTypes) == list, 'meshTypes must be a list'
|
||||
|
||||
# if we just provide one expected order, repeat it for each mesh type
|
||||
if type(self.expectedOrders) == float or type(self.expectedOrders) == int:
|
||||
self.expectedOrders = [self.expectedOrders for i in self.meshTypes]
|
||||
|
||||
assert type(self.expectedOrders) == list, 'expectedOrders must be a list'
|
||||
assert len(self.expectedOrders) == len(self.meshTypes), 'expectedOrders must have the same length as the meshTypes'
|
||||
|
||||
for ii_meshType, meshType in enumerate(self.meshTypes):
|
||||
self._meshType = meshType
|
||||
self._expectedOrder = self.expectedOrders[ii_meshType]
|
||||
|
||||
order = []
|
||||
err_old = 0.
|
||||
max_h_old = 0.
|
||||
for ii, nc in enumerate(self.meshSizes):
|
||||
max_h = self.setupMesh(nc)
|
||||
err = self.getError()
|
||||
if ii == 0:
|
||||
print ''
|
||||
print 'Testing convergence on ' + self.M._meshType + ' of: ' + self.name
|
||||
print '_____________________________________________'
|
||||
print ' h | error | e(i-1)/e(i) | order'
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
|
||||
print '%4i | %8.2e |' % (nc, err)
|
||||
else:
|
||||
order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
|
||||
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
|
||||
err_old = err
|
||||
max_h_old = max_h
|
||||
print '---------------------------------------------'
|
||||
passTest = np.mean(np.array(order)) > self.tolerance*self._expectedOrder
|
||||
if passTest:
|
||||
print ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!'][np.random.randint(5)]
|
||||
else:
|
||||
print 'Failed to pass test on ' + self._meshType + '.'
|
||||
print ''
|
||||
self.assertTrue(passTest)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,103 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
|
||||
from utils import ndgrid
|
||||
|
||||
|
||||
class BasicLOMTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
a = np.array([1, 1, 1])
|
||||
b = np.array([1, 2])
|
||||
c = np.array([1, 4])
|
||||
gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
|
||||
X, Y = ndgrid(gridIt([a, b]), vector=False)
|
||||
self.TM2 = TensorMesh([a, b])
|
||||
self.LOM2 = LogicallyOrthogonalMesh([X, Y])
|
||||
X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
|
||||
self.TM3 = TensorMesh([a, b, c])
|
||||
self.LOM3 = LogicallyOrthogonalMesh([X, Y, Z])
|
||||
|
||||
def test_area_3D(self):
|
||||
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
|
||||
self.assertTrue(np.all(self.LOM3.area == test_area))
|
||||
|
||||
def test_vol_3D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
|
||||
np.testing.assert_almost_equal(self.LOM3.vol, test_vol)
|
||||
self.assertTrue(True) # Pass if you get past the assertion.
|
||||
|
||||
def test_vol_2D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2])
|
||||
t1 = np.all(self.LOM2.vol == test_vol)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_3D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
|
||||
t1 = np.all(self.LOM3.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_2D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
|
||||
t1 = np.all(self.LOM2.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_tangents(self):
|
||||
T = self.LOM2.tangents
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nE[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nE[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nE[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nE[1])))
|
||||
|
||||
T = self.LOM3.tangents
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nE[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nE[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nE[0])))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nE[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nE[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nE[1])))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nE[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nE[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nE[2])))
|
||||
|
||||
def test_normals(self):
|
||||
N = self.LOM2.normals
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nF[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nF[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nF[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nF[1])))
|
||||
|
||||
N = self.LOM3.normals
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nF[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nF[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nF[0])))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nF[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nF[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nF[1])))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nF[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nF[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nF[2])))
|
||||
|
||||
def test_grid(self):
|
||||
self.assertTrue(np.all(self.LOM2.gridFx == self.TM2.gridFx))
|
||||
self.assertTrue(np.all(self.LOM2.gridFy == self.TM2.gridFy))
|
||||
self.assertTrue(np.all(self.LOM2.gridEx == self.TM2.gridEx))
|
||||
self.assertTrue(np.all(self.LOM2.gridEy == self.TM2.gridEy))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.gridFx == self.TM3.gridFx))
|
||||
self.assertTrue(np.all(self.LOM3.gridFy == self.TM3.gridFy))
|
||||
self.assertTrue(np.all(self.LOM3.gridFz == self.TM3.gridFz))
|
||||
self.assertTrue(np.all(self.LOM3.gridEx == self.TM3.gridEx))
|
||||
self.assertTrue(np.all(self.LOM3.gridEy == self.TM3.gridEy))
|
||||
self.assertTrue(np.all(self.LOM3.gridEz == self.TM3.gridEz))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,210 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from OrderTest import OrderTest
|
||||
|
||||
|
||||
# MATLAB code:
|
||||
|
||||
# syms x y z
|
||||
|
||||
# ex = x.^2+y.*z;
|
||||
# ey = (z.^2).*x+y.*z;
|
||||
# ez = y.^2+x.*z;
|
||||
|
||||
# e = [ex;ey;ez];
|
||||
|
||||
# sigma1 = x.*y+1;
|
||||
# sigma2 = x.*z+2;
|
||||
# sigma3 = 3+z.*y;
|
||||
# sigma4 = 0.1.*x.*y.*z;
|
||||
# sigma5 = 0.2.*x.*y;
|
||||
# sigma6 = 0.1.*z;
|
||||
|
||||
# S1 = [sigma1,0,0;0,sigma1,0;0,0,sigma1];
|
||||
# S2 = [sigma1,0,0;0,sigma2,0;0,0,sigma3];
|
||||
# S3 = [sigma1,sigma4,sigma5;sigma4,sigma2,sigma6;sigma5,sigma6,sigma3];
|
||||
|
||||
# i1 = int(int(int(e.'*S1*e,x,0,1),y,0,1),z,0,1);
|
||||
# i2 = int(int(int(e.'*S2*e,x,0,1),y,0,1),z,0,1);
|
||||
# i3 = int(int(int(e.'*S3*e,x,0,1),y,0,1),z,0,1);
|
||||
|
||||
|
||||
class TestInnerProducts(OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
|
||||
meshDimension = 3
|
||||
meshSizes = [16, 32]
|
||||
|
||||
def getError(self):
|
||||
|
||||
call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
ex = lambda x, y, z: x**2+y*z
|
||||
ey = lambda x, y, z: (z**2)*x+y*z
|
||||
ez = lambda x, y, z: y**2+x*z
|
||||
|
||||
sigma1 = lambda x, y, z: x*y+1
|
||||
sigma2 = lambda x, y, z: x*z+2
|
||||
sigma3 = lambda x, y, z: 3+z*y
|
||||
sigma4 = lambda x, y, z: 0.1*x*y*z
|
||||
sigma5 = lambda x, y, z: 0.2*x*y
|
||||
sigma6 = lambda x, y, z: 0.1*z
|
||||
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 647./360 # Found using matlab symbolic toolbox.
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 37./12 # Found using matlab symbolic toolbox.
|
||||
elif self.sigmaTest == 6:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),
|
||||
call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
|
||||
analytic = 69881./21600 # Found using matlab symbolic toolbox.
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
||||
Ec = np.vstack((cart(self.M.gridEx),
|
||||
cart(self.M.gridEy),
|
||||
cart(self.M.gridEz)))
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
A = self.M.getEdgeInnerProduct(sigma)
|
||||
numeric = E.T.dot(A.dot(E))
|
||||
elif self.location == 'faces':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
||||
Fc = np.vstack((cart(self.M.gridFx),
|
||||
cart(self.M.gridFy),
|
||||
cart(self.M.gridFz)))
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_edges(self):
|
||||
self.name = "Edge Inner Product - Isotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges(self):
|
||||
self.name = "Edge Inner Product - Anisotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges(self):
|
||||
self.name = "Edge Inner Product - Full Tensor"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 6
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "Face Inner Product - Isotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces(self):
|
||||
self.name = "Face Inner Product - Anisotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_faces(self):
|
||||
self.name = "Face Inner Product - Full Tensor"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 6
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestInnerProducts2D(OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
|
||||
meshDimension = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
|
||||
z = 5 # Because 5 is just such a great number.
|
||||
|
||||
call = lambda fun, xy: fun(xy[:, 0], xy[:, 1])
|
||||
|
||||
ex = lambda x, y: x**2+y*z
|
||||
ey = lambda x, y: (z**2)*x+y*z
|
||||
|
||||
sigma1 = lambda x, y: x*y+1
|
||||
sigma2 = lambda x, y: x*z+2
|
||||
sigma3 = lambda x, y: 3+z*y
|
||||
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 144877./360 # Found using matlab symbolic toolbox. z=5
|
||||
elif self.sigmaTest == 2:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
analytic = 189959./120 # Found using matlab symbolic toolbox. z=5
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 781427./360 # Found using matlab symbolic toolbox. z=5
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
||||
Ec = np.vstack((cart(self.M.gridEx),
|
||||
cart(self.M.gridEy)))
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
A = self.M.getEdgeInnerProduct(sigma)
|
||||
numeric = E.T.dot(A.dot(E))
|
||||
elif self.location == 'faces':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
||||
Fc = np.vstack((cart(self.M.gridFx),
|
||||
cart(self.M.gridFy)))
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_edges(self):
|
||||
self.name = "2D Edge Inner Product - Isotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges(self):
|
||||
self.name = "2D Edge Inner Product - Anisotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 2
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges(self):
|
||||
self.name = "2D Edge Inner Product - Full Tensor"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "2D Face Inner Product - Isotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.orderTest()
|
||||
|
||||
def test_order2_faces(self):
|
||||
self.name = "2D Face Inner Product - Anisotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 2
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces(self):
|
||||
self.name = "2D Face Inner Product - Full Tensor"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,162 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from OrderTest import OrderTest
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
|
||||
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
|
||||
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
|
||||
cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
|
||||
cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
|
||||
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
|
||||
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
|
||||
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
|
||||
|
||||
|
||||
class TestCurl(OrderTest):
|
||||
name = "Curl"
|
||||
meshTypes = MESHTYPES
|
||||
|
||||
def getError(self):
|
||||
# fun: i (cos(y)) + j (cos(z)) + k (cos(x))
|
||||
# sol: i (sin(z)) + j (sin(x)) + k (sin(y))
|
||||
|
||||
funX = lambda x, y, z: np.cos(2*np.pi*y)
|
||||
funY = lambda x, y, z: np.cos(2*np.pi*z)
|
||||
funZ = lambda x, y, z: np.cos(2*np.pi*x)
|
||||
|
||||
solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z)
|
||||
solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)
|
||||
solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y)
|
||||
|
||||
Ec = cartE3(self.M, funX, funY, funZ)
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
|
||||
Fc = cartF3(self.M, solX, solY, solZ)
|
||||
curlE_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
curlE = self.M.edgeCurl.dot(E)
|
||||
if self._meshType == 'rotateLOM':
|
||||
# Really it is the integration we should be caring about:
|
||||
# So, let us look at the l2 norm.
|
||||
err = np.linalg.norm(self.M.area*(curlE - curlE_anal), 2)
|
||||
else:
|
||||
err = np.linalg.norm((curlE - curlE_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv(OrderTest):
|
||||
name = "Face Divergence"
|
||||
meshTypes = MESHTYPES
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: np.sin(2*np.pi*x)
|
||||
fy = lambda x, y, z: np.sin(2*np.pi*y)
|
||||
fz = lambda x, y, z: np.sin(2*np.pi*z)
|
||||
sol = lambda x, y, z: (2*np.pi*np.cos(2*np.pi*x)+2*np.pi*np.cos(2*np.pi*y)+2*np.pi*np.cos(2*np.pi*z))
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_anal = call3(sol, self.M.gridCC)
|
||||
|
||||
if self._meshType == 'rotateLOM':
|
||||
# Really it is the integration we should be caring about:
|
||||
# So, let us look at the l2 norm.
|
||||
err = np.linalg.norm(self.M.vol*(divF-divF_anal), 2)
|
||||
else:
|
||||
err = np.linalg.norm((divF-divF_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv2D(OrderTest):
|
||||
name = "Face Divergence 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: np.sin(2*np.pi*x)
|
||||
fy = lambda x, y: np.sin(2*np.pi*y)
|
||||
sol = lambda x, y: 2*np.pi*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y))
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_anal = call2(sol, self.M.gridCC)
|
||||
|
||||
err = np.linalg.norm((divF-divF_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad(OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
meshTypes = MESHTYPES
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
# i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
solX = lambda x, y, z: -np.sin(x)
|
||||
solY = lambda x, y, z: -np.sin(y)
|
||||
solZ = lambda x, y, z: -np.sin(z)
|
||||
|
||||
phi = call3(fun, self.M.gridN)
|
||||
gradE = self.M.nodalGrad.dot(phi)
|
||||
|
||||
Ec = cartE3(self.M, solX, solY, solZ)
|
||||
gradE_anal = self.M.projectEdgeVector(Ec)
|
||||
|
||||
err = np.linalg.norm((gradE-gradE_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad2D(OrderTest):
|
||||
name = "Nodal Gradient 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fun = lambda x, y: (np.cos(x)+np.cos(y))
|
||||
# i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
solX = lambda x, y: -np.sin(x)
|
||||
solY = lambda x, y: -np.sin(y)
|
||||
|
||||
phi = call2(fun, self.M.gridN)
|
||||
gradE = self.M.nodalGrad.dot(phi)
|
||||
|
||||
Ec = cartE2(self.M, solX, solY)
|
||||
gradE_anal = self.M.projectEdgeVector(Ec)
|
||||
|
||||
err = np.linalg.norm((gradE-gradE_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -57,83 +57,6 @@ class BasicTensorMeshTests(unittest.TestCase):
|
||||
t1 = np.all(self.mesh2.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
class TestCurl(OrderTest):
|
||||
name = "Curl"
|
||||
|
||||
def getError(self):
|
||||
fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x))
|
||||
sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y))
|
||||
|
||||
Ex = fun(self.M.gridEx[:, 1])
|
||||
Ey = fun(self.M.gridEy[:, 2])
|
||||
Ez = fun(self.M.gridEz[:, 0])
|
||||
E = np.concatenate((Ex, Ey, Ez))
|
||||
|
||||
Fx = sol(self.M.gridFx[:, 2])
|
||||
Fy = sol(self.M.gridFy[:, 0])
|
||||
Fz = sol(self.M.gridFz[:, 1])
|
||||
curlE_anal = np.concatenate((Fx, Fy, Fz))
|
||||
|
||||
# Generate DIV matrix
|
||||
CURL = self.M.edgeCurl
|
||||
|
||||
curlE = CURL*E
|
||||
err = np.linalg.norm((curlE-curlE_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv(OrderTest):
|
||||
name = "Face Divergence"
|
||||
|
||||
def getError(self):
|
||||
DIV = self.M.faceDiv
|
||||
|
||||
#Test function
|
||||
fun = lambda x: np.sin(x)
|
||||
Fx = fun(self.M.gridFx[:, 0])
|
||||
Fy = fun(self.M.gridFy[:, 1])
|
||||
Fz = fun(self.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(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2])
|
||||
|
||||
err = np.linalg.norm((divF-divF_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad(OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
|
||||
def getError(self):
|
||||
GRAD = self.M.nodalGrad
|
||||
#Test function
|
||||
fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
|
||||
phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2])
|
||||
gradE = GRAD*phi
|
||||
|
||||
Ex = sol(self.M.gridEx[:, 0])
|
||||
Ey = sol(self.M.gridEy[:, 1])
|
||||
Ez = sol(self.M.gridEz[:, 2])
|
||||
|
||||
gradE_anal = np.concatenate((Ex, Ey, Ez))
|
||||
err = np.linalg.norm((gradE-gradE_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestPoissonEqn(OrderTest):
|
||||
name = "Poisson Equation"
|
||||
@@ -167,5 +90,7 @@ class TestPoissonEqn(OrderTest):
|
||||
self.name = "Poisson Equation - Backward"
|
||||
self.forward = False
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
+48
-10
@@ -2,7 +2,8 @@ import numpy as np
|
||||
import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from utils import mkvc, ndgrid
|
||||
from utils import mkvc, ndgrid, indexCube
|
||||
from sputils import sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, sp
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
@@ -30,10 +31,8 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
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)
|
||||
self.assertTrue(np.all(XY[:, 0] == X1_test))
|
||||
self.assertTrue(np.all(XY[:, 1] == X2_test))
|
||||
|
||||
def test_ndgrid_3D(self):
|
||||
XYZ = ndgrid([self.a, self.b, self.c])
|
||||
@@ -42,12 +41,51 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
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(np.all(XYZ[:, 0] == X1_test))
|
||||
self.assertTrue(np.all(XYZ[:, 1] == X2_test))
|
||||
self.assertTrue(np.all(XYZ[:, 2] == X3_test))
|
||||
|
||||
def test_indexCube_2D(self):
|
||||
nN = np.array([3, 3])
|
||||
self.assertTrue(np.all(indexCube('A', nN) == np.array([0, 1, 3, 4])))
|
||||
self.assertTrue(np.all(indexCube('B', nN) == np.array([3, 4, 6, 7])))
|
||||
self.assertTrue(np.all(indexCube('C', nN) == np.array([4, 5, 7, 8])))
|
||||
self.assertTrue(np.all(indexCube('D', nN) == np.array([1, 2, 4, 5])))
|
||||
|
||||
def test_indexCube_3D(self):
|
||||
nN = np.array([3, 3, 3])
|
||||
self.assertTrue(np.all(indexCube('A', nN) == np.array([0, 1, 3, 4, 9, 10, 12, 13])))
|
||||
self.assertTrue(np.all(indexCube('B', nN) == np.array([3, 4, 6, 7, 12, 13, 15, 16])))
|
||||
self.assertTrue(np.all(indexCube('C', nN) == np.array([4, 5, 7, 8, 13, 14, 16, 17])))
|
||||
self.assertTrue(np.all(indexCube('D', nN) == np.array([1, 2, 4, 5, 10, 11, 13, 14])))
|
||||
self.assertTrue(np.all(indexCube('E', nN) == np.array([9, 10, 12, 13, 18, 19, 21, 22])))
|
||||
self.assertTrue(np.all(indexCube('F', nN) == np.array([12, 13, 15, 16, 21, 22, 24, 25])))
|
||||
self.assertTrue(np.all(indexCube('G', nN) == np.array([13, 14, 16, 17, 22, 23, 25, 26])))
|
||||
self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23])))
|
||||
|
||||
def test_invXXXBlockDiagonal(self):
|
||||
|
||||
a = [np.random.rand(5, 1) for i in range(4)]
|
||||
|
||||
B = inv2X2BlockDiagonal(*a)
|
||||
|
||||
A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]))),
|
||||
sp.hstack((sdiag(a[2]), sdiag(a[3])))))
|
||||
|
||||
Z2 = B*A - sp.eye(10, 10)
|
||||
self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-12)
|
||||
|
||||
a = [np.random.rand(5, 1) for i in range(9)]
|
||||
B = inv3X3BlockDiagonal(*a)
|
||||
|
||||
A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]), sdiag(a[2]))),
|
||||
sp.hstack((sdiag(a[3]), sdiag(a[4]), sdiag(a[5]))),
|
||||
sp.hstack((sdiag(a[6]), sdiag(a[7]), sdiag(a[8])))))
|
||||
|
||||
Z3 = B*A - sp.eye(15, 15)
|
||||
|
||||
self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12)
|
||||
|
||||
self.assertTrue(xtest and ytest and ztest)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
|
||||
+235
-4
@@ -1,10 +1,6 @@
|
||||
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
|
||||
|
||||
@@ -22,6 +18,9 @@ def mkvc(x, numDims=1):
|
||||
> (3, 1, 1)
|
||||
|
||||
"""
|
||||
if type(x) == np.matrix:
|
||||
x = np.array(x)
|
||||
|
||||
assert type(x) == np.ndarray, "Vector must be a numpy array"
|
||||
|
||||
if numDims == 1:
|
||||
@@ -97,3 +96,235 @@ def ndgrid(*args, **kwargs):
|
||||
return np.c_[X1, X2, X3]
|
||||
else:
|
||||
return XYZ[2], XYZ[1], XYZ[0]
|
||||
|
||||
|
||||
def volTetra(xyz, A, B, C, D):
|
||||
"""
|
||||
Returns the volume for tetrahedras volume specified by the indexes A to D.
|
||||
|
||||
|
||||
Input:
|
||||
xyz - X,Y,Z vertex vector
|
||||
A,B,C,D - vert index of the tetrahedra
|
||||
|
||||
Output:
|
||||
V - volume
|
||||
|
||||
Algorithm: http://en.wikipedia.org/wiki/Tetrahedron#Volume
|
||||
|
||||
V = 1/3 A * h
|
||||
|
||||
V = 1/6 | ( a - d ) o ( ( b - d ) X ( c - d ) ) |
|
||||
|
||||
"""
|
||||
|
||||
AD = xyz[A, :] - xyz[D, :]
|
||||
BD = xyz[B, :] - xyz[D, :]
|
||||
CD = xyz[C, :] - xyz[D, :]
|
||||
|
||||
V = (BD[:, 0]*CD[:, 1] - BD[:, 1]*CD[:, 0])*AD[:, 2] - (BD[:, 0]*CD[:, 2] - BD[:, 2]*CD[:, 0])*AD[:, 1] + (BD[:, 1]*CD[:, 2] - BD[:, 2]*CD[:, 1])*AD[:, 0]
|
||||
return V/6
|
||||
|
||||
|
||||
def indexCube(nodes, gridSize, n=None):
|
||||
"""
|
||||
Returns the index of nodes on the mesh.
|
||||
|
||||
|
||||
Input:
|
||||
nodes - string of which nodes to return. e.g. 'ABCD'
|
||||
gridSize - size of the nodal grid
|
||||
n - number of nodes each i,j,k direction: [ni,nj,nk]
|
||||
|
||||
|
||||
Output:
|
||||
index - index in the order asked e.g. 'ABCD' --> (A,B,C,D)
|
||||
|
||||
TWO DIMENSIONS:
|
||||
|
||||
node(i,j) node(i,j+1)
|
||||
A -------------- B
|
||||
| |
|
||||
| cell(i,j) |
|
||||
| I |
|
||||
| |
|
||||
D -------------- C
|
||||
node(i+1,j) node(i+1,j+1)
|
||||
|
||||
|
||||
THREE DIMENSIONS:
|
||||
|
||||
node(i,j,k+1) node(i,j+1,k+1)
|
||||
E --------------- F
|
||||
/| / |
|
||||
/ | / |
|
||||
/ | / |
|
||||
node(i,j,k) node(i,j+1,k)
|
||||
A -------------- B |
|
||||
| H ----------|---- G
|
||||
| /cell(i,j) | /
|
||||
| / I | /
|
||||
| / | /
|
||||
D -------------- C
|
||||
node(i+1,j,k) node(i+1,j+1,k)
|
||||
|
||||
|
||||
@author Rowan Cockett
|
||||
|
||||
Last modified on: 2013/07/26
|
||||
"""
|
||||
|
||||
assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'"
|
||||
assert type(gridSize) == np.ndarray, "Number of nodes must be an ndarray"
|
||||
nodes = nodes.upper()
|
||||
# Make sure that we choose from the possible nodes.
|
||||
possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
|
||||
for node in nodes:
|
||||
assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
|
||||
dim = gridSize.size
|
||||
if n is None:
|
||||
n = gridSize - 1
|
||||
|
||||
if dim == 2:
|
||||
ij = ndgrid(np.arange(n[0]), np.arange(n[1]))
|
||||
i, j = ij[:, 0], ij[:, 1]
|
||||
elif dim == 3:
|
||||
ijk = ndgrid(np.arange(n[0]), np.arange(n[1]), np.arange(n[2]))
|
||||
i, j, k = ijk[:, 0], ijk[:, 1], ijk[:, 2]
|
||||
else:
|
||||
raise Exception('Only 2 and 3 dimensions supported.')
|
||||
|
||||
nodeMap = {'A': [0, 0, 0], 'B': [0, 1, 0], 'C': [1, 1, 0], 'D': [1, 0, 0],
|
||||
'E': [0, 0, 1], 'F': [0, 1, 1], 'G': [1, 1, 1], 'H': [1, 0, 1]}
|
||||
out = ()
|
||||
for node in nodes:
|
||||
shift = nodeMap[node]
|
||||
if dim == 2:
|
||||
out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1]]).flatten(), )
|
||||
elif dim == 3:
|
||||
out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1], k+shift[2]]).flatten(), )
|
||||
|
||||
return out
|
||||
|
||||
|
||||
def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
|
||||
"""
|
||||
function [N] = faceInfo(y,A,B,C,D)
|
||||
|
||||
Returns the averaged normal, area, and edge lengths for a given set of faces.
|
||||
|
||||
If average option is FALSE then N is a cell array {nA,nB,nC,nD}
|
||||
|
||||
|
||||
Input:
|
||||
xyz - X,Y,Z vertex vector
|
||||
A,B,C,D - vert index of the face (counter clockwize)
|
||||
|
||||
Options:
|
||||
average - [true]/false, toggles returning all normals or the average
|
||||
|
||||
Output:
|
||||
N - average face normal or {nA,nB,nC,nD} if average = false
|
||||
area - average face area
|
||||
edgeLengths - exact edge Lengths, 4 column vector [AB, BC, CD, DA]
|
||||
|
||||
see also testFaceNormal testFaceArea
|
||||
|
||||
@author Rowan Cockett
|
||||
|
||||
Last modified on: 2013/07/26
|
||||
|
||||
"""
|
||||
assert type(average) is bool, 'average must be a boolean'
|
||||
assert type(normalizeNormals) is bool, 'normalizeNormals must be a boolean'
|
||||
# compute normal that is pointing away from you.
|
||||
#
|
||||
# A -------A-B------- B
|
||||
# | |
|
||||
# | |
|
||||
# D-A (X) B-C
|
||||
# | |
|
||||
# | |
|
||||
# D -------C-D------- C
|
||||
|
||||
AB = xyz[B, :] - xyz[A, :]
|
||||
BC = xyz[C, :] - xyz[B, :]
|
||||
CD = xyz[D, :] - xyz[C, :]
|
||||
DA = xyz[A, :] - xyz[D, :]
|
||||
|
||||
def cross(X, Y):
|
||||
return np.c_[X[:, 1]*Y[:, 2] - X[:, 2]*Y[:, 1],
|
||||
X[:, 2]*Y[:, 0] - X[:, 0]*Y[:, 2],
|
||||
X[:, 0]*Y[:, 1] - X[:, 1]*Y[:, 0]]
|
||||
|
||||
nA = cross(AB, DA)
|
||||
nB = cross(BC, AB)
|
||||
nC = cross(CD, BC)
|
||||
nD = cross(DA, CD)
|
||||
|
||||
length = lambda x: np.sqrt(x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)
|
||||
normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(x), 2))
|
||||
if average:
|
||||
# average the normals at each vertex.
|
||||
N = (nA + nB + nC + nD)/4 # this is intrinsically weighted by area
|
||||
# normalize
|
||||
N = normalize(N)
|
||||
else:
|
||||
if normalizeNormals:
|
||||
N = [normalize(nA), normalize(nB), normalize(nC), normalize(nD)]
|
||||
else:
|
||||
N = [nA, nB, nC, nD]
|
||||
|
||||
# Area calculation
|
||||
#
|
||||
# Approximate by 4 different triangles, and divide by 2.
|
||||
# Each triangle is one half of the length of the cross product
|
||||
#
|
||||
# So also could be viewed as the average parallelogram.
|
||||
#
|
||||
# WARNING: This does not compute correctly for concave quadrilaterals
|
||||
area = (length(nA)+length(nB)+length(nC)+length(nD))/4
|
||||
|
||||
return N, area
|
||||
|
||||
|
||||
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 = np.array(mult).reshape(len(mult))
|
||||
|
||||
sub = []
|
||||
|
||||
for i in range(0, len(shape)):
|
||||
sub.extend([np.math.floor(ind / mult[i])])
|
||||
ind = ind - (np.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 = np.array(mult).reshape(len(mult), 1)
|
||||
|
||||
idx = np.dot((subs), (mult))
|
||||
return idx
|
||||
|
||||
|
||||
def getSubArray(A, ind):
|
||||
"""subArray"""
|
||||
assert type(ind) == list, "ind must be a list of vectors"
|
||||
assert len(A.shape) == len(ind), "ind must have the same length as the dimension of A"
|
||||
|
||||
if len(A.shape) == 2:
|
||||
return A[ind[0], :][:, ind[1]]
|
||||
elif len(A.shape) == 3:
|
||||
return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
|
||||
else:
|
||||
raise Exception("getSubArray does not support dimension asked.")
|
||||
|
||||
File diff suppressed because one or more lines are too long
Reference in New Issue
Block a user