-
Notifications
You must be signed in to change notification settings - Fork 0
/
problems.py
118 lines (97 loc) · 4.25 KB
/
problems.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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
import numpy as np
from scipy.special import erf, erfc
def smooth_square():
def u_neumann(t, x_hat):
""" evaluates the neumann trace along the lateral boundary. """
return -np.pi * np.exp(-2 * np.pi**2 * t) * np.sin(np.pi * (x_hat % 1))
def u0(xy):
return np.sin(np.pi * xy[0]) * np.sin(np.pi * xy[1])
def M0u0(t, xy):
x = xy[0]
y = xy[1]
pit = np.pi * t
sqrtt = np.sqrt(t)
return (((-(1 / 16)) * (erf((x - 2 * 1j * pit) / (2 * sqrtt)) + erf(
(1 - x + 2 * 1j * pit) /
(2 * sqrtt)) - np.exp(2 * 1j * x * np.pi) * (erf(
(1 - x - 2 * 1j * pit) / (2 * sqrtt)) + erf(
(x + 2 * 1j * pit) / (2 * sqrtt)))) * (erf(
(y - 2 * 1j * pit) / (2 * sqrtt)) + erf(
(1 - y + 2 * 1j * pit) /
(2 * sqrtt)) - np.exp(2 * 1j * y * np.pi) * (erf(
(1 - y - 2 * 1j * pit) / (2 * sqrtt)) + erf(
(y + 2 * 1j * pit) / (2 * sqrtt))))) /
np.exp(1j * np.pi * (x + y - 2 * 1j * pit))).real
return {'u-trace': u_neumann, 'u0': u0, 'M0u0': M0u0}
def smooth_pisquare():
def u_neumann(t, x_hat):
""" evaluates the neumann trace along the lateral boundary. """
return -np.exp(-2 * t) * np.sin((x_hat % np.pi))
def u0(xy):
return np.sin(xy[0]) * np.sin(xy[1])
def M0u0(t, xy):
x = xy[0]
y = xy[1]
return ((-(1 / 16)) * np.exp((-1j) * (x + y) - 2 * t) * (-1 + erf(
(x - 2 * 1j * t) / (2 * np.sqrt(t))) + np.exp(2 * 1j * x) * (-erf(
(x + 2 * 1j * t) / (2 * np.sqrt(t))) + erf(
(x - np.pi + 2 * 1j * t) / (2 * np.sqrt(t)))) + erfc(
(x - np.pi - 2 * 1j * t) / (2 * np.sqrt(t)))) *
(-1 + erf(
(y - 2 * 1j * t) / (2 * np.sqrt(t))) + np.exp(2 * 1j * y) *
(-erf((y + 2 * 1j * t) / (2 * np.sqrt(t))) + erf(
(y - np.pi + 2 * 1j * t) / (2 * np.sqrt(t)))) + erfc(
(y - np.pi - 2 * 1j * t) / (2 * np.sqrt(t))))).real
return {'u-trace': u_neumann, 'u0': u0, 'M0u0': M0u0}
def singular_square():
def M0u0(t, xy):
a = xy[0]
b = xy[1]
return (1 / 4) * (erf(
(1 - a) / (2 * np.sqrt(t))) + erf(a / (2 * np.sqrt(t)))) * (erf(
(1 - b) / (2 * np.sqrt(t))) + erf(b / (2 * np.sqrt(t))))
return {'u0': lambda xy: 1, 'M0u0': M0u0}
def singular_lshape():
def M0u0(t, xy):
a = xy[0]
b = xy[1]
return (1 / 4) * ((erf((1 - a) / (2 * np.sqrt(t))) + erf(
(1 + a) / (2 * np.sqrt(t)))) * (erf(
(1 - b) /
(2 * np.sqrt(t))) + erf(b / (2 * np.sqrt(t)))) + (erf(
(1 - a) / (2 * np.sqrt(t))) + erf(a / (2 * np.sqrt(t)))) *
(-erf(b / (2 * np.sqrt(t))) + erf(
(1 + b) / (2 * np.sqrt(t)))))
return {'u0': lambda xy: 1, 'M0u0': M0u0}
def problem_helper(problem, domain):
assert problem in ['Smooth', 'Dirichlet', 'Singular', 'MildSingular']
assert domain in ['UnitSquare', 'PiSquare', 'LShape', 'Circle']
result = {}
if problem == 'Smooth':
if domain == 'UnitSquare':
result.update(smooth_square())
elif domain == 'PiSquare':
result.update(smooth_pisquare())
else:
print('Invalid domain for smooth:', domain)
assert False
elif problem == 'Singular':
if domain == 'UnitSquare':
result.update(singular_square())
elif domain == 'LShape':
result.update(singular_lshape())
else:
print('Invalid domain for singular:', domain)
assert False
elif problem == 'Dirichlet':
result['g-linform'] = lambda elems: np.array(
[elem.h_t * elem.h_x for elem in elems])
result['g'] = lambda t, xy: 1
elif problem == 'MildSingular':
result['g-linform'] = lambda elems: np.array([
1 / 3 * elem.h_x *
(elem.time_interval[1]**3 - elem.time_interval[0]**3)
for elem in elems
])
result['g'] = lambda t, xy: t**2
return result