Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1266:ad31d9c4cec2 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add displacement BC in StaggeredAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Thu, 30 Apr 2020 20:55:26 -0700
parents 6696b216b982
children e61083f178be
comparison
equal deleted inserted replaced
1265:6696b216b982 1266:ad31d9c4cec2
147 ed_l{i} = ops{i}.e_dual_l; 147 ed_l{i} = ops{i}.e_dual_l;
148 ed_r{i} = ops{i}.e_dual_r; 148 ed_r{i} = ops{i}.e_dual_r;
149 end 149 end
150 150
151 % Borrowing constants 151 % Borrowing constants
152 for i = 1:dim 152 % for i = 1:dim
153 obj.h11{i} = ops{i}.H_dual(1,1); 153 % obj.h11{i} = ops{i}.H_dual(1,1);
154 end 154 % end
155 obj.h11{1}{1} = ops{1}.H_dual(1,1);
156 obj.h11{1}{2} = ops{2}.H_primal(1,1);
157 obj.h11{2}{1} = ops{1}.H_primal(1,1);
158 obj.h11{2}{2} = ops{2}.H_dual(1,1);
155 159
156 %---- Grid layout ------- 160 %---- Grid layout -------
157 % gu1 = xp o yp; 161 % gu1 = xp o yp;
158 % gu2 = xd o yd; 162 % gu2 = xd o yd;
159 % gs1 = xd o yp; 163 % gs1 = xd o yp;
359 T_s{b,c} = cell(dim, dim); 363 T_s{b,c} = cell(dim, dim);
360 T_n{b,c} = cell(dim, dim); 364 T_n{b,c} = cell(dim, dim);
361 365
362 for i = 1:dim 366 for i = 1:dim
363 for j = 1:dim 367 for j = 1:dim
364 T_w{b,c}{i,j} = sparse(m_we, N_u{c}); 368 T_w{b,c}{i,j} = sparse(N_u{c}, m_we);
365 T_e{b,c}{i,j} = sparse(m_we, N_u{c}); 369 T_e{b,c}{i,j} = sparse(N_u{c}, m_we);
366 T_s{b,c}{i,j} = sparse(m_sn, N_u{c}); 370 T_s{b,c}{i,j} = sparse(N_u{c}, m_sn);
367 T_n{b,c}{i,j} = sparse(m_sn, N_u{c}); 371 T_n{b,c}{i,j} = sparse(N_u{c}, m_sn);
368 end 372 end
369 end 373 end
370 end 374 end
371 end 375 end
372 376
395 tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)'; 399 tau_s{b}{j} = tau_s{b}{j} + (e_s_s{b}'*n_s(i)*sigma_b_ij)';
396 tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)'; 400 tau_n{b}{j} = tau_n{b}{j} + (e_n_s{b}'*n_n(i)*sigma_b_ij)';
397 401
398 S_bc_ijl = C{b}{i,j,k,l}*D1_u2s{b,c}{k}; 402 S_bc_ijl = C{b}{i,j,k,l}*D1_u2s{b,c}{k};
399 403
400 % T_w{b,c}{i,l} = T_w{b,c}{i,l} + (e_w_s{b}'*n_w(i)*S_bc_ijl)'; 404 T_w{b,c}{j,l} = T_w{b,c}{j,l} + (e_w_s{b}'*n_w(i)*S_bc_ijl)';
401 % T_e{b,c}{i,l} = T_e{b,c}{i,l} + (e_e_s{b}'*n_e(i)*S_bc_ijl)'; 405 T_e{b,c}{j,l} = T_e{b,c}{j,l} + (e_e_s{b}'*n_e(i)*S_bc_ijl)';
402 % T_s{b,c}{i,l} = T_s{b,c}{i,l} + (e_s_s{b}'*n_s(i)*S_bc_ijl)'; 406 T_s{b,c}{j,l} = T_s{b,c}{j,l} + (e_s_s{b}'*n_s(i)*S_bc_ijl)';
403 % T_n{b,c}{i,l} = T_n{b,c}{i,l} + (e_n_s{b}'*n_n(i)*S_bc_ijl)'; 407 T_n{b,c}{j,l} = T_n{b,c}{j,l} + (e_n_s{b}'*n_n(i)*S_bc_ijl)';
404 end 408 end
405 end 409 end
406 end 410 end
407 end 411 end
408 end 412 end
504 508
505 e_u = obj.getBoundaryOperatorForScalarField('e_u', boundary); 509 e_u = obj.getBoundaryOperatorForScalarField('e_u', boundary);
506 e_s = obj.getBoundaryOperatorForScalarField('e_s', boundary); 510 e_s = obj.getBoundaryOperatorForScalarField('e_s', boundary);
507 tau = obj.getBoundaryOperator('tau', boundary); 511 tau = obj.getBoundaryOperator('tau', boundary);
508 T = obj.getBoundaryTractionOperator(boundary); 512 T = obj.getBoundaryTractionOperator(boundary);
509 % h11 = obj.getBorrowing(boundary);
510 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); 513 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
511 nu = obj.getNormal(boundary); 514 nu = obj.getNormal(boundary);
512 515
513 U = obj.U; 516 U = obj.U;
514 G = obj.G; 517 G = obj.G;
537 540
538 % Preallocate 541 % Preallocate
539 [~, col] = size(tau{1}{1}); 542 [~, col] = size(tau{1}{1});
540 closure = sparse(m_tot, m_tot); 543 closure = sparse(m_tot, m_tot);
541 penalty = cell(nGrids, 1); 544 penalty = cell(nGrids, 1);
545 for a = 1:nGrids
546 [~, col] = size(e_u{a});
547 penalty{a} = sparse(m_tot, col);
548 end
542 549
543 j = comp; 550 j = comp;
544 switch type 551 switch type
545 552
546 % Dirichlet boundary condition 553 % Dirichlet boundary condition
557 % Z: symmetric part of penalty 564 % Z: symmetric part of penalty
558 % X = Y + Z. 565 % X = Y + Z.
559 566
560 % Nonsymmetric part goes on all components to 567 % Nonsymmetric part goes on all components to
561 % yield traction in discrete energy rate 568 % yield traction in discrete energy rate
562 for a = 1:nGrids 569 for c = 1:nGrids
563 for b = 1:nGrids 570 for m = 1:numel(gridCombos)
571 gc = gridCombos{m};
572 a = gc{1};
573 b = gc{2};
574
564 for i = 1:dim 575 for i = 1:dim
565 Y = T{a,b}{j,i}'; 576 Y = T{a,c}{j,i}';
566 X = e_s{a}*Y; 577 closure = closure + G{c}*U{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(e_u{b}'*U{b}{j}'*G{b}') ));
567 closure = closure + U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a}*(e_u{b}'*U{j}'*G{b}' ); 578 penalty{b} = penalty{b} - G{c}*U{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}) );
568 penalty = penalty - U{i}*RHOi*Hi*X'*e_s{a}*H_gamma{a};
569 end 579 end
570 end 580 end
571 end 581 end
572 582
573 % Symmetric part only required on components with displacement BC. 583 % Symmetric part only required on components with displacement BC.
574 % (Otherwise it's not symmetric.) 584 % (Otherwise it's not symmetric.)
575 for i = dComps 585 for m = 1:numel(gridCombos)
576 Z = sparse(m_tot, m_tot); 586 gc = gridCombos{m};
577 for l = 1:dim 587 a = gc{1};
578 for k = 1:dim 588 b = gc{2};
579 Z = Z + nu(l)*C{l,i,k,j}*nu(k); 589
590 h11 = obj.getBorrowing(b, boundary);
591
592 for i = dComps
593 Z = 0*C{b}{1,1,1,1};
594 for l = 1:dim
595 for k = 1:dim
596 Z = Z + nu(l)*C{b}{l,i,k,j}*nu(k);
597 end
580 end 598 end
599 Z = -tuning*dim/h11*Z;
600 X = e_s{b}'*Z*e_s{b};
601 closure = closure + G{a}*U{a}{i}*((RHO{a}*H{a})\(e_u{a}*X'*H_gamma{b}*(e_u{a}'*U{a}{j}'*G{a}' ) ));
602 penalty{a} = penalty{a} - G{a}*U{a}{i}*((RHO{a}*H{a})\(e_u{a}*X'*H_gamma{b} ));
581 end 603 end
582 Z = -tuning*dim/h11*Z;
583 X = Z;
584 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
585 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
586 end 604 end
587 605
588 % Free boundary condition 606 % Free boundary condition
589 case {'F','f','Free','free','traction','Traction','t','T'} 607 case {'F','f','Free','free','traction','Traction','t','T'}
590 for m = 1:numel(gridCombos) 608 for m = 1:numel(gridCombos)
762 end 780 end
763 end 781 end
764 782
765 % Returns h11 for the boundary specified by the string boundary. 783 % Returns h11 for the boundary specified by the string boundary.
766 % op -- string 784 % op -- string
767 function h11 = getBorrowing(obj, boundary) 785 function h11 = getBorrowing(obj, stressGrid, boundary)
768 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 786 assertIsMember(boundary, {'w', 'e', 's', 'n'})
769 787
770 switch boundary 788 switch boundary
771 case {'w','e'} 789 case {'w','e'}
772 h11 = obj.h11{1}; 790 h11 = obj.h11{stressGrid}{1};
773 case {'s', 'n'} 791 case {'s', 'n'}
774 h11 = obj.h11{2}; 792 h11 = obj.h11{stressGrid}{2};
775 end 793 end
776 end 794 end
777 795
778 % Returns the outward unit normal vector for the boundary specified by the string boundary. 796 % Returns the outward unit normal vector for the boundary specified by the string boundary.
779 function nu = getNormal(obj, boundary) 797 function nu = getNormal(obj, boundary)