% Creates a unifom mesh of ex by ey elements 
% with length of lx and ly in the x and y directions.
% The origin (lower left) can be optionally specified as x0, y0.
function [m] = uniform_mesh(ex, ey, lx, ly, x0, y0)
    % if origin is not specified, then set it to zero.
    if nargin < 5, x0 = 0; end
    if nargin < 6, y0 = 0; end       
        
    m.num_nodes = (ex+1)*(ey+1);
    % Nodal reference coordinates.
    m.X    = zeros(2, m.num_nodes);   
    for j=1:ey+1
        for i=1:ex+1
            m.X(:,i+(j-1)*(ex+1)) = [(i-1)*lx/ex + x0; (j-1)*ly/ey + y0];
        end
    end
    m.num_elements = ex*ey;
    m.connectivity = zeros(4, m.num_elements);
    for j=1:ey
        for i=1:ex            
            % first node in element
            n0 = i+(j-1)*(ex+1);            
            m.connectivity(:,i+(j-1)*ex) = [n0; n0+1; n0+ex+2; n0+1+ex];            
        end
    end
        
    % 4 point Gaussian quadrature rule.
    m.quad_points = [-1, 1, 1,-1;
                     -1,-1, 1, 1] / sqrt(3);
    m.quad_points = [m.quad_points; 1,1,1,1];

    % Call as mesh.draw() to draw this mesh.
    m.draw =       @() plot_mesh(m);
    % Call as mesh.draw_nodal(f) to plot mesh colored by nodal value f.
    m.draw_nodal = @(f) plot_nodal(m,f);
end

function [] = plot_mesh(m)
    p.vertices = m.X';
    p.faces = m.connectivity';
    p.facecolor = 'none';
    patch(p);
end

function[] = plot_nodal(m, f)
    p.vertices = m.X';
    p.faces = m.connectivity';
    p.facecolor = 'interp';
    p.facevertexcdata = f;
    patch(p);    
end

