% Constraints and sensitivities
    constraints = mean(densities)-final_volume_fraction;
    constraints_sensitivities = fea.H*ones(size(fea.mesh.incidence,1),1)/size(fea.mesh.incidence,1);