function Groundoverload_2111238(file_in)
% 本程序采用常应变三角形单元计算地面在局部超载作用下的变形和应力
% 全局变量
%      gNode ------------- 节点坐标
%      gElement ---------- 单元定义
%      gMaterial --------- 材料性质
%      gBC1 -------------- 约束条件
%      gDF --------------- 分布力
%      gK ---------------- 整体刚度矩阵
%      gDelta ------------ 整体节点坐标
%      gNodeStress ------- 节点应力
%      gElementStress ---- 整体初应力
%      gElementStrain ---- 整体应变

    global       gDelta    

    if nargin < 1
        file_in = 'Groundoverload.dat' ;
    end
    
    % 检查文件是否存在
    if exist( file_in ) == 0
         fprintf( '错误：文件 %s 不存在\n', file_in )  
         fprintf( '程序终止\n' ) 
        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 ) == 1
            answer = 'no' ;
        end
        
        answer = lower( answer ) ;
        if answer ~= 'y' | answer ~= 'yes' 
             fprintf( '请使用另外的文件名，或备份已有的文件\n' )  ;
             fprintf( '程序终止\n' ) ;
            return ; 
        end
    end

    % 建立有限元模型并求解，保存结果
    FemModel( file_in ) ;          % 定义有限元模型
    SolveModel ;                   % 求解有限元模型
    SaveResults( file_out ) ;      % 保存计算结果
    
     DisplayModel;
     PlotDisplacement( 2 );
     PlotDisplacementContour( 2, 10, 'r' );
   
    % 计算结束
     fprintf( '荷载中心点竖直位移 %f m\n', full( gDelta(170,1) ) ) ;
     fprintf( '计算正常结束，结果保存在文件 %s 中\n', file_out )  ;
return ;

function FemModel(filename)
%  定义有限元模型
%      该函数定义平面问题的有限元模型数据：
%        gNode ------- 节点定义
%        gElement ---- 单元定义
%        gMaterial --- 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩
%        gBC1 -------- 约束条件
%        gDF --------- 分布力

    global gNode gElement gMaterial gBC1 gDF
    
    % 打开文件
    fid = fopen( filename, 'r' ) ;
    
    % 读取节点坐标
    node_number = fscanf( fid, '%d', 1 ) ;
    gNode = zeros( node_number, 2 ) ;
    for i=1:node_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gNode( i, : ) = fscanf( fid, '%f', [1, 2] ) ;
    end
    
    % 读取单元定义
    element_number = fscanf( fid, '%d', 1 ) ;
    gElement = zeros( element_number, 4 ) ;
    for i=1:element_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gElement( i, : ) = fscanf( fid, '%d', [1, 4] ) ;
    end
	
    % 读取荷载信息
    load_number = fscanf( fid, '%d', 1 ) ;
    gDF = zeros( load_number, 4 ) ;
    for i=1:load_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gDF( i, : ) = fscanf( fid, '%f', [1, 4] ) ;
    end
	
    % 读取材料信息
    material_number = fscanf( fid, '%d', 1 ) ;
    gMaterial = zeros( material_number, 4 ) ;
    for i=1:material_number
        dummy = fscanf( fid, '%d', 1 ) ;
        gMaterial( i, : ) = fscanf( fid, '%f', [1,4] ) ;
    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
    
    % 关闭文件
    fclose( fid ) ;
return

function SolveModel
%  求解有限元模型
%      该函数求解有限元模型，过程如下
%        1. 计算单元刚度矩阵，集成整体刚度矩阵
%        2. 计算单元的等效节点力，集成整体节点力向量
%        3. 处理约束条件，修改整体刚度矩阵和节点力向量
%        4. 求解方程组，得到整体节点位移向量
%        5. 计算单元应力和节点应力

    global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF gElementStrain gFE

    %% step 1. 定义整体刚度矩阵和节点力向量
    [node_number,dummy] = size( gNode ) ;
    gK = sparse( node_number * 2, node_number * 2 ) ;  %系统整体刚度矩阵
    gFE = sparse( node_number * 2, 1 ) ;               %整体内力向量
    f = sparse( node_number * 2, 1 ) ;                 %节点等效荷载向量
    
    %% step 2. 定义初始全局应力应变矩阵
    [element_number, dunmmy] = size( gElement ) ;
    gElementStrain = zeros( element_number, 3) ;       %整体应变矩阵
    gElementStress = zeros( element_number, 3) ;       %整体应力矩阵
    
	  %% step 3. 计算地面超载产生的等效节点力
    [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) )*20;% 乘以10是改变力的大小
        i = gElement( gDF(idf,1), 1 ) ;
        j = gElement( gDF(idf,1), 2 ) ;
        m = gElement( gDF(idf,1), 3 ) ;
        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 ) ;
        f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + enf( 5:6 ) ;
    end 
 
    %% step 5 切线刚度法迭代 
    gDelta=zeros(node_number * 2,1);  %整体节点位移列阵 取初值delta0=0    
    js=0;                             %迭代次数统计指标
    while true
        gK=zeros( node_number * 2, node_number * 2 ); %整体刚度矩阵
        gFE=zeros(node_number * 2,1);                 %整体等效节点内力列阵
        for ie=1:1:element_number       
            delta = NodeDe( ie,gDelta); %整体节点位移列阵中提取单元的节点位移
            eps = MatrixB( ie ) * delta; %几何方程：公式2求epsilon0                      
            sigma0 = unlinerD(ie,eps)* (eps-gElementStrain(ie,:)')+gElementStress(ie,:)';%g公式3求sigmma0
            h  = gMaterial( gElement(ie, 4), 3 ) ;
            B = MatrixB( ie );
            area = ElementArea( ie );
            kt = transpose(B)*unlinerD(ie,eps)*B*h*abs(area) ; %公式4   
            AssembleStiffnessMatrix( ie, kt ) ;
            for ii = 1:1:3
                gElementStrain(ie,ii) = eps(ii);
                gElementStress(ie,ii) = sigma0(ii);
            end
            FE = elementforce( ie ,gElementStress) ;% 单元内力列阵
            AssembleFE( ie, FE ) ; % 整体内力列阵gFE
        end
        Phi=( f-gFE );
%% step 4. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法
        [bc_number,dummy] = size( gBC1 ) ;
        for ibc=1:1:bc_number
            n = gBC1(ibc, 1 ) ;
            d = gBC1(ibc, 2 ) ;
            m = (n-1)*2 + d ;
            Phi(m) = gBC1(ibc, 3)* gK(m,m) * 1e20 ;
            gK(m,m) = gK(m,m) * 1e20 ;
        end
        dd = gK \ Phi;  %公式5 位移增量
        gDelta=dd+gDelta;
        conv=norm(Phi)/norm(f);
        fprintf('迭代次数：%d   不平衡力/总外力：%e \n',js,conv)
        if js>100 || conv<1e-8            %判断是否达到精度要求
            break
        else        
        js=js+1;
        end
    end
   
    
    %% step 6. 计算节点应力(采用绕节点加权平均)
    gNodeStress = zeros( node_number, 6 ) ;       
    for i=1:node_number                         
        S = zeros( 1, 3 ) ;                         
        A = 0 ;
        for ie=1:1:element_number
            for k=1:1:3
                if i == gElement( ie, k ) 
                    area= ElementArea( ie ) ;
                    S = S + gElementStress(ie,1:3 ) * area ;
                    A = A + area ;
                    break ;
                end
            end
        end
        gNodeStress(i,1:3) = S / A ;
        gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
        gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
        gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
    end
    
return

  
function B = MatrixB( ie )
%  计算单元的应变矩阵B
%  输入参数:
%     ie ----  单元号
%  返回值:
%     B  ----  单元应变矩阵
    global gNode gElement
    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    xm = gNode( gElement( ie, 3 ), 1 ) ;
    ym = gNode( gElement( ie, 3 ), 2 ) ;
    ai = xj*ym - xm*yj ;
    aj = xm*yi - xi*ym ;
    am = xi*yj - xj*yi ;
    bi = yj - ym ;
    bj = ym - yi ;
    bm = yi - yj ;
    ci = -(xj-xm) ;
    cj = -(xm-xi) ;
    cm = -(xi-xj) ;
    area = abs((ai+aj+am)/2) ;
    B = [bi  0 bj  0 bm  0
          0 ci  0 cj  0 cm
         ci bi cj bj cm bm] ;
    B = B/2/area ;
return

function area = ElementArea( ie )
%  计算单元面积
%  输入参数:
%     ie ----  单元号
%  返回值:
%     area  ----  单元面积
    global gNode gElement  

    xi = gNode( gElement( ie, 1 ), 1 ) ;
    yi = gNode( gElement( ie, 1 ), 2 ) ;
    xj = gNode( gElement( ie, 2 ), 1 ) ;
    yj = gNode( gElement( ie, 2 ), 2 ) ;
    xm = gNode( gElement( ie, 3 ), 1 ) ;
    ym = gNode( gElement( ie, 3 ), 2 ) ;
    ai = xj*ym - xm*yj ;
    aj = xm*yi - xi*ym ;
    am = xi*yj - xj*yi ;
    area = abs((ai+aj+am)/2) ;
return

function k = StiffnessMatrix( ie )
%  计算单元刚度矩阵
%  输入参数:
%     ie ----  单元号
%  返回值:
%     k  ----  单元刚度矩阵

    global  gElement gMaterial 
    k = zeros( 6, 6 ) ;
    h  = gMaterial( gElement(ie, 4), 3 ) ;
    B = MatrixB( ie );
    area = ElementArea( ie );
    E  = gMaterial( gElement(ie, 4), 1 ) ;
    mu = gMaterial( gElement(ie, 4), 2 ) ;
    D = [ 1-mu    mu      0
           mu    1-mu     0
            0      0   (1-2*mu)/2] ;
    D = D*E/(1-2*mu)/(1+mu) ;
    k = transpose(B)*D*B*h*abs(area) ;    
return

function FE = elementforce( ie ,gElementStress) 
%  计算单元力矩阵
%  输入参数:
%     ie ----  单元号
%  返回值:
%     FE  ----  单元力

    global  gElement gMaterial 
    FE = zeros( 6, 6 ) ;
    thickness = gMaterial( gElement(ie, 4), 3 ) ;
    area = ElementArea(ie);
    B = MatrixB( ie );
   
    sigma0 = [gElementStress(ie,1);gElementStress(ie,2);gElementStress(ie,3)];
    FE = transpose(B)*sigma0*thickness*abs(area) ;  
return

function FEE = AssembleFE(ie, FE)
    global gElement gFE
    
    for i = 1:1:3
        gFE(gElement(ie,i)*2-1)=gFE(gElement(ie,i)*2-1)+FE(i*2-1);
        gFE(gElement(ie,i)*2)=gFE(gElement(ie,i)*2)+FE(i*2);
    end
    FEE=0;
return

function D = unlinerD (ie,eps)
%计算非线性弹性D矩阵
%  输入参数:
%     ie ----  单元号
%  返回值:
%     D  ----  D矩阵
    global  gElement gMaterial
    E  = gMaterial( gElement(ie, 4), 1 ) ;
    mu = gMaterial( gElement(ie, 4), 2 ) ;
    D = [ 1-mu    mu      0
           mu    1-mu     0
            0      0   (1-2*mu)/2] ;
    epsx = eps(1);
    epsy = eps(2);
    D = D*E*(1-100*epsx^2-100*epsy^2)/(1-2*mu)/(1+mu) ;
return

function delta = NodeDe( ie ,gDelta)
% 计算单元节点位移列阵
   global gElement  
    delta = zeros( 6, 1 ) ;   
    for j=1:1:3
        delta( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
        delta( 2*j   ) = gDelta( 2*gElement( ie, j )   ) ;
    end
return

function AssembleStiffnessMatrix( ie, k )
%  把单元刚度矩阵集成到整体刚度矩阵
%  输入参数:
%      ie  --- 单元号
%      k   --- 单元刚度矩阵
%  返回值:
%      无
    global gElement gK
    for i=1:1:3
        for j=1:1:3
            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, aa, bb, ps )
%   计算线性分布荷载的等效节点力
%   输入参数:
%      ie  -----  单元号
%      aa  -----  终点节点号
%      bb  -----  起点节点号
%      ps  -----  荷载分布集度值
%   返回值:
%      enf -----  等效节点力向量
    global gElement gNode gMaterial
    enf = zeros(6,1) ;
    h   = gMaterial( gElement( ie, 4 ), 3 ) ;
    xj  = gNode( aa, 1 ) ;
    yj  = gNode( aa, 2 ) ;
    xi  = gNode( bb, 1 ) ;
    yi  = gNode( bb, 2 ) ;
    f1 = h*ps*(yi-yj)/2 ;
    f2 = h*ps*(xj-xi)/2 ;
	f3 = h*ps*(yi-yj)/2 ;
	f4 = h*ps*(xj-xi)/2 ;
	
%   弥补前处理中的节点乱序
	switch ie
	case 446
	enf = [f1; f2; f3; f4; 0; 0] ;
	case {492, 449, 456, 458, 460, 462}
	enf = [0; 0; f1; f2; f3; f4] ;
	case 437
	enf = [f3; f4; 0; 0; f1; f2] ;
	end
return

function SaveResults( file_out )
%  显示计算结果
%  输入参数：
%     file_out  --- 保存结果文件
%  返回值：
%     无

    global gNode gElement gMaterial gBC1 gDelta gNodeStress gElementStress gK
    save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', 'gDelta', 'gNodeStress', 'gElementStress', 'gK' ) ;
return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function DisplayModel
%  用图形方式显示有限元模型
%  输入参数：
%     无
%  返回值：
%     无
    global gNode gElement gMaterial gBC1 
    
    figure ;
    axis equal ;
    axis off ;
    set( gcf, 'NumberTitle', 'off' ) ;
    set( gcf, 'Name', '有限元模型' ) ;
    
    % 根据不同的材料，显示单元颜色
    [element_number, dummy] = size( gElement ) ;
    material_color = [ 'r','g','b','c','m','y','w','k'] ;
    for i=1:element_number
        x1 = gNode( gElement( i, 1 ), 1 ) ;
        x2 = gNode( gElement( i, 2 ), 1 ) ;
        x3 = gNode( gElement( i, 3 ), 1 ) ;
				x4 = gNode( gElement( i, 4 ), 1 ) ;
        y1 = gNode( gElement( i, 1 ), 2 ) ;
        y2 = gNode( gElement( i, 2 ), 2 ) ;
        y3 = gNode( gElement( i, 3 ), 2 ) ;
				y4 = gNode( gElement( i, 4 ), 2 ) ;
        color_index = mod( 1, length( material_color ) ) ; 
        if color_index == 0 
            color_index = length( material_color ) ;
        end
        patch( [x1;x2;x3], [y1;y2;y3], material_color( color_index ) ) ;
    end
    
    % 显示边界条件
    DisplayBC( 'blue' ) ;
    
return

function DisplayBC( color )
%  用图形方式显示有限元模型的边界条件
%  输入参数：
%     color  ----  边界条件的颜色
%  返回值：
%     无
    global gNode gBC1
    
    % 确定边界条件的大小
    xmin = min( gNode(:,1) ) ;
    xmax = max( gNode(:,1) ) ;
    factor = ( xmax - xmin ) / 25 ;
    
    [bc1_number,dummy] = size( gBC1 ) ;
    dBCSize = factor ;
    for i=1:bc1_number
        if( gBC1( i, 2 ) == 1 )  %  x方向约束
            x0 = gNode( gBC1( i, 1 ), 1 ) ;
            y0 = gNode( gBC1( i, 1 ), 2 ) ;
            x1 = x0 - dBCSize ;
            y1 = y0 + dBCSize/2 ;
            x2 = x1 ;
            y2 = y0 - dBCSize/2 ;
            hLine = line( [x0 x1 x2 x0], [y0 y1 y2 y0] ) ;
            set( hLine, 'Color', color ) ;                
            
            xCenter = x1 - dBCSize/6 ;
            yCenter = y0 + dBCSize/4 ;
            radius = dBCSize/6 ;
            theta=0:pi/6:2*pi ;
            x = radius * cos( theta ) ;
            y = radius * sin( theta ) ;
            hLine = line( x+xCenter, y+yCenter ) ;
            set( hLine, 'Color', color ) ;
            
            hLine = line( x+xCenter, y+yCenter-dBCSize/2 ) ;
            set( hLine, 'Color', color ) ;          
            
            x0 = x0 - dBCSize - dBCSize/3 ;                   
            y0 = y0 + dBCSize/2 ;
            x1 = x0 ;
            y1 = y0 - dBCSize ;
            hLine = line( [x0, x1], [y0, y1] ) ;    
            set( hLine, 'Color', color ) ;          
            
            x = [x0 x0-dBCSize/6] ;
            y = [y0 y0-dBCSize/6] ;
            hLine = line( x, y ) ;
            set( hLine, 'Color', color ) ;
            for j=1:1:4
                hLine = line( x, y - dBCSize/4*j ); 
                set( hLine, 'Color', color ) ;
            end
        else                      %  y方向约束
            x0 = gNode( gBC1( i, 1 ), 1 ) ;
            y0 = gNode( gBC1( i, 1 ), 2 ) ;
            x1 = x0 - dBCSize/2 ;
            y1 = y0 - dBCSize ;
            x2 = x1 + dBCSize ;
            y2 = y1 ;
            hLine = line( [x0 x1 x2 x0], [y0 y1 y2 y0] ) ; 
            set( hLine, 'Color', color ) ;      
            
            xCenter = x0 - dBCSize/4 ;
            yCenter = y1 - dBCSize/6 ;
            radius = dBCSize/6 ;
            theta=0:pi/6:2*pi ;
            x = radius * cos( theta ) ;
            y = radius * sin( theta ) ;
            hLine = line( x+xCenter, y+yCenter ) ;
            set( hLine, 'Color', color ) ; 
            
            hLine = line( x+xCenter+dBCSize/2, y+yCenter ) ;
            set( hLine, 'Color', color ) ; 
            
            hLine = line( [x1, x1+dBCSize], [y1-dBCSize/3, y1-dBCSize/3] ) ;
            set( hLine, 'Color', color ) ; 
            
            x = [x1 x1-dBCSize/6] ;
            y = [y1-dBCSize/3 y1-dBCSize/2] ;
            hLine = line( x, y ) ;
            set( hLine, 'Color', color ) ;
            for j=1:1:4
                hLine = line( x+dBCSize/4*j, y ); 
                set( hLine, 'Color', color ) ;
            end
        end
    end
return

function PlotDisplacement( iDisp )
%  显示位移云图
%  输入参数：
%     iDisp --- 位移分量指示，它可以是下面的值
%                 1  --  x方向位移
%                 2  --  y方向位移
%  返回值：
%     无
    global gNode gElement gDelta
    
    switch iDisp
    case 1
        title = ' x 方向位移' ;
    case 2
        title = ' y 方向位移' ;
    otherwise
         fprintf( '位移分量错误\n' )  ;
        return ;
    end
    
    figure ;
    axis equal ;
    axis off ;
    set( gcf, 'NumberTitle', 'off' ) ;
    set( gcf, 'Name', title ) ;
    
    dispMin = min( gDelta( iDisp:2:length(gDelta) ) ) ;
    dispMax = max( gDelta( iDisp:2:length(gDelta) ) ) ;
    caxis( [dispMin, dispMax] ) ;
    colormap( 'jet' ) ;

    [element_number, dummy] = size( gElement ) ;
    for ie=1:1:element_number
        x = [ gNode( gElement( ie, 1 ), 1 ) ;
              gNode( gElement( ie, 2 ), 1 ) ;
							gNode( gElement( ie, 3 ), 1 ) ] ;
        y = [ gNode( gElement( ie, 1 ), 2 ) ;
              gNode( gElement( ie, 2 ), 2 ) ;
							gNode( gElement( ie, 3 ), 2 ) ] ;
        c = [ gDelta( ( gElement( ie, 1 ) - 1 ) * 2 + iDisp ) ;
              gDelta( ( gElement( ie, 2 ) - 1 ) * 2 + iDisp ) ;
							gDelta( ( gElement( ie, 3 ) - 1 ) * 2 + iDisp )] ;
        set( patch( x, y, c ), 'EdgeColor', 'interp' ) ;
    end
    
    yTick = dispMin:(dispMax-dispMin)/10:dispMax ;
    Label = cell( 1, length(yTick) ); 
    for i=1:length(yTick)
        Label{i} = sprintf( '%.2e', yTick(i) ) ;
    end
    set( colorbar( 'vert' ), 'YTick', yTick, 'YTickLabelMode', 'Manual', 'YTickLabel', Label ) ;
    
    PlotDisplacementContour( iDisp, 10, 'white' ) ;
return

function PlotDisplacementContour( iDisp, nContour, color )
%  显示位移等值线
%  输入参数：
%     iDisp ----- 位移分量指示，它可以是下面的值
%                 1  --  x方向位移
%                 2  --  y方向位移
%     nContour -- 等值线的条数
%     color  ---- 等值线颜色
%  返回值：
%     无
    global gNode gElement gDelta
   
    [element_number, dummy] = size( gElement ) ;
    [node_number, dummy] = size( gNode ) ;
    
    dispMin = min( gDelta( iDisp:2:length(gDelta) ) ) ;
    dispMax = max( gDelta( iDisp:2:length(gDelta) ) ) ;
    dispDelta = (dispMax-dispMin)/( nContour+1 ) ;
    dispContour = dispMin+dispDelta : dispDelta : dispMax - dispDelta ;
    for ie=1:1:element_number
        x = [ gNode( gElement( ie, 1 ), 1 ) ;
              gNode( gElement( ie, 2 ), 1 ) ;
							gNode( gElement( ie, 3 ), 1 ) ] ;
        y = [ gNode( gElement( ie, 1 ), 2 ) ;
              gNode( gElement( ie, 2 ), 2 ) ;
							gNode( gElement( ie, 3 ), 2 )] ;
        s = [ gDelta( ( gElement( ie, 1 ) - 1 ) * 2 + iDisp ) ;
              gDelta( ( gElement( ie, 2 ) - 1 ) * 2 + iDisp ) ;
							gDelta( ( gElement( ie, 3 ) - 1 ) * 2 + iDisp )  ] ;
        for is = 1:1:nContour  
            [smax, ismax] = max( s ) ;
            [smin, ismin] = min( s ) ;
            if dispContour(is) > smax || dispContour(is) < smin
                continue ;
            end
            
            x1 = x(ismin) + ( dispContour(is)- smin ) / (smax-smin) * ( x(ismax) - x(ismin) ) ;
            y1 = y(ismin) + ( dispContour(is)- smin ) / (smax-smin) * ( y(ismax) - y(ismin) ) ;
            
            for ismed=1:1:4
                if ismed ~= ismax && ismed ~= ismin
                    break ;
                end
            end
            
            if dispContour(is) < s( ismed )
                x2 = x(ismin) + (dispContour(is)-smin)/(s(ismed)-smin)*(x(ismed)-x(ismin)) ;
                y2 = y(ismin) + (dispContour(is)-smin)/(s(ismed)-smin)*(y(ismed)-y(ismin)) ;
            else
                x2 = x(ismed) + (dispContour(is)-s(ismed))/(smax-s(ismed))*(x(ismax)-x(ismed)) ;
                y2 = y(ismed) + (dispContour(is)-s(ismed))/(smax-s(ismed))*(y(ismax)-y(ismed)) ;
            end
          
            set( line( [x1;x2], [y1;y2] ), 'color', color ) ;
        end
    end
return