Mercurial > repos > public > sbplib
changeset 679:247b58a4dbe8 feature/poroelastic
Add support for Dirichlet and Traction BC on different components at the same boundary. Remove some unused variables and improve comments.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 05 Feb 2018 14:45:26 -0800 |
| parents | 06676c40e77f |
| children | cd1a76c38565 |
| files | +scheme/elasticVariable.m |
| diffstat | 1 files changed, 39 insertions(+), 44 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/elasticVariable.m Mon Feb 05 11:06:15 2018 -0800 +++ b/+scheme/elasticVariable.m Mon Feb 05 14:45:26 2018 -0800 @@ -3,7 +3,6 @@ % Discretizes the elastic wave equation: % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i - properties m % Number of points in each direction, possibly a vector h % Grid spacing @@ -16,7 +15,7 @@ % Diagonal matrices for varible coefficients LAMBDA % Variable coefficient, related to dilation MU % Shear modulus, variable coefficient - RHO, RHOi % Density + RHO, RHOi % Density, variable D % Total operator D1 % First derivatives @@ -31,6 +30,7 @@ H, Hi % 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 @@ -38,20 +38,12 @@ H_boundary % Boundary inner products - LAMBDA_boundary_l % Variable coefficients at boundaries - LAMBDA_boundary_r - MU_boundary_l - MU_boundary_r - % Kroneckered norms and coefficients RHOi_kron Hi_kron end methods - % Implements the shear part of the elastic wave equation, i.e. - % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i - % where a = mu. function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) default_arg('opSet',@sbp.D2Variable); @@ -81,6 +73,7 @@ beta = ops{1}.borrowing.R.delta_D; obj.H11 = ops{1}.borrowing.H11; obj.phi = beta/obj.H11; + obj.gamma = ops{1}.borrowing.M.d1; I = cell(dim,1); D1 = cell(dim,1); @@ -168,18 +161,6 @@ obj.H_boundary{1} = H{2}; obj.H_boundary{2} = H{1}; - % Boundary coefficient matrices - obj.LAMBDA_boundary_l = cell(dim,1); - obj.LAMBDA_boundary_r = cell(dim,1); - obj.MU_boundary_l = cell(dim,1); - obj.MU_boundary_r = cell(dim,1); - for i = 1:dim - obj.LAMBDA_boundary_l{i} = obj.e_l{i}'*LAMBDA*obj.e_l{i}; - obj.LAMBDA_boundary_r{i} = obj.e_r{i}'*LAMBDA*obj.e_r{i}; - obj.MU_boundary_l{i} = obj.e_l{i}'*MU*obj.e_l{i}; - obj.MU_boundary_r{i} = obj.e_r{i}'*MU*obj.e_r{i}; - end - % E{i}^T picks out component i. E = cell(dim,1); I = speye(m_tot,m_tot); @@ -277,12 +258,12 @@ % Closure functions return the operators applied to the own domain to close the boundary % 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'. - % type is a string specifying the type of boundary condition if there are several. + % type is a cell array of strings specifying the type of boundary condition for each 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. function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) - default_arg('type','free'); + default_arg('type',{'free','free'}); default_arg('parameter', []); % j is the coordinate direction of the boundary @@ -310,46 +291,60 @@ MU = obj.MU; RHOi = obj.RHOi; - phi = obj.phi; - H11 = obj.H11; - h = obj.h; dim = obj.dim; m_tot = obj.grid.N(); RHOi_kron = obj.RHOi_kron; Hi_kron = obj.Hi_kron; + % Preallocate closure = sparse(dim*m_tot, dim*m_tot); penalty = cell(dim,1); - switch type + for k = 1:dim + penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j)); + end + + % Loop over components that we (potentially) have different BC on + for k = 1:dim + switch type{k} + % Dirichlet boundary condition - case {'D','d','dirichlet'} - error('not implemented') + case {'D','d','dirichlet','Dirichlet'} + tuning = 1.2; phi = obj.phi; + h = obj.h(j); + h11 = obj.H11*h; + gamma = obj.gamma; - sigma = tuning * obj.dim/(H11*h(j)) +... - tuning * 1/(H11*h(j)*phi); + a_lambda = dim/h11 + 1/(h11*phi); + a_mu_i = 2/(gamma*h); + a_mu_ij = 2/h11 + 1/(h11*phi); - closure = - sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma*e{j}'*E{j}' ... - + nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma*e{j}'*E{j}'; + 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 ); - penalty = + sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma ... - - nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma; + % 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{j}*H_gamma*(e{j}'*E{k}' ); + penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; + end % Free boundary condition - case {'F','f','Free','free'} - - % Loop over components of traction - for i = 1:dim - closure = closure - E{i}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{i} ); - penalty{i} = E{i}*RHOi*Hi*e{j}*H_gamma; - end - + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); + penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; % Unknown boundary condition otherwise error('No such boundary condition: type = %s',type); + end end end