      SUBROUTINE etc5noe(coor,car,iopt,U,SIGMA)
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C  BUT : CALCUL DU TABLEAU DES CONTRAINTES ELEMENTAIRES AU CENTRE DE
C  ---   L ELEMENT PAR LA FORMULE :
C             SIGMA(I)=2.*E1*(D(I)-.5*TR(D)*ID)+2.*RR*TR(D)*ID
C        POUR L ELEMENT QUADRILATERAL A 4 SOUS TRIANGLES EN CROIX
C              ---   VERSION PLANE  ---
C  N' EST ECRIT CI DESSOUS QUE LA PARTIE DE CALCUL DES CONTRAINTES
C  ELASTIQUES POUR UN MATERIAU HOMOGENE POUR L'OPTION NTHELA=0
C  MANQUE: (1)   LE TRAITEMENT DES OPTIONS NTHELA DIFFERENTES DE ZERO
C          (2)   LE CALCUL POUR DES MATERIAUX NON HOMOGENES
C
c in : coor(noe,ndim) : coordones des 4 sommets.
c      car(2): caracteristiques des materiaux
c              car(1) = young
c              car(2) = poisson
c      iopt = 1 isotrope Deformations Planes
c           = 2 isotrope Contraintes  Planes
C out: sigma(3)  CONTRAINTES SXX,SYY,SXY
C ....................................................................
      implicit none
      integer N(3,4), i, i1, np, j, k, ii, i2, iopt
      DOUBLE PRECISION coor(5,2),car(2),U(2,5),young,poison,DIV1,RR,
     +                 E1,X(5),Y(5),BI1(3),BI2(3),D11(10),D22(10),
     +                 D12(10),DIV(10),S,B11,B12,X31,Y31,X42,Y42,
     +                 DSIGMA(3,10),sigma(3)
      DATA N/1,2,5,2,3,5,3,4,5,4,1,5/
C
      YOUNG  = car(1)
      POISON = car(2)
      E1 = YOUNG*.5D0/(1.D0+POISON)
      IF(iopt .eq. 1) then
         RR=E1*.5/(1.-2.*poison)
      ELSEIF(iopt .eq. 2 )then
         RR=E1*.5*(1.+poison)/(1.-poison)
      ELSE
         call mexErrMsgTxt('CHOIX CONTRAINTES OU DEFORMATIONS PLANES ?')
         RETURN
      ENDIF
C
C     CALCUL DES COORDONNEES DES NOEUDS
C
      DO 1 I = 1,4
        X(I) = coor(I,1)
        Y(I) = coor(I,2)
 1    continue
      X31=X(3)-X(1)
      Y31=Y(3)-Y(1)
      X42=X(4)-X(2)
      Y42=Y(4)-Y(2)
      DIV1=1./(Y31*X42-Y42*X31)
      X(5)=(X(1)*Y31*X42-X(2)*Y42*X31+(Y(2)-Y(1))*X31*X42)*DIV1
      Y(5)=(Y(2)*Y31*X42-Y(1)*Y42*X31-(X(2)-X(1))*Y31*Y42)*DIV1
C
C     INITIALISATION
C
      DIV1=0.
      DO 100 I1=1,10
         DIV(I1)=0.
         D11(I1)=0.
         D12(I1)=0.
         D22(I1)=0.
 100  continue
C
C     BOUCLE SUR LES SOUS-ELEMENTS
C
      DO 10 NP=1,4
         I=N(1,NP)
         J=N(2,NP)
         K=N(3,NP)
C
C        CALCUL DES GRADIENTS BIJ DES COORDONNEES BARYCENTRIQUES
C        ET DE S=2*AIRE DU SOUS ELEMENT
C
         S   = (X(I)-X(J))*(Y(I)-Y(K))-(X(I)-X(K))*(Y(I)-Y(J))
         BI1(1) = (Y(J)-Y(K))
         BI1(2) = (Y(K)-Y(I))
         BI1(3) = (Y(I)-Y(J))
         BI2(1) = (X(K)-X(J))
         BI2(2) = (X(I)-X(K))
         BI2(3) = (X(J)-X(I))
C
C        ACTUALISATION DES DEFORMATIONS
         DIV1 = DIV1 + S
         DO 101 II=1,3
            I1=2*N(II,NP)-1
            I2=I1+1
            B11=BI1(II)*.5
            B12=BI2(II)*.5
            D11(I1)=D11(I1)+B11
            D12(I1)=B12+D12(I1)
            D22(I1)=-B11+D22(I1)
            DIV(I1)=B11*2.+DIV(I1)
            D11(I2)=-B12+D11(I2)
            D12(I2)=B11+D12(I2)
            D22(I2)=B12+D22(I2)
            DIV(I2)=DIV(I2)+B12*2.
  101    CONTINUE
   10 CONTINUE
C
C     CALCUL DE SIGMA AU CENTRE (DANS L'ORDRE SXX,SYY,SXY)
C     SIGMA XX => DSIGMA(1,I)
C     SIGMA YY => DSIGMA(2,I)
C     SIGMA XY => DSIGMA(3,I)
C
      DIV1 = 2./DIV1
      DO 102 J=1,10
         DSIGMA(1,J)=(E1*(D11(J))+RR*DIV(J))*DIV1
         DSIGMA(2,J)=E1*D12(J)*DIV1
         DSIGMA(3,J)=(E1*(D22(J))+RR*DIV(J))*DIV1
  102 CONTINUE
C     [      11]   [1 4 7 10 13 16 19 22 25 28]   [          ]
C     [SIGMA 22] = [2 5 8 11 14 17 20 23 26 29] * [u_solution]
C     [      12]   [3 6 9 12 15 18 21 24 27 30]   [          ]
C             3*1                3*9                     10*1
C
      sigma(1) = DSIGMA(1,1)*U(1,1) + DSIGMA(1, 2)*U(2,1)
     +         + DSIGMA(1,3)*U(1,2) + DSIGMA(1, 4)*U(2,2)
     +         + DSIGMA(1,5)*U(1,3) + DSIGMA(1, 6)*U(2,3)
     +         + DSIGMA(1,7)*U(1,4) + DSIGMA(1, 8)*U(2,4)
     +         + DSIGMA(1,9)*U(1,5) + DSIGMA(1,10)*U(2,5)
      sigma(2) = DSIGMA(2,1)*U(1,1) + DSIGMA(2, 2)*U(2,1)
     +         + DSIGMA(2,3)*U(1,2) + DSIGMA(2, 4)*U(2,2)
     +         + DSIGMA(2,5)*U(1,3) + DSIGMA(2, 6)*U(2,3)
     +         + DSIGMA(2,7)*U(1,4) + DSIGMA(2, 8)*U(2,4)
     +         + DSIGMA(2,9)*U(1,5) + DSIGMA(2,10)*U(2,5)
      sigma(3) = DSIGMA(3,1)*U(1,1) + DSIGMA(3, 2)*U(2,1)
     +         + DSIGMA(3,3)*U(1,2) + DSIGMA(3, 4)*U(2,2)
     +         + DSIGMA(3,5)*U(1,3) + DSIGMA(3, 6)*U(2,3)
     +         + DSIGMA(3,7)*U(1,4) + DSIGMA(3, 8)*U(2,4)
     +         + DSIGMA(3,9)*U(1,5) + DSIGMA(3,10)*U(2,5)
C
       END
