-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathlaplace_biem_analytic_test
99 lines (78 loc) · 2.61 KB
/
laplace_biem_analytic_test
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
94
95
96
97
98
99
import colorcet as cc
import itertools as it
import matplotlib.pyplot as plt; plt.ion()
import numpy as np
N = 200
a = 2
b = 1
nxgrid, nygrid = 100, 100
xmin, xmax, ymin, ymax = -2.2, 2.2, -2.2, 2.2
theta = np.linspace(0, 2*np.pi, N + 1, endpoint=True)
x = a*np.cos(theta)
y = b*np.sin(theta)
h = 1 / (N) * 2 * np.pi
theta = theta[:-1]
x = a*np.cos(theta)
y = b*np.sin(theta)
tx = -a*np.sin(theta)
ty = b*np.cos(theta)
ax = -a*np.cos(theta)
ay = -b*np.sin(theta)
nx = ax/np.sqrt(ax**2 + ay**2)
ny = ax/np.sqrt(ax**2 + ay**2)
w = h*np.sqrt(tx**2 + ty**2)
f = (x**2 - y**2) / (2*np.pi)
K = np.empty((N, N), dtype=np.float64)
K[...] = np.nan
for i, j in it.product(range(N), range(N)):
if i == j:
kappa = (tx[j]*ax[j] - ty[j]*ay[j])/(tx[j]**2 + ty[j]**2)**1.5
kappa = 1
K[i, j] = -kappa/(4*np.pi) * w[j]
else:
dx = x[j] - x[i]
dy = y[j] - y[i]
dot = nx[j]*dx + ny[j]*dy
r = np.sqrt(dx**2 + dy**2)
K[i, j] = dot/(2*np.pi*r**2)
K[i, j] *= w[j]
sigma = np.linalg.solve(-np.eye(N)/2 + K, f)
X, Y = np.meshgrid(np.linspace(xmin, xmax, nxgrid), np.linspace(ymin, ymax, nygrid), indexing='xy')
u = np.empty((nxgrid, nygrid), dtype=np.float64)
for i, j in it.product(range(nxgrid), range(nygrid)):
if (X[i, j]/a)**2 + (Y[i ,j]/b)**2 >= 1:
u[i, j] = np.nan
continue
u[i, j] = 0
for k in range(N):
dx = X[i, j] - x[k]
dy = Y[i, j] - y[k]
dot = nx[k]*dx + ny[k]*dy
r = np.sqrt(dx**2 + dy**2)
u[i, j] -= dot/(2*np.pi*r**2)*w[k]*sigma[k]
U = (X**2 - Y**2) / (2*np.pi)
E = u - U
Uvmax = abs(U[np.isfinite(E)]).max()
Uvmin = -Uvmax
Evmax = np.nanmax(abs(E))
Evmin = -Evmax
# Set up the figure and subplots
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(13, 4))
# Numerical solution subplot
num_sol_ax = axes[0]
img1 = num_sol_ax.imshow(u, extent=[xmin, xmax, ymin, ymax], cmap=cc.cm.gouldian, vmax=Uvmax, vmin=Uvmin)
num_sol_ax.set_title('Numerical Solution')
fig.colorbar(img1, ax=num_sol_ax)
# Analytical solution subplot
analytic_sol_ax = axes[1]
img2 = analytic_sol_ax.imshow(U, extent=[xmin, xmax, ymin, ymax], cmap=cc.cm.gouldian, vmax=Uvmax, vmin=Uvmin)
analytic_sol_ax.set_title('Analytical Solution')
fig.colorbar(img2, ax=analytic_sol_ax)
# Error subplot
error_ax = axes[2]
img3 = error_ax.imshow(E, extent=[xmin, xmax, ymin, ymax], cmap=cc.cm.coolwarm, vmin=Evmin, vmax=Evmax)
error_ax.set_title('Error')
fig.colorbar(img3, ax=error_ax)
# Adjust layout and display the plot
plt.tight_layout()
plt.show()