sbplib: +scheme/Elastic2dCurvilinear.m comparison

comparison +scheme/Elastic2dCurvilinear.m @ 1122:b8bd54332763 feature/poroelastic

Bugfix Elastic2dCurve

author	Martin Almquist <malmquist@stanford.edu>
date	Fri, 10 May 2019 15:26:47 -0700
parents	1f6b2fb69225
children	60c875c18de3

comparison

equal deleted inserted replaced

-:07d0caf915e4
+:b8bd54332763
 classdef Elastic2dCurvilinear < scheme.Scheme
 % Discretizes the elastic wave equation in curvilinear coordinates.
 %
 % Untransformed equation:
 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
 %
 % Transformed equation:
 % J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j
 %                + dk J b_jk mu b_il dl u_j
 %                + dk J b_jk mu b_jl dl u_i
 % opSet should be cell array of opSets, one per dimension. This
 % is useful if we have periodic BC in one direction.
 properties
 m % Number of points in each direction, possibly a vector
 gamma % Borrowing constant for d1 from M
 H11 % First element of H
 e_l, e_r
 d1_l, d1_r % Normal derivatives at the boundary
 E % E{i}^T picks out component i
 H_boundary_l, H_boundary_r % Boundary inner products
 % Kroneckered norms and coefficients
 RHOi_kron
 Ji_kron, J_kron
 xmax = max(ops{1}.x);
 ymax = max(ops{2}.x);
 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
 x_xi = D1Metric{1}*x;
 x_eta = D1Metric{2}*x;
 y_xi = D1Metric{1}*y;
 y_eta = D1Metric{2}*y;
 J = x_xi.*y_eta - x_eta.*y_xi;
 b = cell(dim,dim);
 b{1,1} = y_eta./J;
-b{1,2} = -x_eta./J;
+b{2,1} = -x_eta./J;
-b{2,1} = -y_xi./J;
+b{1,2} = -y_xi./J;
 b{2,2} = x_xi./J;
 % Scale factors for boundary integrals
 beta = cell(dim,1);
 beta{1} = sqrt(x_eta.^2 + y_eta.^2);
 T_l{j}{i,k} = sparse(m_tot,m_tot);
 T_r{j}{i,k} = sparse(m_tot,m_tot);
 for m = 1:dim
 for l = 1:dim
 T_l{j}{i,k} = T_l{j}{i,k} + ...
 -d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
 -d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ...
 -d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m});
 T_r{j}{i,k} = T_r{j}{i,k} + ...
 d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
 d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ...
 d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m});
 end
 end
 T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k};
 T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k};
 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
 end
 d = @kroneckerDelta;  % Kronecker delta
 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
 + d(i,j)* a_mu_i*MU ...
 + db(i,j)*a_mu_ij*MU );
 % Loop over components that Dirichlet penalties end up on
 for i = 1:dim
 C = T{k,i};
 A = -d(i,k)*alpha(i,j);
 B = A + C;
 closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' );
 penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma;
 end
 % Free boundary condition
 case {'F','f','Free','free','traction','Traction','t','T'}
 closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} );
 penalty = penalty + E{k}*RHOi*Ji*Hi*e*H_gamma;
 % Unknown boundary condition
 otherwise
 error('No such boundary condition: type = %s',type);
 % Operators and quantities that correspond to the own domain only
 Hi = obj.Hi;
 RHOi = obj.RHOi;
 dim = obj.dim;
 %--- Other operators ----
 m_tot_u = obj.grid.N();
 E = obj.E;
 LAMBDA_u = obj.LAMBDA;
 MU_u = obj.MU;
 LAMBDA_v = neighbour_scheme.LAMBDA;
 MU_v = neighbour_scheme.MU;
 lambda_v = e_v'*LAMBDA_v*e_v;
 mu_v = e_v'*MU_v*e_v;
 %-------------------------
 % Borrowing constants
 phi_u = obj.phi{j};
 h_u = obj.h(j);
 h11_u = obj.H11{j}*h_u;
 gamma_u = obj.gamma{j};
 h_v = neighbour_scheme.h(j_v);
 h11_v = neighbour_scheme.H11{j_v}*h_v;
 gamma_v = neighbour_scheme.gamma{j_v};
 % E > sum_i 1/(2*alpha_ij)*(tau_i)^2
 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)
 th1 = h11/(2*dim);
 th2 = h11*phi/2;
 th3 = h*gamma;
 a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3);
 a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3);
 alpha_ij = mu*(2/h11 + 1/(phi*h11));
 end
 [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
 [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);
 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
 d = @kroneckerDelta;  % Kronecker delta
 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
 % Loop over components that we have interface conditions on
 for k = 1:dim
 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
 end
 end
 end
 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
 function [j, nj] = get_boundary_number(obj, boundary)
 varargout{i} = obj.d1_l{j};
 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
 varargout{i} = obj.d1_r{j};
 end
 case 'H'
 switch boundary
 case {'w','W','west','West','s','S','south','South'}
 varargout{i} = obj.H_boundary_l{j};
 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
 varargout{i} = obj.H_boundary_r{j};
 end
 switch boundary
 case {'w','W','west','West','s','S','south','South'}
 varargout{i} = obj.tau_l{j};
 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
 varargout{i} = obj.tau_r{j};
 end
 otherwise
 error(['No such operator: operator = ' op{i}]);
 end
 end

Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinear.m @ 1122:b8bd54332763 feature/poroelastic