% Draw a piewise-linear function in two dimensions on a given triangulation. % Due to a bug in plot2svg, it can't export 3D pictures well. % Then, I have to take the 3D picture, and rotate and project it manually to 2D. % That makes the code more complicated.  function main()  %  read the triangulation from the data at the end of the code    dummy_arg = 0;    node=get_nodes(dummy_arg);      [np, k]=size(node);    ele=get_triangles(dummy_arg);   [nt, k]=size(ele);     % the function whose piecewise-linear approximation will be graphed    f=inline('0.07*(22-8*x^2-10*y^2)+0.14');     % will keep here the triangles to plot and their colors    P  = zeros(3*nt, 3);    C  = zeros(nt, 3);        % iterate through triangles, save the coordinates of all the triangles    alpha=0.3;    for i=1:nt;       u=ele(i,2); v=ele(i, 3); w=ele(i, 4);       y1=node(u, 2); x1=node(u, 3); f1=f(x1, y1);       y2=node(v, 2); x2=node(v, 3); f2=f(x2, y2);       y3=node(w, 2); x3=node(w, 3); f3=f(x3, y3);        % the color of the given triangle is chosen randomly       color = alpha*rand(1, 3)+(1-alpha)*[1 1 1];        % store the triangle and its color for the future       m = 3*i - 2;       P(m+0, 1) = x1; P(m+0, 2) = y1; P(m+0, 3) = f1;       P(m+1, 1) = x2; P(m+1, 2) = y2; P(m+1, 3) = f2;       P(m+2, 1) = x3; P(m+2, 2) = y3; P(m+2, 3) = f3;        C(i, :) = color;    end  % the "base", the domain of the piecewise linear function    P0 = P; P0(:, 3) = 0*P0(:, 3);  %  Do a rotation in 3D, then plot the projections onto the xy-plane. %  This has to be done by hand since plot2svg has trouble saving 3D graphics    a = pi/2.5; b = 0; c = 0;    Q = do_rotate(P, a, b, c);    Q0 = do_rotate(P0, a, b, c);     % sort the triangles by the third coordinate of the center of gravity (after the rotation)    R = zeros(nt, 2);    for i=1:nt       m = 3*i-2;       z1=Q(m, 3);       z2=Q(m+1, 3);       z3=Q(m+2, 3);        R(i, 1) = (z1+z2+z3)/3;       R(i, 2) = i;    end    R = sortrows(R, 1);     % plot the projection of the rotated figure and the base shape    clf; hold on; axis equal; axis off;    lw = 0.5; black = [0, 0, 0]; white = [1, 1, 1];    for i = 1:nt        j = R(i, 2);       m = 3*j-2;       fill([Q(m, 1), Q(m+1, 1) Q(m+2, 1)], [Q(m, 2), Q(m+1, 2) Q(m+2, 2)], C(i, :));              fill([Q0(m, 1), Q0(m+1, 1) Q0(m+2, 1)], [Q0(m, 2), Q0(m+1, 2) Q0(m+2, 2)], white);        plot([Q(m, 1), Q(m+1, 1) Q(m+2, 1), Q(m, 1)], [Q(m, 2), Q(m+1, 2) Q(m+2, 2), Q(m, 2)], ... 	   'linewidth', lw, 'color', black);        plot([Q0(m, 1), Q0(m+1, 1) Q0(m+2, 1), Q0(m, 1)], [Q0(m, 2), Q0(m+1, 2) Q0(m+2, 2), Q0(m, 2)], ... 	   'linewidth', lw, 'color', black);    end     % a small fix to avoid a bug with the bounding box when exporting    small = 0.1;    Sx = min(min(Q(:, 1)), min(Q0(:, 1)))-small;    Lx = max(max(Q(:, 1)), max(Q0(:, 1)))+small;    Sy = min(min(Q(:, 2)), min(Q0(:, 2)))-small;    Ly = max(max(Q(:, 2)), max(Q0(:, 2)))+small;    plot(Lx, Ly, '*', 'color', 0.99*white);    plot(Sx, Sy, '*', 'color', 0.99*white);    axis([Sx-small Lx+small, Sy-small, Ly+small])  %  export as eps and svg %  saveas(gcf, 'piecewise_linear2D_proj.eps', 'psc2')    plot2svg('piecewise_linear2D_proj.svg')     function node = get_nodes (dummy_arg)     node =[1 1 0 2 0.913545 0.406737 3 0.669131 0.743145 4 0.309017 0.951057 5 -0.104528 0.994522 6 -0.5 0.866025 7 -0.809017 0.587785 8 -0.978148 0.207912 9 -0.978148 -0.207912 10 -0.809017 -0.587785 11 -0.5 -0.866025 12 -0.104528 -0.994522 13 0.309017 -0.951057 14 0.669131 -0.743145 15 0.913545 -0.406737 16 -0.161265 -0.179103 17 0.313878 0.228046 18 -0.314083 0.348825 19 0.40037 -0.290886 20 0.0609951 -0.58033 21 0.0617879 0.587873 22 -0.587046 1.34875e-16];      function ele = get_triangles(dummy_arg)     ele=[1 10 11 16 2 16 18 22 3 10 22 9 4 10 16 22 5 11 12 20 6 7 22 18 7 21 3 4 8 8 9 22 9 8 22 7 10 1 19 15 11 20 13 14 12 6 18 21 13 6 21 5 14 19 1 17 15 19 16 20 16 11 20 16 17 2 17 1 18 16 17 18 19 6 7 18 20 17 16 19 21 21 4 5 22 3 17 2 23 17 3 21 24 20 12 13 25 19 20 14 26 18 17 21 27 14 15 19];  function Q = do_rotate(P, a, b, c)     M = [1, 0, 0; 0, cos(a), sin(a); 0 -sin(a), cos(a)]*[cos(b), 0, -sin(b); 0, 1, 0; sin(b), 0, cos(b)]...        *[cos(c), sin(c), 0; -sin(c), cos(c), 0; 0, 0, 1];     [m, n] = size(P);     Q = 0*P;     for i=1:m        X = P(i, :)';       X = M*X;              Q(i, 1) = X(1);       Q(i, 2) = X(2);       Q(i, 3) = X(3);    end