Mercurial > repos > public > sbplib
view +parametrization/Ti.m @ 774:66eb4a2bbb72 feature/grids
Remove default scaling of the system.
The scaling doens't seem to help actual solutions. One example that fails in the flexural code.
With large timesteps the solutions seems to blow up. One particular example is profilePresentation
on the tdb_presentation_figures branch with k = 0.0005
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 18 Jul 2018 15:42:52 -0700 |
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