Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 498:324c927d8b1d feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
chnaged sbp interfacein 1d among many things
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 03 Mar 2017 16:19:41 +0100
parents 4905446f165e
children f1465e6aeb26
comparison
equal deleted inserted replaced
497:4905446f165e 498:324c927d8b1d
132 else 132 else
133 ti = parametrization.Ti.points(obj.p{1}(t),obj.p{2}(t),obj.p{3}(t),obj.p{4}(t)); 133 ti = parametrization.Ti.points(obj.p{1}(t),obj.p{2}(t),obj.p{3}(t),obj.p{4}(t));
134 ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t)); 134 ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t));
135 135
136 lcoords=points(obj.grid); 136 lcoords=points(obj.grid);
137 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 137 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
138 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 138 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
139 x = reshape(obj.xm,obj.m_tot,1); 139 x = reshape(obj.xm,obj.m_tot,1);
140 y = reshape(obj.ym,obj.m_tot,1); 140 y = reshape(obj.ym,obj.m_tot,1);
141 obj.x = x; 141 obj.x = x;
142 obj.y = y; 142 obj.y = y;
143 143
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] = 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,I,Ji] = obj.get_boundary_ops(boundary);
202 202
203 a_t = spdiag(coeff_t); 203 a_t = @(t) spdiag(coeff_t(t));
204 a_n = spdiag(coeff_n); 204 a_n = @(t) spdiag(coeff_n(t));
205 F = (s * a_n *d_n' + s * a_t*d_t')'; 205 Ji = spdiag(Ji);
206
207 F = @(t)(s * (d_n * Ji)' + s * a_t(t) *d_t')';
206 tau1 = 1; 208 tau1 = 1;
207 a = spdiag(g); 209 a = @(t)spdiag(g(t));
208 tau2 = (-1*s*a - abs(a))/4; 210 tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4;
209 211
210 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e; 212 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
211 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2; 213 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t);
212 214
213 closure = @(t) obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e'; 215 closure = @(t) obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e';
214 penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'+ obj.Ji*penalty_parameter_2(t)*e'; 216 penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'- obj.Ji*penalty_parameter_2(t)*e';
215 217
216 end 218 end
217 219
218 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 220 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
219 % u denotes the solution in the own domain 221 % u denotes the solution in the own domain
220 % v denotes the solution in the neighbour domain 222 % v denotes the solution in the neighbour domain
221 [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] = 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, I_u,JI_u] = obj.get_boundary_ops(boundary);
222 [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] = 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, I_v,JI_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
223 225
224 a_n_u = spdiag(coeff_n_u); 226 a_n_u = @(t) spdiag(coeff_n_u(t));
225 a_t_u = spdiag(coeff_t_u); 227 a_t_u = @(t) spdiag(coeff_t_u(t));
226 a_n_v = spdiag(coeff_n_v); 228 a_n_v = @(t) spdiag(coeff_n_v(t));
227 a_t_v = spdiag(coeff_t_v); 229 a_t_v = @(t) spdiag(coeff_t_v(t));
228 230
229 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; 231 Ji_u = spdiag(JI_u);
230 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; 232 Ji_v = spdiag(JI_v);
231 233
232 a = spdiag(gamm_u); 234 F_u = @(t)((d_n_u*Ji_u)' + a_t_u(t)*d_t_u')';
233 235 F_v = @(t)((d_n_v*Ji_v)' + a_t_v(t)*d_t_v')';
234 u = obj; 236
235 v = neighbour_scheme; 237 a = @(t)spdiag(gamm_u(t));
236 238
237 tau = -1/2*1i; 239 tau = s_u*1*1i/2;
238 sig = 1/2*1i; 240 sig = -s_u*1*1i/2;
239 gamm = (-1*s_u*a)/2; 241 gamm = @(t) (-s_u*a(t))/2;
240 242
241 penalty_parameter_1 =@(t) halfnorm_inv_u_n*(e_u*tau + sig*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); 243 penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u);
242 penalty_parameter_2 =@(t) halfnorm_inv_u_n * e_u * (-sig + gamm ); 244 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig );
243 245 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) );
244 closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u'); 246
245 penalty =@(t) obj.Ji*obj.c^2 * (-penalty_parameter_1(t)*e_v' + penalty_parameter_2(t)*F_v'); 247
246 end 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');
247 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');
248 250 end
249 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) 251
252
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)
250 254
251 % gridMatrix = zeros(obj.m(2),obj.m(1)); 255 % gridMatrix = zeros(obj.m(2),obj.m(1));
252 % gridMatrix(:) = 1:numel(gridMatrix); 256 % gridMatrix(:) = 1:numel(gridMatrix);
253 257
254 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 258 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
259 d_n = obj.du_w; 263 d_n = obj.du_w;
260 d_t = obj.dv_w; 264 d_t = obj.dv_w;
261 s = -1; 265 s = -1;
262 266
263 I = ind(1,:); 267 I = ind(1,:);
264 coeff_t = obj.a12(I); 268 coeff_t = @(t)obj.a12(I);
265 coeff_n = obj.a11(I); 269 coeff_n =@(t) obj.a11(I);
266 g=obj.g_1(I); 270 g = @(t)obj.g_1(I);
267 case 'e' 271 case 'e'
268 e = obj.e_e; 272 e = obj.e_e;
269 d_n = obj.du_e; 273 d_n = obj.du_e;
270 d_t = obj.dv_e; 274 d_t = obj.dv_e;
271 s = 1; 275 s = 1;
272 276
273 I = ind(end,:); 277 I = ind(end,:);
274 coeff_t = obj.a12(I); 278 coeff_t = @(t)obj.a12(I);
275 coeff_n = obj.a11(I); 279 coeff_n = @(t)obj.a11(I);
276 g=obj.g_1(I); 280 g = @(t)obj.g_1(I);
277 case 's' 281 case 's'
278 e = obj.e_s; 282 e = obj.e_s;
279 d_n = obj.dv_s; 283 d_n = obj.dv_s;
280 d_t = obj.du_s; 284 d_t = obj.du_s;
281 s = -1; 285 s = -1;
282 286
283 I = ind(:,1)'; 287 I = ind(:,1)';
284 coeff_t = obj.a12(I); 288 coeff_t = @(t)obj.a12(I);
285 coeff_n = obj.a11(I); 289 coeff_n = @(t)obj.a11(I);
286 g=obj.g_2(I); 290 g = @(t)obj.g_2(I);
287 case 'n' 291 case 'n'
288 e = obj.e_n; 292 e = obj.e_n;
289 d_n = obj.dv_n; 293 d_n = obj.dv_n;
290 d_t = obj.du_n; 294 d_t = obj.du_n;
291 s = 1; 295 s = 1;
292 296
293 I = ind(:,end)'; 297 I = ind(:,end)';
294 coeff_t = obj.a12(I); 298 coeff_t = @(t)obj.a12(I);
295 coeff_n = obj.a11(I); 299 coeff_n = @(t)obj.a11(I);
296 g=obj.g_2(I); 300 g = @(t)obj.g_2(I);
297 otherwise 301 otherwise
298 error('No such boundary: boundary = %s',boundary); 302 error('No such boundary: boundary = %s',boundary);
299 end 303 end
300 304
301 switch boundary 305 switch boundary
307 case {'s','n'} 311 case {'s','n'}
308 halfnorm_inv_n = obj.Hiv; 312 halfnorm_inv_n = obj.Hiv;
309 halfnorm_inv_t = obj.Hiu; 313 halfnorm_inv_t = obj.Hiu;
310 halfnorm_t = obj.Hu; 314 halfnorm_t = obj.Hu;
311 end 315 end
316 fis = diag(obj.Ji);
317 Ji = fis(I);
312 end 318 end
313 319
314 function N = size(obj) 320 function N = size(obj)
315 N = prod(obj.m); 321 N = prod(obj.m);
316 end 322 end
317 end 323 end
324 methods(Static)
325 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
326 % and bound_v of scheme schm_v.
327 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
328 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
329 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
330 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
331 end
332 end
318 end 333 end