-
Notifications
You must be signed in to change notification settings - Fork 0
/
April24_H1_optimal.py
130 lines (64 loc) · 2.23 KB
/
April24_H1_optimal.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
123
124
125
126
127
128
129
130
# In this file I have the codes needed to compute
# the mesh density function that minimizes the H1 norm
from numerical_integration import *
def interior_u_xx(Uc, Ub, Ua, xc, xb, xa):
"""
Uc : U_i+1
Ua : U_i-1
xc, xb, xa : corresponding grid points
"""
U0 = Ua # U_i-1
U1 = Ub # Ui
U2 = Uc # U_i+1
h1 = xb - xa
h2 = xc - xb
# compute the derivative u_xx
u_xx_top = h1 * Uc + h2 * Ua - (h1 + h2) * Ub
u_xx_bottom = (1/2)* (h1*h2**2 + h2 * h1**2)
u_xx = u_xx_top / u_xx_bottom
return u_xx
def endpoint_u_xx(Ua, Ub, Uc, xa, xb, xc):
"""
Uc = U_i+2 or U_i-2
Ub : U_i
Ua : Ui+1 or U_i-1
xc, xb, xa : corresponding grid points
"""
hb = xb - xa
hc = xc - xb
term_a = 2 / (hb**2 + hb*hc)
term_b = -2 / (hb*hc)
term_c = 2 / (hc**2 + hb*hc)
u_xx_a = term_a * Ua
u_xx_b = term_b * Ub
u_xx_c = term_c * Uc
u_xx = u_xx_a + u_xx_b + u_xx_c
return u_xx
# code to generate M values on a given grid
def M_calc_optimal_H1(U, grid):
"""
U : Vector of Approximations
grid : mesh
"""
u_dpr = [] # empty list to save M values
u_dpr.append(endpoint_u_xx(U[0][0], U[1][0], U[2][0], grid[0][0], grid[1][0], grid[2][0])) # add M value at first grid point using forward approx
for i in range(len(grid)-2): # add M values for interior nodes
# set values of U and y to pass to mesh density function
U1 = U[i+1][0]
U0 = U[i][0]
U2 = U[i+2][0]
x0 = grid[i][0]
x1 = grid[i+1][0]
x2 = grid[i+2][0]
val = interior_u_xx(U2, U1, U0, x2, x1, x0) # compute mesh density
u_dpr.append(val)
u_dpr.append(endpoint_u_xx(U[-3][0], U[-2][0], U[-1][0], grid[-3][0], grid[-2][0], grid[-1][0])) # add M value for last grid point
u_dpr = np.array([u_dpr]).T
# compute the desired integral
u_dpr_power = abs(u_dpr)**(2/3)
# int_u_dpr_power = trap_nonuni(u_dpr_power, grid)
int_u_dpr_power = int_constant(u_dpr_power, grid)
alpha_int = ((1 / (grid[-1][0] - grid[0][0])) * int_u_dpr_power )**(3)
# compute the mesh density function
rho = ( 1 + (1/alpha_int) * (abs(u_dpr))**2)**(1/3)
return rho