Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 970:23d9ca6755be feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryQuadrature in Elastic2dVariable. Rename get_boundary_operator -> getBoundaryOperator. Add operators in getBoundaryOperator, with full size so that they work with multiblock.DiffOp.getBoundaryOperator.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 25 Dec 2018 07:23:38 +0100
parents a4ad90b37998
children 104f0af001e0
comparison
equal deleted inserted replaced
969:adae8063ea2f 970:23d9ca6755be
365 comp = bc{1}; 365 comp = bc{1};
366 type = bc{2}; 366 type = bc{2};
367 367
368 % j is the coordinate direction of the boundary 368 % j is the coordinate direction of the boundary
369 j = obj.get_boundary_number(boundary); 369 j = obj.get_boundary_number(boundary);
370 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 370 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
371 371
372 372
373 E = obj.E; 373 E = obj.E;
374 Hi = obj.Hi; 374 Hi = obj.Hi;
375 LAMBDA = obj.LAMBDA; 375 LAMBDA = obj.LAMBDA;
387 switch type 387 switch type
388 388
389 % Dirichlet boundary condition 389 % Dirichlet boundary condition
390 case {'D','d','dirichlet','Dirichlet'} 390 case {'D','d','dirichlet','Dirichlet'}
391 391
392 alpha = obj.get_boundary_operator('alpha', boundary); 392 alpha = obj.getBoundaryOperator('alpha', boundary);
393 393
394 % Loop over components that Dirichlet penalties end up on 394 % Loop over components that Dirichlet penalties end up on
395 for i = 1:dim 395 for i = 1:dim
396 C = transpose(T{k,i}); 396 C = transpose(T{k,i});
397 A = -e*transpose(alpha{i,k}); 397 A = -e*transpose(alpha{i,k});
441 % j is the coordinate direction of the boundary 441 % j is the coordinate direction of the boundary
442 j = obj.get_boundary_number(boundary); 442 j = obj.get_boundary_number(boundary);
443 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary); 443 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
444 444
445 % Get boundary operators 445 % Get boundary operators
446 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 446 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
447 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary); 447 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary);
448 448
449 % Operators and quantities that correspond to the own domain only 449 % Operators and quantities that correspond to the own domain only
450 Hi = obj.Hi; 450 Hi = obj.Hi;
451 RHOi = obj.RHOi; 451 RHOi = obj.RHOi;
452 dim = obj.dim; 452 dim = obj.dim;
534 end 534 end
535 end 535 end
536 536
537 % Returns the boundary operator op for the boundary specified by the string boundary. 537 % Returns the boundary operator op for the boundary specified by the string boundary.
538 % op: may be a cell array of strings 538 % op: may be a cell array of strings
539 function [varargout] = get_boundary_operator(obj, op, boundary) 539 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
540 function [varargout] = getBoundaryOperator(obj, op, boundary)
540 541
541 switch boundary 542 switch boundary
542 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 543 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
543 j = 1; 544 j = 1;
544 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 545 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
558 case {'w','W','west','West','s','S','south','South'} 559 case {'w','W','west','West','s','S','south','South'}
559 varargout{k} = obj.e_l{j}; 560 varargout{k} = obj.e_l{j};
560 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 561 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
561 varargout{k} = obj.e_r{j}; 562 varargout{k} = obj.e_r{j};
562 end 563 end
564
565 case 'e_tot'
566 e = obj.getBoundaryOperator('e', boundary);
567 I_dim = speye(obj.dim, obj.dim);
568 varargout{k} = kron(e, I_dim);
569
563 case 'd' 570 case 'd'
564 switch boundary 571 switch boundary
565 case {'w','W','west','West','s','S','south','South'} 572 case {'w','W','west','West','s','S','south','South'}
566 varargout{k} = obj.d1_l{j}; 573 varargout{k} = obj.d1_l{j};
567 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 574 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
568 varargout{k} = obj.d1_r{j}; 575 varargout{k} = obj.d1_r{j};
569 end 576 end
577
570 case 'T' 578 case 'T'
571 switch boundary 579 switch boundary
572 case {'w','W','west','West','s','S','south','South'} 580 case {'w','W','west','West','s','S','south','South'}
573 varargout{k} = obj.T_l{j}; 581 varargout{k} = obj.T_l{j};
574 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 582 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
575 varargout{k} = obj.T_r{j}; 583 varargout{k} = obj.T_r{j};
576 end 584 end
585
577 case 'tau' 586 case 'tau'
578 switch boundary 587 switch boundary
579 case {'w','W','west','West','s','S','south','South'} 588 case {'w','W','west','West','s','S','south','South'}
580 varargout{k} = obj.tau_l{j}; 589 varargout{k} = obj.tau_l{j};
581 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 590 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
582 varargout{k} = obj.tau_r{j}; 591 varargout{k} = obj.tau_r{j};
583 end 592 end
593
594 case 'tau_tot'
595 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
596
597 I_dim = speye(obj.dim, obj.dim);
598 e_tot = kron(e, I_dim);
599 E = obj.E;
600 tau_tot = (e_tot'*E{1}*e*tau{1}')';
601 for i = 2:obj.dim
602 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
603 end
604 varargout{k} = tau_tot;
605
584 case 'H' 606 case 'H'
585 varargout{k} = obj.H_boundary{j}; 607 varargout{k} = obj.H_boundary{j};
608
586 case 'alpha' 609 case 'alpha'
587 % alpha = alpha(i,j) is the penalty strength for displacement BC. 610 % alpha = alpha(i,j) is the penalty strength for displacement BC.
588 switch boundary 611 e = obj.getBoundaryOperator('e', boundary);
589 case {'w','W','west','West','s','S','south','South'}
590 e = obj.e_l{j};
591 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
592 e = obj.e_r{j};
593 end
594 612
595 tuning = 1.2; 613 tuning = 1.2;
596 LAMBDA = obj.LAMBDA; 614 LAMBDA = obj.LAMBDA;
597 MU = obj.MU; 615 MU = obj.MU;
598 616
617 alpha{i,l} = d(i,l)*alpha_func(i,j)*e; 635 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
618 end 636 end
619 end 637 end
620 638
621 varargout{k} = alpha; 639 varargout{k} = alpha;
640
641 case 'alpha_tot'
642 % alpha = alpha(i,j) is the penalty strength for displacement BC.
643 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
644 E = obj.E;
645 [m, n] = size(alpha{1,1});
646 alpha_tot = sparse(m*obj.dim, n*obj.dim);
647 for i = 1:obj.dim
648 for l = 1:obj.dim
649 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
650 end
651 end
652 varargout{k} = alpha_tot;
653
622 otherwise 654 otherwise
623 error(['No such operator: operator = ' op{k}]); 655 error(['No such operator: operator = ' op{k}]);
624 end 656 end
625 end 657 end
626 658
627 end 659 end
628 660
661 function H = getBoundaryQuadrature(obj, boundary)
662 switch boundary
663 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
664 j = 1;
665 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
666 j = 2;
667 otherwise
668 error('No such boundary: boundary = %s',boundary);
669 end
670 H = obj.H_boundary{j};
671 I_dim = speye(obj.dim, obj.dim);
672 H = kron(H, I_dim);
673 end
674
629 function N = size(obj) 675 function N = size(obj)
630 N = obj.dim*prod(obj.m); 676 N = obj.dim*prod(obj.m);
631 end 677 end
632 end 678 end
633 end 679 end