shell单元的matlab代码

function [ek]=shellek（dyhmjdzbjdzb1dybh）
%----------------------------------------------------------------------------------------------   
% dyhm: node connectivity for each element  
% jdzb: coordinate of the nodes in the top surface of element
% jdzb1: coordinate of the nodes in the bottom surface of element
% dybh: the number of element whose stiffness matrix is calculated in this time of invoking. 
temp=0.577350292; % value needed by Gauss integral
gaosi=[ 1  1  1;
        1 -1  1;
       -1  1  1;
       -1 -1  1;
        1  1 -1;
        1 -1 -1;
       -1  1 -1;
       -1 -1 -1];
d=[]; % elastic matrix in 7.7-25. five order square matrix. 
% Its value is determined by the character of material and isn’t given in this example.

ek=zeros（40）;
    for i=1:8
        for j=1:3
            v3i（ij）=jdzb1（dybh（dyhmi）j）-jdzb（dybh（dyhmi）j）;
        end
    end
    tempi=[1 0 0];
    for i=1:8
       temp1=cross（v3i（i:）‘tempi）;
       xv2i（i:）=temp1/norm（temp1）;
       temp1=cross（xv2i（i:）v3i（i:））;
       xv1i（i:）=temp1/norm（temp1）;
    end
    xv3i=v3i/t;
    
    for gaosii=1:8  % begin of Gauss integral
    zb1=temp*gaosi（gaosii1）;
    zb2=temp*gaosi（gaosii2）;
    zb3=temp*gaosi（gaosii3）;

    ni（1）=1/4*（1+zb1）*（1+zb2）*（zb1+zb2-1）;
    ni（2）=1/4*（1-zb1）*（1+zb2）*（-zb1+zb2-1）;
    ni（3）=1/4*（1-zb1）*（1-zb2）*（-zb1-zb2-1）;
    ni（4）=1/4*（1+zb1）*（1-zb2）*（zb1-zb2-1）;
    ni（5）=1/2*（1-zb1^2）*（1+zb2）;
    ni（6）=1/2*（1-zb1）*（1-zb2^2）;
    ni（7）=1/2*（1-zb1^2）*（1-zb2）;
    ni（8）=1/2*（1+zb1）*（1-zb2^2）;
             v3=ni（1）*v3i（1:）‘+ni（2）*v3i（2:）‘+ni（3）*v3i（3:）‘+ni（4）*v3i（4:）‘+ni（5）*v3i（5:）‘+ni（6）*v3i（6:）‘+ni（7）*v3i（7:）‘+ni（8）*v3i（8:）‘;
    temp1=cross（v3tempi‘）;
    tn2=sqrt（temp1（1）^2+temp1（2）^2+temp1（3）^2）;
    xv2=temp1/tn2;
    temp1=cross（xv2v3）;
    tn1=sqrt（temp1（1）^2+temp1（2）^2+temp1（3）^2）;
    xv1=temp1/tn1;
    
    tn3=sqrt（v3（1）^2+v3（2）^2+v3（3）^2）;
    xv3=v3/tn3;
    theta=[xv1xv2xv3];
    
    for i=1:8
       for j=1:3
          zmtemp（ij）=（jdzb1（dybh（dyhmi）j）+jdzb（dybh（dyhmi）j））/2;
       end
    end
    
    
    pni（11）=1/4*（1+zb2）*（zb1+zb2-1）+（1/4+1/4*zb1）*（1+zb2）;
    pni（12）=（1/4+1/4*zb1）*（zb1+zb2-1）+（1/4+1/4*zb1）*（1+zb2）;
    pni（13）=0;
    pni（21）=-1/4*（1+zb2）*（-zb1+zb2-1）-（1/4-1/4*zb1）*（1+zb2）;
    pni（22）=（1/4-1/4*zb1）*（-zb1+zb2-1）+（1/4-1/4*zb1）*（1+zb2）;
    pni（23）=0;
    pni（31）=-1/4*（1-zb2）*（-zb1-zb2-1）-（1/4-1/4*zb1）*（1-zb2）;
    pni（32）=-（1/4-1/4*zb1）*（-zb1-zb2-1）-（1/4-1/4*zb1）*（1-zb2）;
    pni（33）=0;
    pni（41）=1/4*（1-zb2）*（zb1-zb2-1）+（1/4+1/4*zb1）*（1-zb2）;
    pni（42）=-（1/4+1/4*zb1）*（zb1-zb2-1）-（1/4+1/4*zb1）*（1-zb2）;
    pni（43）=0;
    pni（51）=-zb1*（1+zb2）;
    pni（52）=1/2-1/2*zb1^2;
    pni（53）=0;
    pni（61）=-1/2+1/2*zb2^2;
    pni（62）=-2*（1/2-1/2*zb1）*zb2;
    pni（63）=0;
    pni（71）=-zb1*（1-zb2）;
    pni（72）=-1/2+1/2*zb1^2;
    pni（73）=0;
    pni（81）=1/2-1/2*zb2^2;
    pni（82）=-2*（