corrected bugs

This commit is contained in:
ehaber99
2013-08-02 01:51:10 -07:00
parent dcae030d81
commit ec4358d872
2 changed files with 15 additions and 11 deletions
+7 -5
View File
@@ -169,14 +169,14 @@ def getEdgeInnerProduct(mesh, sigma):
# Cell volume
v = np.sqrt(mesh.vol)
v3 = np.r_[v, v, v]
#V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry
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 = 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
A = 0.125*A
P = sp.vstack((sdiag(v3)*P000,sdiag(v3)*P001,sdiag(v3)*P010,sdiag(v3)*P011,
sdiag(v3)*P100,sdiag(v3)*P101,sdiag(v3)*P110,sdiag(v3)*P111))
A = 0.125* (P.T*Sigma*P)
return A, P
@@ -186,4 +186,6 @@ if __name__ == '__main__':
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)
A = getFaceInnerProduct(mesh,mu)
B = getEdgeInnerProduct(mesh,mu)
+8 -6
View File
@@ -1,6 +1,7 @@
from scipy import sparse as sp
import numpy as np
def interpmat(x,y,z,xr,yr,zr):
#
# This function does local linear interpolation
@@ -23,7 +24,8 @@ def interpmat(x,y,z,xr,yr,zr):
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))
im = np.argmin(abs(xr[i]-x))
print i,im
if xr[i] - x[im] >= 0: # Point on the left
ind_x[0] = im; ind_x[1] = im+1
else: # Point on the right
@@ -33,7 +35,7 @@ def interpmat(x,y,z,xr,yr,zr):
dx[0] = xr[i] - x[ind_x[0]]
dx[1] = x[ind_x[1]] - xr[i]
im = np.amin(abs(yr[i] - y))
im = np.argmin(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
@@ -43,7 +45,7 @@ def interpmat(x,y,z,xr,yr,zr):
dy[0] = yr[i] - y[ind_y[0]]
dy[1] = y[ind_y[1]] - yr[i];
im = np.amin(abs(zr[i] - z));
im = np.argmin(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
@@ -80,9 +82,9 @@ def interpmat(x,y,z,xr,yr,zr):
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])
x = np.array([1.1, 2.1, 3.6, 4.9])
y = np.array([1.2, 2.2, 3.3, 4.9, 5.6])
z = np.array([0.8, 1.7, 4.9, 6.5])
xr = np.array([2.5,3.2])
yr = np.array([2.4,3.6])