xs=diff(x,s);

xt=diff(x,t);

ys=diff(y,s);

yt=diff(y,t);

J=xs*yt-ys*xt;

for i=1:8

    Ns(i)=diff(N(i),s);

    Nt(i)=diff(N(i),t);

end

for i=1:8

    g(i)=yt*Ns(i)-ys*Nt(i);

    h(i)=xs*Nt(i)-xt*Ns(i);

end

B=sym(zeros(3,16));

for i=1:8

    B(1,2*i-1)=B(1,2*i-1)+g(i);

    B(2,2*i)=B(2,2*i)+h(i);

    B(3,2*i-1)=B(3,2*i-1)+h(i);

    B(3,2*i)=B(3,2*i)+g(i);

end

Bnew=simplify(B);

Jnew=simplify(J);

BD=transpose(Bnew)*D*Bnew/Jnew;

r=int(int(BD,t,-1,1),s,-1,1);

z=hou*r;

k=double(z);

〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓

function K=kk(no,k,K)

%  no表示相应单元在整体坐标系中的编号，k为该单元的刚度矩阵

%  K返回结构的整体刚度矩阵

for i=1:8

    for j=1:8

        K(2*no(i)-1,2*no(j)-1)=K(2*no(i)-1,2*no(j)-1)+k(2*i-1,2*j-1);

        K(2*no(i)-1,2*no(j))=K(2*no(i)-1,2*no(j))+k(2*i-1,2*j);

        K(2*no(i),2*no(j)-1)=K(2*no(i),2*no(j)-1)+k(2*i,2*j-1);

        K(2*no(i),2*no(j))=K(2*no(i),2*no(j))+k(2*i,2*j);

    end

end