Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1108:5ec23b9bf360 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Wed, 10 Apr 2019 11:00:27 -0700
parents 867307f4d80f 84200bbae101
children 01d28cfafe7c
comparison
equal deleted inserted replaced
1087:867307f4d80f 1108:5ec23b9bf360
276 % neighbour_boundary is a string specifying which boundary to interface to. 276 % neighbour_boundary is a string specifying which boundary to interface to.
277 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 277 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
278 default_arg('type','neumann'); 278 default_arg('type','neumann');
279 default_arg('parameter', []); 279 default_arg('parameter', []);
280 280
281 e = obj.getBoundaryOperator('e', boundary); 281 e = obj.getBoundaryOperator('e', boundary);
282 d = obj.getBoundaryOperator('d', boundary); 282 d = obj.getBoundaryOperator('d', boundary);
283 H_b = obj.getBoundaryQuadrature(boundary); 283 H_b = obj.getBoundaryQuadrature(boundary);
284 s_b = obj.getBoundaryScaling(boundary); 284 s_b = obj.getBoundaryScaling(boundary);
285 [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary); 285 [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary);
286 m = obj.getBoundaryNumber(boundary); 286 m = obj.getBoundaryNumber(boundary);
287 287
288 K = obj.K; 288 K = obj.K;
289 J = obj.J; 289 J = obj.J;
290 Hi = obj.Hi; 290 Hi = obj.Hi;
291 a = obj.a; 291 a = obj.a;
356 % v denotes the solution in the neighbour domain 356 % v denotes the solution in the neighbour domain
357 u = obj; 357 u = obj;
358 v = neighbour_scheme; 358 v = neighbour_scheme;
359 359
360 % Boundary operators, u 360 % Boundary operators, u
361 e_u = u.getBoundaryOperator('e', boundary); 361 e_u = u.getBoundaryOperator('e', boundary);
362 d_u = u.getBoundaryOperator('d', boundary); 362 d_u = u.getBoundaryOperator('d', boundary);
363 gamm_u = u.getBoundaryBorrowing(boundary); 363 gamm_u = u.getBoundaryBorrowing(boundary);
364 s_b_u = u.getBoundaryScaling(boundary); 364 s_b_u = u.getBoundaryScaling(boundary);
365 [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary); 365 [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary);
366 m_u = u.getBoundaryNumber(boundary); 366 m_u = u.getBoundaryNumber(boundary);
367 367
368 % Coefficients, u 368 % Coefficients, u
369 K_u = u.K; 369 K_u = u.K;
370 J_u = u.J; 370 J_u = u.J;
371 b_b_u = e_u'*u.b*e_u; 371 b_b_u = e_u'*u.b*e_u;
372 372
373 % Boundary operators, v 373 % Boundary operators, v
374 e_v = v.getBoundaryOperator('e', neighbour_boundary); 374 e_v = v.getBoundaryOperator('e', neighbour_boundary);
375 d_v = v.getBoundaryOperator('d', neighbour_boundary); 375 d_v = v.getBoundaryOperator('d', neighbour_boundary);
376 gamm_v = v.getBoundaryBorrowing(neighbour_boundary); 376 gamm_v = v.getBoundaryBorrowing(neighbour_boundary);
377 s_b_v = v.getBoundaryScaling(neighbour_boundary); 377 s_b_v = v.getBoundaryScaling(neighbour_boundary);
378 [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary); 378 [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary);
379 m_v = v.getBoundaryNumber(neighbour_boundary); 379 m_v = v.getBoundaryNumber(neighbour_boundary);
380 380
381 % Coefficients, v 381 % Coefficients, v
382 K_v = v.K; 382 K_v = v.K;
429 tuning = type.tuning; 429 tuning = type.tuning;
430 430
431 431
432 % u denotes the solution in the own domain 432 % u denotes the solution in the own domain
433 % v denotes the solution in the neighbour domain 433 % v denotes the solution in the neighbour domain
434 e_u = obj.getBoundaryOperator('e', boundary); 434 e_u = obj.getBoundaryOperator('e', boundary);
435 d_u = obj.getBoundaryOperator('d', boundary); 435 d_u = obj.getBoundaryOperator('d', boundary);
436 H_b_u = obj.getBoundaryQuadrature(boundary); 436 H_b_u = obj.getBoundaryQuadrature(boundary);
437 I_u = obj.getBoundaryIndices(boundary); 437 I_u = obj.getBoundaryIndices(boundary);
438 gamm_u = obj.getBoundaryBorrowing(boundary); 438 gamm_u = obj.getBoundaryBorrowing(boundary);
439 439
440 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); 440 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
441 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary); 441 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
442 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); 442 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
443 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); 443 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
444 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); 444 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
445 445
446 446
447 % Find the number of grid points along the interface 447 % Find the number of grid points along the interface
448 m_u = size(e_u, 2); 448 m_u = size(e_u, 2);
528 end 528 end
529 529
530 % Returns the indices of the boundary points in the grid matrix 530 % Returns the indices of the boundary points in the grid matrix
531 % boundary -- string 531 % boundary -- string
532 function I = getBoundaryIndices(obj, boundary) 532 function I = getBoundaryIndices(obj, boundary)
533 assertIsMember(boundary, {'w', 'e', 's', 'n'})
534
533 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 535 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
534 switch boundary 536 switch boundary
535 case 'w' 537 case 'w'
536 I = ind(1,:); 538 I = ind(1,:);
537 case 'e' 539 case 'e'
538 I = ind(end,:); 540 I = ind(end,:);
539 case 's' 541 case 's'
540 I = ind(:,1)'; 542 I = ind(:,1)';
541 case 'n' 543 case 'n'
542 I = ind(:,end)'; 544 I = ind(:,end)';
543 otherwise
544 error('No such boundary: boundary = %s',boundary);
545 end 545 end
546 end 546 end
547 547
548 % Returns borrowing constant gamma 548 % Returns borrowing constant gamma
549 % boundary -- string 549 % boundary -- string
550 function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary) 550 function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary)
551 assertIsMember(boundary, {'w', 'e', 's', 'n'})
552
551 switch boundary 553 switch boundary
552 case {'w','e'} 554 case {'w','e'}
553 theta_H = obj.theta_H_u; 555 theta_H = obj.theta_H_u;
554 theta_M = obj.theta_M_u; 556 theta_M = obj.theta_M_u;
555 theta_R = obj.theta_R_u; 557 theta_R = obj.theta_R_u;
556 case {'s','n'} 558 case {'s','n'}
557 theta_H = obj.theta_H_v; 559 theta_H = obj.theta_H_v;
558 theta_M = obj.theta_M_v; 560 theta_M = obj.theta_M_v;
559 theta_R = obj.theta_R_v; 561 theta_R = obj.theta_R_v;
560 otherwise
561 error('No such boundary: boundary = %s',boundary);
562 end 562 end
563 end 563 end
564 564
565 function N = size(obj) 565 function N = size(obj)
566 N = prod(obj.m); 566 N = prod(obj.m);