Mercurial > repos > public > sbplib
annotate +scheme/Schrodinger2dCurve.m @ 494:c6d03f951a7f feature/quantumTriangles
Meged with default
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Thu, 23 Feb 2017 09:32:30 +0100 |
| parents | 6b8297f66c91 |
| children | 4905446f165e |
| rev | line source |
|---|---|
| 490 | 1 classdef Schrodinger2dCurve < scheme.Scheme |
| 2 properties | |
| 3 m % Number of points in each direction, possibly a vector | |
| 4 h % Grid spacing | |
| 5 | |
| 6 grid | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
7 xm, ym |
| 490 | 8 |
| 9 order % Order accuracy for the approximation | |
| 10 | |
| 11 D % non-stabalized scheme operator | |
| 12 M % Derivative norm | |
| 13 H % Discrete norm | |
| 14 Hi | |
| 15 H_u, H_v % Norms in the x and y directions | |
| 16 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 17 Hi_u, Hi_v | |
| 18 Hiu, Hiv | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
19 D1_v, D1_u |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
20 D2_v, D2_u |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
21 Du, Dv |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
22 x,y |
| 490 | 23 |
| 24 e_w, e_e, e_s, e_n | |
| 25 du_w, dv_w | |
| 26 du_e, dv_e | |
| 27 du_s, dv_s | |
| 28 du_n, dv_n | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
29 g_1, g_2 |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
30 c |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
31 a11, a12, a22 |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
32 m_tot, m_u, m_v |
| 490 | 33 p,p_tau |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
34 Ji |
| 490 | 35 end |
| 36 | |
| 37 methods | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
38 function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau) |
| 490 | 39 default_arg('opSet',@sbp.D2Variable); |
| 40 default_arg('c', 1); | |
| 41 | |
| 42 obj.p=p; | |
| 43 obj.p_tau=p_tau; | |
| 44 obj.c=1; | |
| 45 | |
| 46 m = g.size(); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
47 obj.m_u = m(1); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
48 obj.m_v = m(2); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
49 obj.m_tot = g.N(); |
| 490 | 50 |
| 51 h = g.scaling(); | |
| 52 h_u = h(1); | |
| 53 h_v = h(2); | |
| 54 | |
| 55 % Operators | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
56 ops_u = opSet(obj.m_u, {0, 1}, order); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
57 ops_v = opSet(obj.m_v, {0, 1}, order); |
| 490 | 58 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
59 I_u = speye(obj.m_u); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
60 I_v = speye(obj.m_v); |
| 490 | 61 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
62 obj.D1_u = ops_u.D1; |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
63 obj.D2_u = ops_u.D2; |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
64 |
| 490 | 65 H_u = ops_u.H; |
| 66 Hi_u = ops_u.HI; | |
| 67 e_l_u = ops_u.e_l; | |
| 68 e_r_u = ops_u.e_r; | |
| 69 d1_l_u = ops_u.d1_l; | |
| 70 d1_r_u = ops_u.d1_r; | |
| 71 | |
| 72 obj.D1_v = ops_v.D1; | |
| 73 obj.D2_v = ops_v.D2; | |
| 74 H_v = ops_v.H; | |
| 75 Hi_v = ops_v.HI; | |
| 76 e_l_v = ops_v.e_l; | |
| 77 e_r_v = ops_v.e_r; | |
| 78 d1_l_v = ops_v.d1_l; | |
| 79 d1_r_v = ops_v.d1_r; | |
| 80 | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
81 obj.Du = kr(obj.D1_u,I_v); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
82 obj.Dv = kr(I_u,obj.D1_v); |
| 490 | 83 |
| 84 obj.H = kr(H_u,H_v); | |
| 85 obj.Hi = kr(Hi_u,Hi_v); | |
| 86 obj.Hu = kr(H_u,I_v); | |
| 87 obj.Hv = kr(I_u,H_v); | |
| 88 obj.Hiu = kr(Hi_u,I_v); | |
| 89 obj.Hiv = kr(I_u,Hi_v); | |
| 90 | |
| 91 obj.e_w = kr(e_l_u,I_v); | |
| 92 obj.e_e = kr(e_r_u,I_v); | |
| 93 obj.e_s = kr(I_u,e_l_v); | |
| 94 obj.e_n = kr(I_u,e_r_v); | |
| 95 obj.du_w = kr(d1_l_u,I_v); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
96 obj.dv_w = (obj.e_w'*obj.Dv)'; |
| 490 | 97 obj.du_e = kr(d1_r_u,I_v); |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
98 obj.dv_e = (obj.e_e'*obj.Dv)'; |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
99 obj.du_s = (obj.e_s'*obj.Du)'; |
| 490 | 100 obj.dv_s = kr(I_u,d1_l_v); |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
101 obj.du_n = (obj.e_n'*obj.Du)'; |
| 490 | 102 obj.dv_n = kr(I_u,d1_r_v); |
| 103 | |
| 104 % obj.x_u = x_u; | |
| 105 % obj.x_v = x_v; | |
| 106 % obj.y_u = y_u; | |
| 107 % obj.y_v = y_v; | |
| 108 | |
| 109 obj.m = m; | |
| 110 obj.h = [h_u h_v]; | |
| 111 obj.order = order; | |
| 112 obj.grid = g; | |
| 113 | |
| 114 | |
| 115 end | |
| 116 | |
| 117 | |
| 118 function [D ]= d_fun(obj,t) | |
| 119 % Metric derivatives | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
120 ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t)); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
121 ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3(t),obj.p_tau.p4(t)); |
| 490 | 122 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
123 lcoords=points(obj.grid); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
124 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
125 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
126 x = reshape(obj.xm,obj.m_tot,1); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
127 y = reshape(obj.ym,obj.m_tot,1); |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
128 obj.x=x; |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
129 obj.y=y; |
| 490 | 130 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
131 x_tau = reshape(x_tau,obj.m_tot,1); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
132 y_tau = reshape(y_tau,obj.m_tot,1); |
| 490 | 133 |
| 134 x_u = obj.Du*x; | |
| 135 x_v = obj.Dv*x; | |
| 136 y_u = obj.Du*y; | |
| 137 y_v = obj.Dv*y; | |
| 138 | |
| 139 J = x_u.*y_v - x_v.*y_u; | |
| 140 a11 = 1./J.* (x_v.^2 + y_v.^2); | |
| 141 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); | |
| 142 a22 = 1./J .* (x_u.^2 + y_u.^2); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
143 |
| 493 | 144 obj.a11 = a11; |
| 145 obj.a12 = a12; | |
| 146 obj.a22 = a22; | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
147 |
| 490 | 148 % Assemble full operators |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
149 L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); |
| 490 | 150 Duv = obj.Du*L_12*obj.Dv; |
| 151 Dvu = obj.Dv*L_12*obj.Du; | |
| 152 | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
153 Duu = sparse(obj.m_tot); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
154 Dvv = sparse(obj.m_tot); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
155 ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); |
| 490 | 156 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
157 |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
158 for i = 1:obj.m_v |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
159 D = obj.D2_u(a11(ind(:,i))); |
| 490 | 160 p = ind(:,i); |
| 161 Duu(p,p) = D; | |
| 162 end | |
| 163 | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
164 for i = 1:obj.m_u |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
165 D = obj.D2_v(a22(ind(i,:))); |
| 490 | 166 p = ind(i,:); |
| 167 Dvv(p,p) = D; | |
| 168 end | |
| 169 | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
170 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
171 obj.Ji=Ji; |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
172 obj.g_1 = x_v.*y_tau-x_tau.*y_v; |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
173 obj.g_2 = x_tau.*y_u - y_tau.*x_u; |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
174 |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
175 b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot); |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
176 b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot); |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
177 |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
178 b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot); |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
179 b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot); |
| 490 | 180 |
| 181 %Add the flux splitting | |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
182 % D = Ji*(-b1*obj.Du -b2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
183 D = Ji*(-1/2*(b1*obj.Du-b1_u+obj.Du*b1) - 1/2*(b2*obj.Dv - b2_v +obj.Dv*b2) + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); |
| 490 | 184 |
| 185 % obj.gamm_u = h_u*ops_u.borrowing.M.d1; | |
| 186 % obj.gamm_v = h_v*ops_v.borrowing.M.d1; | |
| 187 | |
| 188 end | |
| 189 | |
| 190 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 191 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 192 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 193 % type is a string specifying the type of boundary condition if there are several. | |
| 194 % data is a function returning the data that should be applied at the boundary. | |
| 195 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 196 % neighbour_boundary is a string specifying which boundary to interface to. | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
197 function [closure, penalty] = boundary_condition(obj, boundary) |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
198 [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
199 |
| 490 | 200 a_t = spdiag(coeff_t); |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
201 F = (s * d_n' + s * a_t*d_t')'; |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
202 tau1 = 1; |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
203 a = spdiag(g); |
| 490 | 204 tau2 = (-1*s*a - abs(a))/4; |
| 205 | |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
206 penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e; |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
207 penalty_parameter_2 = halfnorm_inv_n*e*tau2; |
| 490 | 208 |
|
492
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
209 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' + obj.Ji* penalty_parameter_2*e'; |
|
6b95a894cbd7
fixed a bug in the metric coefficients but something is wrong at the boundaries
Ylva Rydin <ylva.rydin@telia.com>
parents:
491
diff
changeset
|
210 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1*e'+ obj.Ji*penalty_parameter_2*e'; |
| 490 | 211 |
| 212 end | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
213 |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
214 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) |
|
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
215 end |
| 490 | 216 |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
217 function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) |
| 490 | 218 |
| 219 % gridMatrix = zeros(obj.m(2),obj.m(1)); | |
| 220 % gridMatrix(:) = 1:numel(gridMatrix); | |
| 221 | |
| 222 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); | |
| 223 | |
| 224 switch boundary | |
| 225 case 'w' | |
| 226 e = obj.e_w; | |
| 227 d_n = obj.du_w; | |
| 228 d_t = obj.dv_w; | |
| 229 s = -1; | |
| 230 | |
| 231 I = ind(1,:); | |
| 232 coeff_t = obj.a12(I); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
233 g=obj.g_1(I); |
| 490 | 234 case 'e' |
| 235 e = obj.e_e; | |
| 236 d_n = obj.du_e; | |
| 237 d_t = obj.dv_e; | |
| 238 s = 1; | |
| 239 | |
| 240 I = ind(end,:); | |
| 241 coeff_t = obj.a12(I); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
242 g=obj.g_1(I); |
| 490 | 243 case 's' |
| 244 e = obj.e_s; | |
| 245 d_n = obj.dv_s; | |
| 246 d_t = obj.du_s; | |
| 247 s = -1; | |
| 248 | |
| 249 I = ind(:,1)'; | |
| 250 coeff_t = obj.a12(I); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
251 g=obj.g_2(I); |
| 490 | 252 case 'n' |
| 253 e = obj.e_n; | |
| 254 d_n = obj.dv_n; | |
| 255 d_t = obj.du_n; | |
| 256 s = 1; | |
| 257 | |
| 258 I = ind(:,end)'; | |
| 259 coeff_t = obj.a12(I); | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
260 g=obj.g_2(I); |
| 490 | 261 otherwise |
| 262 error('No such boundary: boundary = %s',boundary); | |
| 263 end | |
| 264 | |
| 265 switch boundary | |
| 266 case {'w','e'} | |
| 267 halfnorm_inv_n = obj.Hiu; | |
| 268 halfnorm_inv_t = obj.Hiv; | |
| 269 halfnorm_t = obj.Hv; | |
|
491
26125cfefe11
Schrodinger 2d runs without exploding with fixed coordinates
Ylva Rydin <ylva.rydin@telia.com>
parents:
490
diff
changeset
|
270 |
| 490 | 271 case {'s','n'} |
| 272 halfnorm_inv_n = obj.Hiv; | |
| 273 halfnorm_inv_t = obj.Hiu; | |
| 274 halfnorm_t = obj.Hu; | |
| 275 end | |
| 276 end | |
| 277 | |
| 278 function N = size(obj) | |
| 279 N = prod(obj.m); | |
| 280 end | |
| 281 | |
| 282 | |
| 283 end | |
| 284 end |