tree.m

来自「四面体有限元程序,单元刚度的计算及总刚度矩征集成」· M 代码 · 共 25 行

theta=pi/6;
u =[0,0;0,1];  
subplot(2,3,1)
plot(u(:,1),u(:,2))
axis([-0.5,0.5,0,1])
for n=1:5
  uuu=[];
 for I=0:(length(u)/2-1)
       p1=(u(2*I+1,:)*2+u(2*I+2,:))/3 ;  
       p2=(u(2*I+1,:)+u(2*I+2,:)*2)/3;
       lp=[cos(theta),-sin(theta);sin(theta),cos(theta)]*...
       [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3 ;
       lp=p1+lp';
       rp=[cos(theta),sin(theta);-sin(theta),cos(theta)]*...
       [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3;
       rp=p2+rp';
       uu=[u(2*I+1,:);p1;p1;lp;p1;p2;p2;rp;p2;u(2*I+2,:)];
       uuu=[uuu;uu];
    end   
    u=[uuu];  
    subplot(2,3,n+1)
    plot(u(:,1),u(:,2))
    axis([-0.5,0.5,0,1])
  end