Mercurial > repos > public > sbplib
changeset 855:5751262b323b feature/poroelastic
Add 1D quadrature matrices as property in Elastic2dVariable.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 16 Sep 2018 18:00:14 -0700 |
| parents | bd43256f0c41 |
| children | ee4cfb37534d |
| files | +scheme/Elastic2dVariable.m |
| diffstat | 1 files changed, 21 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dVariable.m Fri Sep 14 10:41:00 2018 -0700 +++ b/+scheme/Elastic2dVariable.m Sun Sep 16 18:00:14 2018 -0700 @@ -1,7 +1,7 @@ classdef Elastic2dVariable < scheme.Scheme % Discretizes the elastic wave equation: -% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i +% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -30,14 +30,14 @@ T_l, T_r tau_l, tau_r - H, Hi % Inner products + H, Hi, H_1D % Inner products phi % Borrowing constant for (d1 - e^T*D1) from R 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 % Boundary inner products % Kroneckered norms and coefficients @@ -171,6 +171,7 @@ obj.H_boundary = cell(dim,1); obj.H_boundary{1} = H{2}; obj.H_boundary{2} = H{1}; + obj.H_1D = {H{1}, H{2}}; % E{i}^T picks out component i. E = cell(dim,1); @@ -230,14 +231,14 @@ tau_l{j}{i} = sparse(m_tot,dim*m_tot); tau_r{j}{i} = sparse(m_tot,dim*m_tot); for k = 1:dim - T_l{j}{i,k} = ... + T_l{j}{i,k} = ... -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... - -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... + -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... -d(i,k)*MU*e_l{j}*d1_l{j}'; - T_r{j}{i,k} = ... + T_r{j}{i,k} = ... d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... - +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... + +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... +d(i,k)*MU*e_r{j}*d1_r{j}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; @@ -270,7 +271,7 @@ % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition - % on the first component. + % on the first component. % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. @@ -317,20 +318,20 @@ 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 ); + + 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*B'*e*H_gamma*(e'*E{k}' ); + closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; - end + end % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); + closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; % Unknown boundary condition @@ -357,7 +358,7 @@ Hi = obj.Hi; RHOi = obj.RHOi; dim = obj.dim; - + %--- Other operators ---- m_tot_u = obj.grid.N(); E = obj.E; @@ -373,7 +374,7 @@ 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); @@ -386,7 +387,7 @@ 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) + function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) th1 = h11/(2*dim); th2 = h11*phi/2; th3 = h*gamma; @@ -398,7 +399,7 @@ 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); + [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; @@ -420,9 +421,9 @@ % 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 + 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 @@ -494,7 +495,7 @@ varargout{i} = obj.tau_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} varargout{i} = obj.tau_r{j}; - end + end otherwise error(['No such operator: operator = ' op{i}]); end