function exam8_2
% 本程序為第八章的第二個算例，采用平面梁單元計算兩鉸拋物線拱的在初始條件下
%  自由振動，并對時程曲線結果進行FFT變換，求得的頻率可與exam8_1.m的結果進行
%  比較，以驗證本程序的可靠性
%      輸入參數： 無
%      輸出結果： 位移的時程曲線及其頻譜特性圖 

    PlaneFrameModel ;             % 定義有限元模型
    SolveModel ;                  % 求解有限元模型
    SaveResults('exam8_2.mat') ;  % 保存計算結果
    DisplayResults ;              % 顯示計算結果
return ;

function PlaneFrameModel
%  定義平面桿系的有限元模型
%  輸入參數：
%      無
%  返回值：
%      無
%  說明：
%      該函數定義平面桿系的有限元模型數據：
%        gNode -------- 節點定義
%        gElement ----- 單元定義
%        gMaterial ---- 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩
%        gBC1 --------- 約束條件
%        gDeltaT ------ 時間步長
%        gTimeEnd ----- 計算結束時刻
%        gDisp -------- 位移時程響應
%        gVelo -------- 速度時程響應
%        gAcce -------- 加速度時程響應

    global gNode gElement gMaterial gBC1 gDeltaT gTimeEnd gDisp gVelo gAcce

    % 給定拋物線拱的幾何特征
    L = 60 ;               %  計算跨徑(m)     
    f = 7.5 ;              %  計算矢高(m)
    
    n = 100 ;              %  單元數目
    x = -L/2:L/n:L/2 ;     %  結點的x坐標
    a = f/L^2*4 ;
    y = - a * x.^2 ;       %  結點的y坐標

    % 節點坐標
    gNode = [x'  y'] ;
    
    % 單元定義
    gElement = zeros( n, 3 ) ;
    for i=1:n
        gElement( i, : ) = [ i, i+1, 1 ] ;
    end
    
    % 材料性質 
    %           彈性模量   抗彎慣性矩   截面積   密度
    gMaterial = [2.06e11,  0.03622,   0.0815,  1435.2/0.0815];   %  材料 1

    % 第一類約束條件
    %     節點號   自由度號    約束值
    gBC1 = [ 1,        1,        0.0
             1,        2,        0.0
             n+1,      1,        0.0
             n+1,      2,        0.0] ;
     
    gDeltaT = 0.01 ;
    gTimeEnd = 4096*gDeltaT  ;    % 計算時間為載荷通過所需時間的兩倍
    timestep = floor(gTimeEnd/gDeltaT) ;

    % 定義位移，速度和加速度
    gDisp = zeros( (n+1)*3, timestep ) ;
    gVelo = zeros( (n+1)*3, timestep ) ;
    gAcce = zeros( (n+1)*3, timestep ) ;
    
    % 初始條件
    gDisp(:,1) = zeros( (n+1)*3, 1 ) ;
    gVelo(:,1) = ones( (n+1)*3, 1 ) ;
return

function SolveModel
%  求解有限元模型
%  輸入參數：
%     無
%  返回值：
%     無
%  說明：
%      該函數求解有限元模型，過程如下
%        1. 計算單元的剛度和質量矩陣，集成整體剛度和質量矩陣
%        2. 用Newmark法計算時程響應

    global gNode gElement gMaterial gBC1 gK gM gDeltaT gTimeEnd gDisp gVelo gAcce

    % step1. 定義整體剛度矩陣和節點力向量
    [node_number,dummy] = size( gNode ) ;
    gK = sparse( node_number * 3, node_number * 3 ) ;
    gM = sparse( node_number * 3, node_number * 3 ) ;

    % step2. 計算單元剛度和質量矩陣，并集成到整體剛度和質量矩陣中
    [element_number,dummy] = size( gElement ) ;
    for ie=1:1:element_number
        k = StiffnessMatrix( ie ) ;
        m = MassMatrix( ie ) ; 
        AssembleGlobalMatrix( ie, k, m ) ;
    end

    % step3. 計算時程響應(Newmark法)
    % step3.1 初始計算
    gama = 0.5 ;
    beta = 0.25 ;
    C = zeros( size( gK ) ) ;
    [N,N] = size( gK ) ;
    alpha0 = 1/beta/gDeltaT^2 ;
    alpha1 = gama/beta/gDeltaT ;
    alpha2 = 1/beta/gDeltaT ;
    alpha3 = 1/2/beta - 1 ;
    alpha4 = gama/beta - 1 ;
    alpha5 = gDeltaT/2*(gama/beta-2) ;
    alpha6 = gDeltaT*(1-gama) ;
    alpha7 = gama*gDeltaT ;
    K1 = gK + alpha0*gM + alpha1*C;
    timestep = floor(gTimeEnd/gDeltaT) ;
    
    % step3.2 對K1進行邊界條件處理
    [bc1_number,dummy] = size( gBC1 ) ;
    K1im = zeros(N,bc1_number) ;
    for ibc=1:1:bc1_number
        n = gBC1(ibc, 1 ) ;
        d = gBC1(ibc, 2 ) ;
        m = (n-1)*3 + d ;
        K1im(:,ibc) = K1(:,m) ; 
        K1(:,m) = zeros( node_number*3, 1 ) ;
        K1(m,:) = zeros( 1, node_number*3 ) ;
        K1(m,m) = 1.0 ;
    end
    [KL,KU] = lu(K1) ;   % 進行三角分解，節省后面的求解時間
    
    % step3.3 計算初始加速度
    gAcce(:,1) = gM\(-gK*gDisp(:,1)-C*gVelo(:,1)) ;
    
    % step3.4 對每一個時間步計算
    for i=2:1:timestep
        if mod(i,100) == 0
            fprintf( '當前時間步：%d\n', i ) ;
        end
        f1 = gM*(alpha0*gDisp(:,i-1)+alpha2*gVelo(:,i-1)+alpha3*gAcce(:,i-1)) ...
                  + C*(alpha1*gDisp(:,i-1)+alpha4*gVelo(:,i-1)+alpha5*gAcce(:,i-1)) ;
        % 對f1進行邊界條件處理
        [bc1_number,dummy] = size( gBC1 ) ;
        for ibc=1:1:bc1_number
            n = gBC1(ibc, 1 ) ;
            d = gBC1(ibc, 2 ) ;
            m = (n-1)*3 + d ;
            f1 = f1 - gBC1(ibc,3) * K1im(:,ibc) ;
            f1(m) = gBC1(ibc,3) ;
        end
        y = KL\f1 ;
        gDisp(:,i) = KU\y ;
        gAcce(:,i) = alpha0*(gDisp(:,i)-gDisp(:,i-1)) - alpha2*gVelo(:,i-1) - alpha3*gAcce(:,i-1) ;
        gVelo(:,i) = gVelo(:,i-1) + alpha6*gAcce(:,i-1) + alpha7*gAcce(:,i) ;
    end
return

function k = StiffnessMatrix( ie )
%  計算單元剛度矩陣
%  輸入參數:
%     ie -------  單元號
%  返回值:
%     k  ----  整體坐標系下的剛度矩陣
    global gNode gElement gMaterial
    k = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    I = gMaterial( gElement(ie, 3), 2 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    k = [  E*A/L           0          0 -E*A/L           0          0
               0  12*E*I/L^3  6*E*I/L^2      0 -12*E*I/L^3  6*E*I/L^2
               0   6*E*I/L^2    4*E*I/L      0  -6*E*I/L^2    2*E*I/L
          -E*A/L           0          0  E*A/L           0          0
               0 -12*E*I/L^3 -6*E*I/L^2      0  12*E*I/L^3 -6*E*I/L^2
               0   6*E*I/L^2    2*E*I/L      0  -6*E*I/L^2    4*E*I/L] ;
    T = TransformMatrix( ie ) ;
    k = T*k*transpose(T) ;
return

function m = MassMatrix( ie )
%  計算單元質量矩陣
%  輸入參數:
%     ie -------  單元號
%  返回值:
%     m  ----  整體坐標系下的質量矩陣
    global gNode gElement gMaterial
    m = zeros( 6, 6 ) ;
    E = gMaterial( gElement(ie, 3), 1 ) ;
    A = gMaterial( gElement(ie, 3), 3 ) ;
    ro = gMaterial( gElement(ie, 3 ), 4 ) ;
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
    m = ro*A*L/420*[140      0      0   70      0      0
                      0    156   22*L    0     54   -13*L
                      0   22*L  4*L^2    0   13*L  -3*L^2
                     70      0      0  140      0       0 
                      0     54   13*L    0    156   -22*L
                      0  -13*L -3*L^2    0  -22*L  4*L^2 ] ;
    T = TransformMatrix( ie ) ;
    m = T*m*transpose(T) ;
return

function AssembleGlobalMatrix( ie, ke, me )
%  把單元剛度和質量矩陣集成到整體剛度矩陣
%  輸入參數:
%      ie  --- 單元號
%      ke  --- 單元剛度矩陣
%      me  --- 單元質量矩陣
%  返回值:
%      無
    global gElement gK gM
    for i=1:1:2
        for j=1:1:2
            for p=1:1:3
                for q =1:1:3
                    m = (i-1)*3+p ;
                    n = (j-1)*3+q ;
                    M = (gElement(ie,i)-1)*3+p ;
                    N = (gElement(ie,j)-1)*3+q ;
                    gK(M,N) = gK(M,N) + ke(m,n) ;
                    gM(M,N) = gM(M,N) + me(m,n) ;
                end
            end
        end
    end
return

function T = TransformMatrix( ie )
%  計算單元的坐標轉換矩陣( 局部坐標 -> 整體坐標 )
%  輸入參數
%      ie  ----- 節點號
%  返回值
%      T ------- 從局部坐標到整體坐標的坐標轉換矩陣
    global gElement gNode
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
    c = (xj-xi)/L ;
    s = (yj-yi)/L ;
    T=[ c  -s   0   0   0   0
        s   c   0   0   0   0
        0   0   1   0   0   0
        0   0   0   c  -s   0
        0   0   0   s   c   0
        0   0   0   0   0   1] ;
return

function SaveResults( file_out )
%  保存計算結果
%  輸入參數：
%     無
%  返回值：
%     無
    global gNode gElement gMaterial gBC1 gDeltaT gTimeEnd gLoad gLoadVelo gDisp gVelo gAcce
    save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1',  ...
          'gDeltaT', 'gTimeEnd', 'gLoad', 'gLoadVelo', 'gDisp', 'gVelo', 'gAcce' ) ;
return

function DisplayResults
%  顯示計算結果
%  輸入參數：
%     無
%  返回值：
%     無

    global gNode gElement gMaterial gBC1 gDisp gVelo gAcce gDeltaT gTimeEnd

    % 繪制時程曲線
    [node_number,dummy] = size(gNode) ;
    t = 0:gDeltaT:gTimeEnd-gDeltaT ;
    d = gDisp((floor(node_number/4)*3)+2,:) ;
    subplot(2,1,1) ;
    plot( t, d ) ;
    title( 'L/4處撓度時程曲線' ) ;
    xlabel( '時間(s)') ;
    ylabel( '撓度(m)' ) ;
    
    % 對時程曲線進行FFT變換，獲取頻譜特性
    fd = fft( d ) ;
    df = 1/gTimeEnd ;
    f = (0:length(d)-1)*df ;
    subplot(2,1,2);
    plot(f,abs(fd)) ;
    set(gca,'xlim',[2,10]) ;
    title( 'L/4處撓度的頻譜圖' ) ;
    xlabel( '頻率(Hz)') ;
    ylabel( '幅值' ) ;
    
    % 標注頻率峰值
    fifi1 = diff(abs(fd));
    n = length(fifi1) ;
    d1 = fifi1(1:n-1);
    d2 = fifi1(2:n) ;
    indmax = find( d1.*d2<0 & d1>0 )+1;
    for i=1:length(indmax)
        if f(indmax(i)) > 10 
            break ;
        end
        text( f(indmax(i)+2), abs(fd(indmax(i)))*0.9, sprintf('f=%.3f',f(indmax(i))));
    end
return