%  ex1.m
% 
%  Example 1: 
%  Generate a triangular mesh for 
%  the swirling velocity filed on the unit square 
%  (proposed by R.LeVeque)
% 
%  Jiangguo (James) Liu
%  ColoState, Jan.-Jun. 2007

clear all;
close all;

pderect([0 1 0 1]);

%  Click icon "mesh"
%  Click icon "refine mesh"
%  Click menu "export mesh"

disp('Export mesh data, then continue...');
pause;

NumElts = size(t,2);
NumNds = size(p,2);

ConnNode = sparse(NumNds, NumNds);
vrtx = zeros(1,3);
for i=1:NumElts
   vrtx(1:3) = t(1:3,i);
   vrtx(4) = vrtx(1);
   for j=1:3
      ConnNode(vrtx(j), vrtx(j+1)) = i;
   end
end

A = triu(ConnNode,1);
B = (tril(ConnNode,-1))';
C = A+B;
NumEdges = nnz(C);
EdgeInfo = zeros(NumEdges,5);

k = 0;
for i=1:NumNds
   for j=(i+1):NumNds
      if (C(i,j)~=0)
         k = k+1;
         EdgeInfo(k,:) = [i,j,A(i,j),B(i,j),sign(A(i,j)*B(i,j))];
      end
   end
end

% figure(1);
% hold on;
% 
% x = zeros(3,1);
% y = zeros(3,1);
% 
% x = zeros(3,1);
% y = zeros(3,1);
% 
% for j=1:NumElts
%    for i=1:3
%       x(i) = p(1,t(i,j));
%       y(i) = p(2,t(i,j));
%    end
%    fill(x,y,'w');
% end
% 
% % for i=1:10
% %    plot([i*0.1 i*0.1], [0 1], 'b-');
% %    plot([0 1], [i*0.1 i*0.1], 'b-');
% % end
% 
% for j=1:NumElts
%    for i=1:3
%       x(i) = p(1,t(i,j));
%       y(i) = p(2,t(i,j));
%    end
%    fill(x,y,'w');
%    xc = (x(1)+x(2)+x(3))/3;
%    yc = (y(1)+y(2)+y(3))/3;
%    
% %    %  This is funny
% %    if j<=9
% %       xc = xc;
% %    end
% %    if j>=10 & j<=99 
% %       xc = xc-0.00625;
% %    end
% %    if j>=100 & j<=999
% %       xc = xc-0.0125;
% %    end
% %    text(xc,yc,num2str(j));
% end
% 
% print -f1 -deps ex1.eps;

%  Write the triangular mesh data to a file
fp = fopen('TrigMesh.dat','wt');
fprintf(fp, '%7d  %7d  %7d\n', NumElts, NumNds, NumEdges);
fprintf(fp, '%7d  %7d  %7d\n', 1, 1, 1);
for i=1:NumElts
   fprintf(fp, '%7d  %2d  %7d  %7d  %7d\n', i, 3, t(1:3,i));
end
for i=1:NumNds
   fprintf(fp,' %7d   %f   %f\n', i, p(:,i));
end
for i=1:NumEdges
   fprintf(fp, '%7d  %2d  %7d  %7d  %7d  %7d\n', i, EdgeInfo(i,5), EdgeInfo(i,1:4));
end
fclose(fp);

disp('Press any key to close all windows...');
pause;
close all;

return;