Mercurial > repos > public > sbplib
changeset 219:f66513508c75 feature/beams
Merged with feature/grids.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 27 Jun 2016 13:26:02 +0200 |
| parents | fc07ebc49412 (diff) da058ce66876 (current diff) |
| children | 5df8d20281fe |
| files | |
| diffstat | 10 files changed, 204 insertions(+), 59 deletions(-) [+] |
line wrap: on
line diff
--- a/+blockmatrix/isBlockmatrixTest.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+blockmatrix/isBlockmatrixTest.m Mon Jun 27 13:26:02 2016 +0200 @@ -52,6 +52,10 @@ }, true % Empty blocks allowed. }, + { + blockmatrix.zero({[1 2 3],[2 3]}), + true % A zero block matrix is a block matrix + }, }; for i = 1:length(cases)
--- a/+blockmatrix/zero.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+blockmatrix/zero.m Mon Jun 27 13:26:02 2016 +0200 @@ -1,5 +1,9 @@ % Creates a block matrix according to the division with zeros everywhere. function bm = zero(div) + if ~blockmatrix.isDivision(div); + error('div is not a valid division'); + end + n = div{1}; m = div{2};
--- a/+multiblock/DiffOp.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+multiblock/DiffOp.m Mon Jun 27 13:26:02 2016 +0200 @@ -64,11 +64,11 @@ [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; + 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; + D{j,i} = ji; end end obj.D = blockmatrix.toMatrix(D); @@ -108,12 +108,12 @@ % 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: 1s or 3w - function [closure, penalty] = boundary_condition(obj,boundary,type,data) - I = boundary(1) - name = boundary(2:end); + function [closure, penalty] = boundary_condition(obj,boundary,type) + I = boundary{1}; + name = boundary{2}; % Get the closure and penaly matrices - [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type, data); + [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); % Expand to matrix for full domain. div = obj.blockmatrixDiv; @@ -126,7 +126,7 @@ % Convert to matrix closure = blockmatrix.toMatrix(closure); - penalty = blockmatrix.toMatrix(closure); + penalty = blockmatrix.toMatrix(penalty); end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
--- a/+multiblock/Grid.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+multiblock/Grid.m Mon Jun 27 13:26:02 2016 +0200 @@ -62,9 +62,9 @@ nBlocks = length(obj.grids); % Collect number of points in each block - N = cell(1,nBlocks); + N = zeros(1,nBlocks); for i = 1:nBlocks - N{i} = obj.grids{i}.N(); + N(i) = obj.grids{i}.N(); end gfs = mat2cell(gf, N, 1);
--- a/+multiblock/stitchSchemes.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+multiblock/stitchSchemes.m Mon Jun 27 13:26:02 2016 +0200 @@ -11,24 +11,24 @@ % Output parameters are cell arrays and cell matrices. % % Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) -function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) +function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound) default_arg('schmParam',[]); - n_blocks = numel(blocks); + n_blocks = numel(grids); % Creating Schemes for i = 1:n_blocks if isempty(schmParam); - schms{i} = schmHand(ms{i},blocks{i},order,[]); + schms{i} = schmHand(grids{i},order,[]); elseif ~iscell(schmParam) param = schmParam(i); - schms{i} = schmHand(ms{i},blocks{i},order,param); + schms{i} = schmHand(grids{i},order,param); else param = schmParam{i}; if iscell(param) - schms{i} = schmHand(ms{i},blocks{i},order,param{:}); + schms{i} = schmHand(grids{i},order,param{:}); else - schms{i} = schmHand(ms{i},blocks{i},order,param); + schms{i} = schmHand(grids{i},order,param); end end
--- a/+noname/convergence.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+noname/convergence.m Mon Jun 27 13:26:02 2016 +0200 @@ -1,5 +1,5 @@ % Reference is either a key or a function handle -function [q, e, h, runtime] = convergence(filename, errorFunc, reference, method, order, m, T) +function [q, e, h, runtime] = convergence(filename, errorFunc, reference, name, order, m, T) default_arg('errorFunc', @scheme.error1d); sf = SolutionFile(filename); @@ -7,7 +7,7 @@ % Generate convergence, error, and efficiency plots for each search key with more than one entry. for i = 1:length(m) - key.method = method; + key.name = name; key.order = order; key.m = m(i); key.T = T; @@ -30,18 +30,11 @@ % Get the reference solution vector if isa(reference,'function_handle'); - x = v_repr.x; + x = v_repr.grid.Points(); v_ref = reference(x,T); else % Downsample the reference solution - x = v_repr.x; - x_ref = reference.x; - - [~,I] = ismember(x,x_ref,'rows'); - if any(I == 0) - error('Solution and reference solution seem to be on different grids.'); - end - v_ref = reference.v(I); + v_ref = reference.grid.restrictFunc(reference.v, v_repr.grid); end e = errorFunc(discr,v, v_ref);
--- a/+noname/plotSolutions.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+noname/plotSolutions.m Mon Jun 27 13:26:02 2016 +0200 @@ -10,7 +10,7 @@ key = sf.keys{i}; entry = sf.get(key); - method = key.method; + name = key.name; order = key.order; m = key.m; T = key.T; @@ -23,11 +23,11 @@ update(repr); if save_figures - figname = sprintf('%s_%s_o%d_m%d_T%d',figename_prefix,method,order,m,i); + figname = sprintf('%s_%s_o%d_m%d_T%d',figename_prefix,name,order,m,i); fprintf('Saving figure to ''%s''\n',figname); saveeps(hand,figname); end end -end \ No newline at end of file +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Beam.m Mon Jun 27 13:26:02 2016 +0200 @@ -0,0 +1,160 @@ +classdef Beam < scheme.Scheme + properties + order % Order accuracy for the approximation + grid + + D % non-stabalized scheme operator + alpha + + H % Discrete norm + Hi + + e_l, e_r + d1_l, d1_r + d2_l, d2_r + d3_l, d3_r + gamm + delt + end + + methods + function obj = Beam(grid, order, alpha, opsGen) + default_arg('alpha', 1); + default_arg('opsGen', @sbp.Higher); + + if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 1 + error('Grid must be 1d cartesian'); + end + + obj.grid = grid; + obj.order = order; + obj.alpha = alpha; + + m = grid.m; + h = grid.scaling(); + + ops = opsGen(m, h, order); + + I = speye(m); + + D4 = sparse(ops.derivatives.D4); + obj.H = sparse(ops.norms.H); + obj.Hi = sparse(ops.norms.HI); + obj.e_l = sparse(ops.boundary.e_1); + obj.e_r = sparse(ops.boundary.e_m); + obj.d1_l = sparse(ops.boundary.S_1); + obj.d1_r = sparse(ops.boundary.S_m); + obj.d2_l = sparse(ops.boundary.S2_1); + obj.d2_r = sparse(ops.boundary.S2_m); + obj.d3_l = sparse(ops.boundary.S3_1); + obj.d3_r = sparse(ops.boundary.S3_m); + + obj.D = alpha*D4; + + obj.gamm = h*ops.borrowing.N.S2/2; + obj.delt = h^3*ops.borrowing.N.S3/2; + 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. + % 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_e, penalty_d] = boundary_condition(obj,boundary,type) + default_arg('type','dn'); + + [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary); + gamm = obj.gamm; + delt = obj.delt; + + switch type + case {'dn'} % Dirichlet-neumann boundary condition + alpha = obj.alpha; + + % tau1 < -alpha^2/gamma + tuning = 1.1; + + tau1 = tuning * alpha/delt; + tau4 = s*alpha; + + sig2 = tuning * alpha/gamm; + sig3 = -s*alpha; + + tau = tau1*e+tau4*d3; + sig = sig2*d1+sig3*d2; + + closure = obj.Hi*(tau*e' + sig*d1'); + + penalty_e = obj.Hi*tau; + penalty_d = obj.Hi*sig; + otherwise % Unknown, boundary condition + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary); + [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + gamm_u = obj.gamm; + delt_u = obj.delt; + + gamm_v = neighbour_scheme.gamm; + delt_v = neighbour_scheme.delt; + + tuning = 2; + + alpha_u = obj.alpha; + alpha_v = neighbour_scheme.alpha; + + tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; + % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; + tau4 = s_u*alpha_u/2; + + sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning; + sig3 = -s_u*alpha_u/2; + + phi2 = s_u*1/2; + + psi1 = -s_u*1/2; + + tau = tau1*e_u + tau4*d3_u; + sig = sig2*d1_u + sig3*d2_u ; + phi = phi2*d1_u ; + psi = psi1*e_u ; + + closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); + penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); + end + + % Returns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d1 = obj.d1_l; + d2 = obj.d2_l; + d3 = obj.d3_l; + s = -1; + case 'r' + e = obj.e_r; + d1 = obj.d1_r; + d2 = obj.d2_r; + d3 = obj.d3_r; + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = prod(obj.m); + end + + end +end
--- a/+scheme/Beam2d.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+scheme/Beam2d.m Mon Jun 27 13:26:02 2016 +0200 @@ -1,10 +1,6 @@ classdef Beam2d < scheme.Scheme properties - m % Number of points in each direction, possibly a vector - N % Number of points total - h % Grid spacing - u,v % Grid - x,y % Values of x and y for each grid point + grid order % Order accuracy for the approximation D % non-stabalized scheme operator @@ -27,21 +23,23 @@ methods function obj = Beam2d(m,lim,order,alpha,opsGen) + default_arg('alpha',1); default_arg('opsGen',@sbp.Higher); - default_arg('a',1); - if length(m) == 1 - m = [m m]; + if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 2 + error('Grid must be 2d cartesian'); end - m_x = m(1); - m_y = m(2); + obj.grid = grid; + obj.alpha = alpha; + obj.order = order; - xlim = lim{1}; - ylim = lim{2}; + m_x = grid.m(1); + m_y = grid.m(2); - [x, h_x] = util.get_grid(xlim{:},m_x); - [y, h_y] = util.get_grid(ylim{:},m_y); + h = grid.scaling(); + h_x = h(1); + h_y = h(2); ops_x = opsGen(m_x,h_x,order); ops_y = opsGen(m_y,h_y,order); @@ -49,9 +47,6 @@ I_x = speye(m_x); I_y = speye(m_y); - - - D4_x = sparse(ops_x.derivatives.D4); H_x = sparse(ops_x.norms.H); Hi_x = sparse(ops_x.norms.HI); @@ -105,16 +100,7 @@ obj.d3_s = kr(I_x,d3_l_y); obj.d3_n = kr(I_x,d3_r_y); - obj.m = m; - obj.h = [h_x h_y]; - obj.order = order; - - obj.alpha = alpha; obj.D = alpha*D4; - obj.u = x; - obj.v = y; - obj.x = kr(x,ones(m_y,1)); - obj.y = kr(ones(m_x,1),y); obj.gamm_x = h_x*ops_x.borrowing.N.S2/2; obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2;
--- a/+scheme/Scheme.m Mon Jun 27 13:24:59 2016 +0200 +++ b/+scheme/Scheme.m Mon Jun 27 13:26:02 2016 +0200 @@ -20,13 +20,11 @@ % '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. - [closure, penalty] = boundary_condition(obj,boundary,type,data) + [closure, penalty] = boundary_condition(obj,boundary,type) [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) % Returns the number of degrees of freedom. @@ -42,4 +40,4 @@ [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); end end -end \ No newline at end of file +end