snow.m

四面体有限元程序,单元刚度的计算及总刚度矩征集成

clear all
new=[0.5+(sqrt(3)/2*i),-0.5+(sqrt(3)/2*i),0,0.5+(sqrt(3)/2*i)];
subplot(2,3,1);
plot(new)
axis equal
for kkk=1:5;
   old=new;
   n=length(old)-1;
   for jj=0:n-1;
      diff=(old(jj+2)-old(jj+1))/3;
      new(4*jj+1)=old(jj+1);
      new(4*jj+2)=old(jj+1)+diff;
      new(4*jj+3)=new(4*jj+2)+diff*((1-sqrt(3)*i)/2);
      new(4*jj+4)=old(jj+1)+2*diff;
   end
   new(4*n+1)=old(n+1);
   subplot(2,3,kkk+1)
   plot(new)
   axis equal;
end