Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1062:e512714fb890 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/getBoundaryOp
author Martin Almquist <malmquist@stanford.edu>
date Mon, 14 Jan 2019 18:14:44 -0800
parents 78db023a7fe3
children 8d73fcdb07a5 1341c6cea9c1
comparison
equal deleted inserted replaced
988:a72038b1f709 1062:e512714fb890
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 end
441
442 % Returns square boundary quadrature matrix, of dimension
443 % corresponding to the number of boundary points
383 % 444 %
384 % I -- the indices of the boundary points in the grid matrix 445 % boundary -- string
385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) 446 function H_b = getBoundaryQuadrature(obj, boundary)
386 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
387 447
388 switch boundary 448 switch boundary
389 case 'w' 449 case 'w'
390 e = obj.e_w;
391 d = obj.d_w;
392 H_b = obj.H_w; 450 H_b = obj.H_w;
451 case 'e'
452 H_b = obj.H_e;
453 case 's'
454 H_b = obj.H_s;
455 case 'n'
456 H_b = obj.H_n;
457 otherwise
458 error('No such boundary: boundary = %s',boundary);
459 end
460 end
461
462 % Returns the indices of the boundary points in the grid matrix
463 % boundary -- string
464 function I = getBoundaryIndices(obj, boundary)
465 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
466 switch boundary
467 case 'w'
393 I = ind(1,:); 468 I = ind(1,:);
394 case 'e' 469 case 'e'
395 e = obj.e_e;
396 d = obj.d_e;
397 H_b = obj.H_e;
398 I = ind(end,:); 470 I = ind(end,:);
399 case 's' 471 case 's'
400 e = obj.e_s;
401 d = obj.d_s;
402 H_b = obj.H_s;
403 I = ind(:,1)'; 472 I = ind(:,1)';
404 case 'n' 473 case 'n'
405 e = obj.e_n;
406 d = obj.d_n;
407 H_b = obj.H_n;
408 I = ind(:,end)'; 474 I = ind(:,end)';
409 otherwise 475 otherwise
410 error('No such boundary: boundary = %s',boundary); 476 error('No such boundary: boundary = %s',boundary);
411 end 477 end
412 478 end
479
480 % Returns borrowing constant gamma
481 % boundary -- string
482 function gamm = getBoundaryBorrowing(obj, boundary)
413 switch boundary 483 switch boundary
414 case {'w','e'} 484 case {'w','e'}
415 gamm = obj.gamm_u; 485 gamm = obj.gamm_u;
416 case {'s','n'} 486 case {'s','n'}
417 gamm = obj.gamm_v; 487 gamm = obj.gamm_v;
488 otherwise
489 error('No such boundary: boundary = %s',boundary);
418 end 490 end
419 end 491 end
420 492
421 function N = size(obj) 493 function N = size(obj)
422 N = prod(obj.m); 494 N = prod(obj.m);