Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropicIncompatible.m @ 1347:ac54767ae1fb feature/poroelastic tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interface, not fully compatible.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 30 Apr 2024 14:58:35 +0200
parents 464e6f65c2c6
children
comparison
equal deleted inserted replaced
1346:464e6f65c2c6 1347:ac54767ae1fb
530 % Fields: 530 % Fields:
531 % -- tuning: penalty strength, defaults to 1.0 531 % -- tuning: penalty strength, defaults to 1.0
532 % -- interpolation: type of interpolation, default 'none' 532 % -- interpolation: type of interpolation, default 'none'
533 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 533 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
534 534
535 defaultType.tuning = 1.0; 535 defaultType.tuning = 1.2;
536 defaultType.interpolation = 'none'; 536 defaultType.interpolation = 'none';
537 default_struct('type', defaultType); 537 default_struct('type', defaultType);
538 538
539 switch type.interpolation 539 switch type.interpolation
540 case {'none', ''} 540 case {'none', ''}
556 v = neighbour_scheme; 556 v = neighbour_scheme;
557 557
558 % Operators, u side 558 % Operators, u side
559 e_u = u.getBoundaryOperatorForScalarField('e', boundary); 559 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
560 tau_u = u.getBoundaryOperator('tau', boundary); 560 tau_u = u.getBoundaryOperator('tau', boundary);
561 h11_u = u.getBorrowing(boundary); 561 [h11_u, th_R_u] = u.getBorrowing(boundary);
562 nu_u = u.getNormal(boundary); 562 nu_u = u.getNormal(boundary);
563 beta_u = u.computeBorrowingFromC(boundary);
563 564
564 E_u = u.E; 565 E_u = u.E;
565 C_u = u.C; 566 C_u = u.C;
566 m_tot_u = u.grid.N(); 567 m_tot_u = u.grid.N();
568 beta_inv_u = spdiag(1./beta_u);
567 569
568 % Operators, v side 570 % Operators, v side
569 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary); 571 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
570 tau_v = v.getBoundaryOperator('tau', neighbour_boundary); 572 tau_v = v.getBoundaryOperator('tau', neighbour_boundary);
571 h11_v = v.getBorrowing(neighbour_boundary); 573 [h11_v, th_R_v] = v.getBorrowing(neighbour_boundary);
572 nu_v = v.getNormal(neighbour_boundary); 574 nu_v = v.getNormal(neighbour_boundary);
575 beta_v = v.computeBorrowingFromC(neighbour_boundary);
573 576
574 E_v = v.E; 577 E_v = v.E;
575 C_v = v.C; 578 C_v = v.C;
576 m_tot_v = v.grid.N(); 579 m_tot_v = v.grid.N();
580 beta_inv_v = spdiag(1./beta_v);
581
582 % Compensate at corners
583 [~, N] = size(e_u);
584 corner_weight = ones(N, 1);
585 corner_weight(1) = 2;
586 corner_weight(end) = 2;
587 corner_weight = spdiag(corner_weight);
577 588
578 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings 589 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
579 flipFlag = false; 590 flipFlag = false;
580 e_v_flip = e_v; 591 e_v_flip = e_v;
581 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ... 592 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
634 for k = 1:dim 645 for k = 1:dim
635 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u; 646 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
636 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v; 647 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
637 end 648 end
638 end 649 end
639 650 Z_u = (corner_weight/(4*h11_u) + beta_inv_u/(4*th_R_u))*Z_u;
651 Z_v = (corner_weight/(4*h11_v) + beta_inv_v/(4*th_R_v))*Z_v;
640 if flipFlag 652 if flipFlag
641 Z_v = rot90(Z_v,2); 653 Z_v = rot90(Z_v,2);
642 end 654 end
643 655 Z = -tuning*(Z_u + Z_v);
644 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
645 X = Y + Z*e_u'; 656 X = Y + Z*e_u';
646 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; 657 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
647 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}'; 658 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
648 659
649 end 660 end