function NQM(fid_out)
%  ------  計算單元桿端力
global  element  delta  nfe  nfe_number nf  element_number %（定義全局變量）
for ie=1:1:element_number
    t=trans(ie);
    k=StiffnessMatrix(ie);
    edelta=elementdelta(ie,delta,element);
    locedelta=transpose(t)*edelta;
    N=k*locedelta;
    for j=1:1:nfe_number
        ee=nfe(j,1);
        if ee==ie 
            FF=efix(j);
            N=N+transpose(FF);
        end
    end
    fprintf( fid_out, 'element = %5d, Ni = %.8e, Qi = %.8e, Mi = %.8e ,Nj = %.8e, Qj = %.8e, Mj = %.8e\n', ...
   ie, N(1), N(2) , N(3) ,N(4), N(5) , N(6));  % output the node displacement
end