Mercurial > repos > public > sbplib
view +parametrization/Ti.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Thu, 17 Jan 2019 10:25:06 +0100 |
parents | edb1d60b0b77 |
children |
line wrap: on
line source
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) if isrow(u) && isrow(v) flipped = false; else flipped = true; u = u'; v = v'; end x1 = g1(u); x2 = g2(v); x3 = g3(1-u); x4 = g4(1-v); o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ... -((1-u).*(1-v).*A(1,:) + u.*(1-v).*B(1,:) + u.*v.*C(1,:) + (1-u).*v.*D(1,:)); 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,:)); if ~flipped o = [o1;o2]; else o = [o1'; o2']; end end obj.S = @S_fun; end % Does this funciton make sense? % Should it always be eval? 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 % Evaluate S for each pair of u and v, % Return same shape as u function [x, y] = eval(obj, u, v) x = zeros(size(u)); y = zeros(size(u)); for i = 1:numel(u) p = obj.S(u(i), 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.border = obj.plot(2,2); set(h.border,'Color',[0.8500 0.3250 0.0980]); set(h.border,'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 = parametrization.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 = parametrization.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 = parametrization.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 = parametrization.Ti(new_gs{:}); end end methods(Static) function obj = points(p1, p2, p3, p4) g1 = parametrization.Curve.line(p1,p2); g2 = parametrization.Curve.line(p2,p3); g3 = parametrization.Curve.line(p3,p4); g4 = parametrization.Curve.line(p4,p1); obj = parametrization.Ti(g1,g2,g3,g4); end function obj = rectangle(a, b) p1 = a; p2 = [b(1), a(2)]; p3 = b; p4 = [a(1), b(2)]; obj = parametrization.Ti.points(p1,p2,p3,p4); end % Like the constructor but allows inputing line curves as 2-cell arrays: % example: parametrization.Ti.linesAndCurves(g1, g2, {a, b} g4) function obj = linesAndCurves(C1, C2, C3, C4) C = {C1, C2, C3, C4}; c = cell(1,4); for i = 1:4 if ~iscell(C{i}) c{i} = C{i}; else c{i} = parametrization.Curve.line(C{i}{:}); end end obj = parametrization.Ti(c{:}); 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); parametrization.place_label(pc,str); parametrization.place_label(pw,'w'); parametrization.place_label(pe,'e'); parametrization.place_label(ps,'s'); parametrization.place_label(pn,'n'); end end end end