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*（