Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1048:adbb80e60b10 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up Elastic2dVariable (partially), Utux, and Utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 22 Jan 2019 11:12:23 -0800
parents a2fcc4cf2298
children 19d821ddc108
comparison
equal deleted inserted replaced
1047:6bc55a773e7c 1048:adbb80e60b10
456 error('Non-conforming interfaces not implemented yet.'); 456 error('Non-conforming interfaces not implemented yet.');
457 end 457 end
458 458
459 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 459 % 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) 460 function [j, nj] = get_boundary_number(obj, boundary)
461 assertIsMember(boundary, {'w', 'e', 's', 'n'})
461 462
462 switch boundary 463 switch boundary
463 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 464 case {'w', 'e'}
464 j = 1; 465 j = 1;
465 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 466 case {'s', 'n'}
466 j = 2; 467 j = 2;
467 otherwise
468 error('No such boundary: boundary = %s',boundary);
469 end 468 end
470 469
471 switch boundary 470 switch boundary
472 case {'w','W','west','West','s','S','south','South'} 471 case {'w', 's'}
473 nj = -1; 472 nj = -1;
474 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 473 case {'e', 'n'}
475 nj = 1; 474 nj = 1;
476 end 475 end
477 end 476 end
478 477
479 % Returns the boundary operator op for the boundary specified by the string boundary. 478 % Returns the boundary operator op for the boundary specified by the string boundary.
480 % op: may be a cell array of strings 479 % op -- string
481 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() 480 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
482 function [varargout] = getBoundaryOperator(obj, op, boundary) 481 function [varargout] = getBoundaryOperator(obj, op, boundary)
482 assertIsMember(boundary, {'w', 'e', 's', 'n'})
483 assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'})
483 484
484 switch boundary 485 switch boundary
485 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 486 case {'w', 'e'}
486 j = 1; 487 j = 1;
487 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 488 case {'s', 'n'}
488 j = 2; 489 j = 2;
489 otherwise 490 end
490 error('No such boundary: boundary = %s',boundary); 491
491 end 492 switch op
492 493 case 'e'
493 if ~iscell(op) 494 switch boundary
494 op = {op}; 495 case {'w', 's'}
495 end 496 o = obj.e_l{j};
496 497 case {'e', 'n'}
497 for k = 1:length(op) 498 o = obj.e_r{j};
498 switch op{k} 499 end
499 case 'e' 500
500 switch boundary 501 case 'e_tot'
501 case {'w','W','west','West','s','S','south','South'} 502 e = obj.getBoundaryOperator('e', boundary);
502 varargout{k} = obj.e_l{j}; 503 I_dim = speye(obj.dim, obj.dim);
503 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 504 o = kron(e, I_dim);
504 varargout{k} = obj.e_r{j}; 505
506 case 'd'
507 switch boundary
508 case {'w', 's'}
509 o = obj.d1_l{j};
510 case {'e', 'n'}
511 o = obj.d1_r{j};
512 end
513
514 case 'T'
515 switch boundary
516 case {'w', 's'}
517 o = obj.T_l{j};
518 case {'e', 'n'}
519 o = obj.T_r{j};
520 end
521
522 case 'tau'
523 switch boundary
524 case {'w', 's'}
525 o = obj.tau_l{j};
526 case {'e', 'n'}
527 o = obj.tau_r{j};
528 end
529
530 case 'tau_tot'
531 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
532
533 I_dim = speye(obj.dim, obj.dim);
534 e_tot = kron(e, I_dim);
535 E = obj.E;
536 tau_tot = (e_tot'*E{1}*e*tau{1}')';
537 for i = 2:obj.dim
538 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
539 end
540 o = tau_tot;
541
542 case 'H'
543 o = obj.H_boundary{j};
544
545 case 'alpha'
546 % alpha = alpha(i,j) is the penalty strength for displacement BC.
547 e = obj.getBoundaryOperator('e', boundary);
548
549 LAMBDA = obj.LAMBDA;
550 MU = obj.MU;
551
552 dim = obj.dim;
553 theta_R = obj.theta_R{j};
554 theta_H = obj.theta_H{j};
555 theta_M = obj.theta_M{j};
556
557 a_lambda = dim/theta_H + 1/theta_R;
558 a_mu_i = 2/theta_M;
559 a_mu_ij = 2/theta_H + 1/theta_R;
560
561 d = @kroneckerDelta; % Kronecker delta
562 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
563 alpha = cell(obj.dim, obj.dim);
564
565 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
566 + d(i,j)* a_mu_i*MU ...
567 + db(i,j)*a_mu_ij*MU;
568 for i = 1:obj.dim
569 for l = 1:obj.dim
570 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
505 end 571 end
506 572 end
507 case 'e_tot' 573
508 e = obj.getBoundaryOperator('e', boundary); 574 o = alpha;
509 I_dim = speye(obj.dim, obj.dim); 575
510 varargout{k} = kron(e, I_dim); 576 case 'alpha_tot'
511 577 % alpha = alpha(i,j) is the penalty strength for displacement BC.
512 case 'd' 578 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
513 switch boundary 579 E = obj.E;
514 case {'w','W','west','West','s','S','south','South'} 580 [m, n] = size(alpha{1,1});
515 varargout{k} = obj.d1_l{j}; 581 alpha_tot = sparse(m*obj.dim, n*obj.dim);
516 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 582 for i = 1:obj.dim
517 varargout{k} = obj.d1_r{j}; 583 for l = 1:obj.dim
584 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
518 end 585 end
519 586 end
520 case 'T' 587 o = alpha_tot;
521 switch boundary 588 end
522 case {'w','W','west','West','s','S','south','South'} 589
523 varargout{k} = obj.T_l{j}; 590 end
524 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 591
525 varargout{k} = obj.T_r{j}; 592 % Returns square boundary quadrature matrix, of dimension
526 end 593 % corresponding to the number of boundary points
527 594 %
528 case 'tau' 595 % 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) 596 function H = getBoundaryQuadrature(obj, boundary)
597 assertIsMember(boundary, {'w', 'e', 's', 'n'})
598
603 switch boundary 599 switch boundary
604 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 600 case {'w','e'}
605 j = 1; 601 j = 1;
606 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 602 case {'s','n'}
607 j = 2; 603 j = 2;
608 otherwise
609 error('No such boundary: boundary = %s',boundary);
610 end 604 end
611 H = obj.H_boundary{j}; 605 H = obj.H_boundary{j};
612 I_dim = speye(obj.dim, obj.dim); 606 I_dim = speye(obj.dim, obj.dim);
613 H = kron(H, I_dim); 607 H = kron(H, I_dim);
614 end 608 end