Mercurial > repos > public > sbplib
changeset 968:a4ad90b37998 feature/poroelastic
Merge with default.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 23 Dec 2018 14:39:31 +0100 |
| parents | 368a2773f78b (current diff) 2412f407749a (diff) |
| children | adae8063ea2f |
| files | +multiblock/DiffOp.m +sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m +sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m +sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m +sbp/InterpAWW.m +sbp/InterpMC.m +scheme/Elastic2dVariable.m +scheme/Wave.m .hgtags |
| diffstat | 71 files changed, 2022 insertions(+), 1059 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Nodes.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,47 @@ +classdef Nodes < grid.Grid + properties + coords + end + + methods + % Creates a grid with one point for each row in coords. + % The dimension equals the number of columns in coords. + function obj = Nodes(coords) + obj.coords = coords; + end + + function o = N(obj) + o = size(obj.coords, 1); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = size(obj.coords, 2); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = obj.coords; + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + error('Not implemented'); + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + error('Not implemented'); + end + + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. + function bs = getBoundaryNames(obj) + error('Not implemented'); + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + error('Not implemented'); + end + end +end
--- a/+multiblock/+domain/Circle.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+multiblock/+domain/Circle.m Sun Dec 23 14:39:31 2018 +0100 @@ -65,10 +65,10 @@ 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.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);
--- a/+multiblock/DiffOp.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+multiblock/DiffOp.m Sun Dec 23 14:39:31 2018 +0100 @@ -10,7 +10,7 @@ end methods - function obj = DiffOp(doHand, g, order, doParam) + function obj = DiffOp(doHand, g, order, doParam, intfTypes) % 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, ...) @@ -24,14 +24,17 @@ % corresponding function handle as extra parameters: % doHand(..., doParam{i}{:}) Otherwise doParam is sent as % extra parameters to all doHand: doHand(..., doParam{:}) + % + % intfTypes (optional) -- nBlocks x nBlocks cell array of types for + % every interface. default_arg('doParam', []) + default_arg('intfTypes', cell(g.nBlocks(), g.nBlocks()) ); [getHand, getParam] = parseInput(doHand, g, doParam); + obj.order = order; nBlocks = g.nBlocks(); - obj.order = order; - % Create the diffOps for each block obj.diffOps = cell(1, nBlocks); for i = 1:nBlocks @@ -70,12 +73,11 @@ continue end - - [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}, intfTypes{i,j}); D{i,i} = D{i,i} + ii; D{i,j} = D{i,j} + ij; - [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}, intfTypes{i,j}); D{j,j} = D{j,j} + jj; D{j,i} = D{j,i} + ji; end @@ -269,33 +271,8 @@ [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); % Expand to matrix for full domain. - div = obj.blockmatrixDiv; - if ~iscell(blockClosure) - temp = blockmatrix.zero(div); - temp{I,I} = blockClosure; - closure = blockmatrix.toMatrix(temp); - else - for i = 1:length(blockClosure) - temp = blockmatrix.zero(div); - temp{I,I} = blockClosure{i}; - closure{i} = blockmatrix.toMatrix(temp); - end - end - - if ~iscell(blockPenalty) - div{2} = size(blockPenalty, 2); % Penalty is a column vector - p = blockmatrix.zero(div); - p{I} = blockPenalty; - penalty = blockmatrix.toMatrix(p); - else - % TODO: used by beam equation, should be eliminated. SHould only set one BC per call - for i = 1:length(blockPenalty) - div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector - p = blockmatrix.zero(div); - p{I} = blockPenalty{i}; - penalty{i} = blockmatrix.toMatrix(p); - end - end + closure = multiblock.local2globalClosure(blockClosure, obj.blockmatrixDiv, I); + penalty = multiblock.local2globalPenalty(blockPenalty, obj.blockmatrixDiv, I); end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/local2globalClosure.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,10 @@ +% Takes the block-local closures and turns it into a global closure +% local -- The local closure +% div -- block matrix division for the diffOp +% I -- Index of blockmatrix block +function closure = local2globalClosure(local, div, I) + closure_bm = blockmatrix.zero(div); + closure_bm{I,I} = local; + + closure = blockmatrix.toMatrix(closure_bm); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/local2globalPenalty.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,11 @@ +% Takes the block-local penalty and turns it into a global penalty +% local -- The local penalty +% div -- block matrix division for the diffOp +% I -- Index of blockmatrix block +function penalty = local2globalPenalty(local, div, I) + penaltyDiv = {div{1}, size(local,2)}; + penalty_bm = blockmatrix.zero(penaltyDiv); + penalty_bm{I,1} = local; + + penalty = blockmatrix.toMatrix(penalty_bm); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/setAllInterfaceTypes.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,13 @@ +% Create interface configuration with a single type for all interfaces +% g -- multiblock grid +% type -- type for all interfaces +function intfTypes = setAllInterfaceTypes(g, type) + intfTypes = cell(g.nBlocks(), g.nBlocks()); + for i = 1:g.nBlocks() + for j = 1:g.nBlocks() + if ~isempty(g.connections{i,j}) + intfTypes{i,j} = type; + end + end + end +end
--- a/+sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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/intOpOP_orders_2to2_ratio2to1.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOpOP_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.intOpOP_orders_2to2_ratio_2to1_accC2F1_accF2C2; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpOP_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/intOpOP_orders_2to2_ratio_2to1_accC2F1_accF2C2.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_2to2_ratio_2to1_accC2F2_accF2C1.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_4to4_ratio2to1.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOpOP_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.intOpOP_orders_4to4_ratio_2to1_accC2F2_accF2C3; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpOP_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/intOpOP_orders_4to4_ratio_2to1_accC2F2_accF2C3.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_4to4_ratio_2to1_accC2F3_accF2C2.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_6to6_ratio2to1.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOpOP_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.intOpOP_orders_6to6_ratio_2to1_accC2F3_accF2C4; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpOP_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/intOpOP_orders_6to6_ratio_2to1_accC2F3_accF2C4.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_6to6_ratio_2to1_accC2F4_accF2C3.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_8to8_ratio2to1.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOpOP_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.intOpOP_orders_8to8_ratio_2to1_accC2F4_accF2C5; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpOP_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/intOpOP_orders_8to8_ratio_2to1_accC2F4_accF2C5.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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/intOpOP_orders_8to8_ratio_2to1_accC2F5_accF2C4.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpOP_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
--- a/+sbp/InterpAWW.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,68 +0,0 @@ -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
--- a/+sbp/InterpMC.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -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
--- a/+sbp/InterpOps.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+sbp/InterpOps.m Sun Dec 23 14:39:31 2018 +0100 @@ -1,9 +1,7 @@ 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" - + Iu2v % Interpolation operator(s) from "u" to "v" + Iv2u % Interpolation operator(s) from "v" to "u" end methods (Abstract)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpOpsMC.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,78 @@ +% Interpolation operators by Mattsson and Carpenter (MC), see +% Mattsson and Carpenter, +% "Stable and Accurate Interpolatino Operators for High-Order Multiblock Finite DIfference Methods", +% https://epubs.siam.org/doi/pdf/10.1137/090750068 +% +% Let ^* denote the adjoint. These operators satsify +% +% Iuv2 = Iv2u^* +% +% Both Iu2v and Iv2u have p:th order accuracy, if the interior stencil is +% of order 2p. +% +% This approach leads to a reduction of the convergence rate by one order for +% PDEs with 2nd derivatives in space, as compared to conforming interfaces. +% To obtain full convergence rate, use the order-preserving (OP) operators in +% InterpOpsOP.m +classdef InterpOpsMC < sbp.InterpOps + properties + + % Structs of interpolation operators, fields .good and .bad + % Here .good and .bad are the same, but this makes them fit in the + % OP (order-preserving) framework. + Iu2v + Iv2u + end + + methods + % m_u, m_v -- number of grid points along the interface + % order_u, order_v -- order of accuracy in the different blocks + function obj = InterpOpsMC(m_u, m_v, order_u, order_v) + + assert(order_u == order_v,... + 'InterpOpsMC: Different orders of accuracy not available'); + + switch order_u + case 2 + intOpSet = @sbp.implementations.intOpMC_orders_2to2_ratio2to1; + case 4 + intOpSet = @sbp.implementations.intOpMC_orders_4to4_ratio2to1; + case 6 + intOpSet = @sbp.implementations.intOpMC_orders_6to6_ratio2to1; + case 8 + intOpSet = @sbp.implementations.intOpMC_orders_8to8_ratio2to1; + otherwise + error('InterpOpsMC: Order of accuracy %d not available.', order_u); + end + + Iu2v = struct; + Iv2u = struct; + + if (m_u-1)/(m_v-1) == 2 + % Block u is fine, v is coarse + m_C = m_v; + [Iv2u.good, Iu2v.bad] = intOpSet(m_C); + Iv2u.bad = Iv2u.good; + Iu2v.good = Iu2v.bad; + + elseif (m_v-1)/(m_u-1) == 2 + % Block v is fine, u is coarse + m_C = m_u; + [Iu2v.good, Iv2u.bad] = intOpSet(m_C); + Iu2v.bad = Iu2v.good; + Iv2u.good = Iv2u.bad; + else + error('InterpOpsMC: Interpolation operators for grid ratio %f have not yet been constructed', (m_u-1)/(m_v-1)); + end + + obj.Iu2v = Iu2v; + obj.Iv2u = Iv2u; + + end + + function str = string(obj) + str = [class(obj)]; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpOpsOP.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,76 @@ +% Order-preserving (OP) interpolation operators, see +% Almquist, Wang, Werpers, +% "Order-Preserving Interpolation for Summation-by-Parts Operators +% at Non-Conforming Interfaces", https://arxiv.org/abs/1806.01931 +% +% Let ^* denote the adjoint. These operators satsify +% +% Iuv2.good = Iv2u.bad^* +% Iv2u.good = Iu2v.bad^* +% +% The .bad operators have the same order of accuracy as the operators +% by Mattsson and Carpenter (MC) in InterpOpsMC, i.e. order p, +% if the interior stencil is order 2p. The .good operators are +% one order more accurate, i.e. order p+1. +% +% For PDEs of second order in space, the OP operators allow for the same +% convergence rate as with conforming interfaces, which is an improvement +% by one order compared what is possible with the MC operators. +classdef InterpOpsOP < sbp.InterpOps + properties + + % Structs of interpolation operators, fields .good and .bad + Iu2v + Iv2u + end + + methods + % m_u, m_v -- number of grid points along the interface + % order_u, order_v -- order of accuracy in the different blocks + function obj = InterpOpsOP(m_u, m_v, order_u, order_v) + + assert(order_u == order_v,... + 'InterpOpsOP: Different orders of accuracy not available'); + + switch order_u + case 2 + intOpSet = @sbp.implementations.intOpOP_orders_2to2_ratio2to1; + case 4 + intOpSet = @sbp.implementations.intOpOP_orders_4to4_ratio2to1; + case 6 + intOpSet = @sbp.implementations.intOpOP_orders_6to6_ratio2to1; + case 8 + intOpSet = @sbp.implementations.intOpOP_orders_8to8_ratio2to1; + otherwise + error('InterpOpsOP: Order of accuracy %d not available.', order_u); + end + + Iu2v = struct; + Iv2u = struct; + + if (m_u-1)/(m_v-1) == 2 + % Block u is fine, v is coarse + m_C = m_v; + [Iv2u.good, Iu2v.bad] = intOpSet(m_C, 1, 'C2F'); + [Iv2u.bad, Iu2v.good] = intOpSet(m_C, 1, 'F2C'); + + elseif (m_v-1)/(m_u-1) == 2 + % Block v is fine, u is coarse + m_C = m_u; + [Iu2v.good, Iv2u.bad] = intOpSet(m_C, 1, 'C2F'); + [Iu2v.bad, Iv2u.good] = intOpSet(m_C, 1, 'F2C'); + else + error('InterpOpsOP: Interpolation operators for grid ratio %f have not yet been constructed', (m_u-1)/(m_v-1)); + end + + obj.Iu2v = Iu2v; + obj.Iv2u = Iv2u; + + end + + function str = string(obj) + str = [class(obj)]; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/+bc/closureSetup.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,25 @@ +% Setup closure and penalty matrices for several boundary conditions at once. +% Each bc is a struct with the fields +% * type -- Type of boundary condition +% * boundary -- Boundary identifier +% * data -- A function_handle for a function which provides boundary data.(see below) +% Also takes S_sign which modifies the sign of the penalty function, [-1,1] +% Returns a closure matrix and a penalty matrices for each boundary condition. +% +% 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, penalties] = closureSetup(diffOp, bcs) + scheme.bc.verifyFormat(bcs, diffOp); + + % Setup storage arrays + closure = spzeros(size(diffOp)); + penalties = cell(1, length(bcs)); + + % Collect closures and penalties + for i = 1:length(bcs) + [localClosure, penalties{i}] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type); + closure = closure + localClosure; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/+bc/forcingSetup.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,86 @@ +% Setup the forcing function for the given boundary conditions and data. +% Each bc is a struct with the fields +% * type -- Type of boundary condition +% * boundary -- Boundary identifier +% * data -- A function_handle for a function which provides boundary data.(see below) +% S_sign allows changing the sign of the function to put on different sides in the system of ODEs. +% default is 1, which the same side as the diffOp. +% Returns 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 S = forcingSetup(diffOp, penalties, bcs, S_sign) + default_arg('S_sign', 1); + + assertType(bcs, 'cell'); + assertIsMember(S_sign, [1, -1]); + + scheme.bc.verifyFormat(bcs, diffOp); + + [gridData, symbolicData] = parseAndSortData(bcs, penalties, diffOp); + + % Setup penalty function + O = spzeros(size(diffOp),1); + function v = S_fun(t) + v = O; + 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 + +% Go through a cell array of boundary condition specifications and return cell arrays +% of structs for grid and symbolic data. +function [gridData, symbolicData] = parseAndSortData(bcs, penalties, diffOp) + gridData = {}; + symbolicData = {}; + for i = 1:length(bcs) + [ok, isSymbolic, data] = parseData(bcs{i}, penalties{i}, diffOp.grid); + + if ~ok + continue % There was no data + end + + if isSymbolic + symbolicData{end+1} = data; + else + gridData{end+1} = data; + end + end +end + +function [ok, isSymbolic, dataStruct] = parseData(bc, penalty, grid) + if ~isfield(bc,'data') || isempty(bc.data) + isSymbolic = []; + dataStruct = struct(); + ok = false; + return + end + ok = true; + + nArg = nargin(bc.data); + + if nArg > 1 + % Symbolic data + isSymbolic = true; + coord = grid.getBoundary(bc.boundary); + dataStruct.penalty = penalty; + dataStruct.func = bc.data; + dataStruct.coords = num2cell(coord, 1); + else + % Grid data + isSymbolic = false; + dataStruct.penalty = penalty; + dataStruct.func = bc.data; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/+bc/verifyFormat.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,31 @@ +% Errors with a more or less detailed error message if there is a problem with the bc specification +function verifyBcFormat(bcs, diffOp) + assertType(bcs, 'cell'); + 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 + + % Find dimension of boundary + b = diffOp.grid.getBoundary(bcs{i}.boundary); + dim = size(b,2); + + % Assert that the data function has a valid number of input arguments + if ~(nargin(bcs{i}.data) == 1 || 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 + + if nargin(bcs{i}.data) == 1 + % Grid data (only function of time) + % Assert that the data has the correct dimension + assertSize(bcs{i}.data(0), 1, size(b,1)); + end + end +end
--- a/+scheme/Beam.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Beam.m Sun Dec 23 14:39:31 2018 +0100 @@ -19,7 +19,7 @@ alphaII alphaIII - opt + opt % TODO: Get rid of this and use the interface type instead end methods @@ -170,7 +170,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
--- a/+scheme/Beam2d.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Beam2d.m Sun Dec 23 14:39:31 2018 +0100 @@ -161,7 +161,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
--- a/+scheme/Euler1d.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Euler1d.m Sun Dec 23 14:39:31 2018 +0100 @@ -446,7 +446,7 @@ closure = @closure_fun; end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) error('NOT DONE') % u denotes the solution in the own domain % v denotes the solution in the neighbour domain
--- a/+scheme/Hypsyst2d.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Hypsyst2d.m Sun Dec 23 14:39:31 2018 +0100 @@ -6,10 +6,10 @@ x,y % Grid X,Y % Values of x and y for each grid point order % Order accuracy for the approximation - + D % non-stabalized scheme operator A, B, E %Coefficient matrices - + H % Discrete norm % Norms in the x and y directions Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -17,65 +17,65 @@ e_w, e_e, e_s, e_n params %parameters for the coeficient matrice end - + methods %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu function obj = Hypsyst2d(m, lim, order, A, B, E, params) default_arg('E', []) xlim = lim{1}; ylim = lim{2}; - + if length(m) == 1 m = [m m]; end - + obj.A=A; obj.B=B; obj.E=E; - + m_x = m(1); m_y = m(2); obj.params = params; - + ops_x = sbp.D2Standard(m_x,xlim,order); ops_y = sbp.D2Standard(m_y,ylim,order); - + obj.x = ops_x.x; obj.y = ops_y.x; - + obj.X = kr(obj.x,ones(m_y,1)); obj.Y = kr(ones(m_x,1),obj.y); - + Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y); Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y); Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y); - + obj.n = length(A(obj.params,0,0)); - + I_n = eye(obj.n);I_x = speye(m_x); obj.I_x = I_x; I_y = speye(m_y); obj.I_y = I_y; - - + + D1_x = kr(I_n, ops_x.D1, I_y); obj.Hxi = kr(I_n, ops_x.HI, I_y); D1_y = kr(I_n, I_x, ops_y.D1); obj.Hyi = kr(I_n, I_x, ops_y.HI); - + obj.e_w = kr(I_n, ops_x.e_l, I_y); obj.e_e = kr(I_n, ops_x.e_r, I_y); obj.e_s = kr(I_n, I_x, ops_y.e_l); obj.e_n = kr(I_n, I_x, ops_y.e_r); - + obj.m = m; obj.h = [ops_x.h ops_y.h]; obj.order = order; - + obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated; - + end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -92,18 +92,18 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + error('Not implemented'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y) params = obj.params; - + if isa(mat,'function_handle') [rows,cols] = size(mat(params,0,0)); matVec = mat(params,X',Y'); @@ -116,7 +116,7 @@ cols = cols/side; end ret = cell(rows,cols); - + for ii = 1:rows for jj=1:cols ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); @@ -124,13 +124,13 @@ end ret = cell2mat(ret); end - + %Characteristic boundary conditions function [closure, penalty] = boundary_condition_char(obj,boundary) params = obj.params; x = obj.x; y = obj.y; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -164,7 +164,7 @@ pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -180,13 +180,13 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,boundary,L) params = obj.params; x = obj.x; y = obj.y; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -218,14 +218,14 @@ boundPos = 'r'; Hi = obj.Hyi; [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); - L = obj.evaluateCoefficientMatrix(L,x,y(end)); + L = obj.evaluateCoefficientMatrix(L,x,y(end)); side = max(length(x)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -233,7 +233,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -243,7 +243,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -251,13 +251,13 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] % D = [ d0 ] % [ d-] - % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y) params = obj.params; syms xs ys @@ -265,12 +265,12 @@ Vi = inv(V); xs = x; ys = y; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret = zeros(obj.n,side*obj.n); - + for ii = 1:obj.n for jj = 1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -278,7 +278,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -286,16 +286,16 @@ Vi = obj.evaluateCoefficientMatrix(Vi,x,y); D = obj.evaluateCoefficientMatrix(D,x,y); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; end - + end end \ No newline at end of file
--- a/+scheme/Hypsyst2dCurve.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Hypsyst2dCurve.m Sun Dec 23 14:39:31 2018 +0100 @@ -4,19 +4,19 @@ n % size of system h % Grid spacing X,Y % Values of x and y for each grid point - + J, Ji % Jacobaian and inverse Jacobian xi,eta Xi,Eta - + A,B X_eta, Y_eta X_xi,Y_xi order % Order accuracy for the approximation - + D % non-stabalized scheme operator Ahat, Bhat, E - + H % Discrete norm Hxii,Hetai % Kroneckerd norms in xi and eta. I_xi,I_eta, I_N, onesN @@ -24,93 +24,93 @@ index_w, index_e,index_s,index_n params % Parameters for the coeficient matrice end - - + + methods % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) default_arg('E', []) xilim = {0 1}; etalim = {0 1}; - + if length(m) == 1 m = [m m]; end obj.params = params; obj.A=A; obj.B=B; - + obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); obj.E=@(params,x,y,~,~)E(params,x,y); - + m_xi = m(1); m_eta = m(2); m_tot=m_xi*m_eta; - + ops_xi = sbp.D2Standard(m_xi,xilim,order); ops_eta = sbp.D2Standard(m_eta,etalim,order); - + obj.xi = ops_xi.x; obj.eta = ops_eta.x; - + obj.Xi = kr(obj.xi,ones(m_eta,1)); obj.Eta = kr(ones(m_xi,1),obj.eta); - + obj.n = length(A(obj.params,0,0)); obj.onesN=ones(obj.n); - + obj.index_w=1:m_eta; - obj.index_e=(m_tot-m_e - + obj.index_e=(m_tot-m_e + metric_termsta+1):m_tot; obj.index_s=1:m_eta:(m_tot-m_eta+1); obj.index_n=(m_eta):m_eta:m_tot; - + I_n = eye(obj.n); I_xi = speye(m_xi); obj.I_xi = I_xi; I_eta = speye(m_eta); obj.I_eta = I_eta; - + D1_xi = kr(I_n, ops_xi.D1, I_eta); obj.Hxii = kr(I_n, ops_xi.HI, I_eta); D1_eta = kr(I_n, I_xi, ops_eta.D1); obj.Hetai = kr(I_n, I_xi, ops_eta.HI); - + obj.e_w = kr(I_n, ops_xi.e_l, I_eta); obj.e_e = kr(I_n, ops_xi.e_r, I_eta); obj.e_s = kr(I_n, I_xi, ops_eta.e_l); - obj.e_n = kr(I_n, I_xi, - + obj.e_n = kr(I_n, I_xi, + metric_termsops_eta.e_r); - + [X,Y] = ti.map(obj.xi,obj.eta); - + [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); - + obj.X = reshape(X,m_tot,1); obj.Y = reshape(Y,m_tot,1); obj.X_xi = reshape(x_xi,m_tot,1); obj.Y_xi = reshape(y_xi,m_tot,1); obj.X_eta = reshape(x_eta,m_tot,1); obj.Y_eta = reshape(y_eta,m_tot,1); - + Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); - + obj.m = m; obj.h = [ops_xi.h ops_eta.h]; obj.order = order; obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi; obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); - + obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; - + end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'. @@ -127,18 +127,18 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions) - error('An interface function does not exist yet'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + error('Not implemented'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) params = obj.params; - + if isa(mat,'function_handle') [rows,cols] = size(mat(params,0,0,0,0)); x_ = kr(obj.onesN,x_); @@ -152,7 +152,7 @@ side = max(length(X),length(Y)); cols = cols/side; end - + ret = cell(rows,cols); for ii = 1:rows for jj = 1:cols @@ -161,7 +161,7 @@ end ret = cell2mat(ret); end - + %Characteristic boundary conditions function [closure, penalty] = boundary_condition_char(obj,boundary) params = obj.params; @@ -169,7 +169,7 @@ Y = obj.Y; xi = obj.xi; eta = obj.eta; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -200,11 +200,11 @@ [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); side = max(length(xi)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -218,10 +218,10 @@ Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); closure = Hi*e_*V*tau*Vi_minus*e_'; penalty = -Hi*e_*V*tau*Vi_minus; - end + end end - - + + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,boundary,L) params = obj.params; @@ -229,7 +229,7 @@ Y = obj.Y; xi = obj.xi; eta = obj.eta; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -237,7 +237,7 @@ boundPos = 'l'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_w)); xi_x = Ji*obj.Y_eta(obj.index_w); @@ -250,7 +250,7 @@ boundPos = 'r'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_e)); xi_x = Ji*obj.Y_eta(obj.index_e); @@ -263,7 +263,7 @@ boundPos = 'l'; Hi = obj.Hetai; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_s)); eta_x = Ji*obj.Y_xi(obj.index_s); @@ -276,7 +276,7 @@ boundPos = 'r'; Hi = obj.Hetai; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_n)); eta_x = Ji*obj.Y_xi(obj.index_n); @@ -284,11 +284,11 @@ L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); side = max(length(xi)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -296,7 +296,7 @@ Vi_minus = Vi(pos+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -306,7 +306,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -314,7 +314,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -329,22 +329,22 @@ else xs_ = 0; end - + if(sum(abs(y_))~= 0) syms ys_; else ys_ = 0; end - + [V, D] = eig(mat(params,xs,ys,xs_,ys_)); Vi = inv(V); syms xs ys xs_ ys_ - + xs = x; ys = y; xs_ = x_; ys_ = y_; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); @@ -356,7 +356,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -364,11 +364,11 @@ D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_); Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- a/+scheme/Hypsyst3d.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Hypsyst3d.m Sun Dec 23 14:39:31 2018 +0100 @@ -7,11 +7,11 @@ X, Y, Z% Values of x and y for each grid point Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces order % Order accuracy for the approximation - + D % non-stabalized scheme operator A, B, C, E % Symbolic coefficient matrices Aevaluated,Bevaluated,Cevaluated, Eevaluated - + H % Discrete norm Hx, Hy, Hz % Norms in the x, y and z directions Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -19,8 +19,8 @@ e_w, e_e, e_s, e_n, e_b, e_t params % Parameters for the coeficient matrice end - - + + methods % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) @@ -28,11 +28,11 @@ xlim = lim{1}; ylim = lim{2}; zlim = lim{3}; - + if length(m) == 1 m = [m m m]; end - + obj.A = A; obj.B = B; obj.C = C; @@ -41,7 +41,7 @@ m_y = m(2); m_z = m(3); obj.params = params; - + switch operator case 'upwind' ops_x = sbp.D1Upwind(m_x,xlim,order); @@ -52,29 +52,29 @@ ops_y = sbp.D2Standard(m_y,ylim,order); ops_z = sbp.D2Standard(m_z,zlim,order); end - + obj.x = ops_x.x; obj.y = ops_y.x; obj.z = ops_z.x; - + obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1)); obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); - + obj.Yx = kr(obj.y,ones(m_z,1)); obj.Zx = kr(ones(m_y,1),obj.z); obj.Xy = kr(obj.x,ones(m_z,1)); obj.Zy = kr(ones(m_x,1),obj.z); obj.Xz = kr(obj.x,ones(m_y,1)); obj.Yz = kr(ones(m_z,1),obj.y); - + obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); - + obj.n = length(A(obj.params,0,0,0)); - + I_n = speye(obj.n); I_x = speye(m_x); obj.I_x = I_x; @@ -83,31 +83,31 @@ I_z = speye(m_z); obj.I_z = I_z; I_N = kr(I_n,I_x,I_y,I_z); - + obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); obj.Hx = ops_x.H; obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); obj.Hy = ops_y.H; obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); obj.Hz = ops_z.H; - + obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); - + obj.m = m; obj.h = [ops_x.h ops_y.h ops_x.h]; obj.order = order; - + switch operator case 'upwind' alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end))))); alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end))))); alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end))))); - + Ap = (obj.Aevaluated+alphaA*I_N)/2; Am = (obj.Aevaluated-alphaA*I_N)/2; Dpx = kr(I_n, ops_x.Dp, I_y,I_z); @@ -116,7 +116,7 @@ temp = Ap*Dmx; obj.D = obj.D-temp; clear Ap Am Dpx Dmx - + Bp = (obj.Bevaluated+alphaB*I_N)/2; Bm = (obj.Bevaluated-alphaB*I_N)/2; Dpy = kr(I_n, I_x, ops_y.Dp,I_z); @@ -126,20 +126,20 @@ temp = Bp*Dmy; obj.D = obj.D-temp; clear Bp Bm Dpy Dmy - - + + Cp = (obj.Cevaluated+alphaC*I_N)/2; Cm = (obj.Cevaluated-alphaC*I_N)/2; Dpz = kr(I_n, I_x, I_y,ops_z.Dp); Dmz = kr(I_n, I_x, I_y,ops_z.Dm); - + temp = Cm*Dpz; obj.D = obj.D-temp; temp = Cp*Dmz; obj.D = obj.D-temp; clear Cp Cm Dpz Dmz obj.D = obj.D-obj.Eevaluated; - + case 'standard' D1_x = kr(I_n, ops_x.D1, I_y,I_z); D1_y = kr(I_n, I_x, ops_y.D1,I_z); @@ -149,7 +149,7 @@ error('Opperator not supported'); end end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -167,15 +167,15 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + error('Not implemented'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) params = obj.params; side = max(length(X),length(Y)); @@ -189,7 +189,7 @@ side = max(length(X),length(Y)); cols = cols/side; end - + ret = cell(rows,cols); for ii = 1:rows for jj = 1:cols @@ -198,10 +198,10 @@ end ret = cell2mat(ret); end - + function [BM] = boundary_matrices(obj,boundary) params = obj.params; - + switch boundary case {'w','W','west'} BM.e_ = obj.e_w; @@ -248,7 +248,7 @@ end BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - + % Characteristic bouyndary consitions function [closure, penalty]=boundary_condition_char(obj,BM) side = BM.side; @@ -260,7 +260,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -276,9 +276,9 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) - function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) + function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) side = BM.side; pos = BM.pos; neg = BM.neg; @@ -288,7 +288,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch boundary case {'w','W','west'} L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); @@ -303,7 +303,7 @@ case {'t','T','top'} L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); end - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -311,7 +311,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -321,7 +321,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -329,7 +329,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -344,13 +344,13 @@ xs = x; ys = y; zs = z; - - + + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret= zeros(obj.n,side*obj.n); - + for ii=1:obj.n for jj=1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -358,7 +358,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -366,11 +366,11 @@ Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z); D = obj.evaluateCoefficientMatrix(D,x,y,z); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- a/+scheme/Hypsyst3dCurve.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Hypsyst3dCurve.m Sun Dec 23 14:39:31 2018 +0100 @@ -5,22 +5,22 @@ h % Grid spacing X, Y, Z% Values of x and y for each grid point Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces - + xi,eta,zeta Xi, Eta, Zeta - + Eta_xi, Zeta_xi, Xi_eta, Zeta_eta, Xi_zeta, Eta_zeta % Metric terms X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta % Metric terms - + order % Order accuracy for the approximation - + D % non-stabalized scheme operator Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices A, B, C, E % Symbolic coeffiecient matrices - + J, Ji % JAcobian and inverse Jacobian - + H % Discrete norm % Norms in the x, y and z directions Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -30,14 +30,14 @@ index_w, index_e,index_s,index_n, index_b, index_t params %parameters for the coeficient matrice end - - + + methods function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator) xilim ={0 1}; etalim = {0 1}; zetalim = {0 1}; - + if length(m) == 1 m = [m m m]; end @@ -47,11 +47,11 @@ m_tot = m_xi*m_eta*m_zeta; obj.params = params; obj.n = length(A(obj,0,0,0)); - + obj.m = m; obj.order = order; obj.onesN = ones(obj.n); - + switch operator case 'upwind' ops_xi = sbp.D1Upwind(m_xi,xilim,order); @@ -64,21 +64,21 @@ otherwise error('Operator not available') end - + obj.xi = ops_xi.x; obj.eta = ops_eta.x; obj.zeta = ops_zeta.x; - + obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1)); obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); - - + + [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); obj.X = X; obj.Y = Y; obj.Z = Z; - + I_n = eye(obj.n); I_xi = speye(m_xi); obj.I_xi = I_xi; @@ -86,19 +86,19 @@ obj.I_eta = I_eta; I_zeta = speye(m_zeta); obj.I_zeta = I_zeta; - + I_N=kr(I_n,I_xi,I_eta,I_zeta); - + O_xi = ones(m_xi,1); O_eta = ones(m_eta,1); O_zeta = ones(m_zeta,1); - - + + obj.Hxi = ops_xi.H; obj.Heta = ops_eta.H; obj.Hzeta = ops_zeta.H; obj.h = [ops_xi.h ops_eta.h ops_zeta.h]; - + switch operator case 'upwind' D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta); @@ -109,11 +109,11 @@ D1_eta = kr(I_xi, ops_eta.D1,I_zeta); D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); end - + obj.A = A; obj.B = B; obj.C = C; - + obj.X_xi = D1_xi*X; obj.X_eta = D1_eta*X; obj.X_zeta = D1_zeta*X; @@ -123,55 +123,55 @@ obj.Z_xi = D1_xi*Z; obj.Z_eta = D1_eta*Z; obj.Z_zeta = D1_zeta*Z; - + obj.Ahat = @transform_coefficient_matrix; obj.Bhat = @transform_coefficient_matrix; obj.Chat = @transform_coefficient_matrix; obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); - + obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); - + switch operator case 'upwind' clear D1_xi D1_eta D1_zeta alphaA = max(abs(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))))); alphaB = max(abs(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))))); alphaC = max(abs(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))))); - + Ap = (obj.Aevaluated+alphaA*I_N)/2; Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); diffSum = -Ap*Dmxi; clear Ap Dmxi - + Am = (obj.Aevaluated-alphaA*I_N)/2; - + obj.Aevaluated = []; Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); temp = Am*Dpxi; diffSum = diffSum-temp; clear Am Dpxi - + Bp = (obj.Bevaluated+alphaB*I_N)/2; Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); temp = Bp*Dmeta; diffSum = diffSum-temp; clear Bp Dmeta - + Bm = (obj.Bevaluated-alphaB*I_N)/2; obj.Bevaluated = []; Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); temp = Bm*Dpeta; diffSum = diffSum-temp; clear Bm Dpeta - + Cp = (obj.Cevaluated+alphaC*I_N)/2; Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); temp = Cp*Dmzeta; diffSum = diffSum-temp; clear Cp Dmzeta - + Cm = (obj.Cevaluated-alphaC*I_N)/2; clear I_N obj.Cevaluated = []; @@ -179,72 +179,72 @@ temp = Cm*Dpzeta; diffSum = diffSum-temp; clear Cm Dpzeta temp - + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; - + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); - + obj.D = obj.Ji*diffSum-obj.Eevaluated; - + case 'standard' D1_xi = kr(I_n,D1_xi); D1_eta = kr(I_n,D1_eta); D1_zeta = kr(I_n,D1_zeta); - + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; - + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); - + obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; otherwise error('Operator not supported') end - + obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); - + obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1); obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1); obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1); obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1); obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1); - + obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta); obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta); obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta); obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta); obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l); obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r); - + obj.Eta_xi = kr(obj.eta,ones(m_xi,1)); obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta); obj.Xi_eta = kr(obj.xi,ones(m_zeta,1)); obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta); obj.Xi_zeta = kr(obj.xi,ones(m_eta,1)); - obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); + obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); end - + function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); ret = ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); ret = ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); end - - + + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -253,7 +253,7 @@ function [closure, penalty] = boundary_condition(obj,boundary,type,L) default_arg('type','char'); BM = boundary_matrices(obj,boundary); - + switch type case{'c','char'} [closure,penalty] = boundary_condition_char(obj,BM); @@ -263,15 +263,15 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + error('Not implemented'); end - + function N = size(obj) N = obj.m; end - + % Evaluates the symbolic Coeffiecient matrix mat function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) params = obj.params; @@ -294,7 +294,7 @@ end matVec(abs(matVec)<10^(-10)) = 0; ret = cell(rows,cols); - + for ii = 1:rows for jj = 1:cols ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); @@ -302,7 +302,7 @@ end ret = cell2mat(ret); end - + function [BM] = boundary_matrices(obj,boundary) params = obj.params; BM.boundary = boundary; @@ -385,7 +385,7 @@ BM.side = sum(BM.index); BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - + % Characteristic boundary condition function [closure, penalty] = boundary_condition_char(obj,BM) side = BM.side; @@ -397,7 +397,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -413,7 +413,7 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) side = BM.side; @@ -426,7 +426,7 @@ D = BM.D; e_ = BM.e_; index = BM.index; - + switch BM.boundary case{'b','B','bottom'} Ji_vec = diag(obj.Ji); @@ -434,10 +434,10 @@ Zeta_x = Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); Zeta_y = Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); Zeta_z = Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); - + L = obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); end - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -445,7 +445,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -455,7 +455,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -463,7 +463,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -478,38 +478,38 @@ else x_1s = 0; end - + if(sum(abs(x_2))>eps) syms x_2s; else x_2s = 0; end - - + + if(sum(abs(y_1))>eps) syms y_1s else y_1s = 0; end - + if(sum(abs(y_2))>eps) syms y_2s; else y_2s = 0; end - + if(sum(abs(z_1))>eps) syms z_1s else z_1s = 0; end - + if(sum(abs(z_2))>eps) syms z_2s; else z_2s = 0; end - + syms xs ys zs [V, D] = eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); Vi = inv(V); @@ -522,12 +522,12 @@ y_2s = y_2; z_1s = z_1; z_2s = z_2; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret = zeros(obj.n,side*obj.n); - + for ii=1:obj.n for jj=1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -535,7 +535,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -543,11 +543,11 @@ D = obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); Vi = obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Laplace1d.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,146 @@ +classdef Laplace1d < scheme.Scheme + properties + grid + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + H % Discrete norm + M % Derivative norm + a + + D2 + Hi + e_l + e_r + d_l + d_r + gamm + end + + methods + function obj = Laplace1d(grid, order, a) + default_arg('a', 1); + + assertType(grid, 'grid.Cartesian'); + + ops = sbp.D2Standard(grid.size(), grid.lim{1}, order); + + obj.D2 = sparse(ops.D2); + obj.H = sparse(ops.H); + obj.Hi = sparse(ops.HI); + obj.M = sparse(ops.M); + obj.e_l = sparse(ops.e_l); + obj.e_r = sparse(ops.e_r); + obj.d_l = -sparse(ops.d1_l); + obj.d_r = sparse(ops.d1_r); + + + obj.grid = grid; + obj.order = order; + + obj.a = a; + obj.D = a*obj.D2; + + obj.gamm = grid.h*ops.borrowing.M.S; + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % 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); + + [e,d,s] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','dirichlet'} + tuning = 1.1; + tau1 = -tuning/obj.gamm; + tau2 = 1; + + tau = tau1*e + tau2*d; + + closure = obj.a*obj.Hi*tau*e'; + penalty = obj.a*obj.Hi*tau; + + % Neumann boundary condition + case {'N','neumann'} + tau = -e; + + closure = obj.a*obj.Hi*tau*d'; + penalty = -obj.a*obj.Hi*tau; + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + + [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + + a_u = obj.a; + a_v = neighbour_scheme.a; + + gamm_u = obj.gamm; + gamm_v = neighbour_scheme.gamm; + + tuning = 1.1; + + tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning; + tau2 = 1/2*a_u; + sig1 = -1/2; + sig2 = 0; + + tau = tau1*e_u + tau2*d_u; + sig = sig1*e_u + sig2*d_u; + + closure = obj.Hi*( tau*e_u' + sig*a_u*d_u'); + penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d,s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d = obj.d_l; + s = -1; + case 'r' + e = obj.e_r; + d = obj.d_r; + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = obj.grid.size(); + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + 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
--- a/+scheme/LaplaceCurvilinear.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/LaplaceCurvilinear.m Sun Dec 23 14:39:31 2018 +0100 @@ -38,6 +38,7 @@ du_n, dv_n gamm_u, gamm_v lambda + end methods @@ -53,7 +54,11 @@ error('Not implemented yet') end - assert(isa(g, 'grid.Curvilinear')) + % assert(isa(g, 'grid.Curvilinear')) + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end m = g.size(); m_u = m(1); @@ -268,11 +273,31 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.2 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.2; + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + tuning = type.tuning; + % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - tuning = 1.2; - % tuning = 20.2; [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); @@ -298,15 +323,66 @@ penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); end - % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + % TODO: Make this work for curvilinear grids + warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); + + % User can request special interpolation operators by specifying type.interpOpSet + default_field(type, 'interpOpSet', @sbp.InterpOpsOP); + interpOpSet = type.interpOpSet; + tuning = type.tuning; + + + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); + [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + % Find the number of grid points along the interface + m_u = size(e_u, 2); + m_v = size(e_v, 2); + + Hi = obj.Hi; + a = obj.a; + + u = obj; + v = neighbour_scheme; + + b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; + b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; + b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; + b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; + + tau_u = -1./(4*b1_u) -1./(4*b2_u); + tau_v = -1./(4*b1_v) -1./(4*b2_v); + + tau_u = tuning * spdiag(tau_u); + tau_v = tuning * spdiag(tau_v); + beta_u = tau_v; + + % Build interpolation operators + intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); + Iu2v = intOps.Iu2v; + Iv2u = intOps.Iv2u; + + closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... + a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... + a*1/2*Hi*d_u*H_b_u*e_u' + ... + -a*1/2*Hi*e_u*H_b_u*d_u'; + + penalty = -a*Hi*e_u*tau_u*H_b_u*Iv2u.good*e_v' + ... + -a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*e_v' + ... + -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... + -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; + + end + + % Returns the boundary ops and sign for the boundary specified by the string boundary. % The right boundary is considered the positive boundary % - % I -- the indecies of the boundary points in the grid matrix + % I -- the indices of the boundary points in the grid matrix function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) - - % gridMatrix = zeros(obj.m(2),obj.m(1)); - % gridMatrix(:) = 1:numel(gridMatrix); - ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); switch boundary
--- a/+scheme/Scheme.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Scheme.m Sun Dec 23 14:39:31 2018 +0100 @@ -26,7 +26,10 @@ % interface to. % penalty may be a cell array if there are several penalties with different weights [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition - [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % type -- sets the type of interface, could be a string or a struct or something else + % depending on the particular scheme implementation + [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % TODO: op = getBoundaryOperator()?? % makes sense to have it available through a method instead of random properties
--- a/+scheme/Schrodinger.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Schrodinger.m Sun Dec 23 14:39:31 2018 +0100 @@ -90,7 +90,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2d.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,326 @@ +classdef Schrodinger2d < scheme.Scheme + +% Discretizes the Laplacian with constant coefficent, +% in the Schrödinger equation way (i.e., the discretization matrix is not necessarily +% definite) +% u_t = a*i*Laplace u +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrix for variable coefficients + a % Constant coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2 + + H, Hi % Inner products + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + e_w, e_e, e_s, e_n + d_w, d_e, d_s, d_n + + H_boundary % Boundary inner products + + end + + methods + + function obj = Schrodinger2d(g ,order, a, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('a',1); + dim = 2; + + assertType(g, 'grid.Cartesian'); + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end + + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + lim = {xlim, ylim}; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + % Constant coeff D2 + for i = 1:dim + D2{i} = D2{i}(ones(m(i),1)); + end + + %====== Assemble full operators ======== + obj.D1 = cell(dim,1); + obj.D2 = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 + obj.D2{1} = kron(D2{1},I{2}); + obj.D2{2} = kron(I{1},D2{2}); + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{2} = H{1}; + + % Differentiation matrix D (without SAT) + D2 = obj.D2; + D = sparse(m_tot,m_tot); + for j = 1:dim + D = D + a*1i*D2{j}; + end + obj.D = D; + %=========================================% + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + obj.a = a; + obj.e_w = obj.e_l{1}; + obj.e_e = obj.e_r{1}; + obj.e_s = obj.e_l{2}; + obj.e_n = obj.e_r{2}; + obj.d_w = obj.d1_l{1}; + obj.d_e = obj.d1_r{1}; + obj.d_s = obj.d1_l{2}; + obj.d_n = obj.d1_r{2}; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators 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. + % 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, parameter) + default_arg('type','Neumann'); + default_arg('parameter', []); + + % j is the coordinate direction of the boundary + % nj: outward unit normal component. + % nj = -1 for west, south, bottom boundaries + % nj = 1 for east, north, top boundaries + [j, nj] = obj.get_boundary_number(boundary); + switch nj + case 1 + e = obj.e_r; + d = obj.d1_r; + case -1 + e = obj.e_l; + d = obj.d1_l; + end + + Hi = obj.Hi; + H_gamma = obj.H_boundary{j}; + a = e{j}'*obj.a*e{j}; + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); + penalty = -nj*Hi*d{j}*a*1i*H_gamma; + + % Free boundary condition + case {'N','n','neumann','Neumann'} + closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*a*1i*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + % type Struct that specifies the interface coupling. + % Fields: + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + + % Get boundary operators + [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e, d, H_gamma] = obj.get_boundary_ops(boundary); + Hi = obj.Hi; + a = obj.a; + + % Get outward unit normal component + [~, n] = obj.get_boundary_number(boundary); + + Hi = obj.Hi; + sigma = -n*1i*a/2; + tau = -n*(1i*a)'/2; + + closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; + penalty = -tau*Hi*d*H_gamma*e_neighbour' ... + -sigma*Hi*e*H_gamma*d_neighbour'; + + end + + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + % User can request special interpolation operators by specifying type.interpOpSet + default_field(type, 'interpOpSet', @sbp.InterpOpsOP); + interpOpSet = type.interpOpSet; + + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary); + Hi = obj.Hi; + a = obj.a; + + % Get outward unit normal component + [~, n] = obj.get_boundary_number(boundary); + + % Find the number of grid points along the interface + m_u = size(e_u, 2); + m_v = size(e_v, 2); + + % Build interpolation operators + intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); + Iu2v = intOps.Iu2v; + Iv2u = intOps.Iv2u; + + sigma = -n*1i*a/2; + tau = -n*(1i*a)'/2; + + closure = tau*Hi*d_u*H_gamma*e_u' + sigma*Hi*e_u*H_gamma*d_u'; + penalty = -tau*Hi*d_u*H_gamma*Iv2u.good*e_v' ... + -sigma*Hi*e_u*H_gamma*Iv2u.bad*d_v'; + + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [j, nj] = get_boundary_number(obj, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + end + end + + % Returns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e, d, H_b] = get_boundary_ops(obj, boundary) + + switch boundary + case 'w' + e = obj.e_w; + d = obj.d_w; + H_b = obj.H_boundary{1}; + case 'e' + e = obj.e_e; + d = obj.d_e; + H_b = obj.H_boundary{1}; + case 's' + e = obj.e_s; + d = obj.d_s; + H_b = obj.H_boundary{2}; + case 'n' + e = obj.e_n; + d = obj.d_n; + H_b = obj.H_boundary{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- a/+scheme/Utux.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Utux.m Sun Dec 23 14:39:31 2018 +0100 @@ -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 @@ -16,42 +16,30 @@ end - 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); - - - switch operator - 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; - otherwise - error('Unvalid operator') - end - obj.x=ops.x; + methods + function obj = Utux(g, order, opSet) + default_arg('opSet',@sbp.D2Standard); - + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + + ops = opSet(m, xlim, order); + obj.D1 = ops.D1; + + 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,19 +49,29 @@ % 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); - tau =-1*obj.e_l; - closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; - + 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 = -obj.Hi*tau; + end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + 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) N = obj.m; end @@ -81,9 +79,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 Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,308 @@ +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]; + % Can either be a constant vector or a cell array of function handles. + + 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 + end + + + methods + function obj = Utux2d(g ,order, opSet, a) + + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + + assertType(g, 'grid.Cartesian'); + if iscell(a) + a1 = grid.evalOn(g, a{1}); + a2 = grid.evalOn(g, a{2}); + a = {spdiag(a1), spdiag(a2)}; + else + a = {a(1), a(2)}; + end + + 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.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 + + % type Struct that specifies the interface coupling. + % Fields: + % -- couplingType String, type of interface coupling + % % Default: 'upwind'. Other: 'centered' + % -- interpolation: type of interpolation, default 'none' + % -- interpolationDamping: damping on upstream and downstream sides, when using interpolation. + % Default {0,0} gives zero damping. + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.couplingType = 'upwind'; + defaultType.interpolation = 'none'; + defaultType.interpolationDamping = {0,0}; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + couplingType = type.couplingType; + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + end + + switch couplingType + + % 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 ' couplingType ' is not available.']) + end + + 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*e_neighbour'; + 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*e_neighbour'; + 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*e_neighbour'; + 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*e_neighbour'; + end + + end + + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + % User can request special interpolation operators by specifying type.interpOpSet + default_field(type, 'interpOpSet', @sbp.InterpOpsOP); + + interpOpSet = type.interpOpSet; + couplingType = type.couplingType; + interpolationDamping = type.interpolationDamping; + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + end + + switch couplingType + + % 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 ' couplingType ' is not available.']) + end + + int_damp_us = interpolationDamping{1}; + int_damp_ds = interpolationDamping{2}; + + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + % Find the number of grid points along the interface + switch boundary + case {'w','e'} + m_u = obj.m(2); + case {'s','n'} + m_u = obj.m(1); + end + m_v = size(e_neighbour, 2); + + % Build interpolation operators + intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order); + Iu2v = intOps.Iu2v; + Iv2u = intOps.Iv2u; + + I_local2neighbour_ds = intOps.Iu2v.bad; + I_local2neighbour_us = intOps.Iu2v.good; + I_neighbour2local_ds = intOps.Iv2u.good; + I_neighbour2local_us = intOps.Iv2u.bad; + + 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*obj.e_w' - obj.Hi*beta*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*obj.e_e' - obj.Hi*beta*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*obj.e_s' - obj.Hi*beta*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*obj.e_n' - obj.Hi*beta*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
--- a/+scheme/Wave.m Sun Dec 23 14:06:26 2018 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,175 +0,0 @@ -classdef Wave < scheme.Scheme - properties - m % Number of points in each direction, possibly a vector - h % Grid spacing - x % Grid - order % Order accuracy for the approximation - - D % non-stabalized scheme operator - H % Discrete norm - M % Derivative norm - alpha - - D2 - Hi - e_l - e_r - d1_l - d1_r - gamm - end - - methods - function obj = Wave(m,xlim,order,alpha) - default_arg('a',1); - [x, h] = util.get_grid(xlim{:},m); - - ops = sbp.Ordinary(m,h,order); - - obj.D2 = sparse(ops.derivatives.D2); - obj.H = sparse(ops.norms.H); - obj.Hi = sparse(ops.norms.HI); - obj.M = sparse(ops.norms.M); - obj.e_l = sparse(ops.boundary.e_1); - obj.e_r = sparse(ops.boundary.e_m); - obj.d1_l = sparse(ops.boundary.S_1); - obj.d1_r = sparse(ops.boundary.S_m); - - - obj.m = m; - obj.h = h; - obj.order = order; - - obj.alpha = alpha; - obj.D = alpha*obj.D2; - obj.x = x; - - obj.gamm = h*ops.borrowing.M.S; - - end - - - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % 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); - - [e,d,s] = obj.get_boundary_ops(boundary); - - switch type - % Dirichlet boundary condition - case {'D','dirichlet'} - alpha = obj.alpha; - - % tau1 < -alpha^2/gamma - tuning = 1.1; - tau1 = -tuning*alpha/obj.gamm; - tau2 = s*alpha; - - p = tau1*e + tau2*d; - - closure = obj.Hi*p*e'; - - pp = obj.Hi*p; - switch class(data) - case 'double' - penalty = pp*data; - case 'function_handle' - penalty = @(t)pp*data(t); - otherwise - error('Wierd data argument!') - end - - - % Neumann boundary condition - case {'N','neumann'} - alpha = obj.alpha; - tau1 = -s*alpha; - tau2 = 0; - tau = tau1*e + tau2*d; - - closure = obj.Hi*tau*d'; - - pp = obj.Hi*tau; - switch class(data) - case 'double' - penalty = pp*data; - case 'function_handle' - penalty = @(t)pp*data(t); - otherwise - error('Wierd data argument!') - end - - % Unknown, boundary condition - otherwise - error('No such boundary condition: type = %s',type); - end - end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - % u denotes the solution in the own domain - % v denotes the solution in the neighbour domain - [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); - [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - - tuning = 1.1; - - alpha_u = obj.alpha; - alpha_v = neighbour_scheme.alpha; - - gamm_u = obj.gamm; - gamm_v = neighbour_scheme.gamm; - - % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) - - tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; - tau2 = s_u*1/2*alpha_u; - sig1 = s_u*(-1/2); - sig2 = 0; - - tau = tau1*e_u + tau2*d_u; - sig = sig1*e_u + sig2*d_u; - - closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u'); - penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v'); - end - - % Ruturns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive boundary - function [e,d,s] = get_boundary_ops(obj,boundary) - switch boundary - case 'l' - e = obj.e_l; - d = obj.d1_l; - s = -1; - case 'r' - e = obj.e_r; - d = obj.d1_r; - s = 1; - otherwise - error('No such boundary: boundary = %s',boundary); - end - end - - function N = size(obj) - N = obj.m; - end - - end - - methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') - 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
--- a/+scheme/Wave2d.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Wave2d.m Sun Dec 23 14:39:31 2018 +0100 @@ -158,7 +158,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
--- a/+scheme/Wave2dCurve.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/Wave2dCurve.m Sun Dec 23 14:39:31 2018 +0100 @@ -243,7 +243,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain tuning = 1.2;
--- a/+scheme/bcSetup.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+scheme/bcSetup.m Sun Dec 23 14:39:31 2018 +0100 @@ -1,10 +1,9 @@ -% function [closure, S] = bcSetup(diffOp, bc) % Takes a diffOp and a cell array of boundary condition definitions. % Each bc is a struct with the fields % * type -- Type of boundary condition % * boundary -- Boundary identifier % * data -- A function_handle for a function which provides boundary data.(see below) -% Also takes S_sign which modifies the sign of S, [-1,1] +% Also takes S_sign which modifies the sign of the penalty function, [-1,1] % Returns a closure matrix and a forcing function S. % % The boundary data function can either be a function of time or a function of time and space coordinates. @@ -16,97 +15,6 @@ 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)); - gridData = {}; - symbolicData = {}; - - % Collect closures, penalties and data - for i = 1:length(bcs) - [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type); - closure = closure + localClosure; - - [ok, isSymbolic, data] = parseData(bcs{i}, penalty, diffOp.grid); - - if ~ok - % There was no data - continue - end - - if isSymbolic - 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(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; + [closure, penalties] = scheme.bc.closureSetup(diffOp, bcs); + S = scheme.bc.forcingSetup(diffOp, penalties, bcs, S_sign); 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, isSymbolic, dataStruct] = parseData(bc, penalty, grid) - if ~isfield(bc,'data') || isempty(bc.data) - isSymbolic = []; - dataStruct = struct(); - ok = false; - return - end - ok = true; - - nArg = nargin(bc.data); - - if nArg > 1 - % Symbolic data - isSymbolic = true; - coord = grid.getBoundary(bc.boundary); - dataStruct.penalty = penalty; - dataStruct.func = bc.data; - dataStruct.coords = num2cell(coord, 1); - else - % Grid data - isSymbolic = false; - dataStruct.penalty = penalty; - dataStruct.func = bcs{i}.data; - end -end
--- a/+time/CdiffImplicit.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+time/CdiffImplicit.m Sun Dec 23 14:39:31 2018 +0100 @@ -17,6 +17,8 @@ % u(t0) = f1 % u_t(t0) = f2 % starting at time t0 with timestep k + + % TODO: Fix order of matrices function obj = CdiffImplicit(A, B, C, G, f1, f2, k, t0) default_arg('A', []); default_arg('C', []);
--- a/+util/ReplaceableString.m Sun Dec 23 14:06:26 2018 +0100 +++ b/+util/ReplaceableString.m Sun Dec 23 14:39:31 2018 +0100 @@ -58,3 +58,5 @@ function b = padStr(a, n) b = sprintf('%-*s', n, a); end + +% TODO: Add a debug mode which prints without replacing?
--- a/.hgtags Sun Dec 23 14:06:26 2018 +0100 +++ b/.hgtags Sun Dec 23 14:39:31 2018 +0100 @@ -1,4 +1,5 @@ 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1 0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 +b723495cdb2f96314d7b3f0aa79723a7dc088c7d v0.2 08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0 08f3ffe63f484d02abce8df4df61e826f568193f Heimisson2018
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/LICENSE.txt Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,25 @@ +MIT License + +Copyright (c) +2015-2018 Jonatan Werpers +2015-2018 Martin Almquist +2016-2018 Ylva Rydin +2018 Vidar Stiernström + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE.
--- a/Map.m Sun Dec 23 14:06:26 2018 +0100 +++ b/Map.m Sun Dec 23 14:39:31 2018 +0100 @@ -59,6 +59,9 @@ function v = subsref(obj, S) switch S(1).type case '()' + if length(S.subs) > 1 + error('sbplib:Map:multipleKeys', 'Multiple dimensions are not supported. Use a cell array as a key instead.'); + end k = S.subs{1}; try v = get(obj, k); @@ -81,6 +84,9 @@ function obj = subsasgn(obj, S, v); switch S(1).type case '()' + if length(S.subs) > 1 + error('sbplib:Map:multipleKeys', 'Multiple dimensions are not supported. Use a cell array as a key instead.'); + end k = S.subs{1}; set(obj, k, v); otherwise
--- a/MapTest.m Sun Dec 23 14:06:26 2018 +0100 +++ b/MapTest.m Sun Dec 23 14:39:31 2018 +0100 @@ -12,6 +12,21 @@ }; end +function testMultiKey(testCase) + map = Map + + function setMultiKey() + map(1,2) = 1; + end + + function getMultiKey() + v = map(1,2); + end + + testCase.verifyError(@setMultiKey,'sbplib:Map:multipleKeys') + testCase.verifyError(@getMultiKey,'sbplib:Map:multipleKeys') +end + function testSetAndGet(testCase) keyValuePairs = getKeyValuePairs();
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.md Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,2 @@ +# SBPLIB +sbplib is a library of primitives and help functions for working with summation-by-parts finite differences in Matlab. To use sbplib download the code and add the sbplib folder to the matlab path.
--- a/TextTable.m Sun Dec 23 14:06:26 2018 +0100 +++ b/TextTable.m Sun Dec 23 14:39:31 2018 +0100 @@ -41,6 +41,14 @@ obj.fmtArray{i,j} = fmt; end + function formatGrid(obj, I, J, fmt) + for i = I + for j = J + obj.fmtArray{i,j} = fmt; + end + end + end + function formatRow(obj, i, fmt) obj.fmtArray(i,:) = {fmt}; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/arrowAnnotation.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,23 @@ +% Draw an arrow from p1 to p2, with text attached +function [h] = arrowAnnotation(p1,p2) + ah = gca; + xl = ah.XLim(1); + xr = ah.XLim(2); + + yl = ah.YLim(1); + yr = ah.YLim(2); + + dx = xr - xl; + dy = yr - yl; + + s = [ + ah.Position(1) + (p1(1) - xl)/dx*ah.Position(3), + ah.Position(1) + (p2(1) - xl)/dx*ah.Position(3), + ]; + t = [ + ah.Position(2) + (p1(2) - yl)/dy*ah.Position(4), + ah.Position(2) + (p2(2) - yl)/dy*ah.Position(4), + ]; + + h = annotation('arrow', s, t); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dealStruct.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,18 @@ +function varargout = dealStruct(s, fields) + default_arg('fields', []); + + if isempty(fields) + out = dealFields(s, fieldnames(s)); + varargout = out(1:nargout); + else + assert(nargout == length(fields), 'Number of output arguements must match the number of fieldnames provided'); + varargout = dealFields(s, fields); + end +end + +function out = dealFields(s, fields) + out = cell(1, length(fields)); + for i = 1:length(fields) + out{i} = s.(fields{i}); + end +end
--- a/default_field.m Sun Dec 23 14:06:26 2018 +0100 +++ b/default_field.m Sun Dec 23 14:39:31 2018 +0100 @@ -1,7 +1,7 @@ function default_field(s, f, val) - if isfield(s,f) + if isfield(s,f) && ~isempty(s.(f)) return end s.(f) = val; assignin('caller', inputname(1),s); -end \ No newline at end of file +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/hgRevision.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,8 @@ +% Returns the short mercurial revision Id. +% ok is false if there are uncommited changes. +function [revId, ok] = hgRevision() + [~, s] = system('hg id -i'); + revId = strtrim(s); + + ok = s(end) ~= '+'; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/structArray.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,16 @@ +% % Usage example: +% c = structArray({'a','b'}, { +% 1, 2; +% 3, 4; +% }); + +function c = structArray(fields, values) + assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal'); + c = struct(); + + for i = 1:size(values, 1) + for j = 1:length(fields) + c(i).(fields{j}) = values{i,j}; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/structCellArray.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,16 @@ +% % Usage example: +% c = structCellArray({'a','b'}, { +% 1, 2; +% 3, 4; +% }); + +function c = structCellArray(fields, values) + assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal'); + c = cell(1, size(values, 1)); + + for i = 1:size(values, 1) + for j = 1:length(fields) + c{i}.(fields{j}) = values{i,j}; + end + end +end