fprintf('输入有误，请重新输入结点个数ne(ne为不小于1的正整数) \n ');

    ne=input('');

end

K1=zeros(2*nn);                   % K1为结构整体刚度矩阵

for i=1:ne

    fprintf('请输入第%d个单元的结点坐标co(形式为[x1 y1;x2 y2;……;x8 y8]);\n',i);

    co=input('');

    fprintf('请输入该单元的结点编号no(形式为[1 2……7 8]) \n ');

    no=input('');

    k1=k(co,E,NU,D,hou);          % 调用单元刚度矩阵函数

    K1=kk(no,k1,K1);              % 将单元刚度矩阵装配到整体刚度矩阵中

end

fprintf('未修正位移边界条件的整体刚度矩阵为K1=\n');

disp(K1);

fprintf('请输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn],约束端的外荷载按0输入) \n ');

F=input('');

while size(F)~=2*nn

    fprintf('输入有误。请重新输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn]) \n ');

    F=input('');

end

                                  %  以下采用对角元素乘大数法修正位移边界条件

fprintf('请输入x方向受到位移约束的结点编号及其约束位移值(形式为[node1 u1;node2 u2;……;noden un]) \n ');