Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 561:3a13916f8ff0 feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Switch to using boundary ops for normal derivative
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Aug 2017 13:10:57 +0200
parents 6132c52bf923
children 11d8d6ccbcd7
comparison
equal deleted inserted replaced
560:6132c52bf923 561:3a13916f8ff0
221 % neighbour_boundary is a string specifying which boundary to interface to. 221 % neighbour_boundary is a string specifying which boundary to interface to.
222 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 222 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
223 default_arg('type','neumann'); 223 default_arg('type','neumann');
224 default_arg('parameter', []); 224 default_arg('parameter', []);
225 225
226 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); 226 [e, d, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
227 switch type 227 switch type
228 % Dirichlet boundary condition 228 % Dirichlet boundary condition
229 case {'D','d','dirichlet'} 229 case {'D','d','dirichlet'}
230 % v denotes the solution in the neighbour domain
231 tuning = 1.2; 230 tuning = 1.2;
232 % tuning = 20.2; 231 % tuning = 20.2;
233 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); 232 [e, F, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
234
235 A_n = spdiag(a_n);
236 A_t = spdiag(a_t);
237
238 F = s*(A_n*d_n' + A_t*d_t')';
239 233
240 u = obj; 234 u = obj;
241 235
242 b1 = gamm*u.lambda./u.a11.^2; 236 b1 = gamm*u.lambda./u.a11.^2;
243 b2 = gamm*u.lambda./u.a22.^2; 237 b2 = gamm*u.lambda./u.a22.^2;
252 penalty = -obj.Ji*obj.a * penalty_parameter_1; 246 penalty = -obj.Ji*obj.a * penalty_parameter_1;
253 247
254 248
255 % Neumann boundary condition 249 % Neumann boundary condition
256 case {'N','n','neumann'} 250 case {'N','n','neumann'}
257 A_n = spdiag(a_n);
258 A_t = spdiag(a_t);
259 d = s*(A_n * d_n' + A_t*d_t')';
260
261 tau1 = -1; 251 tau1 = -1;
262 tau2 = 0; 252 tau2 = 0;
263 tau = s*obj.a*obj.Ji*(tau1*e + tau2*d); 253 tau = s*obj.a*obj.Ji*(tau1*e + tau2*d);
264 254
265 closure = halfnorm_inv*tau*d'; 255 closure = halfnorm_inv*tau*d';
267 257
268 % Characteristic boundary condition 258 % Characteristic boundary condition
269 case {'characteristic', 'char', 'c'} 259 case {'characteristic', 'char', 'c'}
270 default_arg('parameter', 1); 260 default_arg('parameter', 1);
271 beta = parameter; 261 beta = parameter;
272
273 A_n = spdiag(a_n);
274 A_t = spdiag(a_t);
275 d = s*(A_n * d_n' + A_t*d_t')'; % outward facing normal derivative
276 262
277 tau = -obj.a * 1/beta*obj.Ji*e; 263 tau = -obj.a * 1/beta*obj.Ji*e;
278 264
279 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; 265 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
280 closure{2} = halfnorm_inv*tau*beta*d'; 266 closure{2} = halfnorm_inv*tau*beta*d';
289 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 275 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
290 % u denotes the solution in the own domain 276 % u denotes the solution in the own domain
291 % v denotes the solution in the neighbour domain 277 % v denotes the solution in the neighbour domain
292 tuning = 1.2; 278 tuning = 1.2;
293 % tuning = 20.2; 279 % tuning = 20.2;
294 [e_u, d_n_u, d_t_u, a_n_u, a_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); 280 [e_u, F_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
295 [e_v, d_n_v, d_t_v, a_n_v, a_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); 281 [e_v, F_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);
296
297 A_n_u = spdiag(a_n_u);
298 A_t_u = spdiag(a_t_u);
299 A_n_v = spdiag(a_n_v);
300 A_t_v = spdiag(a_t_v);
301
302 F_u = s_u*(A_n_u * d_n_u' + A_t_u*d_t_u')';
303 F_v = s_v*(A_n_v * d_n_v' + A_t_v*d_t_v')';
304 282
305 u = obj; 283 u = obj;
306 v = neighbour_scheme; 284 v = neighbour_scheme;
307 285
308 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 286 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
325 303
326 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 304 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
327 % The right boundary is considered the positive boundary 305 % The right boundary is considered the positive boundary
328 % 306 %
329 % I -- the indecies of the boundary points in the grid matrix 307 % I -- the indecies of the boundary points in the grid matrix
330 function [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary) 308 function [e, d, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
331 309
332 % gridMatrix = zeros(obj.m(2),obj.m(1)); 310 % gridMatrix = zeros(obj.m(2),obj.m(1));
333 % gridMatrix(:) = 1:numel(gridMatrix); 311 % gridMatrix(:) = 1:numel(gridMatrix);
334 312
335 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 313 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
336 314
337 switch boundary 315 switch boundary
338 case 'w' 316 case 'w'
339 e = obj.e_w; 317 e = obj.e_w;
340 d_n = obj.du_w; 318 d = obj.d_w;
341 d_t = obj.dv_w;
342 s = -1; 319 s = -1;
343 320
344 I = ind(1,:); 321 I = ind(1,:);
345 a_n = obj.a11(I);
346 a_t = obj.a12(I);
347 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); 322 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
348 case 'e' 323 case 'e'
349 e = obj.e_e; 324 e = obj.e_e;
350 d_n = obj.du_e; 325 d = obj.d_e;
351 d_t = obj.dv_e;
352 s = 1; 326 s = 1;
353 327
354 I = ind(end,:); 328 I = ind(end,:);
355 a_n = obj.a11(I);
356 a_t = obj.a12(I);
357 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); 329 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
358 case 's' 330 case 's'
359 e = obj.e_s; 331 e = obj.e_s;
360 d_n = obj.dv_s; 332 d = obj.d_s;
361 d_t = obj.du_s;
362 s = -1; 333 s = -1;
363 334
364 I = ind(:,1)'; 335 I = ind(:,1)';
365 a_n = obj.a22(I);
366 a_t = obj.a12(I);
367 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); 336 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
368 case 'n' 337 case 'n'
369 e = obj.e_n; 338 e = obj.e_n;
370 d_n = obj.dv_n; 339 d = obj.d_n;
371 d_t = obj.du_n;
372 s = 1; 340 s = 1;
373 341
374 I = ind(:,end)'; 342 I = ind(:,end)';
375 a_n = obj.a22(I);
376 a_t = obj.a12(I);
377 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); 343 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
378 otherwise 344 otherwise
379 error('No such boundary: boundary = %s',boundary); 345 error('No such boundary: boundary = %s',boundary);
380 end 346 end
381 347