function [xc, yc]=contmesh(L,meshele,meshnod,fieldnod,fieldcontours);
%  contmesh(L,meshele,meshnod,fieldnod,fieldcontours);
%  plots contours of triangular grid 
%  fieldnod is data defined on nodes
%
%  speedall = sqrt(uall.*uall + vall.*vall);
% contmesh(L,meshele,meshnod,speed,[0:.25: 5])

fieldnod=fieldnod(:);



x = meshnod(:,2);
y = meshnod(:,3);

field3 = meshele(:,2:4);
field3(:)=fieldnod(meshele(:,2:4));


x3 = meshele(:,2:4);
x3(:)=x(meshele(:,2:4));
y3 = meshele(:,2:4);
y3(:)=y(meshele(:,2:4));
hold off
clf
fill(x3',y3',field3')
axis(L)
kids = get(gca,'children');
for i = 1:length(kids)
%disp([int2str(i),'   ' ,get(kids(i),'type')])
   if strcmp(get(kids(i),'type'),'patch')
%disp(['visible setter'])
     set(kids(i),'edgecolor','none')
   end
end
%get(kids(100))
%figure(gcf)


%caxis([minfield maxfield]);

%axis('equal')
%axis(L)

nc = length(fieldcontours);
xc = NaN*ones(nc*length(x3),3);
yc = NaN*ones(nc*length(x3),3);
k = 0;

for inod = 1:length(x3)
c1 = fieldnod(meshele(inod,2));
c2 = fieldnod(meshele(inod,3));
c3 = fieldnod(meshele(inod,4));
x1 = x(meshele(inod,2));
x2 = x(meshele(inod,3));
x3 = x(meshele(inod,4));
y1 = y(meshele(inod,2));
y2 = y(meshele(inod,3));
y3 = y(meshele(inod,4));

for ic = 1:nc
c = fieldcontours(ic);
if ( ( c > c1 & c<c2 ) | ( c < c1 & c>c2 ) )
% contour passes between c1 .. c2
if ( ( c > c2 & c<c3 ) | ( c < c2 & c>c3 ) )
% contour passes between c1 .. c2  and between c2..c3
x12 = x1 + ( x2-x1 )*(c-c1)/(c2-c1);
y12 = y1 + ( y2-y1 )*(c-c1)/(c2-c1);
x23 = x2 + ( x3-x2 )*(c-c2)/(c3-c2);
y23 = y2 + ( y3-y2 )*(c-c2)/(c3-c2);
k = k+1;
xc(k,1:2) = [ x12 x23];
yc(k,1:2) = [ y12 y23];
%disp([' y12,y23 ', num2str(c1),' ',num2str(c),' ',num2str(c2),' '])
%disp([' ....... ', num2str(c2),' ',num2str(c),' ',num2str(c3),' '])
else
% contour passes between c1 .. c2  and between c1..c3
x12 = x1 + ( x2-x1 )*(c-c1)/(c2-c1);
y12 = y1 + ( y2-y1 )*(c-c1)/(c2-c1);
x13 = x1 + ( x3-x1 )*(c-c1)/(c3-c1);
y13 = y1 + ( y3-y1 )*(c-c1)/(c3-c1);
k = k+1;
xc(k,1:2) = [ x12 x13];
yc(k,1:2) = [ y12 y13];
%disp([' y12,y13 ', num2str(c1),' ',num2str(c),' ',num2str(c2),' '])
%disp([' ....... ', num2str(c1),' ',num2str(c),' ',num2str(c3),' '])
end

end % close if c1...c2
if ( ( c > c2 & c<c3 ) | ( c < c2 & c>c3 ) ) & ...
   ( ( c > c1 & c<c3 ) | ( c < c1 & c>c3 ) )
% contour passes between c2 .. c3  and therefore must between c1..c3
x13 = x1 + ( x3-x1 )*(c-c1)/(c3-c1);
y13 = y1 + ( y3-y1 )*(c-c1)/(c3-c1);
x23 = x2 + ( x3-x2 )*(c-c2)/(c3-c2);
y23 = y2 + ( y3-y2 )*(c-c2)/(c3-c2);
k = k+1;
xc(k,1:2) = [ x13 x23];
yc(k,1:2) = [ y13 y23];
%disp([' y23,y13 ', num2str(c1),' ',num2str(c),' ',num2str(c3),' '])
%disp([' ....... ', num2str(c2),' ',num2str(c),' ',num2str(c3),' '])
end

end % close loop on contour values
end % close loop on elements


xc = xc(1:k,:);
yc = yc(1:k,:);

hold on
plot(xc',yc','k')
hold off
