Mercurial > repos > public > sbplib
diff +multiblock/DiffOp.m @ 498:324c927d8b1d feature/quantumTriangles
chnaged sbp interfacein 1d among many things
author | Ylva Rydin <ylva.rydin@telia.com> |
---|---|
date | Fri, 03 Mar 2017 16:19:41 +0100 |
parents | 30ff8879162e |
children | 111fcbcff2e9 |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Sat Feb 25 12:44:01 2017 +0100 +++ b/+multiblock/DiffOp.m Fri Mar 03 16:19:41 2017 +0100 @@ -5,12 +5,12 @@ diffOps D H - + blockmatrixDiv end - + methods - function obj = DiffOp(doHand, grid, order, doParam) + function obj = DiffOp(doHand, grid, order, doParam,timeDep) % 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, ...) @@ -25,13 +25,13 @@ % doHand(..., doParam{i}{:}) Otherwise doParam is sent as % extra parameters to all doHand: doHand(..., doParam{:}) default_arg('doParam', []) - + [getHand, getParam] = parseInput(doHand, grid, doParam); - + nBlocks = grid.nBlocks(); - + obj.order = order; - + % Create the diffOps for each block obj.diffOps = cell(1, nBlocks); for i = 1:nBlocks @@ -42,48 +42,73 @@ 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 - D{i,i} = obj.diffOps{i}.D; - end - - for i = 1:nBlocks - for j = 1:nBlocks - intf = grid.connections{i,j}; - if isempty(intf) - continue + switch timeDep + case {'n','no','N','No'} + obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + D = blockmatrix.zero(obj.blockmatrixDiv); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; end - - - [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); - D{i,i} = D{i,i} + ii; - D{i,j} = D{i,j} + ij; - - [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); - D{j,j} = D{j,j} + jj; - D{j,i} = D{j,i} + ji; - 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} = D{i,i} + ii; + D{i,j} = D{i,j} + ij; + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = D{j,j} + jj; + D{j,i} = D{j,i} + ji; + end + end + obj.D = blockmatrix.toMatrix(D); + case {'y','yes','Y','Yes'} + 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} = ij; + + [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} = ji; + end + end + obj.D = D; end - obj.D = blockmatrix.toMatrix(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') @@ -91,41 +116,41 @@ 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 - + function op = getBoundaryOperator(obj, op, boundary) if iscell(boundary) localOpName = [op '_' 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; @@ -135,7 +160,7 @@ % Boundary är en sträng med en boundary group i. 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 @@ -143,10 +168,10 @@ function [closure, penalty] = boundary_condition(obj, boundary, type) I = boundary{1}; name = boundary{2}; - + % Get the closure and penaly matrices [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); - + % Expand to matrix for full domain. div = obj.blockmatrixDiv; if ~iscell(blockClosure) @@ -160,7 +185,7 @@ closure{i} = blockmatrix.toMatrix(temp); end end - + div{2} = size(blockPenalty, 2); % Penalty is a column vector if ~iscell(blockPenalty) p = blockmatrix.zero(div); @@ -174,11 +199,11 @@ end end end - + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - + end - + % Size returns the number of degrees of freedom function N = size(obj) N = 0;