Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 982:a0b3161e44f3 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add the following methods in LaplaceCurvilinear: getBoundaryOperator, getBoundaryQuadrature, getBoundaryBorrowing. Remove get_boundary_ops. Make interface and boundary_condition methods use the new methods.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 08 Jan 2019 11:51:24 +0100
parents 706d1c2b4199
children 78db023a7fe3
comparison
equal deleted inserted replaced
981:a2fcc4cf2298 982:a0b3161e44f3
236 % neighbour_boundary is a string specifying which boundary to interface to. 236 % neighbour_boundary is a string specifying which boundary to interface to.
237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
238 default_arg('type','neumann'); 238 default_arg('type','neumann');
239 default_arg('parameter', []); 239 default_arg('parameter', []);
240 240
241 [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); 241 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
242 H_b = obj.getBoundaryQuadrature(boundary);
243 gamm = obj.getBoundaryBorrowing(boundary);
242 switch type 244 switch type
243 % Dirichlet boundary condition 245 % Dirichlet boundary condition
244 case {'D','d','dirichlet'} 246 case {'D','d','dirichlet'}
245 tuning = 1.2; 247 tuning = 1.2;
246 % tuning = 20.2; 248 % tuning = 20.2;
296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 298 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
297 tuning = type.tuning; 299 tuning = type.tuning;
298 300
299 % u denotes the solution in the own domain 301 % u denotes the solution in the own domain
300 % v denotes the solution in the neighbour domain 302 % v denotes the solution in the neighbour domain
301 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 303 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
302 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 304 H_b_u = obj.getBoundaryQuadrature(boundary);
305 I_u = obj.getBoundaryIndices(boundary);
306 gamm_u = obj.getBoundaryBorrowing(boundary);
307
308 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
309 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
310 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
311 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
303 312
304 u = obj; 313 u = obj;
305 v = neighbour_scheme; 314 v = neighbour_scheme;
306 315
307 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 316 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
334 tuning = type.tuning; 343 tuning = type.tuning;
335 344
336 345
337 % u denotes the solution in the own domain 346 % u denotes the solution in the own domain
338 % v denotes the solution in the neighbour domain 347 % v denotes the solution in the neighbour domain
339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 348 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 349 H_b_u = obj.getBoundaryQuadrature(boundary);
350 I_u = obj.getBoundaryIndices(boundary);
351 gamm_u = obj.getBoundaryBorrowing(boundary);
352
353 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
354 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
355 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
356 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
357
341 358
342 % Find the number of grid points along the interface 359 % Find the number of grid points along the interface
343 m_u = size(e_u, 2); 360 m_u = size(e_u, 2);
344 m_v = size(e_v, 2); 361 m_v = size(e_v, 2);
345 362
376 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... 393 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
377 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; 394 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
378 395
379 end 396 end
380 397
381 % Returns the boundary ops and sign for the boundary specified by the string boundary. 398 % Returns the boundary operator op for the boundary specified by the string boundary.
382 % The right boundary is considered the positive boundary 399 % op -- string or a cell array of strings
400 % boundary -- string
401 function varargout = getBoundaryOperator(obj, op, boundary)
402
403 if ~iscell(op)
404 op = {op};
405 end
406
407 for i = 1:numel(op)
408 switch op{i}
409 case 'e'
410 switch boundary
411 case 'w'
412 e = obj.e_w;
413 case 'e'
414 e = obj.e_e;
415 case 's'
416 e = obj.e_s;
417 case 'n'
418 e = obj.e_n;
419 otherwise
420 error('No such boundary: boundary = %s',boundary);
421 end
422 varargout{i} = e;
423
424 case 'd'
425 switch boundary
426 case 'w'
427 d = obj.d_w;
428 case 'e'
429 d = obj.d_e;
430 case 's'
431 d = obj.d_s;
432 case 'n'
433 d = obj.d_n;
434 otherwise
435 error('No such boundary: boundary = %s',boundary);
436 end
437 varargout{i} = d;
438 end
439 end
440
441 end
442
443 % Returns square boundary quadrature matrix, of dimension
444 % corresponding to the number of boundary points
383 % 445 %
384 % I -- the indices of the boundary points in the grid matrix 446 % boundary -- string
385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) 447 function H_b = getBoundaryQuadrature(obj, boundary)
386 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
387 448
388 switch boundary 449 switch boundary
389 case 'w' 450 case 'w'
390 e = obj.e_w;
391 d = obj.d_w;
392 H_b = obj.H_w; 451 H_b = obj.H_w;
452 case 'e'
453 H_b = obj.H_e;
454 case 's'
455 H_b = obj.H_s;
456 case 'n'
457 H_b = obj.H_n;
458 otherwise
459 error('No such boundary: boundary = %s',boundary);
460 end
461 end
462
463 % Returns the indices of the boundary points in the grid matrix
464 % boundary -- string
465 function I = getBoundaryIndices(obj, boundary)
466 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
467 switch boundary
468 case 'w'
393 I = ind(1,:); 469 I = ind(1,:);
394 case 'e' 470 case 'e'
395 e = obj.e_e;
396 d = obj.d_e;
397 H_b = obj.H_e;
398 I = ind(end,:); 471 I = ind(end,:);
399 case 's' 472 case 's'
400 e = obj.e_s;
401 d = obj.d_s;
402 H_b = obj.H_s;
403 I = ind(:,1)'; 473 I = ind(:,1)';
404 case 'n' 474 case 'n'
405 e = obj.e_n;
406 d = obj.d_n;
407 H_b = obj.H_n;
408 I = ind(:,end)'; 475 I = ind(:,end)';
409 otherwise 476 otherwise
410 error('No such boundary: boundary = %s',boundary); 477 error('No such boundary: boundary = %s',boundary);
411 end 478 end
412 479 end
480
481 % Returns borrowing constant gamma
482 % boundary -- string
483 function gamm = getBoundaryBorrowing(obj, boundary)
413 switch boundary 484 switch boundary
414 case {'w','e'} 485 case {'w','e'}
415 gamm = obj.gamm_u; 486 gamm = obj.gamm_u;
416 case {'s','n'} 487 case {'s','n'}
417 gamm = obj.gamm_v; 488 gamm = obj.gamm_v;
489 otherwise
490 error('No such boundary: boundary = %s',boundary);
418 end 491 end
419 end 492 end
420 493
421 function N = size(obj) 494 function N = size(obj)
422 N = prod(obj.m); 495 N = prod(obj.m);