function CreateMeshTriangle(name,xy,area,options)
#[...] = CreateMeshTriangle(...)
# generate a mesh using triangle
#
# CreateMeshTriangle(name,xy,area,options)
#  name the base filename: the file name.poly will be generated
#       then triangle will generate files name.1.* with the mesh
#
#  xy   vector containing the coordinates of the nodes forming the 
#       outer boundary. Currently no holes can be generated. The format is
#       [x1,y1,t1;x2,y2,t2;...;xn,yn,tn] where
#       xi x-coordinate of node i
#       yi x-coordinate of node i
#       ti boundary marker for segment from node i to node i+1
#         bi=1  Dirichlet condition
#         bi=2  Neumann condition
#       the last given node will be connected to the first given node 
#       to create a closed curve
#
#  area  a scalar given the maximal area of the triangles to be generated
#
#  options additional options to be used when calling triangle
#      the options "pa" and the area will be added automatically
#      default options are "q" resp. "qpa"
#      to supress the verbose information use "Q"
#  The information can then be read and used by 
#     [nodes,elem,edges]=ReadMeshTriangle("name.1");

    if ((nargin<3)|(nargin>4)) 
  usage("CreateMeshTriangle(name,xy,area,options)");
endif

opt="q";
if (nargin==4) opt=options; endif

n=length(xy);

[fid,msg]=fopen(strcat(name,".poly"),"w");
fprintf(fid,"# nodes\n");
fprintf(fid,"%i 2 0 1\n",n);
fprintf(fid,"%i %g %g %d\n",1,xy(1,1),xy(1,2),min(xy(1,3),xy(n,3)));
for k=2:n
  fprintf(fid,"%i %g %g %d\n",k,xy(k,1),xy(k,2),min(xy(k,3),xy(k-1,3)));
endfor

fprintf(fid,"# segments\n");
fprintf(fid,"%i 1\n",n);
for k=1:n-1
  fprintf(fid,"%i %i %i %i\n",k,k,k+1,xy(k,3));
endfor
fprintf(fid,"%i %i %i %i\n\n",n,n,1,xy(n,3));

fprintf(fid,"# holes\n");
fprintf(fid,"0\n",n);

fclose(fid);

system(strcat("triangle -",opt,sprintf("pa%f  ",area),name,".poly"));
system(strcat("CuthillMcKee -s0 -g ",name,".1"));

endfunction




