Mercurial > repos > public > sbplib
comparison +scheme/Elastic2dVariableAnisotropic.m @ 1252:8fc2f9a4c882 feature/poroelastic
Add (temporary) fix for e-s , w-n, and x-x couplings
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 07 Jan 2020 12:49:09 -0800 |
| parents | 7f427275bc9c |
| children | 89dad61cad22 |
comparison
equal
deleted
inserted
replaced
| 1227:02dfe3a56742 | 1252:8fc2f9a4c882 |
|---|---|
| 503 | 503 |
| 504 E_v = v.E; | 504 E_v = v.E; |
| 505 C_v = v.C; | 505 C_v = v.C; |
| 506 m_tot_v = v.grid.N(); | 506 m_tot_v = v.grid.N(); |
| 507 | 507 |
| 508 % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings | |
| 509 flipFlag = false; | |
| 510 e_v_flip = e_v; | |
| 511 if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ... | |
| 512 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ... | |
| 513 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ... | |
| 514 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ... | |
| 515 (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ... | |
| 516 (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ... | |
| 517 (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ... | |
| 518 (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e')) | |
| 519 | |
| 520 flipFlag = true; | |
| 521 e_v_flip = fliplr(e_v); | |
| 522 | |
| 523 t1 = tau_v(:,1:2:end-1); | |
| 524 t2 = tau_v(:,2:2:end); | |
| 525 | |
| 526 t1 = fliplr(t1); | |
| 527 t2 = fliplr(t2); | |
| 528 | |
| 529 tau_v(:,1:2:end-1) = t1; | |
| 530 tau_v(:,2:2:end) = t2; | |
| 531 end | |
| 532 | |
| 508 % Operators that are only required for own domain | 533 % Operators that are only required for own domain |
| 509 Hi = u.Hi_kron; | 534 Hi = u.Hi_kron; |
| 510 RHOi = u.RHOi_kron; | 535 RHOi = u.RHOi_kron; |
| 511 e_kron = u.getBoundaryOperator('e', boundary); | 536 e_kron = u.getBoundaryOperator('e', boundary); |
| 512 T_u = u.getBoundaryTractionOperator(boundary); | 537 T_u = u.getBoundaryTractionOperator(boundary); |
| 539 for k = 1:dim | 564 for k = 1:dim |
| 540 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u; | 565 Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u; |
| 541 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v; | 566 Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v; |
| 542 end | 567 end |
| 543 end | 568 end |
| 569 | |
| 570 if flipFlag | |
| 571 Z_v = rot90(Z_v,2); | |
| 572 end | |
| 573 | |
| 544 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); | 574 Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); |
| 545 X = Y + Z*e_u'; | 575 X = Y + Z*e_u'; |
| 546 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; | 576 closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; |
| 547 penalty = penalty - E_u{i}*X'*H_gamma*e_v'*E_v{j}'; | 577 penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}'; |
| 578 | |
| 548 end | 579 end |
| 549 end | 580 end |
| 550 | 581 |
| 551 % ---- Continuity of traction ------ | 582 % ---- Continuity of traction ------ |
| 552 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u'; | 583 closure = closure - 1/2*e_kron*H_gamma_kron*tau_u'; |