Mercurial > repos > public > sbplib
changeset 717:8e4274ee6dd8 feature/utux2D
Merge with feature/poroelastic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sat, 03 Mar 2018 14:58:21 -0800 |
| parents | 2d85f17a8aec (diff) 60eb7f46d8d9 (current diff) |
| children | 71aa5828cbbf |
| files | diffSymfun.m |
| diffstat | 26 files changed, 1451 insertions(+), 127 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Circle.m Sat Mar 03 14:58:21 2018 -0800 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Rectangle.m Sat Mar 03 14:58:21 2018 -0800 @@ -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
--- a/+multiblock/Def.m Thu Feb 15 16:38:36 2018 -0800 +++ /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 - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DefCurvilinear.m Sat Mar 03 14:58:21 2018 -0800 @@ -0,0 +1,101 @@ +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 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Definition.m Sat Mar 03 14:58:21 2018 -0800 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:56:32 + +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.288511407918266e-1,-6.442557039591332e-2,6.855399328409322e-1,1.572300335795339e-1,-2.438404903851266e-1,7.469202451925633e-1,-8.434338091984989e-2,2.921716904599249e-1,1.379279767221774e-2,-1.318963988361088e-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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:58:05 + +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,-2.991278410589404e-33]; +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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:52:19 + +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([6.641823072821646e-1,-9.461353301006098e-2,6.490928427434416e-2,6.933503739604814e-1,3.504421573787056e-1,-6.890799169004223e-2,-1.977985424802404e-1,5.011015965140368e-1,-1.405177377247462e-1,-2.597998113910025e-1,3.313888743609167e-1,2.533930978221618e-1,7.441130208605595e-2,-8.922642832077801e-2,7.452614450107663e-1,-2.038477449219133e-2,-8.90962147907633e-3,2.928129961309942e-1,2.178962135739788e-2,2.094432837031309e-2,-1.443651496959125e-1,2.813056939226545e-2,-1.144388744639026e-2,-3.684372837980203e-2,1.788671569293071e-2,-1.216141243321431e-2,4.26819573669726e-2,-3.285698393466837e-2,1.840487209794448e-2,-1.153648202068284e-2,1.108922252680633e-2,-5.926946032396911e-3,3.11230890594676e-3],[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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:53:51 + +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,-4.646656476696072e-17,1.330722210499341e-17,1.77e2./2.72e2,3.0./8.0,-5.9e1./6.88e2,1.125919117647019e-2,3.546742584745763e-1,1.335392441860465e-2,-4.9e1./5.44e2,2.335805084745763e-1,3.204941860465116e-1,-1.194852941176471e-2,1.350635593220374e-2,5.308866279069767e-1,-9.0./5.44e2,-1.11228813559322e-2,2.943313953488372e-1,-2.803308823529412e-2,2.105402542372881e-2,-1.580668604651122e-2,-9.0./5.44e2,1.430084745762712e-2,-5.450581395348837e-2,-9.191176470588592e-4,7.944915254237005e-4,-5.450581395348837e-3,1.838235294117539e-3,-1.588983050847458e-3,2.180232558139535e-3,2.297794117649301e-4,-1.986228813559669e-4,2.725290697674479e-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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:48:32 + +stencil_F2C = [-5.957031249998899e-3,3.0./5.12e2,3.574218749999339e-2,-2.5e1./5.12e2,-8.935546874998348e-2,7.5e1./2.56e2,6.19140624999978e-1,7.5e1./2.56e2,-8.935546874998348e-2,-2.5e1./5.12e2,3.574218749999339e-2,3.0./5.12e2,-5.957031249998899e-3]; +if nargout > 1 + BC_F2C = reshape([5.233890618131316e-1,-1.594459192645162e-2,3.532688727636746e-2,8.021234957689218e-1,3.906832051735622e-1,-1.732229516322391e-1,-8.876626864832848e-2,2.900912357966253e-1,-1.600356115708854e-1,-1.04402537502747e-1,2.346179009198368e-1,6.091329306528954e-1,1.56127552270282e-1,-1.168382445709634e-1,1.040987887887253,-2.387061036980731e-1,1.363504965974362e-1,1.082611928255256e-1,-5.745310654053959e-1,3.977694249198524e-1,-9.376217911401942e-1,3.554518646054659e-2,5.609157787390911e-5,-1.564625311232018e-1,6.107907974027045e-1,-3.786608696698134e-1,7.02265869951073e-1,2.054294270642539e-1,-1.302300112378258e-1,2.47894140769089e-1,-2.657085326191222e-1,1.568761445933419e-1,-2.632906518005064e-1,-1.228047556139637e-1,7.182193248980263e-2,-1.129123824234603e-1,5.81125878040401e-2,-3.466364400804729e-2,5.683680203337197e-2,1.781337575097064e-2,-1.030572999042705e-2,1.577197067502767e-2,-1.4438021157683e-2,8.391316308372733e-3,-1.295341173690281e-2,-1.54801862773853e-3,8.794183904937343e-4,-1.2989614072299e-3,1.573818938200423e-3,-8.940753636683842e-4,1.320610764016837e-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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:49:13 + +stencil_F2C = [1.074218750000866e-3,3.0./5.12e2,-6.445312500005198e-3,-2.5e1./5.12e2,1.6113281250013e-2,7.5e1./2.56e2,4.785156249999827e-1,7.5e1./2.56e2,1.6113281250013e-2,-2.5e1./5.12e2,-6.445312500005198e-3,3.0./5.12e2,1.074218750000866e-3]; +if nargout > 1 + BC_F2C = reshape([1.0./2.0,1.257360194105895e-16,-6.179984855389442e-17,3.122529826246464e-18,6.876076086160102e-1,1.5e1./3.2e1,-3.461879841386957e-1,3.502577439820879e-2,3.099721991997423e-3,2.228547003766734e-1,9.36365299936193e-3,-3.157910466040993e-3,-1.057891432064606e-1,2.355631659452257e-1,6.070385985337517e-1,-4.847496617839159e-2,-4.809232246873087e-3,5.154694068408974e-3,7.049223649022339e-1,1.01674934980955e-2,3.460117842882258e-2,-4.99662149858484e-2,5.210650763094795e-1,2.233592875303227e-1,-2.239024779745777e-3,4.254668119952138e-4,9.213872590841121e-3,3.910524351558031e-1,-9.048711215107312e-2,8.597375629526316e-2,-3.468219752858717e-1,3.250682992395968e-1,-1.309107466664176e-1,1.199249722305999e-1,-4.071552384959323e-1,1.474417123938745e-1,-3.028638773902806e-2,3.046016554982085e-2,-1.081315128181481e-1,1.120705238850529e-2,7.977722870035704e-2,-6.772545887787698e-2,2.002704188375896e-1,-5.427183680490626e-2,6.524845570188395e-2,-5.637610037043228e-2,1.711235764016971e-1,-4.20364690240717e-2,4.639548341457899e-3,-4.055884445499834e-3,1.258630924936465e-2,-2.776451559295031e-3,-1.023784366070731e-2,8.817108288936754e-3,-2.659664906860907e-2,7.14445034054856e-3,-1.435466025443095e-3,1.240274208150634e-3,-3.764718680840621e-3,1.02871629501858e-3,1.032012418491475e-3,-8.794183904932091e-4,2.597922814458923e-3,-6.571159498039887e-4,1.892022767233693e-4,-1.612267049239958e-4,4.762858493182703e-4,-1.204712574642361e-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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:35:35 + +stencil_F2C = [1.324608212162265e-3,-1.220703125e-3,-1.059686569729812e-2,1.1962890625e-2,3.708902994054343e-2,-5.9814453125e-2,-7.417805988108686e-2,2.99072265625e-1,5.927225748513586e-1,2.99072265625e-1,-7.417805988108686e-2,-5.9814453125e-2,3.708902994054343e-2,1.1962890625e-2,-1.059686569729812e-2,-1.220703125e-3,1.324608212162265e-3]; +if nargout > 1 + BC_F2C = reshape([4.992887925847094e-1,5.498475469512003e-4,-4.887770641026507e-3,4.665521105138738e-4,-5.08176255953527e-4,9.371526406387121e-1,3.692200393998361e-1,-4.121716266360099e-2,-5.793634963116702e-2,9.198688438966184e-2,-1.944260863292151e-2,1.005123109345939e-1,-1.599307393413331e-1,1.840524510920726e-2,-3.030570567820816e-2,-3.679491052410095e-1,4.675912786945224e-1,7.693828786649279e-1,7.991701233784946e-2,-1.116312262557717e-1,2.200741834857246e-2,-2.268766613354886e-2,1.092003215561338,-4.60408097516935e-2,1.1739885344787e-1,-2.838913563581289e-1,2.003650598970398e-1,-7.263431119606933e-1,4.24704116325088e-1,-4.205790124298841e-1,-3.073510269746697e-1,2.518096165802759e-1,-2.465892549251015,5.412211006532582e-1,-5.069079992942043e-1,-1.031849105852084e-2,-5.890344569485665e-4,3.339752549138614e-2,1.067378101911543e-1,9.080077690027688e-1,1.009716377260241,-8.057493032708039e-1,7.552997687263325,-8.037294231858354e-1,2.428377173577377,5.320914256727015e-5,6.650778016823549e-2,-1.530110402831001,2.984406952457918e-1,-1.991084960210322e-1,-7.06914823978699e-1,6.727706602793345e-1,-7.843914264173437,1.107103572717113,-2.43929976717732,5.604791245128845e-1,-4.425298387359292e-1,3.967766340253329,-3.522369040702429e-1,3.122097817337138e-2,-4.049040744944211e-1,1.757666832053552e-1,2.891078390407337e-1,-3.322243688079751e-1,1.135931681420319,-1.441917590296272e-1,5.871333497354273e-2,2.140714079363586e-1,-1.469694367395001e-1,5.07970154947088e-1,1.015274187485206e-1,-2.054517120288494e-2,-5.67835102181511e-1,1.703408698715292e-1,-4.592057814970091e-1,2.273402578242697e-1,-1.432526269706648e-1,8.529823291916478e-1,-1.712633233737925e-2,-1.192523672270192e-1,-1.148999373309336e-1,7.490058655859051e-2,-4.81232298672226e-1,1.655673895947216e-2,5.219724538051045e-2,4.249355033834381e-3,-9.573018626796837e-3,1.660158192640221e-1,-2.810038374360505e-2,5.562029743544568e-2,2.839191358593554e-2,-1.730855226292001e-2,9.362690670625201e-2,3.103571996021832e-4,-1.975964029518344e-2,-1.554790540676067e-2,1.299770214218259e-2,-1.279218656827878e-1,1.403531594655547e-2,-1.878922540621182e-2,-2.813561320718192e-2,2.083750432170367e-2,-1.736115084499458e-1,1.485151842449431e-2,-1.284855833421381e-2,4.465740308545228e-3,-3.527224685683412e-3,3.228400709661693e-2,-3.214671576338129e-3,3.743039692474543e-3,3.800396827092577e-3,-2.857121549390628e-3,2.438937990996417e-2,-2.183643594689368e-3,2.116332353479732e-3,5.074056896907446e-3,-3.922846806579957e-3,3.487143946867525e-2,-3.328581653207827e-3,3.625546051653038e-3],[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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 7.1. +% 11-Jul-2017 14:44:32 + +stencil_F2C = [-2.564929780864058e-4,-1.220703125e-3,2.051943824691247e-3,1.1962890625e-2,-7.181803386419364e-3,-5.9814453125e-2,1.436360677283873e-2,2.99072265625e-1,4.820454915339516e-1,2.99072265625e-1,1.436360677283873e-2,-5.9814453125e-2,-7.181803386419364e-3,1.1962890625e-2,2.051943824691247e-3,-1.220703125e-3,-2.564929780864058e-4]; +if nargout > 1 + BC_F2C = reshape([3.84889400178856e-1,1.112426550601988e-1,-1.318495369488137,1.887814023562389e-1,-4.11247635928463e-1,2.655099795914844e-2,1.028610922995518,2.364211968142455e-1,2.059181912003073,-4.340414328006023e-1,1.075497633269875,-7.409800036086324e-2,7.37838637890288e-2,1.306633125422226e-2,8.451322719188257e-1,-1.210055485680529e-1,2.636024797025984e-1,-1.701872129868639e-2,-4.635230469768725e-1,6.090754537930958e-1,-1.489759296753512,4.8674037213944e-1,-1.179464341635484,7.944500441079773e-2,-7.260695877364954e-2,7.032243845450474e-2,-3.47295282501933e-2,1.20392464519063e-1,-2.651360013733321e-1,1.767343065829761e-2,-2.096982321858537e-2,4.637607722376962e-3,9.018016569036531e-1,2.906245591263948e-1,-3.442417302951498e-1,2.096340897620274e-2,-2.874137204380543e-3,2.384228040892216e-3,-2.026040077382159e-2,2.573962718116355e-1,1.054557409115141e-2,-3.708993639747955e-3,1.992191138619206e-2,-1.470164138115057e-2,2.155634694638224e-2,1.340780622608816e-1,7.925832983348031e-1,-4.875561424378325e-2,-2.350572186645021e-4,3.076716550868981e-4,-8.141271258942936e-3,3.404241369471586e-3,1.185477795471203,8.70934314005432e-3,-1.262976611986137e-2,8.032247391372543e-3,1.316324784899156e-3,-4.086253158542394e-2,8.015278252366783e-1,2.175869750865401e-1,4.007238901530734e-4,3.864273806524533e-4,-1.508600083692204e-2,3.153609486577761e-3,3.042286208824603e-3,3.82047294861268e-1,4.397753338251498e-2,-2.563454837101723e-2,1.41743994040605e-3,7.608105232124472e-2,-5.977718736752081e-1,3.20050013977378e-1,6.802094242804826e-2,-4.194020001568932e-2,7.25315420916367e-2,9.186386275429062e-2,-6.573062309340446e-1,1.385772992079256e-1,1.775827558199515e-2,-1.127406881785768e-2,1.811068993199809e-2,2.791591769899442e-2,-1.928604231232766e-1,6.083089919559692e-3,-3.221756216941285e-2,1.711212226954443e-2,2.901092273409817e-2,-5.378021102520328e-2,2.992550025960758e-1,-4.666662872498935e-2,-3.465064705533973e-2,2.133105869707025e-2,-4.756361192707918e-2,-3.854942718050275e-2,2.54971564582897e-1,-3.334541871752157e-2,-1.305394187273881e-3,4.611114655140197e-4,8.244321726072152e-3,-4.278000168944293e-3,2.105663050394826e-2,-1.821301944614694e-3,5.395937659538634e-3,-2.783571877605622e-3,-8.060157292374778e-3,1.022596015716632e-2,-5.706697079870192e-2,7.93956741906313e-3,1.155407428597162e-4,2.111811854597317e-4,-7.982283350430192e-3,2.363443780394039e-3,-9.897869059620054e-3,1.174758160329132e-3,-1.714734530349539e-3,1.229820838718235e-3,-7.4133437108253e-3,-5.029932466958482e-4,6.964815291000956e-3,-1.280740277254815e-3,-2.947086158334019e-4,2.091452223113523e-4,-1.20812830748558e-3,-1.054095401829176e-4,1.280105034963581e-3,-2.33676862216909e-4,-4.642133749083303e-4,4.859542003098669e-4,-6.379322438624295e-3,1.046455654592667e-3,-2.762766278515372e-3,2.407540045947593e-4,-1.891185255825058e-4,1.906094692056257e-4,-2.389371525102409e-3,3.700689438675072e-4,-9.076857917948225e-4,7.171015964436335e-5,8.001159359856251e-4,-7.732304515235331e-4,9.164656932421384e-3,-1.312190264641192e-3,2.858518569557242e-3,-1.845518711255185e-4],[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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -0,0 +1,67 @@ +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) + + 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 Sat Mar 03 14:58:21 2018 -0800 @@ -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 Sat Mar 03 14:58:21 2018 -0800 @@ -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/Utux.m Thu Feb 15 16:38:36 2018 -0800 +++ b/+scheme/Utux.m Sat Mar 03 14:58:21 2018 -0800 @@ -2,7 +2,7 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing - x % Grid + grid % Grid order % Order accuracy for the approximation H % Discrete norm @@ -17,41 +17,40 @@ methods - function obj = Utux(m,xlim,order,operator) - default_arg('a',1); - - %Old operators - % [x, h] = util.get_grid(xlim{:},m); - %ops = sbp.Ordinary(m,h,order); - - + function obj = Utux(g ,order, operator) + default_arg('operator','Standard'); + + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + switch operator - case 'NonEquidistant' - ops = sbp.D1Nonequidistant(m,xlim,order); - obj.D1 = ops.D1; +% case 'NonEquidistant' +% ops = sbp.D1Nonequidistant(m,xlim,order); +% obj.D1 = ops.D1; case 'Standard' ops = sbp.D2Standard(m,xlim,order); obj.D1 = ops.D1; - case 'Upwind' - ops = sbp.D1Upwind(m,xlim,order); - obj.D1 = ops.Dm; +% case 'Upwind' +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; otherwise error('Unvalid operator') end - obj.x=ops.x; + + obj.grid = g; - obj.H = ops.H; obj.Hi = ops.HI; obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.D=obj.D1; + obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -61,17 +60,27 @@ % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj,boundary,type,data) - default_arg('type','neumann'); - default_arg('data',0); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; + penalty = -obj.Hi*tau; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + end function N = size(obj) @@ -81,9 +90,9 @@ end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Utux2D.m Sat Mar 03 14:58:21 2018 -0800 @@ -0,0 +1,277 @@ +classdef Utux2D < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + grid % Grid + order % Order accuracy for the approximation + v0 % Initial data + + a % Wave speed a = [a1, a2]; + + H % Discrete norm + H_x, H_y % Norms in the x and y directions + Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms + + % Derivatives + Dx, Dy + + % Boundary operators + e_w, e_e, e_s, e_n + + D % Total discrete operator + + % String, type of interface coupling + % Default: 'upwind' + % Other: 'centered' + coupling_type + + % String, type of interpolation operators + % Default: 'AWW' (Almquist Wang Werpers) + % Other: 'MC' (Mattsson Carpenter) + interpolation_type + + + % Cell array, damping on upwstream and downstream sides. + interpolation_damping + + end + + + methods + function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) + + default_arg('interpolation_damping',{0,0}); + default_arg('interpolation_type','AWW'); + default_arg('coupling_type','upwind'); + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + assert(isa(g, 'grid.Cartesian')) + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + obj.grid = g; + + % Operator sets + ops_x = opSet(m_x, xlim, order); + ops_y = opSet(m_y, ylim, order); + Ix = speye(m_x); + Iy = speye(m_y); + + % Norms + Hx = ops_x.H; + Hy = ops_y.H; + Hxi = ops_x.HI; + Hyi = ops_y.HI; + + obj.H_x = Hx; + obj.H_y = Hy; + obj.H = kron(Hx,Hy); + obj.Hi = kron(Hxi,Hyi); + obj.Hx = kron(Hx,Iy); + obj.Hy = kron(Ix,Hy); + obj.Hxi = kron(Hxi,Iy); + obj.Hyi = kron(Ix,Hyi); + + % Derivatives + Dx = ops_x.D1; + Dy = ops_y.D1; + obj.Dx = kron(Dx,Iy); + obj.Dy = kron(Ix,Dy); + + % Boundary operators + obj.e_w = kr(ops_x.e_l, Iy); + obj.e_e = kr(ops_x.e_r, Iy); + obj.e_s = kr(Ix, ops_y.e_l); + obj.e_n = kr(Ix, ops_y.e_r); + + obj.m = m; + obj.h = [ops_x.h ops_y.h]; + obj.order = order; + obj.a = a; + obj.coupling_type = coupling_type; + obj.interpolation_type = interpolation_type; + obj.interpolation_damping = interpolation_damping; + obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); + + end + % Closure functions return the opertors applied to the own domain to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); + + sigma = -1; % Scalar penalty parameter + switch boundary + case {'w','W','west','West'} + tau = sigma*obj.a(1)*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + + case {'s','S','south','South'} + tau = sigma*obj.a(2)*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + end + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + m_neighbour = neighbour_scheme.m(2); + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + m_neighbour = neighbour_scheme.m(2); + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + m_neighbour = neighbour_scheme.m(1); + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + m_neighbour = neighbour_scheme.m(1); + end + + switch obj.coupling_type + + % Upwind coupling (energy dissipation) + case 'upwind' + sigma_ds = -1; %"Downstream" penalty + sigma_us = 0; %"Upstream" penalty + + % Energy-preserving coupling (no energy dissipation) + case 'centered' + sigma_ds = -1/2; %"Downstream" penalty + sigma_us = 1/2; %"Upstream" penalty + + otherwise + error(['Interface coupling type ' coupling_type ' is not available.']) + end + + % Check grid ratio for interpolation + switch boundary + case {'w','W','west','West','e','E','east','East'} + m = obj.m(2); + case {'s','S','south','South','n','N','north','North'} + m = obj.m(1); + end + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_us = interpOpSet.IF2C; + I_neighbour2local_ds = interpOpSet.IF2C; + I_local2neighbour_us = interpOpSet.IC2F; + I_local2neighbour_ds = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_neighbour2local_us = interpOpSet.IC2F; + I_neighbour2local_ds = interpOpSet.IC2F; + I_local2neighbour_us = interpOpSet.IF2C; + I_local2neighbour_ds = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_us = interpOpSetC2F.IF2C; + I_neighbour2local_ds = interpOpSetF2C.IF2C; + I_local2neighbour_us = interpOpSetC2F.IC2F; + I_local2neighbour_ds = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_us = interpOpSetF2C.IC2F; + I_neighbour2local_ds = interpOpSetC2F.IC2F; + I_local2neighbour_us = interpOpSetF2C.IF2C; + I_local2neighbour_ds = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_us = speye(m,m); + I_neighbour2local_ds = speye(m,m); + end + + int_damp_us = obj.interpolation_damping{1}; + int_damp_ds = obj.interpolation_damping{2}; + + I = speye(m,m); + I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us; + I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds; + + + switch boundary + case {'w','W','west','West'} + tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a(1)... + *obj.e_w*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; + case {'e','E','east','East'} + tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; + closure = obj.Hi*tau*obj.e_e'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a(1)... + *obj.e_e*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; + case {'s','S','south','South'} + tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a(2)... + *obj.e_s*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; + case {'n','N','north','North'} + tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; + closure = obj.Hi*tau*obj.e_n'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a(2)... + *obj.e_n*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; + end + + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file