Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 972:104f0af001e0 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up Elastic2dVariable.interfaceStandard
author Martin Almquist <malmquist@stanford.edu>
date Tue, 25 Dec 2018 08:33:35 +0100
parents 23d9ca6755be
children f5a99cca4654
comparison
equal deleted inserted replaced
970:23d9ca6755be 972:104f0af001e0
436 436
437 % u denotes the solution in the own domain 437 % u denotes the solution in the own domain
438 % v denotes the solution in the neighbour domain 438 % v denotes the solution in the neighbour domain
439 % Operators without subscripts are from the own domain. 439 % Operators without subscripts are from the own domain.
440 440
441 % j is the coordinate direction of the boundary
442 j = obj.get_boundary_number(boundary);
443 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
444
445 % Get boundary operators 441 % Get boundary operators
446 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); 442 [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary);
447 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary); 443 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary);
444 H_gamma = obj.getBoundaryQuadrature(boundary);
448 445
449 % Operators and quantities that correspond to the own domain only 446 % Operators and quantities that correspond to the own domain only
450 Hi = obj.Hi; 447 Hi = obj.Hi_kron;
451 RHOi = obj.RHOi; 448 RHOi = obj.RHOi_kron;
452 dim = obj.dim; 449
453 450 % Penalty strength operators
454 %--- Other operators ---- 451 alpha_u = 1/4*obj.getBoundaryOperator('alpha_tot', boundary);
455 m_tot_u = obj.grid.N(); 452 alpha_v = 1/4*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary);
456 E = obj.E; 453
457 LAMBDA_u = obj.LAMBDA; 454 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
458 MU_u = obj.MU; 455 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
459 lambda_u = e'*LAMBDA_u*e; 456
460 mu_u = e'*MU_u*e; 457 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
461 458 penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
462 m_tot_v = neighbour_scheme.grid.N(); 459
463 E_v = neighbour_scheme.E; 460 closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
464 LAMBDA_v = neighbour_scheme.LAMBDA; 461 penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
465 MU_v = neighbour_scheme.MU; 462
466 lambda_v = e_v'*LAMBDA_v*e_v;
467 mu_v = e_v'*MU_v*e_v;
468 %-------------------------
469
470 % Borrowing constants
471 theta_R_u = obj.theta_R{j};
472 theta_H_u = obj.theta_H{j};
473 theta_M_u = obj.theta_M{j};
474
475 theta_R_v = neighbour_scheme.theta_R{j_v};
476 theta_H_v = neighbour_scheme.theta_H{j_v};
477 theta_M_v = neighbour_scheme.theta_M{j_v};
478
479 function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu)
480 alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M);
481 alpha_ij = mu/(2*th_H) + mu/(4*th_R);
482 end
483
484 [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u);
485 [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v);
486 sigma_ii = alpha_ii_u + alpha_ii_v;
487 sigma_ij = alpha_ij_u + alpha_ij_v;
488
489 d = @kroneckerDelta; % Kronecker delta
490 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
491 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
492
493 % Preallocate
494 closure = sparse(dim*m_tot_u, dim*m_tot_u);
495 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
496
497 % Loop over components that penalties end up on
498 for i = 1:dim
499 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
500 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
501
502 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau{i}';
503 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau_v{i}';
504
505 % Loop over components that we have interface conditions on
506 for k = 1:dim
507 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}';
508 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}';
509 end
510 end
511 end 463 end
512 464
513 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 465 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
514 error('Non-conforming interfaces not implemented yet.'); 466 error('Non-conforming interfaces not implemented yet.');
515 end 467 end