function mesh = mesh_graded_bricks(x, y, z, p, h, g)
%MESH_GRADED_BRICKS Summary of this function goes here
%   Detailed explanation goes here

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

% Edges along xyz with brick index
e_xyz = [kron((1 : size(h,1)).', ones(3,1)), p(e(:,1), :), p(e(:,2), :)];

e_x = e_xyz(e_xyz(:,2) ~= e_xyz(:,5), [1 2 5]);
e_y = e_xyz(e_xyz(:,3) ~= e_xyz(:,6), [1 3 6]);
e_z = e_xyz(e_xyz(:,4) ~= e_xyz(:,7), [1 4 7]);

[e_ux, ~, ie_ux] = unique(e_x(:,2:3), 'rows');
[e_uy, ~, ie_uy] = unique(e_y(:,2:3), 'rows');
[e_uz, ~, ie_uz] = unique(e_z(:,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

lines_z = cell(size(e_uz,1), 1);
for i_z = 1 : size(e_uz, 1)
    lines_z{i_z} = mesh_graded_line( ...
        z(e_uz(i_z,:), 1), z(e_uz(i_z,:), 2), g(3));
end

% Mesh the hexes
bricks = cell(size(h, 1), 1);
for i = 1 : size(h, 1)
    bricks{i} = create_brick( ...
        lines_x{ie_ux(i)}, lines_y{ie_uy(i)}, lines_z{ie_uz(i)});
end

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

end

