Mercurial > repos > public > sbplib
changeset 958:72cd29107a9a feature/poroelastic
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 05 Dec 2018 18:58:10 -0800 |
| parents | e30aaa4a3e09 |
| children | c226fb8c2b8a |
| files | +multiblock/DiffOp.m +scheme/Elastic2dVariable.m |
| diffstat | 2 files changed, 98 insertions(+), 27 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Wed Nov 28 14:04:31 2018 -0800 +++ b/+multiblock/DiffOp.m Wed Dec 05 18:58:10 2018 -0800 @@ -170,6 +170,50 @@ end end + % Get a boundary cell of operators, specified by opName for the given boundary/BoundaryGroup + function opCell = getBoundaryCellOperator(obj, opName, boundary, blockmatrixDiv) + default_arg('blockmatrixDiv', obj.blockmatrixDiv); + + % Get size of cell + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2}); + case 'multiblock.BoundaryGroup' + blockId = boundary{1}{1}; + localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{1}{2}); + otherwise + error('Unknown boundary indentifier') + end + + % Loop over cell elements and build the boundary operator in each cell + opCell = cell(size(localCell)); + for i = 1:numel(opCell) + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localOpCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2}); + localOp = localOpCell{i}; + + div = {blockmatrixDiv, size(localOp,2)} + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + opCell{i} = op; + + case 'multiblock.BoundaryGroup' + op = sparse(size(obj.D,1),0); + for j = 1:length(boundary) + localCell = obj.getBoundaryCellOperator(opName, boundary{j}, blockmatrixDiv); + op = [op, localCell{i}]; + end + opCell{i} = op; + otherwise + error('Unknown boundary indentifier') + end + end + end + function op = getBoundaryQuadrature(obj, boundary) opName = 'H'; switch class(boundary)
--- a/+scheme/Elastic2dVariable.m Wed Nov 28 14:04:31 2018 -0800 +++ b/+scheme/Elastic2dVariable.m Wed Dec 05 18:58:10 2018 -0800 @@ -238,20 +238,28 @@ tau_l{j} = cell(dim,1); tau_r{j} = cell(dim,1); + LAMBDA_l = e_l{j}'*LAMBDA*e_l{j}; + LAMBDA_r = e_r{j}'*LAMBDA*e_r{j}; + MU_l = e_l{j}'*MU*e_l{j}; + MU_r = e_r{j}'*MU*e_r{j}; + + [~, n_l] = size(e_l{j}); + [~, n_r] = size(e_r{j}); + % Loop over components for i = 1:dim - tau_l{j}{i} = sparse(m_tot,dim*m_tot); - tau_r{j}{i} = sparse(m_tot,dim*m_tot); + tau_l{j}{i} = sparse(n_l, dim*m_tot); + tau_r{j}{i} = sparse(n_r, dim*m_tot); for k = 1:dim 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(i,k)*MU*e_l{j}*d1_l{j}'; + -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{k}' + db(i,k)*e_l{j}'*D1{k})... + -d(j,k)*MU_l*(d(i,j)*d1_l{i}' + db(i,j)*e_l{j}'*D1{i})... + -d(i,k)*MU_l*d1_l{j}'; 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(i,k)*MU*e_r{j}*d1_r{j}'; + d(i,j)*LAMBDA_r*(d(i,k)*d1_r{k}' + db(i,k)*e_r{j}'*D1{k})... + +d(j,k)*MU_r*(d(i,j)*d1_r{i}' + db(i,j)*e_r{j}'*D1{i})... + +d(i,k)*MU_r*d1_r{j}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; @@ -259,6 +267,19 @@ end end + + % Transpose T and tau to match boundary operator convention + for i = 1:dim + for j = 1:dim + tau_l{i}{j} = transpose(tau_l{i}{j}); + tau_r{i}{j} = transpose(tau_r{i}{j}); + for k = 1:dim + T_l{i}{j,k} = transpose(T_l{i}{j,k}); + T_r{i}{j,k} = transpose(T_r{i}{j,k}); + end + end + end + obj.T_l = T_l; obj.T_r = T_r; obj.tau_l = tau_l; @@ -344,16 +365,16 @@ % Loop over components that Dirichlet penalties end up on for i = 1:dim - C = T{k,i}; + C = transpose(T{k,i}); A = -d(i,k)*alpha(i,j); - B = A + C; + B = A + e*C; closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 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*tau{k}'; penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; % Unknown boundary condition @@ -434,8 +455,8 @@ % 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}'; + closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}'; + penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}'; end end end @@ -477,37 +498,37 @@ op = {op}; end - for i = 1:length(op) - switch op{i} + for k = 1:length(op) + switch op{k} case 'e' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.e_l{j}; + varargout{k} = obj.e_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.e_r{j}; + varargout{k} = obj.e_r{j}; end case 'd' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.d1_l{j}; + varargout{k} = obj.d1_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.d1_r{j}; + varargout{k} = obj.d1_r{j}; end case 'H' - varargout{i} = obj.H_boundary{j}; + varargout{k} = obj.H_boundary{j}; case 'T' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.T_l{j}; + varargout{k} = obj.T_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.T_r{j}; + varargout{k} = obj.T_r{j}; end case 'tau' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.tau_l{j}; + varargout{k} = obj.tau_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.tau_r{j}; + varargout{k} = obj.tau_r{j}; end case 'alpha' % alpha = alpha(i,j) is the penalty strength for displacement BC. @@ -526,13 +547,19 @@ 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 ... + alpha = cell(obj.dim, obj.dim); + + for i = 1:obj.dim + for j = 1:obj.dim + alpha{i,j} = tuning*( d(i,j)* a_lambda*LAMBDA ... + d(i,j)* a_mu_i*MU ... + db(i,j)*a_mu_ij*MU ); + end + end - varargout{i} = alpha; + varargout{k} = alpha; otherwise - error(['No such operator: operator = ' op{i}]); + error(['No such operator: operator = ' op{k}]); end end