Mercurial > repos > public > sbplib
diff +multiblock/DiffOp.m @ 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 | 7d4f57725192 |
| children | 262b52c3f268 |
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)