function exam6_2
% 本程序為第六章的第二個算例，采用20結點六面體等參元計算端部受橫向力的懸臂梁
%  輸入參數： 
%       無

    file_in = 'exam6_2.dat' ;
    
    % 檢查文件是否存在
    if exist( file_in ) == 0
        disp( sprintf( '錯誤：文件 %s 不存在', file_in ) )
        disp( sprintf( '程序終止' ) )
        return ;
    end
    
    % 根據輸入文件名稱生成輸出文件名稱
    [path_str,name_str,ext_str] = fileparts( file_in ) ;
    ext_str_out = '.mat' ;
    file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
    
    % 檢查輸出文件是否存在
    if exist( file_out ) ~= 0 
        answer = input( sprintf( '文件 %s 已經存在，是否覆蓋? ( Yes / [No] ):  ', file_out ), 's' ) ;
        if length( answer ) == 0
            answer = 'no' ;
        end
        
        answer = lower( answer ) ;
        if answer ~= 'y' | answer ~= 'yes' 
            disp( sprintf( '請使用另外的文件名，或備份已有的文件' ) ) ;
            disp( sprintf( '程序終止' ) );
            return ; 
        end
    end

    FemModel( file_in ) ;       % 定義有限元模型
    SolveModel ;                % 求解有限元模型
    SaveResults( file_out ) ;   % 保存計算結果
    DisplayResults ;            % 把計算結果顯示出來 
return

function FemModel( file_in )
%  定義有限元模型
%  輸入參數：
%      file_in ---- 有限模型數據文件
%  返回值：
%      無

%  全局變量及名稱
%      gNode ------ 節點坐標
%      gElement --- 單元定義
%      gMaterial -- 材料性質
%      gBC1 ------- 第一類約束條件
%      gNF -------- 集中力
    global gNode gElement gMaterial gBC1 gNF
    
    % 打開文件
    fid = fopen( file_in, 'r' ) ;
    
    % 節點坐標
    node_number = fscanf( fid, '%d', 1 ) ;
    gNode = zeros( node_number, 3 ) ;
    for i = 1:node_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gNode( i, : ) = fscanf( fid, '%f', [1 3] ) ;
    end
    
    % 單元定義
    element_number = fscanf( fid, '%d', 1 ) ;
    gElement = zeros( element_number, 21 ) ;
    for i = 1:element_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gElement( i, : ) = fscanf( fid, '%d', [1 21] ) ;
    end
    
    % 材料性質
    material_number = fscanf( fid, '%d', 1 ) ;
    gMaterial = zeros( material_number, 2 ) ;
    for i = 1:material_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gMaterial( i, : ) = fscanf( fid, '%f', [1 2] ) ;
    end
    
    % 第一類約束條件
    bc1_number = fscanf( fid, '%d', 1 ) ;
    gBC1 = zeros( bc1_number, 3 ) ;
    for i=1:bc1_number
        gBC1(i,1) = fscanf( fid, '%d', 1 ) ;  % 結點號
        gBC1(i,2) = fscanf( fid, '%d', 1 ) ;  % 自由度號
        gBC1(i,3) = fscanf( fid, '%f', 1 ) ;  % 約束值
    end
    
    % 集中力
    nf_number = fscanf( fid, '%d', 1 ) ;
    gNF = zeros( nf_number, 3 ) ;
    for i=1:nf_number
        gNF(i,1) = fscanf( fid, '%d', 1 ) ;  % 結點號
        gNF(i,2) = fscanf( fid, '%d', 1 ) ;  % 自由度號
        gNF(i,3) = fscanf( fid, '%f', 1 ) ;  % 集中力大小
    end
    
    % 關閉文件
    fclose( fid ) ;
return

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

%  全局變量及名稱
%      gNode ----------- 節點坐標
%      gElement -------- 單元定義
%      gMaterial ------- 材料性質
%      gBC1 ------------ 第一類約束條件
%      gNF ------------- 集中力
%      gK -------------- 整體剛度矩陣
%      gDelta ---------- 結點位移向量
%      gElementStress -- 每個單元的結點應力
    global gNode gElement gMaterial gBC1 gNF gK gDelta gElementStress

    % step1. 定義整體剛度矩陣和節點力向量
    [node_number,dummy] = size( gNode ) ;
    %gK = sparse( node_number * 3, node_number * 3 ) ;   % 由于稀疏矩陣的存儲速度很慢，因此改成普通矩陣運算
    %f = sparse( node_number * 3, 1 ) ;                  % 如果計算機內存不夠，就恢復這兩句，并注釋掉下面兩句。
    gK = zeros( node_number*3,node_number*3) ;
     f = zeros( node_number*3,1 ) ;
     
    % step2. 計算單元剛度矩陣，并集成到整體剛度矩陣中
    [element_number,dummy] = size( gElement ) ;
    for ie=1:1:element_number
        disp( sprintf(  '計算單元剛度矩陣并集成，當前單元: %d', ie  ) ) ;
        k = StiffnessMatrix( ie ) ;
        AssembleStiffnessMatrix( ie, k ) ;
    end
    
    spy(gK) ;
    nz=nnz(gK) ;
    fprintf( 'non-zeros elements/total elmements = %g%%',nz/node_number/node_number/9*100) ;
    pause ;
    
    % step3. 把集中力集成到整體節點力向量中
    [nf_number, dummy] = size( gNF ) ;
    for idf = 1:1:nf_number
        node  = gNF(idf, 1) ;
        dim   = gNF(idf, 2) ;
        force = gNF(idf, 3) ;
        f( (node-1)*3+dim ) = f( (node-1)*3+dim ) + force ;
    end
    
    % step4. 處理第一類約束條件，修改剛度矩陣和節點力向量。采用乘大數法
    [bc1_number,dummy] = size( gBC1 ) ;
    for ibc=1:1:bc1_number
        n = gBC1(ibc, 1 ) ;
        d = gBC1(ibc, 2 ) ;
        m = (n-1)*3 + d ;
        f(m,:) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
        gK(m,m) = gK(m,m) * 1e15 ;
    end
    
    % step5. 求解方程組，得到節點位移向量
    gDelta = gK \ f ;
    
    % step6. 計算每個單元的結點應力
    gElementStress = zeros( element_number, 20, 6 ) ;
    delta = zeros( 60, 1 ) ;
    x = [-1, 1, 1, -1, -1, 1, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, -1, 1, 1, -1];
    e = [-1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0, -1, -1, 1, 1];
    z = [-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0];
    for ie = 1:element_number
        fprintf( '計算單元應力，當前單元號：%d\n', ie ) ;
        for n=1:20
            B = MatrixB( ie, x(n), e(n), z(n) ) ;
            D = MatrixD( ie ) ;
            delta(1:3:58) = gDelta( gElement(ie,1:20)*3-2) ;
            delta(2:3:59) = gDelta( gElement(ie,1:20)*3-1) ;
            delta(3:3:60) = gDelta( gElement(ie,1:20)*3) ;
            sigma = D*B*delta ;
            gElementStress( ie, n, : ) = sigma ;
        end
    end
return

function k = StiffnessMatrix( ie )   
%  計算平面應變等參數單元的剛度矩陣
%  輸入參數：
%     ie -- 單元號
%  返回值：
%     k  -- 單元剛度矩陣
%  說明：
%     用高斯積分法求解平面等參數單元的剛度矩陣

    k = zeros( 60, 60 ) ;
    D = MatrixD( ie ) ;
    %3 x 3 x 3 高斯積分點和權系數
    %x = [-0.774596669241483,              0.0,0.774596669241483] ;   
    %w = [ 0.555555555555556,0.888888888888889,0.555555555555556] ;  
    %for i=1:1:length(x)
    %    for j=1:1:length(x)
    %        for m=1:1:length(x)
    %            B = MatrixB( ie, x(i), x(j), x(m) ) ;
    %            J = Jacobi( ie, x(i), x(j), x(m) ) ;
    %            k = k + w(i)*w(j)*w(m)*transpose(B)*D*B*det(J) ;   
    %        end
    %    end
    %end
    
    %14點 irons 積分, 詳細參見以下文獻 
    % Irons B. M., 1971 Quadrature rules for brick based finite elements, 
    %   International Journal for Numerical Methods in Engineering, 
    %   3(2):293-294.
    b0 = 320/361 ;
    c0 = 121/361 ;
    b  = sqrt( 19/30 ) ;
    c  = sqrt( 19/33 ) ;
    x1 = [-b,  0,  0
           b,  0,  0
           0, -b,  0
           0,  b,  0
           0,  0, -b
           0,  0,  b ] ;
    x2 = [-c, -c, -c
           c,  c,  c
           c, -c, -c
          -c,  c, -c
          -c, -c,  c
           c,  c, -c
          -c,  c,  c
           c, -c,  c ] ;
    for i=1:6
        B = MatrixB( ie, x1(i,1), x1(i,2), x1(i,3) ) ;
        J = Jacobi( ie, x1(i,1), x1(i,2), x1(i,3) ) ;
        k = k + b0*transpose(B)*D*B*det(J) ;
    end
    for i=1:8
        B = MatrixB( ie, x2(i,1), x2(i,2), x2(i,3) ) ;
        J = Jacobi( ie, x2(i,1), x2(i,2), x2(i,3) ) ;
        k = k + c0*transpose(B)*D*B*det(J) ;
    end
return

function D = MatrixD( ie )
%  計算單元的彈性矩陣D
%  輸入參數：
%     ie --------- 單元號
%  返回值：
%     D  --------- 彈性矩陣D

%  全局變量及名稱
%      gElement -------- 單元定義
%      gMaterial ------- 材料性質
    global gElement gMaterial 

    E  = gMaterial( gElement(ie, 21), 1 ) ;   % 彈性模量
    mu = gMaterial( gElement(ie, 21), 2 ) ;   % 泊松比 
    A1 = mu/(1-mu) ;                        
    A2 = (1-2*mu)/2/(1-mu) ;                
    A3 = E*(1-mu)/(1+mu)/(1-2*mu) ;         
    D = A3*[  1  A1  A1   0   0   0
             A1   1  A1   0   0   0
             A1  A1   1   0   0   0
              0   0   0  A2   0   0
              0   0   0   0  A2   0
              0   0   0   0   0  A2] ;
return

function B = MatrixB( ie, xi, eta, zeta )
%  計算單元的應變矩陣B
%  輸入參數：
%     ie --------------- 單元號
%     xi,eta, zeta ----- 局部坐標  
%  返回值：
%     B  --------------- 在局部坐標處的應變矩陣B

    [N_x, N_y, N_z] = N_xyz( ie, xi, eta, zeta ); 
    B = zeros( 6, 60 ) ;
    for i=1:1:20
        B(:,(3*i-2):3*i) = [ N_x(i)       0       0
                                  0  N_y(i)       0
                                  0       0  N_z(i)
                                  0  N_z(i)  N_y(i)
                             N_z(i)      0   N_x(i)
                             N_y(i)  N_x(i)       0 ] ;
    end
return

function [N_x, N_y, N_z] = N_xyz( ie, xi, eta, zeta )
%  計算形函數對整體坐標的導數