|
1 |
| -import numpy as np |
| 1 | +import numpy as np |
2 | 2 | from . import jacobi
|
3 | 3 |
|
4 | 4 | # The defalut configurations for the base Jacobi parameters.
|
5 |
| -alpha = (-0.5,-0.5) |
| 5 | +alpha = (-0.5, -0.5) |
6 | 6 |
|
7 |
| -def quadrature(Nmax,alpha=alpha,**kw): |
8 |
| - return jacobi.quadrature(Nmax,alpha[0],alpha[1],**kw) |
9 | 7 |
|
10 |
| -def trial_functions(Nmax,z,alpha=alpha): |
| 8 | +def quadrature(Nmax, alpha=alpha, **kw): |
| 9 | + return jacobi.quadrature(Nmax, alpha[0], alpha[1], **kw) |
11 | 10 |
|
12 |
| - init = 1/np.sqrt(jacobi.mass(alpha[0],alpha[1])) + 0*z |
13 |
| - return jacobi.recursion(Nmax,alpha[0],alpha[1],z,init) |
14 | 11 |
|
15 |
| -def operator(dimension,op,Nmax,k,ell,radii,pad=0,alpha=alpha): |
| 12 | +def trial_functions(Nmax, z, alpha=alpha): |
| 13 | + |
| 14 | + init = 1 / np.sqrt(jacobi.mass(alpha[0], alpha[1])) + 0 * z |
| 15 | + return jacobi.recursion(Nmax, alpha[0], alpha[1], z, init) |
| 16 | + |
| 17 | + |
| 18 | +def operator(dimension, op, Nmax, k, ell, radii, pad=0, alpha=alpha): |
16 | 19 | # Pad the matrices by a safe amount before outputting the correct size.
|
17 |
| - |
18 |
| - if radii[1] <= radii[0]: raise ValueError('Inner radius must be greater than outer radius.') |
19 |
| - |
20 |
| - gapwidth = radii[1] - radii[0] |
21 |
| - aspectratio = (radii[1] + radii[0])/gapwidth |
22 |
| - |
23 |
| - a, b, N = k+alpha[0], k+alpha[1], Nmax + pad |
24 |
| - |
| 20 | + |
| 21 | + if radii[1] <= radii[0]: |
| 22 | + raise ValueError('Inner radius must be greater than outer radius.') |
| 23 | + |
| 24 | + gapwidth = radii[1] - radii[0] |
| 25 | + aspectratio = (radii[1] + radii[0]) / gapwidth |
| 26 | + |
| 27 | + a, b, N = k + alpha[0], k + alpha[1], Nmax + pad |
| 28 | + |
25 | 29 | # zeros
|
26 |
| - if (op == '0'): return jacobi.operator('0',N,a,b) |
27 |
| - |
| 30 | + if op == '0': |
| 31 | + return jacobi.operator('0', N, a, b) |
| 32 | + |
28 | 33 | # identity
|
29 |
| - if (op == 'I'): return jacobi.operator('I',N,a,b) |
30 |
| - |
31 |
| - Z = aspectratio*jacobi.operator('I',N+2,a,b) |
32 |
| - Z += jacobi.operator('J',N+2,a,b) |
33 |
| - |
| 34 | + if op == 'I': |
| 35 | + return jacobi.operator('I', N, a, b) |
| 36 | + |
| 37 | + Z = aspectratio * jacobi.operator('I', N + 2, a, b) |
| 38 | + Z += jacobi.operator('J', N + 2, a, b) |
| 39 | + |
34 | 40 | # r multiplication
|
35 |
| - if op == 'R': return (gapwidth/2)*Z[:N+1,:N+1] |
36 |
| - |
37 |
| - E = jacobi.operator('A+',N+2,a,b+1) @ jacobi.operator('B+',N+2,a,b) |
38 |
| - |
| 41 | + if op == 'R': |
| 42 | + return (gapwidth / 2) * Z[: N + 1, : N + 1] |
| 43 | + |
| 44 | + E = jacobi.operator('A+', N + 2, a, b + 1) @ jacobi.operator('B+', N + 2, a, b) |
| 45 | + |
39 | 46 | # conversion
|
40 |
| - if op == 'E': return 0.5*(E @ Z)[:N+1,:N+1] |
41 |
| - |
42 |
| - D = jacobi.operator('D+',N+2,a,b) @ Z |
43 |
| - |
| 47 | + if op == 'E': |
| 48 | + return 0.5 * (E @ Z)[: N + 1, : N + 1] |
| 49 | + |
| 50 | + D = jacobi.operator('D+', N + 2, a, b) @ Z |
| 51 | + |
44 | 52 | # derivatives
|
45 |
| - if op == 'D+': return (D - (ell+k+1 )*E)[:N+1,:N+1]/gapwidth |
46 |
| - if op == 'D-': return (D + (ell-k+dimension-3)*E)[:N+1,:N+1]/gapwidth |
47 |
| - |
48 |
| - # restriction |
49 |
| - if op == 'r=Ri': return jacobi.operator('z=-1',N,a,b) |
50 |
| - if op == 'r=Ro': return jacobi.operator('z=+1',N,a,b) |
| 53 | + if op == 'D+': |
| 54 | + return (D - (ell + k + 1) * E)[: N + 1, : N + 1] / gapwidth |
| 55 | + if op == 'D-': |
| 56 | + return (D + (ell - k + dimension - 3) * E)[: N + 1, : N + 1] / gapwidth |
51 | 57 |
|
| 58 | + # restriction |
| 59 | + if op == 'r=Ri': |
| 60 | + return jacobi.operator('z=-1', N, a, b) |
| 61 | + if op == 'r=Ro': |
| 62 | + return jacobi.operator('z=+1', N, a, b) |
0 commit comments