Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1060:e40899094f20 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Elastic2dVariable: Clean up alpha in getBoundaryOperator. Add alpha1 and alpha2 as boundary operators.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 26 Jan 2019 16:56:35 -0800
parents d8ab528f10f6
children 07d0caf915e4
comparison
equal deleted inserted replaced
1059:d8ab528f10f6 1060:e40899094f20
547 547
548 % Returns the boundary operator op for the boundary specified by the string boundary. 548 % Returns the boundary operator op for the boundary specified by the string boundary.
549 % op -- string 549 % op -- string
550 function o = getBoundaryOperator(obj, op, boundary) 550 function o = getBoundaryOperator(obj, op, boundary)
551 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 551 assertIsMember(boundary, {'w', 'e', 's', 'n'})
552 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha'}) 552 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha', 'alpha1', 'alpha2'})
553 553
554 switch op 554 switch op
555 555
556 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'} 556 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
557 o = obj.([op, '_', boundary]); 557 o = obj.([op, '_', boundary]);
558 558
559 % Yields vector-valued penalty strength given displacement BC on all components
559 case 'alpha' 560 case 'alpha'
560 % alpha = alpha(i,j) is the penalty strength for displacement BC. 561 e = obj.getBoundaryOperator('e', boundary);
561 e = obj.getBoundaryOperatorForScalarField('e', boundary); 562 e_scalar = obj.getBoundaryOperatorForScalarField('e', boundary);
562 e_tot = obj.getBoundaryOperator('e', boundary); 563 alpha_scalar = obj.getBoundaryOperatorForScalarField('alpha', boundary);
563 alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
564 E = obj.E; 564 E = obj.E;
565 [m, n] = size(alpha{1,1}); 565 [m, n] = size(alpha_scalar{1,1});
566 alpha_tot = sparse(m*obj.dim, n*obj.dim); 566 alpha = sparse(m*obj.dim, n*obj.dim);
567 for i = 1:obj.dim 567 for i = 1:obj.dim
568 for l = 1:obj.dim 568 for l = 1:obj.dim
569 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; 569 alpha = alpha + (e'*E{i}*e_scalar*alpha_scalar{i,l}'*E{l}')';
570 end 570 end
571 end 571 end
572 o = alpha_tot; 572 o = alpha;
573
574 % Yields penalty strength for component 1 given displacement BC on all components
575 case 'alpha1'
576 alpha = obj.getBoundaryOperator('alpha', boundary);
577 e = obj.getBoundaryOperator('e', boundary);
578 e1 = obj.getBoundaryOperator('e1', boundary);
579
580 alpha1 = (e1'*e*alpha')';
581 o = alpha1;
582
583 % Yields penalty strength for component 2 given displacement BC on all components
584 case 'alpha2'
585 alpha = obj.getBoundaryOperator('alpha', boundary);
586 e = obj.getBoundaryOperator('e', boundary);
587 e2 = obj.getBoundaryOperator('e2', boundary);
588
589 alpha2 = (e2'*e*alpha')';
590 o = alpha2;
573 end 591 end
574 592
575 end 593 end
576 594
577 % Returns the boundary operator op for the boundary specified by the string boundary. 595 % Returns the boundary operator op for the boundary specified by the string boundary.
584 602
585 case 'e' 603 case 'e'
586 o = obj.(['e_scalar', '_', boundary]); 604 o = obj.(['e_scalar', '_', boundary]);
587 605
588 case 'alpha' 606 case 'alpha'
589 % alpha = alpha(i,j) is the penalty strength for displacement BC. 607
608 % alpha{i,j} is the penalty strength on component i due to
609 % displacement BC for component j.
590 e = obj.getBoundaryOperatorForScalarField('e', boundary); 610 e = obj.getBoundaryOperatorForScalarField('e', boundary);
591 611
592 LAMBDA = obj.LAMBDA; 612 LAMBDA = obj.LAMBDA;
593 MU = obj.MU; 613 MU = obj.MU;
594 dim = obj.dim; 614 dim = obj.dim;
608 a_mu_i = 2/theta_M; 628 a_mu_i = 2/theta_M;
609 a_mu_ij = 2/theta_H + 1/theta_R; 629 a_mu_ij = 2/theta_H + 1/theta_R;
610 630
611 d = @kroneckerDelta; % Kronecker delta 631 d = @kroneckerDelta; % Kronecker delta
612 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 632 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
613 alpha = cell(obj.dim, obj.dim);
614 633
615 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ... 634 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
616 + d(i,j)* a_mu_i*MU ... 635 + d(i,j)* a_mu_i*MU ...
617 + db(i,j)*a_mu_ij*MU; 636 + db(i,j)*a_mu_ij*MU;
637
638 alpha = cell(obj.dim, obj.dim);
618 for i = 1:obj.dim 639 for i = 1:obj.dim
619 for j = 1:obj.dim 640 for j = 1:obj.dim
620 alpha{i,j} = d(i,j)*alpha_func(i,k)*e; 641 alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
621 end 642 end
622 end 643 end
623
624 o = alpha; 644 o = alpha;
625 end 645 end
626 646
627 end 647 end
628 648
629 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary. 649 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
650 % Formula: tau_i = T_ij u_j
630 % op -- string 651 % op -- string
631 function T = getBoundaryTractionOperator(obj, boundary) 652 function T = getBoundaryTractionOperator(obj, boundary)
632 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 653 assertIsMember(boundary, {'w', 'e', 's', 'n'})
633 654
634 T = obj.(['T', '_', boundary]); 655 T = obj.(['T', '_', boundary]);