In [1]:
import rectmesh
import numpy as np
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg
import matplotlib.pyplot as plt
In [41]:
%matplotlib inline
import numpy as np
n1 = 100
n2 = 100
x = np.linspace(0,np.pi,n1)
y = np.linspace(0,np.pi,n2)
# x = np.sort(np.random.random(n1))
# y = np.sort(np.random.random(n2))
# x[0] = 0
# x[-1] = 1.
# y[0] = 0
# y[-1] = 1.
nodes = np.array([[0, 0],
[0, n2-1],
[n1/2 - 1, n2-1],
[n1/2 - 1,n2/2-1],
[n1 - 1, n2/2-1],
[n1-1, 0],
[0, 0]])
boundary = [['N', 1], ['D', 1], ['N', 1], ['N', 1], ['D', 1], ['D', 1]]
# nodes = np.array([[0,0], [0,n2-1], [n1-1,n2-1], [n1-1,0], [0, 0]])
# boundary = [['D', 1], ['N', 1], ['N', 1], ['D', 1]]
a = rectmesh.rectmesh(meshx = x, meshy = y, nodes = nodes, boundary = boundary)
a.create_boundary()
a.create_order()
a.plot()
In [42]:
# Solving \Delta u = f
def build_system(mesh, rhs):
A = sp.sparse.lil_matrix((mesh.mask_nnz, mesh.mask_nnz))
x = mesh.x
y = mesh.y
order = mesh.order
order_inv = dict((v,k) for k, v in order.iteritems())
f = np.zeros(mesh.mask_nnz)
#dx_f = x[1:mesh.size[0]-1] - x[:mesh.size[0]-2]
for k in xrange(mesh.num_inside):
[i,j] = order_inv[k]
dx_f = x[i + 1] - x[i]
dx_b = x[i] - x[i - 1]
dx_c = x[i + 1] - x[i - 1]
dy_f = y[j + 1] - y[j]
dy_b = y[j] - y[j - 1]
dy_c = y[j + 1] - y[j - 1]
k_if_jc = order[(i+1, j)]
k_ib_jc = order[(i-1, j)]
k_ic_jf = order[(i, j+1)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_f/dx_b - 2./dy_f/dy_b
A[k, k_if_jc] = 2./dx_f/dx_c
A[k, k_ib_jc] = 2./dx_b/dx_c
A[k, k_ic_jf] = 2./dy_f/dy_c
A[k, k_ic_jb] = 2./dy_b/dy_c
f[k] = rhs[i,j]
mask_ext = np.zeros((mesh.mesh_size[0]+2, mesh.mesh_size[1]+2))
mask_ext[1:mesh.mesh_size[0]+1, 1:mesh.mesh_size[1]+1] = mesh.mask
for k in xrange(mesh.num_inside, mesh.num_inside + mesh.num_neumann):
[i,j] = order_inv[k]
m_if = mask_ext[i+2, j+1]
m_ib = mask_ext[i, j+1]
m_jf = mask_ext[i+1, j+2]
m_jb = mask_ext[i+1, j]
if mask_ext[i+1, j+1] == 3.:
if m_ib==0:
dx_f = x[i + 1] - x[i]
dy_f = y[j + 1] - y[j]
dy_b = y[j] - y[j - 1]
dy_c = y[j + 1] - y[j - 1]
k_if_jc = order[(i+1, j)]
k_ic_jf = order[(i, j+1)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_f**2 - 2./dy_f/dy_b
A[k, k_if_jc] = 2./dx_f**2
A[k, k_ic_jf] = 2./dy_f/dy_c
A[k, k_ic_jb] = 2./dy_b/dy_c
f[k] = rhs[i,j]
elif m_if==0:
dx_b = x[i] - x[i - 1]
dy_f = y[j + 1] - y[j]
dy_b = y[j] - y[j - 1]
dy_c = y[j + 1] - y[j - 1]
k_ib_jc = order[(i-1, j)]
k_ic_jf = order[(i, j+1)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_b**2 - 2./dy_f/dy_b
A[k, k_ib_jc] = 2./dx_b**2
A[k, k_ic_jf] = 2./dy_f/dy_c
A[k, k_ic_jb] = 2./dy_b/dy_c
f[k] = rhs[i,j]
elif m_jb==0:
dx_f = x[i + 1] - x[i]
dx_b = x[i] - x[i - 1]
dx_c = x[i + 1] - x[i - 1]
dy_f = y[j + 1] - y[j]
k_if_jc = order[(i+1, j)]
k_ib_jc = order[(i-1, j)]
k_ic_jf = order[(i, j+1)]
A[k, k] = -2./dx_f/dx_b - 2./dy_f**2
A[k, k_if_jc] = 2./dx_f/dx_c
A[k, k_ib_jc] = 2./dx_b/dx_c
A[k, k_ic_jf] = 2./dy_f**2
f[k] = rhs[i,j]
elif m_jf==0:
dx_f = x[i + 1] - x[i]
dx_b = x[i] - x[i - 1]
dx_c = x[i + 1] - x[i - 1]
dy_b = y[j] - y[j - 1]
k_if_jc = order[(i+1, j)]
k_ib_jc = order[(i-1, j)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_f/dx_b - 2./dy_b**2
A[k, k_if_jc] = 2./dx_f/dx_c
A[k, k_ib_jc] = 2./dx_b/dx_c
A[k, k_ic_jb] = 2./dy_b**2
f[k] = rhs[i,j]
else:
raise Exception('Problems with mask')
elif mask_ext[i+1, j+1] == 4.:
if m_ib == 0 and m_jb == 0:
dx_f = x[i + 1] - x[i]
dy_f = y[j + 1] - y[j]
k_if_jc = order[(i+1, j)]
k_ic_jf = order[(i, j+1)]
A[k, k] = -2./dx_f**2 - 2./dy_f**2
A[k, k_if_jc] = 2./dx_f**2
A[k, k_ic_jf] = 2./dy_f**2
f[k] = rhs[i,j]
elif m_if == 0 and m_jb == 0:
dx_b = x[i] - x[i - 1]
dy_f = y[j + 1] - y[j]
k_ib_jc = order[(i-1, j)]
k_ic_jf = order[(i, j+1)]
A[k, k] = -2./dx_b**2 - 2./dy_f**2
A[k, k_ib_jc] = 2./dx_b**2
A[k, k_ic_jf] = 2./dy_f**2
f[k] = rhs[i,j]
elif m_ib == 0 and m_jf == 0:
dx_f = x[i + 1] - x[i]
dy_b = y[j] - y[j - 1]
k_if_jc = order[(i+1, j)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_f**2 - 2./dy_b**2
A[k, k_if_jc] = 2./dx_f**2
A[k, k_ic_jb] = 2./dy_b**2
f[k] = rhs[i,j]
elif m_if == 0 and m_jf == 0:
dx_b = x[i] - x[i - 1]
dy_b = y[j] - y[j - 1]
k_ib_jc = order[(i-1, j)]
k_ic_jb = order[(i, j-1)]
A[k, k] = -2./dx_b**2 - 2./dy_b**2
A[k, k_ib_jc] = 2./dx_b**2
A[k, k_ic_jb] = 2./dy_b**2
f[k] = rhs[i,j]
else:
raise Exception('Problems with mask')
else:
raise Exception('Problems with mask')
for k in xrange(mesh.num_inside + mesh.num_neumann, mesh.num_inside + mesh.num_neumann + mesh.num_dirichlet):
A[k,k] = 1.
f[k] = mesh.dirichlet_values[k]
return A, f
def tomeshorder(mesh, sol):
size = mesh.mesh_size
sol2mesh = np.zeros(size)
for i in xrange(size[0]):
for j in xrange(size[1]):
try:
sol2mesh[i, j] = sol[mesh.order[i,j]]
except:
sol2mesh[i, j] = 0.
return sol2mesh
In [43]:
# Testing on [0, pi]x[0, pi] square with solution sin x cos y + 1
Y, X = np.meshgrid(x,y)
rhs = np.zeros((len(y), len(x)))
#rhs = - np.sin(X/2)* np.sin(Y/2)/2
rhs[:, len(y)/4] = -len(x)**2
A, f = build_system(a, rhs.T)
A = A.tocsc()
In [44]:
lu = sp.sparse.linalg.splu(A)
sol = lu.solve(f)
In [45]:
sol2mesh = tomeshorder(a,sol)
In [46]:
plt.imshow(sol2mesh)
plt.colorbar()
Out[46]:
In [26]:
sol_exact = np.sin(X/2)* np.sin(Y/2) + 1
np.linalg.norm(sol_exact - sol2mesh)/np.linalg.norm(sol_exact)
print (sol2mesh[n1-1,n2-1] - sol_exact[n1-1,n2-1])/sol_exact[n1-1,n2-1]
In [103]:
7.3137488292270522e-06/1.8093163949806047e-06
Out[103]:
In [114]:
1.04897125286e-05/2.59611505138e-06
Out[114]:
In [ ]: