%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-13 16:57:28
 % @LastEditTime: 2022-01-13 17:12:54
%--------------------------
function [rho, cpnts] = init_trimesh_density(opt_load_rho, optDesign, vF, v0, nelx, nely, microx, poly)

    if opt_load_rho==0
        e1 = 1; nu1 = 0.3;
    
        triMesh = mexDT(poly, 10);
        nodes = triMesh{1};
        faces = triMesh{2};

        globalx = nelx * microx;
        globaly = nely * microx;
        maxx = globalx / 2; minx = -maxx;
        maxy = globaly / 2; miny = -maxy;          
        
        switch optDesign
        case 'cantilever'
            [~,fixeddofs,fext,~] = designDomain(optDesign,globalx,globaly, vF, microx, 1);

%             rho = top88(globalx,globaly,v0, e1, nu1, fext, fixeddofs);
            load("data/mat/rho_"+optDesign+".mat");
%             rho = flipud(rho);

            [cx, cy] = meshgrid(minx:1:maxx-1, miny:1:maxy-1);
            cpnts = [cx(:), cy(:)] + 0.5;

        case 'Lshape'
            [~, fixeddofs, fext, ~] = design_domain_nonregular(nodes, ...
                optDesign, vF, minx, maxx, miny, maxy);

            rho = top88_nonregular(v0, e1, nu1, fext, fixeddofs, nodes, faces);            
            cpnts = mapping_rho(rho, nodes, faces);
        end        

        save("data/mat/rho_"+optDesign+".mat", "rho", "cpnts");
    else
        load("data/mat/rho_"+optDesign+".mat");
    end
    
end


function cpnts = mapping_rho(rho, nodes, faces)

    x = nodes(:,1);
    y = nodes(:,2);

    cx = mean(x(faces), 2); % [m,1]
    cy = mean(y(faces), 2); % [m,1]

    cpnts = [cx, cy];
end