Mercurial > repos > public > sbplib
changeset 0:48b6fb693025
Initial commit.
line wrap: on
line diff
diff -r 000000000000 -r 48b6fb693025 +anim/animate.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/animate.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,33 @@ +% Calls F(t) repeatedly +% Should there be a Fsetup and a F, two function, to allow creating a plot and then updating it? +% F takes the time to generate the frame for and returns the actual time for the generated frame. +% t = F(t_r) is a function that paints a frame for time t. t is the closest time <=t_r +% it will be called for increasnig t. + +%Todo: make it catch up and produce warnings if it lags behind? Instead of just requesting the next target time +function animate(F, t, tend, time_modifier , frame_rate) + if ~exist('time_modifier') + time_modifier = 1; + end + + if ~exist('frame_rate') + frame_rate = 30; + end + + frame_time = 1/frame_rate; + dt = frame_time*time_modifier; + + animation_start = tic(); + t = F(t); + while t < tend + t = F(t + dt); + + t_left = t/time_modifier-toc(animation_start); + + pause(t_left) + end + time_to_animate = toc(animation_start); + expected_time = tend/time_modifier; + fprintf('Time to animate: %.3f\n', time_to_animate) + fprintf('Expected time : %.3f\n', expected_time) +end
diff -r 000000000000 -r 48b6fb693025 +anim/func.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/func.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,24 @@ +% Animates a function F(x,t) for x in lim from t to tend +function func(xrange,yrange,t,tend, F) + x = linspace(xrange(1),xrange(2), 200); + + fig_handle = figure; + plot_handle = plot(x,F(x,t)); + xlim(xrange) + ylim(yrange) + axis_handle = gca; + + + function t = G(t) + set(plot_handle,'YData',F(x,t)) + title(axis_handle,sprintf('T=%.3f',t)); + drawnow + end + + anim.animate(@G, t, tend) +end + + + + +
diff -r 000000000000 -r 48b6fb693025 +anim/setup_1d_plot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_1d_plot.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,29 @@ +function [update_data,figure_handle,plot_handles] = setup_1d_plot(x,y_lim,yfun) + default_arg('yfun',{@(y)y}); + + figure_handle = figure; + plot_handles(1) = plot(x,0*x); + hold on + for i = 2:length(yfun) + plot_handles(i) = plot(x,0*x); + end + hold off + + axis_handle = gca; + + xlabel('x') + ylabel('y') + xlim([x(1) x(end)]); + ylim(y_lim); + + function update(t,varargin) + if ishandle(figure_handle) && ishandle(axis_handle) + for i = 1:length(yfun) + set(plot_handles(i),'YData',yfun{i}(varargin{:})); + end + title(axis_handle,sprintf('T=%.3f',t)); + drawnow + end + end + update_data = @update; +end
diff -r 000000000000 -r 48b6fb693025 +anim/setup_2d_plot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_2d_plot.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,68 @@ +function [update_data,figure_handle] = setup_2d_plot(x,y,clim,zfun) + default_arg('zfun',@(z)z); + + Z = zeros(length(x),length(y)); + figure_handle = figure; + plot_handle = imagesc(x,y,Z); + axis_handle = gca; + + xlabel('x') + ylabel('y') + xlim([x(1) x(end)]); + ylim([y(1) y(end)]); + caxis(clim); + % axis vis3d + colormap(parula(256)) + colorbar + + function update(t,z) + Z = zfun(z); + % Z = reshape(zfun(z),length(x),length(y)); + if ishandle(plot_handle) && ishandle(axis_handle) + set(plot_handle,'CData',Z) + title(axis_handle,sprintf('T=%.3f',t)); + drawnow + end + end + update_data = @update; +end + + +% TODO +% This function is for squre grid. +% This function is for 2d image. +% Make one for 3d surface +% Make one for curvilinear grids using pcolor + + +function [update_data,figure_handle,plot_handles] = setup_1d_plot(x,y_lim,yfun) + + figure_handle = figure; + plot_handles(1) = plot(x,0*x); + hold on + for i = 2:length(yfun) + plot_handles(i) = plot(x,0*x); + end + hold off + + axis_handle = gca; + + xlabel('x') + ylabel('y') + xlim([x(1) x(end)]); + ylim(y_lim); + + function update(t,varargin) + if ishandle(figure_handle) && ishandle(axis_handle) + for i = 1:length(yfun) + set(plot_handles(i),'YData',yfun{i}(varargin{:})); + end + title(axis_handle,sprintf('T=%.3f',t)); + drawnow + end + end + update_data = @update; +end + + +
diff -r 000000000000 -r 48b6fb693025 +anim/setup_fig_mov.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_fig_mov.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,11 @@ +function save_frame = setup_fig_mov(figure_handle,dirname) + + fresh_dir(dirname) + + frame_nr = 0; + function save_frame_F() + saveas(figure_handle,sprintf('%s/%s%04d.png',dirname,'fig',frame_nr)); + frame_nr = frame_nr + 1; + end + save_frame = @save_frame_F; +end
diff -r 000000000000 -r 48b6fb693025 +anim/setup_img_mov.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_img_mov.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,18 @@ +function save_frame = setup_img_mov(dirname,filename_prefix,ulim,map,F) + default_arg('map',parula(256)); + default_arg('F',@(x)x); + + fresh_dir(dirname); + + n = size(map,1); + + frame_nr = 0; + function save_frame_F(Z) + filename = sprintf('%s/%s%04d.png',dirname,filename_prefix, frame_nr); + z = uint8(F(Z)/ulim(2)*n); + imwrite(z,map,filename); + frame_nr = frame_nr + 1; + end + + save_frame = @save_frame_F; +end
diff -r 000000000000 -r 48b6fb693025 +anim/setup_surf_plot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_surf_plot.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,30 @@ +function [update_data,figure_handle] = setup_2d_plot(x,y,z_lim,zfun) + default_arg('zfun',@(z)z); + + Z = zeros(length(y),length(x)); + figure_handle = figure; + plot_handle = surf(x,y,Z); + plot_handle.LineStyle = 'none'; + axis_handle = gca; + + xlabel('x') + ylabel('y') + xlim([x(1) x(end)]); + ylim([y(1) y(end)]); + zlim(z_lim); + caxis(z_lim); + % axis vis3d + % colormap(parula(256)) + % colorbar + + function update(t,z) + Z = zfun(z); + % Z = reshape(zfun(z),length(x),length(y)); + if ishandle(plot_handle) && ishandle(axis_handle) + set(plot_handle,'ZData',Z) + title(axis_handle,sprintf('T=%.3f',t)); + drawnow + end + end + update_data = @update; +end
diff -r 000000000000 -r 48b6fb693025 +anim/setup_time_quantity_plot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/setup_time_quantity_plot.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,28 @@ +function [update_data, plot_handles] = setup_time_quantity_plot(yfun) + default_arg('yfun',@(y)y); + + t = []; + for i = 1:length(yfun) + plot_handles(i) = line(0,0); + plot_handles(i).XData = []; + plot_handles(i).YData = []; + quantities{i} = []; + end + + axis_handle = gca; + legend() + + + function update(t_now,varargin) + if ishandle(axis_handle) + t = [t t_now]; + for j = 1:length(yfun) + quantities{j} = [quantities{j} yfun{j}(varargin{:})]; + plot_handles(j).XData = t; + plot_handles(j).YData = quantities{j}; + end + drawnow + end + end + update_data = @update; +end
diff -r 000000000000 -r 48b6fb693025 +draw/labelPoint.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+draw/labelPoint.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,21 @@ +% Takes a variable number of points as inputs and plots them with a label to the current figure. +% label_pt(p) +function label_pt(varargin) + for i = 1:length(varargin) + try + placelabel(varargin{i},inputname(i)); + catch e + error('Could not place label for input %d, ''%s''. Not a valid point',i,inputname(i)) + end + end +end + +function placelabel(pt,str) + x = pt(1); + y = pt(2); + h = line(x,y); + h.Marker = '.'; + h = text(x,y,str); + h.HorizontalAlignment = 'center'; + h.VerticalAlignment = 'bottom'; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +draw/point.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+draw/point.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,5 @@ +function l = point(p) + l = line(p(1,:),p(2,:)); + l.LineStyle = 'none'; + l.Marker = '.'; +end
diff -r 000000000000 -r 48b6fb693025 +draw/prompt_bezier.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+draw/prompt_bezier.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,30 @@ +function [C,h] = prompt_bezier(s,varargin) + default_arg('s',[]) + if ~isempty(s) + fprintf(s,varargin{:}); + end + + fprintf('# Draw bezier curve\n') + a = draw.prompt_point('Enter starting point\n'); + p = draw.point(a); + p.Color = Color.green; + p.MarkerSize = 24; + b = draw.prompt_point('Enter stopping point\n'); + p = draw.point(b); + p.Color = Color.red; + p.MarkerSize = 24; + c1 = draw.prompt_point('Enter control point 1\n'); + p = draw.point(c1); + p.Color = Color.yellow; + p.MarkerSize = 16; + c2 = draw.prompt_point('Enter control point 2\n'); + p = draw.point(c2); + p.Color = Color.yellow; + p.MarkerSize = 16; + + C = grid.Curve.bezier(a,c1,c2,b); + fprintf('C = grid.Curve.bezier([%.3g; %.3g],[%.3g; %.3g],[%.3g; %.3g],[%.3g; %.3g]);\n',a,c1,c2,b) + h = C.plot(); + uistack(h,'bottom'); + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +draw/prompt_line.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+draw/prompt_line.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,19 @@ +function [C,h] = prompt_line(s,varargin) + default_arg('s',[]) + if ~isempty(s) + fprintf(s,varargin{:}); + end + + a = draw.prompt_point('Enter starting point\n'); + p = draw.point(a); + p.Color = Color.green; + p.MarkerSize = 24; + b = draw.prompt_point('Enter stopping point\n'); + p = draw.point(b); + p.Color = Color.red; + p.MarkerSize = 24; + + C = grid.Curve.line(a,b); + h = C.plot(); + uistack(h,'bottom'); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +draw/prompt_point.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+draw/prompt_point.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,23 @@ +function [p, button] = prompt_point(s,varargin) + default_arg('s',[]) + + set(gcf,'Pointer','crosshair') + + if ~isempty(s) + fprintf(s,varargin{:}); + end + + a = gca; + + function get_point(src,event) + cp = a.CurrentPoint; + p = cp(1,1:2)'; + a.ButtonDownFcn = []; + end + + a.ButtonDownFcn = @get_point; + waitfor(a,'ButtonDownFcn', []) + + set(gcf,'Pointer','arrow') + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/Curve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Curve.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,347 @@ +classdef Curve + properties + g + gp + transformation + end + methods + %TODO: + % Errors or FD if there is no derivative function added. + + % Concatenation of curves + % Subsections of curves + % Stretching of curve paramter + % Curve to cell array of linesegments + + % Should accept derivative of curve. + function obj = Curve(g,gp,transformation) + default_arg('gp',[]); + default_arg('transformation',[]); + p_test = g(0); + assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); + + if ~isempty(transformation) + transformation.base_g = g; + transformation.base_gp = gp; + [g,gp] = grid.Curve.transform_g(g,gp,transformation); + end + + obj.g = g; + obj.gp = gp; + obj.transformation = transformation; + end + + % Made up length calculation!! Science this before actual use!! + % Calculates the length of the curve. Makes sure the longet segment used is shorter than h_max. + function [L,t] = curve_length(C,h_max) + default_arg('h_max',0.001); + g = C.g; + h = 0.1; + m = 1/h+1; + t = linspace(0,1,m); + + [p,d] = get_d(t,g); + + while any(d>h_max) + I = find(d>h_max); + + % plot(p(1,:),p(2,:),'.') + % waitforbuttonpress + + new_t = []; + for i = I + new_t(end +1) = (t(i)+t(i+1))/2; + end + t = [t new_t]; + t = sort(t); + + [p,d] = get_d(t,g); + end + + L = sum(d); + + function [p,d] = get_d(t,g) + n = length(t); + + p = zeros(2,n); + for i = 1:n + p(:,i) = g(t(i)); + end + + d = zeros(1,n-1); + for i = 1:n-1 + d(i) = norm(p(:,i) - p(:,i+1)); + end + end + end + + function n = normal(obj,t) + deriv = obj.gp(t); + normalization = sqrt(sum(deriv.^2,1)); + n = [-deriv(2,:)./normalization; deriv(1,:)./normalization]; + end + + + % Plots a curve g(t) for 0<t<1, using n points. Retruns a handle h to the plotted curve. + % h = plot_curve(g,n) + function h = plot(obj,n) + default_arg('n',100); + + t = linspace(0,1,n); + + p = obj.g(t); + + h = line(p(1,:),p(2,:)); + end + + function h= plot_normals(obj,l,n,m) + default_arg('l',0.1); + default_arg('n',10); + default_arg('m',100); + t_n = linspace(0,1,n); + + normals = obj.normal(t_n)*l; + + n0 = obj.g(t_n); + n1 = n0 + normals; + + h = line([n0(1,:); n1(1,:)],[n0(2,:); n1(2,:)]); + set(h,'Color',Color.red); + obj.plot(m); + end + + function h= show(obj,name) + p = obj.g(1/2); + n = obj.normal(1/2); + p = p + n*0.1; + + % Add arrow + + h = text(p(1),p(2),name); + h.HorizontalAlignment = 'center'; + h.VerticalAlignment = 'middle'; + + obj.plot(); + end + % Shows curve with name and arrow for direction. + + + % how to make it work for methods without returns + function p = subsref(obj,S) + %Should i add error checking here? + %Maybe if you want performance you fetch obj.g and then use that + switch S(1).type + case '()' + p = obj.g(S.subs{1}); + % case '.' + + % p = obj.(S.subs); + otherwise + p = builtin('subsref',obj,S); + % error() + end + end + + + + + %% TRANSFORMATION OF A CURVE + function D = reverse(C) + % g = C.g; + % gp = C.gp; + % D = grid.Curve(@(t)g(1-t),@(t)-gp(1-t)); + D = C.transform([],[],-1); + end + + function D = transform(C,A,b,flip) + default_arg('A',[1 0; 0 1]); + default_arg('b',[0; 0]); + default_arg('flip',1); + if isempty(C.transformation) + g = C.g; + gp = C.gp; + transformation.A = A; + transformation.b = b; + transformation.flip = flip; + else + g = C.transformation.base_g; + gp = C.transformation.base_gp; + A_old = C.transformation.A; + b_old = C.transformation.b; + flip_old = C.transformation.flip; + + transformation.A = A*A_old; + transformation.b = A*b_old + b; + transformation.flip = flip*flip_old; + end + + D = grid.Curve(g,gp,transformation); + + end + + function D = translate(C,a) + g = C.g; + gp = C.gp; + + % function v = g_fun(t) + % x = g(t); + % v(1,:) = x(1,:)+a(1); + % v(2,:) = x(2,:)+a(2); + % end + + % D = grid.Curve(@g_fun,gp); + + D = C.transform([],a); + end + + function D = mirror(C, a, b) + assert_size(a,[2,1]); + assert_size(b,[2,1]); + + g = C.g; + gp = C.gp; + + l = b-a; + lx = l(1); + ly = l(2); + + + % fprintf('Singular?\n') + + A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l); + + % function v = g_fun(t) + % % v = a + A*(g(t)-a) + % x = g(t); + + % ax1 = x(1,:)-a(1); + % ax2 = x(2,:)-a(2); + % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; + % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; + % end + + % function v = gp_fun(t) + % v = A*gp(t); + % end + + % D = grid.Curve(@g_fun,@gp_fun); + + % g = A(g-a)+a = Ag - Aa + a; + b = - A*a + a; + D = C.transform(A,b); + + end + + function D = rotate(C,a,rad) + assert_size(a, [2,1]); + assert_size(rad, [1,1]); + g = C.g; + gp = C.gp; + + + A = [cos(rad) -sin(rad); sin(rad) cos(rad)]; + + % function v = g_fun(t) + % % v = a + A*(g(t)-a) + % x = g(t); + + % ax1 = x(1,:)-a(1); + % ax2 = x(2,:)-a(2); + % v(1,:) = a(1)+A(1,:)*[ax1;ax2]; + % v(2,:) = a(2)+A(2,:)*[ax1;ax2]; + % end + + % function v = gp_fun(t) + % v = A*gp(t); + % end + + % D = grid.Curve(@g_fun,@gp_fun); + + + % g = A(g-a)+a = Ag - Aa + a; + b = - A*a + a; + D = C.transform(A,b); + end + end + + methods (Static) + function obj = line(p1, p2) + + function v = g_fun(t) + v(1,:) = p1(1) + t.*(p2(1)-p1(1)); + v(2,:) = p1(2) + t.*(p2(2)-p1(2)); + end + g = @g_fun; + + obj = grid.Curve(g); + end + + function obj = circle(c,r,phi) + default_arg('phi',[0; 2*pi]) + default_arg('c',[0; 0]) + default_arg('r',1) + + function v = g_fun(t) + w = phi(1)+t*(phi(2)-phi(1)); + v(1,:) = c(1) + r*cos(w); + v(2,:) = c(2) + r*sin(w); + end + + function v = g_fun_deriv(t) + w = phi(1)+t*(phi(2)-phi(1)); + v(1,:) = -(phi(2)-phi(1))*r*sin(w); + v(2,:) = (phi(2)-phi(1))*r*cos(w); + end + + obj = grid.Curve(@g_fun,@g_fun_deriv); + end + + function obj = bezier(p0, p1, p2, p3) + function v = g_fun(t) + v(1,:) = (1-t).^3*p0(1) + 3*(1-t).^2.*t*p1(1) + 3*(1-t).*t.^2*p2(1) + t.^3*p3(1); + v(2,:) = (1-t).^3*p0(2) + 3*(1-t).^2.*t*p1(2) + 3*(1-t).*t.^2*p2(2) + t.^3*p3(2); + end + + function v = g_fun_deriv(t) + v(1,:) = 3*(1-t).^2*(p1(1)-p0(1)) + 6*(1-t).*t*(p2(1)-p1(1)) + 3*t.^2*(p3(1)-p2(1)); + v(2,:) = 3*(1-t).^2*(p1(2)-p0(2)) + 6*(1-t).*t*(p2(2)-p1(2)) + 3*t.^2*(p3(2)-p2(2)); + end + + obj = grid.Curve(@g_fun,@g_fun_deriv); + end + + + function [g_out,gp_out] = transform_g(g,gp,tr) + A = tr.A; + b = tr.b; + flip = tr.flip; + + function v = g_fun_noflip(t) + % v = A*g + b + x = g(t); + + v(1,:) = A(1,:)*x+b(1); + v(2,:) = A(2,:)*x+b(2); + end + + function v = g_fun_flip(t) + % v = A*g + b + x = g(1-t); + + v(1,:) = A(1,:)*x+b(1); + v(2,:) = A(2,:)*x+b(2); + end + + + switch flip + case 1 + g_out = @g_fun_noflip; + gp_out = @(t)A*gp(t); + case -1 + g_out = @g_fun_flip; + gp_out = @(t)-A*gp(1-t); + end + end + + end +end
diff -r 000000000000 -r 48b6fb693025 +grid/Ti.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Ti.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,201 @@ +classdef Ti + properties + gs % {4}Curve + S % FunctionHandle(u,v) + end + + methods + % TODO function to label boundary names. + % function to find largest and smallest delta h in the grid. Maybe shouldnt live here + function obj = Ti(C1,C2,C3,C4) + obj.gs = {C1,C2,C3,C4}; + + g1 = C1.g; + g2 = C2.g; + g3 = C3.g; + g4 = C4.g; + + A = g1(0); + B = g2(0); + C = g3(0); + D = g4(0); + + function o = S_fun(u,v) + 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,:)); + 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,:)); + + o = [o1;o2]; + end + + obj.S = @S_fun; + end + + function [X,Y] = map(obj,u,v) + default_arg('v',u); + + if isscalar(u) + u = linspace(0,1,u); + end + + if isscalar(v) + v = linspace(0,1,v); + end + + S = obj.S; + + nu = length(u); + nv = length(v); + + X = zeros(nv,nu); + Y = zeros(nv,nu); + + u = rowVector(u); + v = rowVector(v); + + for i = 1:nv + p = S(u,v(i)); + X(i,:) = p(1,:); + Y(i,:) = p(2,:); + end + 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
diff -r 000000000000 -r 48b6fb693025 +grid/equal_step_size.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equal_step_size.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,4 @@ +% Calculates M so that the stepsize m/M is as close to n/M as possible +function M = equal_step_size(n,N,m) + M = round((m*(N-1)+n)/n); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/concat_curve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/concat_curve.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,12 @@ +% Concatenate two curves g1 and g2 intop +% g = concat_curve(g1,g2) +function g = concat_curve(g1,g2) + function v = g_fun(t) + if t < 1/2 + v = g1(2*t); + else + v = g2(2*t-1); + end + end + g = @g_fun; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/curve_discretise.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/curve_discretise.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,97 @@ +% Discretises the curve g with the smallest number of points such that all segments +% are shorter than h. If do_plot is true the points of the discretisation and +% the normals of the curve in those points are plotted. +% +% [t,p,d] = curve_discretise(g,h,do_plot) +% +% t is a vector of input values to g. +% p is a cector of points. +% d are the length of the segments. +function [t,p,d] = curve_discretise(g,h,do_plot) + default_arg('do_plot',false) + + n = 10; + + [t,p,d] = curve_discretise_n(g,n); + + % ni = 0; + while any(d>h) + [t,p,d] = curve_discretise_n(g,n); + n = ceil(n*d(1)/h); + % ni = ni+1; + end + + % nj = 0; + while all(d<h) + [t,p,d] = curve_discretise_n(g,n); + n = n-1; + % nj = nj+1; + end + [t,p,d] = curve_discretise_n(g,n+1); + + % fprintf('ni = %d, nj = %d\n',ni,nj); + + if do_plot + fprintf('n:%d max: %f min: %f\n', n, max(d),min(d)); + p = grid.map_curve(g,t); + figure + show(g,t,h); + end + +end + +function [t,p,d] = curve_discretise_n(g,n) + t = linspace(0,1,n); + t = equalize_d(g,t); + d = D(g,t); + p = grid.map_curve(g,t); +end + +function d = D(g,t) + p = grid.map_curve(g,t); + + d = zeros(1,length(t)-1); + for i = 1:length(d) + d(i) = norm(p(:,i) - p(:,i+1)); + end +end + +function t = equalize_d(g,t) + d = D(g,t); + v = d-mean(d); + while any(abs(v)>0.01*mean(d)) + dt = t(2:end)-t(1:end-1); + t(2:end) = t(2:end) - cumsum(dt.*v./d); + + t = t/t(end); + d = D(g,t); + v = d-mean(d); + end +end + + +function show(g,t,hh) + p = grid.map_curve(g,t); + + + + h = grid.plot_curve(g); + h.LineWidth = 2; + axis equal + hold on + h = plot(p(1,:),p(2,:),'.'); + h.Color = [0.8500 0.3250 0.0980]; + h.MarkerSize = 24; + hold off + + n = grid.curve_normals(g,t); + hold on + for i = 1:length(t) + p0 = p(:,i); + p1 = p0 + hh*n(:,i); + l = [p0, p1]; + h = plot(l(1,:),l(2,:)); + h.Color = [0.8500 0.3250 0.0980]; + end + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/curve_interp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/curve_interp.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,58 @@ +% Create a cubic spline from the points (x,y) using periodic conditions. +% g = curve_interp(x,y) +function g = curve_interp(x,y) + default_arg('x',[0 2 2 1 1 0]) + default_arg('y',[0 0 2 2 1 1]) + % solve for xp and yp + + % x(t) = at^4 + bt^2+ct+d + + % a = xp1 -2x1 + 2x0 + xp0 + % b = 3x1 -xp1 - 3x0 + 2xp0 + % c = xp0 + % d = x0 + + assert(length(x) == length(y)) + n = length(x); + A = spdiags(ones(n,1)*[2, 8, 2],-1:1,n,n); + A(n,1) = 2; + A(1,n) = 2; + + bx = zeros(n,1); + for i = 2:n-1 + bx(i) = -6*x(i-1)+6*x(i+1); + end + bx(1) = -6*x(n)+6*x(2); + bx(n) = -6*x(n-1)+6*x(1); + + by = zeros(n,1); + for i = 2:n-1 + by(i) = -6*y(i-1)+6*y(i+1); + end + by(1) = -6*y(n)+6*y(2); + by(n) = -6*y(n-1)+6*y(1); + + + xp = A\bx; + yp = A\by; + + x(end+1) = x(1); + y(end+1) = y(1); + + xp(end+1) = xp(1); + yp(end+1) = yp(1); + + function v = g_fun(t) + t = mod(t,1); + i = mod(floor(t*n),n) + 1; + t = t * n -(i-1); + X = (2*x(i)-2*x(i+1)+xp(i)+xp(i+1))*t.^3 + (-3*x(i)+3*x(i+1)-2*xp(i)-xp(i+1))*t.^2 + (xp(i))*t + x(i); + Y = (2*y(i)-2*y(i+1)+yp(i)+yp(i+1))*t.^3 + (-3*y(i)+3*y(i+1)-2*yp(i)-yp(i+1))*t.^2 + (yp(i))*t + y(i); + v = [X;Y]; + end + + g = @g_fun; +end + + +
diff -r 000000000000 -r 48b6fb693025 +grid/old/max_h.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/max_h.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,33 @@ +function [d_max, i1_max, j1_max, i2_max, j2_max] = max_h(X,Y) + ni = size(X,1); + nj = size(X,2); + d_max = 0; + + i1_max = 0; + j1_max = 0; + i2_max = 0; + j2_max = 0; + + D = {[0,-1],[1,0],[0,1],[-1,0]}; + + for i = 1:ni + for j = 1:nj + p1 = [X(i,j); Y(i,j)]; + for k = 1:length(D) + i2 = i+D{k}(1); + j2 = j+D{k}(2); + if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj + p2 = [X(i2,j2); Y(i2,j2)]; + d = norm(p2-p1); + if d > d_max; + d_max = d; + i1_max = i; + j1_max = j; + i2_max = i2; + j2_max = j2; + end + end + end + end + end +end
diff -r 000000000000 -r 48b6fb693025 +grid/old/min_h.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/min_h.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,34 @@ +function [d_min, i1_min, j1_min, i2_min, j2_min] = min_h(X,Y) + ni = size(X,1); + nj = size(X,2); + d_min = norm([X(1,1);Y(1,1)] - [X(ni,nj);Y(ni,nj)]); + + i1_min = 0; + j1_min = 0; + i2_min = 0; + j2_min = 0; + + D = {[-1,-1],[0,-1],[1,-1],[1,0],[1,1],[0,1],[-1,1],[-1,0]}; + % D = {[0,-1],[1,0],[0,1],[-1,0]}; + + for i = 1:ni + for j = 1:nj + p1 = [X(i,j); Y(i,j)]; + for k = 1:length(D) + i2 = i+D{k}(1); + j2 = j+D{k}(2); + if i2 >= 1 && i2 <= ni && j2 >= 1 && j2 <= nj + p2 = [X(i2,j2); Y(i2,j2)]; + d = norm(p2-p1); + if d < d_min; + d_min = d; + i1_min = i; + j1_min = j; + i2_min = i2; + j2_min = j2; + end + end + end + end + end +end
diff -r 000000000000 -r 48b6fb693025 +grid/old/plot_shape.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/plot_shape.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,12 @@ +% Plot a shape using n points. Returns cell array of plot handles. +% hs = plot_shape(s,n) +function hs = plot_shape(s,n) + default_arg('n',100); + + hs = {}; + hold on + for i = 1:length(s) + hs{end+1} = grid.plot_curve(s{i},n); + end + hold off +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/shape.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/shape.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,4 @@ +% Creates a shape from a number of curves. A shape is a cell array of curves. +function s = shape(varargin); + s = varargin; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/shape_discretise.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/shape_discretise.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,8 @@ +% Discretises a shape such that points on the curves are no further than h appart. +function p = shape_discretise(s,h) + p = []; + for i = 1:length(s) + [~,pt] = grid.curve_discretise(s{i},h); + p = [p, pt]; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/shape_linesegments.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/shape_linesegments.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +% Converts a shape into a cell array of linesegments shorter than h. +function l = shape_linesegments(s,h) + l = {}; + + for i = 1:length(s) + t = grid.curve_discretise(s{i},h); + l = [l, grid.curve_linesegments(s{i},t)]; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/triang_interp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/triang_interp.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,132 @@ +classdef triang_interp + properties + g1, g2 ,g3 % Curves encirling the tirangle in the positive direction. + A,B,C % The corners of the triangle + Sa, Sb, Sc % Mappings from square with different sides collapsed + end + + methods + function o = triang_interp(g1,g2,g3) + o.g1 = g1; + o.g2 = g2; + o.g3 = g3; + o.A = g1(0); + o.B = g2(0); + o.C = g3(0); + o.Sa = grid.triang_interp.square_to_triangle_interp(g2,g3,g1); + o.Sb = grid.triang_interp.square_to_triangle_interp(g3,g1,g2); + o.Sc = grid.triang_interp.square_to_triangle_interp(g1,g2,g3); + end + + + function show(o,N) + % Show the mapped meridians of the triangle. + % Might be used for the barycentric coordinates. + ma = @(t)o.Sa(1/2,1-t); + mb = @(t)o.Sb(1/2,1-t); + mc = @(t)o.Sc(1/2,1-t); + + na = @(t)o.Sa(t,1/2); + + ka = @(t)(o.g1(1-t)+o.g2(t))/2; + + h = grid.plot_curve(ma); + h.Color = Color.blue; + h = grid.plot_curve(mb); + h.Color = Color.blue; + h = grid.plot_curve(mc); + h.Color = Color.blue; + + h = grid.plot_curve(na); + h.Color = Color.red; + + h = grid.plot_curve(ka); + h.Color = Color.red; + + [a(1),a(2)] = ma(1/3); + [b(1),b(2)] = mb(1/3); + [c(1),c(2)] = mc(1/3); + + d = ka(1-1/3); + + + grid.label_pt(a,b,c,d); + + + % t = linspace(0,1,N); + % for i = 1:N + % sa = @(s)o.Sa(s,t(i)); + % sb = @(s)o.Sb(s,t(i)); + % sc = @(s)o.Sc(s,t(i)); + + % h = grid.plot_curve(sa); + % h.Color = Color.blue; + % h = grid.plot_curve(sb); + % h.Color = Color.blue; + % h = grid.plot_curve(sc); + % h.Color = Color.blue; + % end + + h = grid.plot_curve(o.g1); + h.LineWidth = 2; + h.Color = Color.red; + + h = grid.plot_curve(o.g2); + h.LineWidth = 2; + h.Color = Color.red; + + h = grid.plot_curve(o.g3); + h.LineWidth = 2; + h.Color = Color.red; + + end + + + end + + methods(Static) + % Makes a mapping from the unit square to a triangle by collapsing + % one of the sides of the squares to a corner on the triangle + % The collapsed side is mapped to the corner oposite to g1. + % This is done such that for S(s,t), S(s,1) = g1(s) + function S = square_to_triangle_interp(g1,g2,g3) + corner = grid.line_segment(g3(0),g3(0)); + S = grid.transfinite_interp(corner,g3,f(g1),f(g2)) + + % Function to flip a curve + function h = f(g) + h = @(t)g(1-t); + end + end + end + +end + +% % Return a mapping from u.v to x,y of the domain encircled by g1 g2 g3 in the the positive direction. created be using transfinite interpolation. +% function S = triang_interp(g1,g2,g3) +% A = g1(0) +% B = g2(0) +% C = g3(0) + +% function [x,y] = S_fun(u,v) +% w = sqrt((u-1)^2+v^2)/sqrt(2); % Parameter for g3 +% v = v*(1-u-v)*g1(u) + u*(1-u-v)*g2(v) + u*v*g3(w) ... +% +(1-u)*(1-v)*A+u*(1-v)*B + (1-u)*v*C; +% x = v(1); +% y = v(2); +% end +% S = @S_fun; +% end + + + +% function subsref(obj,S) +% if ~all(isnumeric(S.subs{:})) +% error('Only supports calling object with number') +% end +% if numel(S.subs{:}) > 1 +% disp('You''ve called the object with more than one argument'); +% else +% disp(['You called the object with argument = ',num2str(S.subs{:})]); +% end +% end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +grid/old/triang_interp_pts.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/triang_interp_pts.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,12 @@ +% Creates a transfinite interpolation from connecting the four points wiht straight lines. +function [S, g1, g2, g3] = triang_interp_pts(p1,p2,p3) + if size(p1) ~= [2 1] + error('p1 is strange!'); + end + + g1 = @(t)(p1 + t*(p2-p1)); + g2 = @(t)(p2 + t*(p3-p2)); + g3 = @(t)(p3 + t*(p1-p3)); + + S = grid.triang_interp(g1,g2,g3); +end
diff -r 000000000000 -r 48b6fb693025 +grid/old/triang_map.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/triang_map.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,29 @@ +% Creates a grid [X,Y] from the mapping function S at points in vectors u,v +function [X, Y] = traing_map(S,u,v) + error('not done') + if nargin == 2 + v = u; + end + + if isscalar(u) + u = linspace(0,1,u); + end + + if isscalar(v) + v = linspace(0,1,v); + end + + nu = length(u); + nv = length(v); + + X = zeros(nu,nv); + Y = zeros(nu,nv); + + for i = 1:nu + for j = 1:nv + [x,y] = S(u(i),v(j)); + X(i,j) = x; + Y(i,j) = y; + end + end +end
diff -r 000000000000 -r 48b6fb693025 +grid/old/triang_plot_interp.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/triang_plot_interp.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,114 @@ +% Plots a transfinite interpolation in x,y space using nu and nv curves along u and v axes. + + + + + + +% Plots a interp of a triangle where one the interpolation is from a square +% with one side collapsed to +function h = triang_plot_interp_kindaworking(S,n) + u = linspace(0,1,n); + v = linspace(0,1,n); + + m = 100; + m = 20; + + Xl_curves = cell(n,1); + Xr_curves = cell(n,1); + Y_curves = cell(n,1); + + + function u = wierdness(v,d,N) + if N == 0 + u = 0; + else + u = N*d./(1-v); + end + end + + + %Y curves + t = linspace(0,1,m); + for i = 1:n + x = []; y = []; + for j = 1:length(t) + [x(j),y(j)] = S(t(j),v(i)); + end + Y_curves{i} = [x', y']; + end + + + % Right and left X curves + t = linspace(0,1,m); + d = u(2); + for i = 1:n + xl = []; yl = []; + xr = []; yr = []; + N = i-1; + t = linspace(0,1-N*d,m); + for j = 1:length(t) + w = wierdness(t(j),d,N); + [xr(j),yr(j)] = S(w,t(j)); + [xl(j),yl(j)] = S(1-w,t(j)); + end + Xl_curves{i} = [xl', yl']; + Xr_curves{i} = [xr', yr']; + end + + for i = 1:n-1 + line(Xl_curves{i}(:,1),Xl_curves{i}(:,2)) + line(Xr_curves{i}(:,1),Xr_curves{i}(:,2)) + line(Y_curves{i}(:,1),Y_curves{i}(:,2)) + end +end + + + + +function h = triang_plot_interp_nonworking(S,n) + + u = linspace(0,1,n); + v = linspace(0,1,n); + + m = 100; + + X_curves = cell(n-1,1); + Y_curves = cell(n-1,1); + K_curves = cell(n-1,1); + + + t = linspace(0,1,m); + for i = 1:n-1 + x = []; y = []; + for j = find(t+u(i) <= 1) + [x(j),y(j)] = S(u(i),t(j)); + end + X_curves{i} = [x', y']; + end + + for i = 1:n-1 + x = []; y = []; + for j = find(t+v(i) <= 1) + [x(j),y(j)] = S(t(j),v(i)); + end + Y_curves{i} = [x', y']; + end + + for i = 2:n + x = []; y = []; + for j = find(t<u(i)) + [x(j),y(j)] = S(t(j), u(i)-t(j)); + end + K_curves{i-1} = [x', y']; + end + + for i = 1:n-1 + line(X_curves{i}(:,1),X_curves{i}(:,2)) + line(Y_curves{i}(:,1),Y_curves{i}(:,2)) + line(K_curves{i}(:,1),K_curves{i}(:,2)) + end + + h = -1; + % h = plot(X_curves{:},Y_curves{:},K_curves{:}); +end
diff -r 000000000000 -r 48b6fb693025 +grid/old/triang_show.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/old/triang_show.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,23 @@ +% Show a grid as a matlab figure. +function triang_show(S,n) + + ih = ishold(); + + hold on + h_grid = grid.triang_plot_interp(S,n); + h_bord = grid.triang_plot_interp(S,2); + + set(h_grid,'Color',[0 0.4470 0.7410]); + set(h_bord,'Color',[0.8500 0.3250 0.0980]); + set(h_bord,'LineWidth',2); + + % axis auto + % axis equal + % axis square + + if ih + hold on + else + hold off + end +end
diff -r 000000000000 -r 48b6fb693025 +grid/place_label.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/place_label.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,10 @@ +function place_label(pt,str,horzAl,vertAl) + default_arg('horzAl','center'); + default_arg('vertAl', 'middle'); + + x = pt(1); + y = pt(2); + h = text(x,y,str); + h.HorizontalAlignment = horzAl; + h.VerticalAlignment = vertAl; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +multiblock/equal_step_size.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/equal_step_size.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,4 @@ +% Choose the number of points on m, M, so that the step size is equal to that for N points on n. +function [M] = equal_step_size(n,N,m) + M = round((m*(N-1)+n)/n); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +multiblock/gridVector1d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/gridVector1d.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,8 @@ +function x = gridVector1d(schms) + n = length(schms); + x = cell(n,1); + + for i = 1:n + x{i} = schms{i}.x; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +multiblock/gridVector2d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/gridVector2d.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,10 @@ +function [x,y] = gridVector2d(schms) + n = length(schms); + x = cell(n,1); + y = cell(n,1); + + for i = 1:n + x{i} = schms{i}.x; + y{i} = schms{i}.y; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +multiblock/multiblockgrid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/multiblockgrid.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,70 @@ +% Creates a multi block square grid with defined boundary conditions. +% x,y defines the grid lines. Rember to think of the indexing as a matrix. Order matters! +% bc is a struct defining the boundary conditions on each side of the block. +% bc.w = {'dn',[function or value]} +function [block,conn,bound,ms] = multiblockgrid(x,y,mx,my,bc) + n = length(y)-1; % number of blocks in the y direction. + m = length(x)-1; % number of blocks in the x direction. + N = n*m; % number of blocks + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + if ~issorted(flip(y)) + error('The elements of y seem to be in the wrong order'); + end + % y = sort(y,'descend'); + + % Dimensions of blocks and number of points + block = cell(n,m); + for i = 1:n + for j = 1:m + block{i,j} = { + {x(j),x(j+1)}, {y(i+1),y(i)}; + }; + + ms{i,j} = [mx(i),my(j)]; + end + end + + % Interface couplings + conn = cell(N,N); + for i = 1:n + for j = 1:m + I = flat_index(n,i,j); + if i < n + J = flat_index(n,i+1,j); + conn{I,J} = {'s','n'}; + end + + if j < m + J = flat_index(n,i,j+1); + conn{I,J} = {'e','w'}; + end + end + end + + + % Boundary conditions + bound = cell(n,m); + for i = 1:n + if isfield(bc,'w') + bound{i,1}.w = bc.w; + end + + if isfield(bc,'e') + bound{i,n}.e = bc.e; + end + end + + for j = 1:m + if isfield(bc,'n') + bound{1,j}.n = bc.n; + end + + if isfield(bc,'s') + bound{m,j}.s = bc.s; + end + end +end +
diff -r 000000000000 -r 48b6fb693025 +multiblock/solutionVector2cell.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/solutionVector2cell.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,15 @@ +function v_cell = solutionVector2cell(schms,v,components) + if nargin == 2 + components = 1; + end + + N = length(schms); + v_cell = cell(N,1); + + i_start = 1; + for i =1:N + i_end = i_start + schms{i}.size()*components - 1; + v_cell{i} = v(i_start:i_end); + i_start = i_end+1; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +multiblock/stitchSchemes.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/stitchSchemes.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,85 @@ +% Stitch schemes together given connection matrix and BC vector. +% schmHand - function_handle to a Scheme constructor +% order - order of accuracy +% schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays +% blocks - block definitions, On whatever form the scheme expects. +% ms - grid points in each direction for each block. Ex {[10,10], [10, 20]} +% conn - connection matrix +% bound - boundary condition vector, array of structs with fields w,e,s,n +% each field with a parameter array that is sent to schm.boundary_condition +% +% Output parameters are cell arrays and cell matrices. +function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) + + n_blocks = numel(blocks); + + % Creating Schemes + for i = 1:n_blocks + if ~iscell(schmParam{i}) + param = schmParam(i); + else + param = schmParam{i}; + end + + % class(schmParam) + % class(ms) + % class(blocks) + % class(schmParam{i}) + % class(ms) + + schms{i} = schmHand(ms{i},blocks{i},order,param{:}); + end + + + % Total norm + H = cell(n_blocks,n_blocks); + for i = 1:n_blocks + H{i,i} = schms{i}.H; + end + + + + + + %% Total system matrix + + % Differentiation terms + D = cell(n_blocks,n_blocks); + for i = 1:n_blocks + D{i,i} = schms{i}.D; + end + + % Boundary penalty terms + for i = 1:n_blocks + if ~isstruct(bound{i}) + continue + end + + fn = fieldnames(bound{i}); + for j = 1:length(fn); + bc = bound{i}.(fn{j}); + if isempty(bc) + continue + end + + t = schms{i}.boundary_condition(fn{j},bc{:}); + D{i,i} = D{i,i}+t; + end + end + + % Interface penalty terms + for i = 1:n_blocks + for j = 1:n_blocks + intf = conn{i,j}; + if isempty(intf) + continue + end + + [uu,uv,vv,vu] = noname.Scheme.interface_coupling(schms{i},intf{1},schms{j},intf{2}); + D{i,i} = D{i,i} + uu; + D{i,j} = uv; + D{j,j} = D{j,j} + vv; + D{j,i} = vu; + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +noname/Discretization.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/Discretization.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,42 @@ +classdef Discretization < handle + properties (Abstract) + name %Short description + description %Longer description + order %Order of accuracy + end + + methods (Abstract) + % Prints some info about the discretisation + printInfo() + + % Return the number of DOF + n = size(obj) + + % Returns a timestepper for integrating the discretisation in time + % method is a string that states which timestepping method should be used. + % The implementation should switch on the string and deliver + % the appropriate timestepper. It should also provide a default value. + % time_align is a time that the timesteps should align with so that for some + % integer number of timesteps we end up exactly on time_align + ts = getTimestepper(obj,method,time_align) %% ??? + + % Sets up movie recording to a given file. + % saveFrame is a function_handle with no inputs that records the current state + % as a frame in the moive. + saveFrame = setupMov(obj, file) + + % Sets up a plot of the discretisation + % update is a function_handle accepting a timestepper that updates the plot to the + % state of the timestepper + [update,hand] = setupPlot(obj) + + end + + methods(Abstract,Static) + % Compare two functions u and v in the discrete l2 norm. + e = compareSolutions(u, v) + + % Compare the functions u to the analytical function g in the discrete l2 norm. + e = compareSolutionsAnalytical(u, g) + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +noname/animate.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+noname/animate.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,61 @@ +% animate(dirname,discretization,Tend, time_modifier,time_method) +% +% Example: +% animate('',discr,tend) +% animate('my_mov',discr,tend,time_mod,time_method) + +function hand = animate(dirname,discretization,Tend, time_modifier,time_method) + makemovies = ~strcmp(dirname,''); + if makemovies + dirname = ['mov/' dirname]; + end + + default_arg('Tend',5*60); + default_arg('time_modifier',1); + default_arg('time_method',[]); + + + fprintf('Animating: %s\n',discretization.name); + fprintf('Tend : %.2f\n',Tend); + fprintf('order : %d\n',discretization.order); + fprintf('m : %d\n',size(discretization)); + + fprintf('Creating time discretization'); + tic + ts = discretization.getTimestepper(time_method); + fprintf(' - done %fs\n', toc()) + + [update, figure_handle] = discretization.setupPlot(); + + if makemovies + save_frame = anim.setup_fig_mov(figure_handle,dirname); + end + + + % Initialize loop + str = ''; + % Loop function + function next_t = G(next_t) + ts.evolve(next_t); + update(ts); + % waitforbuttonpress + if makemovies + save_frame(); + end + % pause(0.1) + str = util.replace_string(str,'t = %.2f',ts.t); + + end + update(ts); + + fprintf('Using time step k = %.6f\n',ts.k) + fprintf('System size: %d\n',size(discretization)) + waitforbuttonpress + anim.animate(@G,0,Tend,time_modifier) + str = util.replace_string(str,''); + + % if makemovies + % fprintf('Generating movies...\n') + % system(sprintf('bash make_movie.sh %s',dirname)); + % end +end
diff -r 000000000000 -r 48b6fb693025 +scheme/Beam2d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Beam2d.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,259 @@ +classdef SchmBeam2d < noname.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 + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + alpha + + H % Discrete norm + Hi + H_x, H_y % Norms in the x and y directions + Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_x, Hi_y + Hix, Hiy + e_w, e_e, e_s, e_n + d1_w, d1_e, d1_s, d1_n + d2_w, d2_e, d2_s, d2_n + d3_w, d3_e, d3_s, d3_n + gamm_x, gamm_y + delt_x, delt_y + end + + methods + function obj = SchmBeam2d(m,lim,order,alpha,opsGen) + default_arg('opsGen',@sbp.Higher); + default_arg('a',1); + + if length(m) == 1 + m = [m m]; + end + + m_x = m(1); + m_y = m(2); + + xlim = lim{1}; + ylim = lim{2}; + + [x, h_x] = util.get_grid(xlim{:},m_x); + [y, h_y] = util.get_grid(ylim{:},m_y); + + ops_x = opsGen(m_x,h_x,order); + ops_y = opsGen(m_y,h_y,order); + + 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); + e_l_x = sparse(ops_x.boundary.e_1); + e_r_x = sparse(ops_x.boundary.e_m); + d1_l_x = sparse(ops_x.boundary.S_1); + d1_r_x = sparse(ops_x.boundary.S_m); + d2_l_x = sparse(ops_x.boundary.S2_1); + d2_r_x = sparse(ops_x.boundary.S2_m); + d3_l_x = sparse(ops_x.boundary.S3_1); + d3_r_x = sparse(ops_x.boundary.S3_m); + + D4_y = sparse(ops_y.derivatives.D4); + H_y = sparse(ops_y.norms.H); + Hi_y = sparse(ops_y.norms.HI); + e_l_y = sparse(ops_y.boundary.e_1); + e_r_y = sparse(ops_y.boundary.e_m); + d1_l_y = sparse(ops_y.boundary.S_1); + d1_r_y = sparse(ops_y.boundary.S_m); + d2_l_y = sparse(ops_y.boundary.S2_1); + d2_r_y = sparse(ops_y.boundary.S2_m); + d3_l_y = sparse(ops_y.boundary.S3_1); + d3_r_y = sparse(ops_y.boundary.S3_m); + + + D4 = kr(D4_x, I_y) + kr(I_x, D4_y); + + % Norms + obj.H = kr(H_x,H_y); + obj.Hx = kr(H_x,I_x); + obj.Hy = kr(I_x,H_y); + obj.Hix = kr(Hi_x,I_y); + obj.Hiy = kr(I_x,Hi_y); + obj.Hi = kr(Hi_x,Hi_y); + + % Boundary operators + obj.e_w = kr(e_l_x,I_y); + obj.e_e = kr(e_r_x,I_y); + obj.e_s = kr(I_x,e_l_y); + obj.e_n = kr(I_x,e_r_y); + obj.d1_w = kr(d1_l_x,I_y); + obj.d1_e = kr(d1_r_x,I_y); + obj.d1_s = kr(I_x,d1_l_y); + obj.d1_n = kr(I_x,d1_r_y); + obj.d2_w = kr(d2_l_x,I_y); + obj.d2_e = kr(d2_r_x,I_y); + obj.d2_s = kr(I_x,d2_l_y); + obj.d2_n = kr(I_x,d2_r_y); + obj.d3_w = kr(d3_l_x,I_y); + obj.d3_e = kr(d3_r_x,I_y); + 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; + + obj.gamm_y = h_y*ops_y.borrowing.N.S2/2; + obj.delt_y = h_y^3*ops_y.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. + % 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_e,penalty_d] = boundary_condition(obj,boundary,type,data) + default_arg('type','dn'); + default_arg('data',0); + + [e,d1,d2,d3,s,gamm,delt,halfnorm_inv] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet-neumann boundary condition + case {'dn'} + 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 = halfnorm_inv*(tau*e' + sig*d1'); + + pp_e = halfnorm_inv*tau; + pp_d = halfnorm_inv*sig; + switch class(data) + case 'double' + penalty_e = pp_e*data; + penalty_d = pp_d*data; + case 'function_handle' + penalty_e = @(t)pp_e*data(t); + penalty_d = @(t)pp_d*data(t); + otherwise + error('Wierd data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); + [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + 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 = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); + penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) + switch boundary + case 'w' + e = obj.e_w; + d1 = obj.d1_w; + d2 = obj.d2_w; + d3 = obj.d3_w; + s = -1; + gamm = obj.gamm_x; + delt = obj.delt_x; + halfnorm_inv = obj.Hix; + case 'e' + e = obj.e_e; + d1 = obj.d1_e; + d2 = obj.d2_e; + d3 = obj.d3_e; + s = 1; + gamm = obj.gamm_x; + delt = obj.delt_x; + halfnorm_inv = obj.Hix; + case 's' + e = obj.e_s; + d1 = obj.d1_s; + d2 = obj.d2_s; + d3 = obj.d3_s; + s = -1; + gamm = obj.gamm_y; + delt = obj.delt_y; + halfnorm_inv = obj.Hiy; + case 'n' + e = obj.e_n; + d1 = obj.d1_n; + d2 = obj.d2_n; + d3 = obj.d3_n; + s = 1; + gamm = obj.gamm_y; + delt = obj.delt_y; + halfnorm_inv = obj.Hiy; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = prod(obj.m); + end + + end +end
diff -r 000000000000 -r 48b6fb693025 +scheme/Euler1d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Euler1d.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,248 @@ +classdef SchmBeam2d < noname.Scheme + properties + m % Number of points in each direction, possibly a vector + N % Number of points total + h % Grid spacing + u % Grid values + x % Values of x and y for each + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + alpha + + H % Discrete norm + Hi + e_l, e_r + + end + + methods + function obj = SchmBeam2d(m,xlim,order,gamma,opsGen) + default_arg('opsGen',@sbp.Ordinary); + default_arg('gamma', 1.4); + + [x, h] = util.get_grid(xlim{:},m_x); + + ops = opsGen(m_x,h_x,order); + + I_x = speye(m); + I_3 = speye(3); + + D1 = sparse(ops.derivatives.D1); + H = sparse(ops.norms.H); + Hi = sparse(ops.norms.HI); + e_l = sparse(ops.boundary.e_1); + e_r = sparse(ops.boundary.e_m); + + D1 = kr(D1, I_3); + + % Norms + obj.H = kr(H,I_3); + + % Boundary operators + obj.e_l = kr(e_l,I_3); + obj.e_r = kr(e_r,I_3); + + obj.m = m; + obj.h = h; + obj.order = order; + + + % Man har Q_t+F_x=0 i 1D Euler, där + % q=[rho, rho*u, e]^T + % F=[rho*u, rho*u^2+p, (e+p)*u] ^T + % p=(gamma-1)*(e-rho/2*u^2); + + + %Solving on form q_t + F_x = 0 + function o = F(q) + o = [q(2); q(2).^2/q(1) + p(q); (q(3)+p(q))*q(2)/q(1)]; + end + + % Equation of state + function o = p(q) + o = (gamma-1)*(q(3)-q(2).^2/q(1)/2); + end + + + % R = + % [sqrt(2*(gamma-1))*rho , rho , rho ; + % sqrt(2*(gamma-1))*rho*u , rho*(u+c) , rho*(u-c) ; + % sqrt(2*(gamma-1))*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c, e+(gamma-1)*(e-rho*u^2/2)-rho*u*c]); + function o = R(q) + rho = q(1); + u = q(2)/q(1); + e = q(3); + + sqrt2gamm = sqrt(2*(gamma-1)); + + o = [ + sqrt2gamm*rho , rho , rho ; + sqrt2gamm*rho*u , rho*(u+c) , rho*(u-c) ; + sqrt2gamm*rho*u^2/2, e+(gamma-1)*(e-rho*u^2/2)+rho*u*c , e+(gamma-1)*(e-rho*u^2/2)-rho*u*c + ]; + end + + function o = Fx(q) + o = zeros(size(q)); + for i = 1:3:3*m + o(i:i+2) = F(q(i:i+2)); + end + end + + + + % A=R*Lambda*inv(R), där Lambda=diag(u, u+c, u-c) (c är ljudhastigheten) + % c^2=gamma*p/rho + % function o = A(rho,u,e) + % end + + + obj.D = @Fx; + obj.u = x; + obj.x = kr(x,ones(3,1)); + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % 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, alpha,data) + default_arg('alpha',0); + default_arg('data',0); + + % Boundary condition on form + % w_in = w_out + g, where g is data + + [e,s] = obj.get_boundary_ops(boundary); + + tuning = 1; % ????????????????????????? + + tau = R(q)*lambda(q)*tuning; % SHOULD THIS BE abs(lambda)????? + + function closure_fun(q,t) + q_b = e * q; + end + + function penalty_fun(q,t) + end + + + + + + % tau1 < -alpha^2/gamma + + 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 = halfnorm_inv*(tau*e' + sig*d1'); + + pp_e = halfnorm_inv*tau; + pp_d = halfnorm_inv*sig; + switch class(data) + case 'double' + penalty_e = pp_e*data; + penalty_d = pp_d*data; + case 'function_handle' + penalty_e = @(t)pp_e*data(t); + penalty_d = @(t)pp_d*data(t); + otherwise + error('Wierd data argument!') + 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,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); + [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + 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 = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); + penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) + switch boundary + case 'w' + e = obj.e_w; + d1 = obj.d1_w; + d2 = obj.d2_w; + d3 = obj.d3_w; + s = -1; + gamm = obj.gamm_x; + delt = obj.delt_x; + halfnorm_inv = obj.Hix; + case 'e' + e = obj.e_e; + d1 = obj.d1_e; + d2 = obj.d2_e; + d3 = obj.d3_e; + s = 1; + gamm = obj.gamm_x; + delt = obj.delt_x; + halfnorm_inv = obj.Hix; + case 's' + e = obj.e_s; + d1 = obj.d1_s; + d2 = obj.d2_s; + d3 = obj.d3_s; + s = -1; + gamm = obj.gamm_y; + delt = obj.delt_y; + halfnorm_inv = obj.Hiy; + case 'n' + e = obj.e_n; + d1 = obj.d1_n; + d2 = obj.d2_n; + d3 = obj.d3_n; + s = 1; + gamm = obj.gamm_y; + delt = obj.delt_y; + halfnorm_inv = obj.Hiy; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = prod(obj.m); + end + + end +end
diff -r 000000000000 -r 48b6fb693025 +scheme/Scheme.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Scheme.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,45 @@ +% Start with all matrix returns. When that works see how we should generalize to non-matrix stuff/nonlinear +classdef Scheme < handle + properties (Abstract) + m % Number of points in each direction, possibly a vector + order % Order accuracy for the approximation + + % vectors u,v,w depending on dim that gives were gridpoints are in each dimension + % vectors x,y,z containing the x,y,z values corresponding to each grid point + % matrices X,Y,Z with point coordinates as multi dimensional vectors + + D % non-stabalized scheme operator + H % Discrete norm + + % Should also containg: + % the grid points used + % the grid spacing + end + + methods (Abstract) + % 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. + m = boundary_condition(obj,boundary,type,data) + m = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + N = size(obj) % Returns the number of degrees of freedom. + + + + % Static function to calculate the coupling of two domains with an interface ??? + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +scheme/Schrodinger.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,144 @@ +classdef Schrodinger < noname.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + x % Grid + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + H % Discrete norm + M % Derivative norm + alpha + + Hi + e_l + e_r + d1_l + d1_r + gamm + end + + methods + % Solving SE in the form u_t = i*u_xx -i*V; + function obj = Schrodinger(m,xlim,order,V) + default_arg('V',0); + + [x, h] = util.get_grid(xlim{:},m); + + ops = sbp.Ordinary(m,h,order); + + obj.D2 = sparse(ops.derivatives.D2); + obj.H = sparse(ops.norms.H); + obj.Hi = sparse(ops.norms.HI); + obj.M = sparse(ops.norms.M); + obj.e_l = sparse(ops.boundary.e_1); + obj.e_r = sparse(ops.boundary.e_m); + obj.d1_l = sparse(ops.boundary.S_1); + obj.d1_r = sparse(ops.boundary.S_m); + + + if isa(V,'function_handle') + V_vec = V(x); + else + V_vec = x*0 + V; + end + + V_mat = spdiags(V,0,m,m); + + + D = 1i * obj.D2 - 1i * V; + + obj.m = m; + obj.h = h; + obj.order = order; + + obj.D = alpha*obj.D2; + obj.x = x; + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','neumann'); + default_arg('data',0); + + [e,d,s] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','d','dirichlet'} + tau = -1i* s * d; + + closure = obj.Hi*tau*e'; + + pp = obj.Hi*p; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Wierd data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + a = s/2 * 1i ; + b = - a'; + + tau = b*d_u; + sig = a*e_u; + + closure = obj.Hi * (tau*e_u' + sig*d_u'); + penalty = 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 + function [e,d,s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d = obj.d1_l; + s = -1; + case 'r' + e = obj.e_r; + d = obj.d1_r; + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +scheme/Wave.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Wave.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,175 @@ +classdef SchmWave < noname.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + x % Grid + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + H % Discrete norm + M % Derivative norm + alpha + + D2 + Hi + e_l + e_r + d1_l + d1_r + gamm + end + + methods + function obj = SchmWave(m,xlim,order,alpha) + default_arg('a',1); + [x, h] = util.get_grid(xlim{:},m); + + ops = sbp.Ordinary(m,h,order); + + obj.D2 = sparse(ops.derivatives.D2); + obj.H = sparse(ops.norms.H); + obj.Hi = sparse(ops.norms.HI); + obj.M = sparse(ops.norms.M); + obj.e_l = sparse(ops.boundary.e_1); + obj.e_r = sparse(ops.boundary.e_m); + obj.d1_l = sparse(ops.boundary.S_1); + obj.d1_r = sparse(ops.boundary.S_m); + + + obj.m = m; + obj.h = h; + obj.order = order; + + obj.alpha = alpha; + obj.D = alpha*obj.D2; + obj.x = x; + + obj.gamm = h*ops.borrowing.M.S; + + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','neumann'); + default_arg('data',0); + + [e,d,s] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','dirichlet'} + alpha = obj.alpha; + + % tau1 < -alpha^2/gamma + tuning = 1.1; + tau1 = -tuning*alpha/obj.gamm; + tau2 = s*alpha; + + p = tau1*e + tau2*d; + + closure = obj.Hi*p*e'; + + pp = obj.Hi*p; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Wierd data argument!') + end + + + % Neumann boundary condition + case {'N','neumann'} + alpha = obj.alpha; + tau1 = -s*alpha; + tau2 = 0; + tau = tau1*e + tau2*d; + + closure = obj.Hi*tau*d'; + + pp = obj.Hi*tau; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Wierd data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + tuning = 1.1; + + alpha_u = obj.alpha; + alpha_v = neighbour_scheme.alpha; + + gamm_u = obj.gamm; + gamm_v = neighbour_scheme.gamm; + + % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) + + tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; + tau2 = s_u*1/2*alpha_u; + sig1 = s_u*(-1/2); + sig2 = 0; + + tau = tau1*e_u + tau2*d_u; + sig = sig1*e_u + sig2*d_u; + + closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u'); + penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d,s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d = obj.d1_l; + s = -1; + case 'r' + e = obj.e_r; + d = obj.d1_r; + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +scheme/Wave2d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Wave2d.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,231 @@ +classdef SchmWave2d < noname.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + x,y % Grid + X,Y % Values of x and y for each grid point + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + alpha + + H % Discrete norm + Hi + H_x, H_y % Norms in the x and y directions + Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_x, Hi_y + Hix, Hiy + e_w, e_e, e_s, e_n + d1_w, d1_e, d1_s, d1_n + gamm_x, gamm_y + end + + methods + function obj = SchmWave2d(m,xlim,ylim,order,alpha) + default_arg('a',1); + + if length(m) == 1 + m = [m m]; + end + + m_x = m(1); + m_y = m(2); + + [x, h_x] = util.get_grid(xlim{:},m_x); + [y, h_y] = util.get_grid(ylim{:},m_y); + + ops_x = sbp.Ordinary(m_x,h_x,order); + ops_y = sbp.Ordinary(m_y,h_y,order); + + I_x = speye(m_x); + I_y = speye(m_y); + + D2_x = sparse(ops_x.derivatives.D2); + H_x = sparse(ops_x.norms.H); + Hi_x = sparse(ops_x.norms.HI); + M_x = sparse(ops_x.norms.M); + e_l_x = sparse(ops_x.boundary.e_1); + e_r_x = sparse(ops_x.boundary.e_m); + d1_l_x = sparse(ops_x.boundary.S_1); + d1_r_x = sparse(ops_x.boundary.S_m); + + D2_y = sparse(ops_y.derivatives.D2); + H_y = sparse(ops_y.norms.H); + Hi_y = sparse(ops_y.norms.HI); + M_y = sparse(ops_y.norms.M); + e_l_y = sparse(ops_y.boundary.e_1); + e_r_y = sparse(ops_y.boundary.e_m); + d1_l_y = sparse(ops_y.boundary.S_1); + d1_r_y = sparse(ops_y.boundary.S_m); + + D2 = kr(D2_x, I_y) + kr(I_x, D2_y); + obj.H = kr(H_x,H_y); + obj.Hx = kr(H_x,I_y); + obj.Hy = kr(I_x,H_y); + obj.Hix = kr(Hi_x,I_y); + obj.Hiy = kr(I_x,Hi_y); + obj.Hi = kr(Hi_x,Hi_y); + obj.M = kr(M_x,H_y)+kr(H_x,M_y); + obj.e_w = kr(e_l_x,I_y); + obj.e_e = kr(e_r_x,I_y); + obj.e_s = kr(I_x,e_l_y); + obj.e_n = kr(I_x,e_r_y); + obj.d1_w = kr(d1_l_x,I_y); + obj.d1_e = kr(d1_r_x,I_y); + obj.d1_s = kr(I_x,d1_l_y); + obj.d1_n = kr(I_x,d1_r_y); + + obj.m = m; + obj.h = [h_x h_y]; + obj.order = order; + + obj.alpha = alpha; + obj.D = alpha*D2; + obj.x = x; + obj.y = 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.M.S; + obj.gamm_y = h_y*ops_y.borrowing.M.S; + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','neumann'); + default_arg('data',0); + + [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','d','dirichlet'} + alpha = obj.alpha; + + % tau1 < -alpha^2/gamma + tuning = 1.1; + tau1 = -tuning*alpha/gamm; + tau2 = s*alpha; + + p = tau1*e + tau2*d; + + closure = halfnorm_inv*p*e'; + + pp = halfnorm_inv*p; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Wierd data argument!') + end + + + % Neumann boundary condition + case {'N','n','neumann'} + alpha = obj.alpha; + tau1 = -s*alpha; + tau2 = 0; + tau = tau1*e + tau2*d; + + closure = halfnorm_inv*tau*d'; + + pp = halfnorm_inv*tau; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Wierd data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + tuning = 1.1; + + alpha_u = obj.alpha; + alpha_v = neighbour_scheme.alpha; + + % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) + + tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; + tau2 = s_u*1/2*alpha_u; + sig1 = s_u*(-1/2); + sig2 = 0; + + tau = tau1*e_u + tau2*d_u; + sig = sig1*e_u + sig2*d_u; + + closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); + penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) + switch boundary + case 'w' + e = obj.e_w; + d = obj.d1_w; + s = -1; + gamm = obj.gamm_x; + halfnorm_inv = obj.Hix; + case 'e' + e = obj.e_e; + d = obj.d1_e; + s = 1; + gamm = obj.gamm_x; + halfnorm_inv = obj.Hix; + case 's' + e = obj.e_s; + d = obj.d1_s; + s = -1; + gamm = obj.gamm_y; + halfnorm_inv = obj.Hiy; + case 'n' + e = obj.e_n; + d = obj.d1_n; + s = 1; + gamm = obj.gamm_y; + halfnorm_inv = obj.Hiy; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = prod(obj.m); + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +scheme/Wave2dCurve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Wave2dCurve.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,333 @@ +classdef Wave2dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + u,v % Grid + x,y % Values of x and y for each grid point + X,Y % Grid point locations as matrices + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + c + J, Ji + a11, a12, a22 + + H % Discrete norm + Hi + H_u, H_v % Norms in the x and y directions + Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_u, Hi_v + Hiu, Hiv + e_w, e_e, e_s, e_n + du_w, dv_w + du_e, dv_e + du_s, dv_s + du_n, dv_n + gamm_u, gamm_v + end + + methods + function obj = Wave2dCurve(m,ti,order,c,opSet) + default_arg('opSet',@sbp.Variable); + + if length(m) == 1 + m = [m m]; + end + + m_u = m(1); + m_v = m(2); + m_tot = m_u*m_v; + + [u, h_u] = util.get_grid(0, 1, m_u); + [v, h_v] = util.get_grid(0, 1, m_v); + + + % Operators + ops_u = opSet(m_u,h_u,order); + ops_v = opSet(m_v,h_v,order); + + I_u = speye(m_u); + I_v = speye(m_v); + + D1_u = sparse(ops_u.derivatives.D1); + D2_u = ops_u.derivatives.D2; + H_u = sparse(ops_u.norms.H); + Hi_u = sparse(ops_u.norms.HI); + % M_u = sparse(ops_u.norms.M); + e_l_u = sparse(ops_u.boundary.e_1); + e_r_u = sparse(ops_u.boundary.e_m); + d1_l_u = sparse(ops_u.boundary.S_1); + d1_r_u = sparse(ops_u.boundary.S_m); + + D1_v = sparse(ops_v.derivatives.D1); + D2_v = ops_v.derivatives.D2; + H_v = sparse(ops_v.norms.H); + Hi_v = sparse(ops_v.norms.HI); + % M_v = sparse(ops_v.norms.M); + e_l_v = sparse(ops_v.boundary.e_1); + e_r_v = sparse(ops_v.boundary.e_m); + d1_l_v = sparse(ops_v.boundary.S_1); + d1_r_v = sparse(ops_v.boundary.S_m); + + + % Metric derivatives + [X,Y] = ti.map(u,v); + + [x_u,x_v] = gridDerivatives(X,D1_u,D1_v); + [y_u,y_v] = gridDerivatives(Y,D1_u,D1_v); + + + + J = x_u.*y_v - x_v.*y_u; + a11 = 1./J .* (x_v.^2 + y_v.^2); %% GÖR SOM MATRISER + 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)); + + dof_order = reshape(1:m_u*m_v,m_v,m_u); + + Duu = sparse(m_tot); + Dvv = sparse(m_tot); + + for i = 1:m_v + D = D2_u(a11(i,:)); + p = dof_order(i,:); + Duu(p,p) = D; + end + + for i = 1:m_u + D = D2_v(a22(:,i)); + p = dof_order(:,i); + Dvv(p,p) = D; + end + + L_12 = spdiags(a12(:),0,m_tot,m_tot); + Du = kr(D1_u,I_v); + Dv = kr(I_u,D1_v); + + Duv = Du*L_12*Dv; + Dvu = Dv*L_12*Du; + + + + obj.H = kr(H_u,H_v); + obj.Hi = kr(Hi_u,Hi_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); + + % obj.M = kr(M_u,H_v)+kr(H_u,M_v); + obj.e_w = kr(e_l_u,I_v); + obj.e_e = kr(e_r_u,I_v); + obj.e_s = kr(I_u,e_l_v); + obj.e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (obj.e_w'*Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (obj.e_e'*Dv)'; + obj.du_s = (obj.e_s'*Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (obj.e_n'*Du)'; + obj.dv_n = kr(I_u,d1_r_v); + + obj.m = m; + obj.h = [h_u h_v]; + obj.order = order; + + + obj.c = c; + obj.J = spdiags(J(:),0,m_tot,m_tot); + obj.Ji = spdiags(1./J(:),0,m_tot,m_tot); + obj.a11 = a11; + obj.a12 = a12; + obj.a22 = a22; + obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); + obj.u = u; + obj.v = v; + obj.X = X; + obj.Y = Y; + obj.x = X(:); + obj.y = Y(:); + + obj.gamm_u = h_u*ops_u.borrowing.M.S; + obj.gamm_v = h_v*ops_v.borrowing.M.S; + end + + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','neumann'); + default_arg('data',0); + + [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','d','dirichlet'} + error('not implemented') + alpha = obj.alpha; + + % tau1 < -alpha^2/gamma + tuning = 1.1; + tau1 = -tuning*alpha/gamm; + tau2 = s*alpha; + + p = tau1*e + tau2*d; + + closure = halfnorm_inv*p*e'; + + pp = halfnorm_inv*p; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Weird data argument!') + end + + + % Neumann boundary condition + case {'N','n','neumann'} + c = obj.c; + + + a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); + a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); + d = (a_n * d_n' + a_t*d_t')'; + + tau1 = -s; + tau2 = 0; + tau = c.^2 * obj.Ji*(tau1*e + tau2*d); + + closure = halfnorm_inv*tau*d'; + + pp = halfnorm_inv*tau; + switch class(data) + case 'double' + penalty = pp*data; + case 'function_handle' + penalty = @(t)pp*data(t); + otherwise + error('Weird data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary); + [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(boundary); + + F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u'; + F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v'; + + tau = 111111111111111111111111111111; + sig1 = 1/2; + sig2 = -1/2; + + penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u; + penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; + + tuning = 1.2; + + alpha_u = obj.alpha; + alpha_v = neighbour_scheme.alpha; + + % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) + + tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; + tau2 = s_u*1/2*alpha_u; + sig1 = s_u*(-1/2); + sig2 = 0; + + tau = tau1*e_u + tau2*d_u; + sig = sig1*e_u + sig2*d_u; + + closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); + penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = get_boundary_ops(obj,boundary) + switch boundary + case 'w' + e = obj.e_w; + d_n = obj.du_w; + d_t = obj.dv_w; + s = -1; + + coeff_n = obj.a11(:,1); + coeff_t = obj.a12(:,1); + case 'e' + e = obj.e_e; + d_n = obj.du_e; + d_t = obj.dv_e; + s = 1; + + coeff_n = obj.a11(:,end); + coeff_t = obj.a12(:,end); + case 's' + e = obj.e_s; + d_n = obj.dv_s; + d_t = obj.du_s; + s = -1; + + coeff_n = obj.a22(1,:)'; + coeff_t = obj.a12(1,:)'; + case 'n' + e = obj.e_n; + d_n = obj.dv_n; + d_t = obj.du_n; + s = 1; + + coeff_n = obj.a22(end,:)'; + coeff_t = obj.a12(end,:)'; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','e'} + halfnorm_inv_n = obj.Hiu; + halfnorm_inv_t = obj.Hiv; + halfnorm_t = obj.Hv; + gamm = obj.gamm_u; + case {'s','n'} + halfnorm_inv_n = obj.Hiv; + halfnorm_inv_t = obj.Hiu; + halfnorm_t = obj.Hu; + gamm = obj.gamm_v; + end + end + + function N = size(obj) + N = prod(obj.m); + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/+cdiff/cdiff.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+cdiff/cdiff.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,16 @@ +% Takes a step of +% v_tt = Dv+Ev_t+S +% +% 1/k^2 * (v_next - 2v + v_prev) = Dv + E 1/(2k)(v_next - v_prev) + S +% +function [v_next, v] = cdiff(v, v_prev, k, D, E, S) + % 1/k^2 * (v_next - 2v + v_prev) = Dv + E 1/(2k)(v_next - v_prev) + S + % ekv to + % A v_next = B v + C v_prev + S + I = speye(size(D)); + A = 1/k^2 * I - 1/(2*k)*E; + B = 2/k^2 * I + D; + C = -1/k^2 * I - 1/(2*k)*E; + + v_next = A\(B*v + C*v_prev + S); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/+rk4/get_rk4_time_step.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk4/get_rk4_time_step.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,21 @@ +% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values. +function k = get_rk4_time_step(lambda,l_type) + default_arg('l_type','complex') + + rad = abs(lambda); + if strcmp(l_type,'real') + % Real eigenvalue + % kl > -2.7852 + k = 2.7852/rad; + + elseif strcmp(l_type,'imag') + % Imaginary eigenvalue + % |kl| < 2.8284 + k = 2.8284/rad; + elseif strcmp(l_type,'complex') + % |kl| < 2.5 + k = 2.5/rad; + else + error('l_type must be one of ''real'',''imag'' or ''complex''.') + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/+rk4/rk4_stability.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk4/rk4_stability.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,58 @@ +function rk_stability() + ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); + circ = @(z)(abs(z)); + + + % contour(X,Y,z) + ax = [-4 2 -3 3]; + % hold on + fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2]) + hold on + r = 2.6; + fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r') + hold off + % contour(X,Y,z,[1,1],'b') + axis(ax) + title('4th order Runge-Kutta stability region') + xlabel('Re') + ylabel('Im') + axis equal + grid on + box on + hold off + % surf(X,Y,z) + + + rk4roots() +end + +function fcontour(f,levels,x_lim,y_lim,opt) + default_arg('opt','b') + x = linspace(x_lim(1),x_lim(2)); + y = linspace(y_lim(1),y_lim(2)); + [X,Y] = meshgrid(x,y); + mu = X+ 1i*Y; + + z = f(mu); + + contour(X,Y,z,levels,opt) + +end + + +function rk4roots() + ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); + % Roots for real evalues: + F = @(x)(abs(ruku4(x))-1); + real_x = fzero(F,-3); + + % Roots for imaginary evalues: + F = @(x)(abs(ruku4(1i*x))-1); + imag_x1 = fzero(F,-3); + imag_x2 = fzero(F,3); + + + fprintf('Real x = %f\n',real_x) + fprintf('Imag x = %f\n',imag_x1) + fprintf('Imag x = %f\n',imag_x2) +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/+rk4/rungekutta_4.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk4/rungekutta_4.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,10 @@ +% Takes one time step of size k using the rungekutta method +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = rungekutta_4(v, t , k, F) + k1 = F(v ,t ); + k2 = F(v+0.5*k*k1,t+0.5*k); + k3 = F(v+0.5*k*k2,t+0.5*k); + k4 = F(v+ k*k3,t+ k); + v = v + (1/6)*(k1+2*(k2+k3)+k4)*k; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/+rk4/rungekutta_6.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk4/rungekutta_6.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,31 @@ +% Takes one time step of size k using the rungekutta method +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = rungekutta_6(v, t , k, F) + s = 7 + k = zeros(length(v),s) + a = zeros(7,6); + c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; + b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12]; + a = [ + 0, 0, 0, 0, 0, 0; + 4/7, 0, 0, 0, 0, 0; + 115/112, -5/16, 0, 0, 0, 0; + 589/630, 5/18, -16/45, 0, 0, 0; + 229/1200 - 29/6000*sqrt(5), 119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5), 0, 0; + 71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4, 0; + -49/480 + 43/160*sqrt(5), -425/96 + 51/32*sqrt(5), 52/15 - 4/5*sqrt(5), -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5); + ] + + for i = 1:s + u = v + for j = 1: i-1 + u = u + h*a(i,j) * k(:,j) + end + k(:,i) = F(t+c(i)*k,u) + end + + for i = 1:s + v = v + k*b(i)*k(:,i) + end +end
diff -r 000000000000 -r 48b6fb693025 +time/Cdiff.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Cdiff.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,52 @@ +classdef Cdiff < time.Timestepper + properties + D + E + S + k + t + v + v_prev + n + end + + + methods + function obj = Cdiff(D, E, S, k, t0, v, v_prev) + m = size(D,1); + default_arg('E',sparse(m,m)); + default_arg('S',sparse(m,1)); + + if ~(issparse(D) && issparse(E) && issparse(S)) + warning('One of the matrices D, E, S is not sparse!') + print_issparse(D) + print_issparse(E) + print_issparse(S) + end + + obj.D = D; + obj.E = E; + obj.S = S; + obj.k = k; + obj.t = t0+k; + obj.v = v; + obj.v_prev = v_prev; + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function [vt,t] = getVt(obj) + vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) + t = obj.t; + end + + function obj = step(obj) + [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D, obj.E, obj.S); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/CdiffNonlin.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/CdiffNonlin.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,52 @@ +classdef Cdiff < time.Timestepper + properties + D + E + S + k + t + v + v_prev + n + end + + + methods + function obj = Cdiff(D, E, S, k, t0, v, v_prev) + m = size(D,1); + default_arg('E',sparse(m,m)); + default_arg('S',sparse(m,1)); + + if ~(issparse(D) && issparse(E) && issparse(S)) + warning('One of the matrices D, E, S is not sparse!') + print_issparse(D) + print_issparse(E) + print_issparse(S) + end + + obj.D = D; + obj.E = E; + obj.S = S; + obj.k = k; + obj.t = t0+k; + obj.v = v; + obj.v_prev = v_prev; + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function [vt,t] = getVt(obj) + vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) + t = obj.t; + end + + function obj = step(obj) + [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D, obj.E, obj.S); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/Ode45.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Ode45.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,67 @@ +classdef Ode45 < time.Timestepper + properties + F + k + t + w + m + D + E + S + M + C + n + end + + + methods + function obj = Ode45(D, E, S, k, t0, v0, v0t) + obj.D = D; + obj.E = E; + obj.S = S; + obj.m = length(v0); + + I = speye(obj.m); + O = sparse(obj.m,obj.m); + obj.M = [O, I; D, E*I]; % Multiply with I to allow 0 as input. + + if S == 0 + obj.C = zeros(2*obj.m,1); + else + obj.C = [zeros(obj.m,1), S]; + end + + obj.k = k; + obj.t = t0; + obj.w = [v0; v0t]; + + obj.F = @(w,t)(obj.M*w + obj.C); + end + + function [v,t] = getV(obj) + v = obj.w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + vt = obj.w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + [t,w] = ode45(@(t,w)(obj.F(w,t)),[obj.t obj.t+obj.k],obj.w); + + obj.t = t(end); + obj.w = w(end,:)'; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/Rk4SecondOrderNonlin.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rk4SecondOrderNonlin.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,74 @@ +classdef Rk4SecondOrderNonlin < time.Timestepper + properties + F + k + t + w + m + + D + E + S + + n + end + + + methods + function obj = Rk4SecondOrderNonlin(D, E, S, k, t0, v0, v0t) + + if S == 0 + S = @(v,t)0; + end + + if E == 0 + E = @(v,t)0; + end + + obj.k = k; + obj.t = t0; + obj.w = [v0; v0t]; + + m = length(v0); + function wt = F(w,t) + v = w(1:m); + vt = w(m+1:end); + + % Def: w = [v; vt] + wt(1:m,1) = vt; + wt(m+1:2*m,1) = D(v)*v + E(v)*vt + S(v,t); + + end + + obj.F = @F; + obj.D = D; + obj.E = E; + obj.S = S; + obj.m = m; + end + + function [v,t] = getV(obj) + v = obj.w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + vt = obj.w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/Rungekutta4.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rungekutta4.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,50 @@ +classdef Rungekutta4 < time.Timestepper + properties + D + S + F + k + t + v + m + n + end + + + methods + function obj = Rungekutta4(D, S, k, t0, v0) + obj.D = D; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.m = length(v0); + + if S == 0 + obj.S = zeros(obj.m,1); + else + obj.S = S; + end + + obj.F = @(v,t)(obj.D*v + obj.S); + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/Rungekutta4SecondOrder.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rungekutta4SecondOrder.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,65 @@ +classdef Rungekutta4SecondOrder < time.Timestepper + properties + F + k + t + w + m + D + E + S + M + C + n + end + + + methods + function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) + obj.D = D; + obj.E = E; + obj.S = S; + obj.m = length(v0); + + I = speye(obj.m); + O = sparse(obj.m,obj.m); + obj.M = [O, I; D, E*I]; % Multiply with I to allow 0 as input. + + if S == 0 + obj.C = zeros(2*obj.m,1); + else + obj.C = [zeros(obj.m,1), S]; + end + + obj.k = k; + obj.t = t0; + obj.w = [v0; v0t]; + + obj.F = @(w,t)(obj.M*w + obj.C); + end + + function [v,t] = getV(obj) + v = obj.w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + vt = obj.w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +time/Timestepper.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Timestepper.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,81 @@ +classdef Timestepper < handle + properties (Abstract) + t + k + n + end + + methods (Abstract) + [v,t] = getV(obj) + obj = step(obj) + end + + + methods + function [v,t] = stepN(obj,n,progress_bar) + + if progress_bar && n > 1500 + n1000 = floor(n/1000); + + s = util.replace_string('',' %d %%',0); + for i=1:n + obj.step(); + if mod(i,n1000) == 0 + s = util.replace_string(s,' %.2f %%',i/n*100); + end + end + s = util.replace_string(s,''); + else + for i=1:n + obj.step(); + end + end + v = obj.getV; + t = obj.t; + end + + function [v,t] = evolve(obj, tend, progress_bar) + default_arg('progress_bar',false) + if progress_bar + obj.evolve_with_progress(tend); + else + obj.evolve_without_progress(tend); + end + v = obj.getV; + t = obj.t; + end + + function evolve_with_progress(obj, tend) + dt = tend-obj.t; + n = floor(dt/obj.k); + n1000 = floor(n/1000); + + i = 0; + s = util.replace_string('',' %d %%',0); + while obj.t < tend - obj.k/100 + obj.step(); + + i = i + 1; + if mod(i,n1000) == 0 + s = util.replace_string(s,' %.2f %%',i/n*100); + end + end + + % if t < tend + % v = rk4.rungekutta_4(v, t, tend-t,F); + % end + + s = util.replace_string(s,''); + end + + function evolve_without_progress(obj, tend) + while obj.t < tend - obj.k/100 + obj.step(); + end + + % if t < tend + % v = rk4.rungekutta_4(v, t, tend-t,F); + % end + end + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/calc_borrowing.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/calc_borrowing.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,71 @@ + +m = 30; +h = 1; + + +%% 4th order non-compatible +[H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4(m,h); +S1 = S_1*S_1' + S_m*S_m'; +S2 = S2_1*S2_1' + S2_m*S2_m'; +S3 = S3_1*S3_1' + S3_m*S3_m'; + +alpha_II = util.matrixborrow(M4, h*S2 ); +alpha_III = util.matrixborrow(M4, h^3*S3); +fprintf('4th order non-compatible\n') +fprintf('alpha_II: %.10f\n',alpha_II) +fprintf('alpha_III: %.10f\n',alpha_III) +fprintf('\n') + + +%% 6th order non-compatible +[H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6(m,h); +S1 = S_1*S_1' + S_m*S_m'; +S2 = S2_1*S2_1' + S2_m*S2_m'; +S3 = S3_1*S3_1' + S3_m*S3_m'; + +alpha_II = util.matrixborrow(M4, h*S2 ); +alpha_III = util.matrixborrow(M4, h^3*S3); +fprintf('6th order non-compatible\n') +fprintf('alpha_II: %.10f\n',alpha_II) +fprintf('alpha_III: %.10f\n',alpha_III) +fprintf('\n') + + +%% 2nd order compatible +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher2_compatible(m,h); +S1 = S_1*S_1' + S_m*S_m'; +S2 = S2_1*S2_1' + S2_m*S2_m'; +S3 = S3_1*S3_1' + S3_m*S3_m'; + +alpha_II = util.matrixborrow(M4, h*S2 ); +alpha_III = util.matrixborrow(M4, h^3*S3); +fprintf('2nd order compatible\n') +fprintf('alpha_II: %.10f\n',alpha_II) +fprintf('alpha_III: %.10f\n',alpha_III) +fprintf('\n') + + +%% 4th order compatible +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4_compatible(m,h); +S1 = S_1*S_1' + S_m*S_m'; +S2 = S2_1*S2_1' + S2_m*S2_m'; +S3 = S3_1*S3_1' + S3_m*S3_m'; + +alpha_II = util.matrixborrow(M4, h*S2 ); +alpha_III = util.matrixborrow(M4, h^3*S3); +fprintf('4th order compatible\n') +fprintf('alpha_II: %.10f\n',alpha_II) +fprintf('alpha_III: %.10f\n',alpha_III) +fprintf('\n') + +%% 6th order compatible +[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6_compatible(m,h); +S1 = S_1*S_1' + S_m*S_m'; +S2 = S2_1*S2_1' + S2_m*S2_m'; +S3 = S3_1*S3_1' + S3_m*S3_m'; + +alpha_II = util.matrixborrow(M4, h*S2 ); +alpha_III = util.matrixborrow(M4, h^3*S3); +fprintf('6th order compatible\n') +fprintf('alpha_II: %.10f\n',alpha_II) +fprintf('alpha_III: %.10f\n',alpha_III)
diff -r 000000000000 -r 48b6fb693025 +util/get_convergence.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/get_convergence.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,5 @@ +function q = get_convergence(e,h) + for i = 1:length(e)-1 + q(i) = log(e(i)/e(i+1))/log(h(i+1)/h(i)); + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/get_grid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/get_grid.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,5 @@ +function [x,h] = get_grid(a,b,m) + % TODO: allow the interval to be a vector + x = linspace(a,b,m)'; + h = (b-a)/(m-1); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/matrixborrow.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/matrixborrow.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,15 @@ +function alpha = matrixborrow(A,B) + + function y = neg(x) + y = min(eig(A-x*B)); + end + + + x0 = 1; + while neg(x0) > -1 + x0 = 2*x0; + end + + opt = optimset('Display','none'); + alpha = fsolve(@neg, x0,opt); +end
diff -r 000000000000 -r 48b6fb693025 +util/plotf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/plotf.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,11 @@ +function plotf(f,a,b,opt) + if ~exist('opt') + opt = ''; + end + + x = linspace(a,b); + for i = 1:length(x) + y(i) = f(x(i)); + end + plot(x,y,opt) +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/plotfi.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/plotfi.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,12 @@ +% runs f for each value in v and plots the result. +function plotfv(f,v,opt) + + if ~exist('opt') + opt = ''; + end + + for i = 1:length(v) + y(i) = f(v(i)); + end + plot(v,y,opt) +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/replace_string.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/replace_string.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +% Replaces string old on the console with new. +% new may be a format string like the one accepted by fprintf +% old has to be the last string printed. +function s = replace_string(old,new,varargin) + reverseStr = repmat(sprintf('\b'), 1, length(old)); + blankStr = repmat(sprintf(' ') , 1, length(old)); + fprintf([reverseStr,blankStr,reverseStr, new],varargin{:}); + s = sprintf(new, varargin{:}); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 +util/semidiscrete2function.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/semidiscrete2function.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,13 @@ +function F = semidiscrete2function(D,S) + % TODO: Find threshhold when sparse makes sense. + % TODO: Why does having the spase block slow down convergence.m (or does it?) + if size(D,1) > 500 + D = sparse(D); + end + + if ~isa(S,'function_handle') + F = @(w,t)(D*w + S); + else + F = @(w,t)(D*w + S(t)); + end +end
diff -r 000000000000 -r 48b6fb693025 +util/tt2t.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/tt2t.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,28 @@ +% Converts a semidiscretized PDE from second order in time to first order in time. +% v_tt = Dv + S +% becomes +% w_t = Mw + C +% where +% w = [ v ; v_t] +% and +% M = [0 I; +% D 0]; +% C = [0; +% S]; +function [M,C] = tt2t(D,S) + default_arg('S',sparse(size(D))) + time_dependent_bc = isa(S,'function_handle'); + + I = eye(size(D)); + O = zeros(size(D)); + + M = [O I; + D O]; + + if ~time_dependent_bc + C = [zeros(size(S)); S]; + else + o = zeros(size(S(0))); + C = @(t)([o; S(t)]); + end +end
diff -r 000000000000 -r 48b6fb693025 +util/unique_filename.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+util/unique_filename.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +% Creates a unique filename on the form base_name.x.suffix where x is some number. +function filename = unique_filename(base_name, suffix) + filename = strcat(base_name,suffix); + i = 1; + while exist(filename) + filename = strcat(base_name,'.',num2str(i),suffix); + i = i+1; + end +end
diff -r 000000000000 -r 48b6fb693025 Color.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Color.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,11 @@ +classdef Color + properties(Constant) + blue = [0.000 0.447 0.741]; + red = [0.850 0.325 0.098]; + yellow = [0.929 0.694 0.125]; + purple = [0.494 0.184 0.556]; + green = [0.466 0.674 0.188]; + lightblue = [0.301 0.745 0.933]; + darkred = [0.635 0.078 0.184]; + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 assert_size.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assert_size.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,16 @@ +% Assert that array A has the size s. +function assert_size(A,s) + errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), format_vector(s), format_vector(size(A))); + assert(all(size(A) == s),errmsg); +end + +function str = format_vector(a) + l = length(a); + str = sprintf('[%d',a(1)); + + for i = 2:l + str = [str sprintf(', %d',a(i))]; + end + + str = [str ']']; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 cell2sparse.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cell2sparse.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,48 @@ +function A = cell2sparse(C) + n = row_height(C); + m = col_width(C); + + N = sum(n); + M = sum(m); + + A = sparse(N,M); + + n_ind = [0 cumsum(n)]; + m_ind = [0 cumsum(m)]; + + for i = 1:size(C,1) + for j = 1:size(C,2) + if ~has_matrix(C{i,j}) + continue + end + A(n_ind(i)+1:n_ind(i+1),m_ind(j)+1:m_ind(j+1)) = C{i,j}; + end + end + +end + +function m = col_width(C) + for j = 1:size(C,2) + for i = 1:size(C,1) + if ~has_matrix(C{i,j}) + continue + end + m(j) = size(C{i,j},2); + end + end +end + +function n = row_height(C) + for i = 1:size(C,1) + for j = 1:size(C,2) + if ~has_matrix(C{i,j}) + continue + end + n(i) = size(C{i,j},1); + end + end +end + +function b = has_matrix(c) + b = ~(isempty(c) || (numel(c)==1 && c == 0)); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 default_arg.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/default_arg.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,5 @@ +function default_arg(s, val) + if evalin('caller',sprintf('~exist(''%s'',''var'') || isempty(%s)',s,s)) + assignin('caller',s,val) + end +end
diff -r 000000000000 -r 48b6fb693025 find_elements.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/find_elements.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,17 @@ +% I = find_elements(a,b) +% Finds the index of elements a in b. +% a and b have to be in the same order. +function I = find_elements(a,b) + I = []; + + j = 1; + for i = 1:length(a) + while b(j) ~= a(i) + j = j + 1; + end + I(end+1) = j; + j = j+1; + end + + assert(length(I) == length(a),'Expected %d but got %d elements',length(a),length(I)) +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 flat_index.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/flat_index.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,3 @@ +function I = flat_index(n,i,j) + I = i + (j - 1)*n; +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 fresh_dir.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fresh_dir.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,12 @@ +function fresh_dir(dirname) + if isdir(dirname) + fprintf('Directory %s already exists.\n',dirname); + yN = input('Delete its content and continue? [y/N]: ','s'); + if strcmp(yN,'y') + system(['rm -rf ' dirname]); + else + error('Can''t use directory.') + end + end + mkdir(dirname); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 gridDerivatives.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gridDerivatives.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,15 @@ +% Calculates derivatives in all directions of function. +% Uses multi-dim vector form. +% gridDerivatives(F,Dx,Dy) +% gridDerivatives(F,Dx,Dy,Dz) +function varargout = gridDerivatives(F, varargin) + assert(length(varargin) == ndims(F)); + + switch ndims(F) + case 2 + varargout{1} = (varargin{1}*F')'; + varargout{2} = varargin{2}*F; + otherwise + error('Not implemented for ndims(F) = %d',ndims(F)); + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 kr.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/kr.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +function A = kr(varargin) + n = nargin; + + A = 1; + + for i = 1:n + A = kron(A,varargin{i}); + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 printSize.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/printSize.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,4 @@ +function printSize(A) + s = size(A); + fprintf('%8s has size: [%d, %d]\n',inputname(1),s(1),s(2)); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 print_issparse.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/print_issparse.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +function print_issparse(A) + b = issparse(A); + if b + s = 'true'; + else + s = 'false'; + end + fprintf('%8s is sparse: %s\n',inputname(1),s); +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 rowVector.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rowVector.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,9 @@ +function r = rowVector(v) + if size(v,1) == 1 + r = v + elseif size(v,2) == 1 + r = v'; + else + error('v is not a matrix'); + end +end \ No newline at end of file
diff -r 000000000000 -r 48b6fb693025 saveeps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/saveeps.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,60 @@ +% Saves a figure to an .eos file with corrected bounding box. +function saveeps(handle,filename) + if length(filename) < 4 || ~strcmp(filename(end-3:end),'.eps') + filename = [filename '.eps']; + end + + handle_units = handle.Units; % Save the current units to be able to restore + + % Copy size of figure in centimeters to a place where saveas will honor it + handle.Units = 'centimeters'; + handle.PaperUnits = 'centimeters'; + handle.PaperPosition(3:4) = handle.Position(3:4); + + % Save as a bugged eps file. + saveas(handle,filename,'epsc'); + + handle.Units = handle_units; % Restore the old units + + % Correct the buggy eps file + correct_stupid_matlab_bug(filename); +end + +% Corrects the format of an eps file so that the bounding box is defined at the top of the file instead of +% at the bottom. +function correct_stupid_matlab_bug(filename) + contents = fileread(filename); + lines = strsplit(contents,'\n'); + + % Find the line + pagel = findPrefix(lines,'%%Pages:'); + boundl = findPrefix(lines,'%%BoundingBox:'); + + if ~(length(pagel) == 2 && length(boundl) == 2) + error('Undexpected number of found lines'); + end + + if ~(strcmp(lines{pagel(1)},'%%Pages: (atend)') && strcmp(lines{boundl(1)},'%%BoundingBox: (atend)')) + error('Does the file really contain the error?'); + end + + % Overwrite the nasty lines with the nice ones. + lines{pagel(1)} = lines{pagel(2)}; + lines{boundl(1)} = lines{boundl(2)}; + + % Delete the duplicates + lines(pagel(2)) = []; + lines(boundl(2)) = []; + + + %Rewrite the file + contents = strjoin(lines,'\n'); + + fh = fopen(filename,'w'); + fprintf(fh, '%s',contents); + fclose(fh); +end + +function I = findPrefix(lines, prefix) + I = find(strncmp(lines,prefix,length(prefix))); +end
diff -r 000000000000 -r 48b6fb693025 sortd.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sortd.m Thu Sep 17 10:12:50 2015 +0200 @@ -0,0 +1,11 @@ +% Sorts x with respect to the function d +% d is a function that accepts an element of x and returns a scalar. +function y = sortd(d,x) + v = []; + for i = 1:length(x) + v(i) = d(x(i)); + end + + [~,I] = sort(v); + y = x(I); +end \ No newline at end of file