Mercurial > repos > public > sbplib
changeset 592:4422c4476650 feature/utux2D
Merge with feature/grids
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Mon, 11 Sep 2017 14:17:15 +0200 |
| parents | 39554f2de783 (current diff) 97b9a0023d38 (diff) |
| children | a2ddaccf5fd1 |
| files | +multiblock/BoundaryGroupTest.m textTable.m |
| diffstat | 47 files changed, 1454 insertions(+), 335 deletions(-) [+] |
line wrap: on
line diff
--- a/+blockmatrix/fromMatrixTest.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/fromMatrixTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -46,6 +46,20 @@ [23 5; 4 6; 10 12; 11 18], [7 14 16; 13 20 22; 19 21 3; 25 2 9]; }; }, + { + { + magic(3), + {[1 0 2],[1 2 0]} + }, + mat2cell(magic(3),[1 0 2],[1 2 0]) + }, + { + { + zeros(0,1), + {0,1}, + }, + {zeros(0,1)} + }, }; for i = 1:length(cases) out = convertToFull(blockmatrix.fromMatrix(cases{i}{1}{:}));
--- a/+blockmatrix/getDivision.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/getDivision.m Mon Sep 11 14:17:15 2017 +0200 @@ -16,7 +16,7 @@ m = zeros(1,size(C,2)); for j = 1:size(C,2) for i = 1:size(C,1) - if isempty(C{i,j}) + if isNullMatrix(C{i,j}) continue end m(j) = size(C{i,j},2); @@ -28,10 +28,15 @@ n = zeros(1,size(C,1)); for i = 1:size(C,1) for j = 1:size(C,2) - if isempty(C{i,j}) + if isNullMatrix(C{i,j}) continue end n(i) = size(C{i,j},1); end end -end \ No newline at end of file +end + +function b = isNullMatrix(A) + [n, m] = size(A); + b = n == 0 && m == 0; +end
--- a/+blockmatrix/getDivisionTest.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/getDivisionTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -56,6 +56,18 @@ }, {[2 1], 2} }, + { + {zeros(3,0)}, + {3, 0}, + }, + { + {zeros(3,0), zeros(3,0)}, + {3, [0, 0]}, + }, + { + {zeros(3,0); zeros(2,0)}, + {[3 2],0}, + }, }; for i = 1:length(cases)
--- a/+blockmatrix/isBlockmatrix.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/isBlockmatrix.m Mon Sep 11 14:17:15 2017 +0200 @@ -4,7 +4,7 @@ return end - % Make sure all blocks are numerica matrices + % Make sure all blocks are numerical matrices for i = 1:length(bm) if ~isnumeric(bm{i}) b = false;
--- a/+blockmatrix/isDivision.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/isDivision.m Mon Sep 11 14:17:15 2017 +0200 @@ -21,11 +21,11 @@ function b = isDivisionVector(v) if isempty(v) - b = false; + b = true; return end - if any(v <= 0) + if any(v < 0) b = false; return end
--- a/+blockmatrix/isDivisionTest.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/isDivisionTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -4,14 +4,16 @@ function testIsDivision(testCase) cases = { + {[1 2] ,false}, % Must be a cell array + {{[1 2 3]} ,false}, % Must have two vectors + {{[],[]}, true} % No blocks is a valid blockmatrix + {{[1 2],[]} ,true}, + {{[],[1 2]} ,true}, {{[2 2 2],[1 2]} ,true}, - {{[1 2],[1 0]} ,false}, - {{[0 2],[1 1]} ,false}, - {{[1 2],[]} ,false}, + {{[1 2],[1 0]} ,true}, + {{[0 2],[1 1]} ,true}, {{[1 2],[1]} ,true}, {{[1 2],[1], [1 2 3]} ,false}, - {{[1 2 3]} ,false}, - {[1 2] ,false}, }; for i = 1:length(cases)
--- a/+blockmatrix/toMatrixTest.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/toMatrixTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -51,6 +51,22 @@ [2 2 1; 2 1 2], }, + { + {zeros(0,0)}, + [], + }, + { + {zeros(3,0), zeros(3,0)}, + zeros(3,0), + }, + { + {zeros(3,0); zeros(2,0)}, + zeros(5,0), + }, + { + {zeros(0,3), zeros(0,2)}, + zeros(0,5), + }, }; for i = 1:length(cases)
--- a/+blockmatrix/zero.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/zero.m Mon Sep 11 14:17:15 2017 +0200 @@ -1,6 +1,6 @@ % Creates a block matrix according to the division with zeros everywhere. function bm = zero(div) - if ~blockmatrix.isDivision(div); + if ~blockmatrix.isDivision(div) error('div is not a valid division'); end
--- a/+blockmatrix/zeroTest.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+blockmatrix/zeroTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -32,6 +32,22 @@ {[1 2],[2 1]}, {[0 0],[0];[0 0; 0 0],[0; 0]}; }, + { + {[3],[0]}, + {zeros(3,0)}, + }, + + { + {[0],[3]}, + {zeros(0,3)}, + }, + { + {[0 2],[0 3]}, + { + zeros(0,0), zeros(0,3); + zeros(2,0), zeros(2,3); + }, + }, }; for i = 1:length(cases)
--- a/+grid/Cartesian.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+grid/Cartesian.m Mon Sep 11 14:17:15 2017 +0200 @@ -15,6 +15,8 @@ obj.d = length(varargin); for i = 1:obj.d + assert(isvector(varargin{i}), 'Coordinate inputs must be a vectors.') + obj.x{i} = varargin{i}; obj.m(i) = length(varargin{i}); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Empty.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,53 @@ +classdef Empty < grid.Grid & grid.Structured + properties + dim + end + + methods + function obj = Empty(D) + obj.dim = D; + end + % n returns the number of points in the grid + function o = N(obj) + o = 0; + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.dim; + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = sparse(0,obj.dim); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + error('Restrict does not make sense for an empty grid') + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + error('Project does not make sense for an empty grid') + end + + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. + function bs = getBoundaryNames(obj) + bs = {}; + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = sparse(0,obj.dim-1); + end + + function h = scaling(obj) + h = 1; + end + + function s = size(obj) + s = zeros(1, obj.dim); + end + end +end \ No newline at end of file
--- a/+grid/Grid.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+grid/Grid.m Mon Sep 11 14:17:15 2017 +0200 @@ -16,7 +16,7 @@ % Projects the grid function gf on obj to the grid g. gf = projectFunc(obj, gf, g) - % Return the names of all boundaries in this grid. + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. bs = getBoundaryNames(obj) % Return coordinates for the given boundary
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryIdentifier.txt Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,5 @@ +A grid.boundaryIdentifier identifies a boundary of a grid. +For a Cartesian grid it is simply 's', 'n', 'w', 'e'. +For a multiblock grid it might be something like {1, 's'}. +For some other grid it will be up to the developer to chose +a good way to identify boundaries.
--- a/+grid/evalOn.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+grid/evalOn.m Mon Sep 11 14:17:15 2017 +0200 @@ -27,7 +27,11 @@ x = num2cell(g.points()); % Find the number of components - x0 = x(1,:); + if size(x,1) ~= 0 + x0 = x(1,:); + else + x0 = num2cell(ones(1,size(x,2))); + end f0 = func(x0{:}); k = length(f0);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/gridFunction.txt Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,1 @@ +TODO: Add documentation for gridfunctions here \ No newline at end of file
--- a/+multiblock/BoundaryGroup.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+multiblock/BoundaryGroup.m Mon Sep 11 14:17:15 2017 +0200 @@ -1,30 +1,10 @@ % BoundaryGroup defines a boundary grouping in a multiblock grid. -classdef BoundaryGroup - properties - blockIDs - names - end - +% It workds like a cell array and collects boundary identifiers +% Within the multiblock package a BoundaryGroup is a valid boundary identifier as well. +classdef BoundaryGroup < Cell methods - function obj = BoundaryGroup(varargin) - % Input arguemnts are arbitrary number or 1x2 cell arrays - % representing each boundary in the group. - % The 1st element of the cell array is an integer defining which grid it belongs to. - % The 2nd element of the cell array is the name of the boundary within the block. - % - % Ex: - % bg = multiblock.BoundaryGroup({1,'n'},{1,'s'},{2,'s'}) - - - obj.blockIDs = []; - obj.names = {}; - for i = 1:length(varargin) - if ~iscell(varargin{i}) || ~all(size(varargin{i}) == [1 2]) - error('multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', 'Inputs must be 1x2 cell arrays'); - end - obj.blockIDs(i) = varargin{i}{1}; - obj.names{i} = varargin{i}{2}; - end + function obj = BoundaryGroup(data) + obj = obj@Cell(data); end function display(obj, name) @@ -33,19 +13,7 @@ disp([name, ' =']) disp(' ') - if length(obj.names) == 1 - fprintf(' {}\n\n') - return - end - - fprintf(' {') - - fprintf('%d:%s', obj.blockIDs(1), obj.names{1}) - for i = 2:length(obj.names) - fprintf(', %d:%s', obj.blockIDs(i), obj.names{i}); - end - - fprintf('}\n\n') + fprintf(' BoundaryGroup%s\n\n', toString(obj.data)); end end -end \ No newline at end of file +end
--- a/+multiblock/BoundaryGroupTest.m Mon Sep 11 14:12:54 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -function tests = BoundaryGroupTest() - tests = functiontests(localfunctions); -end - -function testCreation(testCase) - in = {{3,'n'},{2,'hoho'},{1,'s'}}; - - blockIDs = [3 2 1]; - names = {'n', 'hoho', 's'}; - - bg = multiblock.BoundaryGroup(in{:}); - testCase.verifyEqual(bg.blockIDs, blockIDs); - testCase.verifyEqual(bg.names, names); -end - -function testInputError(testCase) - in = { - {'n', 's'}, - {{3,'n'},{2,2,'hoho'},{1,'s'}}, - }; - - out = { - 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', - 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', - }; - - for i = 1:length(in) - testCase.verifyError(@()multiblock.BoundaryGroup(in{i}{:}), out{i}); - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Contour.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,65 @@ +classdef Contour < handle + properties + grid + contours + nContours + + ZData + CData + + end + + methods + function obj = Contour(g, gf, nContours) + obj.grid = g; + obj.nContours = nContours; + + coords = obj.grid.points(); + X = obj.grid.funcToPlotMatrices(coords(:,1)); + Y = obj.grid.funcToPlotMatrices(coords(:,2)); + + V = obj.grid.funcToPlotMatrices(gf); + + + holdState = ishold(); + hold on + + contours = {1, obj.grid.nBlocks}; + for i = 1:obj.grid.nBlocks + [~, contours{i}] = contour(X{i}, Y{i}, V{i},obj.nContours); + contours{i}.LevelList = contours{1}.LevelList; + end + + if holdState == false + hold off + end + + obj.contours = [contours{:}]; + + obj.ZData = gf; + obj.CData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.contours, propertyName, propertyValue); + end + + function obj = set.ZData(obj, gf) + obj.ZData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.contours(i).ZData = V{i}; + end + end + + function obj = set.CData(obj, gf) + obj.CData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.contours(i).CData = V{i}; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Def.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,101 @@ +classdef Def + properties + nBlocks + blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation + blockNames + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + end + + methods + % Defines a multiblock setup for transfinite interpolation blocks + % TODO: How to bring in plotting of points? + function obj = Def(blockMaps, connections, boundaryGroups, blockNames) + default_arg('boundaryGroups', struct()); + default_arg('blockNames',{}); + + nBlocks = length(blockMaps); + + obj.nBlocks = nBlocks; + + obj.blockMaps = blockMaps; + + assert(all(size(connections) == [nBlocks, nBlocks])); + obj.connections = connections; + + + if isempty(blockNames) + obj.blockNames = cell(1, nBlocks); + for i = 1:length(blockMaps) + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == nBlocks); + obj.blockNames = blockNames; + end + + obj.boundaryGroups = boundaryGroups; + end + + function g = getGrid(obj, varargin) + ms = obj.getGridSizes(varargin{:}); + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(2,2); + end + axis equal + return + end + + if gridLines + ms = obj.getGridSizes(varargin{:}); + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockMaps{i}, labels{i}); + end + + axis equal + end + end + + methods (Abstract) + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + ms = getGridSizes(obj, varargin) + % end + end + +end + +
--- a/+multiblock/DiffOp.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+multiblock/DiffOp.m Mon Sep 11 14:17:15 2017 +0200 @@ -120,27 +120,54 @@ 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); + % 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 - else - % Boundary är en sträng med en boundary group i. + 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'} + % 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' + [n,m] = size(obj.D); + closure = sparse(n,m); + penalty = sparse(n,0); + for i = 1:length(boundary) + [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type); + closure = closure + closurePart; + penalty = [penalty, penaltyPart]; + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function [closure, penalty] = singleBoundaryCondition(obj, boundary, type) I = boundary{1}; name = boundary{2}; @@ -177,7 +204,7 @@ end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - + error('not implemented') end % Size returns the number of degrees of freedom
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/EmptyGrid.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,47 @@ +classdef EmptyGrid < grid.Empty & multiblock.Grid + methods + function obj = EmptyGrid(D) + obj@multiblock.Grid({},cell(0,0)); + obj@grid.Empty(D); + end + + % n returns the number of points in the grid + function o = N(obj) + o = N@grid.Empty(obj); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = D@grid.Empty(obj); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = points@grid.Empty(obj); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gf = restrictFunc@grid.Empty(obj); + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + gf = projectFunc@grid.Empty(obj); + end + + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. + function bs = getBoundaryNames(obj) + bs = getBoundaryNames@grid.Empty(obj); + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = getBoundary@grid.Empty(name); + end + + function s = size(obj) + s = size@multiblock.Grid(obj); + end + end +end
--- a/+multiblock/Grid.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+multiblock/Grid.m Mon Sep 11 14:17:15 2017 +0200 @@ -9,16 +9,20 @@ % General multiblock grid methods - - % grids -- cell array of N grids - % connections -- NxN upper triangular cell matrix. connections{i,j} - % specifies the connection between block i and j. If - % it's empty there is no connection otherwise it's a 2 - % -cell-vector with strings naming the boundaries to be - % connected. (inverted coupling?) - %% Should we have boundary groups at all? maybe it can be handled in a - %% cleaner way outside of the class. Maybe the grid class doesn't need to care at all. All boundaryGroup interaction is in DiffOp? + % grids -- cell array of N grids + % connections -- NxN upper triangular cell matrix. connections{i,j} + % specifies the connection between block i and j. If + % it's empty there is no connection otherwise it's a 2 + % -cell-vector with strings naming the boundaries to be + % connected. (inverted coupling?) + % boundaryGroups -- A struct of BoundaryGroups. The field names of the + % struct are the names of each boundary group. + % The boundary groups can be used to collect block + % boundaries into physical boundaries to simplify + % getting boundary operators and setting boundary conditions function obj = Grid(grids, connections, boundaryGroups) + default_arg('boundaryGroups', struct()); + assertType(grids, 'cell') obj.grids = grids; obj.connections = connections; @@ -27,7 +31,7 @@ obj.nPoints = obj.nPoints + grids{i}.N(); end - % if iscell(boundaryGroups) + obj.boundaryGroups = boundaryGroups; end function n = size(obj) @@ -42,7 +46,7 @@ % Ns returns the number of points in each sub grid as a vector function o = Ns(obj) ns = zeros(1,obj.nBlocks); - for i = 1:obj.nBlocks; + for i = 1:obj.nBlocks ns(i) = obj.grids{i}.N(); end o = ns; @@ -59,7 +63,7 @@ % points returns a n x d matrix containing the coordinates for all points. function X = points(obj) - X = []; + X = sparse(0,obj.D()); for i = 1:length(obj.grids) X = [X; obj.grids{i}.points]; end @@ -76,12 +80,19 @@ N(i) = obj.grids{i}.N(); end - gfs = mat2cell(gf, N, 1); + gfs = blockmatrix.fromMatrix(gf, {N,1}); end + % TODO: Split op? + % Should the method to split an operator be moved here instead of being in multiblock.DiffOp? + % Converts a gridfunction to a set of plot matrices % Takes a grid function and and a structured grid. function F = funcToPlotMatrices(obj, gf) + % TODO: This method should problably not be here. + % The funcToPlotMatrix uses .size poperty of the grids + % Which doesn't always exist for all types of grids. + % It's only valid for structured grids gfs = obj.splitFunc(gf); F = cell(1, obj.nBlocks()); @@ -121,13 +132,43 @@ end + % Find all non interface boundaries of all blocks. + % Return their grid.boundaryIdentifiers in a cell array. function bs = getBoundaryNames(obj) - bs = []; + bs = {}; + for i = 1:obj.nBlocks() + candidates = obj.grids{i}.getBoundaryNames(); + for j = 1:obj.nBlocks() + if ~isempty(obj.connections{i,j}) + candidates = setdiff(candidates, obj.connections{i,j}{1}); + end + + if ~isempty(obj.connections{j,i}) + candidates = setdiff(candidates, obj.connections{j,i}{2}); + end + end + + for k = 1:length(candidates) + bs{end+1} = {i, candidates{k}}; + end + end end - % Return coordinates for the given boundary - function b = getBoundary(obj, name) - b = []; + % Return coordinates for the given boundary/boundaryGroup + function b = getBoundary(obj, boundary) + switch class(boundary) + case 'cell' + I = boundary{1}; + name = boundary{2}; + b = obj.grids{I}.getBoundary(name); + case 'multiblock.BoundaryGroup' + b = sparse(0,obj.D()); + for i = 1:length(boundary) + b = [b; obj.getBoundary(boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + end end end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Line.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,47 @@ +classdef Line < handle + properties + grid + lines + + YData + end + + methods + function obj = Line(g, gf) + assertType(g, 'multiblock.Grid') + obj.grid = g; + + X = obj.grid.splitFunc(obj.grid.points()); + Y = obj.grid.splitFunc(gf); + + holdState = ishold(); + hold on + + lines = cell(1, obj.grid.nBlocks); + for i = 1:obj.grid.nBlocks + lines{i} = plot(X{i}, Y{i}); + end + + if holdState == false + hold off + end + + obj.lines = [lines{:}]; + + obj.YData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.lines, propertyName, propertyValue); + end + + function obj = set.YData(obj, gf) + obj.YData = gf; + + Y = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.lines(i).YData = Y{i}; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Surface.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,62 @@ +classdef Surface < handle + properties + grid + surfs + + ZData + CData + + end + + methods + function obj = Surface(g, gf) + obj.grid = g; + + coords = obj.grid.points(); + X = obj.grid.funcToPlotMatrices(coords(:,1)); + Y = obj.grid.funcToPlotMatrices(coords(:,2)); + + V = obj.grid.funcToPlotMatrices(gf); + + + holdState = ishold(); + hold on + + surfs = cell(1, obj.grid.nBlocks); + for i = 1:obj.grid.nBlocks + surfs{i} = surf(X{i}, Y{i}, V{i}); + end + + if holdState == false + hold off + end + + obj.surfs = [surfs{:}]; + + obj.ZData = gf; + obj.CData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.surfs, propertyName, propertyValue); + end + + function obj = set.ZData(obj, gf) + obj.ZData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.surfs(i).ZData = V{i}; + end + end + + function obj = set.CData(obj, gf) + obj.CData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.surfs(i).CData = V{i}; + end + end + end +end
--- a/+parametrization/Ti.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+parametrization/Ti.m Mon Sep 11 14:17:15 2017 +0200 @@ -182,6 +182,15 @@ obj = parametrization.Ti(g1,g2,g3,g4); end + function obj = rectangle(a, b) + p1 = a; + p2 = [b(1), a(2)]; + p3 = b; + p4 = [a(1), b(2)]; + + obj = parametrization.Ti.points(p1,p2,p3,p4); + end + % Like the constructor but allows inputing line curves as 2-cell arrays: % example: parametrization.Ti.linesAndCurves(g1, g2, {a, b} g4) function obj = linesAndCurves(C1, C2, C3, C4)
--- a/+scheme/Beam.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+scheme/Beam.m Mon Sep 11 14:17:15 2017 +0200 @@ -90,6 +90,9 @@ gamm = obj.gamm; delt = obj.delt; + + % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ... + switch type case {'dn', 'clamped'} % Dirichlet-neumann boundary condition alpha = obj.alpha;
--- a/+scheme/LaplaceCurvilinear.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+scheme/LaplaceCurvilinear.m Mon Sep 11 14:17:15 2017 +0200 @@ -7,42 +7,48 @@ order % Order accuracy for the approximation - D % non-stabalized scheme operator - M % Derivative norm - a,b + a,b % Parameters of the operator + + + % Inner products and operators for physical coordinates + D % Laplace operator + H, Hi % Inner product + e_w, e_e, e_s, e_n + d_w, d_e, d_s, d_n % Normal derivatives at the boundary + H_w, H_e, H_s, H_n % Boundary inner products + Dx, Dy % Physical derivatives + M % Gradient inner product + + % Metric coefficients J, Ji a11, a12, a22 + x_u + x_v + y_u + y_v - H % Discrete norm - Hi + % Inner product and operators for logical coordinates H_u, H_v % Norms in the x and y directions + Hi_u, Hi_v Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. - Hi_u, Hi_v Hiu, Hiv - e_w, e_e, e_s, e_n du_w, dv_w du_e, dv_e du_s, dv_s du_n, dv_n gamm_u, gamm_v lambda - - Dx, Dy % Physical derivatives - - x_u - x_v - y_u - y_v end methods % Implements a*div(b*grad(u)) as a SBP scheme + % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) + function obj = LaplaceCurvilinear(g ,order, a, b, opSet) default_arg('opSet',@sbp.D2Variable); default_arg('a', 1); default_arg('b', 1); - if b ~=1 error('Not implemented yet') end @@ -58,7 +64,8 @@ h_u = h(1); h_v = h(2); - % Operators + + % 1D operators ops_u = opSet(m_u, {0, 1}, order); ops_v = opSet(m_v, {0, 1}, order); @@ -83,10 +90,30 @@ d1_l_v = ops_v.d1_l; d1_r_v = ops_v.d1_r; + + % Logical operators Du = kr(D1_u,I_v); Dv = kr(I_u,D1_v); + obj.Hu = kr(H_u,I_v); + obj.Hv = kr(I_u,H_v); + obj.Hiu = kr(Hi_u,I_v); + obj.Hiv = kr(I_u,Hi_v); - % Metric derivatives + e_w = kr(e_l_u,I_v); + e_e = kr(e_r_u,I_v); + e_s = kr(I_u,e_l_v); + e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (e_w'*Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (e_e'*Dv)'; + obj.du_s = (e_s'*Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (e_n'*Du)'; + obj.dv_n = kr(I_u,d1_r_v); + + + % Metric coefficients coords = g.points(); x = coords(:,1); y = coords(:,2); @@ -102,8 +129,14 @@ a22 = 1./J .* (x_u.^2 + y_u.^2); lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); + obj.x_u = x_u; + obj.x_v = x_v; + obj.y_u = y_u; + obj.y_v = y_v; + + % Assemble full operators - L_12 = spdiags(a12, 0, m_tot, m_tot); + L_12 = spdiag(a12); Duv = Du*L_12*Dv; Dvu = Dv*L_12*Du; @@ -123,31 +156,55 @@ Dvv(p,p) = D; end - obj.H = kr(H_u,H_v); - obj.Hi = kr(Hi_u,Hi_v); - obj.Hu = kr(H_u,I_v); - obj.Hv = kr(I_u,H_v); - obj.Hiu = kr(Hi_u,I_v); - obj.Hiv = kr(I_u,Hi_v); + + % Physical operators + obj.J = spdiag(J); + obj.Ji = spdiag(1./J); + + obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); + obj.H = obj.J*kr(H_u,H_v); + obj.Hi = obj.Ji*kr(Hi_u,Hi_v); + + obj.e_w = e_w; + obj.e_e = e_e; + obj.e_s = e_s; + obj.e_n = e_n; + + %% normal derivatives + I_w = ind(1,:); + I_e = ind(end,:); + I_s = ind(:,1); + I_n = ind(:,end); - obj.e_w = kr(e_l_u,I_v); - obj.e_e = kr(e_r_u,I_v); - obj.e_s = kr(I_u,e_l_v); - obj.e_n = kr(I_u,e_r_v); - obj.du_w = kr(d1_l_u,I_v); - obj.dv_w = (obj.e_w'*Dv)'; - obj.du_e = kr(d1_r_u,I_v); - obj.dv_e = (obj.e_e'*Dv)'; - obj.du_s = (obj.e_s'*Du)'; - obj.dv_s = kr(I_u,d1_l_v); - obj.du_n = (obj.e_n'*Du)'; - obj.dv_n = kr(I_u,d1_r_v); + a11_w = spdiag(a11(I_w)); + a12_w = spdiag(a12(I_w)); + a11_e = spdiag(a11(I_e)); + a12_e = spdiag(a12(I_e)); + a22_s = spdiag(a22(I_s)); + a12_s = spdiag(a12(I_s)); + a22_n = spdiag(a22(I_n)); + a12_n = spdiag(a12(I_n)); + + s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); + s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); + s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); + s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); - obj.x_u = x_u; - obj.x_v = x_v; - obj.y_u = y_u; - obj.y_v = y_v; + obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))'; + obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))'; + obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))'; + obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))'; + + obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; + obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; + %% Boundary inner products + obj.H_w = H_v*spdiag(s_w); + obj.H_e = H_v*spdiag(s_e); + obj.H_s = H_u*spdiag(s_s); + obj.H_n = H_u*spdiag(s_n); + + % Misc. obj.m = m; obj.h = [h_u h_v]; obj.order = order; @@ -155,17 +212,11 @@ obj.a = a; obj.b = b; - obj.J = spdiags(J, 0, m_tot, m_tot); - obj.Ji = spdiags(1./J, 0, m_tot, m_tot); obj.a11 = a11; obj.a12 = a12; obj.a22 = a22; - obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); obj.lambda = lambda; - obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; - obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; - obj.gamm_u = h_u*ops_u.borrowing.M.d1; obj.gamm_v = h_v*ops_v.borrowing.M.d1; end @@ -182,62 +233,34 @@ default_arg('type','neumann'); default_arg('parameter', []); - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); + [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); switch type % Dirichlet boundary condition case {'D','d','dirichlet'} - % v denotes the solution in the neighbour domain tuning = 1.2; % tuning = 20.2; - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); - a_n = spdiag(coeff_n); - a_t = spdiag(coeff_t); - - F = (s * a_n * d_n' + s * a_t*d_t')'; - - u = obj; + b1 = gamm*obj.lambda./obj.a11.^2; + b2 = gamm*obj.lambda./obj.a22.^2; - b1 = gamm*u.lambda./u.a11.^2; - b2 = gamm*u.lambda./u.a22.^2; + tau1 = tuning * spdiag(-1./b1 - 1./b2); + tau2 = 1; - tau = -1./b1 - 1./b2; - tau = tuning * spdiag(tau); - sig1 = 1; + tau = (tau1*e + tau2*d)*H_b; - penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; - - closure = obj.Ji*obj.a * penalty_parameter_1*e'; - penalty = -obj.Ji*obj.a * penalty_parameter_1; + closure = obj.a*obj.Hi*tau*e'; + penalty = -obj.a*obj.Hi*tau; % Neumann boundary condition case {'N','n','neumann'} - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = (a_n * d_n' + a_t*d_t')'; - - tau1 = -s; + tau1 = -1; tau2 = 0; - tau = obj.a * obj.Ji*(tau1*e + tau2*d); - - closure = halfnorm_inv*tau*d'; - penalty = -halfnorm_inv*tau; + tau = (tau1*e + tau2*d)*H_b; - % Characteristic boundary condition - case {'characteristic', 'char', 'c'} - default_arg('parameter', 1); - beta = parameter; + closure = obj.a*obj.Hi*tau*d'; + penalty = -obj.a*obj.Hi*tau; - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative - - tau = -obj.a * 1/beta*obj.Ji*e; - - closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; - closure{2} = halfnorm_inv*tau*beta*d'; - penalty = -halfnorm_inv*tau; % Unknown, boundary condition otherwise @@ -250,16 +273,8 @@ % v denotes the solution in the neighbour domain tuning = 1.2; % tuning = 20.2; - [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); - [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - - a_n_u = spdiag(coeff_n_u); - a_t_u = spdiag(coeff_t_u); - a_n_v = spdiag(coeff_n_v); - a_t_v = spdiag(coeff_t_v); - - F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; - F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; + [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); + [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); u = obj; v = neighbour_scheme; @@ -269,24 +284,25 @@ b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; - tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - tau = tuning * spdiag(tau); - sig1 = 1/2; - sig2 = -1/2; + tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); + tau1 = tuning * spdiag(tau1); + tau2 = 1/2; - penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); - penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; + sig1 = -1/2; + sig2 = 0; + tau = (e_u*tau1 + tau2*d_u)*H_b_u; + sig = (sig1*e_u + sig2*d_u)*H_b_u; - closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); - penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); + closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); + penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary. % The right boundary is considered the positive boundary % % I -- the indecies of the boundary points in the grid matrix - function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary) + function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) % gridMatrix = zeros(obj.m(2),obj.m(1)); % gridMatrix(:) = 1:numel(gridMatrix); @@ -296,58 +312,32 @@ switch boundary case 'w' e = obj.e_w; - d_n = obj.du_w; - d_t = obj.dv_w; - s = -1; - + d = obj.d_w; + H_b = obj.H_w; I = ind(1,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); case 'e' e = obj.e_e; - d_n = obj.du_e; - d_t = obj.dv_e; - s = 1; - + d = obj.d_e; + H_b = obj.H_e; I = ind(end,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); case 's' e = obj.e_s; - d_n = obj.dv_s; - d_t = obj.du_s; - s = -1; - + d = obj.d_s; + H_b = obj.H_s; I = ind(:,1)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); case 'n' e = obj.e_n; - d_n = obj.dv_n; - d_t = obj.du_n; - s = 1; - + d = obj.d_n; + H_b = obj.H_n; I = ind(:,end)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); otherwise error('No such boundary: boundary = %s',boundary); end switch boundary case {'w','e'} - halfnorm_inv_n = obj.Hiu; - halfnorm_inv_t = obj.Hiv; - halfnorm_t = obj.Hv; gamm = obj.gamm_u; case {'s','n'} - halfnorm_inv_n = obj.Hiv; - halfnorm_inv_t = obj.Hiu; - halfnorm_t = obj.Hu; gamm = obj.gamm_v; end end @@ -355,7 +345,5 @@ function N = size(obj) N = prod(obj.m); end - - end -end \ No newline at end of file +end
--- a/+scheme/Scheme.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+scheme/Scheme.m Mon Sep 11 14:17:15 2017 +0200 @@ -25,9 +25,12 @@ % neighbour_boundary is a string specifying which boundary to % interface to. % penalty may be a cell array if there are several penalties with different weights - [closure, penalty] = boundary_condition(obj,boundary,type) + [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % TODO: op = getBoundaryOperator()?? + % makes sense to have it available through a method instead of random properties + % Returns the number of degrees of freedom. N = size(obj) end
--- a/+time/Rungekutta4SecondOrder.m Mon Sep 11 14:12:54 2017 +0200 +++ b/+time/Rungekutta4SecondOrder.m Mon Sep 11 14:17:15 2017 +0200 @@ -15,7 +15,19 @@ methods - % Solves u_tt = Du + Eu_t + S + % Solves u_tt = Du + Eu_t + S by + % Rewriting on first order form: + % w_t = M*w + C(t) + % where + % M = [ + % 0, I; + % D, E; + % ] + % and + % C(t) = [ + % 0; + % S(t) + % ] % D, E, S can either all be constants or all be function handles, % They can also be omitted by setting them equal to the empty matrix. function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) @@ -41,23 +53,15 @@ S = @(t)S; end - I = speye(obj.m); - O = sparse(obj.m,obj.m); - - obj.M = @(t)[ - O, I; - D(t), E(t); - ]; - obj.C = @(t)[ - zeros(obj.m,1); - S(t); - ]; - obj.k = k; obj.t = t0; obj.w = [v0; v0t]; - obj.F = @(w,t)(obj.M(t)*w + obj.C(t)); + % Avoid matrix formulation because it is VERY slow + obj.F = @(w,t)[ + w(obj.m+1:end); + D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); + ]; else default_arg('D', sparse(obj.m, obj.m));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Cell.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,92 @@ +% Cell is a reimplementation of matlabs cell array with the benefit that it is subclassable +% It might be used for giving a typename to a cellarray to increase readability of the code. +classdef Cell + properties + data + end + + methods + function obj = Cell(data) + default_arg('data', {}); + if ~iscell(data) + error('Input argument to Cell must be a cell array'); + end + + obj.data = data; + end + + function str = toString(obj) + str = sprintf('%s%s', class(obj), toString(obj.data)); + end + + function s = size(A) + s = size(A.data); + end + + function b = isempty(A) + b = prod(size(A)) == 0; + end + + function l = length(A) + l = length(A.data); + end + + function ind = end(A,k,n) + ind = builtin('end',A.data, k, n); + end + + function B = transpose(A) + b = A.data.'; + B = callConstructor(A, b); + end + + function B = ctranspose(A) + b = A.data'; + B = callConstructor(A, b); + end + + function A = subsasgn(A, S, B) + a = subsasgn(A.data, S, B); + A = callConstructor(A, a); + end + + function B = subsref(A, S) + switch S(1).type + case '()' + b = subsref(A.data, S(1)); + B = callConstructor(A, b); + if length(S) > 1 + B = subsref(B,S(2:end)); + end + case '{}' + B = subsref(A.data, S); + case '.' + B = builtin('subsref',A, S); + otherwise + error('unreachable'); + end + end + + function C = horzcat(varargin) + dataArray = cell(1, length(varargin)); + + for i = 1:length(varargin) + dataArray{i} = varargin{i}.data; + end + + c = horzcat(dataArray{:}); + C = callConstructor(varargin{1}, c); + end + + function C = vertcat(varargin) + dataArray = cell(1, length(varargin)); + + for i = 1:length(varargin) + dataArray{i} = varargin{i}.data; + end + + c = vertcat(dataArray{:}); + C = callConstructor(varargin{1}, c); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CellTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,279 @@ +function tests = CellTest() + tests = functiontests(localfunctions); +end + +function testSize(testCase) + cases = { + {{}, [0, 0]}, + {{1}, [1, 1]}, + {{1, 2}, [1, 2]}, + {{1; 2}, [2, 1]}, + {{1, 2; 3, 4}, [2,2]}, + }; + + for i = 1:length(cases) + A = Cell(cases{i}{1}); + expected = cases{i}{2}; + + testCase.verifyEqual(size(A),expected); + end +end + +function testLength(testCase) + cases = { + {{}, 0}, + {{1}, 1}, + {{1, 2}, 2}, + {{1; 2}, 2}, + {{1, 2; 3, 4}, 2}, + }; + + for i = 1:length(cases) + A = Cell(cases{i}{1}); + expected = cases{i}{2}; + + testCase.verifyEqual(length(A),expected); + end +end + +function testIsEmpty(testCase) + cases = { + {cell(0,0), true}, + {cell(1,0), true}, + {cell(0,1), true}, + {cell(1,1), false}, + }; + + for i = 1:length(cases) + A = Cell(cases{i}{1}); + expected = cases{i}{2}; + testCase.verifyEqual(isempty(A),expected); + end +end + +function testTranspose(testCase) + testCase.verifyEqual(Cell({1i, 2}).', Cell({1i; 2})); + testCase.verifyEqual(Cell({1i; 2}).', Cell({1i, 2})); +end + +function testCtranspose(testCase) + testCase.verifyEqual(Cell({1i, 2})', Cell({1i; 2})); + testCase.verifyEqual(Cell({1i; 2})', Cell({1i, 2})); +end + +function testRoundIndexWithProperty(testCase) + A = Cell({3,2,1}); + + result = A([1,3]).data; + testCase.verifyEqual(result, {3, 1}); +end + +function testSubAssignmentRound(testCase) + cases = { + % { + % lArray, + % index, + % rhs, + % expectedResult + % }, + { + {}, + 1, + {'a'}, + {'a'}, + }, + { + {1}, + 1, + {'a'}, + {'a'}, + }, + { + {1,2,3}, + 2, + {'a'}, + {1,'a',3}, + }, + { + {1,2,3}, + 2, + [], + {1,3}, + }, + }; + + for i = 1:length(cases) + lArray = Cell(cases{i}{1}); + index = cases{i}{2}; + rhs = cases{i}{3}; + expectedResult = cases{i}{4}; + + lArray(index) = rhs; + + testCase.verifyEqual(lArray.data, expectedResult) + end +end + +function testSubAssignmentCurly(testCase) + cases = { + % { + % lArray, + % index, + % rhs, + % expectedResult + % }, + { + {}, + 1, + 'a', + {'a'}, + }, + { + {1}, + 1, + 'a', + {'a'}, + }, + { + {1,2,3}, + 2, + 'a', + {1,'a',3}, + }, + }; + + for i = 1:length(cases) + lArray = Cell(cases{i}{1}); + index = cases{i}{2}; + rhs = cases{i}{3}; + expectedResult = cases{i}{4}; + + lArray{index} = rhs; + + testCase.verifyEqual(lArray.data, expectedResult) + end +end + +function testIndexreferenceRound(testCase) + cases = { + % { + % array, + % index, + % roundResult + % }, + { + {1,2,3}, + 1, + {1}, + }, + { + {1,3,2}, + 2, + {3}, + }, + { + {1,3,2}, + [1 3], + {1, 2}, + }, + }; + + + for i = 1:length(cases) + array = Cell(cases{i}{1}); + index = cases{i}{2}; + expected = cases{i}{3}; + + result = array(index); + + testCase.verifyTrue(isa(result, 'Cell')); + testCase.verifyEqual(result.data, expected); + end +end + +function testEndIndexing(testCase) + C = Cell({1,2,3}); + + testCase.verifyEqual(C(end), Cell({3})); + testCase.verifyEqual(C{end}, 3); +end + +function testColonIndexing(testCase) + C = Cell({1, 2, 3}); + D = Cell({1; 2; 3}); + + testCase.verifyEqual(C(:), D); + testCase.verifyEqual(D(:), D); +end + +function testIndexreferenceCurly(testCase) + cases = { + % { + % array, + % index, + % curlyResult + % }, + { + {1,2,3}, + 1, + 1 + }, + { + {1,3,2}, + 2, + 3 + }, + }; + + + for i = 1:length(cases) + array = Cell(cases{i}{1}); + index = cases{i}{2}; + expected = cases{i}{3}; + + result = array{index}; + + testCase.verifyEqual(result, expected); + end +end + +function testConcat(testCase) + cases = { + {{},{}}, + {{1},{}}, + {{},{1}}, + {{1},{2}}, + {{1, 2},{3, 4}}, + {{1; 2},{3; 4}}, + }; + + horzCat = { + {}, + {1}, + {1}, + {1,2}, + {1, 2, 3, 4}, + {1, 3; 2, 4}, + }; + + vertCat = { + {}, + {1}, + {1}, + {1; 2}, + {1, 2; 3, 4}, + {1; 2; 3; 4}, + }; + + for i = 1:length(cases) + A = Cell(cases{i}{1}); + B = Cell(cases{i}{2}); + + C_horz = [A, B]; + C_vert = [A; B]; + + testCase.verifyEqual(C_horz.data, horzCat{i}); + testCase.verifyEqual(C_vert.data, vertCat{i}); + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/TextTable.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,122 @@ +classdef TextTable < handle + properties + data + fmtArray + vertDiv + horzDiv + + nCols + nRows + end + + methods + function obj = TextTable(data, vertDiv, horzDiv); + default_arg('vertDiv', []); + default_arg('horzDiv', []); + + + obj.data = data; + obj.vertDiv = vertDiv; + obj.horzDiv = horzDiv; + + [obj.nRows, obj.nCols] = size(data); + obj.fmtArray = cell(size(data)); + obj.formatAll('%s'); + + end + + function formatAll(obj, fmt) + obj.fmtArray(:,:) = {fmt}; + end + + function formatCell(obj, i, j, fmt) + obj.fmtArray{i,j} = fmt; + end + + function formatRow(obj, i, fmt) + obj.fmtArray(i,:) = {fmt}; + end + + function formatColumn(obj, j, fmt) + obj.fmtArray(:,j) = {fmt}; + end + + function widths = getWidths(obj) + strArray = obj.getStringArray(); + + widths = zeros(1,obj.nCols); + for j = 1:obj.nCols + for i = 1:obj.nRows + widths(j) = max(widths(j), length(strArray{i,j})); + end + end + end + + function str = toString(obj) + strArray = obj.getStringArray(); + widths = obj.getWidths(); + + str = ''; + + % First horzDiv + if ismember(0, obj.horzDiv) + str = [str, obj.getHorzDiv(widths)]; + end + + for i = 1:obj.nRows + str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv, obj.horzDiv)]; + + % Interior horzDiv + if ismember(i, obj.horzDiv) + str = [str, obj.getHorzDiv(widths)]; + end + end + end + + function str = getHorzDiv(obj, widths) + str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv, obj.horzDiv); + str(find(' ' == str)) = '-'; + str(find('|' == str)) = '+'; + end + + function strArray = getStringArray(obj) + strArray = cell(size(obj.data)); + + for i = 1:obj.nRows + for j = 1:obj.nCols + strArray{i,j} = sprintf(obj.fmtArray{i,j}, obj.data{i,j}); + end + end + end + end + + methods (Static) + function str = rowToString(strs, widths, vertDiv, horzDiv) + % First vertDiv + if ismember(0, vertDiv) + str = '| '; + else + str = ' '; + end + + % Interior cols + for j = 1:length(strs) - 1 + str = [str, sprintf('%*s ', widths(j), strs{j})]; + + % Interior vertDiv + if ismember(j, vertDiv) + str = [str, '| ']; + end + end + + % Last col + str = [str, sprintf('%*s ', widths(end), strs{end})]; + + if ismember(length(strs), vertDiv) + str = [str, '|']; + end + + str = [str, sprintf('\n')]; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertSymbolic.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,3 @@ +function assertSymbolic(s) + assert(logical(simplify(s))); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertType.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,11 @@ +function assertType(obj, type) + if ~iscell(type) + if ~isa(obj, type) + error('sbplib:assertType:wrongType', '"%s" must have type "%s", found "%s"', inputname(1), type, class(obj)); + end + else + if ~isAnyOf(obj, type) + error('sbplib:assertType:wrongType', '"%s" must be one of the types %s, found "%s"', inputname(1), toString(type), class(obj)); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/callConstructor.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,6 @@ +% Calls the constructor of an object. +% Might be usefull to call the constructor of a subclass object in the superclass +function obj = callConstructor(subclassObj, varargin) + fun = str2func(class(subclassObj)); + obj = fun(varargin{:}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gaussian.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,3 @@ +function z = gaussian(x,x0,d) + z = exp(-norm(x-x0).^2/d^2); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/isAnyOf.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,11 @@ +% Returns true of obj is any of he types in the cell array types +% b = isAnyOf(obj, types) +function b = isAnyOf(obj, types) + for i = 1:length(types) + if isa(obj, types{i}); + b = true; + return + end + end + b = false; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/isEquidistant.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,19 @@ +% Tests if consecutive elements of vector v are euidistant +function b = isEquidistant(v) + if length(v) < 2 + error('sbplib:isEquidistant:inputTooShort', 'Input vector is too short'); + end + + tol = 1e-8; + + d = v(2:end) - v(1:end-1); + err = abs(d - d(1)); + + relErr = err./abs(d); + + I_zero = find(d < tol); + + relErr(I_zero) = err(I_zero); + + b = all(relErr < tol); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/isEquidistantTest.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,31 @@ +function tests = isEquidistantTest() + tests = functiontests(localfunctions); +end + +function testTooShortInput(testCase) + testCase.verifyError(@()isEquidistant([]), 'sbplib:isEquidistant:inputTooShort') +end + +function testCorrectOutput(testCase) + cases = { + % {input, expected}, + {[0,0,0,0,0], true}, + {[1,1,1,1,1], true}, + {[1,2,3,4,5], true}, + {[1,3,4,5], false}, + {[1,2,3,5], false}, + {[1,2,4,5], false}, + {linspace(0,pi, 3), true}, + {linspace(0,1, 4), true}, + {linspace(0,1, 4123), true}, + {linspace(0,sin(1), 123), true}, + }; + + for i = 1:length(cases) + input = cases{i}{1}; + expected = cases{i}{2}; + result = isEquidistant(input); + + testCase.verifyEqual(result,expected); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/logsurf.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,19 @@ +function [srfHandle, cbHandle] = logsurf(X,Y,Z, lim) + absLogZ = log10(abs(Z)); + srfHandle = surf(X,Y,absLogZ); + + cbHandle = colorbar(); + colormap(hot(256)); + ah = gca(); + ah.CLim = lim; + + oldTickLabels = cbHandle.TickLabels; + + newTickLabels = {}; + + for i = 1:length(oldTickLabels) + newTickLabels{i} = sprintf('10^{%s}',oldTickLabels{i}); + end + + cbHandle.TickLabels = newTickLabels; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/skewPart.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,4 @@ +% Returns the skew of A +function S = skewPart(A, tol) + S = 1/2*(A - A'); +end
--- a/spdiag.m Mon Sep 11 14:12:54 2017 +0200 +++ b/spdiag.m Mon Sep 11 14:17:15 2017 +0200 @@ -1,5 +1,10 @@ function A = spdiag(a,i) default_arg('i',0); + + if isrow(a) + a = a'; + end + n = length(a)-abs(i); A = spdiags(a,i,n,n); end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spzeros.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,18 @@ +function S = spzeros(varargin) + switch length(varargin) + case 2 + S = sparse(varargin{1}, varargin{2}); + case 1 + v = varargin{1}; + switch length(v) + case 1 + S = sparse(v,v); + case 2 + S = sparse(v(1), v(2)); + otherwise + error('Input must be either one integer, two integers or a vector with two integers'); + end + otherwise + error('Too many input arguments.'); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/surfSym.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,54 @@ +function h = surfSym(X,Y,Z) + if isvector(X) && isvector(Y) + [X,Y] = meshgrid(X,Y); + end + + V_nodes = [X(:), Y(:), Z(:)]; + + X_centers = neighbourMean(X); + Y_centers = neighbourMean(Y); + Z_centers = neighbourMean(Z); + V_centers = [X_centers(:), Y_centers(:), Z_centers(:)]; + + + N = prod(size(X)); + nodeIndecies = reshape(1:N, size(X)); + centerIndecies = reshape(N+(1:prod(size(X)-[1,1])), size(X) - [1,1]); + + % figure() + % h = line(V_nodes(:,1),V_nodes(:,2),V_nodes(:,3)); + % h.LineStyle = 'none'; + % h.Marker = '.'; + % h = line(V_centers(:,1),V_centers(:,2),V_centers(:,3)); + % h.LineStyle = 'none'; + % h.Marker = '.'; + % h.Color = Color.red; + % axis equal + + + I_0 = nodeIndecies(1:end-1, 1:end-1); + I_1 = nodeIndecies(2:end, 1:end-1); + I_2 = nodeIndecies(2:end, 2:end); + I_3 = nodeIndecies(1:end-1, 2:end); + I_C = centerIndecies; + + S.Vertices = [ + V_nodes; + V_centers; + ]; + + S.Faces = [ + I_0(:), I_1(:), I_C(:); + I_1(:), I_2(:), I_C(:); + I_2(:), I_3(:), I_C(:); + I_3(:), I_0(:), I_C(:); + ]; + + % figure() + h = patch(S, 'FaceVertexCData', S.Vertices(:,3),'FaceColor','flat'); +end + +% Calculate the mean of four neighbours around a patch +function M = neighbourMean(A) + M = (A(1:end-1, 1:end-1) + A(2:end, 1:end-1) + A(1:end-1, 2:end) + A(2:end, 2:end))/4; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/symmetricPart.m Mon Sep 11 14:17:15 2017 +0200 @@ -0,0 +1,4 @@ +% Returns the symmetric of A +function S = symmetricPart(A, tol) + S = 1/2*(A + A'); +end
--- a/textTable.m Mon Sep 11 14:12:54 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,64 +0,0 @@ - -% data -- cell array of numbers -% leftColstrings -- cell array of strings, for left column -% topRowStrings -- cell array of strings, for top row -% dataFormat -- (optional) format specifier, e.g. '%.2f' -function textTable(data, leftColStrings, topRowStrings, dataFormat) - - default_arg('dataFormat','%.2f') - - nRows = length(leftColStrings); - nCols = length(topRowStrings); - [m,n] = size(data); - - if(m ~= nRows || n ~=nCols) - error('Data dimensions must match labels'); - end - - % Find column widths - headerLength = 0; - for i = 1:nCols - headerLength = max(headerLength, length(topRowStrings{i} )); - end - - dataLength = 0; - for i = 1:nRows - for j = 1:nCols - temp = length(sprintf(dataFormat, data{i,j})); - dataLength = max(dataLength, temp); - end - end - dataLength = length(sprintf(dataFormat, data{1,1})); - - colWidth = max(headerLength,dataLength); - - % Print headers - fprintf(' %*s |',colWidth,'') - for i = 1:nCols - fprintf(' %-*s |', colWidth, topRowStrings{i}); - end - fprintf('\n'); - - % Print divider - m_dev = repmat('-',1,colWidth); - column_dev = repmat('-',1,colWidth); - fprintf('-%s-+',m_dev); - for i = 1:nCols - fprintf('-%s-+', column_dev); - end - fprintf('\n'); - - - % Print each row - dataFormat = ['%*' dataFormat(2:end)]; - for i = 1:nRows - fprintf(' %*s |',colWidth,leftColStrings{i}); - for j = 1:nCols - fprintf([' ' dataFormat ' |'], colWidth, data{i,j}); - end - fprintf('\n'); - end - - fprintf('\n'); - -end \ No newline at end of file
--- a/toString.m Mon Sep 11 14:12:54 2017 +0200 +++ b/toString.m Mon Sep 11 14:17:15 2017 +0200 @@ -51,6 +51,14 @@ end function str = struct2string(s) + if isscalar(s) + str = structScalar2string(s); + else + str = structArray2string(s); + end +end + +function str = structScalar2string(s) fn = fieldnames(s); if length(fn) == 0 @@ -68,3 +76,32 @@ str = [str sprintf('%s: %s}',fn{end}, value2string(value))]; end +function str = structArray2string(s) + fn = fieldnames(s); + + if length(fn) == 0 + str = '{}'; + return + end + + stringArray = cell(length(s)+1, length(fn)+1); + + stringArray(1,2:end) = fn; + + for i = 1:length(s) + stringArray{i+1,1} = i; + for j = 1:length(fn) + valueStr = value2string(s(i).(fn{j})); + stringArray{i+1,j+1} = valueStr; + end + end + + tt = TextTable(stringArray); + tt.fmtArray(2:end, 1) = {'%d'}; + tt.vertDiv = [1]; + tt.horzDiv = [1]; + str = tt.toString(); +end + + +