Mercurial > repos > public > sbplib
changeset 200:ef41fde95ac4 feature/beams
Merged feature/grids into feature/beams.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 13 Jun 2016 16:59:02 +0200 |
| parents | 419ec303e97d (current diff) d18096820ed4 (diff) |
| children | 3c4ffbfbfb84 |
| files | +grid/Multiblock.m +multiblock/gridVector1d.m +multiblock/gridVector2d.m +multiblock/solutionVector2cell.m +scheme/Scheme.m |
| diffstat | 31 files changed, 938 insertions(+), 96 deletions(-) [+] |
line wrap: on
line diff
diff -r 419ec303e97d -r ef41fde95ac4 +anim/animate.m --- a/+anim/animate.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+anim/animate.m Mon Jun 13 16:59:02 2016 +0200 @@ -19,6 +19,7 @@ dTau_target = 1/target_frame_rate; % Real time between frames rs = util.ReplaceableString(); + rs.appendFormat(' t: %d\n'); rs.appendFormat(' tau: %d\n'); rs.appendFormat(' target tau: %d\n'); rs.appendFormat(' Target fps: %.2f\n'); @@ -59,7 +60,7 @@ % Update information about this frame tau = toc(animation_start); - rs.updateParam(tau, targetTau, 1/dTau_target, 1/dTau, time_modifier_bound, time_modifier); + rs.updateParam(t, tau, targetTau, 1/dTau_target, 1/dTau, time_modifier_bound, time_modifier); end
diff -r 419ec303e97d -r ef41fde95ac4 +anim/setup_1d_plot.m --- a/+anim/setup_1d_plot.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+anim/setup_1d_plot.m Mon Jun 13 16:59:02 2016 +0200 @@ -9,8 +9,9 @@ % be passed to functions in yfun. % plot_handles - Array of plot_handles. One for each yfun. % axis_handle - Handle to the axis plotted to. -function [update_data, plot_handles, axis_handle] = setup_1d_plot(x,yfun) +function [update_data, plot_handles, axis_handle] = setup_1d_plot(x,yfun,show_title) default_arg('yfun',{@(y)y}); + default_arg('show_title', true) if isa(yfun,'function_handle') yfun = {yfun}; @@ -24,8 +25,6 @@ axis_handle = gca; - xlabel('x') - ylabel('y') xlim([x(1) x(end)]); function update(t,varargin) @@ -33,7 +32,10 @@ for i = 1:length(yfun) set(plot_handles(i),'YData',yfun{i}(varargin{:})); end - title(axis_handle,sprintf('T=%.3f',t)); + + if show_title + title(axis_handle,sprintf('T=%.3f',t)); + end end end update_data = @update;
diff -r 419ec303e97d -r ef41fde95ac4 +grid/Cartesian.m --- a/+grid/Cartesian.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/Cartesian.m Mon Jun 13 16:59:02 2016 +0200 @@ -122,5 +122,68 @@ function gf = projectFunc(obj, gf, g) error('grid:Cartesian:NotImplemented') end + + % Return the names of all boundaries in this grid. + function bs = getBoundaryNames(obj) + switch obj.D() + case 1 + bs = {'l', 'r'}; + case 2 + bs = {'w', 'e', 's', 'n'}; + case 3 + bs = {'w', 'e', 's', 'n', 'd', 'u'}; + otherwise + error('not implemented'); + end + end + + % Return coordinates for the given boundary + function X = getBoundary(obj, name) + % In what dimension is the boundary? + switch name + case {'l', 'r', 'w', 'e'} + D = 1; + case {'s', 'n'} + D = 2; + case {'d', 'u'} + D = 3; + otherwise + error('not implemented'); + end + + % At what index is the boundary? + switch name + case {'l', 'w', 's', 'd'} + index = 1; + case {'r', 'e', 'n', 'u'} + index = obj.m(D); + otherwise + error('not implemented'); + end + + + + I = cell(1, obj.d); + for i = 1:obj.d + if i == D + I{i} = index; + else + I{i} = ':'; + end + end + + % Calculate size of result: + m = obj.m; + m(D) = []; + N = prod(m); + + X = zeros(N, obj.d); + + coordMat = obj.matrices(); + for i = 1:length(coordMat) + Xtemp = coordMat{i}(I{:}); + X(:,i) = reshapeRowMaj(Xtemp, [N,1]); + end + end end end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +grid/CartesianTest.m --- a/+grid/CartesianTest.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/CartesianTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -194,3 +194,125 @@ testCase.verifyEqual(g.scaling(),[2 1]); end + + +function testGetBoundaryNames(testCase) + in = { + {[1 2 3]}, + {[1 2 3], [4 5]}, + {[1 2 3], [4 5], [6 7 8]}, + }; + + out = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.getBoundaryNames(), out{i}); + end +end + +function testGetBoundary(testCase) + grids = { + {[1 2 3]}, + {[1 2 3], [4 5]}, + {[1 2 3], [4 5], [6 7 8]}, + }; + + boundaries = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + + % 1d + out{1,1} = 1; + out{1,2} = 3; + + % 2d + out{2,1} = [ + 1,4; + 1,5; + ]; + out{2,2} = [ + 3,4; + 3,5; + ]; + out{2,3} = [ + 1,4; + 2,4; + 3,4; + ]; + out{2,4} = [ + 1,5; + 2,5; + 3,5; + ]; + + % 3d + out{3,1} = [ + 1,4,6; + 1,4,7; + 1,4,8; + 1,5,6; + 1,5,7; + 1,5,8; + ]; + out{3,2} = [ + 3,4,6; + 3,4,7; + 3,4,8; + 3,5,6; + 3,5,7; + 3,5,8; + ]; + out{3,3} = [ + 1,4,6; + 1,4,7; + 1,4,8; + 2,4,6; + 2,4,7; + 2,4,8; + 3,4,6; + 3,4,7; + 3,4,8; + ]; + out{3,4} = [ + 1,5,6; + 1,5,7; + 1,5,8; + 2,5,6; + 2,5,7; + 2,5,8; + 3,5,6; + 3,5,7; + 3,5,8; + ]; + out{3,5} = [ + 1,4,6; + 1,5,6; + 2,4,6; + 2,5,6; + 3,4,6; + 3,5,6; + ]; + out{3,6} = [ + 1,4,8; + 1,5,8; + 2,4,8; + 2,5,8; + 3,4,8; + 3,5,8; + ]; + + for ig = 1:length(grids) + g = grid.Cartesian(grids{ig}{:}); + for ib = 1:length(boundaries{ig}) + testCase.verifyEqual(g.getBoundary(boundaries{ig}{ib}), out{ig,ib}); + end + end +end
diff -r 419ec303e97d -r ef41fde95ac4 +grid/Curvilinear.m --- a/+grid/Curvilinear.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/Curvilinear.m Mon Jun 13 16:59:02 2016 +0200 @@ -80,6 +80,64 @@ end h = obj.logic.h; end + + % Return the names of all boundaries in this grid. + function bs = getBoundaryNames(obj) + bs = obj.logic.getBoundaryNames(); + end + + % Return coordinates for the given boundary + function X = getBoundary(obj, name) + % In what dimension is the boundary? + switch name + case {'l', 'r', 'w', 'e'} + D = 1; + case {'s', 'n'} + D = 2; + case {'d', 'u'} + D = 3; + otherwise + error('not implemented'); + end + + % At what index is the boundary? + switch name + case {'l', 'w', 's', 'd'} + index = 1; + case {'r', 'e', 'n', 'u'} + index = obj.logic.m(D); + otherwise + error('not implemented'); + end + + + + I = cell(1, obj.D); + for i = 1:obj.D + if i == D + I{i} = index; + else + I{i} = ':'; + end + end + + % Calculate size of result: + m = obj.logic.m; + m(D) = []; + N = prod(m); + + X = zeros(N, obj.D); + + p = obj.points; + for i = 1:obj.D() + coordMat{i} = reshapeRowMaj(p(:,i), obj.logic.m); + end + + for i = 1:length(coordMat) + Xtemp = coordMat{i}(I{:}); + X(:,i) = reshapeRowMaj(Xtemp, [N,1]); + end + end end end
diff -r 419ec303e97d -r ef41fde95ac4 +grid/CurvilinearTest.m --- a/+grid/CurvilinearTest.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/CurvilinearTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -67,7 +67,7 @@ end function testMappingInputError(testCase) - testCase.assumeFail(); + testCase.verifyFail(); end function testScaling(testCase) @@ -80,3 +80,52 @@ testCase.verifyEqual(g.scaling(),[2 1]); end +function testGetBoundaryNames(testCase) + in = { + {{1:10}, @(x) exp(x)}, + {{1:10,1:6}, @(x,y) [exp(x+y); exp(x-y)]}, + {{1:10,1:5,1:7}, @(x,y,z)[exp(x+y+z); exp(x-y-z); 2+x+y-z]}, + }; + + out = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + for i = 1:length(in) + g = grid.Curvilinear(in{i}{2},in{i}{1}{:}); + testCase.verifyEqual(g.getBoundaryNames(), out{i}); + end +end + +function testGetBoundary(testCase) + grids = { + {{1:10}, @(x) exp(x)}, + {{1:10,1:6}, @(x,y) [exp(x+y); exp(x-y)]}, + {{1:10,1:5,1:7}, @(x,y,z)[exp(x+y+z); exp(x-y-z); 2+x+y-z]}, + }; + + boundaries = { + {'l', 'r'}, + {'w', 'e', 's', 'n'}, + {'w', 'e', 's', 'n', 'd', 'u'}, + }; + + + for ig = 1:length(grids) + g = grid.Curvilinear(grids{ig}{2},grids{ig}{1}{:}); + + logicalGrid = grid.Cartesian(grids{ig}{1}{:}); + + for ib = 1:length(boundaries{ig}) + + logicalBoundary = logicalGrid.getBoundary(boundaries{ig}{ib}); + + x = num2cell(logicalBoundary',2); + expectedBoundary = grids{ig}{2}(x{:})'; + testCase.verifyEqual(g.getBoundary(boundaries{ig}{ib}), expectedBoundary); + end + end +end +
diff -r 419ec303e97d -r ef41fde95ac4 +grid/Grid.m --- a/+grid/Grid.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/Grid.m Mon Jun 13 16:59:02 2016 +0200 @@ -15,13 +15,11 @@ % 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. + bs = getBoundaryNames(obj) + + % Return coordinates for the given boundary + b = getBoundary(obj, name) end end - - -%% Should it be able to return a cell size aswell? For an equidistant grid this would be know -%% for other grids the constructor would have to make something up. -%% For example the grid.Cartesian constructor would take a h (1 x d) vector as an in parameter. - - -%Should define boundaries somehow, look in stitchSchemes. \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +grid/Multiblock.m --- a/+grid/Multiblock.m Mon Feb 29 15:40:34 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -classdef Multiblock < grid.Grid - % General multiblock grid - methods (Abstract) - % NBlocks returns the number of blocks in the grid. - o = NBlocks(obj); - - % Grid returns the ith grid in the multiblockgrid - gs = grid(obj,i); - - % Grids returns a cell array of all the grids in the multiblock grid. - gs = grids(obj); - - % Split a grid function on obj to a cell array of grid function on each block - gf = splitFunc(gf) - end -end - - -% Should define boundaries and connections between grids. \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +grid/bspline.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/bspline.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,47 @@ +% Calculates a D dimensional p-order bspline at t with knots T and control points P. +% T = [t0 t1 t2 ... tm] is a 1 x (m+1) vector with non-decresing elements and t0 = 0 tm = 1. +% P = [P0 P1 P2 ... Pn] is a D x (n+1) matrix. + +% knots p+1 to m-p-1 are the internal knots + +% Implemented from: http://mathworld.wolfram.com/B-Spline.html +function C = bspline(t,p,P,T) + m = length(T) - 1; + n = size(P,2) - 1; + D = size(P,1); + + assert(p == m - n - 1); + + C = zeros(D,length(t)); + + for i = 0:n + for k = 1:D + C(k,:) = C(k,:) + P(k,1+i)*B(i,p,t,T); + end + end + + % Curve not defined for t = 1 ? Ugly fix: + I = find(t == 1); + C(:,I) = repmat(P(:,end),[1,length(I)]); +end + +function o = B(i, j, t, T) + if j == 0 + o = T(1+i) <= t & t < T(1+i+1); + return + end + + if T(1+i+j)-T(1+i) ~= 0 + a = (t-T(1+i))/(T(1+i+j)-T(1+i)); + else + a = t*0; + end + + if T(1+i+j+1)-T(1+i+1) ~= 0 + b = (T(1+i+j+1)-t)/(T(1+i+j+1)-T(1+i+1)); + else + b = t*0; + end + + o = a.*B(i, j-1, t, T) + b.*B(i+1, j-1, t, T); +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +grid/equidistantCurvilinearTest.m --- a/+grid/equidistantCurvilinearTest.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+grid/equidistantCurvilinearTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -3,5 +3,5 @@ end function testNoTests(testCase) - testCase.assumeFail(); + testCase.verifyFail(); end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/BoundaryGroup.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/BoundaryGroup.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,51 @@ +% BoundaryGroup defines a boundary grouping in a multiblock grid. +classdef BoundaryGroup + properties + blockIDs + names + end + + 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 + end + + function display(obj, name) + + disp(' ') + 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') + end + end +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/BoundaryGroupTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/BoundaryGroupTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,30 @@ +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
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/DiffOp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOp.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,121 @@ +classdef DiffOp < scheme.Scheme + properties + grid + order + diffOps + D + H + end + + methods + function obj = DiffOp(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); + + 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); + obj.diffOps{i} = h(grid.grid{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 = cell2sparse(H); + + + % Build the differentiation matrix + D = cell(nBlocks, nBlocks); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; + end + + for i = 1:nBlocks + for j = i:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + [ii, ij] = obj.diffOps{i}.inteface_coupling(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}.inteface_coupling(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 = cell2sparse(D); + + % end + + 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){}; + elseif iscell(doParam) && length(doParam) == grid.nBlocks() + getParam = @(i)doParam{i}; + else + getParam = @(i)doParam; + end + 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 m = boundary_condition(obj,boundary,type,data) + + end + + function m = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + end + + + function N = size(obj) + N = 0; + for i = 1:length(diffOps) + N = N + diffOps.size(); + end + end + end +end
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/DiffOpTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,20 @@ +function tests = DiffOpTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + g = multiblock.Grid({},{}); + doHand = @(grid,order)[]; + order = 0; + do = multiblock.DiffOp(doHand, g, order); +end + +function testMissing(testCase) + testCase.verifyFail(); +end + + +% function do = mockDiffOp() +% do.H = 1; +% do.D = 1; +% end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/Grid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Grid.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,111 @@ +classdef Grid < grid.Grid + properties + grids + connections + boundaryGroups + + nPoints + end + + % 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. + function obj = Grid(grids, connections, boundaryGroups) + obj.grids = grids; + obj.connections = connections; + + obj.nPoints = 0; + for i = 1:length(grids) + obj.nPoints = obj.nPoints + grids{i}.N(); + end + + % if iscell(boundaryGroups) + end + + function n = size(obj) + n = length(obj.grids); + end + + % n returns the number of points in the grid + function o = N(obj) + o = obj.nPoints; + end + + function n = nBlocks(obj) + n = length(obj.grids); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.grids{1}.D(); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = []; + for i = 1:length(obj.grids) + X = [X; obj.grids{i}.points]; + end + end + + % Split a grid function on obj to a cell array of grid function on each block + function gfs = splitFunc(obj, gf) + nComponents = length(gf)/obj.nPoints; + nBlocks = length(obj.grids); + + % Collect number of points in each block + N = cell(1,nBlocks); + for i = 1:nBlocks + N{i} = obj.grids{i}.N(); + end + + gfs = mat2cell(gf, N, 1); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gfs = obj.splitFunc(gf); + + for i = 1:length(obj.grids) + gfs{i} = obj.grids{i}.restrictFunc(gfs{i}, g.grids{i}); + end + + gf = cell2mat(gfs); + end + + % Projects the grid function gf on obj to the grid g. + function o = projectFunc(obj, gf, g) + error('not implemented') + + p = g.points(); + o = zeros(length(p),1); + for i = 1:length(p) + I = whatGrid(p(i)); + o(i) = obj.grids{I}.projectFunc(gf, p(i)); + end + + + function I = whatGrid(p) + % Find what grid a point lies on + end + + end + + function bs = getBoundaryNames(obj) + bs = []; + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = []; + end + end +end
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/GridTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/GridTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,37 @@ +function tests = GridTest() + tests = functiontests(localfunctions); +end + +function testCreation(testCase) + g = multiblock.Grid({},{}); +end + +function testMissing(testCase) + testCase.verifyFail(); +end + +function testGetBoundaryNames(testCase) + [grids, conn] = prepareAdjecentBlocks(); + + mbg = multiblock.Grid(grids, conn, multiblock.BoundaryGroup({1,'w'},{2,'w'}) ); + + testCase.verifyFail(); +end + +function testGetBoundary(testCase) + [grids, conn] = prepareAdjecentBlocks(); + + mbg = multiblock.Grid(grids, conn, multiblock.BoundaryGroup({1,'w'},{2,'w'}) ); + testCase.verifyFail(); +end + + +function [grids, conn] = prepareAdjecentBlocks() + grids = { + grid.Cartesian([0 1 2], [3 4 5]); + grid.Cartesian([1 2], [10 20]); + }; + + conn = cell(2,2); + conn{1, 2} = {'s','n'}; +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/gridVector1d.m --- a/+multiblock/gridVector1d.m Mon Feb 29 15:40:34 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function x = gridVector1d(schms) - n = length(schms); - x = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - end -end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/gridVector2d.m --- a/+multiblock/gridVector2d.m Mon Feb 29 15:40:34 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -function [x,y] = gridVector2d(schms) - n = length(schms); - x = cell(n,1); - y = cell(n,1); - - for i = 1:n - x{i} = schms{i}.x; - y{i} = schms{i}.y; - end -end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +multiblock/solutionVector2cell.m --- a/+multiblock/solutionVector2cell.m Mon Feb 29 15:40:34 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,15 +0,0 @@ -function v_cell = solutionVector2cell(schms,v,components) - if nargin == 2 - components = 1; - end - - N = length(schms); - v_cell = cell(N,1); - - i_start = 1; - for i =1:N - i_end = i_start + schms{i}.size()*components - 1; - v_cell{i} = v(i_start:i_end); - i_start = i_end+1; - end -end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 +scheme/Scheme.m --- a/+scheme/Scheme.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+scheme/Scheme.m Mon Jun 13 16:59:02 2016 +0200 @@ -1,40 +1,45 @@ -% Start with all matrix returns. When that works see how we should generalize to non-matrix stuff/nonlinear +% Start with all matrix returns. When that works see how we should generalize +% to non-matrix stuff/nonlinear classdef Scheme < handle properties (Abstract) order % Order accuracy for the approximation - % vectors u,v,w depending on dim that gives were gridpoints are in each dimension - % vectors x,y,z containing the x,y,z values corresponding to each grid point - % matrices X,Y,Z with point coordinates as multi dimensional vectors + grid D % non-stabalized scheme operator H % Discrete norm - - % Should also containg: - % the grid points used - % the grid spacing end methods (Abstract) - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % neighbour_scheme is an instance of Scheme that should be interfaced to. - % neighbour_boundary is a string specifying which boundary to interface to. - [closure, penalty] = boundary_condition(obj,boundary,type) + % Closure functions return the opertors applied to the own doamin to + % close the boundary Penalty functions return the opertors to force + % the solution. In the case of an interface it returns the operator + % applied to the other doamin. In some cases the penalty return value + % can be ommited and the closure function take care of both parts. + % boundary is a string specifying the boundary e.g. + % 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of + % boundary condition if there are several. + % data is a function returning the data that + % should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be + % interfaced to. + % neighbour_boundary is a string specifying which boundary to + % interface to. + [closure, penalty] = boundary_condition(obj,boundary,type,data) [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - N = size(obj) % Returns the number of degrees of freedom. + % Returns the number of degrees of freedom. + N = size(obj) end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % Calculates the matrcis need for the inteface coupling between + % boundary bound_u of scheme schm_u and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); end end -end \ No newline at end of file +end
diff -r 419ec303e97d -r ef41fde95ac4 +util/calc_borrowing.m --- a/+util/calc_borrowing.m Mon Feb 29 15:40:34 2016 +0100 +++ b/+util/calc_borrowing.m Mon Jun 13 16:59:02 2016 +0200 @@ -32,7 +32,7 @@ %% 2nd order compatible -[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher2_compatible(m,h); +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible2(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; @@ -46,7 +46,7 @@ %% 4th order compatible -[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4_compatible(m,h); +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible4(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; @@ -59,7 +59,7 @@ fprintf('\n') %% 6th order compatible -[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6_compatible(m,h); +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible6(m,h); S1 = S_1*S_1' + S_m*S_m'; S2 = S2_1*S2_1' + S2_m*S2_m'; S3 = S3_1*S3_1' + S3_m*S3_m'; @@ -69,3 +69,27 @@ fprintf('6th order compatible\n') fprintf('alpha_II: %.10f\n',alpha_II) fprintf('alpha_III: %.10f\n',alpha_III) +fprintf('\n')3 + + + + + +% Ordinary + +for order = [2 4 6 8 10] + op = sbp.Ordinary(m,h, order); + + S_1 = op.boundary.S_1; + S_m = op.boundary.S_m; + + M = op.norms.M; + + S1 = S_1*S_1' + S_m*S_m'; + alpha = util.matrixborrow(M, h*S1); + fprintf('%dth order Ordinary\n', order) + fprintf('alpha: %.10f\n', alpha) + fprintf('\n') +end + +
diff -r 419ec303e97d -r ef41fde95ac4 getVarname.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/getVarname.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,7 @@ +function names = getVarname(varargin) + names = cell(size(varargin)); + + for i = 1:numel(varargin) + names{i} = inputname(i); + end +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 labelTask.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/labelTask.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,12 @@ +function o = labelTask(in) + + switch class(in) + case 'char' + fprintf(in); + o = tic(); + case 'uint64' + o = toc(in); + fprintf(' - done %fs\n', o); + otherwise + error('Unknow input type: %s', class(in)) + end
diff -r 419ec303e97d -r ef41fde95ac4 prof.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prof.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,10 @@ +function prof(f) + profile on + try + f(); + profile viewer + catch e + fprintf(2, '\n%s', getReport(e)); + profile clear + end +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 saveFigurePosition.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/saveFigurePosition.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,6 @@ +function saveFigurePosition() + defaultPosition = get(0,'defaultfigureposition'); + f = gcf; + defaultPosition(1:2) = f.Position(1:2); + set(0,'defaultfigureposition',defaultPosition); +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 savepng.m --- a/savepng.m Mon Feb 29 15:40:34 2016 +0100 +++ b/savepng.m Mon Jun 13 16:59:02 2016 +0200 @@ -1,4 +1,5 @@ function savepng(h, filename, dpi) + default_arg('dpi', 300) print(h,filename,'-dpng',sprintf('-r%d',dpi)); % Let print size in inches as input parameter % Smaller boundingbox
diff -r 419ec303e97d -r ef41fde95ac4 setFontSize.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/setFontSize.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,7 @@ +% setFontSize(fig, pts) +% Sets all fontsizes within a figure to size pts +% The default value for pts is 16. +function setFontSize(fig, pts) + default_arg('pts', 16); + set(findall(fig,'-property','FontSize'),'FontSize',pts); +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 sparse2cell.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sparse2cell.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,28 @@ +% sparse2cell breaks a sparse matrix up into a cell matrix of sparse matrices. +% any zero submatrix creates a empty cell in the cell matrix. +% A -- NxM sparse matrix +% d1, d2 -- vectors of sub matrix sizes for each dimensions. Must have sum(di) == Ni. +% Example: +% C = sparse2cell(A,[5 10], [10 5]) +function C = sparse2cell(A, d1, d2) + [n, m] = size(A); + if n ~= sum(d1) || m ~= sum(d2) + error('sparse2cell:NonMatchingDim','The elements of d1 and d2 must sum to N and M.'); + end + + C = cell(length(d1), length(d2)); + I = 1; + for i = 1:length(d1) + J = 1; + for j = 1:length(d2) + Asub = A(I:(I + d1(i)-1), J:(J + d2(j)-1)); + if nnz(Asub) == 0 + C{i,j} = []; + else + C{i,j} = Asub; + end + J = J + d2(j); + end + I = I + d1(i); + end +end
diff -r 419ec303e97d -r ef41fde95ac4 sparse2cellTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sparse2cellTest.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,43 @@ +function tests = sparse2cellTest() + tests = functiontests(localfunctions); +end + +function testErrorNonMatchingDim(testCase) + in = { + {magic(5), [1 2 3], [4]}, + {magic(5), [1 1 1 1 1 1], [5]}, + {magic(5), [5], [1 1 1 1 1 1]}, + {ones(4,2),[2 3],[2]}, + {ones(4,2),[2 2],[3]}, + }; + + for i = 1:length(in) + testCase.verifyError(@()sparse2cell(in{i}{:}),'sparse2cell:NonMatchingDim'); + end +end + +function testOutput(testCase) + in = {}; + out = {}; + in{1}{1} =[17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; 11 18 25 2 9]; + in{1}{2} = [1 4]; + in{1}{3} = [2 3]; + + out{1} = { + [17 24], [1 8 15]; + [23 5; 4 6; 10 12; 11 18], [7 14 16; 13 20 22; 19 21 3; 25 2 9]; + }; + + in{1}{1} = [17 24 1 8 15; 23 5 0 0 0; 4 6 0 0 0; 10 12 0 0 0; 11 18 0 0 0]; + in{1}{2} = [1 4]; + in{1}{3} = [2 3]; + + out{1} = { + [17 24], [1 8 15]; + [23 5; 4 6; 10 12; 11 18], []; + }; + + for i = 1:length(in) + testCase.verifyEqual(sparse2cell(in{i}{:}), out{i}); + end +end \ No newline at end of file
diff -r 419ec303e97d -r ef41fde95ac4 time.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/time.m Mon Jun 13 16:59:02 2016 +0200 @@ -0,0 +1,34 @@ +function t_out = time(f, n) + default_arg('n',1); + + if n == 1 + t = timeOnce(f); + else + t = timeAvg(f, n); + end + + if nargout == 1 + t_out = t; + else + fprintf('Elapsed time is %.6f seconds.\n', t) + end +end + +function t = timeOnce(f) + s = tic(); + + f(); + + t = toc(s); +end + + +function t = timeAvg(f, n) + s = tic(); + + for i = 1:n + f(); + end + + t = toc(s)/n; +end
diff -r 419ec303e97d -r ef41fde95ac4 toString.m --- a/toString.m Mon Feb 29 15:40:34 2016 +0100 +++ b/toString.m Mon Jun 13 16:59:02 2016 +0200 @@ -28,19 +28,26 @@ end function str = cell2string(c) - len = length(c); - - if len == 0 + if isempty(c) str = '{}'; return end + [n, m] = size(c); + str = '{'; - for i =1:len-1 - str = [str sprintf('%s, ', value2string(c{i}))]; + for i = 1:n-1 + for j = 1:m-1 + str = [str sprintf('%s, ', value2string(c{i,j}))]; + end + str = [str sprintf('%s; ', value2string(c{i,end}))]; end - str = [str sprintf('%s}', value2string(c{end}))]; + + for j = 1:m-1 + str = [str sprintf('%s, ', value2string(c{end,j}))]; + end + str = [str sprintf('%s}', value2string(c{end,end}))]; end function str = struct2string(s)