Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1057:ff274c7404cc feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 25 Jan 2019 15:52:21 -0800
parents b4fa176b4287
children 84933722ec0e
comparison
equal deleted inserted replaced
1056:b4fa176b4287 1057:ff274c7404cc
360 comp = obj.getComponent(comp, boundary); 360 comp = obj.getComponent(comp, boundary);
361 end 361 end
362 362
363 % j is the coordinate direction of the boundary 363 % j is the coordinate direction of the boundary
364 j = obj.get_boundary_number(boundary); 364 j = obj.get_boundary_number(boundary);
365 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); 365 e = obj.getBoundaryOperator('e', boundary);
366 T = obj.getBoundaryOperator('T', boundary);
367 tau = obj.getBoundaryOperator('tau', boundary);
368 H_gamma = obj.getBoundaryOperator('H', boundary);
366 369
367 E = obj.E; 370 E = obj.E;
368 Hi = obj.Hi; 371 Hi = obj.Hi;
369 LAMBDA = obj.LAMBDA; 372 LAMBDA = obj.LAMBDA;
370 MU = obj.MU; 373 MU = obj.MU;
431 % u denotes the solution in the own domain 434 % u denotes the solution in the own domain
432 % v denotes the solution in the neighbour domain 435 % v denotes the solution in the neighbour domain
433 % Operators without subscripts are from the own domain. 436 % Operators without subscripts are from the own domain.
434 437
435 % Get boundary operators 438 % Get boundary operators
436 [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary); 439 e = obj.getBoundaryOperator('e_tot', boundary);
437 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary); 440 tau = obj.getBoundaryOperator('tau_tot', boundary);
441
442 e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary);
443 tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary);
444
438 H_gamma = obj.getBoundaryQuadrature(boundary); 445 H_gamma = obj.getBoundaryQuadrature(boundary);
439 446
440 % Operators and quantities that correspond to the own domain only 447 % Operators and quantities that correspond to the own domain only
441 Hi = obj.Hi_kron; 448 Hi = obj.Hi_kron;
442 RHOi = obj.RHOi_kron; 449 RHOi = obj.RHOi_kron;
503 end 510 end
504 end 511 end
505 end 512 end
506 513
507 % Returns the boundary operator op for the boundary specified by the string boundary. 514 % Returns the boundary operator op for the boundary specified by the string boundary.
508 % op: may be a cell array of strings 515 % op -- string
509 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() 516 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
510 function [varargout] = getBoundaryOperator(obj, op, boundary) 517 function o = getBoundaryOperator(obj, op, boundary)
518 assertIsMember(boundary, {'w', 'e', 's', 'n'})
519 assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'})
511 520
512 switch boundary 521 switch boundary
513 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 522 case {'w', 'e'}
514 j = 1; 523 j = 1;
515 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 524 case {'s', 'n'}
516 j = 2; 525 j = 2;
517 otherwise 526 end
518 error('No such boundary: boundary = %s',boundary); 527
519 end 528 switch op
520 529 case 'e'
521 if ~iscell(op) 530 switch boundary
522 op = {op}; 531 case {'w', 's'}
523 end 532 o = obj.e_l{j};
524 533 case {'e', 'n'}
525 for k = 1:length(op) 534 o = obj.e_r{j};
526 switch op{k} 535 end
527 case 'e' 536
528 switch boundary 537 case 'e_tot'
529 case {'w','W','west','West','s','S','south','South'} 538 e = obj.getBoundaryOperator('e', boundary);
530 varargout{k} = obj.e_l{j}; 539 I_dim = speye(obj.dim, obj.dim);
531 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 540 o = kron(e, I_dim);
532 varargout{k} = obj.e_r{j}; 541
542 case 'd'
543 switch boundary
544 case {'w', 's'}
545 o = obj.d1_l{j};
546 case {'e', 'n'}
547 o = obj.d1_r{j};
548 end
549
550 case 'T'
551 switch boundary
552 case {'w', 's'}
553 o = obj.T_l{j};
554 case {'e', 'n'}
555 o = obj.T_r{j};
556 end
557
558 case 'tau'
559 switch boundary
560 case {'w', 's'}
561 o = obj.tau_l{j};
562 case {'e', 'n'}
563 o = obj.tau_r{j};
564 end
565
566 case 'tau_tot'
567 e = obj.getBoundaryOperator('e', boundary);
568 tau = obj.getBoundaryOperator('tau', boundary);
569
570 I_dim = speye(obj.dim, obj.dim);
571 e_tot = kron(e, I_dim);
572 E = obj.E;
573 tau_tot = (e_tot'*E{1}*e*tau{1}')';
574 for i = 2:obj.dim
575 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
576 end
577 o = tau_tot;
578
579 case 'H'
580 o = obj.H_boundary{j};
581
582 case 'alpha'
583 % alpha = alpha(i,j) is the penalty strength for displacement BC.
584 e = obj.getBoundaryOperator('e', boundary);
585
586 LAMBDA = obj.LAMBDA;
587 MU = obj.MU;
588
589 dim = obj.dim;
590 theta_R = obj.theta_R{j};
591 theta_H = obj.theta_H{j};
592 theta_M = obj.theta_M{j};
593
594 a_lambda = dim/theta_H + 1/theta_R;
595 a_mu_i = 2/theta_M;
596 a_mu_ij = 2/theta_H + 1/theta_R;
597
598 d = @kroneckerDelta; % Kronecker delta
599 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
600 alpha = cell(obj.dim, obj.dim);
601
602 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
603 + d(i,j)* a_mu_i*MU ...
604 + db(i,j)*a_mu_ij*MU;
605 for i = 1:obj.dim
606 for l = 1:obj.dim
607 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
533 end 608 end
534 609 end
535 case 'e_tot' 610
536 e = obj.getBoundaryOperator('e', boundary); 611 o = alpha;
537 I_dim = speye(obj.dim, obj.dim); 612
538 varargout{k} = kron(e, I_dim); 613 case 'alpha_tot'
539 614 % alpha = alpha(i,j) is the penalty strength for displacement BC.
540 case 'd' 615 e = obj.getBoundaryOperator('e', boundary);
541 switch boundary 616 e_tot = obj.getBoundaryOperator('e_tot', boundary);
542 case {'w','W','west','West','s','S','south','South'} 617 alpha = obj.getBoundaryOperator('alpha', boundary);
543 varargout{k} = obj.d1_l{j}; 618 E = obj.E;
544 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 619 [m, n] = size(alpha{1,1});
545 varargout{k} = obj.d1_r{j}; 620 alpha_tot = sparse(m*obj.dim, n*obj.dim);
621 for i = 1:obj.dim
622 for l = 1:obj.dim
623 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
546 end 624 end
547 625 end
548 case 'T' 626 o = alpha_tot;
549 switch boundary
550 case {'w','W','west','West','s','S','south','South'}
551 varargout{k} = obj.T_l{j};
552 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
553 varargout{k} = obj.T_r{j};
554 end
555
556 case 'tau'
557 switch boundary
558 case {'w','W','west','West','s','S','south','South'}
559 varargout{k} = obj.tau_l{j};
560 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
561 varargout{k} = obj.tau_r{j};
562 end
563
564 case 'tau_tot'
565 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
566
567 I_dim = speye(obj.dim, obj.dim);
568 e_tot = kron(e, I_dim);
569 E = obj.E;
570 tau_tot = (e_tot'*E{1}*e*tau{1}')';
571 for i = 2:obj.dim
572 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
573 end
574 varargout{k} = tau_tot;
575
576 case 'H'
577 varargout{k} = obj.H_boundary{j};
578
579 case 'alpha'
580 % alpha = alpha(i,j) is the penalty strength for displacement BC.
581 e = obj.getBoundaryOperator('e', boundary);
582
583 LAMBDA = obj.LAMBDA;
584 MU = obj.MU;
585
586 dim = obj.dim;
587 theta_R = obj.theta_R{j};
588 theta_H = obj.theta_H{j};
589 theta_M = obj.theta_M{j};
590
591 a_lambda = dim/theta_H + 1/theta_R;
592 a_mu_i = 2/theta_M;
593 a_mu_ij = 2/theta_H + 1/theta_R;
594
595 d = @kroneckerDelta; % Kronecker delta
596 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
597 alpha = cell(obj.dim, obj.dim);
598
599 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
600 + d(i,j)* a_mu_i*MU ...
601 + db(i,j)*a_mu_ij*MU;
602 for i = 1:obj.dim
603 for l = 1:obj.dim
604 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
605 end
606 end
607
608 varargout{k} = alpha;
609
610 case 'alpha_tot'
611 % alpha = alpha(i,j) is the penalty strength for displacement BC.
612 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
613 E = obj.E;
614 [m, n] = size(alpha{1,1});
615 alpha_tot = sparse(m*obj.dim, n*obj.dim);
616 for i = 1:obj.dim
617 for l = 1:obj.dim
618 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
619 end
620 end
621 varargout{k} = alpha_tot;
622
623 otherwise
624 error(['No such operator: operator = ' op{k}]);
625 end
626 end 627 end
627 628
628 end 629 end
629 630
630 function H = getBoundaryQuadrature(obj, boundary) 631 function H = getBoundaryQuadrature(obj, boundary)