% Assembles self weight load for CST-DKT elements
function [F,dF] = assmFplate(F,nodes,plates,th,massden)
nplates = size(plates,1);
DF = 0*F;
for i = 1:nplates
    % Element nodes
    n1 = plates(i,1); n2 = plates(i,2); n3 = plates(i,3);
    % Node coordinates
    plate_X = [nodes(n1,1) nodes(n2,1) nodes(n3,1)];
    plate_Y = [nodes(n1,2) nodes(n2,2) nodes(n3,2)];
    % Area
    area = 0.5*(plate_X(1)*(plate_Y(2)-plate_Y(3)) + ...
         plate_X(2)*(plate_Y(3)-plate_Y(1)) + ...
         plate_X(3)*(plate_Y(1)-plate_Y(2)));
    % Total self weight
    fsw = area*th*massden;
    loaddofs = [6*n1-3,6*n2-3,6*n3-3];
    DF(loaddofs,1) = DF(loaddofs,1) - fsw/3;
end
F = F + DF;
dF = DF/th;
end