function exam4_1( m, n )
% 本程序為第四章的第一個算例，采用矩形單元計算純彎梁的變形
%      exam4_1(m,n) 
%  輸入參數： 
%      m ------  x方向單元數目
%      n ------  y方向單元數目

% 定義全局變量
%      gNode ------ 節點坐標
%      gElement --- 單元定義
%      gMaterial -- 材料性質
%      gBC1 ------- 第一類約束條件
%      gDF -------- 分布力
%      gK --------- 整體剛度矩陣
%      gDelta ----- 整體節點坐標
    global gNode gElement gMaterial gBC1 gDF gK gDelta

    if nargin < 1 
        m = 4 ;
        n = 4 ;
    elseif nargin < 2
        n = 4 ;
    end
    
    FemModel(m, n) ;       % 定義有限元模型
    SolveModel ;           % 求解有限元模型
    DisplayResults ;       % 顯示計算結果
return ;

function FemModel(m, n)
%  定義有限元模型
%  輸入參數：
%      m ---  x方向單元數目
%      n ---  y方向單元數目
%  返回值：
%      無
%  說明：
%      該函數定義平面桿系的有限元模型數據：
%        gNode ------- 節點定義
%        gElement ---- 單元定義
%        gMaterial --- 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩
%        gBC --------- 約束條件
%        gDF --------- 分布力

    global gNode gElement gMaterial gBC1 gDF

    length = 0.045 ;    % 計算部分的長度(x方向)
    height = 0.03 ;     % 計算部分的高度(y方向)
    dx = length/m ;     % 矩形單元的寬度
    dy = height/n ;     % 矩形單元的高度
    p = 10e6 ;          % 彎曲應力10MPa
    
    % 節點坐標
    gNode = zeros( (m+1)*(n+1), 2 ) ;
    for i=1:n+1
        for j=1:m+1
            k = (i-1)*(m+1)+j ;        % 節點號
            xk = (j-1)*dx ;            % 節點的x坐標
            yk = (i-1)*dy ;            % 節點的y坐標
            gNode(k,:) = [xk, yk] ;
        end
    end
     
    % 單元定義
    gElement = zeros( m*n, 4 ) ;
    for i=1:n
        for j=1:m
            k = (i-1)*m+j ;             % 單元號
            n1 = (i-1)*(m+1)+j ;        % 第一個節點號
            n2 = (i-1)*(m+1)+j+1 ;      % 第二個節點號
            n3 = i*(m+1)+j+1 ;          % 第三個節點號
            n4 = i*(m+1)+j ;            % 第四個節點號
            gElement(k,:) = [n1, n2, n3, n4] ;
        end
    end

    % 材料性質 
    %           彈性模量    泊松比   厚度
    gMaterial = [2.0e11,    0.3,     0.01] ;   %  材料 1

    % 第一類約束條件
    gBC1 = zeros( m+1+n+1, 3 ) ;
    for j=1:m+1
        gBC1(j,:) = [j, 1, 0.0] ;
    end
    for i=2:n+1
        gBC1(m+1+i,:) = [(i-1)*(m+1)+1, 1, 0.0] ; 
    end
    gBC1(m+1+1,:) = [1, 2, 0.0] ;

    % 分布載荷（線性分布）
    gDF = zeros( n, 5 ) ;
    for i=1:n
        k = i*m ;
        gDF(i,:) = [ k, 2, (i-1)*p/n, i*p/n, 1] ;
    end
return

function SolveModel
%  求解有限元模型
%  輸入參數：
%     無
%  返回值：
%     無
%  說明：
%      該函數求解有限元模型，過程如下
%        1. 計算單元剛度矩陣，集成整體剛度矩陣
%        2. 計算單元的等效節點力，集成整體節點力向量
%        3. 處理約束條件，修改整體剛度矩陣和節點力向量
%        4. 求解方程組，得到整體節點位移向量

    global gNode gElement gMaterial gBC1 gDF gK gDelta

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

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

    % step3. 計算分布力的等效節點力，并集成到整體節點力向量中
    [df_number, dummy] = size( gDF ) ;
    for idf = 1:1:df_number
        enf = EquivalentNodeForce( gDF(idf,1), gDF(idf,2), gDF(idf,3), gDF(idf,4), gDF(idf,5) ) ;
        ielem = gDF(idf,1) ;
        iedge = gDF(idf,2) ;
        i = gElement( ielem, iedge ) ;
        if iedge < 4 
            j = gElement( ielem, iedge+1 ) ;
        else
            j = gElement( ielem, 1 );
        end
        
        f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + enf( 1:2 ) ;
        f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + enf( 3:4 ) ;
    end
  
    % step4. 處理約束條件，修改剛度矩陣和節點力向量。采用乘大數法
    [bc_number,dummy] = size( gBC1 ) ;
    for ibc=1:1:bc_number
        n = gBC1(ibc, 1 ) ;
        d = gBC1(ibc, 2 ) ;
        m = (n-1)*2 + d ;
        f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
        gK(m,m) = gK(m,m) * 1e15 ;
    end

    % step 5. 求解方程組，得到節點位移向量
    gDelta = gK \ f ;
return


function k = StiffnessMatrix( ie )
%  計算單元剛度矩陣
%  輸入參數:
%     ie ----  單元號
%  返回值:
%     k  ----  單元剛度矩陣

global gNode gElement gMaterial 
    k = zeros( 8, 8 ) ;
    E  = gMaterial( 1 ) ;
    mu = gMaterial( 2 ) ;
    h  = gMaterial( 3 ) ;
    x1 = gNode( gElement( ie, 1 ), 1 ) ;
    y1 = gNode( gElement( ie, 1 ), 2 ) ;
    x3 = gNode( gElement( ie, 3 ), 1 ) ;
    y3 = gNode( gElement( ie, 3 ), 2 ) ;
    a = (x3-x1)/2 ;
    b = (y3-y1)/2 ;
    xi  = [ -1, 1, 1, -1 ] ;
    eta = [ -1, -1, 1, 1 ] ;
    for i=1:1:4
        for j=1:1:4
            k( (i-1)*2 + 1, (j-1)*2 + 1 ) = b/a*xi(i)*xi(j)*(1+1/3*eta(i)*eta(j)) + ...
                (1-mu)/2*a/b*eta(i)*eta(j)*(1+1/3*xi(i)*xi(j)); 
            k( (i-1)*2 + 1, (j-1)*2 + 2 ) = mu*eta(j)*xi(i) + (1-mu)/2*xi(j)*eta(i) ;
            k( (i-1)*2 + 2, (j-1)*2 + 1 ) = mu*xi(j)*eta(i) + (1-mu)/2*eta(j)*xi(i) ;
            k( (i-1)*2 + 2, (j-1)*2 + 2 ) = a/b*eta(i)*eta(j)*(1+1/3*xi(i)*xi(j)) + ...
                (1-mu)/2*b/a*xi(i)*xi(j)*(1+1/3*eta(i)*eta(j)) ;
        end
    end
    eh = E*h/4/(1-mu^2) ;
    k = eh*k ;
return

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

function enf = EquivalentNodeForce( ie, iedge, p1, p2, idof )
%   計算線性分布荷載的等效節點力
%   輸入參數:
%      ie  -----  單元號
%      ieged ---  作用的邊號
%      p1  -----  第一個節點上的分布力集度值
%      p2  -----  第二個節點上的分布力集度值
%      idof  ---  分布力的方向
%                  1 ---  x方向
%                  2 ---  y方向
%   返回值:
%      enf -----  等效節點力向量
    global gElement gNode gMaterial
    h = gMaterial( 3 ) ;
    x1 = gNode( gElement( ie, 1 ), 1 ) ;
    y1 = gNode( gElement( ie, 1 ), 2 ) ;
    x3 = gNode( gElement( ie, 3 ), 1 ) ;
    y3 = gNode( gElement( ie, 3 ), 2 ) ;
    a = ( x3 - x1 ) / 2 ;
    b = ( y3 - y1 ) / 2 ;
    if iedge == 1 | iedge == 3
        length = a ;        
    else
        length = b ;
    end
    f1 = h*length*(2*p1+p2)/3 ;
    f2 = h*length*(p1+2*p2)/3 ;
    if idof == 1
        enf = [f1; 0; f2; 0] ;
    else
        enf = [0; f1; 0; f2] ;
    end
return


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

    global gNode gDelta
    
    fprintf( '節點位移\n' ) ; 
    fprintf( '  節點號         x方向位移               y方向位移\n' ) ; 
    [node_number,dummy] = size( gNode ) ;
    for i=1:node_number
        fprintf(  '%6d       %16.8e        %16.8e\n',...
                  i, gDelta((i-1)*2+1), gDelta((i-1)*2+2) ) ; 
    end
return