Mercurial > repos > public > sbplib
changeset 676:9926efb39330 feature/poroelastic
Clean up dilation BC code.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 18 Jan 2018 14:36:59 -0800 |
| parents | 90bf651abc7c |
| children | eeaf9a00e304 |
| files | +scheme/elasticDilationVariable.m |
| diffstat | 1 files changed, 26 insertions(+), 39 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/elasticDilationVariable.m Thu Jan 18 13:36:56 2018 -0800 +++ b/+scheme/elasticDilationVariable.m Thu Jan 18 14:36:59 2018 -0800 @@ -9,7 +9,7 @@ order % Order accuracy for the approximation A % Variable coefficient lambda of the operator (as diagonal matrix here) - RHO % Density (as diagonal matrix here) + RHO, RHOi % Density (as diagonal matrix here) D % Total operator D1 % First derivatives @@ -26,6 +26,10 @@ A_boundary_l % Variable coefficient at boundaries A_boundary_r % + + % Kroneckered norms and coefficients + RHOi_kron + Hi_kron end methods @@ -87,7 +91,7 @@ obj.A = A; RHO = spdiag(rho); obj.RHO = RHO; - + obj.RHOi = inv(RHO); obj.D1 = cell(dim,1); obj.D2 = cell(dim,1); @@ -183,6 +187,11 @@ end obj.Div = Div; + % Kroneckered norms and coefficients + I_dim = speye(dim); + obj.RHOi_kron = kron(obj.RHOi, I_dim); + obj.Hi_kron = kron(obj.Hi, I_dim); + % Misc. obj.m = m; obj.h = h; @@ -195,7 +204,6 @@ % 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. - % Here penalty{i,j} enforces data component j on solution component i % 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. % data is a function returning the data that should be applied at the boundary. @@ -205,9 +213,6 @@ default_arg('type','free'); default_arg('parameter', []); - delta = @kroneckerDelta; % Kronecker delta - delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta - % j is the coordinate direction of the boundary % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries @@ -226,57 +231,39 @@ Hi = obj.Hi; H_gamma = obj.H_boundary{j}; A = obj.A; - RHO = obj.RHO; - D1 = obj.D1; + RHOi = obj.RHOi; phi = obj.phi; H11 = obj.H11; h = obj.h; + RHOi_kron = obj.RHOi_kron; + Hi_kron = obj.Hi_kron; + + % Divergence operator + Div = obj.Div{j}; + switch type % Dirichlet boundary condition case {'D','d','dirichlet'} tuning = 1.2; phi = obj.phi; - % Initialize with zeros - m_tot = obj.grid.N(); - closure = sparse(m_tot*obj.dim, m_tot*obj.dim); - penalty = 0*E{j}*e{j}; - - % 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); + sigma = tuning * obj.dim/(H11*h(j)) +... + tuning * 1/(H11*h(j)*phi); - 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}'; + 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}'; - 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; - - end + penalty = + sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma ... + - nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma; % Free boundary condition case {'F','f','Free','free'} - closures = cell(obj.dim,1); - penalties = cell(obj.dim,obj.dim); - % Divergence operator - Div = obj.Div{j}; - - closure = -nj*E{j}*inv(RHO)*Hi*e{j} ... + closure = -nj*E{j}*RHOi*Hi*e{j} ... *H_gamma*e{j}'*A*e{j}*e{j}'*Div; - penalty = nj*E{j}*inv(RHO)*Hi*e{j} ... + penalty = nj*E{j}*RHOi*Hi*e{j} ... *H_gamma*e{j}'*A*e{j};