Mercurial > repos > public > sbplib
changeset 707:0de70ec8bf60 feature/quantumTriangles
merge with feature/optim
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Fri, 10 Nov 2017 14:22:56 +0100 |
| parents | 7c16b5af8d98 (current diff) a95c89436916 (diff) |
| children | acb58769610e |
| files | +grid/Ti3D.m +multiblock/BoundaryGroupTest.m |
| diffstat | 72 files changed, 2186 insertions(+), 526 deletions(-) [+] |
line wrap: on
line diff
--- a/+blockmatrix/fromMatrixTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/fromMatrixTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -46,6 +46,20 @@ [23 5; 4 6; 10 12; 11 18], [7 14 16; 13 20 22; 19 21 3; 25 2 9]; }; }, + { + { + magic(3), + {[1 0 2],[1 2 0]} + }, + mat2cell(magic(3),[1 0 2],[1 2 0]) + }, + { + { + zeros(0,1), + {0,1}, + }, + {zeros(0,1)} + }, }; for i = 1:length(cases) out = convertToFull(blockmatrix.fromMatrix(cases{i}{1}{:}));
--- a/+blockmatrix/getDivision.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/getDivision.m Fri Nov 10 14:22:56 2017 +0100 @@ -16,7 +16,7 @@ m = zeros(1,size(C,2)); for j = 1:size(C,2) for i = 1:size(C,1) - if isempty(C{i,j}) + if isNullMatrix(C{i,j}) continue end m(j) = size(C{i,j},2); @@ -28,10 +28,15 @@ n = zeros(1,size(C,1)); for i = 1:size(C,1) for j = 1:size(C,2) - if isempty(C{i,j}) + if isNullMatrix(C{i,j}) continue end n(i) = size(C{i,j},1); end end -end \ No newline at end of file +end + +function b = isNullMatrix(A) + [n, m] = size(A); + b = n == 0 && m == 0; +end
--- a/+blockmatrix/getDivisionTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/getDivisionTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -56,6 +56,18 @@ }, {[2 1], 2} }, + { + {zeros(3,0)}, + {3, 0}, + }, + { + {zeros(3,0), zeros(3,0)}, + {3, [0, 0]}, + }, + { + {zeros(3,0); zeros(2,0)}, + {[3 2],0}, + }, }; for i = 1:length(cases)
--- a/+blockmatrix/isBlockmatrix.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/isBlockmatrix.m Fri Nov 10 14:22:56 2017 +0100 @@ -4,7 +4,7 @@ return end - % Make sure all blocks are numerica matrices + % Make sure all blocks are numerical matrices for i = 1:length(bm) if ~isnumeric(bm{i}) b = false;
--- a/+blockmatrix/isDivision.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/isDivision.m Fri Nov 10 14:22:56 2017 +0100 @@ -21,11 +21,11 @@ function b = isDivisionVector(v) if isempty(v) - b = false; + b = true; return end - if any(v <= 0) + if any(v < 0) b = false; return end
--- a/+blockmatrix/isDivisionTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/isDivisionTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -4,14 +4,16 @@ function testIsDivision(testCase) cases = { + {[1 2] ,false}, % Must be a cell array + {{[1 2 3]} ,false}, % Must have two vectors + {{[],[]}, true} % No blocks is a valid blockmatrix + {{[1 2],[]} ,true}, + {{[],[1 2]} ,true}, {{[2 2 2],[1 2]} ,true}, - {{[1 2],[1 0]} ,false}, - {{[0 2],[1 1]} ,false}, - {{[1 2],[]} ,false}, + {{[1 2],[1 0]} ,true}, + {{[0 2],[1 1]} ,true}, {{[1 2],[1]} ,true}, {{[1 2],[1], [1 2 3]} ,false}, - {{[1 2 3]} ,false}, - {[1 2] ,false}, }; for i = 1:length(cases)
--- a/+blockmatrix/toMatrix.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/toMatrix.m Fri Nov 10 14:22:56 2017 +0100 @@ -12,16 +12,12 @@ 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}) - continue + bm{i,j} = sparse(n(i),m(j)); end - % TODO: If this ever fails for large matrices. Try cell2mat instead. - A(n_ind(i)+1:n_ind(i+1),m_ind(j)+1:m_ind(j+1)) = bm{i,j}; end end + A = cell2mat(bm); end
--- a/+blockmatrix/toMatrixTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/toMatrixTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -51,6 +51,22 @@ [2 2 1; 2 1 2], }, + { + {zeros(0,0)}, + [], + }, + { + {zeros(3,0), zeros(3,0)}, + zeros(3,0), + }, + { + {zeros(3,0); zeros(2,0)}, + zeros(5,0), + }, + { + {zeros(0,3), zeros(0,2)}, + zeros(0,5), + }, }; for i = 1:length(cases)
--- a/+blockmatrix/zero.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/zero.m Fri Nov 10 14:22:56 2017 +0100 @@ -1,6 +1,6 @@ % Creates a block matrix according to the division with zeros everywhere. function bm = zero(div) - if ~blockmatrix.isDivision(div); + if ~blockmatrix.isDivision(div) error('div is not a valid division'); end
--- a/+blockmatrix/zeroTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+blockmatrix/zeroTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -32,6 +32,22 @@ {[1 2],[2 1]}, {[0 0],[0];[0 0; 0 0],[0; 0]}; }, + { + {[3],[0]}, + {zeros(3,0)}, + }, + + { + {[0],[3]}, + {zeros(0,3)}, + }, + { + {[0 2],[0 3]}, + { + zeros(0,0), zeros(0,3); + zeros(2,0), zeros(2,3); + }, + }, }; for i = 1:length(cases)
--- a/+grid/Cartesian.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+grid/Cartesian.m Fri Nov 10 14:22:56 2017 +0100 @@ -15,6 +15,8 @@ obj.d = length(varargin); for i = 1:obj.d + assert(isvector(varargin{i}), 'Coordinate inputs must be a vectors.') + obj.x{i} = varargin{i}; obj.m(i) = length(varargin{i}); end
--- a/+grid/Curvilinear.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+grid/Curvilinear.m Fri Nov 10 14:22:56 2017 +0100 @@ -19,6 +19,9 @@ % If mapping is a function evaluate it if isa(mapping, 'function_handle') + if nargin(mapping) ~= length(varargin) + error('The dimension of the mapping does not match the dimension of the logical coordinates') + end mapping = grid.evalOn(obj.logic, mapping); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Empty.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,53 @@ +classdef Empty < grid.Grid & grid.Structured + properties + dim + end + + methods + function obj = Empty(D) + obj.dim = D; + end + % n returns the number of points in the grid + function o = N(obj) + o = 0; + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.dim; + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = sparse(0,obj.dim); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + error('Restrict does not make sense for an empty grid') + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + error('Project does not make sense for an empty grid') + end + + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. + function bs = getBoundaryNames(obj) + bs = {}; + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = sparse(0,obj.dim-1); + end + + function h = scaling(obj) + h = 1; + end + + function s = size(obj) + s = zeros(1, obj.dim); + end + end +end \ No newline at end of file
--- a/+grid/Grid.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+grid/Grid.m Fri Nov 10 14:22:56 2017 +0100 @@ -16,7 +16,7 @@ % Projects the grid function gf on obj to the grid g. gf = projectFunc(obj, gf, g) - % Return the names of all boundaries in this grid. + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. bs = getBoundaryNames(obj) % Return coordinates for the given boundary
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/TODO.txt Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,1 @@ +% TODO: Rename grid package. name conflicts with built in function
--- a/+grid/Ti3D.m Thu Oct 05 18:04:23 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,253 +0,0 @@ -classdef Ti3D - properties - gs % {6}Surfaces - V % FunctionHandle(XI,ETA,ZETA) - end - - methods - % TODO write all fancy features for flipping around with the surfaces - % Each surface is defined with an outward facing outward and choosing - % the "corner" where XI=0 if not possible the corner where ETA=0 is choosen - function obj = Ti3D(CW,CE,CS,CN,CB,CT) - obj.gs = {CE,CW,CS,CN,CB,CT}; - - gw = CW.g; - ge = CE.g; - gs = CS.g; - gn = CN.g; - gb = CB.g; - gt = CT.g; - - function o = V_fun(XI,ETA,ZETA) - XI=XI'; - ETA=ETA'; - ZETA=ZETA'; - - one=0*ETA+1; - zero=0*ETA; - - Sw = gw(ETA,(1-ZETA)); - Se = ge((1-ETA),(1-ZETA)); - Ss = gs(XI,ZETA); - Sn = gn((1-XI),(1-ZETA)); - Sb = gb((1-XI),ETA); - St = gt(XI,ETA); - - Ewt = gw(ETA,zero); - Ewb = gw(ETA,one); - Ews = gw(zero,1-ZETA); - Ewn = gw(one,1-ZETA); - Eet = ge(1-ETA,zero); - Eeb = ge(1-ETA,one); - Ees = ge(one,1-ZETA); - Een = ge(zero,1-ZETA); - Enb = gn(1-XI,one); - Ent = gn(1-XI,zero); - Est = gs(XI,one); - Esb = gs(XI,zero); - - Cwbs = gw(zero,one); - Cwbn = gw(one,one); - Cwts = gw(zero,zero); - Cwtn = gw(one,zero); - Cebs = ge(one,one); - Cebn = ge(zero,one); - Cets = ge(one,zero); - Cetn = ge(zero,zero); - - - X1 = (1-XI).*Sw(1,:,:) + XI.*Se(1,:,:); - X2 = (1-ETA).*Ss(1,:,:) + ETA.*Sn(1,:,:); - X3 = (1-ZETA).*Sb(1,:,:) + ZETA.*St(1,:,:); - - X12 = (1-XI).*(1-ETA).*Ews(1,:,:) + (1-XI).*ETA.*Ewn(1,:,:) + XI.*(1-ETA).*Ees(1,:,:) + XI.*ETA.*Een(1,:,:); - X13 = (1-XI).*(1-ZETA).*Ewb(1,:,:) + (1-XI).*ZETA.*Ewt(1,:,:) + XI.*(1-ZETA).*Eeb(1,:,:) + XI.*ZETA.*Eet(1,:,:); - X23 = (1-ETA).*(1-ZETA).*Esb(1,:,:) + (1-ETA).*ZETA.*Est(1,:,:) + ETA.*(1-ZETA).*Enb(1,:,:) + ETA.*ZETA.*Ent(1,:,:); - - X123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(1,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(1,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(1,:,:) + ... - (1-XI).*ETA.*ZETA.*Cwtn(1,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(1,:,:) + XI.*(1-ETA).*ZETA.*Cets(1,:,:) + ... - XI.*ETA.*(1-ZETA).*Cebn(1,:,:) + XI.*ETA.*ZETA.*Cetn(1,:,:); - - X = X1 + X2 + X3 - X12 - X13 - X23 + X123; - - - Y1 = (1-XI).*Sw(2,:,:) + XI.*Se(2,:,:); - Y2 = (1-ETA).*Ss(2,:,:) + ETA.*Sn(2,:,:); - Y3 = (1-ZETA).*Sb(2,:,:) + ZETA.*St(2,:,:); - - Y12 = (1-XI).*(1-ETA).*Ews(2,:,:) + (1-XI).*ETA.*Ewn(2,:,:) + XI.*(1-ETA).*Ees(2,:,:) + XI.*ETA.*Een(2,:,:); - Y13 = (1-XI).*(1-ZETA).*Ewb(2,:,:) + (1-XI).*ZETA.*Ewt(2,:,:) + XI.*(1-ZETA).*Eeb(2,:,:) + XI.*ZETA.*Eet(2,:,:); - Y23 = (1-ETA).*(1-ZETA).*Esb(2,:,:) + (1-ETA).*ZETA.*Est(2,:,:) + ETA.*(1-ZETA).*Enb(2,:,:) + ETA.*ZETA.*Ent(2,:,:); - - Y123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(2,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(2,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(2,:,:) + ... - (1-XI).*ETA.*ZETA.*Cwtn(2,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(2,:,:) + XI.*(1-ETA).*ZETA.*Cets(2,:,:) + ... - XI.*ETA.*(1-ZETA).*Cebn(2,:,:) + XI.*ETA.*ZETA.*Cetn(2,:,:); - - Y = Y1 + Y2 + Y3 - Y12 - Y13 - Y23 + Y123; - - - Z1 = (1-XI).*Sw(3,:,:) + XI.*Se(3,:,:); - Z2 = (1-ETA).*Ss(3,:,:) + ETA.*Sn(3,:,:); - Z3 = (1-ZETA).*Sb(3,:,:) + ZETA.*St(3,:,:); - - Z12 = (1-XI).*(1-ETA).*Ews(3,:,:) + (1-XI).*ETA.*Ewn(3,:,:) + XI.*(1-ETA).*Ees(3,:,:) + XI.*ETA.*Een(3,:,:); - Z13 = (1-XI).*(1-ZETA).*Ewb(3,:,:) + (1-XI).*ZETA.*Ewt(3,:,:) + XI.*(1-ZETA).*Eeb(3,:,:) + XI.*ZETA.*Eet(3,:,:); - Z23 = (1-ETA).*(1-ZETA).*Esb(3,:,:) + (1-ETA).*ZETA.*Est(3,:,:) + ETA.*(1-ZETA).*Enb(3,:,:) + ETA.*ZETA.*Ent(3,:,:); - - Z123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(3,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(3,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(3,:,:) + ... - (1-XI).*ETA.*ZETA.*Cwtn(3,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(3,:,:) + XI.*(1-ETA).*ZETA.*Cets(3,:,:) + ... - XI.*ETA.*(1-ZETA).*Cebn(3,:,:) + XI.*ETA.*ZETA.*Cetn(3,:,:); - - Z = Z1 + Z2 + Z3 - Z12 - Z13 - Z23 + Z123; - o = [X;Y;Z]; - end - - obj.V = @V_fun; - end - - %Should be rewritten so that the input is xi eta zeta - function [X,Y,Z] = map(obj,XI,ETA,ZETA) - - V = obj.V; - - p = V(XI,ETA,ZETA); - X = p(1,:)'; - Y = p(2,:)'; - Z = p(3,:)'; - - end - - % function h = plot(obj,nu,nv) - % S = obj.S; - % - % default_arg('nv',nu) - % - % u = linspace(0,1,nu); - % v = linspace(0,1,nv); - % - % m = 100; - % - % X = zeros(nu+nv,m); - % Y = zeros(nu+nv,m); - % - % - % t = linspace(0,1,m); - % for i = 1:nu - % p = S(u(i),t); - % X(i,:) = p(1,:); - % Y(i,:) = p(2,:); - % end - % - % for i = 1:nv - % p = S(t,v(i)); - % X(i+nu,:) = p(1,:); - % Y(i+nu,:) = p(2,:); - % end - % - % h = line(X',Y'); - % end - % - % - % function h = show(obj,nu,nv) - % default_arg('nv',nu) - % S = obj.S; - % - % if(nu>2 || nv>2) - % h_grid = obj.plot(nu,nv); - % set(h_grid,'Color',[0 0.4470 0.7410]); - % end - % - % h_bord = obj.plot(2,2); - % set(h_bord,'Color',[0.8500 0.3250 0.0980]); - % set(h_bord,'LineWidth',2); - % end - % - % - % % TRANSFORMATIONS - % function ti = translate(obj,a) - % gs = obj.gs; - % - % for i = 1:length(gs) - % new_gs{i} = gs{i}.translate(a); - % end - % - % ti = grid.Ti(new_gs{:}); - % end - % - % % Mirrors the Ti so that the resulting Ti is still left handed. - % % (Corrected by reversing curves and switching e and w) - % function ti = mirror(obj, a, b) - % gs = obj.gs; - % - % new_gs = cell(1,4); - % - % new_gs{1} = gs{1}.mirror(a,b).reverse(); - % new_gs{3} = gs{3}.mirror(a,b).reverse(); - % new_gs{2} = gs{4}.mirror(a,b).reverse(); - % new_gs{4} = gs{2}.mirror(a,b).reverse(); - % - % ti = grid.Ti(new_gs{:}); - % end - % - % function ti = rotate(obj,a,rad) - % gs = obj.gs; - % - % for i = 1:length(gs) - % new_gs{i} = gs{i}.rotate(a,rad); - % end - % - % ti = grid.Ti(new_gs{:}); - % end - % - % function ti = rotate_edges(obj,n); - % new_gs = cell(1,4); - % for i = 0:3 - % new_i = mod(i - n,4); - % new_gs{new_i+1} = obj.gs{i+1}; - % end - % ti = grid.Ti(new_gs{:}); - % end - % end - % - % methods(Static) - % function obj = points(p1, p2, p3, p4) - % g1 = grid.Curve.line(p1,p2); - % g2 = grid.Curve.line(p2,p3); - % g3 = grid.Curve.line(p3,p4); - % g4 = grid.Curve.line(p4,p1); - % - % obj = grid.Ti(g1,g2,g3,g4); - % end - % - % function label(varargin) - % if nargin == 2 && ischar(varargin{2}) - % label_impl(varargin{:}); - % else - % for i = 1:length(varargin) - % label_impl(varargin{i},inputname(i)); - % end - % end - % - % - % function label_impl(ti,str) - % S = ti.S; - % - % pc = S(0.5,0.5); - % - % margin = 0.1; - % pw = S( margin, 0.5); - % pe = S(1-margin, 0.5); - % ps = S( 0.5, margin); - % pn = S( 0.5, 1-margin); - % - % - % ti.show(2,2); - % grid.place_label(pc,str); - % grid.place_label(pw,'w'); - % grid.place_label(pe,'e'); - % grid.place_label(ps,'s'); - % grid.place_label(pn,'n'); - % end - % end - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryIdentifier.txt Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,5 @@ +A grid.boundaryIdentifier identifies a boundary of a grid. +For a Cartesian grid it is simply 's', 'n', 'w', 'e'. +For a multiblock grid it might be something like {1, 's'}. +For some other grid it will be up to the developer to chose +a good way to identify boundaries.
--- a/+grid/evalOn.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+grid/evalOn.m Fri Nov 10 14:22:56 2017 +0100 @@ -7,39 +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 and convert to cell array for easier use as a parameter - x = num2cell(g.points()); - - % Find the number of components - x0 = x(1,:); +% 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); - - if size(f0,2) ~= 1 - error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') - end +end - gf = zeros(g.N*k, 1); - % keyboard - for i = 1:g.N - % (1 + (i-1)*k):(i*k) - % func(x{i,:}) - gf((1 + (i-1)*k):(i*k)) = func(x{i,:}); +% 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 -end \ No newline at end of file +end
--- a/+grid/evalOnTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+grid/evalOnTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -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');
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/gridFunction.txt Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,1 @@ +TODO: Add documentation for gridfunctions here \ No newline at end of file
--- a/+multiblock/BoundaryGroup.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+multiblock/BoundaryGroup.m Fri Nov 10 14:22:56 2017 +0100 @@ -1,30 +1,10 @@ % BoundaryGroup defines a boundary grouping in a multiblock grid. -classdef BoundaryGroup - properties - blockIDs - names - end - +% It workds like a cell array and collects boundary identifiers +% Within the multiblock package a BoundaryGroup is a valid boundary identifier as well. +classdef BoundaryGroup < Cell methods - function obj = BoundaryGroup(varargin) - % Input arguemnts are arbitrary number or 1x2 cell arrays - % representing each boundary in the group. - % The 1st element of the cell array is an integer defining which grid it belongs to. - % The 2nd element of the cell array is the name of the boundary within the block. - % - % Ex: - % bg = multiblock.BoundaryGroup({1,'n'},{1,'s'},{2,'s'}) - - - obj.blockIDs = []; - obj.names = {}; - for i = 1:length(varargin) - if ~iscell(varargin{i}) || ~all(size(varargin{i}) == [1 2]) - error('multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', 'Inputs must be 1x2 cell arrays'); - end - obj.blockIDs(i) = varargin{i}{1}; - obj.names{i} = varargin{i}{2}; - end + function obj = BoundaryGroup(data) + obj = obj@Cell(data); end function display(obj, name) @@ -33,19 +13,7 @@ disp([name, ' =']) disp(' ') - if length(obj.names) == 1 - fprintf(' {}\n\n') - return - end - - fprintf(' {') - - fprintf('%d:%s', obj.blockIDs(1), obj.names{1}) - for i = 2:length(obj.names) - fprintf(', %d:%s', obj.blockIDs(i), obj.names{i}); - end - - fprintf('}\n\n') + fprintf(' BoundaryGroup%s\n\n', toString(obj.data)); end end -end \ No newline at end of file +end
--- a/+multiblock/BoundaryGroupTest.m Thu Oct 05 18:04:23 2017 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -function tests = BoundaryGroupTest() - tests = functiontests(localfunctions); -end - -function testCreation(testCase) - in = {{3,'n'},{2,'hoho'},{1,'s'}}; - - blockIDs = [3 2 1]; - names = {'n', 'hoho', 's'}; - - bg = multiblock.BoundaryGroup(in{:}); - testCase.verifyEqual(bg.blockIDs, blockIDs); - testCase.verifyEqual(bg.names, names); -end - -function testInputError(testCase) - in = { - {'n', 's'}, - {{3,'n'},{2,2,'hoho'},{1,'s'}}, - }; - - out = { - 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', - 'multiblock:BoundaryGroup:BoundaryGroup:InvalidInput', - }; - - for i = 1:length(in) - testCase.verifyError(@()multiblock.BoundaryGroup(in{i}{:}), out{i}); - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Contour.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,65 @@ +classdef Contour < handle + properties + grid + contours + nContours + + ZData + CData + + end + + methods + function obj = Contour(g, gf, nContours) + obj.grid = g; + obj.nContours = nContours; + + coords = obj.grid.points(); + X = obj.grid.funcToPlotMatrices(coords(:,1)); + Y = obj.grid.funcToPlotMatrices(coords(:,2)); + + V = obj.grid.funcToPlotMatrices(gf); + + + holdState = ishold(); + hold on + + contours = {1, obj.grid.nBlocks}; + for i = 1:obj.grid.nBlocks + [~, contours{i}] = contour(X{i}, Y{i}, V{i},obj.nContours); + contours{i}.LevelList = contours{1}.LevelList; + end + + if holdState == false + hold off + end + + obj.contours = [contours{:}]; + + obj.ZData = gf; + obj.CData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.contours, propertyName, propertyValue); + end + + function obj = set.ZData(obj, gf) + obj.ZData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.contours(i).ZData = V{i}; + end + end + + function obj = set.CData(obj, gf) + obj.CData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.contours(i).CData = V{i}; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Def.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,101 @@ +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 + +
--- a/+multiblock/DiffOp.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+multiblock/DiffOp.m Fri Nov 10 14:22:56 2017 +0100 @@ -5,12 +5,12 @@ diffOps D H - + blockmatrixDiv end - + methods - function obj = DiffOp(doHand, grid, order, doParam,timeDep) + function obj = DiffOp(doHand, grid, order, doParam) % doHand -- may either be a function handle or a cell array of % function handles for each grid. The function handle(s) % should be on the form do = doHand(grid, order, ...) @@ -25,13 +25,13 @@ % doHand(..., doParam{i}{:}) Otherwise doParam is sent as % extra parameters to all doHand: doHand(..., doParam{:}) default_arg('doParam', []) - + [getHand, getParam] = parseInput(doHand, grid, doParam); - + nBlocks = grid.nBlocks(); - + obj.order = order; - + % Create the diffOps for each block obj.diffOps = cell(1, nBlocks); for i = 1:nBlocks @@ -42,73 +42,48 @@ end obj.diffOps{i} = h(grid.grids{i}, order, p{:}); end - - + + % Build the norm matrix H = cell(nBlocks, nBlocks); for i = 1:nBlocks H{i,i} = obj.diffOps{i}.H; end obj.H = blockmatrix.toMatrix(H); - - + + % Build the differentiation matrix - switch timeDep - case {'n','no','N','No'} - obj.blockmatrixDiv = {grid.Ns, grid.Ns}; - D = blockmatrix.zero(obj.blockmatrixDiv); - for i = 1:nBlocks - D{i,i} = obj.diffOps{i}.D; + obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + D = blockmatrix.zero(obj.blockmatrixDiv); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; + end + + for i = 1:nBlocks + for j = 1:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue end - - for i = 1:nBlocks - for j = 1:nBlocks - intf = grid.connections{i,j}; - if isempty(intf) - continue - end - - - [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); - 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}); - D{j,j} = D{j,j} + jj; - D{j,i} = D{j,i} + ji; - end - end - obj.D = blockmatrix.toMatrix(D); - case {'y','yes','Y','Yes'} - for i = 1:nBlocks - D{i,i} = @(t)obj.diffOps{i}.D(t); - end - - for i = 1:nBlocks - for j = 1:nBlocks - intf = grid.connections{i,j}; - if isempty(intf) - continue - end - - [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); - D{i,i} = @(t)D{i,i}(t) + ii(t); - D{i,j} = ij; - - [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); - D{j,j} = @(t)D{j,j}(t) + jj(t); - D{j,i} = ji; - end - end - obj.D = D; + + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + 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}); + D{j,j} = D{j,j} + jj; + D{j,i} = D{j,i} + ji; + end end - - + obj.D = blockmatrix.toMatrix(D); + + function [getHand, getParam] = parseInput(doHand, grid, doParam) if ~isa(grid, 'multiblock.Grid') error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); end - + if iscell(doHand) && length(doHand) == grid.nBlocks() getHand = @(i)doHand{i}; elseif isa(doHand, 'function_handle') @@ -116,62 +91,89 @@ else error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); end - + if isempty(doParam) getParam = @(i){}; return end - + if ~iscell(doParam) getParam = @(i)doParam; return end - + % doParam is a non-empty cell-array - + if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) % doParam is a cell-array of cell-arrays getParam = @(i)doParam{i}; return end - + getParam = @(i)doParam; end end - + function ops = splitOp(obj, op) % Splits a matrix operator into a cell-matrix of matrix operators for % each grid. ops = sparse2cell(op, obj.NNN); end - - function op = getBoundaryOperator(obj, op, boundary) - if iscell(boundary) - localOpName = [op '_' boundary{2}]; - blockId = boundary{1}; - localOp = obj.diffOps{blockId}.(localOpName); - - div = {obj.blockmatrixDiv{1}, size(localOp,2)}; - blockOp = blockmatrix.zero(div); - blockOp{blockId,1} = localOp; - op = blockmatrix.toMatrix(blockOp); - return - else - % Boundary är en sträng med en boundary group i. + + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperator(obj, opName, boundary) + switch class(boundary) + case 'cell' + localOpName = [opName '_' boundary{2}]; + blockId = boundary{1}; + localOp = obj.diffOps{blockId}.(localOpName); + + div = {obj.blockmatrixDiv{1}, size(localOp,2)}; + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + return + case 'multiblock.BoundaryGroup' + op = sparse(size(obj.D,1),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperator(opName, boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') end end - + % Creates the closure and penalty matrix for a given boundary condition, % boundary -- the name of the boundary on the form {id,name} where % id is the number of a block and name is the name of a - % boundary of that block example: {1,'s'} or {3,'w'} + % boundary of that block example: {1,'s'} or {3,'w'}. It + % can also be a boundary group function [closure, penalty] = boundary_condition(obj, boundary, type) + switch class(boundary) + case 'cell' + [closure, penalty] = obj.singleBoundaryCondition(boundary, type); + case 'multiblock.BoundaryGroup' + [n,m] = size(obj.D); + closure = sparse(n,m); + penalty = sparse(n,0); + for i = 1:length(boundary) + [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type); + closure = closure + closurePart; + penalty = [penalty, penaltyPart]; + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function [closure, penalty] = singleBoundaryCondition(obj, boundary, type) I = boundary{1}; name = boundary{2}; - + % Get the closure and penaly matrices [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); - + % Expand to matrix for full domain. div = obj.blockmatrixDiv; if ~iscell(blockClosure) @@ -185,25 +187,27 @@ closure{i} = blockmatrix.toMatrix(temp); end end - - div{2} = size(blockPenalty, 2); % Penalty is a column vector + 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 end - + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - + error('not implemented') end - + % Size returns the number of degrees of freedom function N = size(obj) N = 0;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTimeDep.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,210 @@ +classdef DiffOpTimeDep < scheme.Scheme + properties + grid + order + diffOps + D + H + + blockmatrixDiv + end + + methods + function obj = DiffOpTimeDep(doHand, grid, order, doParam) + % doHand -- may either be a function handle or a cell array of + % function handles for each grid. The function handle(s) + % should be on the form do = doHand(grid, order, ...) + % Additional parameters for each doHand may be provided in + % the doParam input. + % grid -- a multiblock grid + % order -- integer specifying the order of accuracy + % doParam -- may either be a cell array or a cell array of cell arrays + % for each block. If it is a cell array with length equal + % to the number of blocks then each element is sent to the + % corresponding function handle as extra parameters: + % doHand(..., doParam{i}{:}) Otherwise doParam is sent as + % extra parameters to all doHand: doHand(..., doParam{:}) + default_arg('doParam', []) + + [getHand, getParam] = parseInput(doHand, grid, doParam); + obj.grid = grid; + nBlocks = grid.nBlocks(); + + obj.order = order; + + % Create the diffOps for each block + obj.diffOps = cell(1, nBlocks); + for i = 1:nBlocks + h = getHand(i); + p = getParam(i); + if ~iscell(p) + p = {p}; + end + obj.diffOps{i} = h(grid.grids{i}, order, p{:}); + end + + + % Build the norm matrix + H = cell(nBlocks, nBlocks); + for i = 1:nBlocks + H{i,i} = obj.diffOps{i}.H; + end + obj.H = blockmatrix.toMatrix(H); + + + % Build the differentiation matrix + obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + % D = blockmatrix.zero(obj.blockmatrixDiv); + + for i = 1:nBlocks + for j = 1:nBlocks + D{i,j} = @(t) 0; + D{j,i} = @(t) 0; + end + end + + + for i = 1:nBlocks + D{i,i} = @(t)obj.diffOps{i}.D(t); + end + + for i = 1:nBlocks + for j = 1:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = @(t) D{i,i}(t) + ii(t); + D{i,j} = @(t) D{i,j}(t) + ij(t); + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = @(t) D{j,j}(t) + jj(t); + D{j,i} = @(t) D{j,i}(t) + ji(t); + end + end + obj.D = D; + + + function [getHand, getParam] = parseInput(doHand, grid, doParam) + if ~isa(grid, 'multiblock.Grid') + error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); + end + + if iscell(doHand) && length(doHand) == grid.nBlocks() + getHand = @(i)doHand{i}; + elseif isa(doHand, 'function_handle') + getHand = @(i)doHand; + else + error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); + end + + if isempty(doParam) + getParam = @(i){}; + return + end + + if ~iscell(doParam) + getParam = @(i)doParam; + return + end + + % doParam is a non-empty cell-array + + if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) + % doParam is a cell-array of cell-arrays + getParam = @(i)doParam{i}; + return + end + + getParam = @(i)doParam; + end + end + + function ops = splitOp(obj, op) + % Splits a matrix operator into a cell-matrix of matrix operators for + % each grid. + ops = sparse2cell(op, obj.NNN); + end + + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperator(obj, opName, boundary) + switch class(boundary) + case 'cell' + localOpName = [opName '_' boundary{2}]; + blockId = boundary{1}; + localOp = obj.diffOps{blockId}.(localOpName); + + div = {obj.blockmatrixDiv{1}, size(localOp,2)}; + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + return + case 'multiblock.BoundaryGroup' + op = sparse(size(obj.D,1),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperator(opName, boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + end + end + + % Creates the closure and penalty matrix for a given boundary condition, + % boundary -- the name of the boundary on the form {id,name} where + % id is the number of a block and name is the name of a + % boundary of that block example: {1,'s'} or {3,'w'}. It + % can also be a boundary group + function [closure, penalty] = boundary_condition(obj, boundary, type) + switch class(boundary) + case 'cell' + [closure, penalty] = obj.singleBoundaryCondition(boundary, type); + case 'multiblock.BoundaryGroup' + nBlocks = obj.grid.nBlocks(); + %[n,m] = size(obj.D); + %closure = sparse(n,m); + %penalty = sparse(n,0); + % closure =@(t)0; + % penalty = @(t)0; + for i = 1:nBlocks + for j = 1:nBlocks + closure{j,i} = @(t)0; + penalty{j,i} = @(t)0; + end + end + + + for i = 1:length(boundary) + [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type); + closure{i,i} = @(t)closure{i,i}(t) + closurePart(t); + penalty{i,i} = @(t)penalty{i,i}(t) + penaltyPart(t); + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function [blockClosure, blockPenalty] = singleBoundaryCondition(obj, boundary, type) + I = boundary{1}; + name = boundary{2}; + % Get the closure and penaly matrices + [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('not implemented') + end + + % Size returns the number of degrees of freedom + function N = size(obj) + N = 0; + for i = 1:length(obj.diffOps) + N = N + obj.diffOps{i}.size(); + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/EmptyGrid.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,47 @@ +classdef EmptyGrid < grid.Empty & multiblock.Grid + methods + function obj = EmptyGrid(D) + obj@multiblock.Grid({},cell(0,0)); + obj@grid.Empty(D); + end + + % n returns the number of points in the grid + function o = N(obj) + o = N@grid.Empty(obj); + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = D@grid.Empty(obj); + end + + % points returns a n x d matrix containing the coordinates for all points. + function X = points(obj) + X = points@grid.Empty(obj); + end + + % Restricts the grid function gf on obj to the subgrid g. + function gf = restrictFunc(obj, gf, g) + gf = restrictFunc@grid.Empty(obj); + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + gf = projectFunc@grid.Empty(obj); + end + + % Return the grid.boundaryIdentifiers of all boundaries in a cell array. + function bs = getBoundaryNames(obj) + bs = getBoundaryNames@grid.Empty(obj); + end + + % Return coordinates for the given boundary + function b = getBoundary(obj, name) + b = getBoundary@grid.Empty(name); + end + + function s = size(obj) + s = size@multiblock.Grid(obj); + end + end +end
--- a/+multiblock/Grid.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+multiblock/Grid.m Fri Nov 10 14:22:56 2017 +0100 @@ -9,16 +9,20 @@ % General multiblock grid methods - - % grids -- cell array of N grids - % connections -- NxN upper triangular cell matrix. connections{i,j} - % specifies the connection between block i and j. If - % it's empty there is no connection otherwise it's a 2 - % -cell-vector with strings naming the boundaries to be - % connected. (inverted coupling?) - %% Should we have boundary groups at all? maybe it can be handled in a - %% cleaner way outside of the class. + % grids -- cell array of N grids + % connections -- NxN upper triangular cell matrix. connections{i,j} + % specifies the connection between block i and j. If + % it's empty there is no connection otherwise it's a 2 + % -cell-vector with strings naming the boundaries to be + % connected. (inverted coupling?) + % boundaryGroups -- A struct of BoundaryGroups. The field names of the + % struct are the names of each boundary group. + % The boundary groups can be used to collect block + % boundaries into physical boundaries to simplify + % getting boundary operators and setting boundary conditions function obj = Grid(grids, connections, boundaryGroups) + default_arg('boundaryGroups', struct()); + assertType(grids, 'cell') obj.grids = grids; obj.connections = connections; @@ -27,7 +31,7 @@ obj.nPoints = obj.nPoints + grids{i}.N(); end - % if iscell(boundaryGroups) + obj.boundaryGroups = boundaryGroups; end function n = size(obj) @@ -42,7 +46,7 @@ % Ns returns the number of points in each sub grid as a vector function o = Ns(obj) ns = zeros(1,obj.nBlocks); - for i = 1:obj.nBlocks; + for i = 1:obj.nBlocks ns(i) = obj.grids{i}.N(); end o = ns; @@ -59,7 +63,7 @@ % points returns a n x d matrix containing the coordinates for all points. function X = points(obj) - X = []; + X = sparse(0,obj.D()); for i = 1:length(obj.grids) X = [X; obj.grids{i}.points]; end @@ -76,9 +80,29 @@ N(i) = obj.grids{i}.N(); end - gfs = mat2cell(gf, N, 1); + gfs = blockmatrix.fromMatrix(gf, {N,1}); end + % TODO: Split op? + % Should the method to split an operator be moved here instead of being in multiblock.DiffOp? + + % Converts a gridfunction to a set of plot matrices + % Takes a grid function and and a structured grid. + function F = funcToPlotMatrices(obj, gf) + % TODO: This method should problably not be here. + % The funcToPlotMatrix uses .size poperty of the grids + % Which doesn't always exist for all types of grids. + % It's only valid for structured grids + gfs = obj.splitFunc(gf); + + F = cell(1, obj.nBlocks()); + + for i = 1:obj.nBlocks() + F{i} = grid.funcToPlotMatrix(obj.grids{i}, gfs{i}); + end + end + + % Restricts the grid function gf on obj to the subgrid g. function gf = restrictFunc(obj, gf, g) gfs = obj.splitFunc(gf); @@ -108,13 +132,43 @@ end + % Find all non interface boundaries of all blocks. + % Return their grid.boundaryIdentifiers in a cell array. function bs = getBoundaryNames(obj) - bs = []; + bs = {}; + for i = 1:obj.nBlocks() + candidates = obj.grids{i}.getBoundaryNames(); + for j = 1:obj.nBlocks() + if ~isempty(obj.connections{i,j}) + candidates = setdiff(candidates, obj.connections{i,j}{1}); + end + + if ~isempty(obj.connections{j,i}) + candidates = setdiff(candidates, obj.connections{j,i}{2}); + end + end + + for k = 1:length(candidates) + bs{end+1} = {i, candidates{k}}; + end + end end - % Return coordinates for the given boundary - function b = getBoundary(obj, name) - b = []; + % Return coordinates for the given boundary/boundaryGroup + function b = getBoundary(obj, boundary) + switch class(boundary) + case 'cell' + I = boundary{1}; + name = boundary{2}; + b = obj.grids{I}.getBoundary(name); + case 'multiblock.BoundaryGroup' + b = sparse(0,obj.D()); + for i = 1:length(boundary) + b = [b; obj.getBoundary(boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + end end end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Line.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,47 @@ +classdef Line < handle + properties + grid + lines + + YData + end + + methods + function obj = Line(g, gf) + assertType(g, 'multiblock.Grid') + obj.grid = g; + + X = obj.grid.splitFunc(obj.grid.points()); + Y = obj.grid.splitFunc(gf); + + holdState = ishold(); + hold on + + lines = cell(1, obj.grid.nBlocks); + for i = 1:obj.grid.nBlocks + lines{i} = plot(X{i}, Y{i}); + end + + if holdState == false + hold off + end + + obj.lines = [lines{:}]; + + obj.YData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.lines, propertyName, propertyValue); + end + + function obj = set.YData(obj, gf) + obj.YData = gf; + + Y = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.lines(i).YData = Y{i}; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Surface.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,62 @@ +classdef Surface < handle + properties + grid + surfs + + ZData + CData + + end + + methods + function obj = Surface(g, gf) + obj.grid = g; + + coords = obj.grid.points(); + X = obj.grid.funcToPlotMatrices(coords(:,1)); + Y = obj.grid.funcToPlotMatrices(coords(:,2)); + + V = obj.grid.funcToPlotMatrices(gf); + + + holdState = ishold(); + hold on + + surfs = cell(1, obj.grid.nBlocks); + for i = 1:obj.grid.nBlocks + surfs{i} = surf(X{i}, Y{i}, V{i}); + end + + if holdState == false + hold off + end + + obj.surfs = [surfs{:}]; + + obj.ZData = gf; + obj.CData = gf; + end + + function set(obj, propertyName, propertyValue) + set(obj.surfs, propertyName, propertyValue); + end + + function obj = set.ZData(obj, gf) + obj.ZData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.surfs(i).ZData = V{i}; + end + end + + function obj = set.CData(obj, gf) + obj.CData = gf; + + V = obj.grid.funcToPlotMatrices(gf); + for i = 1:obj.grid.nBlocks + obj.surfs(i).CData = V{i}; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/Animation.m Fri Nov 10 14:22:56 2017 +0100 @@ -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
--- a/+noname/animate.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+noname/animate.m Fri Nov 10 14:22:56 2017 +0100 @@ -1,10 +1,10 @@ -% hand = noname.animate(discretization, time_modifier, Tend, dirname, opt) +% noname.animate(discretization, time_modifier, Tend, dirname, opt) % % Example: % noname.animate(discr,timemodifier,tend) % noname.animate(discr,1, [tstart tend],'my_mov', opt) -function hand = animate(discretization, time_modifier, Tend, dirname, opt) +function animate(discretization, time_modifier, Tend, dirname, opt) default_arg('time_modifier', 1); default_arg('Tend', Inf); default_arg('dirname', ''); @@ -87,7 +87,8 @@ pause anim.animate(@G, Tstart, Tend, time_modifier); else - while true + pause + while ts.t < Tend ts.step(); sol = discretization.getTimeSnapshot(ts); update(sol);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/calculateErrors.m Fri Nov 10 14:22:56 2017 +0100 @@ -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));
--- a/+parametrization/Ti.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+parametrization/Ti.m Fri Nov 10 14:22:56 2017 +0100 @@ -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; @@ -182,6 +195,15 @@ obj = parametrization.Ti(g1,g2,g3,g4); end + function obj = rectangle(a, b) + p1 = a; + p2 = [b(1), a(2)]; + p3 = b; + p4 = [a(1), b(2)]; + + obj = parametrization.Ti.points(p1,p2,p3,p4); + end + % Like the constructor but allows inputing line curves as 2-cell arrays: % example: parametrization.Ti.linesAndCurves(g1, g2, {a, b} g4) function obj = linesAndCurves(C1, C2, C3, C4)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/Ti3D.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,253 @@ +classdef Ti3D + properties + gs % {6}Surfaces + V % FunctionHandle(XI,ETA,ZETA) + end + + methods + % TODO write all fancy features for flipping around with the surfaces + % Each surface is defined with an outward facing outward and choosing + % the "corner" where XI=0 if not possible the corner where ETA=0 is choosen + function obj = Ti3D(CW,CE,CS,CN,CB,CT) + obj.gs = {CE,CW,CS,CN,CB,CT}; + + gw = CW.g; + ge = CE.g; + gs = CS.g; + gn = CN.g; + gb = CB.g; + gt = CT.g; + + function o = V_fun(XI,ETA,ZETA) + XI=XI'; + ETA=ETA'; + ZETA=ZETA'; + + one=0*ETA+1; + zero=0*ETA; + + Sw = gw(ETA,(1-ZETA)); + Se = ge((1-ETA),(1-ZETA)); + Ss = gs(XI,ZETA); + Sn = gn((1-XI),(1-ZETA)); + Sb = gb((1-XI),ETA); + St = gt(XI,ETA); + + Ewt = gw(ETA,zero); + Ewb = gw(ETA,one); + Ews = gw(zero,1-ZETA); + Ewn = gw(one,1-ZETA); + Eet = ge(1-ETA,zero); + Eeb = ge(1-ETA,one); + Ees = ge(one,1-ZETA); + Een = ge(zero,1-ZETA); + Enb = gn(1-XI,one); + Ent = gn(1-XI,zero); + Est = gs(XI,one); + Esb = gs(XI,zero); + + Cwbs = gw(zero,one); + Cwbn = gw(one,one); + Cwts = gw(zero,zero); + Cwtn = gw(one,zero); + Cebs = ge(one,one); + Cebn = ge(zero,one); + Cets = ge(one,zero); + Cetn = ge(zero,zero); + + + X1 = (1-XI).*Sw(1,:,:) + XI.*Se(1,:,:); + X2 = (1-ETA).*Ss(1,:,:) + ETA.*Sn(1,:,:); + X3 = (1-ZETA).*Sb(1,:,:) + ZETA.*St(1,:,:); + + X12 = (1-XI).*(1-ETA).*Ews(1,:,:) + (1-XI).*ETA.*Ewn(1,:,:) + XI.*(1-ETA).*Ees(1,:,:) + XI.*ETA.*Een(1,:,:); + X13 = (1-XI).*(1-ZETA).*Ewb(1,:,:) + (1-XI).*ZETA.*Ewt(1,:,:) + XI.*(1-ZETA).*Eeb(1,:,:) + XI.*ZETA.*Eet(1,:,:); + X23 = (1-ETA).*(1-ZETA).*Esb(1,:,:) + (1-ETA).*ZETA.*Est(1,:,:) + ETA.*(1-ZETA).*Enb(1,:,:) + ETA.*ZETA.*Ent(1,:,:); + + X123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(1,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(1,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(1,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(1,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(1,:,:) + XI.*(1-ETA).*ZETA.*Cets(1,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(1,:,:) + XI.*ETA.*ZETA.*Cetn(1,:,:); + + X = X1 + X2 + X3 - X12 - X13 - X23 + X123; + + + Y1 = (1-XI).*Sw(2,:,:) + XI.*Se(2,:,:); + Y2 = (1-ETA).*Ss(2,:,:) + ETA.*Sn(2,:,:); + Y3 = (1-ZETA).*Sb(2,:,:) + ZETA.*St(2,:,:); + + Y12 = (1-XI).*(1-ETA).*Ews(2,:,:) + (1-XI).*ETA.*Ewn(2,:,:) + XI.*(1-ETA).*Ees(2,:,:) + XI.*ETA.*Een(2,:,:); + Y13 = (1-XI).*(1-ZETA).*Ewb(2,:,:) + (1-XI).*ZETA.*Ewt(2,:,:) + XI.*(1-ZETA).*Eeb(2,:,:) + XI.*ZETA.*Eet(2,:,:); + Y23 = (1-ETA).*(1-ZETA).*Esb(2,:,:) + (1-ETA).*ZETA.*Est(2,:,:) + ETA.*(1-ZETA).*Enb(2,:,:) + ETA.*ZETA.*Ent(2,:,:); + + Y123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(2,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(2,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(2,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(2,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(2,:,:) + XI.*(1-ETA).*ZETA.*Cets(2,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(2,:,:) + XI.*ETA.*ZETA.*Cetn(2,:,:); + + Y = Y1 + Y2 + Y3 - Y12 - Y13 - Y23 + Y123; + + + Z1 = (1-XI).*Sw(3,:,:) + XI.*Se(3,:,:); + Z2 = (1-ETA).*Ss(3,:,:) + ETA.*Sn(3,:,:); + Z3 = (1-ZETA).*Sb(3,:,:) + ZETA.*St(3,:,:); + + Z12 = (1-XI).*(1-ETA).*Ews(3,:,:) + (1-XI).*ETA.*Ewn(3,:,:) + XI.*(1-ETA).*Ees(3,:,:) + XI.*ETA.*Een(3,:,:); + Z13 = (1-XI).*(1-ZETA).*Ewb(3,:,:) + (1-XI).*ZETA.*Ewt(3,:,:) + XI.*(1-ZETA).*Eeb(3,:,:) + XI.*ZETA.*Eet(3,:,:); + Z23 = (1-ETA).*(1-ZETA).*Esb(3,:,:) + (1-ETA).*ZETA.*Est(3,:,:) + ETA.*(1-ZETA).*Enb(3,:,:) + ETA.*ZETA.*Ent(3,:,:); + + Z123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(3,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(3,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(3,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(3,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(3,:,:) + XI.*(1-ETA).*ZETA.*Cets(3,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(3,:,:) + XI.*ETA.*ZETA.*Cetn(3,:,:); + + Z = Z1 + Z2 + Z3 - Z12 - Z13 - Z23 + Z123; + o = [X;Y;Z]; + end + + obj.V = @V_fun; + end + + %Should be rewritten so that the input is xi eta zeta + function [X,Y,Z] = map(obj,XI,ETA,ZETA) + + V = obj.V; + + p = V(XI,ETA,ZETA); + X = p(1,:)'; + Y = p(2,:)'; + Z = p(3,:)'; + + end + + % function h = plot(obj,nu,nv) + % S = obj.S; + % + % default_arg('nv',nu) + % + % u = linspace(0,1,nu); + % v = linspace(0,1,nv); + % + % m = 100; + % + % X = zeros(nu+nv,m); + % Y = zeros(nu+nv,m); + % + % + % t = linspace(0,1,m); + % for i = 1:nu + % p = S(u(i),t); + % X(i,:) = p(1,:); + % Y(i,:) = p(2,:); + % end + % + % for i = 1:nv + % p = S(t,v(i)); + % X(i+nu,:) = p(1,:); + % Y(i+nu,:) = p(2,:); + % end + % + % h = line(X',Y'); + % end + % + % + % function h = show(obj,nu,nv) + % default_arg('nv',nu) + % S = obj.S; + % + % if(nu>2 || nv>2) + % h_grid = obj.plot(nu,nv); + % set(h_grid,'Color',[0 0.4470 0.7410]); + % end + % + % h_bord = obj.plot(2,2); + % set(h_bord,'Color',[0.8500 0.3250 0.0980]); + % set(h_bord,'LineWidth',2); + % end + % + % + % % TRANSFORMATIONS + % function ti = translate(obj,a) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.translate(a); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % % Mirrors the Ti so that the resulting Ti is still left handed. + % % (Corrected by reversing curves and switching e and w) + % function ti = mirror(obj, a, b) + % gs = obj.gs; + % + % new_gs = cell(1,4); + % + % new_gs{1} = gs{1}.mirror(a,b).reverse(); + % new_gs{3} = gs{3}.mirror(a,b).reverse(); + % new_gs{2} = gs{4}.mirror(a,b).reverse(); + % new_gs{4} = gs{2}.mirror(a,b).reverse(); + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate(obj,a,rad) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.rotate(a,rad); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate_edges(obj,n); + % new_gs = cell(1,4); + % for i = 0:3 + % new_i = mod(i - n,4); + % new_gs{new_i+1} = obj.gs{i+1}; + % end + % ti = grid.Ti(new_gs{:}); + % end + % end + % + % methods(Static) + % function obj = points(p1, p2, p3, p4) + % g1 = grid.Curve.line(p1,p2); + % g2 = grid.Curve.line(p2,p3); + % g3 = grid.Curve.line(p3,p4); + % g4 = grid.Curve.line(p4,p1); + % + % obj = grid.Ti(g1,g2,g3,g4); + % end + % + % function label(varargin) + % if nargin == 2 && ischar(varargin{2}) + % label_impl(varargin{:}); + % else + % for i = 1:length(varargin) + % label_impl(varargin{i},inputname(i)); + % end + % end + % + % + % function label_impl(ti,str) + % S = ti.S; + % + % pc = S(0.5,0.5); + % + % margin = 0.1; + % pw = S( margin, 0.5); + % pe = S(1-margin, 0.5); + % ps = S( 0.5, margin); + % pn = S( 0.5, 1-margin); + % + % + % ti.show(2,2); + % grid.place_label(pc,str); + % grid.place_label(pw,'w'); + % grid.place_label(pe,'e'); + % grid.place_label(ps,'s'); + % grid.place_label(pn,'n'); + % end + % end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/TiTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -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
--- a/+sbp/D1Gauss.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+sbp/D1Gauss.m Fri Nov 10 14:22:56 2017 +0100 @@ -37,5 +37,9 @@ obj.borrowing = []; end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end end end
--- a/+sbp/D2BlockNorm.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+sbp/D2BlockNorm.m Fri Nov 10 14:22:56 2017 +0100 @@ -50,6 +50,10 @@ obj.m = m; end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end end
--- a/+sbp/D4Compatible.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+sbp/D4Compatible.m Fri Nov 10 14:22:56 2017 +0100 @@ -28,8 +28,8 @@ methods function obj = D4Compatible(m,lim,order) - - + + x_l = lim{1}; x_r = lim{2}; L = x_r-x_l; @@ -68,6 +68,10 @@ end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end end
--- a/+sbp/D4Standard.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+sbp/D4Standard.m Fri Nov 10 14:22:56 2017 +0100 @@ -28,7 +28,7 @@ methods function obj = D4Standard(m,lim,order) - + x_l = lim{1}; x_r = lim{2}; L = x_r-x_l; @@ -56,6 +56,10 @@ obj.m = m; end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end end
--- a/+scheme/Beam.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+scheme/Beam.m Fri Nov 10 14:22:56 2017 +0100 @@ -90,6 +90,9 @@ gamm = obj.gamm; delt = obj.delt; + + % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ... + switch type case {'dn', 'clamped'} % Dirichlet-neumann boundary condition alpha = obj.alpha; @@ -123,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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/LaplaceCurvilinear.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,349 @@ +classdef LaplaceCurvilinear < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + + order % Order accuracy for the approximation + + a,b % Parameters of the operator + + + % Inner products and operators for physical coordinates + D % Laplace operator + H, Hi % Inner product + e_w, e_e, e_s, e_n + d_w, d_e, d_s, d_n % Normal derivatives at the boundary + H_w, H_e, H_s, H_n % Boundary inner products + Dx, Dy % Physical derivatives + M % Gradient inner product + + % Metric coefficients + J, Ji + a11, a12, a22 + x_u + x_v + y_u + y_v + + % Inner product and operators for logical coordinates + H_u, H_v % Norms in the x and y directions + Hi_u, Hi_v + Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hiu, Hiv + du_w, dv_w + du_e, dv_e + du_s, dv_s + du_n, dv_n + gamm_u, gamm_v + lambda + 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) + default_arg('opSet',@sbp.D2Variable); + default_arg('a', 1); + default_arg('b', 1); + + if b ~=1 + error('Not implemented yet') + end + + assert(isa(g, 'grid.Curvilinear')) + + m = g.size(); + m_u = m(1); + m_v = m(2); + m_tot = g.N(); + + h = g.scaling(); + h_u = h(1); + h_v = h(2); + + + % 1D operators + ops_u = opSet(m_u, {0, 1}, order); + ops_v = opSet(m_v, {0, 1}, order); + + I_u = speye(m_u); + I_v = speye(m_v); + + D1_u = ops_u.D1; + D2_u = ops_u.D2; + H_u = ops_u.H; + Hi_u = ops_u.HI; + e_l_u = ops_u.e_l; + e_r_u = ops_u.e_r; + d1_l_u = ops_u.d1_l; + d1_r_u = ops_u.d1_r; + + D1_v = ops_v.D1; + D2_v = ops_v.D2; + H_v = ops_v.H; + Hi_v = ops_v.HI; + e_l_v = ops_v.e_l; + e_r_v = ops_v.e_r; + d1_l_v = ops_v.d1_l; + d1_r_v = ops_v.d1_r; + + + % Logical operators + Du = kr(D1_u,I_v); + Dv = kr(I_u,D1_v); + obj.Hu = kr(H_u,I_v); + obj.Hv = kr(I_u,H_v); + obj.Hiu = kr(Hi_u,I_v); + obj.Hiv = kr(I_u,Hi_v); + + e_w = kr(e_l_u,I_v); + e_e = kr(e_r_u,I_v); + e_s = kr(I_u,e_l_v); + e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (e_w'*Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (e_e'*Dv)'; + obj.du_s = (e_s'*Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (e_n'*Du)'; + obj.dv_n = kr(I_u,d1_r_v); + + + % Metric coefficients + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + x_u = Du*x; + x_v = Dv*x; + y_u = Du*y; + y_v = Dv*y; + + J = x_u.*y_v - x_v.*y_u; + a11 = 1./J .* (x_v.^2 + y_v.^2); + a12 = -1./J .* (x_u.*x_v + y_u.*y_v); + a22 = 1./J .* (x_u.^2 + y_u.^2); + lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); + + obj.x_u = x_u; + obj.x_v = x_v; + obj.y_u = y_u; + obj.y_v = y_v; + + + % Assemble full operators + L_12 = spdiag(a12); + Duv = Du*L_12*Dv; + Dvu = Dv*L_12*Du; + + Duu = sparse(m_tot); + Dvv = sparse(m_tot); + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m_v + D = D2_u(a11(ind(:,i))); + p = ind(:,i); + Duu(p,p) = D; + end + + for i = 1:m_u + D = D2_v(a22(ind(i,:))); + p = ind(i,:); + Dvv(p,p) = D; + end + + + % Physical operators + obj.J = spdiag(J); + obj.Ji = spdiag(1./J); + + obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); + obj.H = obj.J*kr(H_u,H_v); + obj.Hi = obj.Ji*kr(Hi_u,Hi_v); + + obj.e_w = e_w; + obj.e_e = e_e; + obj.e_s = e_s; + obj.e_n = e_n; + + %% normal derivatives + I_w = ind(1,:); + I_e = ind(end,:); + I_s = ind(:,1); + I_n = ind(:,end); + + a11_w = spdiag(a11(I_w)); + a12_w = spdiag(a12(I_w)); + a11_e = spdiag(a11(I_e)); + a12_e = spdiag(a12(I_e)); + a22_s = spdiag(a22(I_s)); + a12_s = spdiag(a12(I_s)); + a22_n = spdiag(a22(I_n)); + a12_n = spdiag(a12(I_n)); + + s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); + s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); + s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); + s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); + + obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))'; + obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))'; + obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))'; + obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))'; + + obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; + obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; + + %% Boundary inner products + obj.H_w = H_v*spdiag(s_w); + obj.H_e = H_v*spdiag(s_e); + obj.H_s = H_u*spdiag(s_s); + obj.H_n = H_u*spdiag(s_n); + + % Misc. + obj.m = m; + obj.h = [h_u h_v]; + obj.order = order; + obj.grid = g; + + obj.a = a; + obj.b = b; + obj.a11 = a11; + obj.a12 = a12; + obj.a22 = a22; + obj.lambda = lambda; + + obj.gamm_u = h_u*ops_u.borrowing.M.d1; + obj.gamm_v = h_v*ops_v.borrowing.M.d1; + 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, parameter) + default_arg('type','neumann'); + default_arg('parameter', []); + + [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); + switch type + % Dirichlet boundary condition + case {'D','d','dirichlet'} + tuning = 1.2; + % tuning = 20.2; + + b1 = gamm*obj.lambda./obj.a11.^2; + b2 = gamm*obj.lambda./obj.a22.^2; + + tau1 = tuning * spdiag(-1./b1 - 1./b2); + tau2 = 1; + + tau = (tau1*e + tau2*d)*H_b; + + closure = obj.a*obj.Hi*tau*e'; + penalty = -obj.a*obj.Hi*tau; + + + % Neumann boundary condition + case {'N','n','neumann'} + tau1 = -1; + tau2 = 0; + tau = (tau1*e + tau2*d)*H_b; + + 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) + % 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; + + 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; + + tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); + tau1 = tuning * spdiag(tau1); + tau2 = 1/2; + + sig1 = -1/2; + sig2 = 0; + + tau = (e_u*tau1 + tau2*d_u)*H_b_u; + sig = (sig1*e_u + sig2*d_u)*H_b_u; + + 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. + % The right boundary is considered the positive boundary + % + % I -- the indecies 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 + case 'w' + e = obj.e_w; + d = obj.d_w; + H_b = obj.H_w; + I = ind(1,:); + case 'e' + e = obj.e_e; + d = obj.d_e; + H_b = obj.H_e; + I = ind(end,:); + case 's' + e = obj.e_s; + d = obj.d_s; + H_b = obj.H_s; + I = ind(:,1)'; + case 'n' + e = obj.e_n; + d = obj.d_n; + H_b = obj.H_n; + I = ind(:,end)'; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','e'} + gamm = obj.gamm_u; + case {'s','n'} + gamm = obj.gamm_v; + end + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- a/+scheme/Scheme.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+scheme/Scheme.m Fri Nov 10 14:22:56 2017 +0100 @@ -25,9 +25,12 @@ % neighbour_boundary is a string specifying which boundary to % interface to. % penalty may be a cell array if there are several penalties with different weights - [closure, penalty] = boundary_condition(obj,boundary,type) + [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % TODO: op = getBoundaryOperator()?? + % makes sense to have it available through a method instead of random properties + % Returns the number of degrees of freedom. N = size(obj) end
--- a/+scheme/Schrodinger1dCurve.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+scheme/Schrodinger1dCurve.m Fri Nov 10 14:22:56 2017 +0100 @@ -122,14 +122,7 @@ tau1 = @(t) s * 1i*obj.Ji(p,p)*d; tau2 = @(t) (1*s*obj.a(p,p))/2*e; closure = @(t)obj.Hi*sqrt(obj.Ji)*(tau1(t)*e' + tau2(obj.a)*e')*sqrt(obj.Ji); - - switch class(data) - case 'double' - penalty = @(t) -obj.Hi*sqrt(obj.Ji)*(tau1*data+tau2(obj.a)*data)*sqrt(obj.Ji); - - otherwise - error('Weird data argument!') - end + penalty = @(t) -obj.Hi*sqrt(obj.Ji)*(tau1(t)*e' + tau2(obj.a)*e')*sqrt(obj.Ji); % Unknown, boundary condition otherwise error('No such boundary condition: type = %s',type);
--- a/+scheme/Schrodinger2dCurve.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+scheme/Schrodinger2dCurve.m Fri Nov 10 14:22:56 2017 +0100 @@ -38,6 +38,7 @@ m_tot, m_u, m_v p Ji, J + Hamiltonian end methods @@ -67,7 +68,7 @@ obj.D1_u = ops_u.D1; obj.D2_u = ops_u.D2; - + H_u = ops_u.H; Hi_u = ops_u.HI; e_l_u = ops_u.e_l; @@ -115,10 +116,11 @@ obj.h = [h_u h_v]; obj.order = order; obj.D = @(t)obj.d_fun(t); + obj.Hamiltonian = @(t)obj.h_fun(t); obj.variable_update(0); end - function [D] = d_fun(obj,t) + function [D,n] = d_fun(obj,t) % obj.update_vairables(t); In driscretization? % D = obj.Ji*(-1/2*(obj.b1*obj.Du-obj.b1_u+obj.Du*obj.b1) - % 1/2*(obj.b2*obj.Dv - obj.b2_v +obj.Dv*obj.b2) + @@ -128,6 +130,9 @@ D = sqrt(obj.Ji)*(-1/2*(obj.b1*obj.Du + obj.Du*obj.b1) - 1/2*(obj.b2*obj.Dv + obj.Dv*obj.b2) + 1i*obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV))*sqrt(obj.Ji); end + function [Hamiltonian] = h_fun(obj,t) + Hamiltonian = -sqrt(obj.Ji)*(obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV))*sqrt(obj.Ji); + end function [D ]= variable_update(obj,t) % Metric derivatives @@ -204,7 +209,7 @@ % 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,~) + function [closure, penalty,closureHamiltonian,penaltyHamiltonian] = boundary_condition(obj, boundary,~) [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); a_t = @(t) spdiag(coeff_t(t)); @@ -221,6 +226,9 @@ closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); + + closureHamiltonian = @(t) 1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji); + penaltyHamiltonian = @(t) -1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/TODO.txt Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,1 @@ +% TODO: Rename package and abstract class to diffOp
--- a/+scheme/Wave2dCurve.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+scheme/Wave2dCurve.m Fri Nov 10 14:22:56 2017 +0100 @@ -27,6 +27,8 @@ gamm_u, gamm_v lambda + Dx, Dy % Physical derivatives + x_u x_v y_u @@ -38,6 +40,8 @@ default_arg('opSet',@sbp.D2Variable); default_arg('c', 1); + warning('Use LaplaceCruveilinear instead') + assert(isa(g, 'grid.Curvilinear')) m = g.size(); @@ -153,6 +157,9 @@ obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); obj.lambda = lambda; + obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; + obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; + obj.gamm_u = h_u*ops_u.borrowing.M.d1; obj.gamm_v = h_v*ops_v.borrowing.M.d1; end @@ -225,6 +232,7 @@ tau = -c.^2 * 1/beta*obj.Ji*e; + warning('is this right?! /c?') closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e'; closure{2} = halfnorm_inv*tau*beta*d'; penalty = -halfnorm_inv*tau;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/bcSetup.m Fri Nov 10 14:22:56 2017 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Arnoldi.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,37 @@ +function y = Arnoldi(A,v,t,tol,m) + +n = size (v,1); +V = zeros(n ,m+1); +H = zeros(m+1,m); + +beta = norm(v); +V(:,1) = v/beta; +resnorm = inf; + +j=0; + +while resnorm > tol + j = j+1; + z = A*V(:,j); + for i=1:j + H(i,j) = z'*V(:,i); + z = z - H(i,j)*V(:,i); + end + H(j+1,j) = norm(z); + e1 = zeros(j,1); e1(1) = 1; + ej = zeros(j,1); ej(j) = 1; + s = [0.01, 1/3, 2/3, 1]*t; + for q=1:length(s) + u = expm(-s(q)*H(1:j,1:j))*e1; + beta_j(q) = -H(j+1,j)* (ej'*u); + end + resnorm = norm(beta_j,'inf'); + fprintf('j = %d, resnorm = %.2e\n',j,resnorm); + if resnorm<=toler + break + elseif j==m + disp('warning: no convergence within m steps'); + end + V(:,j+1) = z/H(j+1,j); +end +y = V(:,1:j)*(beta*u);
--- a/+time/+expint/Magnus_mp.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/+expint/Magnus_mp.m Fri Nov 10 14:22:56 2017 +0100 @@ -13,7 +13,7 @@ case 'expm' v = expm(k*Omega)*v; case 'Arnoldi' - v = time.expint.expm_Arnoldi(-Omega,v,k,tol,100); + v = time.expint.expm_Arnoldi(-Omega,v,k,tol,1000); otherwise error('No such matrix exponential evaluation')
--- a/+time/+expint/Magnus_mp.m~ Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/+expint/Magnus_mp.m~ Fri Nov 10 14:22:56 2017 +0100 @@ -4,14 +4,18 @@ function v = Magnus_mp(v,D, t , k,matrixexp,tol) if isa(D,'function_handle') - switch matrixexp - case 'expm' - v = expm(k*D(t +k/2))*v; - case 'Arnol' - v = time.expint.expm_Arnoldi(-D(t +k/2),v,k,toler,100); + Omega = D(t +k/2); else - %v = krylov(k*D,v); - % v = expm(k*D)*v; + Omega = D; end +switch matrixexp + case 'expm' + v = expm(k*Omega)*v; + case 'Arnoldi' + [v, i] = time.expint.expm_Arnoldi(-Omega,v,k,tol,1000); + otherwise + error('No such matrix exponential evaluation') + +end end \ No newline at end of file
--- a/+time/+expint/expm_Arnoldi.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/+expint/expm_Arnoldi.m Fri Nov 10 14:22:56 2017 +0100 @@ -1,5 +1,7 @@ function y = expm_Arnoldi(A,v,t,toler,m) -% + +global iter + % y = expm_Arnoldi(A,v,t,toler,m) % % computes $y = \exp(-t A) v$ @@ -26,8 +28,11 @@ resnorm = inf; j=0; +iter = 0; + while resnorm > toler + iter = iter +1; j = j+1; w = A*V(:,j); for i=1:j
--- a/+time/CdiffImplicit.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/CdiffImplicit.m Fri Nov 10 14:22:56 2017 +0100 @@ -63,7 +63,9 @@ v_prev = f1; I = speye(m); - v = (1/k^2*A)\((1/k^2*A - 1/2*B)*f1 + (1/k*I - 1/2*C)*f2 + 1/2*G(0)); + % v = (1/k^2*A)\((1/k^2*A - 1/2*B)*f1 + (1/k*I - 1/2*C)*f2 + 1/2*G(0)); + v = f1 + k*f2; + if ~issparse(A) || ~issparse(B) || ~issparse(C) error('LU factorization with full pivoting only works for sparse matrices.')
--- a/+time/Rungekutta4SecondOrder.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/Rungekutta4SecondOrder.m Fri Nov 10 14:22:56 2017 +0100 @@ -15,7 +15,19 @@ methods - % Solves u_tt = Du + Eu_t + S + % Solves u_tt = Du + Eu_t + S by + % Rewriting on first order form: + % w_t = M*w + C(t) + % where + % M = [ + % 0, I; + % D, E; + % ] + % and + % C(t) = [ + % 0; + % S(t) + % ] % D, E, S can either all be constants or all be function handles, % They can also be omitted by setting them equal to the empty matrix. function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) @@ -41,23 +53,15 @@ S = @(t)S; end - I = speye(obj.m); - O = sparse(obj.m,obj.m); - - obj.M = @(t)[ - O, I; - D(t), E(t); - ]; - obj.C = @(t)[ - zeros(obj.m,1); - S(t); - ]; - obj.k = k; obj.t = t0; obj.w = [v0; v0t]; - obj.F = @(w,t)(obj.M(t)*w + obj.C(t)); + % Avoid matrix formulation because it is VERY slow + obj.F = @(w,t)[ + w(obj.m+1:end); + D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); + ]; else default_arg('D', sparse(obj.m, obj.m));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeImplicitFormulation.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,130 @@ +classdef SBPInTimeImplicitFormulation < time.Timestepper + % The SBP in time method. + % Implemented for A*v_t = B*v + f(t), v(0) = v0 + properties + A,B + f + + k % total time step. + + blockSize % number of points in each block + N % Number of components + + order + nodes + + M,K % System matrices + L,U,p,q % LU factorization of M + e_T + + % Time state + t + v + n + end + + methods + function obj = SBPInTimeImplicitFormulation(A, B, f, k, t0, v0, TYPE, order, blockSize) + + default_arg('TYPE','gauss'); + default_arg('f',[]); + + if(strcmp(TYPE,'gauss')) + default_arg('order',4) + default_arg('blockSize',4) + else + default_arg('order', 8); + default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE)); + end + + obj.A = A; + obj.B = B; + + if ~isempty(f) + obj.f = f; + else + obj.f = @(t)sparse(length(v0),1); + end + + obj.k = k; + obj.blockSize = blockSize; + obj.N = length(v0); + + obj.n = 0; + obj.t = t0; + + %==== Build the time discretization matrix =====% + switch TYPE + case 'equidistant' + ops = sbp.D2Standard(blockSize,{0,obj.k},order); + case 'optimal' + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order); + case 'minimal' + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal'); + case 'gauss' + ops = sbp.D1Gauss(blockSize,{0,obj.k}); + end + + I = speye(size(A)); + I_t = speye(blockSize,blockSize); + + D1 = kron(ops.D1, I); + HI = kron(ops.HI, I); + e_0 = kron(ops.e_l, I); + e_T = kron(ops.e_r, I); + obj.nodes = ops.x; + + % Convert to form M*w = K*v0 + f(t) + tau = kron(I_t, A) * e_0; + M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B); + + K = HI*tau; + + obj.M = M; + obj.K = K; + obj.e_T = e_T; + + % LU factorization + [obj.L,obj.U,obj.p,obj.q] = lu(obj.M, 'vector'); + + obj.v = v0; + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + RHS = zeros(obj.blockSize*obj.N,1); + + for i = 1:obj.blockSize + RHS((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i)); + end + + RHS = RHS + obj.K*obj.v; + + y = obj.L\RHS(obj.p); + z = obj.U\y; + + w = zeros(size(z)); + w(obj.q) = z; + + obj.v = obj.e_T'*w; + + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + methods(Static) + function N = smallestBlockSize(order,TYPE) + default_arg('TYPE','gauss') + + switch TYPE + case 'gauss' + N = 4; + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeSecondOrderFormImplicit.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,80 @@ +classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper + properties + A, B, C, f + AA, BB, ff + + n + t + k + + firstOrderTimeStepper + end + + methods + % Solves A*u_tt + B*u_t + C*u = f(t) + % A, B can either both be constants or both be function handles, + % They can also be omitted by setting them equal to the empty matrix. + function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize) + default_arg('f', []); + default_arg('TYPE', []); + default_arg('order', []); + default_arg('blockSize',[]); + + m = length(v0); + + default_arg('A', sparse(m, m)); + default_arg('B', sparse(m, m)); + default_arg('C', sparse(m, m)); + + I = speye(m); + O = sparse(m,m); + + % Rewrite to + % AA*w_t = BB*w + ff(t); + + obj.AA = [ + I, O; + O, A; + ]; + obj.BB = [ + O, I; + -C, -B; + ]; + + if ~isempty(f) + obj.ff = @(t)[ + sparse(m,1); + f(t); + ]; + else + obj.ff = @(t) sparse(2*m,1); + end + + w0 = [v0; v0t]; + + obj.k = k; + obj.t = t0; + obj.n = 0; + + obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize); + end + + function [v,t] = getV(obj) + w = obj.firstOrderTimeStepper.getV(); + v = w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + w = obj.firstOrderTimeStepper.getV(); + vt = w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.firstOrderTimeStepper.step(); + obj.t = obj.firstOrderTimeStepper.t; + obj.n = obj.firstOrderTimeStepper.n; + end + end +end
--- a/+time/Timestepper.m Thu Oct 05 18:04:23 2017 +0200 +++ b/+time/Timestepper.m Fri Nov 10 14:22:56 2017 +0100 @@ -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);
--- a/Cell.m Thu Oct 05 18:04:23 2017 +0200 +++ b/Cell.m Fri Nov 10 14:22:56 2017 +0100 @@ -23,6 +23,10 @@ s = size(A.data); end + function b = isempty(A) + b = prod(size(A)) == 0; + end + function l = length(A) l = length(A.data); end
--- a/CellTest.m Thu Oct 05 18:04:23 2017 +0200 +++ b/CellTest.m Fri Nov 10 14:22:56 2017 +0100 @@ -36,6 +36,21 @@ end end +function testIsEmpty(testCase) + cases = { + {cell(0,0), true}, + {cell(1,0), true}, + {cell(0,1), true}, + {cell(1,1), false}, + }; + + for i = 1:length(cases) + A = Cell(cases{i}{1}); + expected = cases{i}{2}; + testCase.verifyEqual(isempty(A),expected); + end +end + function testTranspose(testCase) testCase.verifyEqual(Cell({1i, 2}).', Cell({1i; 2})); testCase.verifyEqual(Cell({1i; 2}).', Cell({1i, 2}));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertNumberOfArguments.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,5 @@ +function assertNumberOfArguments(fun, N) + if nargin(fun) ~= N + error('sbplib:assertNumberOfArguments:wrongNumberOfArguments', '"%s" must have %d, found %d', inputname(1), N, nargin(fun)); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertScalar.m Fri Nov 10 14:22:56 2017 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertSymbolic.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,3 @@ +function assertSymbolic(s) + assert(logical(simplify(s))); +end
--- a/assertType.m Thu Oct 05 18:04:23 2017 +0200 +++ b/assertType.m Fri Nov 10 14:22:56 2017 +0100 @@ -1,5 +1,11 @@ function assertType(obj, type) - if ~isa(obj, type) - error('sbplib:assertType:wrongType', '"%s" must have type "%s", found "%s"', inputname(1), type, class(obj)); + if ~iscell(type) + if ~isa(obj, type) + error('sbplib:assertType:wrongType', '"%s" must have type "%s", found "%s"', inputname(1), type, class(obj)); + end + else + if ~isAnyOf(obj, type) + error('sbplib:assertType:wrongType', '"%s" must be one of the types %s, found "%s"', inputname(1), toString(type), class(obj)); + end end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/four.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,19 @@ +% four returns the fourier transform u_hat of the function u and the frequencies w +function [w, u_hat] = four(x, u) + u_hat = fft(u); + + N = length(x); + L = x(end) - x(1); + + k = shift_k(0:N-1); + + u_hat = fftshift(u_hat); + + dw = 2*pi/L; + w = dw*k; +end + +function k_shifted = shift_k(k) + N = length(k); + k_shifted = [-floor(N/2):-1, 0, 1:ceil(N/2)-1]; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fourInv.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,4 @@ +function u = ifour(u_hat) + u_hat = ifftshift(u_hat); + u = ifft(u_hat); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gaussian.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,3 @@ +function z = gaussian(x,x0,d) + z = exp(-norm(x-x0).^2/d^2); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/isAnyOf.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,11 @@ +% Returns true of obj is any of he types in the cell array types +% b = isAnyOf(obj, types) +function b = isAnyOf(obj, types) + for i = 1:length(types) + if isa(obj, types{i}); + b = true; + return + end + end + b = false; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sbplibLocation.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,4 @@ +function location = sbplibLocation() + scriptname = mfilename('fullpath'); + [location, ~, ~] = fileparts(scriptname); +end
--- a/sbplibVersion.m Thu Oct 05 18:04:23 2017 +0200 +++ b/sbplibVersion.m Fri Nov 10 14:22:56 2017 +0100 @@ -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'; + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/skewPart.m Fri Nov 10 14:22:56 2017 +0100 @@ -0,0 +1,4 @@ +% Returns the skew of A +function S = skewPart(A, tol) + S = 1/2*(A - A'); +end