function [element_nodes_id, element_nodes_x, element_nodes_y] = ...
    find_bounding_nodes(dflow_floc, dflow_beg, dflow_end, poly_elems)
%Finds the bounding DFlow nodes for elements defined in input variable
%   Detailed explanation goes here

% Determine size of input variable and create node matrices 
% *This assumes that 
[r,~] = size(poly_elems);
element_nodes_id = zeros(r,7);
element_nodes_x = zeros(r,8);
element_nodes_y = zeros(r,8);

for j = 1:r
    
    j_elem = poly_elems(j,1);
    j_domain = poly_elems(j,2);
    j_domain_format = sprintf('%04d',int8(j_domain));
    j_file = strcat(dflow_floc,dflow_beg,j_domain_format,dflow_end);
    
    j_elem_ids = ncread(j_file,'FlowElemGlobalNr');
    j_elem_index = find(j_elem_ids(:) == j_elem);
    j_elem_nodes = ncread(j_file,'NetElemNode');
    j_elem_nodes_x = ncread(j_file,'FlowElemContour_x');
    j_elem_nodes_y = ncread(j_file,'FlowElemContour_y');
    
    element_nodes_id(j,1) = j_elem;
    element_nodes_id(j,2) = j_domain;
    element_nodes_x(j,1) = j_elem;
    element_nodes_x(j,2) = j_domain;
    element_nodes_y(j,1) = j_elem;
    element_nodes_y(j,2) = j_domain;
    
    [j_r,~] = size(j_elem_nodes);
    sides = 0;
    
    for k = 1:j_r
        element_nodes_id(j,k+2) = j_elem_nodes(k,j_elem_index);
        
        if element_nodes_id(j,k+2) > 0
            sides = sides + 1;
            element_nodes_x(j,k+3) = j_elem_nodes_x(k,j_elem_index);
            element_nodes_y(j,k+3) = j_elem_nodes_y(k,j_elem_index);
        else
            element_nodes_x(j,k+3) = 0;
            element_nodes_y(j,k+3) = 0;
        end
        
    end
    
    element_nodes_x(j,3) = sides;
    element_nodes_y(j,3) = sides;
    
end

end

