Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 694:1157375c678a feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add oposit jacobian to penalty term
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 19 Sep 2017 09:47:14 +0200
parents f235284e2eb1
children ba0d31ce4121
comparison
equal deleted inserted replaced
693:a5d4e3ced9a5 694:1157375c678a
98 obj.e_e = kr(e_r_u,I_v); 98 obj.e_e = kr(e_r_u,I_v);
99 obj.e_s = kr(I_u,e_l_v); 99 obj.e_s = kr(I_u,e_l_v);
100 obj.e_n = kr(I_u,e_r_v); 100 obj.e_n = kr(I_u,e_r_v);
101 obj.du_w = kr(d1_l_u,I_v); 101 obj.du_w = kr(d1_l_u,I_v);
102 obj.dv_w = (obj.e_w'*obj.Dv)'; 102 obj.dv_w = (obj.e_w'*obj.Dv)';
103 obj.du_e = kr(d1_r_u,I_v); 103 obj.du_e = kr(d1_r_u,I_v);
104 obj.dv_e = (obj.e_e'*obj.Dv)'; 104 obj.dv_e = (obj.e_e'*obj.Dv)';
105 obj.du_s = (obj.e_s'*obj.Du)'; 105 obj.du_s = (obj.e_s'*obj.Du)';
106 obj.dv_s = kr(I_u,d1_l_v); 106 obj.dv_s = kr(I_u,d1_l_v);
107 obj.du_n = (obj.e_n'*obj.Du)'; 107 obj.du_n = (obj.e_n'*obj.Du)';
108 obj.dv_n = kr(I_u,d1_r_v); 108 obj.dv_n = kr(I_u,d1_r_v);
225 end 225 end
226 226
227 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 227 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
228 % u denotes the solution in the own domain 228 % u denotes the solution in the own domain
229 % v denotes the solution in the neighbour domain 229 % v denotes the solution in the neighbour domain
230 [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); 230 [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,Ji_u] = obj.get_boundary_ops(boundary);
231 [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); 231 [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,Ji_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
232 232
233 a_n_u = @(t) spdiag(coeff_n_u(t)); 233 a_n_u = @(t) spdiag(coeff_n_u(t));
234 a_t_u = @(t) spdiag(coeff_t_u(t)); 234 a_t_u = @(t) spdiag(coeff_t_u(t));
235 a_n_v = @(t) spdiag(coeff_n_v(t)); 235 a_n_v = @(t) spdiag(coeff_n_v(t));
236 a_t_v = @(t) spdiag(coeff_t_v(t)); 236 a_t_v = @(t) spdiag(coeff_t_v(t));
248 penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u); 248 penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u);
249 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig ); 249 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig );
250 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) ); 250 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) );
251 251
252 252
253 closure =@(t) sqrt(obj.Ji)*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(obj.Ji); 253 closure =@(t) sqrt(Ji_u)*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(Ji_u);
254 penalty =@(t) sqrt(obj.Ji)*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(obj.Ji); 254 penalty =@(t) sqrt(Ji_v)*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(Ji_v);
255 end 255 end
256 256
257 257
258 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) 258 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)
259 259
316 case {'s','n'} 316 case {'s','n'}
317 halfnorm_inv_n = obj.Hiv; 317 halfnorm_inv_n = obj.Hiv;
318 halfnorm_inv_t = obj.Hiu; 318 halfnorm_inv_t = obj.Hiu;
319 halfnorm_t = obj.Hu; 319 halfnorm_t = obj.Hu;
320 end 320 end
321 Ji = obj.Ji;
321 end 322 end
322 323
323 function N = size(obj) 324 function N = size(obj)
324 N = prod(obj.m); 325 N = prod(obj.m);
325 end 326 end