function mesh = mesh_graded_quads(x, y, p, q, g)
%MESH_GRADED_BRICKS Summary of this function goes here
%   Detailed explanation goes here

% Extract different edges
e = reshape(q(:, [1 2 1 4]).', 2, []).';

% Edges along xy with brick index
e_xy = [kron((1 : size(q,1)).', ones(2,1)), p(e(:,1), :), p(e(:,2), :)];

e_x = e_xy(e_xy(:,2) ~= e_xy(:,4), [1 2 4]);
e_y = e_xy(e_xy(:,3) ~= e_xy(:,5), [1 3 5]);

[e_ux, ~, ie_ux] = unique(e_x(:,2:3), 'rows');
[e_uy, ~, ie_uy] = unique(e_y(:,2:3), 'rows');

% Mesh the unique lines using the specified grading
lines_x = cell(size(e_ux,1), 1);
for i_x = 1 : size(e_ux, 1)
    lines_x{i_x} = mesh_graded_line( ...
        x(e_ux(i_x,:), 1), x(e_ux(i_x,:), 2), g(1));
end

lines_y = cell(size(e_uy,1), 1);
for i_y = 1 : size(e_uy, 1)
    lines_y{i_y} = mesh_graded_line( ...
        y(e_uy(i_y,:), 1), y(e_uy(i_y,:), 2), g(2));
end

% Mesh the hexes
bricks = cell(size(q, 1), 1);
for i = 1 : size(q, 1)
    bricks{i} = create_slab(lines_x{ie_ux(i)}(:), lines_y{ie_uy(i)}(:));
end

% Join the whole mesh
if size(q,1) > 1
    mesh = drop_mesh_IDs(drop_unused_nodes(merge_coincident_nodes( ...
        join_meshes(bricks{:}), min([x(:,2); y(:,2)])/1e6)));
else
    mesh = bricks{1};
end

end

