Mercurial > repos > public > sbplib
changeset 729:aa8cf3851de8 feature/poroelastic
Update multiblock.DiffOp to work for systems.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 20 Apr 2018 16:56:49 -0700 |
| parents | 0aff87f6fb2c |
| children | 9f28cf266f86 |
| files | +multiblock/DiffOp.m +scheme/Elastic2dVariable.m |
| diffstat | 2 files changed, 20 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Fri Apr 20 12:06:38 2018 -0700 +++ b/+multiblock/DiffOp.m Fri Apr 20 16:56:49 2018 -0700 @@ -53,7 +53,11 @@ % Build the differentiation matrix - obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + Ns = zeros(nBlocks,1); + for i = 1:nBlocks + Ns(i) = length(obj.diffOps{i}.D); + end + obj.blockmatrixDiv = {Ns, Ns}; D = blockmatrix.zero(obj.blockmatrixDiv); for i = 1:nBlocks D{i,i} = obj.diffOps{i}.D;
--- a/+scheme/Elastic2dVariable.m Fri Apr 20 12:06:38 2018 -0700 +++ b/+scheme/Elastic2dVariable.m Fri Apr 20 16:56:49 2018 -0700 @@ -64,6 +64,13 @@ h = g.scaling(); lim = g.lim; + if isempty(lim) + x = g.x; + lim = cell(length(x),1); + for i = 1:length(x) + lim{i} = {min(x{i}), max(x{i})}; + end + end % 1D operators ops = cell(dim,1); @@ -270,6 +277,10 @@ default_arg('type',{'free','free'}); default_arg('parameter', []); + if ~iscell(type) + type = {type, type}; + end + % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); @@ -407,11 +418,11 @@ % Loop over components that penalties end up on for i = 1:dim - closure = closure - sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e'*E{i}'; - penalty = penalty + sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e_v'*E_v{i}'; + closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}'; + penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}'; - closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}'; - penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}'; + closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; + penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; % Loop over components that we have interface conditions on for k = 1:dim