function ShowLevelCurves(nodes,elem,u,levels)
%ShowLevelCurves(...)
%  shows level curves of a function on the mesh
%
%ShowLevelCurves(nodes,elem,u,levels);
%nodes, elem contain information about the mesh
%         see ReadMesh() for the description of the format
%u     contains the values of the solution at the nodes
%level is a list of values for which the level curves are drawn
%
%see also ReadMesh, ShowMesh, FEMEquation, FEMSolve, FEMValue


levels=sort(levels);
lmin=levels(1);
lmax=levels(length(levels));

[fid,msg]=fopen("/tmp/leveldata.gnu","w");

for el=1:length(elem)
  values=[u(elem(el,1)),u(elem(el,2)),u(elem(el,3))];
  minval=min(values);
  maxval=max(values);
  locallevel=levels(minval<=levels+eps);      # select the levels to be plotted
  if size(locallevel)>0
    locallevel=locallevel(locallevel<=maxval+eps);
  endif
  for level=locallevel
    points=zeros(3,2);
    npoints=1;
    dl=values-level;
%    if (dl(1)==0) points(npoints,:)=nodes(elem(el,1),1:2);npoints++;endif 
%    if (dl(2)==0) points(npoints,:)=nodes(elem(el,2),1:2);npoints++;endif 
%    if (dl(3)==0) points(npoints,:)=nodes(elem(el,3),1:2);npoints++;endif 

    if (abs(dl(1))<=10*eps) 
      points(npoints,:)=nodes(elem(el,1),1:2);npoints++;endif 
    if (abs(dl(2))<=10*eps)  
      points(npoints,:)=nodes(elem(el,2),1:2);npoints++;endif 
    if (abs(dl(3))<=10*eps) 
      points(npoints,:)=nodes(elem(el,3),1:2);npoints++;endif 

    if (npoints<=2)
      if ((dl(1)*dl(2)) < 0)    # intersection between 1st and 2nd point
	points(npoints,:)= ( dl(2)*nodes(elem(el,1),1:2)...
			    -dl(1)*nodes(elem(el,2),1:2))/(dl(2)-dl(1));
	npoints++;
      endif
      if ((dl(1)*dl(3)) < 0)    # intersection between 1st and 3rd point
	points(npoints,:)= ( dl(3)*nodes(elem(el,1),1:2)...
			    -dl(1)*nodes(elem(el,3),1:2))/(dl(3)-dl(1));
	npoints++;
      endif
      if ((dl(3)*dl(2)) < 0)    # intersection between 2nd and 3rd point
	points(npoints,:)= ( dl(2)*nodes(elem(el,3),1:2)...
			    -dl(3)*nodes(elem(el,2),1:2))/(dl(2)-dl(3));
	npoints++;
      endif
    endif
    for k=1:npoints-1
      fprintf(fid,"%g %g\n",points(k,:));
    endfor
    if npoints>=4 # full triangle  write to file
      fprintf(fid,"%g %g\n",points(1,: )); 
    endif
    fprintf(fid,"\n");   
  endfor # level
endfor # el

fclose(fid);

gplot  "/tmp/leveldata.gnu" with lines 4
endfunction
