Mercurial > repos > public > sbplib
changeset 1057:ff274c7404cc feature/poroelastic
In Elastic2dVariable: copy-paste from feature/getBoundaryOp and fix bugs introduced there.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 25 Jan 2019 15:52:21 -0800 |
| parents | b4fa176b4287 |
| children | 84933722ec0e |
| files | +scheme/Elastic2dVariable.m |
| diffstat | 1 files changed, 100 insertions(+), 99 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dVariable.m Fri Jan 25 15:15:44 2019 -0800 +++ b/+scheme/Elastic2dVariable.m Fri Jan 25 15:52:21 2019 -0800 @@ -362,7 +362,10 @@ % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); - [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); + e = obj.getBoundaryOperator('e', boundary); + T = obj.getBoundaryOperator('T', boundary); + tau = obj.getBoundaryOperator('tau', boundary); + H_gamma = obj.getBoundaryOperator('H', boundary); E = obj.E; Hi = obj.Hi; @@ -433,8 +436,12 @@ % Operators without subscripts are from the own domain. % Get boundary operators - [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary); - [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary); + e = obj.getBoundaryOperator('e_tot', boundary); + tau = obj.getBoundaryOperator('tau_tot', boundary); + + e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary); + tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); % Operators and quantities that correspond to the own domain only @@ -505,124 +512,118 @@ end % Returns the boundary operator op for the boundary specified by the string boundary. - % op: may be a cell array of strings + % op -- string % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() - function [varargout] = getBoundaryOperator(obj, op, boundary) + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'}) switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} + case {'w', 'e'} j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} + case {'s', 'n'} j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); - end - - if ~iscell(op) - op = {op}; end - for k = 1:length(op) - switch op{k} - case 'e' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.e_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.e_r{j}; - end + switch op + case 'e' + switch boundary + case {'w', 's'} + o = obj.e_l{j}; + case {'e', 'n'} + o = obj.e_r{j}; + end - case 'e_tot' - e = obj.getBoundaryOperator('e', boundary); - I_dim = speye(obj.dim, obj.dim); - varargout{k} = kron(e, I_dim); - - case 'd' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.d1_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.d1_r{j}; - end + case 'e_tot' + e = obj.getBoundaryOperator('e', boundary); + I_dim = speye(obj.dim, obj.dim); + o = kron(e, I_dim); - case 'T' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.T_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.T_r{j}; - end + case 'd' + switch boundary + case {'w', 's'} + o = obj.d1_l{j}; + case {'e', 'n'} + o = obj.d1_r{j}; + end - case 'tau' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.tau_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.tau_r{j}; - end + case 'T' + switch boundary + case {'w', 's'} + o = obj.T_l{j}; + case {'e', 'n'} + o = obj.T_r{j}; + end - case 'tau_tot' - [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary); + case 'tau' + switch boundary + case {'w', 's'} + o = obj.tau_l{j}; + case {'e', 'n'} + o = obj.tau_r{j}; + end - I_dim = speye(obj.dim, obj.dim); - e_tot = kron(e, I_dim); - E = obj.E; - tau_tot = (e_tot'*E{1}*e*tau{1}')'; - for i = 2:obj.dim - tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')'; - end - varargout{k} = tau_tot; + case 'tau_tot' + e = obj.getBoundaryOperator('e', boundary); + tau = obj.getBoundaryOperator('tau', boundary); - case 'H' - varargout{k} = obj.H_boundary{j}; + I_dim = speye(obj.dim, obj.dim); + e_tot = kron(e, I_dim); + E = obj.E; + tau_tot = (e_tot'*E{1}*e*tau{1}')'; + for i = 2:obj.dim + tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')'; + end + o = tau_tot; - case 'alpha' - % alpha = alpha(i,j) is the penalty strength for displacement BC. - e = obj.getBoundaryOperator('e', boundary); - - LAMBDA = obj.LAMBDA; - MU = obj.MU; + case 'H' + o = obj.H_boundary{j}; - dim = obj.dim; - theta_R = obj.theta_R{j}; - theta_H = obj.theta_H{j}; - theta_M = obj.theta_M{j}; + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + e = obj.getBoundaryOperator('e', boundary); + + LAMBDA = obj.LAMBDA; + MU = obj.MU; - a_lambda = dim/theta_H + 1/theta_R; - a_mu_i = 2/theta_M; - a_mu_ij = 2/theta_H + 1/theta_R; + dim = obj.dim; + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta - alpha = cell(obj.dim, obj.dim); + a_lambda = dim/theta_H + 1/theta_R; + a_mu_i = 2/theta_M; + a_mu_ij = 2/theta_H + 1/theta_R; - alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ... - + d(i,j)* a_mu_i*MU ... - + db(i,j)*a_mu_ij*MU; - for i = 1:obj.dim - for l = 1:obj.dim - alpha{i,l} = d(i,l)*alpha_func(i,j)*e; - end + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + alpha = cell(obj.dim, obj.dim); + + alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU; + for i = 1:obj.dim + for l = 1:obj.dim + alpha{i,l} = d(i,l)*alpha_func(i,j)*e; end + end - varargout{k} = alpha; + o = alpha; - case 'alpha_tot' - % alpha = alpha(i,j) is the penalty strength for displacement BC. - [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary); - E = obj.E; - [m, n] = size(alpha{1,1}); - alpha_tot = sparse(m*obj.dim, n*obj.dim); - for i = 1:obj.dim - for l = 1:obj.dim - alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; - end + case 'alpha_tot' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + e = obj.getBoundaryOperator('e', boundary); + e_tot = obj.getBoundaryOperator('e_tot', boundary); + alpha = obj.getBoundaryOperator('alpha', boundary); + E = obj.E; + [m, n] = size(alpha{1,1}); + alpha_tot = sparse(m*obj.dim, n*obj.dim); + for i = 1:obj.dim + for l = 1:obj.dim + alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; end - varargout{k} = alpha_tot; - - otherwise - error(['No such operator: operator = ' op{k}]); - end + end + o = alpha_tot; end end