Mercurial > repos > public > sbplib
changeset 675:90bf651abc7c feature/poroelastic
Add Dirichlet BC to dilation. Stable with variable coeff and conv study yields p+1/2 which probably is ok.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Thu, 18 Jan 2018 13:36:56 -0800 |
parents | dd84b8862aa8 |
children | 9926efb39330 |
files | +scheme/elasticDilationVariable.m |
diffstat | 1 files changed, 27 insertions(+), 88 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/elasticDilationVariable.m Tue Jan 16 17:41:55 2018 -0800 +++ b/+scheme/elasticDilationVariable.m Thu Jan 18 13:36:56 2018 -0800 @@ -172,8 +172,10 @@ % Divergence operator for BC Div = cell(dim,1); + % Loop over boundaries for i = 1:dim Div{i} = sparse(m_tot,dim*m_tot); + % Loop over components for j = 1:dim Div{i} = Div{i} + d(i,j)*(obj.e_l{i}*obj.d1_l{i}' + obj.e_r{i}*obj.d1_r{i}')*E{j}' ... + db(i,j)*obj.D1{j}*E{j}'; @@ -234,83 +236,35 @@ switch type % Dirichlet boundary condition case {'D','d','dirichlet'} - error('Not implemented'); tuning = 1.2; phi = obj.phi; - closures = cell(obj.dim,1); - penalties = cell(obj.dim,obj.dim); - % Loop over components - for i = 1:obj.dim - H_gamma_i = obj.H_boundary{i}; - sigma_ij = tuning*delta(i,j)*2/(gamma*h(j)) +... - tuning*delta_b(i,j)*(2/(H11*h(j)) + 2/(H11*h(j)*phi)); - - ci = E{i}*inv(RHO)*nj*Hi*... - ( (e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{i}' + ... - delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ... - ) ... - - sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{i}'; - - cj = E{j}*inv(RHO)*nj*Hi*... - ( delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{i}' ... - ); + % Initialize with zeros + m_tot = obj.grid.N(); + closure = sparse(m_tot*obj.dim, m_tot*obj.dim); + penalty = 0*E{j}*e{j}; - if isempty(closures{i}) - closures{i} = ci; - else - closures{i} = closures{i} + ci; - end - - if isempty(closures{j}) - closures{j} = cj; - else - closures{j} = closures{j} + cj; - end - - pii = - E{i}*inv(RHO)*nj*Hi*... - ( (H_gamma*e{j}'*A*e{j}*d{j}')' + ... - delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ... - ) ... - + sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma; + % Loop over components to put penalties on + for i = 1:obj.dim + sigma_i = tuning * obj.dim/(H11*h(j)) +... + tuning * 1/(H11*h(j)*phi); - pji = - E{j}*inv(RHO)*nj*Hi*... - ( delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' ); - - % Dummies - pij = - 0*E{i}*e{j}; - pjj = - 0*E{j}*e{j}; - - if isempty(penalties{i,i}) - penalties{i,i} = pii; - else - penalties{i,i} = penalties{i,i} + pii; - end - - if isempty(penalties{j,i}) - penalties{j,i} = pji; - else - penalties{j,i} = penalties{j,i} + pji; - end + closure = closure + E{i}*inv(RHO)*nj*Hi*... + ( ... + delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ... + + delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{j}' ... + ) ... + - delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{j}'; - if isempty(penalties{i,j}) - penalties{i,j} = pij; - else - penalties{i,j} = penalties{i,j} + pij; - end + penalty = penalty - ... + E{i}*inv(RHO)*nj*Hi*... + ( ... + delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ... + + delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' ... + ) ... + + delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma; - if isempty(penalties{j,j}) - penalties{j,j} = pjj; - else - penalties{j,j} = penalties{j,j} + pjj; - end end - [rows, cols] = size(closures{1}); - closure = sparse(rows, cols); - for i = 1:obj.dim - closure = closure + closures{i}; - end - penalty = penalties; % Free boundary condition case {'F','f','Free','free'} @@ -320,25 +274,10 @@ % Divergence operator Div = obj.Div{j}; - % Loop over components - %for i = 1:obj.dim - closure = -nj*E{j}*inv(RHO)*Hi*e{j} ... - *H_gamma*e{j}'*A*e{j}*e{j}'*Div; - penalty = nj*E{j}*inv(RHO)*Hi*e{j} ... - *H_gamma*e{j}'*A*e{j}; - %end - % [rows, cols] = size(closures{1}); - % closure = sparse(rows, cols); - % for i = 1:obj.dim - % closure = closure + closures{i}; - % for j = 1:obj.dim - % if i~=j - % [rows cols] = size(penalties{j,j}); - % penalties{i,j} = sparse(rows,cols); - % end - % end - % end - % penalty = penalties; + closure = -nj*E{j}*inv(RHO)*Hi*e{j} ... + *H_gamma*e{j}'*A*e{j}*e{j}'*Div; + penalty = nj*E{j}*inv(RHO)*Hi*e{j} ... + *H_gamma*e{j}'*A*e{j}; % Unknown boundary condition