% Sensitivity of load-dependent compliance w.r.t nodes: beams
function dFdxU = assmdFdXU_beams(U,nodes,edges,b,h,massden)
dFdxU = zeros(size(nodes,1)*2,1);
% Loop over beams
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);
    loaddof = [6*n1-3,6*n2-3];
    b_beam = b(i);
    coeff = 0.5*((x2-x1)^2+(y2-y1)^2)^(-0.5);
    % Derivatives of beam length w.r.t. coors
    dLdx1 = -coeff*2*(x2-x1);
    dLdx2 = coeff*2*(x2-x1);
    dLdy1 = -coeff*2*(y2-y1);
    dLdy2 = coeff*2*(y2-y1);
    % Multiply load by U
    dFU = -b_beam*h*massden/2*sum(U(loaddof,1));
    % Assemble
    dFdxU(2*n1-1) = dFdxU(2*n1-1) + dFU*dLdx1;
    dFdxU(2*n1) = dFdxU(2*n1) + dFU*dLdy1;
    dFdxU(2*n2-1) = dFdxU(2*n2-1) + dFU*dLdx2;
    dFdxU(2*n2) = dFdxU(2*n2) + dFU*dLdy2;
end

