Scilab function

mesh2d - triangulation of n points in the plane

### Calling Sequence

  `[nutr,A] = mesh2d(x,y,[front])`

### Parameters

- `x` : real row array
- `y` : real row array
- `front` : integer row array
- `nutr` : integer matrix
- `A` : sparse 0-1 matrix

### Description

The arrays `x` and `y` are the coordinates of n points in the plane. `mesh2d` returns a matrix `nutr(3,nbt)` of the numbers of the nodes of the `nbt` triangles of the triangulation of the points. It returns also a sparse matrix `A` representing the connections between the nodes (`A(i,j)=1` if `(i,j)` is a side of one of the triangles or `i=j` ). In the case of 3 parameters `front` is the array defining the boundary: it is the array of the indices of the points located on the boundary . The boundary is defined such that the normal to the boundary is oriented towards outside. The boundary is given by its connected components: a component is the part `(i1,i2)` such that `front(i1)=front(i2)` (the external boundary is defined in the counterclockwise way, see the examples below). The error cases are the following: err = 0 if no errors were encountered; err = 3 all nodes are collinear.

If the boundary is given, the other error cases are: err = 2 some points are identical; err = 5 wrong boundary array; err = 6 crossed boundary; err = 7 wrong orientation of the boundary; err = 10 an interior point is on the boundary; err = 8 size limitation; err = 9 crossed boundary; err = 12 some points are identical or size limitation.

### Examples

    // FIRST CASE
    theta=0.025*[1:40]*2.*%pi;
    x=1+cos(theta);
    y=1.+sin(theta);
    theta=0.05*[1:20]*2.*%pi;
    x1=1.3+0.4*cos(theta);
    y1=1.+0.4*sin(theta);
    theta=0.1*[1:10]*2.*%pi;
    x2=0.5+0.2*cos(theta);
    y2=1.+0.2*sin(theta);
    x=[x x1 x2];
    y=[y y1 y2];
    //
    nu=mesh2d(x,y);
    nbt=size(nu,2);
    jj=[nu(1,:)' nu(2,:)';nu(2,:)' nu(3,:)';nu(3,:)' nu(1,:)'];
    as=sparse(jj,ones(size(jj,1),1));
    ast=tril(as+abs(as'-as));
    [jj,v,mn]=spget(ast);
    n=size(x,2);
    g=make_graph('foo',0,n,jj(:,1)',jj(:,2)');
    g('node_x')=300*x;
    g('node_y')=300*y;
    g('default_node_diam')=10;
    show_graph(g)
    // SECOND CASE !!! NEEDS x,y FROM FIRST CASE
    x3=2.*rand(1:200);
    y3=2.*rand(1:200);
    wai=((x3-1).*(x3-1)+(y3-1).*(y3-1));
    ii=find(wai >= .94);
    x3(ii)=[];y3(ii)=[];
    wai=((x3-0.5).*(x3-0.5)+(y3-1).*(y3-1));
    ii=find(wai <= 0.055);
    x3(ii)=[];y3(ii)=[];
    wai=((x3-1.3).*(x3-1.3)+(y3-1).*(y3-1));
    ii=find(wai <= 0.21);
    x3(ii)=[];y3(ii)=[];
    xnew=[x x3];ynew=[y y3];
    fr1=[[1:40] 1];fr2=[[41:60] 41];fr2=fr2($:-1:1);
    fr3=[[61:70] 61];fr3=fr3($:-1:1);
    front=[fr1 fr2 fr3];
    //
    nu=mesh2d(xnew,ynew,front);
    nbt=size(nu,2);
    jj=[nu(1,:)' nu(2,:)';nu(2,:)' nu(3,:)';nu(3,:)' nu(1,:)'];
    as=sparse(jj,ones(size(jj,1),1));
    ast=tril(as+abs(as'-as));
    [jj,v,mn]=spget(ast);
    n=size(xnew,2);
    g=make_graph('foo',0,n,jj(:,1)',jj(:,2)');
    g('node_x')=300*xnew;
    g('node_y')=300*ynew;
    g('default_node_diam')=10;
    show_graph(g)
    // REGULAR CASE !!! NEEDS PREVIOUS CASES FOR x,y,front
    xx=0.1*[1:20];
    yy=xx.*.ones(1,20);
    zz= ones(1,20).*.xx;
    x3=yy;y3=zz;
    wai=((x3-1).*(x3-1)+(y3-1).*(y3-1));
    ii=find(wai >= .94);
    x3(ii)=[];y3(ii)=[];
    wai=((x3-0.5).*(x3-0.5)+(y3-1).*(y3-1));
    ii=find(wai <= 0.055);
    x3(ii)=[];y3(ii)=[];
    wai=((x3-1.3).*(x3-1.3)+(y3-1).*(y3-1));
    ii=find(wai <= 0.21);
    x3(ii)=[];y3(ii)=[];
    xnew=[x x3];ynew=[y y3];
    nu=mesh2d(xnew,ynew,front);
    nbt=size(nu,2);
    jj=[nu(1,:)' nu(2,:)';nu(2,:)' nu(3,:)';nu(3,:)' nu(1,:)'];
    as=sparse(jj,ones(size(jj,1),1));
    ast=tril(as+abs(as'-as));
    [jj,v,mn]=spget(ast);
    n=size(xnew,2);
    g=make_graph('foo',0,n,jj(:,1)',jj(:,2)');
    g('node_x')=300*xnew;
    g('node_y')=300*ynew;
    g('default_node_diam')=3;
    show_graph(g)

    //An example with a random set of points
    function []=test(X,Y)
      Tr=mesh2d(X,Y);
      plot2d(X,Y,[-1,-2,3]);
      [m,n]=size(Tr);
      xpols= matrix(X(Tr),m,n);
      ypols= matrix(Y(Tr),m,n);
      xset("colormap",rand(2*n,3));
      xfpolys(xpols,ypols,[n/4:n/4+n-1]);
    endfunction
    N=1000;xbasc();X=rand(1,N); Y=rand(1,N);
    xset("wdim",700,700);
    test(X,Y);