Mercurial > repos > public > sbplib
changeset 167:15baeb35f74e feature/grids
Merge in changes from default.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 23 Feb 2016 13:25:43 +0100 |
| parents | 7cb97c1988d9 (diff) 344bde2f9d9b (current diff) |
| children | ba1ae5b2c45e |
| files | |
| diffstat | 53 files changed, 1851 insertions(+), 1180 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Cartesian.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,75 @@ +classdef Cartesian < grid.Structured + properties + n % Number of points in the grid + d % Number of dimensions + m % Number of points in each direction + x % Cell array of vectors with node placement for each dimension. + end + + % General d dimensional grid with n points + methods + % Creates a cartesian grid given vectors conatining the coordinates + % in each direction + function obj = Cartesian(varargin) + obj.d = length(varargin); + for i = 1:obj.d + obj.x{i} = varargin{i}; + obj.m(i) = length(varargin{i}); + end + obj.n = prod(obj.m); + if obj.n == 0 + error('grid:Cartesian:EmptyGrid','Input parameter gives an empty grid.') + end + end + % n returns the number of points in the grid + function o = N(obj) + o = obj.n; + end + + % d returns the spatial dimension of the grid + function o = D(obj) + o = obj.d; + end + + % points returns a n x d matrix containing the coordianets for all points. + % points are ordered according to the kronecker product with X*Y*Z + function X = points(obj) + X = zeros(obj.n, obj.d); + + for i = 1:obj.d + if iscolumn(obj.x{i}) + c = obj.x{i}; + else + c = obj.x{i}'; + end + + m_before = prod(obj.m(1:i-1)); + m_after = prod(obj.m(i+1:end)); + + X(:,i) = kr(ones(m_before,1),c,ones(m_after,1)); + end + end + + % matrices returns a cell array with coordinates in matrix form. + % For 2d case these will have to be transposed to work with plotting routines. + function X = matrices(obj) + + if obj.d == 1 % There is no 1d matrix data type in matlab, handle special case + X{1} = reshape(obj.x{1}, [obj.m 1]); + return + end + + X = cell(1,obj.d); + for i = 1:obj.d + s = ones(1,obj.d); + s(i) = obj.m(i); + + t = reshape(obj.x{i},s); + + s = obj.m; + s(i) = 1; + X{i} = repmat(t,s); + end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/CartesianTest.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,98 @@ +function tests = CartesianTest() + tests = functiontests(localfunctions); +end + + +function testWarningEmptyGrid(testCase) + in = { + {[]}, + {[],[1]}, + {[1],[2], []}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.Cartesian(in{i}{:}),'grid:Cartesian:EmptyGrid'); + end +end + +function testN(testCase) + in = { + {[1 2 3]}, + {[1 2 3],[1 2]}, + {[1 2 3],[1 2 3]}, + {[1 2 3],[1 2 3], [1]}, + {[1 2 3],[1 2 3], [1 3 4]}, + }; + + out = [3,6,9,9,27]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.N(),out(i)); + end +end + + +function testD(testCase) + in = { + {[1 2 3]}, + {[1 2 3],[1 2]}, + {[1 2 3],[1 2 3]}, + {[1 2 3],[1 2 3], [1]}, + {[1 2 3],[1 2 3], [1 3 4]}, + }; + + out = [1,2,2,3,3]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.D(),out(i)); + end +end + + +function testPoints(testCase) + in = { + {[1 2]}, + {[1 2],[3 4]}, + {[1 2],[3 4], [5 6]}, + }; + + out = { + [[1; 2]], + [[1; 1; 2; 2],[3; 4; 3; 4]], + [[1; 1; 1; 1; 2; 2; 2; 2],[3; 3; 4; 4; 3; 3; 4; 4],[ 5; 6; 5; 6; 5; 6; 5; 6]], + }; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}); + end +end + +function testMatrices(testCase) + in = { + {[1 2]}, + {[1 2],[3 4]}, + {[1 2],[3 4], [5 6]}, + }; + + out{1}{1} = [1; 2]; + + out{2}{1} = [1, 1; 2, 2]; + out{2}{2} = [3, 4; 3, 4]; + + out{3}{1}(:,:,1) = [1, 1; 2, 2]; + out{3}{1}(:,:,2) = [1, 1; 2, 2]; + + out{3}{2}(:,:,1) = [3, 4; 3, 4]; + out{3}{2}(:,:,2) = [3, 4; 3, 4]; + + out{3}{3}(:,:,1) = [5, 5; 5, 5]; + out{3}{3}(:,:,2) = [6, 6; 6, 6]; + + for i = 1:length(in) + g = grid.Cartesian(in{i}{:}); + testCase.verifyEqual(g.matrices(),out{i}); + end +end \ No newline at end of file
--- a/+grid/Curve.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,388 +0,0 @@ -classdef Curve - properties - g - gp - transformation - end - - methods - %TODO: - % Concatenation of curves - % Subsections of curves - % Stretching of curve parameter - done for arc length. - % Curve to cell array of linesegments - - % Returns a curve object. - % g -- curve parametrization for parameter between 0 and 1 - % gp -- parametrization of curve derivative - 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 - - function n = normal(obj,t) - assert(~isempty(obj.gp),'Curve has no derivative!'); - 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. Returns 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. - - % Length of arc from parameter t0 to t1 (which may be vectors). - % Computed using derivative. - function L = arcLength(obj,t0,t1) - assert(~isempty(obj.gp),'Curve has no derivative!'); - speed = @(t) sqrt(sum(obj.gp(t).^2,1)); - L = integral_vec(speed,t0,t1); - end - - % Creates the arc length parameterization of a curve. - % N -- number of points used to approximate the arclength function - function curve = arcLengthParametrization(obj,N) - default_arg('N',100); - assert(~isempty(obj.gp),'Curve has no derivative!'); - - % Construct arcLength function using splines - tvec = linspace(0,1,N); - arcVec = obj.arcLength(0,tvec); - tFunc = spline(arcVec,tvec); % t as a function of arcLength - L = obj.arcLength(0,1); - arcPar = @(s) tFunc(s*L); - - % New function and derivative - g_new = @(t)obj.g(arcPar(t)); - gp_new = @(t) normalize(obj.gp(arcPar(t))); - curve = grid.Curve(g_new,gp_new); - - end - - % 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) - - % Computes the derivative of g: R -> R^2 using an operator D1 - function gp_out = numericalDerivative(g,D1) - m = length(D1); - t = linspace(0,1,m); - gVec = g(t)'; - gpVec = (D1*gVec)'; - - gp1_fun = spline(t,gpVec(1,:)); - gp2_fun = spline(t,gpVec(2,:)); - gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; - end - - 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 - - - -function g_norm = normalize(g0) - g1 = g0(1,:); - g2 = g0(2,:); - normalization = sqrt(sum(g0.^2,1)); - g_norm = [g1./normalization; g2./normalization]; -end - -function I = integral_vec(f,a,b) -% Wrapper around the built-in function integral that -% handles multiple limits. - - Na = length(a); - Nb = length(b); - assert(Na == 1 || Nb == 1 || Na==Nb,... - 'a and b must have same length, unless one is a scalar.'); - - if(Na>Nb); - I = zeros(size(a)); - for i = 1:Na - I(i) = integral(f,a(i),b); - end - elseif(Nb>Na) - I = zeros(size(b)); - for i = 1:Nb - I(i) = integral(f,a,b(i)); - end - else - I = zeros(size(b)); - for i = 1:Nb - I(i) = integral(f,a(i),b(i)); - end - end -end - -% Returns a function handle to the spline. -function f = spline(tval,fval,spline_order) - default_arg('spline_order',4); - [m,~] = size(tval); - assert(m==1,'Need row vectors.'); - - f_spline = spapi( optknt(tval,spline_order), tval, fval ); - f = @(t) fnval(f_spline,t); -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Curvilinear.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,10 @@ +classdef Curvilinear < grid.Structured + % General grid mapping + methods (Abstract) + % baseGrid returns the domain grid of the mapping. + g = baseGrid(obj); + + % mapping returns the mapped coordinates as a grid.Function + m = mapping(obj); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Grid.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,27 @@ +classdef Grid < handle + % General d dimensional grid with n points + methods (Abstract) + % n returns the number of points in the grid + o = N(obj) + + % d returns the spatial dimension of the grid + o = D(obj) + + % points returns a n x d matrix containing the coordinates for all points. + X = points(obj) + + % Restricts the grid function gf on obj to the subgrid g. + gf = restrictFunc(obj, gf, g) + + % Projects the grid function gf on obj to the grid g. + gf = projectFunc(obj, gf, g) + end +end + + +%% Should it be able to return a cell size aswell? For an equidistant grid this would be know +%% for other grids the constructor would have to make something up. +%% For example the grid.Cartesian constructor would take a h (1 x d) vector as an in parameter. + + +%Should define boundaries somehow, look in stitchSchemes. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Multiblock.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,19 @@ +classdef Multiblock < grid.Grid + % General multiblock grid + methods (Abstract) + % NBlocks returns the number of blocks in the grid. + o = NBlocks(obj); + + % Grid returns the ith grid in the multiblockgrid + gs = grid(obj,i); + + % Grids returns a cell array of all the grids in the multiblock grid. + gs = grids(obj); + + % Split a grid function on obj to a cell array of grid function on each block + gf = splitFunc(gf) + end +end + + +% Should define boundaries and connections between grids. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Structured.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,6 @@ +classdef Structured < grid.Grid + methods (Abstract) + % Returns the size of the grid in each dimension m = [mx my mz ...] + m = size(obj) + end +end
--- a/+grid/Ti.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,201 +0,0 @@ -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
--- a/+grid/equal_step_size.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -% 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistant.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,27 @@ +% Creates a cartesian grid of dimension length(m). +% over the doman xlim, ylim, ... +% Examples: +% g = grid.equidistant([mx, my], xlim, ylim) +% g = grid.equidistant([10, 15], {0,1}, {0,2}) +function g = equidistant(m, varargin) + if length(m) ~= length(varargin) + error('grid:equidistant:NonMatchingParameters','The number of provided dimensions do not match.') + end + + for i = 1:length(m) + if ~iscell(varargin{i}) || numel(varargin{i}) ~= 2 + error('grid:equidistant:InvalidLimits','The limits should be cell arrays with 2 elements.'); + end + + if varargin{i}{1} > varargin{i}{2} + error('grid:equidistant:InvalidLimits','The elements of the limit must be increasing.'); + end + end + + X = {}; + for i = 1:length(m) + X{i} = util.get_grid(varargin{i}{:},m(i)); + end + + g = grid.Cartesian(X{:}); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/equidistantTest.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,50 @@ +function tests = equidistantTest() + tests = functiontests(localfunctions); +end + + +function testErrorInvalidLimits(testCase) + in = { + {10,{1}}, + {10,[0,1]}, + {[10,10],{0,1},{1}}, + {[10,10],{1},{1,0}}, + {10,{1,0}}, + {[10, 5],{1,0}, {0,-1}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.equidistant(in{i}{:}),'grid:equidistant:InvalidLimits',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorNonMatchingParam(testCase) + in = { + {[],{1}}, + {[],{1},{0,1}}, + {[5,5],{0,1},{0,1},{0,1}}, + {[5,5,4],{0,1},{0,1}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.equidistant(in{i}{:}),'grid:equidistant:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + + +function testCompiles(testCase) + in = { + {5, {0,1}}, + {[3 3],{0,1},{0,2}}, + }; + + out = { + [[0; 0.25; 0.5; 0.75; 1]], + [[0; 0; 0; 0.5; 0.5; 0.5; 1; 1; 1;],[0; 1; 2; 0; 1; 2; 0; 1; 2;]], + }; + + for i = 1:length(in) + g = grid.equidistant(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/evalOn.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,37 @@ +% Takes a funciton and evaluates it on a grid to return a grid function in the +% form of a (n*k)x1 vector, where n is the number of grid points and k is the +% number of components of the function. +% g -- Grid to evaluate on. +% func -- Function to evaluate. May be a function handle or a constant. If +% it is a vector value it has to be provided as a column vector, +function gf = evalOn(g, func) + if ~isa(func, 'function_handle') + % We should have a constant. + if size(func,2) ~= 1 + error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') + end + + gf = repmat(func,[g.N, 1]); + return + end + % func should now be a function_handle + + % Get coordinates and convert to cell array for easier use as a parameter + x = g.points(); + X = {}; + for i = 1:size(x, 2) + X{i} = x(:,i); + end + + % Find the number of components + x0 = num2cell(x(1,:)); + f0 = func(x0{:}); + k = length(f0); + + if size(f0,2) ~= 1 + error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector') + end + + gf = func(X{:}); + gf = reshape(reshape(gf, [g.N, k])', [g.N*k, 1]); % Reorder so that componets sits together. +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/evalOnTest.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,118 @@ +function tests = evalOnTest() + tests = functiontests(localfunctions); +end + +function testInputConstant(testCase) + in = { + 0, + 47, + 1, + [1; 2], + }; + + out = { + [0; 0; 0], + [47; 47; 47], + [1; 1; 1], + [1; 2; 1; 2; 1; 2], + }; + + g = getTestGrid('1d'); + + for i = 1:length(in) + gf = grid.evalOn(g,in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + +function testInputScalarFunction1d(testCase) + in = { + @(x)1+x*0, + @(x)x, + @(x)x.*x, + }; + + out = { + [1; 1; 1], + [0; 1; 2], + [0; 1; 4], + }; + + g = getTestGrid('1d'); + + for i = 1:length(in) + gf = grid.evalOn(g,in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + +function testInputScalarFunction2d(testCase) + in = { + @(x,y)1+x*0, + @(x,y)x-y, + @(x,y)x./(1+y), + }; + + out = { + [1; 1; 1; 1; 1; 1; 1; 1; 1], + [0; -1; -2; 1; 0; -1; 2; 1; 0], + [0; 0; 0; 1; 1/2; 1/3; 2; 1; 2/3], + }; + + g = getTestGrid('2d'); + + for i = 1:length(in) + gf = grid.evalOn(g, in{i}); + testCase.verifyEqual(gf, out{i}); + end +end + + +function testInputVectorFunction(testCase) + g = getTestGrid('1d'); + in = @(x)[x; -2*x]; + out = [0; 0; 1; -2; 2; -4]; + + gf = grid.evalOn(g,in); + testCase.verifyEqual(gf, out); + + g = getTestGrid('2d'); + in = @(x,y)[x.^2; -2*y]; + out = [ + 0; 0; + 0; -2; + 0; -4; + 1; 0; + 1; -2; + 1; -4; + 4; 0; + 4; -2; + 4; -4; + ]; + + gf = grid.evalOn(g,in); + testCase.verifyEqual(gf, out); +end + + +function testInputErrorVectorValued(testCase) + in = { + [1,2,3], + @(x,y)[x,-y]; + }; + + g = getTestGrid('2d'); + + for i = 1:length(in) + testCase.verifyError(@()grid.evalOn(g, in{i}),'grid:evalOn:VectorValuedWrongDim',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function g = getTestGrid(d) + switch d + case '1d' + g = grid.equidistant(3,{0,2}); + case '2d' + g = grid.equidistant([3,3],{0,2},{0,2}); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/funcToMatrix.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,5 @@ +% Converts a gridfunction to a matrix +% Takes a grid function and and a structured grid. +function F = funcToMatrix(g, gf) + reshapeRowMaj(gf, g.size()); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/funcToPlotMatrix.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,5 @@ +% Converts a gridfunction to a plot matrix +% Takes a grid function and and a structured grid. +function F = funcToPlotMatrix(g, gf) + reshapeToPlotMatrix(gf, g.size()); +end \ No newline at end of file
--- a/+grid/old/concat_curve.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% 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
--- a/+grid/old/curve_discretise.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -% 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
--- a/+grid/old/curve_interp.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -% 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 - - -
--- a/+grid/old/max_h.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -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
--- a/+grid/old/min_h.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -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
--- a/+grid/old/plot_shape.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% 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
--- a/+grid/old/shape.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -% 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
--- a/+grid/old/shape_discretise.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -% 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
--- a/+grid/old/shape_linesegments.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,9 +0,0 @@ -% 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
--- a/+grid/old/triang_interp.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,132 +0,0 @@ -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
--- a/+grid/old/triang_interp_pts.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -% 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
--- a/+grid/old/triang_map.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,29 +0,0 @@ -% 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
--- a/+grid/old/triang_plot_interp.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,114 +0,0 @@ -% 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
--- a/+grid/old/triang_show.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -% 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
--- a/+grid/place_label.m Tue Feb 23 13:22:00 2016 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/Curve.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,388 @@ +classdef Curve + properties + g + gp + transformation + end + + methods + %TODO: + % Concatenation of curves + % Subsections of curves + % Stretching of curve parameter - done for arc length. + % Curve to cell array of linesegments + + % Returns a curve object. + % g -- curve parametrization for parameter between 0 and 1 + % gp -- parametrization of curve derivative + 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 + + function n = normal(obj,t) + assert(~isempty(obj.gp),'Curve has no derivative!'); + 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. Returns 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. + + % Length of arc from parameter t0 to t1 (which may be vectors). + % Computed using derivative. + function L = arcLength(obj,t0,t1) + assert(~isempty(obj.gp),'Curve has no derivative!'); + speed = @(t) sqrt(sum(obj.gp(t).^2,1)); + L = integral_vec(speed,t0,t1); + end + + % Creates the arc length parameterization of a curve. + % N -- number of points used to approximate the arclength function + function curve = arcLengthParametrization(obj,N) + default_arg('N',100); + assert(~isempty(obj.gp),'Curve has no derivative!'); + + % Construct arcLength function using splines + tvec = linspace(0,1,N); + arcVec = obj.arcLength(0,tvec); + tFunc = spline(arcVec,tvec); % t as a function of arcLength + L = obj.arcLength(0,1); + arcPar = @(s) tFunc(s*L); + + % New function and derivative + g_new = @(t)obj.g(arcPar(t)); + gp_new = @(t) normalize(obj.gp(arcPar(t))); + curve = grid.Curve(g_new,gp_new); + + end + + % 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) + + % Computes the derivative of g: R -> R^2 using an operator D1 + function gp_out = numericalDerivative(g,D1) + m = length(D1); + t = linspace(0,1,m); + gVec = g(t)'; + gpVec = (D1*gVec)'; + + gp1_fun = spline(t,gpVec(1,:)); + gp2_fun = spline(t,gpVec(2,:)); + gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; + end + + 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 + + + +function g_norm = normalize(g0) + g1 = g0(1,:); + g2 = g0(2,:); + normalization = sqrt(sum(g0.^2,1)); + g_norm = [g1./normalization; g2./normalization]; +end + +function I = integral_vec(f,a,b) +% Wrapper around the built-in function integral that +% handles multiple limits. + + Na = length(a); + Nb = length(b); + assert(Na == 1 || Nb == 1 || Na==Nb,... + 'a and b must have same length, unless one is a scalar.'); + + if(Na>Nb); + I = zeros(size(a)); + for i = 1:Na + I(i) = integral(f,a(i),b); + end + elseif(Nb>Na) + I = zeros(size(b)); + for i = 1:Nb + I(i) = integral(f,a,b(i)); + end + else + I = zeros(size(b)); + for i = 1:Nb + I(i) = integral(f,a(i),b(i)); + end + end +end + +% Returns a function handle to the spline. +function f = spline(tval,fval,spline_order) + default_arg('spline_order',4); + [m,~] = size(tval); + assert(m==1,'Need row vectors.'); + + f_spline = spapi( optknt(tval,spline_order), tval, fval ); + f = @(t) fnval(f_spline,t); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/Ti.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/equal_step_size.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/concat_curve.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/curve_discretise.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/curve_interp.m Tue Feb 23 13:25:43 2016 +0100 @@ -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 + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/max_h.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/min_h.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/plot_shape.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape_discretise.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/shape_linesegments.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_interp.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_interp_pts.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_map.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_plot_interp.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/old/triang_show.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/place_label.m Tue Feb 23 13:25:43 2016 +0100 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeRowMaj.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,12 @@ +% Reshapes a matrix as if it was stored in row major order. +function B = reshapeRowMaj(A, m) + D = length(m); + + if D == 1 + m = [m 1]; + D = 2; + end + + % Reshape and reverse order of indecies + B = permute(reshape(permute(A, ndims(A):-1:1), rot90(m,2)), D:-1:1); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeRowMajTest.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,113 @@ +function tests = reshapeRowMajTest() + tests = functiontests(localfunctions); +end + +function test1D(testCase) + in = { + {5,[1; 2; 3; 4; 5]}, + {5,[1 2 3 4 5]}, + }; + out = { + [1; 2; 3; 4; 5], + [1; 2; 3; 4; 5], + }; + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + + +function testIdentity(testCase) + in = { + {[2,2], magic(2)}, + {[3,3], magic(3)}, + {[2,3], [1 2 3; 4 5 6]}, + }; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),in{i}{2}); + end +end + +function test2D(testCase) + in = { + {[2,2],[11; 12; 21; 22]}, + {[3,2],[1 2 3; 4 5 6]}, + {[6 1],[1 2 3; 4 5 6]}, + {[1 6],[1 2 3; 4 5 6]}, + }; + + out{1}(1,1) = 11; + out{1}(1,2) = 12; + out{1}(2,1) = 21; + out{1}(2,2) = 22; + + out{2} = [1 2; 3 4; 5 6]; + out{3} = [1; 2; 3; 4; 5; 6]; + out{4} = [1 2 3 4 5 6]; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + +function test3D(testCase) + in = { + {[2, 2, 2], [111; 112; 121; 122; 211; 212; 221; 222]}, + {[8 1], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[1 8], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[2 4], cat(3,[1 2; 3 4],[5 6; 7 8])}, + {[4 2], cat(3,[1 2; 3 4],[5 6; 7 8])}, + }; + + out{1}(1,1,1) = 111; + out{1}(1,1,2) = 112; + out{1}(1,2,1) = 121; + out{1}(1,2,2) = 122; + out{1}(2,1,1) = 211; + out{1}(2,1,2) = 212; + out{1}(2,2,1) = 221; + out{1}(2,2,2) = 222; + + out{2} = [1; 5; 2; 6; 3; 7; 4; 8]; + out{3} = [1 5 2 6 3 7 4 8]; + out{4} = [1 5 2 6; 3 7 4 8]; + out{5} = [1 5; 2 6; 3 7; 4 8]; + + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end + +function testNonSquare(testCase) + in = { + {[2, 3, 4],[111; 112; 113; 114; 121; 122; 123; 124; 131; 132; 133; 134; 211; 212; 213; 214; 221; 222; 223; 224; 231; 232; 233; 234]}, + }; + out{1}(1,1,1) = 111; + out{1}(1,1,2) = 112; + out{1}(1,1,3) = 113; + out{1}(1,1,4) = 114; + out{1}(1,2,1) = 121; + out{1}(1,2,2) = 122; + out{1}(1,2,3) = 123; + out{1}(1,2,4) = 124; + out{1}(1,3,1) = 131; + out{1}(1,3,2) = 132; + out{1}(1,3,3) = 133; + out{1}(1,3,4) = 134; + out{1}(2,1,1) = 211; + out{1}(2,1,2) = 212; + out{1}(2,1,3) = 213; + out{1}(2,1,4) = 214; + out{1}(2,2,1) = 221; + out{1}(2,2,2) = 222; + out{1}(2,2,3) = 223; + out{1}(2,2,4) = 224; + out{1}(2,3,1) = 231; + out{1}(2,3,2) = 232; + out{1}(2,3,3) = 233; + out{1}(2,3,4) = 234; + for i = 1:length(in) + testCase.verifyEqual(reshapeRowMaj(in{i}{2}, in{i}{1}),out{i}); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeToPlotMatrix.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,17 @@ +% Takes a grid function and reshapes it into a matrix of shape m for plotting. +function F = reshapeToPlotMatrix(gf, m) + D = length(m); + + switch D + case 1 + F = gf; + case 2 + F = reshape(gf, rot90(m,2)); + case 3 + % After the reshape the indecies will be M(z,y,x). Plot need them to be M(y,x,z) + p = [2 3 1]; % Permuation + F = permute(reshape(gf,rot90(m,2)), p); + otherwise + error('reshapeToPlotMatrix:NotImplemented','Grid function to matrix is not implemented for dimension = %d', length(m)); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/reshapeToPlotMatrixTest.m Tue Feb 23 13:25:43 2016 +0100 @@ -0,0 +1,48 @@ +function tests = reshapeToPlotMatrixTest() + tests = functiontests(localfunctions); +end + +function test1D(testCase) + inGf = [1 2 3 4 5]'; + inM = 5; + out = [1 2 3 4 5]'; + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end + +function test2D(testCase) + x = 1:2; + y = 1:3; + + f = @(x,y) x + y*10; + + xx = [1; 1; 1; 2; 2; 2]; + yy = [1; 2; 3; 1; 2; 3]; + inGf = f(xx,yy); + + [X,Y] = meshgrid(x,y); + out = f(X,Y); + + inM = [2, 3]; + + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end + +function test3D(testCase) + x = 1:2; + y = 1:3; + z = 1:4; + + f = @(x,y,z) x + y*10 + z*100; + + xx = [repmat(1, [12, 1]); repmat(2, [12, 1])]; + yy = repmat([1; 1; 1; 1; 2; 2; 2; 2; 3; 3; 3; 3], [2, 1]); + zz = repmat([1; 2; 3; 4], [6, 1]); + inGf = f(xx,yy,zz); + + [X,Y,Z] = meshgrid(x,y,z); + out = f(X,Y,Z); + + inM = [2, 3, 4]; + + testCase.verifyEqual(reshapeToPlotMatrix(inGf, inM),out); +end \ No newline at end of file