Mercurial > repos > public > sbplib
changeset 752:0be9b4d6737b feature/utux2D
Merge with feature/interpolation to get new interpolation operators.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 22 May 2018 13:31:09 -0700 |
| parents | f4595f14d696 (diff) 005a8d071da3 (current diff) |
| children | f891758ad7a4 703183ed8c8b |
| files | |
| diffstat | 38 files changed, 2498 insertions(+), 217 deletions(-) [+] |
line wrap: on
line diff
diff -r 005a8d071da3 -r 0be9b4d6737b +blockmatrix/toMatrix.m --- a/+blockmatrix/toMatrix.m Tue May 22 13:29:47 2018 -0700 +++ b/+blockmatrix/toMatrix.m Tue May 22 13:31:09 2018 -0700 @@ -12,12 +12,9 @@ A = sparse(N,M); - n_ind = [0 cumsum(n)]; - m_ind = [0 cumsum(m)]; - for i = 1:size(bm,1) for j = 1:size(bm,2) - if(isempty(bm{i,j})) + if isempty(bm{i,j}) bm{i,j} = sparse(n(i),m(j)); end end
diff -r 005a8d071da3 -r 0be9b4d6737b +grid/Cartesian.m --- a/+grid/Cartesian.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/Cartesian.m Tue May 22 13:31:09 2018 -0700 @@ -5,6 +5,7 @@ m % Number of points in each direction x % Cell array of vectors with node placement for each dimension. h % Spacing/Scaling + lim % Cell array of left and right boundaries for each dimension. end % General d dimensional grid with n points @@ -27,6 +28,7 @@ end obj.h = []; + obj.lim = []; end % n returns the number of points in the grid function o = N(obj)
diff -r 005a8d071da3 -r 0be9b4d6737b +grid/TODO.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/TODO.txt Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,1 @@ +% TODO: Rename grid package. name conflicts with built in function
diff -r 005a8d071da3 -r 0be9b4d6737b +grid/evalOn.m --- a/+grid/evalOn.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/evalOn.m Tue May 22 13:31:09 2018 -0700 @@ -7,68 +7,34 @@ function gf = evalOn(g, func) if ~isa(func, 'function_handle') % We should have a constant. - if size(func,2) ~= 1 - error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') - end + assert(size(func,2) == 1,'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector'); gf = repmat(func,[g.N, 1]); return end % func should now be a function_handle + assert(g.D == nargin(func),'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') - if g.D ~= nargin(func) - g.D - nargin(func) - error('grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') - end + x = num2cell(g.points(),1); + k = numberOfComponents(func); + gf = func(x{:}); + gf = reorderComponents(gf, k); +end - % Get coordinates - x = g.points(); +% Find the number of vector components of func +function k = numberOfComponents(func) + x0 = num2cell(ones(1,nargin(func))); + f0 = func(x0{:}); + assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector'); + k = length(f0); +end - % Find the number of components - if size(x,1) ~= 0 - x0 = x(1,:); - else - x0 = num2cell(ones(1,size(x,2))); +% Reorder the components of the function to sit together +function gf = reorderComponents(a, k) + N = length(a)/k; + gf = zeros(N*k, 1); + for i = 1:k + gf(i:k:end) = a((i-1)*N + 1 : i*N); end - - dim = length(x0); - % Evaluate f0 = func(x0(1),x0(2),...,x0(dim)); - if(dim == 1) - f0 = func(x0); - else - eval_str = 'f0 = func(x0(1)'; - for i = 2:dim - eval_str = [eval_str, sprintf(',x0(%d)',i)]; - end - eval_str = [eval_str, ');']; - eval(eval_str); - end - - % k = number of components - k = length(f0); - - if size(f0,2) ~= 1 - error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') - end - - % Evaluate gf = func(x(:,1),x(:,2),...,x(:,dim)); - if(dim == 1) - gf = func(x); - else - eval_str = 'gf = func(x(:,1)'; - for i = 2:dim - eval_str = [eval_str, sprintf(',x(:,%d)',i)]; - end - eval_str = [eval_str, ');']; - eval(eval_str); - end - - % Reorganize gf - gf_temp = gf; - gf = zeros(g.N*k, 1); - for i = 1:k - gf(i:k:end) = gf_temp((i-1)*g.N + 1 : i*g.N); - end -end \ No newline at end of file +end
diff -r 005a8d071da3 -r 0be9b4d6737b +grid/evalOnTest.m --- a/+grid/evalOnTest.m Tue May 22 13:29:47 2018 -0700 +++ b/+grid/evalOnTest.m Tue May 22 13:31:09 2018 -0700 @@ -31,7 +31,7 @@ cases = { {getTestGrid('1d'), @(x,y)x-y}, {getTestGrid('2d'), @(x)x }, - } + }; for i = 1:length(cases) g = cases{i}{1}; @@ -111,9 +111,9 @@ function testInputErrorVectorValued(testCase) - in = { + in = { [1,2,3], - @(x,y)[x,-y]; + @(x,y)[x,-y], }; g = getTestGrid('2d');
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/+domain/Circle.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Circle.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,98 @@ +classdef Circle < multiblock.DefCurvilinear + properties + r, c + + hs + r_arc + omega + end + + methods + function obj = Circle(r, c, hs) + default_arg('r', 1); + default_arg('c', [0; 0]); + default_arg('hs', 0.435); + + + % alpha = 0.75; + % hs = alpha*r/sqrt(2); + + % Square should not be a square, it should be an arc. The arc radius + % is chosen so that the three angles of the meshes are all equal. + % This gives that the (half)arc opening angle of should be omega = pi/12 + omega = pi/12; + r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); % = hs* 1/sin(omega) + c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0]; + + cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]); + + c2 = cir(0); + c3 = cir(1); + + s1 = [-hs; -hs]; + s2 = [ hs; -hs]; + s3 = [ hs; hs]; + s4 = [-hs; hs]; + + sp2 = parametrization.Curve.line(s2,c2); + sp3 = parametrization.Curve.line(s3,c3); + + Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]); + Se2 = Se1.rotate(c,pi/2); + Se3 = Se2.rotate(c,pi/2); + Se4 = Se3.rotate(c,pi/2); + + + S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1); + + A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse); + B = A.rotate(c,1*pi/2).rotate_edges(-1); + C = A.rotate(c,2*pi/2).rotate_edges(-1); + D = A.rotate(c,3*pi/2).rotate_edges(0); + + blocks = {S,A,B,C,D}; + blocksNames = {'S','A','B','C','D'}; + + conn = cell(5,5); + conn{1,2} = {'e','w'}; + conn{1,3} = {'n','s'}; + conn{1,4} = {'w','s'}; + conn{1,5} = {'s','w'}; + + conn{2,3} = {'n','e'}; + conn{3,4} = {'w','e'}; + conn{4,5} = {'w','s'}; + conn{5,2} = {'n','s'}; + + boundaryGroups = struct(); + boundaryGroups.E = multiblock.BoundaryGroup({2,'e'}); + boundaryGroups.N = multiblock.BoundaryGroup({3,'n'}); + boundaryGroups.W = multiblock.BoundaryGroup({4,'n'}); + boundaryGroups.S = multiblock.BoundaryGroup({5,'e'}); + boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}}); + + obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames); + + obj.r = r; + obj.c = c; + obj.hs = hs; + obj.r_arc = r_arc; + obj.omega = omega; + end + + function ms = getGridSizes(obj, m) + m_S = m; + + % m_Radial + s = 2*obj.hs; + innerArc = obj.r_arc*obj.omega; + outerArc = obj.r*pi/2; + shortSpoke = obj.r-s/sqrt(2); + x = (1/(2-sqrt(3))-1)*obj.hs; + longSpoke = (obj.r+x)-obj.r_arc; + m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2); + + ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]}; + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/+domain/Rectangle.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Rectangle.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,188 @@ +classdef Rectangle < multiblock.Definition + properties + + blockTi % Transfinite interpolation objects used for plotting + xlims + ylims + blockNames % Cell array of block labels + nBlocks + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + + end + + + methods + % Creates a divided rectangle + % x and y are vectors of boundary and interface positions. + % blockNames: cell array of labels. The id is default. + function obj = Rectangle(x,y,blockNames) + default_arg('blockNames',[]); + + n = length(y)-1; % number of blocks in the y direction. + m = length(x)-1; % number of blocks in the x direction. + N = n*m; % number of blocks + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + if ~issorted(flip(y)) + error('The elements of y seem to be in the wrong order'); + end + + % Dimensions of blocks and number of points + blockTi = cell(N,1); + xlims = cell(N,1); + ylims = cell(N,1); + for i = 1:n + for j = 1:m + p1 = [x(j), y(i+1)]; + p2 = [x(j+1), y(i)]; + I = flat_index(m,j,i); + blockTi{I} = parametrization.Ti.rectangle(p1,p2); + xlims{I} = {x(j), x(j+1)}; + ylims{I} = {y(i+1), y(i)}; + end + end + + % Interface couplings + conn = cell(N,N); + for i = 1:n + for j = 1:m + I = flat_index(m,j,i); + if i < n + J = flat_index(m,j,i+1); + conn{I,J} = {'s','n'}; + end + + if j < m + J = flat_index(m,j+1,i); + conn{I,J} = {'e','w'}; + end + end + end + + % Block names (id number as default) + if isempty(blockNames) + obj.blockNames = cell(1, N); + for i = 1:N + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == N); + obj.blockNames = blockNames; + end + nBlocks = N; + + % Boundary groups + boundaryGroups = struct(); + nx = m; + ny = n; + E = cell(1,ny); + W = cell(1,ny); + S = cell(1,nx); + N = cell(1,nx); + for i = 1:ny + E_id = flat_index(m,nx,i); + W_id = flat_index(m,1,i); + E{i} = {E_id,'e'}; + W{i} = {W_id,'w'}; + end + for j = 1:nx + S_id = flat_index(m,j,ny); + N_id = flat_index(m,j,1); + S{j} = {S_id,'s'}; + N{j} = {N_id,'n'}; + end + boundaryGroups.E = multiblock.BoundaryGroup(E); + boundaryGroups.W = multiblock.BoundaryGroup(W); + boundaryGroups.S = multiblock.BoundaryGroup(S); + boundaryGroups.N = multiblock.BoundaryGroup(N); + boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]); + boundaryGroups.WS = multiblock.BoundaryGroup([W,S]); + boundaryGroups.WN = multiblock.BoundaryGroup([W,N]); + boundaryGroups.ES = multiblock.BoundaryGroup([E,S]); + boundaryGroups.EN = multiblock.BoundaryGroup([E,N]); + + obj.connections = conn; + obj.nBlocks = nBlocks; + obj.boundaryGroups = boundaryGroups; + obj.blockTi = blockTi; + obj.xlims = xlims; + obj.ylims = ylims; + + end + + + % Returns a multiblock.Grid given some parameters + % ms: cell array of [mx, my] vectors + % For same [mx, my] in every block, just input one vector. + function g = getGrid(obj, ms, varargin) + + default_arg('ms',[21,21]) + + % Extend ms if input is a single vector + if (numel(ms) == 2) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockTi{i}.show(2,2); + end + axis equal + return + end + + if gridLines + m = 10; + for i = 1:obj.nBlocks + obj.blockTi{i}.show(m,m); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockTi{i}, labels{i}); + end + + axis equal + end + + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + function ms = getGridSizes(obj, varargin) + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/Def.m --- a/+multiblock/Def.m Tue May 22 13:29:47 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,101 +0,0 @@ -classdef Def - properties - nBlocks - blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation - blockNames - connections % Cell array specifying connections between blocks - boundaryGroups % Structure of boundaryGroups - end - - methods - % Defines a multiblock setup for transfinite interpolation blocks - % TODO: How to bring in plotting of points? - function obj = Def(blockMaps, connections, boundaryGroups, blockNames) - default_arg('boundaryGroups', struct()); - default_arg('blockNames',{}); - - nBlocks = length(blockMaps); - - obj.nBlocks = nBlocks; - - obj.blockMaps = blockMaps; - - assert(all(size(connections) == [nBlocks, nBlocks])); - obj.connections = connections; - - - if isempty(blockNames) - obj.blockNames = cell(1, nBlocks); - for i = 1:length(blockMaps) - obj.blockNames{i} = sprintf('%d', i); - end - else - assert(length(blockNames) == nBlocks); - obj.blockNames = blockNames; - end - - obj.boundaryGroups = boundaryGroups; - end - - function g = getGrid(obj, varargin) - ms = obj.getGridSizes(varargin{:}); - - grids = cell(1, obj.nBlocks); - for i = 1:obj.nBlocks - grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); - end - - g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); - end - - function show(obj, label, gridLines, varargin) - default_arg('label', 'name') - default_arg('gridLines', false); - - if isempty('label') && ~gridLines - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(2,2); - end - axis equal - return - end - - if gridLines - ms = obj.getGridSizes(varargin{:}); - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); - end - end - - - switch label - case 'name' - labels = obj.blockNames; - case 'id' - labels = {}; - for i = 1:obj.nBlocks - labels{i} = num2str(i); - end - otherwise - axis equal - return - end - - for i = 1:obj.nBlocks - parametrization.Ti.label(obj.blockMaps{i}, labels{i}); - end - - axis equal - end - end - - methods (Abstract) - % Returns the grid size of each block in a cell array - % The input parameters are determined by the subclass - ms = getGridSizes(obj, varargin) - % end - end - -end - -
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/DefCurvilinear.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DefCurvilinear.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,101 @@ +classdef DefCurvilinear < multiblock.Definition + properties + nBlocks + blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation + blockNames + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + end + + methods + % Defines a multiblock setup for transfinite interpolation blocks + % TODO: How to bring in plotting of points? + function obj = DefCurvilinear(blockMaps, connections, boundaryGroups, blockNames) + default_arg('boundaryGroups', struct()); + default_arg('blockNames',{}); + + nBlocks = length(blockMaps); + + obj.nBlocks = nBlocks; + + obj.blockMaps = blockMaps; + + assert(all(size(connections) == [nBlocks, nBlocks])); + obj.connections = connections; + + + if isempty(blockNames) + obj.blockNames = cell(1, nBlocks); + for i = 1:length(blockMaps) + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == nBlocks); + obj.blockNames = blockNames; + end + + obj.boundaryGroups = boundaryGroups; + end + + function g = getGrid(obj, varargin) + ms = obj.getGridSizes(varargin{:}); + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(2,2); + end + axis equal + return + end + + if gridLines + ms = obj.getGridSizes(varargin{:}); + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockMaps{i}, labels{i}); + end + + axis equal + end + end + + methods (Abstract) + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + ms = getGridSizes(obj, varargin) + % end + end + +end + +
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/Definition.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Definition.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,11 @@ +classdef Definition + methods (Abstract) + + % Returns a multiblock.Grid given some parameters + g = getGrid(obj, varargin) + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + show(obj, label, gridLines, varargin) + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +multiblock/DiffOp.m --- a/+multiblock/DiffOp.m Tue May 22 13:29:47 2018 -0700 +++ b/+multiblock/DiffOp.m Tue May 22 13:31:09 2018 -0700 @@ -194,6 +194,7 @@ p{I} = blockPenalty; penalty = blockmatrix.toMatrix(p); else + % TODO: used by beam equation, should be eliminated. SHould only set one BC per call for i = 1:length(blockPenalty) div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector p = blockmatrix.zero(div);
diff -r 005a8d071da3 -r 0be9b4d6737b +noname/Animation.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/Animation.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,75 @@ +classdef Animation < handle + properties + timeStepper + representationMaker + updaters + end + + % add input validation + + methods + function obj = Animation(timeStepper, representationMaker, updaters); + obj.timeStepper = timeStepper; + obj.updaters = updaters; + obj.representationMaker = representationMaker; + end + + function update(obj, r) + for i = 1:length(obj.updaters) + obj.updaters{i}(r); + end + drawnow + end + + function run(obj, tEnd, timeModifier, do_pause) + default_arg('do_pause', false) + + function next_t = G(next_t) + obj.timeStepper.evolve(next_t); + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + if do_pause + pause + end + end + + anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier); + end + + function step(obj, tEnd, do_pause) + default_arg('do_pause', false) + + while obj.timeStepper.t < tEnd + obj.timeStepper.step(); + + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + % TODO: Make it never go faster than a certain fram rate + + if do_pause + pause + end + end + end + + function saveMovie(obj, tEnd, timeModifier, figureHandle, dirname) + save_frame = anim.setup_fig_mov(figureHandle, dirname); + + function next_t = G(next_t) + obj.timeStepper.evolve(next_t); + r = obj.representationMaker(obj.timeStepper); + obj.update(r); + + save_frame(); + end + + fprintf('Generating and saving frames to: ..\n') + anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier); + fprintf('Generating movies...\n') + cmd = sprintf('bash %s/+anim/make_movie.sh %s', sbplibLocation(),dirname); + system(cmd); + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +noname/calculateErrors.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/calculateErrors.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,40 @@ +% [discr, trueSolution] = schemeFactory(m) +% where trueSolution should be a timeSnapshot of the true solution a time T +% T is the end time +% m are grid size parameters. +% N are number of timesteps to use for each gird size +% timeOpt are options for the timeStepper +function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt) + assertType(schemeFactory, 'function_handle'); + assertNumberOfArguments(schemeFactory, 1); + assertScalar(T); + assert(length(m) == length(N), 'Vectors m and N must have the same length'); + assertType(errorFun, 'function_handle'); + assertNumberOfArguments(errorFun, 2); + default_arg('timeOpt'); + + e = []; + for i = 1:length(m) + done = timeTask('m = %3d ', m(i)); + + [discr, trueSolution] = schemeFactory(m(i)); + + timeOpt.k = T/N(i); + ts = discr.getTimestepper(timeOpt); + ts.stepTo(N(i), true); + approxSolution = discr.getTimeSnapshot(ts); + + e(i) = errorFun(trueSolution, approxSolution); + + fprintf('e = %.4e', e(i)) + done() + end + fprintf('\n') +end + + +%% Example error function +% u_true = grid.evalOn(dr.grid, @(x,y)trueSolution(T,x,y)); +% err = u_true-u_false; +% e(i) = norm(err)/norm(u_true); +% % e(i) = sqrt(err'*d.H*d.J*err/(u_true'*d.H*d.J*u_true));
diff -r 005a8d071da3 -r 0be9b4d6737b +parametrization/Ti.m --- a/+parametrization/Ti.m Tue May 22 13:29:47 2018 -0700 +++ b/+parametrization/Ti.m Tue May 22 13:31:09 2018 -0700 @@ -21,16 +21,29 @@ D = g4(0); function o = S_fun(u,v) + if isrow(u) && isrow(v) + flipped = false; + else + flipped = true; + u = u'; + v = v'; + end + x1 = g1(u); x2 = g2(v); x3 = g3(1-u); x4 = g4(1-v); + o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ... - -((1-u)*(1-v).*A(1,:) + u*(1-v).*B(1,:) + u*v.*C(1,:) + (1-u)*v.*D(1,:)); + -((1-u).*(1-v).*A(1,:) + u.*(1-v).*B(1,:) + u.*v.*C(1,:) + (1-u).*v.*D(1,:)); o2 = (1-v).*x1(2,:) + u.*x2(2,:) + v.*x3(2,:) + (1-u).*x4(2,:) ... - -((1-u)*(1-v).*A(2,:) + u*(1-v).*B(2,:) + u*v.*C(2,:) + (1-u)*v.*D(2,:)); + -((1-u).*(1-v).*A(2,:) + u.*(1-v).*B(2,:) + u.*v.*C(2,:) + (1-u).*v.*D(2,:)); - o = [o1;o2]; + if ~flipped + o = [o1;o2]; + else + o = [o1'; o2']; + end end obj.S = @S_fun;
diff -r 005a8d071da3 -r 0be9b4d6737b +parametrization/TiTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/TiTest.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,52 @@ +function tests = TiTest() + tests = functiontests(localfunctions); +end + +function testScalarInput(testCase) + ti = getMinimumTi(); + + cases = { + % {u, v, out}, + {0, 0, [1; 2]}, + {0, 1, [1; 4]}, + {1, 0, [3; 2]}, + {1, 1, [3; 4]}, + {0.5, 0.5, [2; 3]}, + }; + + for i = 1:length(cases) + u = cases{i}{1}; + v = cases{i}{2}; + expected = cases{i}{3}; + + testCase.verifyEqual(ti.S(u,v), expected, sprintf('Case: %d',i)); + end +end + +function testRowVectorInput(testCase) + ti = getMinimumTi(); + + u = [0, 0.5, 1]; + v = [0, 0, 0.5]; + expected = [ + 1, 2, 3; + 2, 2, 3; + ]; + + testCase.verifyEqual(ti.S(u,v), expected); +end + +function testColumnvectorInput(testCase) + ti = getMinimumTi(); + + u = [0; 0.5; 1]; + v = [0; 0; 0.5]; + expected = [1; 2; 3; 2; 2; 3]; + + testCase.verifyEqual(ti.S(u,v), expected); +end + + +function ti = getMinimumTi() + ti = parametrization.Ti.rectangle([1; 2], [3; 4]); +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b +sbp/+implementations/d2_variable_periodic_2.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_2.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,50 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<3) + error(['Operator requires at least ' num2str(3) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + vals = [Mm1,M0,Mp1]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b +sbp/+implementations/d2_variable_periodic_4.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_4.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,57 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<5) + error(['Operator requires at least ' num2str(5) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + S = d1_l*d1_l' + d1_r*d1_r'; + + % D1 operator + stencil = [1/12 -2/3 0 2/3 -1/12]; + diags = -2:2; + Q = stripeMatrixPeriodic(stencil, diags, m); + D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); + + + scheme_width = 5; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + % Note: these coefficients are for -M. + Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r); + Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1); + M0 = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2); + Mp1 = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2); + Mp2 = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2); + + vals = -[Mm2,Mm1,M0,Mp1,Mp2]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + + end + D2 = @D2_fun; +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b +sbp/+implementations/d2_variable_periodic_6.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_6.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,58 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<7) + error(['Operator requires at least ' num2str(7) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + + % D1 operator + diags = -3:3; + stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + % D2 operator + scheme_width = 7; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm3 = c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r); + Mm2 = c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2); + Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1); + M0 = c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r); + Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1); + Mp2 = c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2); + Mp3 = c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3); + + vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; + + +end
diff -r 005a8d071da3 -r 0be9b4d6737b +sbp/D2Variable.m --- a/+sbp/D2Variable.m Tue May 22 13:29:47 2018 -0700 +++ b/+sbp/D2Variable.m Tue May 22 13:31:09 2018 -0700 @@ -26,22 +26,39 @@ obj.x = linspace(x_l,x_r,m)'; switch order + + case 6 + + [obj.H, obj.HI, obj.D1, obj.D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + obj.d1_l, obj.d1_r] = ... + sbp.implementations.d4_variable_6(m, obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + case 4 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_4(m,obj.h); obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R case 2 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_2(m,obj.h); obj.borrowing.M.d1 = 0.3636363636; % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R otherwise error('Invalid operator order %d.',order); end - + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, obj.m = m; obj.M = []; end
diff -r 005a8d071da3 -r 0be9b4d6737b +sbp/D2VariablePeriodic.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariablePeriodic.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,71 @@ +classdef D2VariablePeriodic < sbp.OpSet + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariablePeriodic(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/m; + x = linspace(x_l,x_r,m+1)'; + obj.x = x(1:end-1); + + switch order + + case 6 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_6(m,obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + + case 4 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_4(m,obj.h); + obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_2(m,obj.h); + obj.borrowing.M.d1 = 0.3636363636; + % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + + obj.m = m; + obj.M = []; + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + + +end
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Beam.m --- a/+scheme/Beam.m Tue May 22 13:29:47 2018 -0700 +++ b/+scheme/Beam.m Tue May 22 13:31:09 2018 -0700 @@ -126,6 +126,44 @@ penalty{1} = -obj.Hi*tau; penalty{1} = -obj.Hi*sig; + case 'e' + alpha = obj.alpha; + tuning = 1.1; + + tau1 = tuning * alpha/delt; + tau4 = s*alpha; + + tau = tau1*e+tau4*d3; + + closure = obj.Hi*tau*e'; + penalty = -obj.Hi*tau; + case 'd1' + alpha = obj.alpha; + + tuning = 1.1; + + sig2 = tuning * alpha/gamm; + sig3 = -s*alpha; + + sig = sig2*d1+sig3*d2; + + closure = obj.Hi*sig*d1'; + penalty = -obj.Hi*sig; + + case 'd2' + a = obj.alpha; + + tau = s*a*d1; + + closure = obj.Hi*tau*d2'; + penalty = -obj.Hi*tau; + case 'd3' + a = obj.alpha; + + sig = -s*a*e; + + closure = obj.Hi*sig*d3'; + penalty = -obj.Hi*sig; otherwise % Unknown, boundary condition error('No such boundary condition: type = %s',type);
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Elastic2dVariable.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dVariable.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,420 @@ +classdef Elastic2dVariable < scheme.Scheme + +% Discretizes the elastic wave equation: +% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i +% 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 matrices for varible coefficients + LAMBDA % Variable coefficient, related to dilation + MU % Shear modulus, variable coefficient + RHO, RHOi % Density, variable + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_lambda + D2_mu + + % Traction operators used for BC + T_l, T_r + tau_l, tau_r + + H, Hi % Inner products + phi % Borrowing constant for (d1 - e^T*D1) from R + gamma % Borrowing constant for d1 from M + H11 % First element of H + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + E % E{i}^T picks out component i + + H_boundary % Boundary inner products + + % Kroneckered norms and coefficients + RHOi_kron + Hi_kron + end + + methods + + function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('lambda_fun', @(x,y) 0*x+1); + default_arg('mu_fun', @(x,y) 0*x+1); + default_arg('rho_fun', @(x,y) 0*x+1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + + lambda = grid.evalOn(g, lambda_fun); + mu = grid.evalOn(g, mu_fun); + rho = grid.evalOn(g, rho_fun); + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + lim = g.lim; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + end + + % Borrowing constants + for i = 1:dim + beta = ops{i}.borrowing.R.delta_D; + obj.H11{i} = ops{i}.borrowing.H11; + obj.phi{i} = beta/obj.H11{i}; + obj.gamma{i} = ops{i}.borrowing.M.d1; + 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 + + %====== Assemble full operators ======== + LAMBDA = spdiag(lambda); + obj.LAMBDA = LAMBDA; + MU = spdiag(mu); + obj.MU = MU; + RHO = spdiag(rho); + obj.RHO = RHO; + obj.RHOi = inv(RHO); + + obj.D1 = cell(dim,1); + obj.D2_lambda = cell(dim,1); + obj.D2_mu = 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 + for i = 1:dim + obj.D2_lambda{i} = sparse(m_tot); + obj.D2_mu{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_lambda = D2{1}(lambda(ind(:,i))); + D_mu = D2{1}(mu(ind(:,i))); + + p = ind(:,i); + obj.D2_lambda{1}(p,p) = D_lambda; + obj.D2_mu{1}(p,p) = D_mu; + end + + for i = 1:m(1) + D_lambda = D2{2}(lambda(ind(i,:))); + D_mu = D2{2}(mu(ind(i,:))); + + p = ind(i,:); + obj.D2_lambda{2}(p,p) = D_lambda; + obj.D2_mu{2}(p,p) = D_mu; + end + + % 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}; + + % E{i}^T picks out component i. + E = cell(dim,1); + I = speye(m_tot,m_tot); + for i = 1:dim + e = sparse(dim,1); + e(i) = 1; + E{i} = kron(I,e); + end + obj.E = E; + + % Differentiation matrix D (without SAT) + D2_lambda = obj.D2_lambda; + D2_mu = obj.D2_mu; + D1 = obj.D1; + D = sparse(dim*m_tot,dim*m_tot); + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + for i = 1:dim + for j = 1:dim + D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +... + db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ... + ); + D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +... + db(i,j)*D1{j}*MU*D1{i}*E{j}' + ... + D2_mu{j}*E{i}' ... + ); + end + end + obj.D = D; + %=========================================% + + % Numerical traction operators for BC. + % Because d1 =/= e0^T*D1, the numerical tractions are different + % at every boundary. + T_l = cell(dim,1); + T_r = cell(dim,1); + tau_l = cell(dim,1); + tau_r = cell(dim,1); + % tau^{j}_i = sum_k T^{j}_{ik} u_k + + d1_l = obj.d1_l; + d1_r = obj.d1_r; + e_l = obj.e_l; + e_r = obj.e_r; + D1 = obj.D1; + + % Loop over boundaries + for j = 1:dim + T_l{j} = cell(dim,dim); + T_r{j} = cell(dim,dim); + tau_l{j} = cell(dim,1); + tau_r{j} = cell(dim,1); + + % Loop over components + for i = 1:dim + tau_l{j}{i} = sparse(m_tot,dim*m_tot); + tau_r{j}{i} = sparse(m_tot,dim*m_tot); + for k = 1:dim + T_l{j}{i,k} = ... + -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... + -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... + -d(i,k)*MU*e_l{j}*d1_l{j}'; + + T_r{j}{i,k} = ... + d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... + +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... + +d(i,k)*MU*e_r{j}*d1_r{j}'; + + tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; + tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; + end + + end + end + obj.T_l = T_l; + obj.T_r = T_r; + obj.tau_l = tau_l; + obj.tau_r = tau_r; + + % Kroneckered norms and coefficients + I_dim = speye(dim); + obj.RHOi_kron = kron(obj.RHOi, I_dim); + obj.Hi_kron = kron(obj.Hi, I_dim); + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + + 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 cell array of strings specifying the type of boundary condition for each component. + % 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',{'free','free'}); + 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; + tau = obj.tau_r{j}; + T = obj.T_r{j}; + case -1 + e = obj.e_l; + d = obj.d1_l; + tau = obj.tau_l{j}; + T = obj.T_l{j}; + end + + E = obj.E; + Hi = obj.Hi; + H_gamma = obj.H_boundary{j}; + LAMBDA = obj.LAMBDA; + MU = obj.MU; + RHOi = obj.RHOi; + + dim = obj.dim; + m_tot = obj.grid.N(); + + RHOi_kron = obj.RHOi_kron; + Hi_kron = obj.Hi_kron; + + % Preallocate + closure = sparse(dim*m_tot, dim*m_tot); + penalty = cell(dim,1); + for k = 1:dim + penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j)); + end + + % Loop over components that we (potentially) have different BC on + for k = 1:dim + switch type{k} + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + + tuning = 1.2; + phi = obj.phi{j}; + h = obj.h(j); + h11 = obj.H11{j}*h; + gamma = obj.gamma{j}; + + a_lambda = dim/h11 + 1/(h11*phi); + a_mu_i = 2/(gamma*h); + a_mu_ij = 2/h11 + 1/(h11*phi); + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU ); + + % Loop over components that Dirichlet penalties end up on + for i = 1:dim + C = T{k,i}; + A = -d(i,k)*alpha(i,j); + B = A + C; + closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); + penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; + end + + % Free boundary condition + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); + penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + 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 + tuning = 1.2; + % tuning = 20.2; + error('Interface not implemented'); + 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 coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, 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 op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Heat2dVariable.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dVariable.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,262 @@ +classdef Heat2dVariable < scheme.Scheme + +% Discretizes the Laplacian with variable coefficent, +% In the Heat equation way (i.e., the discretization matrix is not necessarily +% symmetric) +% u_t = div * (kappa * grad 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 + KAPPA % Variable coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_kappa + + H, Hi % Inner products + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + + H_boundary % Boundary inner products + + end + + methods + + function obj = Heat2dVariable(g ,order, kappa_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('kappa_fun', @(x,y) 0*x+1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + + kappa = grid.evalOn(g, kappa_fun); + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + lim = g.lim; + + % 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 + + %====== Assemble full operators ======== + KAPPA = spdiag(kappa); + obj.KAPPA = KAPPA; + + obj.D1 = cell(dim,1); + obj.D2_kappa = 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 + for i = 1:dim + obj.D2_kappa{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_kappa = D2{1}(kappa(ind(:,i))); + p = ind(:,i); + obj.D2_kappa{1}(p,p) = D_kappa; + end + + for i = 1:m(1) + D_kappa = D2{2}(kappa(ind(i,:))); + p = ind(i,:); + obj.D2_kappa{2}(p,p) = D_kappa; + end + + % 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_kappa = obj.D2_kappa; + D1 = obj.D1; + D = sparse(m_tot,m_tot); + for i = 1:dim + D = D + D2_kappa{i}; + end + obj.D = D; + %=========================================% + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + + 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}; + KAPPA = obj.KAPPA; + kappa_gamma = e{j}'*KAPPA*e{j}; + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); + penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; + + % Free boundary condition + case {'N','n','neumann','Neumann'} + closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; + + % 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 + error('Interface not implemented'); + 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 coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, 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 op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/LaplaceCurvilinear.m --- a/+scheme/LaplaceCurvilinear.m Tue May 22 13:29:47 2018 -0700 +++ b/+scheme/LaplaceCurvilinear.m Tue May 22 13:31:09 2018 -0700 @@ -38,22 +38,29 @@ du_n, dv_n gamm_u, gamm_v lambda + + interpolation_type end methods % Implements a*div(b*grad(u)) as a SBP scheme % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) - function obj = LaplaceCurvilinear(g ,order, a, b, opSet) + function obj = LaplaceCurvilinear(g ,order, a, b, opSet, interpolation_type) default_arg('opSet',@sbp.D2Variable); default_arg('a', 1); default_arg('b', 1); + default_arg('interpolation_type','AWW'); if b ~=1 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); @@ -209,6 +216,7 @@ obj.h = [h_u h_v]; obj.order = order; obj.grid = g; + obj.interpolation_type = interpolation_type; obj.a = a; obj.b = b; @@ -269,13 +277,70 @@ end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + [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); + Hi = obj.Hi; + a = obj.a; + + m_u = length(H_b_u); + m_v = length(H_b_v); + + grid_ratio = m_u/m_v; + if grid_ratio ~= 1 + + [ms, index] = sort([m_u, m_v]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_v2u_good = interpOpSet.IF2C; + I_v2u_bad = interpOpSet.IF2C; + I_u2v_good = interpOpSet.IC2F; + I_u2v_bad = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_v2u_good = interpOpSet.IC2F; + I_v2u_bad = interpOpSet.IC2F; + I_u2v_good = interpOpSet.IF2C; + I_u2v_bad = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_v2u_good = interpOpSetF2C.IF2C; + I_v2u_bad = interpOpSetC2F.IF2C; + I_u2v_good = interpOpSetC2F.IC2F; + I_u2v_bad = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_v2u_good = interpOpSetC2F.IC2F; + I_v2u_bad = interpOpSetF2C.IC2F; + I_u2v_good = interpOpSetF2C.IF2C; + I_u2v_bad = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_v2u_good = speye(m_u,m_u); + I_v2u_bad = speye(m_u,m_u); + I_u2v_good = speye(m_u,m_u); + I_u2v_bad = speye(m_u,m_u); + end + % 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); - u = obj; v = neighbour_scheme; @@ -284,18 +349,24 @@ 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; - tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - tau1 = tuning * spdiag(tau1); - tau2 = 1/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; - sig1 = -1/2; - sig2 = 0; + closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... + a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_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'; - tau = (e_u*tau1 + tau2*d_u)*H_b_u; - sig = (sig1*e_u + sig2*d_u)*H_b_u; + penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ... + -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ... + -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ... + -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v'; + - closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); - 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.
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Schrodinger2d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2d.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,374 @@ +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 + + interpolation_type % MC or AWW + + end + + methods + + function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type) + default_arg('interpolation_type','AWW'); + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('a',1); + dim = 2; + + assert(isa(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}; + obj.interpolation_type = interpolation_type; + + 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 + + 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 + % Get neighbour boundary operator + + [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary); + [coord, n] = get_boundary_number(obj, boundary); + switch n_nei + case 1 + % North or east boundary + e_neighbour = neighbour_scheme.e_r; + d_neighbour = neighbour_scheme.d1_r; + case -1 + % South or west boundary + e_neighbour = neighbour_scheme.e_l; + d_neighbour = neighbour_scheme.d1_l; + end + + e_neighbour = e_neighbour{coord_nei}; + d_neighbour = d_neighbour{coord_nei}; + H_gamma = obj.H_boundary{coord}; + Hi = obj.Hi; + a = obj.a; + + switch coord_nei + case 1 + m_neighbour = neighbour_scheme.m(2); + case 2 + m_neighbour = neighbour_scheme.m(1); + end + + switch coord + case 1 + m = obj.m(2); + case 2 + m = obj.m(1); + end + + switch n + case 1 + % North or east boundary + e = obj.e_r; + d = obj.d1_r; + case -1 + % South or west boundary + e = obj.e_l; + d = obj.d1_l; + end + e = e{coord}; + d = d{coord}; + + Hi = obj.Hi; + sigma = -n*1i*a/2; + tau = -n*(1i*a)'/2; + + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_e = interpOpSet.IF2C; + I_neighbour2local_d = interpOpSet.IF2C; + I_local2neighbour_e = interpOpSet.IC2F; + I_local2neighbour_d = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_neighbour2local_e = interpOpSet.IC2F; + I_neighbour2local_d = interpOpSet.IC2F; + I_local2neighbour_e = interpOpSet.IF2C; + I_local2neighbour_d = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_e = interpOpSetF2C.IF2C; + I_neighbour2local_d = interpOpSetC2F.IF2C; + I_local2neighbour_e = interpOpSetC2F.IC2F; + I_local2neighbour_d = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_e = interpOpSetC2F.IC2F; + I_neighbour2local_d = interpOpSetF2C.IC2F; + I_local2neighbour_e = interpOpSetF2C.IF2C; + I_local2neighbour_d = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_e = speye(m,m); + I_neighbour2local_d = speye(m,m); + I_local2neighbour_e = speye(m,m); + I_local2neighbour_d = speye(m,m); + end + + closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; + penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ... + -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; + + 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 coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, 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 op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operator = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/TODO.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/TODO.txt Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,1 @@ +% TODO: Rename package and abstract class to diffOp
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Utux.m --- a/+scheme/Utux.m Tue May 22 13:29:47 2018 -0700 +++ b/+scheme/Utux.m Tue May 22 13:31:09 2018 -0700 @@ -2,7 +2,7 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing - x % Grid + grid % Grid order % Order accuracy for the approximation H % Discrete norm @@ -17,41 +17,40 @@ methods - function obj = Utux(m,xlim,order,operator) - default_arg('a',1); - - %Old operators - % [x, h] = util.get_grid(xlim{:},m); - %ops = sbp.Ordinary(m,h,order); - - + function obj = Utux(g ,order, operator) + default_arg('operator','Standard'); + + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + switch operator - case 'NonEquidistant' - ops = sbp.D1Nonequidistant(m,xlim,order); - obj.D1 = ops.D1; +% case 'NonEquidistant' +% ops = sbp.D1Nonequidistant(m,xlim,order); +% obj.D1 = ops.D1; case 'Standard' ops = sbp.D2Standard(m,xlim,order); obj.D1 = ops.D1; - case 'Upwind' - ops = sbp.D1Upwind(m,xlim,order); - obj.D1 = ops.Dm; +% case 'Upwind' +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; otherwise error('Unvalid operator') end - obj.x=ops.x; + + obj.grid = g; - obj.H = ops.H; obj.Hi = ops.HI; obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.D=obj.D1; + obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -61,17 +60,27 @@ % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj,boundary,type,data) - default_arg('type','neumann'); - default_arg('data',0); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; + penalty = -obj.Hi*tau; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + end function N = size(obj) @@ -81,9 +90,9 @@ end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/Utux2D.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Utux2D.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,286 @@ +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 + + % String, type of interface coupling + % Default: 'upwind' + % Other: 'centered' + coupling_type + + % String, type of interpolation operators + % Default: 'AWW' (Almquist Wang Werpers) + % Other: 'MC' (Mattsson Carpenter) + interpolation_type + + + % Cell array, damping on upwstream and downstream sides. + interpolation_damping + + end + + + methods + function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) + + default_arg('interpolation_damping',{0,0}); + default_arg('interpolation_type','AWW'); + default_arg('coupling_type','upwind'); + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + + assert(isa(g, 'grid.Cartesian')) + 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.coupling_type = coupling_type; + obj.interpolation_type = interpolation_type; + obj.interpolation_damping = interpolation_damping; + obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); + + end + % Closure functions return the opertors applied to the own domain to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); + + sigma = -1; % Scalar penalty parameter + switch boundary + case {'w','W','west','West'} + tau = sigma*obj.a{1}*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + + case {'s','S','south','South'} + tau = sigma*obj.a{2}*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + end + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + m_neighbour = neighbour_scheme.m(2); + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + m_neighbour = neighbour_scheme.m(2); + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + m_neighbour = neighbour_scheme.m(1); + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + m_neighbour = neighbour_scheme.m(1); + end + + switch obj.coupling_type + + % Upwind coupling (energy dissipation) + case 'upwind' + sigma_ds = -1; %"Downstream" penalty + sigma_us = 0; %"Upstream" penalty + + % Energy-preserving coupling (no energy dissipation) + case 'centered' + sigma_ds = -1/2; %"Downstream" penalty + sigma_us = 1/2; %"Upstream" penalty + + otherwise + error(['Interface coupling type ' coupling_type ' is not available.']) + end + + % Check grid ratio for interpolation + switch boundary + case {'w','W','west','West','e','E','east','East'} + m = obj.m(2); + case {'s','S','south','South','n','N','north','North'} + m = obj.m(1); + end + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_us = interpOpSet.IF2C; + I_neighbour2local_ds = interpOpSet.IF2C; + I_local2neighbour_us = interpOpSet.IC2F; + I_local2neighbour_ds = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_neighbour2local_us = interpOpSet.IC2F; + I_neighbour2local_ds = interpOpSet.IC2F; + I_local2neighbour_us = interpOpSet.IF2C; + I_local2neighbour_ds = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_us = interpOpSetC2F.IF2C; + I_neighbour2local_ds = interpOpSetF2C.IF2C; + I_local2neighbour_us = interpOpSetC2F.IC2F; + I_local2neighbour_ds = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_us = interpOpSetF2C.IC2F; + I_neighbour2local_ds = interpOpSetC2F.IC2F; + I_local2neighbour_us = interpOpSetF2C.IF2C; + I_local2neighbour_ds = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_us = speye(m,m); + I_neighbour2local_ds = speye(m,m); + end + + int_damp_us = obj.interpolation_damping{1}; + int_damp_ds = obj.interpolation_damping{2}; + + I = speye(m,m); + I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us; + I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds; + + + switch boundary + case {'w','W','west','West'} + tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a{1}... + *obj.e_w*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; + case {'e','E','east','East'} + tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; + closure = obj.Hi*tau*obj.e_e'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a{1}... + *obj.e_e*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; + case {'s','S','south','South'} + tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a{2}... + *obj.e_s*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; + case {'n','N','north','North'} + tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x; + closure = obj.Hi*tau*obj.e_n'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a{2}... + *obj.e_n*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; + end + + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b +scheme/bcSetup.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/bcSetup.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,48 @@ +% 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 with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem +% Also takes S_sign which modifies the sign of S, [-1,1] +% Returns a closure matrix and a forcing function S +function [closure, S] = bcSetup(diffOp, bc, S_sign) + default_arg('S_sign', 1); + assertType(bc, 'cell'); + assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); + + + closure = spzeros(size(diffOp)); + penalties = {}; + dataFunctions = {}; + dataParams = {}; + + for i = 1:length(bc) + assertType(bc{i}, 'struct'); + [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); + closure = closure + localClosure; + + if isempty(bc{i}.data) + continue + end + assertType(bc{i}.data, 'function_handle'); + + coord = diffOp.grid.getBoundary(bc{i}.boundary); + assertNumberOfArguments(bc{i}.data, 1+size(coord,2)); + + penalties{end+1} = penalty; + dataFunctions{end+1} = bc{i}.data; + dataParams{end+1} = num2cell(coord ,1); + end + + O = spzeros(size(diffOp),1); + function v = S_fun(t) + v = O; + for i = 1:length(dataFunctions) + v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:}); + end + + v = S_sign * v; + end + S = @S_fun; +end
diff -r 005a8d071da3 -r 0be9b4d6737b +time/Timestepper.m --- a/+time/Timestepper.m Tue May 22 13:29:47 2018 -0700 +++ b/+time/Timestepper.m Tue May 22 13:31:09 2018 -0700 @@ -62,6 +62,7 @@ function [v, t] = stepTo(obj, n, progress_bar) + assertScalar(n); default_arg('progress_bar',false); [v, t] = obj.stepN(n-obj.n, progress_bar);
diff -r 005a8d071da3 -r 0be9b4d6737b assertScalar.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertScalar.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,5 @@ +function assertScalar(obj) + if ~isscalar(obj) + error('sbplib:assertScalar:notScalar', '"%s" must be scalar, found size "%s"', inputname(1), toString(size(obj))); + end +end
diff -r 005a8d071da3 -r 0be9b4d6737b diffSymfun.m --- a/diffSymfun.m Tue May 22 13:29:47 2018 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -% Differentiates a symbolic function like diff does, but keeps the function as a symfun -function g = diffSymfun(f, varargin) - assertType(f, 'symfun'); - - args = argnames(f); - g = symfun(diff(f,varargin{:}), args); -end
diff -r 005a8d071da3 -r 0be9b4d6737b kroneckerDelta.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/kroneckerDelta.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,6 @@ +function d = kroneckerDelta(i,j) + +d = 0; +if i==j + d = 1; +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b sbplibLocation.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sbplibLocation.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,4 @@ +function location = sbplibLocation() + scriptname = mfilename('fullpath'); + [location, ~, ~] = fileparts(scriptname); +end
diff -r 005a8d071da3 -r 0be9b4d6737b sbplibVersion.m --- a/sbplibVersion.m Tue May 22 13:29:47 2018 -0700 +++ b/sbplibVersion.m Tue May 22 13:31:09 2018 -0700 @@ -1,11 +1,10 @@ % Prints the version and location of the sbplib currently in use. function sbplibVersion() - scriptname = mfilename('fullpath'); - [folder,~,~] = fileparts(scriptname); + location = sbplibLocation(); name = 'sbplib (feature/grids)'; ver = '0.0.x'; fprintf('%s %s\n', name, ver); - fprintf('Running in:\n%s\n',folder); + fprintf('Running in:\n%s\n', location); end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b spdiagVariable.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagVariable.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,17 @@ +function A = spdiagVariable(a,i) + default_arg('i',0); + + if isrow(a) + a = a'; + end + + n = length(a)+abs(i); + + if i > 0 + a = [sparse(i,1); a]; + elseif i < 0 + a = [a; sparse(abs(i),1)]; + end + + A = spdiags(a,i,n,n); +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b spdiagsVariablePeriodic.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagsVariablePeriodic.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,42 @@ +function A = spdiagsVariablePeriodic(vals,diags) + % Creates an m x m periodic discretization matrix. + % vals - m x ndiags matrix of values + % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on + % vals that are not on main diagonal are going to spill over to + % off-diagonal corners. + + default_arg('diags',0); + + [m, ~] = size(vals); + + A = sparse(m,m); + + for i = 1:length(diags) + + d = diags(i); + a = vals(:,i); + + % Sub-diagonals + if d < 0 + a_bulk = a(1+abs(d):end); + a_corner = a(1:1+abs(d)-1); + corner_diag = m-abs(d); + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Super-diagonals + elseif d > 0 + a_bulk = a(1:end-d); + a_corner = a(end-d+1:end); + corner_diag = -m + d; + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Main diagonal + else + A = A + spdiagVariable(a, 0); + end + + end + +end \ No newline at end of file
diff -r 005a8d071da3 -r 0be9b4d6737b stripeMatrixPeriodic.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stripeMatrixPeriodic.m Tue May 22 13:31:09 2018 -0700 @@ -0,0 +1,8 @@ +% Creates a periodic discretization matrix of size n x n +% with the values of val on the diagonals diag. +% A = stripeMatrix(val,diags,n) +function A = stripeMatrixPeriodic(val,diags,n) + + D = ones(n,1)*val; + A = spdiagsVariablePeriodic(D,diags); +end \ No newline at end of file