-
Notifications
You must be signed in to change notification settings - Fork 0
/
solver.py
93 lines (67 loc) · 2.14 KB
/
solver.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.sparse import diags
def solve_poisson(U, Rho=None, Eps=None):
Nr, Nc = U.shape
if Rho is None:
Rho = np.zeros_like(U)
if Eps is None:
Eps = np.ones_like(U)
Epsil = np.r_[np.ones(Nc), Eps.flatten(), 1]
Bound = np.zeros_like(U, dtype=bool)
Bound[0, :] = True
Bound[Nr-1, :] = True
Bound[:, 0] = True
Bound[:, Nc-1] = True
Bound[abs(U) > 1e-7] = True
Bound = Bound.flatten()
s = [-(Epsil[Nc:Nr*Nc+Nc] + Epsil[:Nr*Nc] +
Epsil[1+Nc:Nr*Nc+Nc+1] + Epsil[1:Nr*Nc+1]),
0.5*(Epsil[1:Nr*Nc] + Epsil[1+Nc:Nr*Nc+Nc]),
0.5*(Epsil[1+Nc:Nr*Nc+Nc] + Epsil[1:Nr*Nc]),
0.5*(Epsil[1+Nc:Nr*Nc+1] + Epsil[Nc:Nr*Nc]),
0.5*(Epsil[Nc:Nr*Nc] + Epsil[1+Nc:Nr*Nc+1])]
s[0][Bound] = 1
s[1][Bound[1:Nr*Nc]] = 0
s[2][Bound[:Nr*Nc-1]] = 0
s[3][Bound[Nc:Nr*Nc]] = 0
s[4][Bound[:Nr*Nc-Nc]] = 0
S = diags(s, [0, -1, 1, -Nc, Nc])
S = csr_matrix(S)
b = U.flatten() - Rho.flatten()
phi = spsolve(S, b)
phi = phi.reshape(Nr, Nc)
return phi
def solve_poisson_(U, Rho=None, Eps=None):
N, M = U.shape
assert N == M
if Rho is None:
Rho = np.zeros_like(U)
if Eps is None:
Eps = np.ones_like(U)
Epsil = np.r_[np.ones(N), Eps.flatten(), 1]
Bound = np.zeros_like(U, dtype=bool)
Bound[0, :] = True
Bound[N-1, :] = True
Bound[:, 0] = True
Bound[:, N-1] = True
Bound[abs(U) > 1e-7] = True
Bound = Bound.flatten()
s = [-(Epsil[N:N*N+N] + Epsil[:N*N] +
Epsil[1+N:N*N+N+1] + Epsil[1:N*N+1]),
0.5*(Epsil[1:N*N] + Epsil[1+N:N*N+N]),
0.5*(Epsil[1+N:N*N+N] + Epsil[1:N*N]),
0.5*(Epsil[1+N:N*N+1] + Epsil[N:N*N]),
0.5*(Epsil[N:N*N] + Epsil[1+N:N*N+1])]
s[0][Bound] = 1
s[1][Bound[1:N*N]] = 0
s[2][Bound[:N*N-1]] = 0
s[3][Bound[N:N*N]] = 0
s[4][Bound[:N*N-N]] = 0
S = diags(s, [0, -1, 1, -N, N])
S = csr_matrix(S)
b = U.flatten() - Rho.flatten()
phi = spsolve(S, b)
phi = phi.reshape(N, N)
return phi