Mercurial > repos > public > sbplib
changeset 781:69ab0e69f972 feature/interpolation
Merge with feature/grids
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 24 Jul 2018 20:14:29 -0700 |
| parents | 005a8d071da3 (current diff) e7a6744499fa (diff) |
| children | 3c3280ebabb3 |
| files | +multiblock/Def.m diffSymfun.m |
| diffstat | 44 files changed, 1208 insertions(+), 256 deletions(-) [+] |
line wrap: on
line diff
diff -r 005a8d071da3 -r 69ab0e69f972 +anim/setup_time_quantity_plot.m --- a/+anim/setup_time_quantity_plot.m Tue May 22 13:29:47 2018 -0700 +++ b/+anim/setup_time_quantity_plot.m Tue Jul 24 20:14:29 2018 -0700 @@ -16,7 +16,7 @@ if ishandle(axis_handle) % t = [t t_now]; for j = 1:length(yfun) - addpoints(plot_handles(j),t_now,yfun{j}(varargin{:})); + addpoints(plot_handles(j),t_now,full(yfun{j}(varargin{:}))); end [t,~] = getpoints(plot_handles(1));
diff -r 005a8d071da3 -r 69ab0e69f972 +blockmatrix/toMatrix.m --- a/+blockmatrix/toMatrix.m Tue May 22 13:29:47 2018 -0700 +++ b/+blockmatrix/toMatrix.m Tue Jul 24 20:14:29 2018 -0700 @@ -12,12 +12,9 @@ A = sparse(N,M); - n_ind = [0 cumsum(n)]; - m_ind = [0 cumsum(m)]; - for i = 1:size(bm,1) for j = 1:size(bm,2) - if(isempty(bm{i,j})) + if isempty(bm{i,j}) bm{i,j} = sparse(n(i),m(j)); end end
diff -r 005a8d071da3 -r 69ab0e69f972 +draw/prompt_point.m --- a/+draw/prompt_point.m Tue May 22 13:29:47 2018 -0700 +++ b/+draw/prompt_point.m Tue Jul 24 20:14:29 2018 -0700 @@ -1,22 +1,23 @@ -function [p, button] = prompt_point(s,varargin) +function [p, button] = prompt_point(s, varargin) default_arg('s',[]) set(gcf,'Pointer','crosshair') if ~isempty(s) - fprintf(s,varargin{:}); + fprintf(s, varargin{:}); end - a = gca; + fh = gcf(); + ah = gca(); - function get_point(src,event) - cp = a.CurrentPoint; + function get_point(src, event) + cp = ah.CurrentPoint; p = cp(1,1:2)'; - a.ButtonDownFcn = []; + fh.WindowButtonUpFcn = []; end - a.ButtonDownFcn = @get_point; - waitfor(a,'ButtonDownFcn', []) + fh.WindowButtonUpFcn = @get_point; + waitfor(fh,'WindowButtonUpFcn', []) set(gcf,'Pointer','arrow')
diff -r 005a8d071da3 -r 69ab0e69f972 +grid/Cartesian.m --- a/+grid/Cartesian.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/Cartesian.m Tue Jul 24 20:14:29 2018 -0700 @@ -15,7 +15,7 @@ obj.d = length(varargin); for i = 1:obj.d - assert(isvector(varargin{i}), 'Coordinate inputs must be a vectors.') + assert(isnumeric(varargin{i}), 'Coordinate inputs must be vectors.') obj.x{i} = varargin{i}; obj.m(i) = length(varargin{i});
diff -r 005a8d071da3 -r 69ab0e69f972 +grid/TODO.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/TODO.txt Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,1 @@ +% TODO: Rename grid package. name conflicts with built in function
diff -r 005a8d071da3 -r 69ab0e69f972 +grid/evalOn.m --- a/+grid/evalOn.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/evalOn.m Tue Jul 24 20:14:29 2018 -0700 @@ -7,68 +7,34 @@ function gf = evalOn(g, func) if ~isa(func, 'function_handle') % We should have a constant. - if size(func,2) ~= 1 - error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') - end + assert(size(func,2) == 1,'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector'); gf = repmat(func,[g.N, 1]); return end % func should now be a function_handle + assert(g.D == nargin(func),'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') - if g.D ~= nargin(func) - g.D - nargin(func) - error('grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') - end + x = num2cell(g.points(),1); + k = numberOfComponents(func); + gf = func(x{:}); + gf = reorderComponents(gf, k); +end - % Get coordinates - x = g.points(); +% Find the number of vector components of func +function k = numberOfComponents(func) + x0 = num2cell(ones(1,nargin(func))); + f0 = func(x0{:}); + assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector'); + k = length(f0); +end - % Find the number of components - if size(x,1) ~= 0 - x0 = x(1,:); - else - x0 = num2cell(ones(1,size(x,2))); +% Reorder the components of the function to sit together +function gf = reorderComponents(a, k) + N = length(a)/k; + gf = zeros(N*k, 1); + for i = 1:k + gf(i:k:end) = a((i-1)*N + 1 : i*N); end - - dim = length(x0); - % Evaluate f0 = func(x0(1),x0(2),...,x0(dim)); - if(dim == 1) - f0 = func(x0); - else - eval_str = 'f0 = func(x0(1)'; - for i = 2:dim - eval_str = [eval_str, sprintf(',x0(%d)',i)]; - end - eval_str = [eval_str, ');']; - eval(eval_str); - end - - % k = number of components - k = length(f0); - - if size(f0,2) ~= 1 - error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') - end - - % Evaluate gf = func(x(:,1),x(:,2),...,x(:,dim)); - if(dim == 1) - gf = func(x); - else - eval_str = 'gf = func(x(:,1)'; - for i = 2:dim - eval_str = [eval_str, sprintf(',x(:,%d)',i)]; - end - eval_str = [eval_str, ');']; - eval(eval_str); - end - - % Reorganize gf - gf_temp = gf; - gf = zeros(g.N*k, 1); - for i = 1:k - gf(i:k:end) = gf_temp((i-1)*g.N + 1 : i*g.N); - end -end \ No newline at end of file +end
diff -r 005a8d071da3 -r 69ab0e69f972 +grid/evalOnTest.m --- a/+grid/evalOnTest.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/evalOnTest.m Tue Jul 24 20:14:29 2018 -0700 @@ -31,7 +31,7 @@ cases = { {getTestGrid('1d'), @(x,y)x-y}, {getTestGrid('2d'), @(x)x }, - } + }; for i = 1:length(cases) g = cases{i}{1}; @@ -111,9 +111,9 @@ function testInputErrorVectorValued(testCase) - in = { + in = { [1,2,3], - @(x,y)[x,-y]; + @(x,y)[x,-y], }; g = getTestGrid('2d');
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/+domain/Circle.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Circle.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,98 @@ +classdef Circle < multiblock.DefCurvilinear + properties + r, c + + hs + r_arc + omega + end + + methods + function obj = Circle(r, c, hs) + default_arg('r', 1); + default_arg('c', [0; 0]); + default_arg('hs', 0.435); + + + % alpha = 0.75; + % hs = alpha*r/sqrt(2); + + % Square should not be a square, it should be an arc. The arc radius + % is chosen so that the three angles of the meshes are all equal. + % This gives that the (half)arc opening angle of should be omega = pi/12 + omega = pi/12; + r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); % = hs* 1/sin(omega) + c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0]; + + cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]); + + c2 = cir(0); + c3 = cir(1); + + s1 = [-hs; -hs]; + s2 = [ hs; -hs]; + s3 = [ hs; hs]; + s4 = [-hs; hs]; + + sp2 = parametrization.Curve.line(s2,c2); + sp3 = parametrization.Curve.line(s3,c3); + + Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]); + Se2 = Se1.rotate(c,pi/2); + Se3 = Se2.rotate(c,pi/2); + Se4 = Se3.rotate(c,pi/2); + + + S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1); + + A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse); + B = A.rotate(c,1*pi/2).rotate_edges(-1); + C = A.rotate(c,2*pi/2).rotate_edges(-1); + D = A.rotate(c,3*pi/2).rotate_edges(0); + + blocks = {S,A,B,C,D}; + blocksNames = {'S','A','B','C','D'}; + + conn = cell(5,5); + conn{1,2} = {'e','w'}; + conn{1,3} = {'n','s'}; + conn{1,4} = {'w','s'}; + conn{1,5} = {'s','w'}; + + conn{2,3} = {'n','e'}; + conn{3,4} = {'w','e'}; + conn{4,5} = {'w','s'}; + conn{5,2} = {'n','s'}; + + boundaryGroups = struct(); + boundaryGroups.E = multiblock.BoundaryGroup({2,'e'}); + boundaryGroups.N = multiblock.BoundaryGroup({3,'n'}); + boundaryGroups.W = multiblock.BoundaryGroup({4,'n'}); + boundaryGroups.S = multiblock.BoundaryGroup({5,'e'}); + boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}}); + + obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames); + + obj.r = r; + obj.c = c; + obj.hs = hs; + obj.r_arc = r_arc; + obj.omega = omega; + end + + function ms = getGridSizes(obj, m) + m_S = m; + + % m_Radial + s = 2*obj.hs; + innerArc = obj.r_arc*obj.omega; + outerArc = obj.r*pi/2; + shortSpoke = obj.r-s/sqrt(2); + x = (1/(2-sqrt(3))-1)*obj.hs; + longSpoke = (obj.r+x)-obj.r_arc; + m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2); + + ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]}; + end + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/+domain/Rectangle.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Rectangle.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,188 @@ +classdef Rectangle < multiblock.Definition + properties + + blockTi % Transfinite interpolation objects used for plotting + xlims + ylims + blockNames % Cell array of block labels + nBlocks + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + + end + + + methods + % Creates a divided rectangle + % x and y are vectors of boundary and interface positions. + % blockNames: cell array of labels. The id is default. + function obj = Rectangle(x,y,blockNames) + default_arg('blockNames',[]); + + n = length(y)-1; % number of blocks in the y direction. + m = length(x)-1; % number of blocks in the x direction. + N = n*m; % number of blocks + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + if ~issorted(flip(y)) + error('The elements of y seem to be in the wrong order'); + end + + % Dimensions of blocks and number of points + blockTi = cell(N,1); + xlims = cell(N,1); + ylims = cell(N,1); + for i = 1:n + for j = 1:m + p1 = [x(j), y(i+1)]; + p2 = [x(j+1), y(i)]; + I = flat_index(m,j,i); + blockTi{I} = parametrization.Ti.rectangle(p1,p2); + xlims{I} = {x(j), x(j+1)}; + ylims{I} = {y(i+1), y(i)}; + end + end + + % Interface couplings + conn = cell(N,N); + for i = 1:n + for j = 1:m + I = flat_index(m,j,i); + if i < n + J = flat_index(m,j,i+1); + conn{I,J} = {'s','n'}; + end + + if j < m + J = flat_index(m,j+1,i); + conn{I,J} = {'e','w'}; + end + end + end + + % Block names (id number as default) + if isempty(blockNames) + obj.blockNames = cell(1, N); + for i = 1:N + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == N); + obj.blockNames = blockNames; + end + nBlocks = N; + + % Boundary groups + boundaryGroups = struct(); + nx = m; + ny = n; + E = cell(1,ny); + W = cell(1,ny); + S = cell(1,nx); + N = cell(1,nx); + for i = 1:ny + E_id = flat_index(m,nx,i); + W_id = flat_index(m,1,i); + E{i} = {E_id,'e'}; + W{i} = {W_id,'w'}; + end + for j = 1:nx + S_id = flat_index(m,j,ny); + N_id = flat_index(m,j,1); + S{j} = {S_id,'s'}; + N{j} = {N_id,'n'}; + end + boundaryGroups.E = multiblock.BoundaryGroup(E); + boundaryGroups.W = multiblock.BoundaryGroup(W); + boundaryGroups.S = multiblock.BoundaryGroup(S); + boundaryGroups.N = multiblock.BoundaryGroup(N); + boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]); + boundaryGroups.WS = multiblock.BoundaryGroup([W,S]); + boundaryGroups.WN = multiblock.BoundaryGroup([W,N]); + boundaryGroups.ES = multiblock.BoundaryGroup([E,S]); + boundaryGroups.EN = multiblock.BoundaryGroup([E,N]); + + obj.connections = conn; + obj.nBlocks = nBlocks; + obj.boundaryGroups = boundaryGroups; + obj.blockTi = blockTi; + obj.xlims = xlims; + obj.ylims = ylims; + + end + + + % Returns a multiblock.Grid given some parameters + % ms: cell array of [mx, my] vectors + % For same [mx, my] in every block, just input one vector. + function g = getGrid(obj, ms, varargin) + + default_arg('ms',[21,21]) + + % Extend ms if input is a single vector + if (numel(ms) == 2) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockTi{i}.show(2,2); + end + axis equal + return + end + + if gridLines + m = 10; + for i = 1:obj.nBlocks + obj.blockTi{i}.show(m,m); + 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.blockTi{i}, labels{i}); + end + + axis equal + end + + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + function ms = getGridSizes(obj, varargin) + end + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/Def.m --- a/+multiblock/Def.m Tue May 22 13:29:47 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,101 +0,0 @@ -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 - -
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/DefCurvilinear.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DefCurvilinear.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,102 @@ +classdef DefCurvilinear < multiblock.Definition + 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 = DefCurvilinear(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 h = show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + h = []; + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + h = [h, obj.blockMaps{i}.show(2,2)]; + end + axis equal + return + end + + if gridLines + ms = obj.getGridSizes(varargin{:}); + for i = 1:obj.nBlocks + h = [h, 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 + case 'none' + 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 + +
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/Definition.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Definition.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,11 @@ +classdef Definition + methods (Abstract) + + % Returns a multiblock.Grid given some parameters + g = getGrid(obj, varargin) + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + show(obj, label, gridLines, varargin) + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/DiffOp.m --- a/+multiblock/DiffOp.m Tue May 22 13:29:47 2018 -0700 +++ b/+multiblock/DiffOp.m Tue Jul 24 20:14:29 2018 -0700 @@ -10,13 +10,13 @@ end methods - function obj = DiffOp(doHand, grid, order, doParam) + function obj = DiffOp(doHand, g, 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 + % g -- 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 @@ -26,9 +26,9 @@ % extra parameters to all doHand: doHand(..., doParam{:}) default_arg('doParam', []) - [getHand, getParam] = parseInput(doHand, grid, doParam); + [getHand, getParam] = parseInput(doHand, g, doParam); - nBlocks = grid.nBlocks(); + nBlocks = g.nBlocks(); obj.order = order; @@ -40,7 +40,7 @@ if ~iscell(p) p = {p}; end - obj.diffOps{i} = h(grid.grids{i}, order, p{:}); + obj.diffOps{i} = h(g.grids{i}, order, p{:}); end @@ -53,7 +53,7 @@ % Build the differentiation matrix - obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + obj.blockmatrixDiv = {g.Ns, g.Ns}; D = blockmatrix.zero(obj.blockmatrixDiv); for i = 1:nBlocks D{i,i} = obj.diffOps{i}.D; @@ -61,7 +61,7 @@ for i = 1:nBlocks for j = 1:nBlocks - intf = grid.connections{i,j}; + intf = g.connections{i,j}; if isempty(intf) continue end @@ -77,14 +77,15 @@ end end obj.D = blockmatrix.toMatrix(D); + obj.grid = g; - function [getHand, getParam] = parseInput(doHand, grid, doParam) - if ~isa(grid, 'multiblock.Grid') + function [getHand, getParam] = parseInput(doHand, g, doParam) + if ~isa(g, 'multiblock.Grid') error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); end - if iscell(doHand) && length(doHand) == grid.nBlocks() + if iscell(doHand) && length(doHand) == g.nBlocks() getHand = @(i)doHand{i}; elseif isa(doHand, 'function_handle') getHand = @(i)doHand; @@ -104,7 +105,7 @@ % doParam is a non-empty cell-array - if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) + if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam)) % doParam is a cell-array of cell-arrays getParam = @(i)doParam{i}; return @@ -116,7 +117,7 @@ function ops = splitOp(obj, op) % Splits a matrix operator into a cell-matrix of matrix operators for - % each grid. + % each g. ops = sparse2cell(op, obj.NNN); end @@ -143,6 +144,27 @@ end end + function op = getBoundaryQuadrature(obj, boundary) + opName = 'H'; + switch class(boundary) + case 'cell' + localOpName = [opName '_' boundary{2}]; + blockId = boundary{1}; + op = obj.diffOps{blockId}.(localOpName); + + return + case 'multiblock.BoundaryGroup' + N = length(boundary); + H_bm = cell(N,N); + for i = 1:N + H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i}); + end + op = blockmatrix.toMatrix(H_bm); + 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 @@ -194,6 +216,7 @@ p{I} = blockPenalty; penalty = blockmatrix.toMatrix(p); else + % TODO: used by beam equation, should be eliminated. SHould only set one BC per call for i = 1:length(blockPenalty) div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector p = blockmatrix.zero(div);
diff -r 005a8d071da3 -r 69ab0e69f972 +multiblock/Laplace.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Laplace.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,56 @@ +classdef Laplace < scheme.Scheme + properties + grid + order + mbDiffOp + + D + H + J + end + methods + function obj = Laplace(g, order, a, b, opGen) + default_arg('order', 4); + default_arg('a', 1); + default_arg('b', 1); + default_arg('opGen', @sbp.D4Variable); + + obj.grid = g; + obj.order = order; + obj.mbDiffOp = multiblock.DiffOp(@scheme.LaplaceCurvilinear, obj.grid, order, {a,b,opGen}); + + obj.D = obj.mbDiffOp.D; + obj.J = obj.jacobian(); + obj.H = obj.mbDiffOp.H; + end + + function s = size(obj) + s = size(obj.mbDiffOp); + end + + function J = jacobian(obj) + N = obj.grid.nBlocks; + J = cell(N,N); + + for i = 1:N + J{i,i} = obj.mbDiffOp.diffOps{i}.J; + end + J = blockmatrix.toMatrix(J); + end + + function op = getBoundaryOperator(obj, opName, boundary) + op = getBoundaryOperator(obj.mbDiffOp, opName, boundary); + end + + function op = getBoundaryQuadrature(obj, boundary) + op = getBoundaryQuadrature(obj.mbDiffOp, boundary); + end + + function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition + [closure, penalty] = boundary_condition(obj.mbDiffOp, boundary, type); + end + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('Not implemented') + end + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +noname/Animation.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/Animation.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,75 @@ +classdef Animation < handle + properties + timeStepper + representationMaker + updaters + end + + % add input validation + + methods + function obj = Animation(timeStepper, representationMaker, updaters); + obj.timeStepper = timeStepper; + obj.updaters = updaters; + obj.representationMaker = representationMaker; + end + + function update(obj, r) + for i = 1:length(obj.updaters) + obj.updaters{i}(r); + end + drawnow + end + + function run(obj, tEnd, timeModifier, do_pause) + default_arg('do_pause', false) + + function next_t = G(next_t) + obj.timeStepper.evolve(next_t); + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + if do_pause + pause + end + end + + anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier); + end + + function step(obj, tEnd, do_pause) + default_arg('do_pause', false) + + while obj.timeStepper.t < tEnd + obj.timeStepper.step(); + + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + % TODO: Make it never go faster than a certain fram rate + + if do_pause + pause + end + end + end + + function saveMovie(obj, tEnd, timeModifier, figureHandle, dirname) + save_frame = anim.setup_fig_mov(figureHandle, dirname); + + function next_t = G(next_t) + obj.timeStepper.evolve(next_t); + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + save_frame(); + end + + fprintf('Generating and saving frames to: ..\n') + anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier); + fprintf('Generating movies...\n') + cmd = sprintf('bash %s/+anim/make_movie.sh %s', sbplibLocation(),dirname); + system(cmd); + end + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +noname/calculateErrors.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/calculateErrors.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,53 @@ +% [discr, trueSolution] = schemeFactory(m) +% where trueSolution should be a timeSnapshot of the true solution a time T +% T is the end time +% m are grid size parameters. +% N are number of timesteps to use for each gird size +% timeOpt are options for the timeStepper +% errorFun is a function_handle taking 2 or 3 arguments, errorFun(trueSolution, approxSolution), errorFun(trueSolution, approxSolution, discr) +function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt) + %TODO: Ability to choose paralell or not + assertType(schemeFactory, 'function_handle'); + assertNumberOfArguments(schemeFactory, 1); + assertScalar(T); + assert(length(m) == length(N), 'Vectors m and N must have the same length'); + assertType(errorFun, 'function_handle'); + + if ~ismember(nargin(errorFun), [2,3]) + error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun)); + end + + default_arg('timeOpt', struct()); + + + e = zeros(1,length(m)); + parfor i = 1:length(m) + done = timeTask('m = %3d ', m(i)); + + [discr, trueSolution] = schemeFactory(m(i)); + + timeOptTemp = timeOpt; + timeOptTemp.k = T/N(i); + ts = discr.getTimestepper(timeOptTemp); + ts.stepTo(N(i), true); + approxSolution = discr.getTimeSnapshot(ts); + + switch nargin(errorFun) + case 2 + e(i) = errorFun(trueSolution, approxSolution); + case 3 + e(i) = errorFun(trueSolution, approxSolution, discr); + end + + fprintf('e = %.4e', e(i)) + done() + end + fprintf('\n') +end + + +%% Example error function +% u_true = grid.evalOn(dr.grid, @(x,y)trueSolution(T,x,y)); +% err = u_true-u_false; +% e(i) = norm(err)/norm(u_true); +% % e(i) = sqrt(err'*d.H*d.J*err/(u_true'*d.H*d.J*u_true));
diff -r 005a8d071da3 -r 69ab0e69f972 +parametrization/Curve.m --- a/+parametrization/Curve.m Tue May 22 13:29:47 2018 -0700 +++ b/+parametrization/Curve.m Tue Jul 24 20:14:29 2018 -0700 @@ -181,8 +181,8 @@ end function D = mirror(C, a, b) - assert_size(a,[2,1]); - assert_size(b,[2,1]); + assertSize(a,[2,1]); + assertSize(b,[2,1]); g = C.g; gp = C.gp; @@ -219,8 +219,8 @@ end function D = rotate(C,a,rad) - assert_size(a, [2,1]); - assert_size(rad, [1,1]); + assertSize(a, [2,1]); + assertSize(rad, [1,1]); g = C.g; gp = C.gp;
diff -r 005a8d071da3 -r 69ab0e69f972 +parametrization/Ti.m --- a/+parametrization/Ti.m Tue May 22 13:29:47 2018 -0700 +++ b/+parametrization/Ti.m Tue Jul 24 20:14:29 2018 -0700 @@ -21,16 +21,29 @@ D = g4(0); function o = S_fun(u,v) + if isrow(u) && isrow(v) + flipped = false; + else + flipped = true; + u = u'; + v = v'; + end + x1 = g1(u); x2 = g2(v); x3 = g3(1-u); x4 = g4(1-v); + o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ... - -((1-u)*(1-v).*A(1,:) + u*(1-v).*B(1,:) + u*v.*C(1,:) + (1-u)*v.*D(1,:)); + -((1-u).*(1-v).*A(1,:) + u.*(1-v).*B(1,:) + u.*v.*C(1,:) + (1-u).*v.*D(1,:)); o2 = (1-v).*x1(2,:) + u.*x2(2,:) + v.*x3(2,:) + (1-u).*x4(2,:) ... - -((1-u)*(1-v).*A(2,:) + u*(1-v).*B(2,:) + u*v.*C(2,:) + (1-u)*v.*D(2,:)); + -((1-u).*(1-v).*A(2,:) + u.*(1-v).*B(2,:) + u.*v.*C(2,:) + (1-u).*v.*D(2,:)); - o = [o1;o2]; + if ~flipped + o = [o1;o2]; + else + o = [o1'; o2']; + end end obj.S = @S_fun; @@ -116,13 +129,13 @@ S = obj.S; if(nu>2 || nv>2) - h_grid = obj.plot(nu,nv); - set(h_grid,'Color',[0 0.4470 0.7410]); + h.grid = obj.plot(nu,nv); + set(h.grid,'Color',[0 0.4470 0.7410]); end - h_bord = obj.plot(2,2); - set(h_bord,'Color',[0.8500 0.3250 0.0980]); - set(h_bord,'LineWidth',2); + h.border = obj.plot(2,2); + set(h.border,'Color',[0.8500 0.3250 0.0980]); + set(h.border,'LineWidth',2); end
diff -r 005a8d071da3 -r 69ab0e69f972 +parametrization/TiTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/TiTest.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,52 @@ +function tests = TiTest() + tests = functiontests(localfunctions); +end + +function testScalarInput(testCase) + ti = getMinimumTi(); + + cases = { + % {u, v, out}, + {0, 0, [1; 2]}, + {0, 1, [1; 4]}, + {1, 0, [3; 2]}, + {1, 1, [3; 4]}, + {0.5, 0.5, [2; 3]}, + }; + + for i = 1:length(cases) + u = cases{i}{1}; + v = cases{i}{2}; + expected = cases{i}{3}; + + testCase.verifyEqual(ti.S(u,v), expected, sprintf('Case: %d',i)); + end +end + +function testRowVectorInput(testCase) + ti = getMinimumTi(); + + u = [0, 0.5, 1]; + v = [0, 0, 0.5]; + expected = [ + 1, 2, 3; + 2, 2, 3; + ]; + + testCase.verifyEqual(ti.S(u,v), expected); +end + +function testColumnvectorInput(testCase) + ti = getMinimumTi(); + + u = [0; 0.5; 1]; + v = [0; 0; 0.5]; + expected = [1; 2; 3; 2; 2; 3]; + + testCase.verifyEqual(ti.S(u,v), expected); +end + + +function ti = getMinimumTi() + ti = parametrization.Ti.rectangle([1; 2], [3; 4]); +end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 +scheme/Beam.m --- a/+scheme/Beam.m Tue May 22 13:29:47 2018 -0700 +++ b/+scheme/Beam.m Tue Jul 24 20:14:29 2018 -0700 @@ -126,6 +126,44 @@ penalty{1} = -obj.Hi*tau; penalty{1} = -obj.Hi*sig; + case 'e' + alpha = obj.alpha; + tuning = 1.1; + + tau1 = tuning * alpha/delt; + tau4 = s*alpha; + + tau = tau1*e+tau4*d3; + + closure = obj.Hi*tau*e'; + penalty = -obj.Hi*tau; + case 'd1' + alpha = obj.alpha; + + tuning = 1.1; + + sig2 = tuning * alpha/gamm; + sig3 = -s*alpha; + + sig = sig2*d1+sig3*d2; + + closure = obj.Hi*sig*d1'; + penalty = -obj.Hi*sig; + + case 'd2' + a = obj.alpha; + + tau = s*a*d1; + + closure = obj.Hi*tau*d2'; + penalty = -obj.Hi*tau; + case 'd3' + a = obj.alpha; + + sig = -s*a*e; + + closure = obj.Hi*sig*d3'; + penalty = -obj.Hi*sig; otherwise % Unknown, boundary condition error('No such boundary condition: type = %s',type);
diff -r 005a8d071da3 -r 69ab0e69f972 +scheme/TODO.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/TODO.txt Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,1 @@ +% TODO: Rename package and abstract class to diffOp
diff -r 005a8d071da3 -r 69ab0e69f972 +scheme/bcSetup.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/bcSetup.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,75 @@ +% function [closure, S] = bcSetup(diffOp, bc) +% Takes a diffOp and a cell array of boundary condition definitions. +% Each bc is a struct with the fields +% * type -- Type of boundary condition +% * boundary -- Boundary identifier +% * data -- A function_handle for a function which provides boundary data.(see below) +% Also takes S_sign which modifies the sign of S, [-1,1] +% Returns a closure matrix and a forcing function S. +% +% The boundary data function can either be a function of time or a function of time and space coordinates. +% In the case where it only depends on time it should return the data as grid function for the boundary. +% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain. +% For example in the 2D case: f(t,x,y). +function [closure, S] = bcSetup(diffOp, bc, S_sign) + default_arg('S_sign', 1); + assertType(bc, 'cell'); + assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); + + + % Setup storage arrays + closure = spzeros(size(diffOp)); + gridDataPenalties = {}; + gridDataFunctions = {}; + symbolicDataPenalties = {}; + symbolicDataFunctions = {}; + symbolicDataCoords = {}; + + % Collect closures, penalties and data + for i = 1:length(bc) + assertType(bc{i}, 'struct'); + [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); + closure = closure + localClosure; + + if ~isfield(bc{i},'data') || isempty(bc{i}.data) + % Skip to next loop if there is no data + continue + end + assertType(bc{i}.data, 'function_handle'); + + % Find dimension + dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2); + + if nargin(bc{i}.data) == 1 + % Grid data + boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1]; + assertSize(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size. + gridDataPenalties{end+1} = penalty; + gridDataFunctions{end+1} = bc{i}.data; + elseif nargin(bc{i}.data) == 1+dim + % Symbolic data + coord = diffOp.grid.getBoundary(bc{i}.boundary); + symbolicDataPenalties{end+1} = penalty; + symbolicDataFunctions{end+1} = bc{i}.data; + symbolicDataCoords{end+1} = num2cell(coord ,1); + else + error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + end + end + + % Setup penalty function + O = spzeros(size(diffOp),1); + function v = S_fun(t) + v = O; + for i = 1:length(gridDataFunctions) + v = v + gridDataPenalties{i}*gridDataFunctions{i}(t); + end + + for i = 1:length(symbolicDataFunctions) + v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:}); + end + + v = S_sign * v; + end + S = @S_fun; +end
diff -r 005a8d071da3 -r 69ab0e69f972 +time/SBPInTimeScaled.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeScaled.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,139 @@ +classdef SBPInTimeScaled < time.Timestepper + % The SBP in time method. + % Implemented for A*v_t = B*v + f(t), v(0) = v0 + % The resulting system of equations is + % M*u_next= K*u_prev_end + f + properties + A,B + f + + k % total time step. + + blockSize % number of points in each block + N % Number of components + + order + nodes + + Mtilde,Ktilde % System matrices + L,U,p,q % LU factorization of M + e_T + + scaling + S, Sinv % Scaling matrices + + % Time state + t + vtilde + n + end + + methods + function obj = SBPInTimeScaled(A, B, f, k, t0, v0, scaling, TYPE, order, blockSize) + default_arg('TYPE','gauss'); + default_arg('f',[]); + + if(strcmp(TYPE,'gauss')) + default_arg('order',4) + default_arg('blockSize',4) + else + default_arg('order', 8); + default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE)); + end + + obj.A = A; + obj.B = B; + obj.scaling = scaling; + + if ~isempty(f) + obj.f = f; + else + obj.f = @(t)sparse(length(v0),1); + end + + obj.k = k; + obj.blockSize = blockSize; + obj.N = length(v0); + + obj.n = 0; + obj.t = t0; + + %==== Build the time discretization matrix =====% + switch TYPE + case 'equidistant' + ops = sbp.D2Standard(blockSize,{0,obj.k},order); + case 'optimal' + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order); + case 'minimal' + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal'); + case 'gauss' + ops = sbp.D1Gauss(blockSize,{0,obj.k}); + end + + I = speye(size(A)); + I_t = speye(blockSize,blockSize); + + D1 = kron(ops.D1, I); + HI = kron(ops.HI, I); + e_0 = kron(ops.e_l, I); + e_T = kron(ops.e_r, I); + obj.nodes = ops.x; + + % Convert to form M*w = K*v0 + f(t) + tau = kron(I_t, A) * e_0; + M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B); + + K = HI*tau; + + obj.S = kron(I_t, spdiag(scaling)); + obj.Sinv = kron(I_t, spdiag(1./scaling)); + + obj.Mtilde = obj.Sinv*M*obj.S; + obj.Ktilde = obj.Sinv*K*spdiag(scaling); + obj.e_T = e_T; + + + % LU factorization + [obj.L,obj.U,obj.p,obj.q] = lu(obj.Mtilde, 'vector'); + + obj.vtilde = (1./obj.scaling).*v0; + end + + function [v,t] = getV(obj) + v = obj.scaling.*obj.vtilde; + t = obj.t; + end + + function obj = step(obj) + forcing = zeros(obj.blockSize*obj.N,1); + + for i = 1:obj.blockSize + forcing((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i)); + end + + RHS = obj.Sinv*forcing + obj.Ktilde*obj.vtilde; + + y = obj.L\RHS(obj.p); + z = obj.U\y; + + w = zeros(size(z)); + w(obj.q) = z; + + obj.vtilde = obj.e_T'*w; + + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + methods(Static) + function N = smallestBlockSize(order,TYPE) + default_arg('TYPE','gauss') + + switch TYPE + case 'gauss' + N = 4; + end + end + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 +time/SBPInTimeSecondOrderFormImplicit.m --- a/+time/SBPInTimeSecondOrderFormImplicit.m Tue May 22 13:29:47 2018 -0700 +++ b/+time/SBPInTimeSecondOrderFormImplicit.m Tue Jul 24 20:14:29 2018 -0700 @@ -14,15 +14,16 @@ % Solves A*u_tt + B*u_t + C*u = f(t) % A, B can either both be constants or both be function handles, % They can also be omitted by setting them equal to the empty matrix. - function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize) + function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, do_scaling, TYPE, order, blockSize) default_arg('f', []); default_arg('TYPE', []); default_arg('order', []); default_arg('blockSize',[]); + default_arg('do_scaling', false); m = length(v0); - default_arg('A', sparse(m, m)); + default_arg('A', speye(m, m)); default_arg('B', sparse(m, m)); default_arg('C', sparse(m, m)); @@ -56,7 +57,12 @@ obj.t = t0; obj.n = 0; - obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize); + if do_scaling + scaling = [ones(m,1); sqrt(diag(C))]; + obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize); + else + obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize); + end end function [v,t] = getV(obj)
diff -r 005a8d071da3 -r 69ab0e69f972 +time/Timestepper.m --- a/+time/Timestepper.m Tue May 22 13:29:47 2018 -0700 +++ b/+time/Timestepper.m Tue Jul 24 20:14:29 2018 -0700 @@ -62,6 +62,7 @@ function [v, t] = stepTo(obj, n, progress_bar) + assertScalar(n); default_arg('progress_bar',false); [v, t] = obj.stepN(n-obj.n, progress_bar);
diff -r 005a8d071da3 -r 69ab0e69f972 .hgtags --- a/.hgtags Tue May 22 13:29:47 2018 -0700 +++ b/.hgtags Tue Jul 24 20:14:29 2018 -0700 @@ -1,1 +1,2 @@ 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1 +0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1
diff -r 005a8d071da3 -r 69ab0e69f972 Color.m --- a/Color.m Tue May 22 13:29:47 2018 -0700 +++ b/Color.m Tue Jul 24 20:14:29 2018 -0700 @@ -10,6 +10,10 @@ black = [0.000 0.000 0.000]; white = [1.000 1.000 1.000]; colors = { Color.blue, Color.red, Color.yellow, Color.green, Color.purple, Color.lightblue, Color.darkred, Color.black, Color.white}; + markers = {'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'}; + lineStyles = {'-', '--', ':', '-.'}; + + solidMarkers = {'o', 'square', 'diamond', 'v', 'pentagram', '^', '>', '<', 'hexagram'}; notabilityYellow = [100.0 99.0 22.0 ]/100; notabilityOrange = [97.0 61.0 15.0 ]/100; @@ -34,13 +38,11 @@ methods(Static) function sample() - markers ={'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'}; % Filled and non-filled markers? - lineStyles = {'-', '--', ':', '-.'}; function showMarkers(x0, y0, lx, ly, color, filled) - n = length(markers); + n = length(Color.markers); s = ceil(sqrt(n)); x = linspace(x0, x0 + lx, s); @@ -50,7 +52,7 @@ for i = 1:n lh = line(X(i),Y(i)); - lh.Marker = markers{i}; + lh.Marker = Color.markers{i}; lh.MarkerSize = 12; lh.Color = color; @@ -79,13 +81,13 @@ end function showLines(y0, ly, A, w) - n = length(lineStyles); + n = length(Color.lineStyles); x = linspace(0,1,100); y = linspace(y0, y0+ ly, n); for i = 1:n lh = line(x, y(i) + A*sin(pi*x*w)); lh.LineWidth = 2; - lh.LineStyle = lineStyles{i}; + lh.LineStyle = Color.lineStyles{i}; end end
diff -r 005a8d071da3 -r 69ab0e69f972 TextTable.m --- a/TextTable.m Tue May 22 13:29:47 2018 -0700 +++ b/TextTable.m Tue Jul 24 20:14:29 2018 -0700 @@ -4,28 +4,36 @@ fmtArray vertDiv horzDiv - - nCols - nRows end methods - function obj = TextTable(data, vertDiv, horzDiv); + 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 n = nRows(obj) + n = size(obj.data, 1); + end + + function n = nCols(obj) + n = size(obj.data, 2); + end + + function print(obj) + disp(obj.toString()); + end + function formatAll(obj, fmt) + obj.fmtArray = cell(size(obj.data)); obj.fmtArray(:,:) = {fmt}; end @@ -58,28 +66,31 @@ str = ''; + N = size(strArray, 2); + % First horzDiv - if ismember(0, obj.horzDiv) + if isDiv(0, obj.horzDiv, N); str = [str, obj.getHorzDiv(widths)]; end for i = 1:obj.nRows - str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv, obj.horzDiv)]; + str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv)]; % Interior horzDiv - if ismember(i, obj.horzDiv) + if isDiv(i, obj.horzDiv, N) 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 = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv); str(find(' ' == str)) = '-'; str(find('|' == str)) = '+'; end function strArray = getStringArray(obj) + assert(all(size(obj.data) == size(obj.fmtArray)), 'Sizes of format matrix and data matrix do not match.') strArray = cell(size(obj.data)); for i = 1:obj.nRows @@ -91,32 +102,42 @@ end methods (Static) - function str = rowToString(strs, widths, vertDiv, horzDiv) + function str = rowToString(strs, widths, vertDiv) + N = length(strs); + % First vertDiv - if ismember(0, vertDiv) - str = '| '; + if isDiv(0, vertDiv, N) + prefix = '| '; else - str = ' '; + prefix = ' '; end - % Interior cols - for j = 1:length(strs) - 1 - str = [str, sprintf('%*s ', widths(j), strs{j})]; + % Pad strings + for i = 1:N + strs{i} = sprintf('%*s', widths(i), strs{i}); + end - % Interior vertDiv - if ismember(j, vertDiv) - str = [str, '| ']; + % Column delimiters + delims = cell(1,N-1); + for i = 1:length(delims) + if isDiv(i, vertDiv, N); + delims{i} = '| '; + else + delims{i} = ' '; end end - % Last col - str = [str, sprintf('%*s ', widths(end), strs{end})]; - - if ismember(length(strs), vertDiv) - str = [str, '|']; + if isDiv(N, vertDiv, N); + suffix = '|'; + else + suffix = ''; end - str = [str, sprintf('\n')]; + str = [prefix, strjoin(strs, delims), suffix, sprintf('\n')]; end end +end + +function b = isDiv(i, div, N) + b = ismember(i, div) || ismember(i, N+div+1); end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 assertIsMember.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertIsMember.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,3 @@ +function assertIsMember(v, allowed) + assert(ismember(v, allowed), 'Expected ''%s'' to be in the set %s', inputname(1), toString(allowed)); +end
diff -r 005a8d071da3 -r 69ab0e69f972 assertScalar.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertScalar.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,5 @@ +function assertScalar(obj) + if ~isscalar(obj) + error('sbplib:assertScalar:notScalar', '"%s" must be scalar, found size "%s"', inputname(1), toString(size(obj))); + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 assertSize.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertSize.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,16 @@ +% Assert that array A has the size s. +function assertSize(A,varargin) + if length(varargin) == 1 + s = varargin{1}; + errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), toString(s), toString(size(A))); + assert(all(size(A) == s), errmsg); + elseif length(varargin) == 2 + dim = varargin{1}; + s = varargin{2}; + + errmsg = sprintf('Expected %s to have size %d along dimension %d, got: %d',inputname(1), s, dim, size(A,dim)); + assert(size(A,dim) == s, errmsg); + else + error('Expected 2 or 3 arguments to assertSize()'); + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 assertStructFields.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertStructFields.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,12 @@ +% Assert that the struct s has the all the field names in the cell array fns. +function assertStructFields(s, fns) + assertType(s, 'struct'); + assertType(fns, 'cell'); + + ok = ismember(fns, fieldnames(s)); + if ~all(ok) + str1 = sprintf("'%s' must have the fields %s\n", inputname(1), toString(fns)); + str2 = sprintf("The following fields are missing: %s", toString(fns(~ok))); + error(str1 + str2); + end +end
diff -r 005a8d071da3 -r 69ab0e69f972 assert_size.m --- a/assert_size.m Tue May 22 13:29:47 2018 -0700 +++ b/assert_size.m Tue Jul 24 20:14:29 2018 -0700 @@ -1,16 +1,5 @@ % Assert that array A has the size s. function assert_size(A,s) - errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), format_vector(s), format_vector(size(A))); - assert(all(size(A) == s),errmsg); -end - -function str = format_vector(a) - l = length(a); - str = sprintf('[%d',a(1)); - - for i = 2:l - str = [str sprintf(', %d',a(i))]; - end - - str = [str ']']; + warning('Use assertSize() instead!') + assertSize(A,s); end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 convergencePlot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/convergencePlot.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,55 @@ +function hand = convergencePlot(orders, h, e) + N = length(orders); + + fh = figure(); + ah = axes(); + ah.XScale = 'log'; + ah.YScale = 'log'; + hold on + ph = {}; + phc = {}; + legends = {}; + for i = 1:N + ph{i} = loglog(h{i}, e{i}); + phc{i} = plotConvergenceFit(orders{i}, h{i}, e{i}); + + ph{i}.LineStyle = 'none'; + ph{i}.Marker = Color.solidMarkers{i}; + ph{i}.MarkerSize = 12; + ph{i}.Color = Color.colors{i}; + ph{i}.MarkerFaceColor = Color.colors{i}; + + legends{i} = sprintf('$o = %d$', orders{i}); + end + hold off + + lh = legend([ph{:}], legends); + lh.Interpreter = 'latex'; + lh.Location = 'SouthEast'; + + for i = 1:N + uistack(phc{i}, 'bottom'); + end + + xlabel('$h$', 'interpreter', 'latex') + ylabel('Error', 'interpreter', 'latex') + + % xlim([0.7e-2, 1e-1]) + % ylim([3e-5, 4]) + + grid on + + ah = gca(); + ah.TickLabelInterpreter = 'latex'; + setFontSize(fh); + + % if savePngs + % savepng(fh, 'fig/conv/conv',600) + % end + + hand = struct(); + hand.fig = fh; + hand.data = ph; + hand.fits = phc; + hand.legend = lh; +end
diff -r 005a8d071da3 -r 69ab0e69f972 diffSymfun.m --- a/diffSymfun.m Tue May 22 13:29:47 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -% Differentiates a symbolic function like diff does, but keeps the function as a symfun -function g = diffSymfun(f, varargin) - assertType(f, 'symfun'); - - args = argnames(f); - g = symfun(diff(f,varargin{:}), args); -end
diff -r 005a8d071da3 -r 69ab0e69f972 gaussian.m --- a/gaussian.m Tue May 22 13:29:47 2018 -0700 +++ b/gaussian.m Tue Jul 24 20:14:29 2018 -0700 @@ -1,3 +1,3 @@ function z = gaussian(x,x0,d) - z = exp(-norm(x-x0).^2/d^2); + z = exp(-sum((x-x0).^2,2)/d^2); end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 mononomial.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mononomial.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,17 @@ +% calculate a N-D mononomial with powers k in points x: +% z = x(:,1).^k(1) * x(:,2).^k(2) * ... +function z = mononomial(x, k) + assert(size(x,2) == length(k), 'k must have the same length as the width of x'); + + if any(k < 0) + z = x(:,1)*0; + return + end + + denom = prod(factorial(k)); + + for i = 1:length(k) + x(:,i) = x(:,i).^k(i); + end + z = prod(x,2)/denom; +end
diff -r 005a8d071da3 -r 69ab0e69f972 nextColor.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nextColor.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,5 @@ +function c = nextColor(ah) + default_arg('ah', gca); + + c = ah.ColorOrder(ah.ColorOrderIndex, :); +end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 pointIndex.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pointIndex.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,4 @@ +% Get the index of the points p within the tall array of points ps +function [I, ok] = pointIndex(p, ps) + [ok, I] = ismember(p, ps, 'rows'); +end
diff -r 005a8d071da3 -r 69ab0e69f972 rickerWavelet.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rickerWavelet.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,3 @@ +function y = rickerWavelet(x, x0, A) + y = (1-2*pi^2*A^2*(x-x0).^2).*exp(-pi^2*A^2*(x-x0).^2); +end
diff -r 005a8d071da3 -r 69ab0e69f972 sbplibLocation.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sbplibLocation.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,4 @@ +function location = sbplibLocation() + scriptname = mfilename('fullpath'); + [location, ~, ~] = fileparts(scriptname); +end
diff -r 005a8d071da3 -r 69ab0e69f972 sbplibVersion.m --- a/sbplibVersion.m Tue May 22 13:29:47 2018 -0700 +++ b/sbplibVersion.m Tue Jul 24 20:14:29 2018 -0700 @@ -1,11 +1,10 @@ % Prints the version and location of the sbplib currently in use. function sbplibVersion() - scriptname = mfilename('fullpath'); - [folder,~,~] = fileparts(scriptname); + location = sbplibLocation(); name = 'sbplib (feature/grids)'; ver = '0.0.x'; fprintf('%s %s\n', name, ver); - fprintf('Running in:\n%s\n',folder); + fprintf('Running in:\n%s\n', location); end \ No newline at end of file
diff -r 005a8d071da3 -r 69ab0e69f972 stencilEquation.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stencilEquation.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,13 @@ +% Find the equation for the stencil for d^k/dx^k +function [A,b] = stencilEquation(k, offsets, order) + q = sym('q', [1, length(offsets)]); + + p = 0:(order-1+k); + + v = vandermonde(offsets, p); + vdiff = vandermonde( 0, p-k); + + eq = q*v == vdiff; + + [A,b] = equationsToMatrix(eq, q); +end
diff -r 005a8d071da3 -r 69ab0e69f972 vandermonde.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/vandermonde.m Tue Jul 24 20:14:29 2018 -0700 @@ -0,0 +1,15 @@ +% Create vandermonde matrix for points x and polynomials of order p +% x is a list of N points of size [N,dim], +% p is a list of polynomial orders of size [M, dim]. +% the given mononomials are evaluated and the NxM matrix V is returned. +function V = vandermonde(x, p) + assert(size(x,2) == size(p,2), 'x and p must have the same number of columns') + n = size(x,1); + m = size(p,1); + + for i = 1:m + V(:,i) = mononomial(x, p(i,:)); + end + + assertSize(V,[n,m]); +end