% 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 02 0.913545 0.4067373 0.669131 0.7431454 0.309017 0.9510575 -0.104528 0.9945226 -0.5 0.8660257 -0.809017 0.5877858 -0.978148 0.2079129 -0.978148 -0.20791210 -0.809017 -0.58778511 -0.5 -0.86602512 -0.104528 -0.99452213 0.309017 -0.95105714 0.669131 -0.74314515 0.913545 -0.40673716 -0.161265 -0.17910317 0.313878 0.22804618 -0.314083 0.34882519 0.40037 -0.29088620 0.0609951 -0.5803321 0.0617879 0.58787322 -0.587046 1.34875e-16];   function ele = get_triangles(dummy_arg)   ele=[1 10 11 162 16 18 223 10 22 94 10 16 225 11 12 206 7 22 187 21 3 48 8 9 229 8 22 710 1 19 1511 20 13 1412 6 18 2113 6 21 514 19 1 1715 19 16 2016 11 20 1617 2 17 118 16 17 1819 6 7 1820 17 16 1921 21 4 522 3 17 223 17 3 2124 20 12 1325 19 20 1426 18 17 2127 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