Mercurial > repos > public > sbplib
changeset 705:e6fbdc9ccfc4 feature/optim
Add diffOp time Dep
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Mon, 06 Nov 2017 11:39:59 +0100 |
| parents | 111fcbcff2e9 |
| children | a95c89436916 |
| files | +multiblock/DiffOpTimeDep.m |
| diffstat | 1 files changed, 210 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
diff -r 111fcbcff2e9 -r e6fbdc9ccfc4 +multiblock/DiffOpTimeDep.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTimeDep.m Mon Nov 06 11:39:59 2017 +0100 @@ -0,0 +1,210 @@ +classdef DiffOpTimeDep < scheme.Scheme + properties + grid + order + diffOps + D + H + + blockmatrixDiv + end + + methods + function obj = DiffOpTimeDep(doHand, grid, order, doParam) + % doHand -- may either be a function handle or a cell array of + % function handles for each grid. The function handle(s) + % should be on the form do = doHand(grid, order, ...) + % Additional parameters for each doHand may be provided in + % the doParam input. + % grid -- a multiblock grid + % order -- integer specifying the order of accuracy + % doParam -- may either be a cell array or a cell array of cell arrays + % for each block. If it is a cell array with length equal + % to the number of blocks then each element is sent to the + % corresponding function handle as extra parameters: + % doHand(..., doParam{i}{:}) Otherwise doParam is sent as + % extra parameters to all doHand: doHand(..., doParam{:}) + default_arg('doParam', []) + + [getHand, getParam] = parseInput(doHand, grid, doParam); + obj.grid = grid; + nBlocks = grid.nBlocks(); + + obj.order = order; + + % Create the diffOps for each block + obj.diffOps = cell(1, nBlocks); + for i = 1:nBlocks + h = getHand(i); + p = getParam(i); + if ~iscell(p) + p = {p}; + end + obj.diffOps{i} = h(grid.grids{i}, order, p{:}); + end + + + % Build the norm matrix + H = cell(nBlocks, nBlocks); + for i = 1:nBlocks + H{i,i} = obj.diffOps{i}.H; + end + obj.H = blockmatrix.toMatrix(H); + + + % Build the differentiation matrix + obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + % D = blockmatrix.zero(obj.blockmatrixDiv); + + for i = 1:nBlocks + for j = 1:nBlocks + D{i,j} = @(t) 0; + D{j,i} = @(t) 0; + end + end + + + for i = 1:nBlocks + D{i,i} = @(t)obj.diffOps{i}.D(t); + end + + for i = 1:nBlocks + for j = 1:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = @(t) D{i,i}(t) + ii(t); + D{i,j} = @(t) D{i,j}(t) + ij(t); + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = @(t) D{j,j}(t) + jj(t); + D{j,i} = @(t) D{j,i}(t) + ji(t); + end + end + obj.D = D; + + + function [getHand, getParam] = parseInput(doHand, grid, doParam) + if ~isa(grid, 'multiblock.Grid') + error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); + end + + if iscell(doHand) && length(doHand) == grid.nBlocks() + getHand = @(i)doHand{i}; + elseif isa(doHand, 'function_handle') + getHand = @(i)doHand; + else + error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); + end + + if isempty(doParam) + getParam = @(i){}; + return + end + + if ~iscell(doParam) + getParam = @(i)doParam; + return + end + + % doParam is a non-empty cell-array + + if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) + % doParam is a cell-array of cell-arrays + getParam = @(i)doParam{i}; + return + end + + getParam = @(i)doParam; + end + end + + function ops = splitOp(obj, op) + % Splits a matrix operator into a cell-matrix of matrix operators for + % each grid. + ops = sparse2cell(op, obj.NNN); + end + + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperator(obj, opName, boundary) + switch class(boundary) + case 'cell' + localOpName = [opName '_' boundary{2}]; + blockId = boundary{1}; + localOp = obj.diffOps{blockId}.(localOpName); + + div = {obj.blockmatrixDiv{1}, size(localOp,2)}; + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + return + case 'multiblock.BoundaryGroup' + op = sparse(size(obj.D,1),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperator(opName, boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + end + end + + % Creates the closure and penalty matrix for a given boundary condition, + % boundary -- the name of the boundary on the form {id,name} where + % id is the number of a block and name is the name of a + % boundary of that block example: {1,'s'} or {3,'w'}. It + % can also be a boundary group + function [closure, penalty] = boundary_condition(obj, boundary, type) + switch class(boundary) + case 'cell' + [closure, penalty] = obj.singleBoundaryCondition(boundary, type); + case 'multiblock.BoundaryGroup' + nBlocks = obj.grid.nBlocks(); + %[n,m] = size(obj.D); + %closure = sparse(n,m); + %penalty = sparse(n,0); + % closure =@(t)0; + % penalty = @(t)0; + for i = 1:nBlocks + for j = 1:nBlocks + closure{j,i} = @(t)0; + penalty{j,i} = @(t)0; + end + end + + + for i = 1:length(boundary) + [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type); + closure{i,i} = @(t)closure{i,i}(t) + closurePart(t); + penalty{i,i} = @(t)penalty{i,i}(t) + penaltyPart(t); + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function [blockClosure, blockPenalty] = singleBoundaryCondition(obj, boundary, type) + I = boundary{1}; + name = boundary{2}; + % Get the closure and penaly matrices + [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('not implemented') + end + + % Size returns the number of degrees of freedom + function N = size(obj) + N = 0; + for i = 1:length(obj.diffOps) + N = N + obj.diffOps{i}.size(); + end + end + end +end