Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 499:f1465e6aeb26 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Commit before changing branch
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 06 Mar 2017 17:23:52 +0100
parents 324c927d8b1d
children 83734c26b8e3
comparison
equal deleted inserted replaced
498:324c927d8b1d 499:f1465e6aeb26
196 % type is a string specifying the type of boundary condition if there are several. 196 % type is a string specifying the type of boundary condition if there are several.
197 % data is a function returning the data that should be applied at the boundary. 197 % data is a function returning the data that should be applied at the boundary.
198 % neighbour_scheme is an instance of Scheme that should be interfaced to. 198 % neighbour_scheme is an instance of Scheme that should be interfaced to.
199 % neighbour_boundary is a string specifying which boundary to interface to. 199 % neighbour_boundary is a string specifying which boundary to interface to.
200 function [closure, penalty] = boundary_condition(obj, boundary,~) 200 function [closure, penalty] = boundary_condition(obj, boundary,~)
201 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g,I,Ji] = obj.get_boundary_ops(boundary); 201 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
202 202
203 a_t = @(t) spdiag(coeff_t(t)); 203 a_t = @(t) spdiag(coeff_t(t));
204 a_n = @(t) spdiag(coeff_n(t)); 204 a_n = @(t) spdiag(coeff_n(t));
205 Ji = spdiag(Ji); 205
206 206
207 F = @(t)(s * (d_n * Ji)' + s * a_t(t) *d_t')'; 207 F = @(t)(s * a_n(t)*d_n' + s * a_t(t) *d_t')';
208 tau1 = 1; 208 tau1 = 1;
209 a = @(t)spdiag(g(t)); 209 a = @(t)spdiag(g(t));
210 tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4; 210 tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4;
211 211
212 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e; 212 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
218 end 218 end
219 219
220 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 220 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
221 % u denotes the solution in the own domain 221 % u denotes the solution in the own domain
222 % v denotes the solution in the neighbour domain 222 % v denotes the solution in the neighbour domain
223 [e_u, d_n_u, d_t_u, coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u, I_u,JI_u] = obj.get_boundary_ops(boundary); 223 [e_u, d_n_u, d_t_u, coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u] = obj.get_boundary_ops(boundary);
224 [e_v, d_n_v, d_t_v, coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v, I_v,JI_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 224 [e_v, d_n_v, d_t_v, coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
225 225
226 a_n_u = @(t) spdiag(coeff_n_u(t)); 226 a_n_u = @(t) spdiag(coeff_n_u(t));
227 a_t_u = @(t) spdiag(coeff_t_u(t)); 227 a_t_u = @(t) spdiag(coeff_t_u(t));
228 a_n_v = @(t) spdiag(coeff_n_v(t)); 228 a_n_v = @(t) spdiag(coeff_n_v(t));
229 a_t_v = @(t) spdiag(coeff_t_v(t)); 229 a_t_v = @(t) spdiag(coeff_t_v(t));
230 230
231 Ji_u = spdiag(JI_u); 231
232 Ji_v = spdiag(JI_v); 232 F_u = @(t)(a_n_u(t)*d_n_u' + a_t_u(t)*d_t_u')';
233 233 F_v = @(t)(a_n_v(t)*d_n_v' + a_t_v(t)*d_t_v')';
234 F_u = @(t)((d_n_u*Ji_u)' + a_t_u(t)*d_t_u')';
235 F_v = @(t)((d_n_v*Ji_v)' + a_t_v(t)*d_t_v')';
236 234
237 a = @(t)spdiag(gamm_u(t)); 235 a = @(t)spdiag(gamm_u(t));
238 236
239 tau = s_u*1*1i/2; 237 tau = s_u*1*1i/2;
240 sig = -s_u*1*1i/2; 238 sig = -s_u*1*1i/2;
248 closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u'); 246 closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u');
249 penalty =@(t) obj.Ji*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v'); 247 penalty =@(t) obj.Ji*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v');
250 end 248 end
251 249
252 250
253 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I,Ji] = get_boundary_ops(obj, boundary) 251 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = get_boundary_ops(obj, boundary)
254 252
255 % gridMatrix = zeros(obj.m(2),obj.m(1)); 253 % gridMatrix = zeros(obj.m(2),obj.m(1));
256 % gridMatrix(:) = 1:numel(gridMatrix); 254 % gridMatrix(:) = 1:numel(gridMatrix);
257 255
258 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 256 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
284 d_t = obj.du_s; 282 d_t = obj.du_s;
285 s = -1; 283 s = -1;
286 284
287 I = ind(:,1)'; 285 I = ind(:,1)';
288 coeff_t = @(t)obj.a12(I); 286 coeff_t = @(t)obj.a12(I);
289 coeff_n = @(t)obj.a11(I); 287 coeff_n = @(t)obj.a22(I);
290 g = @(t)obj.g_2(I); 288 g = @(t)obj.g_2(I);
291 case 'n' 289 case 'n'
292 e = obj.e_n; 290 e = obj.e_n;
293 d_n = obj.dv_n; 291 d_n = obj.dv_n;
294 d_t = obj.du_n; 292 d_t = obj.du_n;
295 s = 1; 293 s = 1;
296 294
297 I = ind(:,end)'; 295 I = ind(:,end)';
298 coeff_t = @(t)obj.a12(I); 296 coeff_t = @(t)obj.a12(I);
299 coeff_n = @(t)obj.a11(I); 297 coeff_n = @(t)obj.a22(I);
300 g = @(t)obj.g_2(I); 298 g = @(t)obj.g_2(I);
301 otherwise 299 otherwise
302 error('No such boundary: boundary = %s',boundary); 300 error('No such boundary: boundary = %s',boundary);
303 end 301 end
304 302
311 case {'s','n'} 309 case {'s','n'}
312 halfnorm_inv_n = obj.Hiv; 310 halfnorm_inv_n = obj.Hiv;
313 halfnorm_inv_t = obj.Hiu; 311 halfnorm_inv_t = obj.Hiu;
314 halfnorm_t = obj.Hu; 312 halfnorm_t = obj.Hu;
315 end 313 end
316 fis = diag(obj.Ji);
317 Ji = fis(I);
318 end 314 end
319 315
320 function N = size(obj) 316 function N = size(obj)
321 N = prod(obj.m); 317 N = prod(obj.m);
322 end 318 end