Mercurial > repos > public > sbplib
changeset 1316:bd8e607768ce feature/poroelastic
Add interfaceForcing method in multiblock.DiffOp to facilitate mms forcing
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 26 Jul 2020 20:05:10 -0700 |
| parents | 6c308da9dcbc |
| children | 34997aced843 |
| files | +multiblock/DiffOp.m |
| diffstat | 1 files changed, 51 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Sun Jul 26 17:38:28 2020 -0700 +++ b/+multiblock/DiffOp.m Sun Jul 26 20:05:10 2020 -0700 @@ -209,6 +209,57 @@ error('not implemented') end + function penalties = interfaceForcing(obj, boundary, neighbour_boundary, type) + default_arg('type', []); + switch class(boundary) + case 'cell' + penalties = obj.singleInterfaceForcing(boundary, neighbour_boundary, type); + case 'multiblock.BoundaryGroup' + [n,m] = size(obj.D); + penaltyZero = sparse(n,0); + + for i = 1:length(boundary) + penaltyParts = obj.interfaceForcing(boundary{i}, neighbour_boundary{i}, type); + + % Initalize if this is the first boundary + if i==1 + penalties = cell(numel(penaltyParts), 1); + for j = 1:numel(penaltyParts) + penalties{j} = penaltyZero; + end + end + + for j = 1:numel(penaltyParts) + penalties{j} = [penalties{j}, penaltyParts{j}]; + end + + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function penalties = singleInterfaceForcing(obj, boundary, neighbour_boundary, type) + I = boundary{1}; + b1 = boundary{2}; + + J = neighbour_boundary{1}; + b2 = neighbour_boundary{2}; + + % Get local penalties + [~, ~, blockPenalties] = obj.diffOps{I}.interface(b1, obj.diffOps{J}, b2, type); + % [~, ~, blockPenalties2] = obj.diffOps{J}.interface(b2, b1, obj.diffOps{I}); + + % Expand to matrices for full domain. + n = numel(blockPenalties); + penalties = cell(n, 1); + for i = 1:n + penalties{i} = multiblock.local2globalPenalty(blockPenalties{i}, obj.blockmatrixDiv, I); + % + multiblock.local2globalPenalty(blockPenalties2{i}, obj.blockmatrixDiv, J); + end + end + % Size returns the number of degrees of freedom function N = size(obj) N = 0;