function [xyz,trilist,isbdry] = unitsquare_meshgen(n, type)
%% [xyz,trilist,isbdry] = unitsquare_meshgen(n, type)
%% This function generates a triangular mesh for the unit square
%% with boundary conditions described below.
%% Input arguments: 
%%   n - the number of nodes used will be n^2 arranged even an evenly
%%     spaced grid pattern.
%%   type - type of boundary conditions.
%%             0 = constant 1
%%             1 = linear: z=x+y
%%             2 = sin functions on each edge
%% Output arguments:
%%   xyz - an n^2 by 3 matrix; each row is an (x,y,z) coordinate of
%%      of a node.  The x and y coordinates are in the unit square.
%%      The z-coordinate is zero except on the boundaries.  On the
%%      boundaries, the z coordinate is set according to input argument
%%      type.
%%   trilist - an m-by-3 matrix, where m is the number of triangles.
%%      Each row gives the three indices (integers in 1:n^2) of the
%%      three vertices of the triangle.
%%   isbdry -- an n^2 column vector of type logical (0's and 1's).
%%      0 means the node is not on the boundary; 1 means that it is
%%      on the boundary.

xyz = zeros(n^2,3);
trilist = zeros(2*(n-1)^2,3);
isbdry = logical(zeros(n^2,1));
pos = 1;
for j = 1 : n
  for k = 1 : n
    xyz(pos,:) = [(j - 1) / (n - 1), (k - 1) / (n - 1), 0];
    if j == 1 | j == n | k == 1 | k == n
      isbdry(pos) = 1;
      if type == 0
        xyz(pos,3) = 1;
      elseif type == 1
        xyz(pos,3) = xyz(pos,1) + xyz(pos,2);
      elseif type == 2
        if j == 1 | j == n
          xyz(pos,3) = sin((k - 1) / (n - 1) * pi);
        else
          xyz(pos,3) = -sin((j - 1) / (n - 1) * pi);
        end
      end
    end
    pos = pos + 1;
  end
end

pos = 1;
for j = 1 : n - 1
  for k = 1 : n - 1
    trilist(pos,:) = [(j - 1) * n + k, j * n + k, j * n + k + 1];
    pos = pos + 1;
    trilist(pos,:) = [(j - 1) * n + k, j * n + k + 1, (j - 1) * n + k + 1];
    pos = pos + 1;
  end
end

