% 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