Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 85:643bc513b8b8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 30 Nov 2015 11:41:34 +0100
parents 08fbc284718f
children 53fe4b64f65e
comparison
equal deleted inserted replaced
84:9d514669d372 85:643bc513b8b8
242 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 242 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
243 % u denotes the solution in the own domain 243 % u denotes the solution in the own domain
244 % v denotes the solution in the neighbour domain 244 % v denotes the solution in the neighbour domain
245 tuning = 1.2; 245 tuning = 1.2;
246 % tuning = 20.2; 246 % tuning = 20.2;
247 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary); 247 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
248 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 248 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
249 249
250 a_n_u = spdiag(coeff_n_u); 250 a_n_u = spdiag(coeff_n_u);
251 a_t_u = spdiag(coeff_t_u); 251 a_t_u = spdiag(coeff_t_u);
252 a_n_v = spdiag(coeff_n_v); 252 a_n_v = spdiag(coeff_n_v);
253 a_t_v = spdiag(coeff_t_v); 253 a_t_v = spdiag(coeff_t_v);
256 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; 256 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')';
257 257
258 u = obj; 258 u = obj;
259 v = neighbour_scheme; 259 v = neighbour_scheme;
260 260
261 b1_u = gamm_u*u.lambda./u.a11.^2; 261 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
262 b2_u = gamm_u*u.lambda./u.a22.^2; 262 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
263 b1_v = gamm_v*v.lambda./v.a11.^2; 263 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
264 b2_v = gamm_v*v.lambda./v.a22.^2; 264 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
265
266 265
267 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 266 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
268 tau = tuning * spdiag(tau(:)); 267 tau = tuning * spdiag(tau(:)); % Probably correct until here, see eq 27
269 sig1 = 1/2; 268 sig1 = 1/2;
270 sig2 = -1/2; 269 sig2 = -1/2;
271 270
272 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; 271 % penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; %% This is what is in the paper, but there is an error in dimensions.
272 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); %% Random guess at a fix, should check theory for this.
273 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; 273 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
274 274
275 275
276 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); 276 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
277 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); 277 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v');
278 end 278 end
279 279
280 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 280 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
281 % The right boundary is considered the positive boundary 281 % The right boundary is considered the positive boundary
282 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,a_n, a_t] = get_boundary_ops(obj,boundary) 282 %
283 % I -- the indecies of the boundary points in the grid matrix
284 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I] = get_boundary_ops(obj,boundary)
285
286 gridMatrix = zeros(obj.m(2),obj.m(1));
287 gridMatrix(:) = 1:numel(gridMatrix);
288
283 switch boundary 289 switch boundary
284 case 'w' 290 case 'w'
285 e = obj.e_w; 291 e = obj.e_w;
286 d_n = obj.du_w; 292 d_n = obj.du_w;
287 d_t = obj.dv_w; 293 d_t = obj.dv_w;
288 s = -1; 294 s = -1;
289 295
296 I = gridMatrix(:,1);
290 coeff_n = obj.a11(:,1); 297 coeff_n = obj.a11(:,1);
291 coeff_t = obj.a12(:,1); 298 coeff_t = obj.a12(:,1);
292 case 'e' 299 case 'e'
293 e = obj.e_e; 300 e = obj.e_e;
294 d_n = obj.du_e; 301 d_n = obj.du_e;
295 d_t = obj.dv_e; 302 d_t = obj.dv_e;
296 s = 1; 303 s = 1;
297 304
305 I = gridMatrix(:,end);
298 coeff_n = obj.a11(:,end); 306 coeff_n = obj.a11(:,end);
299 coeff_t = obj.a12(:,end); 307 coeff_t = obj.a12(:,end);
300 case 's' 308 case 's'
301 e = obj.e_s; 309 e = obj.e_s;
302 d_n = obj.dv_s; 310 d_n = obj.dv_s;
303 d_t = obj.du_s; 311 d_t = obj.du_s;
304 s = -1; 312 s = -1;
305 313
314 I = gridMatrix(1,:)';
306 coeff_n = obj.a22(1,:)'; 315 coeff_n = obj.a22(1,:)';
307 coeff_t = obj.a12(1,:)'; 316 coeff_t = obj.a12(1,:)';
308 case 'n' 317 case 'n'
309 e = obj.e_n; 318 e = obj.e_n;
310 d_n = obj.dv_n; 319 d_n = obj.dv_n;
311 d_t = obj.du_n; 320 d_t = obj.du_n;
312 s = 1; 321 s = 1;
313 322
323 I = gridMatrix(end,:)';
314 coeff_n = obj.a22(end,:)'; 324 coeff_n = obj.a22(end,:)';
315 coeff_t = obj.a12(end,:)'; 325 coeff_t = obj.a12(end,:)';
316 otherwise 326 otherwise
317 error('No such boundary: boundary = %s',boundary); 327 error('No such boundary: boundary = %s',boundary);
318 end 328 end