Mercurial > repos > public > sbplib
changeset 797:5cf9fdf4c98f feature/poroelastic
Merge with feature/grids and bugfix bcSetup
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 26 Jul 2018 10:53:05 -0700 |
| parents | aa1ed37a1b56 (current diff) c7c622e26a53 (diff) |
| children | 8c65ef13df89 |
| files | +grid/Cartesian.m +multiblock/DiffOp.m +multiblock/multiblockgrid.m +multiblock/stitchSchemes.m +scheme/bcSetup.m .hgtags |
| diffstat | 48 files changed, 1307 insertions(+), 272 deletions(-) [+] |
line wrap: on
line diff
--- a/+anim/setup_time_quantity_plot.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+anim/setup_time_quantity_plot.m Thu Jul 26 10:53:05 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));
--- a/+draw/prompt_point.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+draw/prompt_point.m Thu Jul 26 10:53:05 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')
--- a/+grid/Cartesian.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+grid/Cartesian.m Thu Jul 26 10:53:05 2018 -0700 @@ -16,7 +16,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});
--- a/+multiblock/DefCurvilinear.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+multiblock/DefCurvilinear.m Thu Jul 26 10:53:05 2018 -0700 @@ -48,13 +48,14 @@ g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); end - function show(obj, label, gridLines, varargin) + 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 - obj.blockMaps{i}.show(2,2); + h = [h, obj.blockMaps{i}.show(2,2)]; end axis equal return @@ -63,7 +64,7 @@ if gridLines ms = obj.getGridSizes(varargin{:}); for i = 1:obj.nBlocks - obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + h = [h, obj.blockMaps{i}.show(ms{i}(1),ms{i}(2))]; end end @@ -76,7 +77,7 @@ for i = 1:obj.nBlocks labels{i} = num2str(i); end - otherwise + case 'none' axis equal return end
--- a/+multiblock/DiffOp.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+multiblock/DiffOp.m Thu Jul 26 10:53:05 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 @@ -65,7 +65,7 @@ for i = 1:nBlocks for j = 1:nBlocks - intf = grid.connections{i,j}; + intf = g.connections{i,j}; if isempty(intf) continue end @@ -81,15 +81,15 @@ end end obj.D = blockmatrix.toMatrix(D); - obj.grid = grid; + 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; @@ -109,7 +109,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 @@ -148,6 +148,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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Laplace.m Thu Jul 26 10:53:05 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
--- a/+multiblock/multiblockgrid.m Wed Jul 25 18:53:07 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -% Creates a multi block square grid with defined boundary conditions. -% x,y defines the grid lines. Rember to think of the indexing as a matrix. Order matters! -% bc is a struct defining the boundary conditions on each side of the block. -% bc.w = {'dn',[function or value]} -function [block,conn,bound,ms] = multiblockgrid(x,y,mx,my,bc) - 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 - % y = sort(y,'descend'); - - % Dimensions of blocks and number of points - block = cell(n,m); - for i = 1:n - for j = 1:m - block{i,j} = { - {x(j),x(j+1)}, {y(i+1),y(i)}; - }; - - ms{i,j} = [mx(i),my(j)]; - end - end - - % Interface couplings - conn = cell(N,N); - for i = 1:n - for j = 1:m - I = flat_index(n,i,j); - if i < n - J = flat_index(n,i+1,j); - conn{I,J} = {'s','n'}; - end - - if j < m - J = flat_index(n,i,j+1); - conn{I,J} = {'e','w'}; - end - end - end - - - % Boundary conditions - bound = cell(n,m); - for i = 1:n - if isfield(bc,'w') - bound{i,1}.w = bc.w; - end - - if isfield(bc,'e') - bound{i,n}.e = bc.e; - end - end - - for j = 1:m - if isfield(bc,'n') - bound{1,j}.n = bc.n; - end - - if isfield(bc,'s') - bound{m,j}.s = bc.s; - end - end -end -
--- a/+multiblock/stitchSchemes.m Wed Jul 25 18:53:07 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,92 +0,0 @@ -% Stitch schemes together given connection matrix and BC vector. -% schmHand - function_handle to a Scheme constructor -% order - order of accuracy -% schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays -% blocks - block definitions, On whatever form the scheme expects. -% ms - grid points in each direction for each block. Ex {[10,10], [10, 20]} -% conn - connection matrix -% bound - boundary condition vector, array of structs with fields w,e,s,n -% each field with a parameter array that is sent to schm.boundary_condition -% -% Output parameters are cell arrays and cell matrices. -% -% Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) -function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound) - default_arg('schmParam',[]); - - n_blocks = numel(grids); - - % Creating Schemes - for i = 1:n_blocks - if isempty(schmParam); - schms{i} = schmHand(grids{i},order,[]); - elseif ~iscell(schmParam) - param = schmParam(i); - schms{i} = schmHand(grids{i},order,param); - else - param = schmParam{i}; - if iscell(param) - schms{i} = schmHand(grids{i},order,param{:}); - else - schms{i} = schmHand(grids{i},order,param); - end - end - - % class(schmParam) - % class(ms) - % class(blocks) - % class(schmParam{i}) - % class(ms) - - - end - - - % Total norm - H = cell(n_blocks,n_blocks); - for i = 1:n_blocks - H{i,i} = schms{i}.H; - end - - %% Total system matrix - - % Differentiation terms - D = cell(n_blocks,n_blocks); - for i = 1:n_blocks - D{i,i} = schms{i}.D; - end - - % Boundary penalty terms - for i = 1:n_blocks - if ~isstruct(bound{i}) - continue - end - - fn = fieldnames(bound{i}); - for j = 1:length(fn); - bc = bound{i}.(fn{j}); - if isempty(bc) - continue - end - - [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); - D{i,i} = D{i,i}+closure; - end - end - - % Interface penalty terms - for i = 1:n_blocks - for j = 1:n_blocks - intf = conn{i,j}; - if isempty(intf) - continue - end - - [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2}); - D{i,i} = D{i,i} + uu; - D{i,j} = uv; - D{j,j} = D{j,j} + vv; - D{j,i} = vu; - end - end -end \ No newline at end of file
--- a/+noname/calculateErrors.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+noname/calculateErrors.m Thu Jul 26 10:53:05 2018 -0700 @@ -4,27 +4,40 @@ % 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'); - assertNumberOfArguments(errorFun, 2); - default_arg('timeOpt'); + + if ~ismember(nargin(errorFun), [2,3]) + error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun)); + end - e = []; - for i = 1:length(m) + 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)); - timeOpt.k = T/N(i); - ts = discr.getTimestepper(timeOpt); + timeOptTemp = timeOpt; + timeOptTemp.k = T/N(i); + ts = discr.getTimestepper(timeOptTemp); ts.stepTo(N(i), true); approxSolution = discr.getTimeSnapshot(ts); - e(i) = errorFun(trueSolution, approxSolution); + 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()
--- a/+parametrization/Curve.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+parametrization/Curve.m Thu Jul 26 10:53:05 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;
--- a/+parametrization/Ti.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+parametrization/Ti.m Thu Jul 26 10:53:05 2018 -0700 @@ -129,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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_2to2_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2 +%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:07 + +stencil_F2C = [-1.0./8.0,1.0./4.0,3.0./4.0,1.0./4.0,-1.0./8.0]; +if nargout > 1 + BC_F2C = reshape([6.288191560529559e-1,-6.440957802647795e-2,6.855471317807184e-1,1.572264341096408e-1,-2.438028498638558e-1,7.469014249319279e-1,-8.431231982626708e-2,2.921561599131335e-1,1.374888185644862e-2,-1.318744409282243e-1],[2,5]); +end +if nargout > 2 + HcU = 1.0./2.0; +end +if nargout > 3 + HfU = 1.0./2.0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1 +%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:08 + +stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0]; +if nargout > 1 + BC_F2C = [1.0./2.0,1.0./2.0]; +end +if nargout > 2 + HcU = 1.0./2.0; +end +if nargout > 3 + HfU = 1.0./2.0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_4to4_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3 +%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:01 + +stencil_F2C = [7.0./2.56e2,-1.0./3.2e1,-7.0./6.4e1,9.0./3.2e1,8.5e1./1.28e2,9.0./3.2e1,-7.0./6.4e1,-1.0./3.2e1,7.0./2.56e2]; +if nargout > 1 + BC_F2C = reshape([7.523257802630956e-2,2.447812262221267e-1,-1.679313063616916e-1,1.290510950666589,6.315723344677289e-3,1.67178747937954e-1,1.982667903557025,-7.554383893379468e-1,7.215271362899867e-1,-2.820478831383137,1.807034411305536,-7.589683751979843e-1,-7.685973268458095e-1,3.965751544535173e-1,4.119789638051451e-1,1.574000556785898e-1,-1.113618964927466e-1,3.631000220124657e-1,1.639694219982991,-9.114074742274862e-1,4.952715520862987e-1,-4.524162151353456e-1,2.65481378213589e-1,-2.268273408674622e-1,3.455365956773523e-1,-2.009765975089827e-1,1.722179564305811e-1,-5.52933488235368e-1,3.18109976271229e-1,-2.171481232558432e-1,1.033835580108027e-1,-5.911351224351343e-2,3.960076712054991e-2],[3,11]); +end +if nargout > 2 + t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2 +%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:05 + +stencil_F2C = [-1.0./2.56e2,-1.0./3.2e1,1.0./6.4e1,9.0./3.2e1,6.1e1./1.28e2,9.0./3.2e1,1.0./6.4e1,-1.0./3.2e1,-1.0./2.56e2]; +if nargout > 1 + BC_F2C = reshape([1.0./2.0,0.0,0.0,1.77e2./2.72e2,3.0./8.0,-5.9e1./6.88e2,1.125919117647059e-2,3.546742584745763e-1,1.335392441860465e-2,-4.9e1./5.44e2,2.335805084745763e-1,3.204941860465116e-1,-1.194852941176471e-2,1.350635593220339e-2,5.308866279069767e-1,-9.0./5.44e2,-1.11228813559322e-2,2.943313953488372e-1,-2.803308823529412e-2,2.105402542372881e-2,-1.580668604651163e-2,-9.0./5.44e2,1.430084745762712e-2,-5.450581395348837e-2,-9.191176470588235e-4,7.944915254237288e-4,-5.450581395348837e-3,1.0./5.44e2,-1.588983050847458e-3,2.180232558139535e-3,2.297794117647059e-4,-1.986228813559322e-4,2.725290697674419e-4],[3,11]); +end +if nargout > 2 + t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_6to6_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4 +%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:47 + +stencil_F2C = [-5.95703125e-3,3.0./5.12e2,3.57421875e-2,-2.5e1./5.12e2,-8.935546875e-2,7.5e1./2.56e2,3.17e2./5.12e2,7.5e1./2.56e2,-8.935546875e-2,-2.5e1./5.12e2,3.57421875e-2,3.0./5.12e2,-5.95703125e-3]; +if nargout > 1 + BC_F2C = reshape([5.233890618131365e-1,-1.594459192645467e-2,3.532688727637403e-2,8.021234957689208e-1,3.90683205173562e-1,-1.73222951632239e-1,-8.87662686483442e-2,2.90091235796637e-1,-1.600356115709148e-1,-1.044025375027475e-1,2.346179009198368e-1,6.091329306528956e-1,1.561275522703128e-1,-1.168382445709856e-1,1.040987887887311,-2.387061036980731e-1,1.363504965974361e-1,1.082611928255256e-1,-5.745310654054326e-1,3.977694249198785e-1,-9.376217911402619e-1,3.554518646054656e-2,5.609157787396987e-5,-1.564625311232018e-1,6.107907974027401e-1,-3.786608696698368e-1,7.02265869951125e-1,2.054294270642538e-1,-1.302300112378257e-1,2.478941407690889e-1,-2.657085326191479e-1,1.568761445933572e-1,-2.632906518005349e-1,-1.228047556139644e-1,7.182193248980271e-2,-1.1291238242346e-1,5.811258780405158e-2,-3.466364400805378e-2,5.683680203338252e-2,1.781337575097077e-2,-1.030572999042704e-2,1.577197067502767e-2,-1.443802115768554e-2,8.391316308374261e-3,-1.29534117369052e-2,-1.548018627738296e-3,8.794183904936319e-4,-1.298961407229805e-3,1.573818938200601e-3,-8.940753636685258e-4,1.320610764016968e-3],[3,17]); +end +if nargout > 2 + t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3 +%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:56 + +stencil_F2C = [1.07421875e-3,3.0./5.12e2,-6.4453125e-3,-2.5e1./5.12e2,1.611328125e-2,7.5e1./2.56e2,2.45e2./5.12e2,7.5e1./2.56e2,1.611328125e-2,-2.5e1./5.12e2,-6.4453125e-3,3.0./5.12e2,1.07421875e-3]; +if nargout > 1 + BC_F2C = reshape([1.0./2.0,0.0,0.0,0.0,6.876076086160158e-1,1.5e1./3.2e1,-3.461879841386942e-1,3.502577439820862e-2,3.099721991995751e-3,2.228547003766753e-1,9.363652999354482e-3,-3.15791046603844e-3,-1.05789143206462e-1,2.35563165945226e-1,6.070385985337514e-1,-4.847496617839149e-2,-4.809232246867902e-3,5.154694068405061e-3,7.049223649022501e-1,1.016749349808733e-2,3.460117842882262e-2,-4.996621498584866e-2,5.210650763094799e-1,2.233592875303228e-1,-2.239024779745769e-3,4.254668119953384e-4,9.213872590833641e-3,3.910524351558127e-1,-9.048711215107334e-2,8.597375629526346e-2,-3.468219752858724e-1,3.250682992395969e-1,-1.309107466664224e-1,1.199249722306043e-1,-4.071552384959425e-1,1.474417123938701e-1,-3.02863877390285e-2,3.046016554982103e-2,-1.081315128181483e-1,1.120705238850532e-2,7.977722870036266e-2,-6.772545887788229e-2,2.00270418837606e-1,-5.427183680490763e-2,6.524845570188292e-2,-5.637610037043203e-2,1.711235764016968e-1,-4.203646902407166e-2,4.639548341453586e-3,-4.055884445496545e-3,1.258630924935448e-2,-2.776451559292779e-3,-1.023784366070774e-2,8.817108288936985e-3,-2.659664906860937e-2,7.14445034054861e-3,-1.435466025441424e-3,1.240274208149505e-3,-3.764718680837329e-3,1.028716295017727e-3,1.032012418492197e-3,-8.794183904936319e-4,2.597922814459609e-3,-6.571159498040679e-4,1.892022767235695e-4,-1.612267049238325e-4,4.76285849317595e-4,-1.204712574640791e-4],[4,17]); +end +if nargout > 2 + t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_8to8_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5 +%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:02 + +stencil_F2C = [1.324608212425595e-3,-1.220703125e-3,-1.059686569940476e-2,1.1962890625e-2,3.708902994791667e-2,-5.9814453125e-2,-7.417805989583333e-2,2.99072265625e-1,5.927225748697917e-1,2.99072265625e-1,-7.417805989583333e-2,-5.9814453125e-2,3.708902994791667e-2,1.1962890625e-2,-1.059686569940476e-2,-1.220703125e-3,1.324608212425595e-3]; +if nargout > 1 + BC_F2C = reshape([2.087365925359584,-1.227221822236823,1.090916714284468e1,-1.041312149906314,1.134214373258162,-1.269686811479259,2.075368250436277,-1.520771409696474e1,1.389750473974189,-1.48485748783569,-1.040579640776707,8.899721508624742e-1,-7.177691661032869,6.882712043721976e-1,-7.599347599660409e-1,-1.26440526850463,1.160658043364166,-5.391509178445742,6.679923135342851e-1,-7.521725463922262e-1,3.752349810676949,-2.906684049793639,2.672876913566556e1,-2.493145660873,2.782825829667436,1.443588084644871e-1,-1.307230271348211e-1,2.216804748506809,1.437719804483573e-1,-1.145830952112369e-1,3.2149320801083e-1,-2.34361439272989e-1,1.855837403850803,1.286992417665136e-1,-5.758238484341108e-2,-1.301103772185027,9.973408313121463e-1,-8.837520107531879,9.534930382604288e-1,-1.429199518808459e-2,-2.384453600787036,1.818345997558567,-1.577342441239303e1,1.422846302346762,3.154159322410927e-3,3.589377915408905e-1,-2.109524979864723e-1,9.363227187483732e-1,6.301238593470628e-2,5.732390257964316e-2,2.793582225299658,-2.033528056927806,1.621324028555783e1,-1.189222718426265,6.189671387692164e-2,5.342131492864642e-1,-4.222231296997605e-1,3.787253511209165,-3.350064226195165e-1,1.245325350855226e-2,-2.184146687275183,1.551331905163865,-1.19387312118926e1,8.349589004472998e-1,-1.35383591594438e-1,-8.014367493869704e-1,5.668415377001197e-1,-4.302842643863972,2.84183319528176e-1,3.83514596570461e-2,9.45794880458861e-1,-6.732638898386457e-1,5.234390273661413,-3.83498639699941e-1,1.440453495443018e-1,5.645065377027613e-1,-4.039220855868618e-1,3.170156487564899,-2.383074337879712e-1,1.216617375214008e-1,-1.937863224145611e-1,1.358891032584724e-1,-1.023378722711107,6.830621267493578e-2,-4.169135139842575e-3,2.270306426219975e-2,-2.38399209865252e-2,2.928388757260058e-1,-4.02060181933195e-2,6.880595419697505e-2,9.638828970166005e-2,-6.98778025332229e-2,5.609317796857676e-1,-4.429526993510467e-2,2.882554468024363e-2,3.582431391759518e-2,-2.671910644265286e-2,2.251335440088556e-1,-1.966486325944329e-2,1.791756577091828e-2,-3.350790201878844e-1,2.581411034605136e-1,-2.283079071439122,2.162064132948073e-1,-2.321676317817421e-1,6.356644993195555e-2,-4.921907141164279e-2,4.384534411523264e-1,-4.198474091936761e-2,4.597203884996652e-2,2.261662512481835e-2,-1.740429407604355e-2,1.537038474929879e-1,-1.452709057733876e-2,1.556101844926456e-2,3.097679325854209e-2,-2.394872918869654e-2,2.128879105995857e-1,-2.032077838507769e-2,2.213372706946953e-2],[5,24]); +end +if nargout > 2 + t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4 +%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:39 + +stencil_F2C = [-2.564929780505952e-4,-1.220703125e-3,2.051943824404762e-3,1.1962890625e-2,-7.181803385416667e-3,-5.9814453125e-2,1.436360677083333e-2,2.99072265625e-1,4.820454915364583e-1,2.99072265625e-1,1.436360677083333e-2,-5.9814453125e-2,-7.181803385416667e-3,1.1962890625e-2,2.051943824404762e-3,-1.220703125e-3,-2.564929780505952e-4]; +if nargout > 1 + BC_F2C = reshape([-1.673327822994457e1,1.665420585653232e1,-1.973927473454954e2,2.826257908914756e1,-6.156813484052936e1,3.974966126696054,2.880639373222355e1,-2.660797294399895e1,3.202303847857831e2,-4.598960974033412e1,1.003152507870138e2,-6.481221721577338,8.061874893005659,-7.706608975219419,9.234201483040626e1,-1.322147613066874e1,2.880209985918895e1,-1.859523138300886,-2.19528557898034e1,2.13763239391502e1,-2.476320412310993e2,3.572924955195535e1,-7.795290920161567e1,5.036097395109802,-1.287197329845376,1.244100160910233,-1.394684346537857e1,2.112320869857597,-4.60441856527933,2.978264879234002e-1,1.935199269380359,-1.88580031845933,2.330807783420723e1,-2.91749306966336,6.644427567385864,-4.302395355488564e-1,-1.148984264037686,1.109982778464842,-1.314799057941425e1,2.137017065505111,-4.084082382904703,2.606486561845403e-1,1.565462612018858,-1.508309105065004,1.772442871757061e1,-2.40060893717649,6.31422909250156,-4.052444871368164e-1,3.179380137695096e-1,-3.071741522636364e-1,3.636264382483117,-5.183996466001354e-1,2.322192602480252,-6.467931586900172e-2,6.277081831671831e-1,-6.107890751034906e-1,7.335850436798604,-1.091016872160133,3.089218691662174,6.988892111841499e-2,1.639756762532723,-1.583883992641009,1.876236237708381e1,-2.685390559486947,5.859855615166555,3.919312696253764e-3,1.611803934956753,-1.540778869185931,1.795955382224485e1,-2.495154500270907,5.003492494640851,-4.157919873785453e-2,-6.213413520816026e-1,6.242581523412814e-1,-7.823535324240331,1.222418137438074,-3.120142966701915,2.975831327541895e-1,2.12230710100161,-2.0451058741188,2.412395225004063e1,-3.423544496956294,7.325914828631551,-4.793445914188063e-1,-4.31809447483846,4.158975059021769,-4.906211424150023e1,6.975058620124194,-1.501259943091044e1,9.418983630154896e-1,-6.138797742037704e-1,5.810969368378602e-1,-6.6821474226667,9.113864351863538e-1,-1.814400062103318,1.002574935832871e-1,-9.860044240994036e-2,9.448685284232805e-2,-1.106188847377502,1.552858946803229e-1,-3.265426230113952e-1,2.062042480203908e-2,-2.745034502159738,2.655226942873441,-3.151193386341271e1,4.520932327069998,-9.883338555855938,6.423439846307804e-1,2.139968493382333,-2.067738462547628,2.450223844197562e1,-3.50699574523559,7.635005983192545,-4.923960464980415e-1,1.036406872394752,-1.002008838408551,1.188339891704641e1,-1.70302259679446,3.715789696036691,-2.407301455417931e-1,1.606428513214277,-1.552524985439685,1.840245870449254e1,-2.635131753653705,5.741508438443495,-3.708346905830083e-1,5.659562678739624e-2,-5.465656201303744e-2,6.471932864664631e-1,-9.253169145609999e-2,2.010909570941804e-1,-1.292046400706276e-2,5.435391241988039e-1,-5.252672844681129e-1,6.225561585258754,-8.91344337215917e-1,1.941632399616401,-1.253426955105431e-1,-1.552116972709218,1.499962759958308,-1.777819805127299e1,2.545472086708348,-5.545140384142844,3.580057322157566e-1],[6,24]); +end +if nargout > 2 + t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_2to2_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,35 @@ +function [IC2F,IF2C] = intOpMC_orders_2to2_ratio2to1(mc) + +mf = 2*(mc-1) + 1; + +stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0]; +stencil_width = length(stencil_F2C); +stencil_halfwidth = (stencil_width-1)/2; + +BC_F2C = [1.0./2.0,1.0./2.0]; + +Hc = speye(mc,mc); +Hc(1,1) = 1/2; +Hc(end,end) = 1/2; + +Hf = 1/2*speye(mf,mf); +Hf(1,1) = 1/4; +Hf(end,end) = 1/4; + +IF2C = sparse(mc,mf); +[BCrows, BCcols] = size(BC_F2C); +IF2C(1:BCrows, 1:BCcols) = BC_F2C; +IF2C(mc-BCrows+1:mc, mf-BCcols+1:mf) = rot90(BC_F2C,2); + +for i = BCrows+1 : mc-BCrows + IF2C(i,(2*i-stencil_halfwidth-1) :(2*i+stencil_halfwidth-1))... + = stencil_F2C; +end + +IC2F = Hf\(IF2C'*Hc); + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_4to4_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,66 @@ +% Marks New interpolation operators +% 4th order to 2nd order accurate (diagonal norm) +% M=9 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_4to4_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); +t1=[ 2047/2176 , 129/1088 , -129/2176 , 0 , 0 , 0 , 0 ,... + 0 ;... + 429/944 , 279/472 , -43/944 , 0 , 0 , 0 , 0 , 0 ;... + 111/2752 , 4913/5504 , 3/32 , -147/5504 ,... + 0 , 0 , 0 , 0 ;... + -103/784 , 549/784 , 387/784 , -1/16 , 0 ,... + 0 , 0 , 0 ;... + -335/3072 , 205/768 , 2365/3072 , 49/512 ,... + -3/128 , 0 , 0 , 0 ;... + -9/256 , 5/256 , 129/256 , 147/256 ,... + -1/16 , 0 , 0 , 0 ;... + 5/192 , -59/1024 , 43/512 , 2695/3072 ,... + 3/32 , -3/128 , 0 , 0 ;... + 23/768 , -37/768 , -43/768 , 147/256 ,... + 9/16 , -1/16 , 0 , 0 ;... + 13/2048 , -11/1024 , -43/2048 , 49/512 ,... + 55/64 , 3/32 , -3/128 , 0 ;... + -1/384 , 1/256 , 0 , -49/768 , 9/16 ,... + 9/16 , -1/16 , 0 ;... + -1/1024 , 3/2048 , 0 , -49/2048 , 3/32 ,... + 55/64 , 3/32 , -3/128]; + + +t2=[-3/128 3/32 55/64 3/32 -3/128]; +t3=[-1/16 9/16 9/16 -1/16]; +I1(1:11,1:8)=t1; +I1(M_F-10:M_F,M_C-7:M_C)=fliplr(flipud(t1)); +I1(12,5:8)=t3; +for i=13:2:M_F-12 + j=(i-1)/2; + I1(i,j-1:j+3)=t2; + I1(i+1,j:j+3)=t3; +end + +% Fine to coarse +I2=zeros(M_C,M_F); + +t1=[ 2047/4352 , 429/544 , 111/2176 , -103/544 ... + , -335/2176 , -27/544 , 5/136 , 23/544 ,... + 39/4352 , -1/272 , -3/2176 ;... + 129/7552 , 279/944 , 4913/15104 , 549/1888 ... + , 205/1888 , 15/1888 , -3/128 , -37/1888 ... + , -33/7552 , 3/1888 , 9/15104 ]; +t2=[-3/256 -1/32 3/64 9/32 55/128 9/32 3/64 -1/32 -3/256]; + +I2(1:2,1:11)=t1; + +I2(M_C-1:M_C,M_F-10:M_F)=fliplr(flipud(t1)); + +for i=3:M_C-2 + j=2*(i-3)+1; + I2(i,j:j+8)=t2; +end + +IC2F = sparse(I1); +IF2C = sparse(I2);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_6to6_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,115 @@ +% Marks New interpolation operators +% 6th order accurate (diagonal norm) +% M=19 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_6to6_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); + +t1= [6854313/6988288 , 401925/6988288 , -401925/6988288 ... + , 133975/6988288 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + 560547/1537664 , 1201479/1537664 , -240439/1537664 ... + , 16077/1537664 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + 203385/5552128 , 1225647/1388032 , 364155/2776064 ... + , -80385/1388032 , 39385/5552128 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -145919/2743808 , 721527/1371904 , 105687/171488 , ... + -25/256 , 23631/2743808 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -178863/4033024 , 1178085/8066048 , ... + 1658587/2016512 , 401925/4033024 , -15/512 , ... + 43801/8066048 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -1668147/11213056 , 4193225/11213056 , ... + 375675/2803264 , 2009625/2803264 , ... + -984625/11213056 , 3/256 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -561187/2949120 , 831521/1474560 , -788801/1474560 ... + , 412643/368640 , 39385/589824 , -43801/1474560 ... + , 5/1024 , 0 , 0 , 0 , 0 , 0 ; ... + 23/1024 , 23/147456 , -43435/221184 , ... + 26795/36864 , 39385/73728 , -43801/442368 , ... + 3/256 , 0 , 0 , 0 , 0 , 0 ; ... + 79379/368640 , -1664707/2949120 , 284431/737280 , ... + 26795/294912 , 606529/737280 , 43801/589824 , ... + -15/512 , 5/1024 , 0 , 0 , 0 , 0 ; ... + 3589/27648 , -2225/6144 , 22939/73728 , ... + -26795/221184 , 39385/73728 , 43801/73728 , ... + -25/256 , 3/256 , 0 , 0 , 0 , 0 ; ... + -720623/14745600 , 10637/92160 , -89513/1474560 , ... + -5359/147456 , 39385/589824 , 3372677/3686400 , ... + 75/1024 , -15/512 , 5/1024 , 0 , 0 , 0 ; ... + -6357/81920 , 55219/276480 , -8707/61440 , ... + 5359/368640 , -39385/442368 , 43801/73728 , ... + 75/128 , -25/256 , 3/256 , 0 , 0 , 0 ; ... + -13315/884736 , 2589/65536 , -479/16384 , ... + 5359/884736 , -7877/294912 , 43801/589824 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 , 0 , 0 ; ... + 8299/737280 , -7043/245760 , 5473/276480 , 0 , ... + 7877/737280 , -43801/442368 , 75/128 , ... + 75/128 , -25/256 , 3/256 , 0 , 0 ; ... + 11027/2949120 , -8461/884736 , 655/98304 , 0 , ... + 7877/1769472 , -43801/1474560 , 75/1024 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 , 0 ; ... + -601/614400 , 601/245760 , -601/368640 , 0 , 0 , ... + 43801/3686400 , -25/256 , 75/128 , ... + 75/128 , -25/256 , 3/256 , 0 ; ... + -601/1474560 , 601/589824 , -601/884736 , 0 , 0 , ... + 43801/8847360 , -15/512 , 75/1024 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 ] ; + + t2= [5/1024 , -15/512 , 75/1024, 231/256 , 75/1024 , -15/512 , 5/1024]; + + t3= [3/256 , -25/256 , 75/128 , 75/128 , -25/256 , 3/256]; + +I1(1:17,1:12)=t1; +I1(M_F-16:M_F,M_C-11:M_C)=fliplr(flipud(t1)); +I1(18,7:12)=t3; +for i=19:2:M_F-18 + j=(i-3)/2; + I1(i,j-1:j+5)=t2; + I1(i+1,j:j+5)=t3; +end + +% Fine to coarse +I2=zeros(M_C,M_F); + +t1=[6854313/13976576 , 2802735/3494144 , ... + 1016925/27953152 , -729595/6988288 , ... + -894315/13976576 , -1668147/6988288 , ... + -8417805/27953152 , 15525/436768 , ... + 1190685/3494144 , 89725/436768 , ... + -2161869/27953152 , -858195/6988288 , ... + -332875/13976576 , 124485/6988288 , ... + 165405/27953152 , -5409/3494144 , -9015/13976576; ... + 80385/12301312 , 1201479/3075328 , 1225647/6150656 ... + , 721527/3075328 , 1178085/24602624 , ... + 838645/6150656 , 60843/300032 , 345/6150656 , ... + -4994121/24602624 , -100125/768832 , ... + 31911/768832 , 55219/768832 , 349515/24602624 , ... + -63387/6150656 , -42305/12301312 , 5409/6150656 ... + , 9015/24602624 ; ... + -80385/5552128 , -240439/1388032 , 364155/5552128 ... + , 105687/173504 , 1658587/2776064 , 75135/694016 ... + , -2366403/5552128 , -217175/1388032 , ... + 853293/2776064 , 344085/1388032 , ... + -268539/5552128 , -78363/694016 , -64665/2776064 ... + , 5473/347008 , 29475/5552128 , -1803/1388032 , ... + -3005/5552128]; + + + + t2=[ 5/2048 , 3/512 , -15/1024 , -25/512 , ... + 75/2048 , 75/256 , 231/512 , 75/256 , ... + 75/2048 , -25/512 , -15/1024 , 3/512 , 5/2048]; + +I2(1:3,1:17)=t1; + +I2(M_C-2:M_C,M_F-16:M_F)=fliplr(flipud(t1)); + +for i=4:M_C-3 + j=2*(i-4)+1; + I2(i,j:j+12)=t2; +end +IC2F = sparse(I1); +IF2C = sparse(I2); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_8to8_ratio2to1.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,53 @@ +% Marks New interpolation operators +% 8th order accurate (diagonal norm) +% M=19 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_8to8_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); + +t1= [0.99340647972821530406e0 0.26374081087138783775e-1 -0.39561121630708175663e-1 0.26374081087138783775e-1 -0.65935202717846959438e-2 0 0 0 0 0 0 0 0 0 0 0; 0.31183959860287729600e0 0.94014160558849081601e0 -0.31646240838273622401e0 0.65141605588490816007e-1 -0.66040139712270400169e-3 0 0 0 0 0 0 0 0 0 0 0; -0.33163591467432411328e-1 0.11092588191512759415e1 -0.10539936193077965253e0 -0.77189144409925778387e-2 0.60418595406382404105e-1 -0.23395546718453703858e-1 0 0 0 0 0 0 0 0 0 0; -0.63951987538041216890e-1 0.56657207530367360435e0 0.56073157416571775151e0 -0.67107298938782711711e-1 0.54915118559238359570e-2 -0.17358748484912632117e-2 0 0 0 0 0 0 0 0 0 0; 0.22970968995770386894e0 -0.84424237330911973043e0 0.20693327723706358111e1 -0.42910115554542331920e0 -0.13191478386190563918e0 0.11675567167691342983e0 -0.10539821288804421121e-1 0 0 0 0 0 0 0 0 0; 0.10195433776615307267e0 -0.45344410547614678949e0 0.11049699103276728065e1 0.26297465812771521534e0 -0.38617448080010680341e-1 0.23925781250000000000e-1 -0.17631339153836246986e-2 0 0 0 0 0 0 0 0 0; -0.33882356049961665258e0 0.12718893449497745294e1 -0.16822328941798626751e1 0.17813590728082409984e1 0.11793036905675500906e0 -0.18266200597575250398e0 0.37689938246258754051e-1 -0.51502644057974593120e-2 0 0 0 0 0 0 0 0; -0.50400650482311136546e0 0.19901277318227816880e1 -0.29708252195230746436e1 0.23722122500767090037e1 0.24457622727717445252e0 -0.15152936311741758892e0 0.21886284536938453771e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0 0; 0.37882732714731866331e-2 0.12115606818363056270e0 -0.53924884156527999826e0 0.88886598075786086512e0 0.27660232216640139169e0 0.33730204543181760125e0 -0.12179633685076087365e0 0.38041730885639608773e-1 -0.47112422807823442564e-2 0 0 0 0 0 0 0; 0.41123781086486062886e0 -0.14418811327145423418e1 0.16793759240044487847e1 -0.57156511000645013503e0 0.24685906775203751929e0 0.76471858554416232639e0 -0.11045284035765094522e0 0.24149101163134162809e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0; 0.36463669929508579110e0 -0.13698743320477484358e1 0.18268133133896437122e1 -0.93074264690533336964e0 0.10888458847499176143e0 0.85685706622610101431e0 0.24359267370152174731e0 -0.13314605809973863070e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0 0; 0.25541802394537278164e0 -0.10497853549910180363e1 0.15921730465359157986e1 -0.96615555845660927855e0 -0.49371813550407503858e-1 0.76471858554416232639e0 0.55226420178825472608e0 -0.12074550581567081404e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0 0; -0.93469656691661424206e-1 0.26634337856674539612e0 -0.17694629839462736800e0 -0.64947940656920645005e-1 -0.54442294237495880713e-1 0.33730204543181760125e0 0.61880473767909428853e0 0.26629211619947726141e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0; -0.49445400002561969522e0 0.18168098496802059221e1 -0.23462477366129557292e1 0.11091140417405116705e1 0.98743627100815007716e-2 -0.15294371710883246528e0 0.55226420178825472608e0 0.60372752907835407022e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0; -0.24633538738293361369e0 0.93021448723433316305e0 -0.12527182680719820488e1 0.63698038058706249354e0 0.15554941210713108775e-1 -0.16865102271590880063e0 0.24359267370152174731e0 0.67646871560981190145e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0; 0.23150628239599695492e0 -0.82682792555928440531e0 0.10237569354829334077e1 -0.45129113643047460593e0 -0.10075880316409694665e-2 0.30588743421766493056e-1 -0.11045284035765094522e0 0.60372752907835407022e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0; 0.19518658723021413499e0 -0.70515100787561674409e0 0.88870680011413432419e0 -0.40458631782898657259e0 -0.19443676513391385969e-2 0.48186006490259657322e-1 -0.12179633685076087365e0 0.26629211619947726141e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0; -0.43403699494753775353e-1 0.15442805550926787092e0 -0.18997683853068679729e0 0.82584189359473172950e-1 0 -0.31213003491598462302e-2 0.22090568071530189043e-1 -0.12074550581567081404e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0; -0.58783367481931040778e-1 0.20994477957758583150e0 -0.25976152530269196017e0 0.11403438083569734975e0 0 -0.60232508112824571652e-2 0.34798953385931678187e-1 -0.13314605809973863070e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0; 0.63390647713259204145e-2 -0.22373993350504987875e-1 0.27185871992161665013e-1 -0.11561529976981026786e-1 0 0 -0.22541395991357335758e-2 0.24149101163134162809e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0; 0.10649583863025939259e-1 -0.37634916628131671889e-1 0.45812371547331598336e-1 -0.19540204529147604911e-1 0 0 -0.43498691732414597734e-2 0.38041730885639608773e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0; -0.45575492476359756866e-3 0.15951422366725914903e-2 -0.19141706840071097884e-2 0.79757111833629574515e-3 0 0 0 -0.24641939962381798784e-2 0.23925781250000000000e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0; -0.87948159845213680356e-3 0.30781855945824788125e-2 -0.36938227134989745750e-2 0.15390927972912394062e-2 0 0 0 -0.47552163607049510966e-2 0.37689938246258754051e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2;]; + + + t2= [-0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.1856585148354920384923865861096125662293072684152233190798249652677391616531107981e82 / 0.2770141785660523684734143913470637620565185365865679672509451726827594519874538713e82 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83;]; + + t3= [-0.5e1 / 0.2048e4 0.49e2 / 0.2048e4 -0.245e3 / 0.2048e4 0.1225e4 / 0.2048e4 0.1225e4 / 0.2048e4 -0.245e3 / 0.2048e4 0.49e2 / 0.2048e4 -0.5e1 / 0.2048e4;]; +I1(1:23,1:16)=t1; +I1(M_F-22:M_F,M_C-15:M_C)=fliplr(flipud(t1)); +I1(24,9:16)=t3; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DEBUGGING + +for i=25:2:M_F-24 + j=(i-3)/2; + I1(i,j-2:j+6)=t2; + I1(i+1,j-1:j+6)=t3; +end + +% Fine to coarse: A +I2=zeros(M_C,M_F); + +t1=[0.49670323986410765203e0 0.80670591743192512511e0 -0.14476656476698073533e-1 -0.19497555250083395132e0 0.16074268813765884450e0 0.22100911221277072755e0 -0.53042418939472250452e0 -0.86254139670456354517e0 0.64231825172866668299e-2 0.69727164011248914487e0 0.61825742343094006838e0 0.43307239695721032895e0 -0.15848187861199173328e0 -0.83836831742920925562e0 -0.41767239062238727391e0 0.39252899650233765024e0 0.33094736964907844809e0 -0.73592865087013788440e-1 -0.99669762780958210515e-1 0.10748160731101219580e-1 0.18056833808814782786e-1 -0.77275234787125894193e-3 -0.14911993994710636483e-2; 0.25487859584304862664e-2 0.47007080279424540800e0 0.93589176759289320118e-1 0.33386224041716468423e0 -0.11418395501435496457e0 -0.18998278818813064037e0 0.38484431284491986603e0 0.65828018436035862232e0 0.39704539050575728856e-1 -0.47252462545568042557e0 -0.44892698918501097924e0 -0.34402935194949605213e0 0.87284452472801643396e-1 0.59539401290875351190e0 0.30484430526276345964e0 -0.27096308216867871783e0 -0.23108785344796321535e0 0.50608234918774219796e-1 0.68801842319351676214e-1 -0.73322707316320135247e-2 -0.12333488857215226757e-1 0.52275043402033040521e-3 0.10087644967132781708e-2; -0.22656973277393401539e-1 -0.93771184228725468159e0 -0.52699680965389826267e-1 0.19581430555012001670e1 0.16586148101136731602e1 0.27435837109255836264e1 -0.30164708462284474624e1 -0.58235016129647971089e1 -0.10472767830023645070e1 0.32615209891555982371e1 0.35478595826728208891e1 0.30921640208240511054e1 -0.34364793368678654581e0 -0.45566547247355317663e1 -0.24329078834671892590e1 0.19882413967858906122e1 0.17259601262271516763e1 -0.36895458453625989435e0 -0.50448363278283993228e0 0.52797763052066775846e-1 0.88972343374038343284e-1 -0.37175165926095403240e-2 -0.71737841052106668688e-2; 0.21626748627943672418e-2 0.27636709246281031633e-1 -0.55259484658035070394e-3 -0.33553649469391355855e-1 -0.49244245327794891233e-1 0.93489376222103426899e-1 0.45734620580442859300e0 0.66579609152388478842e0 0.24716623315221892947e0 -0.15893464065423852815e0 -0.25881084331023040874e0 -0.26865808253702445366e0 -0.18060020510041282529e-1 0.30841043055726240020e0 0.17712461121229459766e0 -0.12549015561536110372e0 -0.11250298507032009025e0 0.22964117700293735857e-1 0.31709446610808019222e-1 -0.32149051440245356510e-2 -0.54335286230388551025e-2 0.22177994574852164638e-3 0.42797426992744435493e-3;]; + + + + t2=[-0.23556211403911721282e-2 -0.12207031250000000000e-2 0.18844969123129377026e-1 0.11962890625000000000e-1 -0.65957391930952819589e-1 -0.59814453125000000000e-1 0.13191478386190563918e0 0.29907226562500000000e0 0.33510652017261795103e0 0.29907226562500000000e0 0.13191478386190563918e0 -0.59814453125000000000e-1 -0.65957391930952819589e-1 0.11962890625000000000e-1 0.18844969123129377026e-1 -0.12207031250000000000e-2 -0.23556211403911721282e-2;]; + + +I2(1:4,1:23)=t1; + +I2(M_C-3:M_C,M_F-22:M_F)=fliplr(flipud(t1)); + +for i=5:M_C-4 + j=2*(i-5)+1; + I2(i,j:j+16)=t2; +end + +IC2F = sparse(I1); +IF2C = sparse(I2); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpAWW.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,68 @@ +classdef InterpAWW < sbp.InterpOps + properties + + % Interpolation operators + IC2F + IF2C + + % Orders used on coarse and fine sides + order_C + order_F + + % Grid points, refinement ratio. + ratio + m_C + m_F + + % Boundary accuracy of IC2F and IF2C. + acc_C2F + acc_F2C + end + + methods + % accOp : String, 'C2F' or 'F2C'. Specifies which of the operators + % should have higher accuracy. + function obj = InterpAWW(m_C,m_F,order_C,order_F,accOp) + assertIsMember(accOp, {'C2F','F2C'}); + + ratio = (m_F-1)/(m_C-1); + h_C = 1; + + assert(order_C == order_F,... + 'Error: Different orders of accuracy not available'); + + switch ratio + case 2 + switch order_C + case 2 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_2to2_ratio2to1(m_C, h_C, accOp); + case 4 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_4to4_ratio2to1(m_C, h_C, accOp); + case 6 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_6to6_ratio2to1(m_C, h_C, accOp); + case 8 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_8to8_ratio2to1(m_C, h_C, accOp); + otherwise + error(['Order ' num2str(order_C) ' not available.']); + end + otherwise + error(['Grid ratio ' num2str(ratio) ' not available']); + end + + obj.IC2F = IC2F; + obj.IF2C = IF2C; + obj.order_C = order_C; + obj.order_F = order_F; + obj.ratio = ratio; + obj.m_C = m_C; + obj.m_F = m_F; + + end + + function str = string(obj) + str = [class(obj) '_orders' num2str(obj.order_F) 'to'... + num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1']; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpMC.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,61 @@ +classdef InterpMC < sbp.InterpOps + properties + + % Interpolation operators + IC2F + IF2C + + % Orders used on coarse and fine sides + order_C + order_F + + % Grid points, refinement ratio. + ratio + m_C + m_F + end + + methods + function obj = InterpMC(m_C,m_F,order_C,order_F) + + ratio = (m_F-1)/(m_C-1); + + assert(order_C == order_F,... + 'Error: Different orders of accuracy not available'); + + switch ratio + case 2 + switch order_C + case 2 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_2to2_ratio2to1(m_C); + case 4 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_4to4_ratio2to1(m_C); + case 6 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_6to6_ratio2to1(m_C); + case 8 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_8to8_ratio2to1(m_C); + otherwise + error(['Order ' num2str(order_C) ' not available.']); + end + otherwise + error(['Grid ratio ' num2str(ratio) ' not available']); + end + + obj.IC2F = IC2F; + obj.IF2C = IF2C; + obj.order_C = order_C; + obj.order_F = order_F; + obj.ratio = ratio; + obj.m_C = m_C; + obj.m_F = m_F; + + + end + + function str = string(obj) + str = [class(obj) '_orders' num2str(obj.order_F) 'to'... + num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1']; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpOps.m Thu Jul 26 10:53:05 2018 -0700 @@ -0,0 +1,14 @@ +classdef (Abstract) InterpOps + properties (Abstract) + % C and F may refer to coarse and fine, but it's not a must. + IC2F % Interpolation operator from "C" to "F" + IF2C % Interpolation operator from "F" to "C" + + end + + methods (Abstract) + % Returns a string representation of the type of operator. + str = string(obj) + end + +end \ No newline at end of file
--- a/+scheme/bcSetup.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+scheme/bcSetup.m Thu Jul 26 10:53:05 2018 -0700 @@ -3,46 +3,108 @@ % Each bc is a struct with the fields % * type -- Type of boundary condition % * boundary -- Boundary identifier -% * data -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem +% * 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 -function [closure, S] = bcSetup(diffOp, bc, S_sign) +% 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, bcs, S_sign) default_arg('S_sign', 1); - assertType(bc, 'cell'); + assertType(bcs, 'cell'); assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); + verifyBcFormat(bcs, diffOp); + % Setup storage arrays closure = spzeros(size(diffOp)); - penalties = {}; - dataFunctions = {}; - dataParams = {}; + gridData = {}; + symbolicData = {}; - for i = 1:length(bc) - assertType(bc{i}, 'struct'); - [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); + % Collect closures, penalties and data + for i = 1:length(bcs) + [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type); closure = closure + localClosure; - if isempty(bc{i}.data) + [ok, isSym, data] = parseData(bcs{i}, penalty, diffOp.grid); + + if ~ok + % There was no data continue end - assertType(bc{i}.data, 'function_handle'); - coord = diffOp.grid.getBoundary(bc{i}.boundary); - assertNumberOfArguments(bc{i}.data, 1+size(coord,2)); - - penalties{end+1} = penalty; - dataFunctions{end+1} = bc{i}.data; - dataParams{end+1} = num2cell(coord ,1); + if isSym + symbolicData{end+1} = data; + else + gridData{end+1} = data; + end end + % Setup penalty function O = spzeros(size(diffOp),1); function v = S_fun(t) v = O; - for i = 1:length(dataFunctions) - v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:}); + for i = 1:length(gridData) + v = v + gridData{i}.penalty*gridData{i}.func(t); + end + + for i = 1:length(symbolicData) + v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:}); end v = S_sign * v; end S = @S_fun; end + +function verifyBcFormat(bcs, diffOp) + for i = 1:length(bcs) + assertType(bcs{i}, 'struct'); + assertStructFields(bcs{i}, {'type', 'boundary'}); + + if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data) + continue + end + + if ~isa(bcs{i}.data, 'function_handle') + error('bcs{%d}.data should be a function of time or a function of time and space',i); + end + + b = diffOp.grid.getBoundary(bcs{i}.boundary); + + dim = size(b,2); + + if nargin(bcs{i}.data) == 1 + % Grid data (only function of time) + assertSize(bcs{i}.data(0), 1, size(b)); + elseif nargin(bcs{i}.data) ~= 1+dim + error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bcs{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + end + end +end + +function [ok, isSym, dataStruct] = parseData(bc, penalty, grid) + if ~isfield(bc,'data') || isempty(bc.data) + ok = false; + return + end + ok = true; + + nArg = nargin(bc.data); + + if nArg > 1 + % Symbolic data + isSym = true; + coord = grid.getBoundary(bc.boundary); + dataStruct.penalty = penalty; + dataStruct.func = bc.data; + dataStruct.coords = num2cell(coord, 1); + else + % Grid data + isSym = false; + dataStruct.penalty = penalty; + dataStruct.func = bcs{i}.data; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeScaled.m Thu Jul 26 10:53:05 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
--- a/+time/SBPInTimeSecondOrderFormImplicit.m Wed Jul 25 18:53:07 2018 -0700 +++ b/+time/SBPInTimeSecondOrderFormImplicit.m Thu Jul 26 10:53:05 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)
--- a/.hgtags Wed Jul 25 18:53:07 2018 -0700 +++ b/.hgtags Thu Jul 26 10:53:05 2018 -0700 @@ -1,2 +1,3 @@ 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1 +0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0
--- a/Color.m Wed Jul 25 18:53:07 2018 -0700 +++ b/Color.m Thu Jul 26 10:53:05 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
--- a/TextTable.m Wed Jul 25 18:53:07 2018 -0700 +++ b/TextTable.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertIsMember.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertSize.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertStructFields.m Thu Jul 26 10:53:05 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
--- a/assert_size.m Wed Jul 25 18:53:07 2018 -0700 +++ b/assert_size.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/convergencePlot.m Thu Jul 26 10:53:05 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
--- a/gaussian.m Wed Jul 25 18:53:07 2018 -0700 +++ b/gaussian.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mononomial.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nextColor.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pointIndex.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rickerWavelet.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stencilEquation.m Thu Jul 26 10:53:05 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/vandermonde.m Thu Jul 26 10:53:05 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