-
Notifications
You must be signed in to change notification settings - Fork 0
/
h_adaptive.py
122 lines (110 loc) · 4.55 KB
/
h_adaptive.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
119
120
121
122
from new_fem import *
from dual_fem import solve_fem as solve_dual_fem
from helpers import new_norm, dual_norm, f_norm, State
from numpy import sqrt
from matplotlib import pyplot as plt
a = 0
b = 1
def draw(state):
# size, solution, dual_solution, nodes, norm, dual, fn, error
xs = []
ys = []
yds = []
nodes = state.nodes
un = state.solution
dual = state.dual_solution
size = state.size
nodes_0 = [0 for i in range(size)]
nodes_y = [un(nodes[i]) for i in range(size)]
nodes_yd = [dual(nodes[i]) for i in range(size)]
size = 1000
h = (b - a) / size
for i in range(size + 1):
x = a + h * i
xs.append(x)
ys.append(un(x))
yds.append(dual(x))
plt.plot(xs, ys, 'b', nodes, nodes_y, 'b^', xs, yds, 'g', nodes, nodes_yd, 'g^', nodes, nodes_0, 'rs')
h = nodes[-1] - nodes[0]
plt.xlim([nodes[0] - 0.05 * h, nodes[-1] + 0.05 * h])
plt.show()
# def h_adaptive_fem(m, sigma, f, alpha, _u, nodes, accuracy, states):
# solution = solve_fem(m=m, sigma=sigma, f=f, alpha=alpha, _u=_u, nodes=nodes)
# dual_solution = solve_dual_fem(m, sigma, f, alpha, _u, nodes)
# new_nodes = []
# straight_norms = []
# dual_norms = []
# f_norms = []
# errors = []
# size = len(nodes) - 1
# print(size)
# fn_full = f_norm(sigma, f, alpha, _u, nodes[0], nodes[-1], nodes[-1])
# for i in range(size):
# norm = new_norm(m, sigma, alpha, solution, nodes[i], nodes[i + 1], nodes[-1])
# dual = dual_norm(m, sigma, alpha, dual_solution, nodes[i], nodes[i + 1], nodes[-1])
# fn = f_norm(sigma, f, alpha, _u, nodes[i], nodes[i + 1], nodes[-1])
# straight_norms.append(norm)
# dual_norms.append(dual)
# new_nodes.append(nodes[i])
# f_norms.append(fn)
# error = abs(fn - (norm + dual))
# errors.append(error)
# # if error > accuracy:
# new_nodes.append((nodes[i] + nodes[i + 1]) / 2)
# new_nodes.append(nodes[-1])
# # fn = f_norm(sigma, f, alpha, nodes[0], nodes[-1], nodes[-1])
# print(fn)
# print(sum(f_norms))
# print(alpha * pow(_u, 2))
# error = sqrt(abs(fn_full - (sum(straight_norms) + sum(dual_norms))) / fn_full)
# new_state = State(size + 1, solution, dual_solution, nodes, straight_norms, dual_norms, fn_full, error, f_norms, errors)
# states.append(new_state)
# # draw(new_state)
# if len(nodes) < len(new_nodes)and len(nodes) < 150:
# return h_adaptive_fem(m, sigma, f, alpha, _u, new_nodes, accuracy, states)
# return states
def h_adaptive_fem(m, beta, sigma, f, alpha, _u, nodes, accuracy, states, even=False):
solution = solve_fem(m=m, beta=beta, sigma=sigma, f=f, alpha=alpha, _u=_u, nodes=nodes)
dual_solution = solve_dual_fem(m, beta, sigma, f, alpha, _u, nodes)
error_solution = solve_error(m, beta, sigma, f, alpha, _u, nodes, solution)
straight_norms = []
error_norms = []
size = len(nodes) - 1
print(size)
dual_norms = []
f_norms = []
norm_full = 0
error_norm_full = 0
errors_with_dual = []
for i in range(size):
e = new_norm(m, beta, sigma, alpha, error_solution, nodes[i], nodes[i + 1], nodes[-1])
error_norms.append(e)
error_norm_full += e
norm = new_norm(m, beta, sigma, alpha, solution, nodes[i], nodes[i + 1], nodes[-1])
straight_norms.append(norm)
norm_full += norm
dual = dual_norm(m, beta, sigma, f, alpha, _u, dual_solution, nodes[i], nodes[i + 1], nodes[-1])
dual_norms.append(dual)
fn = f_norm(sigma, f, alpha, _u, nodes[i], nodes[i + 1], nodes[-1])
f_norms.append(fn)
error = abs(fn - (norm + dual))
errors_with_dual.append(error)
new_nodes = []
errors = []
needs_repeat = False
for i in range(size):
new_nodes.append(nodes[i])
error = sqrt(size * error_norms[i] / (norm_full + error_norm_full))
errors.append(error)
if error > accuracy:
needs_repeat = True
if error > accuracy or even:
new_nodes.append((nodes[i] + nodes[i + 1]) / 2)
new_nodes.append(nodes[-1])
dual = dual_norm(m, beta, sigma, f, alpha, _u, dual_solution, nodes[0], nodes[-1], nodes[-1])
fn = f_norm(sigma, f, alpha, _u, nodes[0], nodes[-1], nodes[-1])
new_state = State(size + 1, solution, dual_solution, nodes, straight_norms, dual_norms, dual, f_norms, error_norms, fn)
states.append(new_state)
if needs_repeat and len(nodes) < 100:
return h_adaptive_fem(m, beta, sigma, f, alpha, _u, new_nodes, accuracy, states, even)
return states