Mercurial > repos > public > sbplib
comparison +multiblock/DiffOp.m @ 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 | adae8063ea2f |
| children | 60c875c18de3 |
comparison
equal
deleted
inserted
replaced
| 1315:6c308da9dcbc | 1316:bd8e607768ce |
|---|---|
| 207 | 207 |
| 208 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | 208 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) |
| 209 error('not implemented') | 209 error('not implemented') |
| 210 end | 210 end |
| 211 | 211 |
| 212 function penalties = interfaceForcing(obj, boundary, neighbour_boundary, type) | |
| 213 default_arg('type', []); | |
| 214 switch class(boundary) | |
| 215 case 'cell' | |
| 216 penalties = obj.singleInterfaceForcing(boundary, neighbour_boundary, type); | |
| 217 case 'multiblock.BoundaryGroup' | |
| 218 [n,m] = size(obj.D); | |
| 219 penaltyZero = sparse(n,0); | |
| 220 | |
| 221 for i = 1:length(boundary) | |
| 222 penaltyParts = obj.interfaceForcing(boundary{i}, neighbour_boundary{i}, type); | |
| 223 | |
| 224 % Initalize if this is the first boundary | |
| 225 if i==1 | |
| 226 penalties = cell(numel(penaltyParts), 1); | |
| 227 for j = 1:numel(penaltyParts) | |
| 228 penalties{j} = penaltyZero; | |
| 229 end | |
| 230 end | |
| 231 | |
| 232 for j = 1:numel(penaltyParts) | |
| 233 penalties{j} = [penalties{j}, penaltyParts{j}]; | |
| 234 end | |
| 235 | |
| 236 end | |
| 237 otherwise | |
| 238 error('Unknown boundary indentifier') | |
| 239 end | |
| 240 | |
| 241 end | |
| 242 | |
| 243 function penalties = singleInterfaceForcing(obj, boundary, neighbour_boundary, type) | |
| 244 I = boundary{1}; | |
| 245 b1 = boundary{2}; | |
| 246 | |
| 247 J = neighbour_boundary{1}; | |
| 248 b2 = neighbour_boundary{2}; | |
| 249 | |
| 250 % Get local penalties | |
| 251 [~, ~, blockPenalties] = obj.diffOps{I}.interface(b1, obj.diffOps{J}, b2, type); | |
| 252 % [~, ~, blockPenalties2] = obj.diffOps{J}.interface(b2, b1, obj.diffOps{I}); | |
| 253 | |
| 254 % Expand to matrices for full domain. | |
| 255 n = numel(blockPenalties); | |
| 256 penalties = cell(n, 1); | |
| 257 for i = 1:n | |
| 258 penalties{i} = multiblock.local2globalPenalty(blockPenalties{i}, obj.blockmatrixDiv, I); | |
| 259 % + multiblock.local2globalPenalty(blockPenalties2{i}, obj.blockmatrixDiv, J); | |
| 260 end | |
| 261 end | |
| 262 | |
| 212 % Size returns the number of degrees of freedom | 263 % Size returns the number of degrees of freedom |
| 213 function N = size(obj) | 264 function N = size(obj) |
| 214 N = 0; | 265 N = 0; |
| 215 for i = 1:length(obj.diffOps) | 266 for i = 1:length(obj.diffOps) |
| 216 N = N + obj.diffOps{i}.size(); | 267 N = N + obj.diffOps{i}.size(); |