Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.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 19d821ddc108
children 60c875c18de3
comparison
equal deleted inserted replaced
1087:867307f4d80f 1108:5ec23b9bf360
427 % u denotes the solution in the own domain 427 % u denotes the solution in the own domain
428 % v denotes the solution in the neighbour domain 428 % v denotes the solution in the neighbour domain
429 % Operators without subscripts are from the own domain. 429 % Operators without subscripts are from the own domain.
430 430
431 % Get boundary operators 431 % Get boundary operators
432 [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary); 432 e = obj.getBoundaryOperator('e_tot', boundary);
433 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary); 433 tau = obj.getBoundaryOperator('tau_tot', boundary);
434
435 e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary);
436 tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary);
437
434 H_gamma = obj.getBoundaryQuadrature(boundary); 438 H_gamma = obj.getBoundaryQuadrature(boundary);
435 439
436 % Operators and quantities that correspond to the own domain only 440 % Operators and quantities that correspond to the own domain only
437 Hi = obj.Hi_kron; 441 Hi = obj.Hi_kron;
438 RHOi = obj.RHOi_kron; 442 RHOi = obj.RHOi_kron;
456 error('Non-conforming interfaces not implemented yet.'); 460 error('Non-conforming interfaces not implemented yet.');
457 end 461 end
458 462
459 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 463 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
460 function [j, nj] = get_boundary_number(obj, boundary) 464 function [j, nj] = get_boundary_number(obj, boundary)
465 assertIsMember(boundary, {'w', 'e', 's', 'n'})
461 466
462 switch boundary 467 switch boundary
463 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 468 case {'w', 'e'}
464 j = 1; 469 j = 1;
465 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 470 case {'s', 'n'}
466 j = 2; 471 j = 2;
467 otherwise
468 error('No such boundary: boundary = %s',boundary);
469 end 472 end
470 473
471 switch boundary 474 switch boundary
472 case {'w','W','west','West','s','S','south','South'} 475 case {'w', 's'}
473 nj = -1; 476 nj = -1;
474 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 477 case {'e', 'n'}
475 nj = 1; 478 nj = 1;
476 end 479 end
477 end 480 end
478 481
479 % Returns the boundary operator op for the boundary specified by the string boundary. 482 % Returns the boundary operator op for the boundary specified by the string boundary.
480 % op: may be a cell array of strings 483 % op -- string
481 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() 484 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
482 function [varargout] = getBoundaryOperator(obj, op, boundary) 485 function [varargout] = getBoundaryOperator(obj, op, boundary)
486 assertIsMember(boundary, {'w', 'e', 's', 'n'})
487 assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'})
483 488
484 switch boundary 489 switch boundary
485 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 490 case {'w', 'e'}
486 j = 1; 491 j = 1;
487 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 492 case {'s', 'n'}
488 j = 2; 493 j = 2;
489 otherwise 494 end
490 error('No such boundary: boundary = %s',boundary); 495
491 end 496 switch op
492 497 case 'e'
493 if ~iscell(op) 498 switch boundary
494 op = {op}; 499 case {'w', 's'}
495 end 500 o = obj.e_l{j};
496 501 case {'e', 'n'}
497 for k = 1:length(op) 502 o = obj.e_r{j};
498 switch op{k} 503 end
499 case 'e' 504
500 switch boundary 505 case 'e_tot'
501 case {'w','W','west','West','s','S','south','South'} 506 e = obj.getBoundaryOperator('e', boundary);
502 varargout{k} = obj.e_l{j}; 507 I_dim = speye(obj.dim, obj.dim);
503 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 508 o = kron(e, I_dim);
504 varargout{k} = obj.e_r{j}; 509
510 case 'd'
511 switch boundary
512 case {'w', 's'}
513 o = obj.d1_l{j};
514 case {'e', 'n'}
515 o = obj.d1_r{j};
516 end
517
518 case 'T'
519 switch boundary
520 case {'w', 's'}
521 o = obj.T_l{j};
522 case {'e', 'n'}
523 o = obj.T_r{j};
524 end
525
526 case 'tau'
527 switch boundary
528 case {'w', 's'}
529 o = obj.tau_l{j};
530 case {'e', 'n'}
531 o = obj.tau_r{j};
532 end
533
534 case 'tau_tot'
535 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
536
537 I_dim = speye(obj.dim, obj.dim);
538 e_tot = kron(e, I_dim);
539 E = obj.E;
540 tau_tot = (e_tot'*E{1}*e*tau{1}')';
541 for i = 2:obj.dim
542 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
543 end
544 o = tau_tot;
545
546 case 'H'
547 o = obj.H_boundary{j};
548
549 case 'alpha'
550 % alpha = alpha(i,j) is the penalty strength for displacement BC.
551 e = obj.getBoundaryOperator('e', boundary);
552
553 LAMBDA = obj.LAMBDA;
554 MU = obj.MU;
555
556 dim = obj.dim;
557 theta_R = obj.theta_R{j};
558 theta_H = obj.theta_H{j};
559 theta_M = obj.theta_M{j};
560
561 a_lambda = dim/theta_H + 1/theta_R;
562 a_mu_i = 2/theta_M;
563 a_mu_ij = 2/theta_H + 1/theta_R;
564
565 d = @kroneckerDelta; % Kronecker delta
566 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
567 alpha = cell(obj.dim, obj.dim);
568
569 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
570 + d(i,j)* a_mu_i*MU ...
571 + db(i,j)*a_mu_ij*MU;
572 for i = 1:obj.dim
573 for l = 1:obj.dim
574 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
505 end 575 end
506 576 end
507 case 'e_tot' 577
508 e = obj.getBoundaryOperator('e', boundary); 578 o = alpha;
509 I_dim = speye(obj.dim, obj.dim); 579
510 varargout{k} = kron(e, I_dim); 580 case 'alpha_tot'
511 581 % alpha = alpha(i,j) is the penalty strength for displacement BC.
512 case 'd' 582 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
513 switch boundary 583 E = obj.E;
514 case {'w','W','west','West','s','S','south','South'} 584 [m, n] = size(alpha{1,1});
515 varargout{k} = obj.d1_l{j}; 585 alpha_tot = sparse(m*obj.dim, n*obj.dim);
516 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 586 for i = 1:obj.dim
517 varargout{k} = obj.d1_r{j}; 587 for l = 1:obj.dim
588 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
518 end 589 end
519 590 end
520 case 'T' 591 o = alpha_tot;
521 switch boundary 592 end
522 case {'w','W','west','West','s','S','south','South'} 593
523 varargout{k} = obj.T_l{j}; 594 end
524 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 595
525 varargout{k} = obj.T_r{j}; 596 % Returns square boundary quadrature matrix, of dimension
526 end 597 % corresponding to the number of boundary points
527 598 %
528 case 'tau' 599 % boundary -- string
529 switch boundary
530 case {'w','W','west','West','s','S','south','South'}
531 varargout{k} = obj.tau_l{j};
532 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
533 varargout{k} = obj.tau_r{j};
534 end
535
536 case 'tau_tot'
537 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
538
539 I_dim = speye(obj.dim, obj.dim);
540 e_tot = kron(e, I_dim);
541 E = obj.E;
542 tau_tot = (e_tot'*E{1}*e*tau{1}')';
543 for i = 2:obj.dim
544 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
545 end
546 varargout{k} = tau_tot;
547
548 case 'H'
549 varargout{k} = obj.H_boundary{j};
550
551 case 'alpha'
552 % alpha = alpha(i,j) is the penalty strength for displacement BC.
553 e = obj.getBoundaryOperator('e', boundary);
554
555 LAMBDA = obj.LAMBDA;
556 MU = obj.MU;
557
558 dim = obj.dim;
559 theta_R = obj.theta_R{j};
560 theta_H = obj.theta_H{j};
561 theta_M = obj.theta_M{j};
562
563 a_lambda = dim/theta_H + 1/theta_R;
564 a_mu_i = 2/theta_M;
565 a_mu_ij = 2/theta_H + 1/theta_R;
566
567 d = @kroneckerDelta; % Kronecker delta
568 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
569 alpha = cell(obj.dim, obj.dim);
570
571 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
572 + d(i,j)* a_mu_i*MU ...
573 + db(i,j)*a_mu_ij*MU;
574 for i = 1:obj.dim
575 for l = 1:obj.dim
576 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
577 end
578 end
579
580 varargout{k} = alpha;
581
582 case 'alpha_tot'
583 % alpha = alpha(i,j) is the penalty strength for displacement BC.
584 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
585 E = obj.E;
586 [m, n] = size(alpha{1,1});
587 alpha_tot = sparse(m*obj.dim, n*obj.dim);
588 for i = 1:obj.dim
589 for l = 1:obj.dim
590 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
591 end
592 end
593 varargout{k} = alpha_tot;
594
595 otherwise
596 error(['No such operator: operator = ' op{k}]);
597 end
598 end
599
600 end
601
602 function H = getBoundaryQuadrature(obj, boundary) 600 function H = getBoundaryQuadrature(obj, boundary)
601 assertIsMember(boundary, {'w', 'e', 's', 'n'})
602
603 switch boundary 603 switch boundary
604 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 604 case {'w','e'}
605 j = 1; 605 j = 1;
606 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 606 case {'s','n'}
607 j = 2; 607 j = 2;
608 otherwise
609 error('No such boundary: boundary = %s',boundary);
610 end 608 end
611 H = obj.H_boundary{j}; 609 H = obj.H_boundary{j};
612 I_dim = speye(obj.dim, obj.dim); 610 I_dim = speye(obj.dim, obj.dim);
613 H = kron(H, I_dim); 611 H = kron(H, I_dim);
614 end 612 end