% Assembles self weight load for beam elements
function [F,dFdh] = assmFbeam(F,nodes,edges,b,h,massden)
DF = 0*F;
for i = 1:size(edges,1)
    n1 = edges(i,1);
    n2 = edges(i,2);
    x1 = nodes(n1,1); y1 = nodes(n1,2);
    x2 = nodes(n2,1); y2 = nodes(n2,2);
    L = ((x2-x1)^2+(y2-y1)^2)^0.5;
    b_beam = b(i);
    % Beam self weight
    fsw = L*b_beam*h*massden;
    loaddof = [6*n1-3,6*n2-3];
    DF(loaddof,1) = DF(loaddof,1) - fsw/2;   
end
F = F + DF;
dFdh = DF/h;
end