%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-03-18 21:53:54
 % @LastEditTime: 2021-05-31 16:54:40
%--------------------------
% get the node_id & ele_id to be removed whose rho = 0
function [ref_nodes_new, def_nodes_new, eles] = separate_rho(rho, edofMat, ref_nodes, def_nodes, density)

%     node_id_saved = [];
%     ele_id_saved  = [];   

    ref_nodes(:, 4) = 1:size(ref_nodes, 1);

%     for ele = 1 : size(edofMat, 1)
%         if rho(ele) == density
%             edof = edofMat(ele, :);
%             nid  = edof(3:3:end) / 3;
%             node_id_saved = [node_id_saved, nid];
%             ele_id_saved  = [ele_id_saved, ele];
%         end
%     end

    edofMat_new = edofMat(rho==density, :);
    node_id_saved = edofMat_new(:, 3:3:end) / 3;
    
    ref_nodes_new = ref_nodes(unique(node_id_saved), :);
    def_nodes_new = def_nodes(unique(node_id_saved), :);
%     edofMat_new   = edofMat(ele_id_saved, :);    

    eles = zeros(size(edofMat_new, 1), 8);

    % now node_id in edofMat_new is still old
    for ele = 1:size(edofMat_new, 1)
        edof = edofMat_new(ele, :);
        nid0 = edof(3:3:end) / 3;

        [nid, ~] = find(ref_nodes_new(:,4)==nid0);
        eles(ele, :) = nid;
    end

    ref_nodes_new = ref_nodes_new(:, 1:3); % remove last column nid
end