function es = ElementStress( ie )
%  计算单元的应力分量
%  输入参数
%      ie  ----- 单元号
%  返回值
%      es ----- 单元应力分量列阵（1×6）： [sx, sy, txy, s1, s2, tmax]
    global gElement gDelta  gMaterial
    
    es = zeros( 1, 6 ) ;   % 单元的应力分量
    de = zeros( 6, 1 ) ;   % 单元节点位移列阵
    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) ;
    B = MatrixB( ie ) ;
    for j=1:1:3
        de( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
        de( 2*j   ) = gDelta( 2*gElement( ie, j )   ) ;
    end
    es(1:3) = D * B * de ;
    es(6) = 0.5*sqrt((es(1)-es(2))^2 + 4*es(3)^2 ) ;
    es(4) = 0.5*(es(1)+es(2)) + es(6) ;
    es(5) = 0.5*(es(1)+es(2)) - es(6) ;
return