Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1066:d64062bed5fb feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up getBoundarOperator/Quadrature as was done on feature/getBoundaryOp. Temporarily re-add lambda so that old interface method works.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 22 Jan 2019 12:53:13 -0800
parents c2bd7f15da48
children 9a858436f8fa
comparison
equal deleted inserted replaced
1065:c2bd7f15da48 1066:d64062bed5fb
42 42
43 % Borrowing constants 43 % Borrowing constants
44 theta_M_u, theta_M_v 44 theta_M_u, theta_M_v
45 theta_R_u, theta_R_v 45 theta_R_u, theta_R_v
46 theta_H_u, theta_H_v 46 theta_H_u, theta_H_v
47
48 % Temporary
49 lambda
50 gamm_u, gamm_v
47 51
48 end 52 end
49 53
50 methods 54 methods
51 % Implements a*div(b*grad(u)) as a SBP scheme 55 % Implements a*div(b*grad(u)) as a SBP scheme
137 141
138 J = x_u.*y_v - x_v.*y_u; 142 J = x_u.*y_v - x_v.*y_u;
139 a11 = 1./J .* (x_v.^2 + y_v.^2); 143 a11 = 1./J .* (x_v.^2 + y_v.^2);
140 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 144 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
141 a22 = 1./J .* (x_u.^2 + y_u.^2); 145 a22 = 1./J .* (x_u.^2 + y_u.^2);
142 % lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); 146 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
143 147
144 K = cell(dim, dim); 148 K = cell(dim, dim);
145 K{1,1} = spdiag(y_v./J); 149 K{1,1} = spdiag(y_v./J);
146 K{1,2} = spdiag(-y_u./J); 150 K{1,2} = spdiag(-y_u./J);
147 K{2,1} = spdiag(-x_v./J); 151 K{2,1} = spdiag(-x_v./J);
232 obj.a = a; 236 obj.a = a;
233 obj.b = b; 237 obj.b = b;
234 obj.a11 = a11; 238 obj.a11 = a11;
235 obj.a12 = a12; 239 obj.a12 = a12;
236 obj.a22 = a22; 240 obj.a22 = a22;
237 % obj.lambda = lambda;
238 obj.s_w = spdiag(s_w); 241 obj.s_w = spdiag(s_w);
239 obj.s_e = spdiag(s_e); 242 obj.s_e = spdiag(s_e);
240 obj.s_s = spdiag(s_s); 243 obj.s_s = spdiag(s_s);
241 obj.s_n = spdiag(s_n); 244 obj.s_n = spdiag(s_n);
242 245
246 obj.theta_R_u = h_u*ops_u.borrowing.R.delta_D; 249 obj.theta_R_u = h_u*ops_u.borrowing.R.delta_D;
247 obj.theta_R_v = h_v*ops_v.borrowing.R.delta_D; 250 obj.theta_R_v = h_v*ops_v.borrowing.R.delta_D;
248 251
249 obj.theta_H_u = h_u*ops_u.borrowing.H11; 252 obj.theta_H_u = h_u*ops_u.borrowing.H11;
250 obj.theta_H_v = h_v*ops_v.borrowing.H11; 253 obj.theta_H_v = h_v*ops_v.borrowing.H11;
254
255 % Temporary
256 obj.lambda = lambda;
257 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
258 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
251 end 259 end
252 260
253 261
254 % Closure functions return the opertors applied to the own doamin to close the boundary 262 % Closure functions return the opertors applied to the own doamin to close the boundary
255 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 263 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
260 % neighbour_boundary is a string specifying which boundary to interface to. 268 % neighbour_boundary is a string specifying which boundary to interface to.
261 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 269 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
262 default_arg('type','neumann'); 270 default_arg('type','neumann');
263 default_arg('parameter', []); 271 default_arg('parameter', []);
264 272
265 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); 273 e = obj.getBoundaryOperator('e', boundary);
274 d = obj.getBoundaryOperator('d', boundary);
266 H_b = obj.getBoundaryQuadrature(boundary); 275 H_b = obj.getBoundaryQuadrature(boundary);
267 s_b = obj.getBoundaryScaling(boundary); 276 s_b = obj.getBoundaryScaling(boundary);
268 [th_H, th_M, th_R] = obj.getBoundaryBorrowing(boundary); 277 [th_H, th_M, th_R] = obj.getBoundaryBorrowing(boundary);
269 m = obj.getBoundaryNumber(boundary); 278 m = obj.getBoundaryNumber(boundary);
270 279
320 % Fields: 329 % Fields:
321 % -- tuning: penalty strength, defaults to 1.2 330 % -- tuning: penalty strength, defaults to 1.2
322 % -- interpolation: type of interpolation, default 'none' 331 % -- interpolation: type of interpolation, default 'none'
323 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 332 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
324 333
325 error('Not implemented') 334 % error('Not implemented')
326 335
327 defaultType.tuning = 1.2; 336 defaultType.tuning = 1.2;
328 defaultType.interpolation = 'none'; 337 defaultType.interpolation = 'none';
329 default_struct('type', defaultType); 338 default_struct('type', defaultType);
330 339
341 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 350 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
342 tuning = type.tuning; 351 tuning = type.tuning;
343 352
344 % u denotes the solution in the own domain 353 % u denotes the solution in the own domain
345 % v denotes the solution in the neighbour domain 354 % v denotes the solution in the neighbour domain
346 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); 355 e_u = obj.getBoundaryOperator('e', boundary);
356 d_u = obj.getBoundaryOperator('d', boundary);
347 H_b_u = obj.getBoundaryQuadrature(boundary); 357 H_b_u = obj.getBoundaryQuadrature(boundary);
348 I_u = obj.getBoundaryIndices(boundary); 358 I_u = obj.getBoundaryIndices(boundary);
349 gamm_u = obj.getBoundaryBorrowing(boundary); 359 gamm_u = obj.getBoundaryBorrowing(boundary);
350 360
351 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); 361 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
362 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
352 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); 363 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
353 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); 364 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
354 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); 365 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
355 366
356 u = obj; 367 u = obj;
386 tuning = type.tuning; 397 tuning = type.tuning;
387 398
388 399
389 % u denotes the solution in the own domain 400 % u denotes the solution in the own domain
390 % v denotes the solution in the neighbour domain 401 % v denotes the solution in the neighbour domain
391 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); 402 e_u = obj.getBoundaryOperator('e', boundary);
403 d_u = obj.getBoundaryOperator('d', boundary);
392 H_b_u = obj.getBoundaryQuadrature(boundary); 404 H_b_u = obj.getBoundaryQuadrature(boundary);
393 I_u = obj.getBoundaryIndices(boundary); 405 I_u = obj.getBoundaryIndices(boundary);
394 gamm_u = obj.getBoundaryBorrowing(boundary); 406 gamm_u = obj.getBoundaryBorrowing(boundary);
395 407
396 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); 408 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
409 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
397 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); 410 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
398 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); 411 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
399 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); 412 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
400 413
401 414
437 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; 450 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
438 451
439 end 452 end
440 453
441 % Returns the boundary operator op for the boundary specified by the string boundary. 454 % Returns the boundary operator op for the boundary specified by the string boundary.
442 % op -- string or a cell array of strings 455 % op -- string
443 % boundary -- string 456 % boundary -- string
444 function varargout = getBoundaryOperator(obj, op, boundary) 457 function o = getBoundaryOperator(obj, op, boundary)
445 458 assertIsMember(op, {'e', 'd'})
446 if ~iscell(op) 459 assertIsMember(boundary, {'w', 'e', 's', 'n'})
447 op = {op}; 460
448 end 461 o = obj.([op, '_', boundary]);
449
450 for i = 1:numel(op)
451 switch op{i}
452 case 'e'
453 switch boundary
454 case 'w'
455 e = obj.e_w;
456 case 'e'
457 e = obj.e_e;
458 case 's'
459 e = obj.e_s;
460 case 'n'
461 e = obj.e_n;
462 otherwise
463 error('No such boundary: boundary = %s',boundary);
464 end
465 varargout{i} = e;
466
467 case 'd'
468 switch boundary
469 case 'w'
470 d = obj.d_w;
471 case 'e'
472 d = obj.d_e;
473 case 's'
474 d = obj.d_s;
475 case 'n'
476 d = obj.d_n;
477 otherwise
478 error('No such boundary: boundary = %s',boundary);
479 end
480 varargout{i} = d;
481 end
482 end
483 end 462 end
484 463
485 % Returns square boundary quadrature matrix, of dimension 464 % Returns square boundary quadrature matrix, of dimension
486 % corresponding to the number of boundary points 465 % corresponding to the number of boundary points
487 % 466 %
488 % boundary -- string 467 % boundary -- string
489 function H_b = getBoundaryQuadrature(obj, boundary) 468 function H_b = getBoundaryQuadrature(obj, boundary)
490 469 assertIsMember(boundary, {'w', 'e', 's', 'n'})
491 switch boundary 470
492 case 'w' 471 H_b = obj.(['H_', boundary]);
493 H_b = obj.H_w;
494 case 'e'
495 H_b = obj.H_e;
496 case 's'
497 H_b = obj.H_s;
498 case 'n'
499 H_b = obj.H_n;
500 otherwise
501 error('No such boundary: boundary = %s',boundary);
502 end
503 end 472 end
504 473
505 % Returns square boundary quadrature scaling matrix, of dimension 474 % Returns square boundary quadrature scaling matrix, of dimension
506 % corresponding to the number of boundary points 475 % corresponding to the number of boundary points
507 % 476 %
508 % boundary -- string 477 % boundary -- string
509 function s_b = getBoundaryScaling(obj, boundary) 478 function s_b = getBoundaryScaling(obj, boundary)
510 479 assertIsMember(boundary, {'w', 'e', 's', 'n'})
511 switch boundary 480
512 case 'w' 481 s_b = obj.(['s_', boundary]);
513 s_b = obj.s_w;
514 case 'e'
515 s_b = obj.s_e;
516 case 's'
517 s_b = obj.s_s;
518 case 'n'
519 s_b = obj.s_n;
520 otherwise
521 error('No such boundary: boundary = %s',boundary);
522 end
523 end 482 end
524 483
525 % Returns the coordinate number corresponding to the boundary 484 % Returns the coordinate number corresponding to the boundary
526 % 485 %
527 % boundary -- string 486 % boundary -- string
528 function m = getBoundaryNumber(obj, boundary) 487 function m = getBoundaryNumber(obj, boundary)
488 assertIsMember(boundary, {'w', 'e', 's', 'n'})
529 489
530 switch boundary 490 switch boundary
531 case {'w', 'e'} 491 case {'w', 'e'}
532 m = 1; 492 m = 1;
533 case {'s', 'n'} 493 case {'s', 'n'}
534 m = 2; 494 m = 2;
535 otherwise
536 error('No such boundary: boundary = %s',boundary);
537 end 495 end
538 end 496 end
539 497
540 % Returns the indices of the boundary points in the grid matrix 498 % Returns the indices of the boundary points in the grid matrix
541 % boundary -- string 499 % boundary -- string