function pml_mesh = project_pml_xyz(bound, L_pml, N_pml, ctrl, tol)
%PROJECT_ Summary of this function goes here
%   Detailed explanation goes here

% CTRL: a - both plus and minus
%       m - minus only
%       p - plus only
%       n - do nothing

if nargin < 5
    tol = 1e-4;
end

if nargin < 4 || isempty(ctrl)
    ctrl = 'aaa';
end

N_dims = 3;

bb_orig = boundingbox(bound);

pml_mesh = create_empty_mesh;


% Go through each dimension
for i_dim = 1 : N_dims
    ctr = ctrl(i_dim);
    if ctr == 'n';
        continue
    end
    % Get the bounding box
    bb = boundingbox(bound);
    
    for d = 1 : 2 % d:1 negative, d:2 positive
        if (d == 2 && ctr == 'p') || (d == 2 && ctr == 'n') || ctr == 'a'
            % Create a box for sectioning the boundary
            box = circshift([bb(d, i_dim) -inf -inf; bb(d, i_dim) inf inf], [0 i_dim-1]);
            box = bsxfun(@plus, box, [-tol; tol]);
            bou = drop_mesh_IDs(drop_unused_nodes(mesh_section(bound, box)));

            % Extrude the boundary
            dir = circshift([(-1)^d 0 0]*L_pml/N_pml, [0 i_dim-1]);
            ext = extrude_mesh(bou, dir, N_pml);
            
            % Join extrusion with existing ones
            pml_mesh = join_meshes(pml_mesh, ext);

            % Update the boundary
            bound = merge_coincident_nodes( ...
                join_meshes(bound, get_boundary(ext)), tol); 
        end
    end
end

% Assign PML element types to extrusion
pml_mesh = drop_mesh_IDs(merge_coincident_nodes(pml_mesh, tol));
pml_mesh.Elements(:, 2) = pml_mesh.Elements(:,2) + 200;

pml_mesh.Properties = [1 3 1 L_pml 0 0];
pml_mesh.PMLData = {1 @pml_sigma_glob @pml_sigma_glob @pml_sigma_glob ...
        bb_orig(:,1) bb_orig(:,2) bb_orig(:,3)};

end
