fem_run.m

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

function [obj,resp] = fem_run (in_data,type_analysis)

% =============================================================================
% FEM analysis
% Anton Zaicenco <[e-mail]>, started: Dec 2004
% updated: May 2006
% =============================================================================

% DATA STRUCTURES:
% obj:
%     - data about the the analyzed object: K, M
% in_data:
%     - input data: nodes, elements, material properties, constrains, etc.
% resp:
%     - structural response: resp.static.D & resp.dynamic.[DYN_D DYN_V DYN_A]

% assembling global stiffness and mass matrices ===============================
% dofN - dof's per node
if in_data.EL(1,2)==0 | in_data.EL(1,2)==1 | in_data.EL(1,2)==2 dofN=3; end; % 2D-beam
if in_data.EL(1,2)==3 dofN = 6; end; % 3D-beam  elements
if in_data.EL(1,2)==4 dofN = 2; end; % CST  elements
if in_data.EL(1,2)==5 dofN = 2; end; % CSQ  elements
if in_data.EL(1,2)==6 dofN = 3; end; % "6" - 3D-BRICK element (8-nodes)
if in_data.EL(1,2)==9 dofN = 3; end; % BCIZ elements

dof_ = size(in_data.ND,1)*dofN; % total dof
[obj.Ksys.Kgl, obj.M]    = K_assembly (in_data);
[obj.Ksys.Kglobal,L,Rch] = K_CON (dofN, dof_, obj, in_data);

% ========================== STATIC ANALYSIS ==================================
if type_analysis(1)==1
    if size(in_data.LOAD_)~=0
        [LVEC] = LOAD_V (dofN, dof_, in_data, L);    % form the load vector
        displacements = inv(obj.Ksys.Kglobal)*LVEC'; % compute displacements
        resp.static.D(L) = displacements;
        SF_ = zeros(size(in_data.ND,1), dofN+1);
        SF_(1:size(in_data.ND,1),1)=[1:size(in_data.ND,1)]';
        for i=1:size(in_data.ND,1)
            SF_(i,dofN-(dofN-1)+1:dofN+1) = [resp.static.D(i*dofN-(dofN-1):i*dofN)];
        end;
        save FEM_STATIC_DISPL.txt SF_ -ascii;
        disp('--- Static Displacements saved to "FEM_STATIC_DISPL.txt": [elem# d1 ... dn]');
        clear displacements
        plot_fem_static(in_data, resp);
    end;
end;
% =============================== END =========================================

% ========================== DYNAMIC ANALYSIS =================================
if type_analysis(2)==1
    if size(in_data.dynam.TIMEH)~=0
        if in_data.EL(1,2)==9 | in_data.EL(1,2)==4 | in_data.EL(1,2)==5 | ...
                in_data.EL(1,2)==3 | in_data.EL(1,2)==6 | in_data.EL(1,2)==0 ...
                | in_data.EL(1,2)==1 | in_data.EL(1,2)==2
            [resp.dynamic.DYN_D, resp.dynamic.DYN_V, resp.dynamic.DYN_A, t] = ...
                dyn_FEM ( in_data, obj, L);
        end;
        zdd    = load(in_data.dynam.TIMEH);
        RESP_ = resp.dynamic.DYN_A(:,in_data.dynam.TIMEHPL(1));
        [FT_T,omega] = FFTsig (zdd(1:(end-1)),RESP_,abs(t(2)-t(1)));
        save FEM_DYN_RESP.txt RESP_ -ascii;
        disp('--- Dynamic response saved to "FEM_DYN_RESP.txt"');
    end;
end;
% =============================== END =========================================

% ========================== MODAL ANALYSIS ===================================
if type_analysis(3)==1
    if size(in_data.dynam.MODA)~=0 & size(in_data.dynam.DAMP_C)~=0
        if in_data.EL(1,2)==4 | in_data.EL(1,2)==5 | in_data.EL(1,2)==3 | in_data.EL(1,2)==6 | ...
                in_data.EL(1,2)==0 | in_data.EL(1,2)==1 | in_data.EL(1,2)==2 | in_data.EL(1,2)==9
            [E_Vec_ok, wn] = modal_analysis ( obj , L, in_data);
            i = size(E_Vec_ok);
            plot_modal_shape(in_data,E_Vec_ok(:,i(2)+1-in_data.dynam.MODA)',dof_,33); % plot mode
            mode_T_sec = 2*pi/wn(in_data.dynam.MODA);
            title(['Mode shape:  '  num2str(in_data.dynam.MODA) '.  T =  ' num2str(mode_T_sec) '  s']);
            disp('--------- MODAL ANALYSIS DONE ---------'); resp.modal='';
        end;
    end;
end;
% =============================== END =========================================