plot_modal_shape.m

结构力学中的有限元例子,包含了7个分类文件夹

function plot_modal_shape ( in_data, D, dof, NF) 

% plots modal shape for CST, CSQ models
%    use with modal analysis
maxX=0; minX=0; maxY=0; minY=0; maxZ=0; minZ=0; labz=0;
if in_data.EL(1,2)==9              % "9" - BCIZ element (bending triangle)
    dof_per_node = 3;
    maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
    maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
    labx = (maxX / 6); laby =  (maxY / 6);
    labx = min(labx,laby); laby = labx;
    SCz = labx*1.2; % Z deflection scale
    labz = SCz/2;

    max_D_z = max(abs(D(1:3:dof(1)))); 
    ND_d = in_data.ND;
    ND_d(:,4) = (D(1:3:dof(1))'./max_D_z)*0.5*SCz;
end;

if in_data.EL(1,2)==0 | in_data.EL(1,2)==1 | in_data.EL(1,2)==2 %  2D-beam element
    n_dof_node = 3;
    maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
    maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
    labx = (maxX / 6); laby = (maxY / 6);
    labx = min([labx laby]); laby = labx;
    max_D_x = max(abs(D(1:n_dof_node:dof(1))));
    max_D_y = max(abs(D(2:n_dof_node:dof(1))));
    ND_d = in_data.ND;
    deN = max([max_D_x  max_D_y]);
    ND_d(:,2) = in_data.ND(:,2)+(D(1:n_dof_node:dof(1))'./deN)*.5*labx;
    ND_d(:,3) = in_data.ND(:,3)+(D(2:n_dof_node:dof(1))'./deN)*.5*laby;
end;
if in_data.EL(1,2)==4 | in_data.EL(1,2)==5 % "4/5" - CST/CSQ element (plain stress) 
    n_dof_node = 2;
    maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
    maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
    labx = (maxX / 6); laby = (maxY / 6);
     labx = min([labx laby]); laby = labx;
    max_D_x = max(abs(D(1:n_dof_node:dof(1))));
    max_D_y = max(abs(D(n_dof_node:n_dof_node:dof(1))));
    ND_d = in_data.ND;
    deN = max([max_D_x  max_D_y]);
    ND_d(:,2) = in_data.ND(:,2)+(D(1:n_dof_node:dof(1))'./deN)*.5*labx;
    ND_d(:,3) = in_data.ND(:,3)+(D(n_dof_node:n_dof_node:dof(1))'./deN)*.5*laby;
end;
if in_data.EL(1,2)==3 % "3" - 3D-beam element
    n_dof_node = 6;
    maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
    maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
    maxZ = max(in_data.ND(:,4)); minZ = min(in_data.ND(:,4));
    labx = (maxX / 6); laby = (maxY / 6); labz = (maxZ / 6);
    labx = min([labx laby labz]); laby = labx;  labz=labx;
    max_D_x = max(abs(D(1:n_dof_node:dof(1))));
    max_D_y = max(abs(D(2:n_dof_node:dof(1))));
    max_D_z = max(abs(D(3:n_dof_node:dof(1))));
    ND_d = in_data.ND;
    deN = max([max_D_x  max_D_y  max_D_z]);
    ND_d(:,2) = in_data.ND(:,2)+(D(1:n_dof_node:dof(1))'./deN)*.5*labx;
    ND_d(:,3) = in_data.ND(:,3)+(D(2:n_dof_node:dof(1))'./deN)*.5*laby;
    ND_d(:,4) = in_data.ND(:,4)+(D(3:n_dof_node:dof(1))'./deN)*.5*labz;
end;
if in_data.EL(1,2)==6              % "6" - BRICK element (8-nodes)
    n_dof_node = 3;
    maxX = max(in_data.ND(:,2)); minX = min(in_data.ND(:,2));
    maxY = max(in_data.ND(:,3)); minY = min(in_data.ND(:,3));
    maxZ = max(in_data.ND(:,4)); minZ = min(in_data.ND(:,4));
    labx = (maxX / 6); laby = (maxY / 6); labz = (maxZ / 6);
    labx = min([labx laby labz]); laby = labx;  labz=labx;
    max_D_x = max(abs(D(1:n_dof_node:dof(1))));
    max_D_y = max(abs(D(2:n_dof_node:dof(1))));
    max_D_z = max(abs(D(3:n_dof_node:dof(1))));
    ND_d = in_data.ND;
    deN = max([max_D_x  max_D_y  max_D_z]);
    ND_d(:,2) = in_data.ND(:,2)+(D(1:n_dof_node:dof(1))'./deN)*0.5*labx;
    ND_d(:,3) = in_data.ND(:,3)+(D(2:n_dof_node:dof(1))'./deN)*0.5*laby;
    ND_d(:,4) = in_data.ND(:,4)+(D(3:n_dof_node:dof(1))'./deN)*0.5*labz;
end;

figure(NF);
if in_data.EL(1,2)==0 | in_data.EL(1,2)==1 | in_data.EL(1,2)==2 % "0/1/2" - 2D-beam element
    plot(ND_d(:,2),ND_d(:,3),'r.'); hold on;
end;

if in_data.EL(1,2)==4 | in_data.EL(1,2)==5 % "4/5" - CST/CSQ element 
    plot(ND_d(:,2),ND_d(:,3),'r.'); hold on;
end;
if in_data.EL(1,2)==3 | in_data.EL(1,2)==6 % "3/6" - 3D-beam/BRICK element
    plot3(ND_d(:,2),ND_d(:,3),ND_d(:,4),'r.'); hold on;
end;

for i=1:size(in_data.EL)
    if in_data.EL(i,2)==9              % "9" - BCIZ element (bending triangle)  
        % plot3(ND_d(:,2),ND_d(:,3),ND_d(:,4),'r.');
        node1 = find(ND_d(:,1)==in_data.EL(i,3));
        node2 = find(ND_d(:,1)==in_data.EL(i,4));
        node3 = find(ND_d(:,1)==in_data.EL(i,5));
        plot3([ND_d(node1,2) ND_d(node2,2) ND_d(node3,2) ND_d(node1,2)], ...
            [ND_d(node1,3) ND_d(node2,3) ND_d(node3,3) ND_d(node1,3)], ...
            [ND_d(node1,4) ND_d(node2,4) ND_d(node3,4) ND_d(node1,4)],'r-');
        axis equal; axis off; view(3); hold on;
        axis([(-1.1*minX) (1.1*maxX) (-1.1*minY) (1.1*maxY) (-SCz)    (SCz)]);
    end;


    if in_data.EL(i,2)==0 | in_data.EL(i,2)==1 | in_data.EL(i,2)==2  % "0/1/2" - 2D-beam element
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        plot([ND_d(node1,2) ND_d(node2,2)], ...
            [ND_d(node1,3) ND_d(node2,3)],'r-');
        axis equal; axis off;
        axis([(minX-labx) (maxX+labx) (minY-laby) (maxY+laby)]); % plot scale
    end;    
    if in_data.EL(i,2)==4  % "4" - CST element (triangle) 
        node1 = find(ND_d(:,1)==in_data.EL(i,3));
        node2 = find(ND_d(:,1)==in_data.EL(i,4));
        node3 = find(ND_d(:,1)==in_data.EL(i,5));
        plot([ND_d(node1,2) ND_d(node2,2) ND_d(node3,2) ND_d(node1,2)], ...
            [ND_d(node1,3) ND_d(node2,3) ND_d(node3,3) ND_d(node1,3)],'r-');
        axis equal; axis off;
        axis([(minX-labx) (maxX+labx) (minY-laby) (maxY+laby) ]); % plot scale
    end;
    if in_data.EL(i,2)==5  % "5" - CSQ element (quadrilateral) 
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
        plot([ND_d(node1,2) ND_d(node2,2) ND_d(node3,2) ND_d(node4,2) ND_d(node1,2)],...
            [ND_d(node1,3) ND_d(node2,3) ND_d(node3,3) ND_d(node4,3) ND_d(node1,3)],'r-');
        axis equal; axis off;
        axis([(minX-labx) (maxX+labx) (minY-laby) (maxY+laby)]); % plot scale
    end;
    if in_data.EL(i,2)==3  % "3" - 3D-beam element
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        plot3([ND_d(node1,2) ND_d(node2,2)], ...
            [ND_d(node1,3) ND_d(node2,3)],[ND_d(node1,4) ND_d(node2,4)],'r-');
        axis equal; axis off;
        axis([(minX-labx) (maxX+labx) (minY-laby) (maxY+laby) (minZ-labz)    (maxZ+labz)]); % plot scale
    end;
    if in_data.EL(1,2)==6              % "6" - BRICK element (8-nodes)
        node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
        node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
        node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
        node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
        node5 = find(in_data.ND(:,1)==in_data.EL(i,7));
        node6 = find(in_data.ND(:,1)==in_data.EL(i,8));
        node7 = find(in_data.ND(:,1)==in_data.EL(i,9));
        node8 = find(in_data.ND(:,1)==in_data.EL(i,10));
        
        plot3([ND_d(node1,2) ND_d(node2,2) ND_d(node3,2) ND_d(node4,2) ND_d(node1,2) ...
                ND_d(node5,2) ND_d(node6,2) ND_d(node7,2) ND_d(node8,2) ND_d(node5,2)],...
            [ND_d(node1,3) ND_d(node2,3) ND_d(node3,3) ND_d(node4,3) ND_d(node1,3) ...
                ND_d(node5,3) ND_d(node6,3) ND_d(node7,3) ND_d(node8,3) ND_d(node5,3)],...
            [ND_d(node1,4) ND_d(node2,4) ND_d(node3,4) ND_d(node4,4) ND_d(node1,4) ...
                ND_d(node5,4) ND_d(node6,4) ND_d(node7,4) ND_d(node8,4) ND_d(node5,4)],'r-','LineWidth',1);    
    
        plot3([ND_d(node2,2) ND_d(node6,2) ND_d(node7,2) ND_d(node3,2) ND_d(node2,2) ...
                ND_d(node1,2) ND_d(node5,2) ND_d(node8,2) ND_d(node4,2) ND_d(node1,2)],...
            [ND_d(node2,3) ND_d(node6,3) ND_d(node7,3) ND_d(node3,3) ND_d(node2,3) ...
                ND_d(node1,3) ND_d(node5,3) ND_d(node8,3) ND_d(node4,3) ND_d(node1,3)],...
            [ND_d(node2,4) ND_d(node6,4) ND_d(node7,4) ND_d(node3,4) ND_d(node2,4) ...
                ND_d(node1,4) ND_d(node5,4) ND_d(node8,4) ND_d(node4,4) ND_d(node1,4)],'r-','LineWidth',1);
        axis equal; axis off;
        axis([(minX-labx) (maxX+labx) (minY-laby) (maxY+laby) (minZ-labz)    (maxZ+labz)]); % plot scale        
    end;
end;
% title(['Mode shape:  '  num2str(in_data.dynam.MODA)  ]); hold on;

% plot constrains ==============================================================
if in_data.EL(1,2)==4 | in_data.EL(1,2)==5 % "4/5" - CST/CSQ element (plain stress) 
    for i=1:size(in_data.CON)
        node_i = find(in_data.ND(:,1)==in_data.CON(i,1));
        if in_data.CON(i,2)==0
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)-labx],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],'b-','LineWidth',3); hold on;
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)-labx],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],'b-','LineWidth',3); hold on;
        end;
        if in_data.CON(i,3)==0
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)-laby],'b-','LineWidth',3); hold on;
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)-laby],'b-','LineWidth',3); hold on;
        end;
    end;
    hold off;
end;

if in_data.EL(i,2)==3  % "3" - 3D-beam element
    maxZ = max(in_data.ND(:,4));
    labz = (maxZ / 8);
    for i=1:size(in_data.CON)
        node_i = find(in_data.ND(:,1)==in_data.CON(i,1));
        if in_data.CON(i,2)==0
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)-labx],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)],'b-','LineWidth',3); hold on;
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)-labx],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)],'b-','LineWidth',3); hold on;
        end;
        if in_data.CON(i,3)==0
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)-laby],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)],'b-','LineWidth',3); hold on;
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)-laby],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)],'b-','LineWidth',3); hold on;
        end;
        if in_data.CON(i,4)==0
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)-labz],'b-','LineWidth',3); hold on;
            plot3([in_data.ND(node_i,2) in_data.ND(node_i,2)],...
                [in_data.ND(node_i,3) in_data.ND(node_i,3)],...
                [in_data.ND(node_i,4) in_data.ND(node_i,4)-labz],'b-','LineWidth',3); hold on;
        end;
    end;
    hold off;
end;

if in_data.EL(i,2)==0 | in_data.EL(i,2)==1 | in_data.EL(i,2)==2  % "0/1/2" - 2D-beam element
    for i=1:size(in_data.CON)
        node_i = find(in_data.ND(:,1)==in_data.CON(i,1));
        if in_data.CON(i,2)==0
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)-labx],[in_data.ND(node_i,3) in_data.ND(node_i,3)],'b-',...
                'LineWidth',3); hold on;
        end;
        if in_data.CON(i,3)==0
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)],[in_data.ND(node_i,3) in_data.ND(node_i,3)-laby],'b-',...
                'LineWidth',3); hold on;
        end;
        if in_data.CON(i,4)==0
            plot([in_data.ND(node_i,2) in_data.ND(node_i,2)],[in_data.ND(node_i,3) in_data.ND(node_i,3)],'bs',...
                'LineWidth',3); hold on;
        end;
    end;
end;