mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 16:32:43 +08:00
Merge branch 'eldadswork' of https://bitbucket.org/rcockett/simpeg into LOM
Conflicts: SimPEG/utils.py
This commit is contained in:
Rowan Cockett
@@ -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)."
|
||||
|
||||
@@ -0,0 +1,187 @@
|
||||
from scipy import sparse as sp
|
||||
from sputils import sdiag
|
||||
from utils import sub2ind, ndgrid, mkvc
|
||||
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):
|
||||
if self._meshType == 'TENSOR':
|
||||
pass
|
||||
elif self._meshType == 'LOM':
|
||||
pass # todo: we should be doing something slightly different here!
|
||||
return getFaceInnerProduct(self, mu)
|
||||
|
||||
def getEdgeInnerProduct(self, sigma):
|
||||
if self._meshType == 'TENSOR':
|
||||
pass
|
||||
elif self._meshType == 'LOM':
|
||||
pass # todo: we should be doing something slightly different here!
|
||||
return getEdgeInnerProduct(self, sigma)
|
||||
|
||||
|
||||
def getFaceInnerProduct(mesh, mu):
|
||||
|
||||
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]
|
||||
|
||||
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()
|
||||
|
||||
return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr()
|
||||
|
||||
# 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 | 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
|
||||
P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
|
||||
P110 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
|
||||
P001 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]])
|
||||
P101 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]])
|
||||
P011 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]])
|
||||
P111 = 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))
|
||||
|
||||
# Cell volume
|
||||
v = np.sqrt(mesh.vol)
|
||||
v3 = np.r_[v, v, v]
|
||||
V = sdiag(v3)*mu*sdiag(v3) # to keep symmetry
|
||||
|
||||
A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
|
||||
|
||||
A = 0.125*A
|
||||
|
||||
return A
|
||||
|
||||
|
||||
def getEdgeInnerProduct(mesh, sigma):
|
||||
|
||||
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]
|
||||
|
||||
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()
|
||||
|
||||
return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
|
||||
|
||||
# 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
|
||||
P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]])
|
||||
P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]])
|
||||
P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]])
|
||||
P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]])
|
||||
P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]])
|
||||
P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]])
|
||||
P111 = 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))
|
||||
|
||||
# Cell volume
|
||||
v = np.sqrt(mesh.vol)
|
||||
v3 = np.r_[v, v, v]
|
||||
V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry
|
||||
|
||||
A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
|
||||
|
||||
A = 0.125*A
|
||||
|
||||
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 = mesh.getFaceInnerProduct(mu)
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -1,87 +0,0 @@
|
||||
from scipy import sparse as sp
|
||||
from sputils import sdiag
|
||||
from utils import sub2ind, ndgrid, mkvc
|
||||
import numpy as np
|
||||
|
||||
|
||||
def getEdgeInnerProduct(mesh, sigma):
|
||||
|
||||
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]
|
||||
|
||||
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()
|
||||
|
||||
return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
|
||||
|
||||
# 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
|
||||
P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]])
|
||||
P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]])
|
||||
P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]])
|
||||
P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]])
|
||||
P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]])
|
||||
P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]])
|
||||
P111 = 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))
|
||||
|
||||
# Cell volume
|
||||
v = np.sqrt(mesh.vol)
|
||||
v3 = np.r_[v, v, v]
|
||||
V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry
|
||||
|
||||
A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
|
||||
|
||||
A = 0.125*A
|
||||
|
||||
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)
|
||||
sigma = np.ones((mesh.nC, 6))
|
||||
A = getEdgeInnerProduct(mesh, sigma)
|
||||
@@ -0,0 +1,92 @@
|
||||
from scipy import sparse as sp
|
||||
import numpy as np
|
||||
|
||||
def interpmat(x,y,z,xr,yr,zr):
|
||||
#
|
||||
# This function does local linear interpolation
|
||||
# computed for each receiver point in turn
|
||||
#
|
||||
# [Q] = linint(x,y,z,xr,yr,zr)
|
||||
# Interpolation matrix
|
||||
#
|
||||
|
||||
nx = np.size(x)
|
||||
ny = np.size(y)
|
||||
nz = np.size(z)
|
||||
|
||||
nps = np.size(xr)
|
||||
|
||||
#Q = spalloc(np,nx*ny*nz,8*np);
|
||||
Q = sp.lil_matrix((nps,nx*ny*nz))
|
||||
ind_x = np.array([0,0])
|
||||
ind_y = np.array([0,0])
|
||||
ind_z = np.array([0,0])
|
||||
dx, dy, dz = np.zeros(2), np.zeros(2), np.zeros(2)
|
||||
for i in range(0, nps):
|
||||
im = np.amin(abs(xr[i]-x))
|
||||
if xr[i] - x[im] >= 0: # Point on the left
|
||||
ind_x[0] = im; ind_x[1] = im+1
|
||||
else: # Point on the right
|
||||
ind_x[0] = im-1; ind_x[1] = im
|
||||
|
||||
|
||||
dx[0] = xr[i] - x[ind_x[0]]
|
||||
dx[1] = x[ind_x[1]] - xr[i]
|
||||
|
||||
im = np.amin(abs(yr[i] - y))
|
||||
if yr[i] - y[im] >= 0: # Point on the left
|
||||
ind_y[0] = im; ind_y[1] = im+1
|
||||
else: # Point on the right
|
||||
ind_y[0] = im-1; ind_y[1] = im
|
||||
|
||||
|
||||
dy[0] = yr[i] - y[ind_y[0]]
|
||||
dy[1] = y[ind_y[1]] - yr[i];
|
||||
|
||||
im = np.amin(abs(zr[i] - z));
|
||||
if zr[i] -z[im] >= 0: # Point on the left
|
||||
ind_z[0] = im; ind_z[1] = im+1
|
||||
else: # Point on the right
|
||||
ind_z[0] = im-1; ind_z[1] = im;
|
||||
|
||||
dz[0] = zr[i] - z[ind_z[0]]; dz[1] = z[ind_z[1]] - zr[i]
|
||||
|
||||
Dx = x[ind_x[1]] - x[ind_x[0]]
|
||||
Dy = y[ind_y[1]] - y[ind_y[0]]
|
||||
Dz = z[ind_z[1]] - z[ind_z[0]]
|
||||
#dv = Dx*Dy*Dz
|
||||
|
||||
# Get the row in the matrix
|
||||
v = np.zeros([nx, ny,nz])
|
||||
|
||||
v[ ind_x[0], ind_y[0], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[0], ind_y[1], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[1], ind_y[0], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[1], ind_y[1], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz)
|
||||
v[ ind_x[0], ind_y[0], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[0], ind_y[1], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[1], ind_y[0], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz)
|
||||
v[ ind_x[1], ind_y[1], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz)
|
||||
|
||||
|
||||
print(np.shape(v.flatten('F')))
|
||||
print(np.shape(Q))
|
||||
|
||||
Q[i,:] = v.flatten('F')
|
||||
|
||||
|
||||
return Q.tocsr()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
x = np.array([1, 2, 3, 4])
|
||||
y = np.array([1, 2, 3, 4, 5])
|
||||
z = np.array([0, 1, 4, 6])
|
||||
|
||||
xr = np.array([2.5,3.2])
|
||||
yr = np.array([2.4,3.6])
|
||||
zr = np.array([2.5,3.9])
|
||||
|
||||
A = interpmat(x,y,z,xr,yr,zr)
|
||||
|
||||
+28
-12
@@ -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:
|
||||
@@ -63,17 +63,32 @@ class OrderTest(unittest.TestCase):
|
||||
|
||||
name = "Order Test"
|
||||
expectedOrder = 2
|
||||
tolerance = 0.85
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
meshType = 'uniformTensorMesh'
|
||||
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
|
||||
|
||||
def getError(self):
|
||||
"""For given h, generate A[h], f and A(f) and return norm of error."""
|
||||
@@ -89,9 +104,9 @@ class OrderTest(unittest.TestCase):
|
||||
"""
|
||||
order = []
|
||||
err_old = 0.
|
||||
nc_old = 0.
|
||||
max_h_old = 0.
|
||||
for ii, nc in enumerate(self.meshSizes):
|
||||
self.setupMesh(nc)
|
||||
max_h = self.setupMesh(nc)
|
||||
err = self.getError()
|
||||
if ii == 0:
|
||||
print ''
|
||||
@@ -101,13 +116,14 @@ class OrderTest(unittest.TestCase):
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
|
||||
print '%4i | %8.2e |' % (nc, err)
|
||||
else:
|
||||
order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc)))
|
||||
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
|
||||
nc_old = nc
|
||||
max_h_old = max_h
|
||||
print '---------------------------------------------'
|
||||
self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order)))
|
||||
|
||||
passTest = np.mean(np.array(order)) > self.tolerance*self.expectedOrder
|
||||
# passTest = len(np.where(np.array(order) > self.tolerance*self.expectedOrder)[0]) > np.floor(0.75*len(order))
|
||||
self.assertTrue(passTest)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -1,15 +1,36 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from OrderTest import OrderTest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from getEdgeInnerProducts import *
|
||||
|
||||
|
||||
class TestEdgeInnerProduct(OrderTest):
|
||||
"""Integrate an edge function over a unit cube domain using edgeInnerProducts."""
|
||||
# MATLAB code:
|
||||
|
||||
name = "Edge Inner Product"
|
||||
# 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."""
|
||||
|
||||
def getError(self):
|
||||
|
||||
@@ -26,22 +47,70 @@ class TestEdgeInnerProduct(OrderTest):
|
||||
sigma5 = lambda x, y, z: 0.2*x*y
|
||||
sigma6 = lambda x, y, z: 0.1*z
|
||||
|
||||
Ex = call(ex, self.M.gridEx)
|
||||
Ey = call(ey, self.M.gridEy)
|
||||
Ez = call(ez, self.M.gridEz)
|
||||
|
||||
E = np.matrix(np.r_[Ex, Ey, Ez]).T
|
||||
Gc = self.M.gridCC
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),
|
||||
call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
|
||||
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':
|
||||
Ex = call(ex, self.M.gridEx)
|
||||
Ey = call(ey, self.M.gridEy)
|
||||
Ez = call(ez, self.M.gridEz)
|
||||
E = np.matrix(np.r_[Ex, Ey, Ez]).T
|
||||
A = self.M.getEdgeInnerProduct(sigma)
|
||||
numeric = E.T*A*E
|
||||
elif self.location == 'faces':
|
||||
Fx = call(ex, self.M.gridFx)
|
||||
Fy = call(ey, self.M.gridFy)
|
||||
Fz = call(ez, self.M.gridFz)
|
||||
F = np.matrix(np.r_[Fx, Fy, Fz]).T
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T*A*F
|
||||
|
||||
A = getEdgeInnerProduct(self.M, sigma)
|
||||
numeric = E.T*A*E
|
||||
analytic = 69881./21600 # Found using matlab symbolic toolbox.
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
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()
|
||||
|
||||
|
||||
|
||||
+10
-11
@@ -1,5 +1,4 @@
|
||||
import numpy as np
|
||||
from numpy import *
|
||||
|
||||
|
||||
def mkvc(x, numDims=1):
|
||||
@@ -283,11 +282,6 @@ def faceInfo(xyz, A, B, C, D, average=True):
|
||||
return N, area, edgeLengths
|
||||
|
||||
|
||||
def getSubArray(A, ind):
|
||||
"""subArray"""
|
||||
return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
|
||||
|
||||
|
||||
def ind2sub(shape, ind):
|
||||
"""From the given shape, returns the subscrips of the given index"""
|
||||
revshp = []
|
||||
@@ -295,13 +289,13 @@ def ind2sub(shape, ind):
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult))
|
||||
mult = np.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])
|
||||
sub.extend([np.math.floor(ind / mult[i])])
|
||||
ind = ind - (np.math.floor(ind/mult[i]) * mult[i])
|
||||
return sub
|
||||
|
||||
|
||||
@@ -311,7 +305,12 @@ def sub2ind(shape, subs):
|
||||
mult = [1]
|
||||
for i in range(0, len(revshp)-1):
|
||||
mult.extend([mult[i]*revshp[i]])
|
||||
mult = array(mult).reshape(len(mult), 1)
|
||||
mult = np.array(mult).reshape(len(mult), 1)
|
||||
|
||||
idx = dot((subs), (mult))
|
||||
idx = np.dot((subs), (mult))
|
||||
return idx
|
||||
|
||||
|
||||
def getSubArray(A, ind):
|
||||
"""subArray"""
|
||||
return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
|
||||
|
||||
Reference in New Issue
Block a user