      SUBROUTINE ets2q1d(coor,fomega,fgamma,pressi,norefs,alpha,
     +                   theta,car,iopt,be)
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C BUT: CALCUL DE SECONDS MEMBRES DE L ELEMENT QUAD 2Q1D
C ---  
c in : 
c   coor(noe,ndim)        : coordones des 4 sommets.
c   fomega(ndim,noe)      : fx, fy aux noeuds
c   fgamma(ndim,2*nbarete): fx, fy aux noeuds de chaque arete.
c   pressi(2*nbarete)     : pression aux noeuds de chaque arete
c   norefs(nbarete,2): norefs(i,1) = 0 si fgamma   = 0 sur arete_i
c                      norefs(i,2) = 0 si pression = 0 sur arete_i
C   theta(4): temperature aux 4 sommets
C   alpha(3): tenseur de dilatation thermique (xx, xy, yy)
C
c   car            : caracteristiques des materiaux
c   iopt = 1 isotrope Deformations Planes
c        = 2 isotrope Contraintes  Planes
c        = sinon anisotrope
c   car(6): caracteristiques des materiaux
c           if(iopt .eq. 1 .or. iopt.eq. 2) then
c             car(1) = young
c             car(2) = poisson
c           else
c             car: E11, E12, E22, E13, E23, E33 avec
c 
C                     E11   E12   E13       
C                           E22   E23
C                                 E33           
c            end if
c
C out: BE(8): second membre.
C  ....................................................................
      IMPLICIT NONE
      integer I , J , IJ, iopt,norefs(4,2)
      double precision coor(4,2),car(6),fomega(2,4),fgamma(2,8),be(8),
     +                 pressi(8),theta(4),ALPHA(3),X21,Y21,X41,Y41,X32,
     +                 Y32,X42,Y42,X43,Y43,X31,Y31,DELT1,DELT2,DELT3,
     +                 DELT4,DEL(4),DELTA,E(6),A(3),YOUNG,
     +              POISSON,C,UNMNU,xji,yji,xnu,ynu,arelon
C
      x21  = coor(2,1) - coor(1,1)
      y21  = coor(2,2) - coor(1,2)
      x31  = coor(3,1) - coor(1,1)
      y31  = coor(3,2) - coor(1,2)
      x32  = coor(3,1) - coor(2,1)
      y32  = coor(3,2) - coor(2,2)
      x41  = coor(4,1) - coor(1,1)
      y41  = coor(4,2) - coor(1,2)
      x42  = coor(4,1) - coor(2,1)
      y42  = coor(4,2) - coor(2,2)
      x43  = coor(4,1) - coor(3,1)
      y43  = coor(4,2) - coor(3,2)
C
      DELT1  = X21*Y41 - X41*Y21
      DELT2  = X21*Y32 - X32*Y21
      DELT3  = X32*Y43 - X43*Y32
      DELT4  = X41*Y43 - X43*Y41
C
      DEL(1)  = DELT1 * .25D+00
      DEL(2)  = DELT2 * .25D+00
      DEL(3)  = DELT3 * .25D+00
      DEL(4)  = DELT4 * .25D+00
C
C
C --  FOMEGA: INTEGRATION AUX 4 SOMMETS (F(0,0)+F(1,0)+F(0,1)+F(1,1)
C
      DO 1 I=1,2
         DO 1 J=1,4
            IJ = 2 * ( J - 1 ) + I
            BE(IJ) = FOMEGA(I,J) * DEL(J)
 1    CONTINUE
C
C --  FGAMMA: INTEGRATION AUX 2 EXTREMITES DES COTES (F(0)+F(1))/2
C
      if(norefs(1,2) .ne. 0) then
        XJI    = coor(2,1) - coor(1,1)
        YJI    = coor(2,2) - coor(1,2)
        ARELON = DSQRT(XJI**2 + YJI**2)
        XNU    =  YJI/ARELON
        YNU    = -XJI/ARELON
      end if
      DELTA = DSQRT(X21*X21 + Y21*Y21) / 2.D+00
      BE(1) = BE(1) + DELTA*(FGAMMA(1,1) - pressi(1) * XNU) 
      BE(3) = BE(3) + DELTA*(FGAMMA(1,2) - pressi(2) * XNU)
      BE(2) = BE(2) + DELTA*(FGAMMA(2,1) - pressi(1) * YNU)
      BE(4) = BE(4) + DELTA*(FGAMMA(2,2) - pressi(2) * YNU)
C
      if(norefs(2,2) .ne. 0) then
        XJI    = coor(3,1) - coor(2,1)
        YJI    = coor(3,2) - coor(2,2)
        ARELON = DSQRT(XJI**2 + YJI**2)
        XNU    =  YJI/ARELON
        YNU    = -XJI/ARELON
      end if
      DELTA = DSQRT(X32*X32 + Y32*Y32) / 2.D+00
      BE(3) = BE(3) + DELTA*(FGAMMA(1,3) - pressi(3) * XNU)
      BE(5) = BE(5) + DELTA*(FGAMMA(1,4) - pressi(4) * XNU)
      BE(4) = BE(4) + DELTA*(FGAMMA(2,3) - pressi(3) * YNU)
      BE(6) = BE(6) + DELTA*(FGAMMA(2,4) - pressi(4) * YNU)
C
      if(norefs(3,2) .ne. 0) then
        XJI    = coor(4,1) - coor(3,1)
        YJI    = coor(4,2) - coor(3,2)
        ARELON = DSQRT(XJI**2 + YJI**2)
        XNU    =  YJI/ARELON
        YNU    = -XJI/ARELON
      end if
      DELTA = DSQRT(X43*X43 + Y43*Y43) / 2.D+00
      BE(5) = BE(5) + DELTA*(FGAMMA(1,5) - pressi(5) * XNU)
      BE(7) = BE(7) + DELTA*(FGAMMA(1,6) - pressi(6) * XNU)
      BE(6) = BE(6) + DELTA*(FGAMMA(2,5) - pressi(5) * YNU)
      BE(8) = BE(8) + DELTA*(FGAMMA(2,6) - pressi(6) * YNU)
C
      if(norefs(4,2) .ne. 0) then
        XJI    = coor(1,1) - coor(4,1)
        YJI    = coor(1,2) - coor(4,2)
        ARELON = DSQRT(XJI**2 + YJI**2)
        XNU    =  YJI/ARELON
        YNU    = -XJI/ARELON
      end if
      DELTA = DSQRT(X41*X41 + Y41*Y41) / 2.D+00
      BE(1) = BE(1) + DELTA*(FGAMMA(1,8) - pressi(8) * XNU)
      BE(7) = BE(7) + DELTA*(FGAMMA(1,7) - pressi(7) * XNU)
      BE(2) = BE(2) + DELTA*(FGAMMA(2,8) - pressi(8) * YNU)
      BE(8) = BE(8) + DELTA*(FGAMMA(2,7) - pressi(7) * YNU)
 2    CONTINUE
C
C     CONTRAINTES THERMIQUES = HOOKE * ALPHA * THETA
C                         3,1     3,3   3,3       3,1
      if(iopt .eq. 1) then
C  --    CONTRAINTES PLANES     (ISOTROPE)     -----
         YOUNG   = car(1)
         POISSON = car(2)
         C   = YOUNG / (1.D+00 - POISSON*POISSON)
         E(1)  = C
         E(2)  = C*POISSON
         E(3)  = C
         E(4)  = 0.D0
         E(5)  = 0.D0
         E(6)  = C * (1.D+00 - POISSON) / 2.D+00
      else if(iopt .eq. 2) then
C  --    DEFORMATIONS PLANES (ISOTROPE)     -----
         YOUNG   = car(1)
         POISSON = car(2)
         UNMNU   = 1.D+00 - POISSON
         C   = YOUNG*UNMNU
         C   = C /( (1.D+00 + POISSON) * (1.D+00 - 2.D+00*POISSON) )
         E(1)  = C
         E(2)  = POISSON*C / UNMNU
         E(3)  = C
         E(4)  = 0.D0
         E(5)  = 0.D0
         E(6)  = (C * (1.D+00 - 2.D+00*POISSON) ) / (2.D+00*UNMNU)
      else
C  --    CAS ANISOTROPE     -----
        do 3 i = 1 , 6
           E(i) = car(i)
 3      continue
      end if
C
C     TENSEUR D'ELASTICITE * TENSEUR THERMIQUE
C
      A(1) = ALPHA(1)*E(1) + ALPHA(3)*E(2)  + ALPHA(2)*E(4)
      A(2) = ALPHA(1)*E(2) + ALPHA(3)*E(3)  + ALPHA(2)*E(5)
      A(3) = ALPHA(1)*E(4) + ALPHA(3)*E(5)  + ALPHA(2)*E(6)
C
      BE(1) = BE(1)
     +     + ( X42 * A(3)*THETA(1) - Y42 * A(1)*THETA(1) ) / 4.
     +     + ( X32 * A(3)*THETA(2) - Y32 * A(1)*THETA(2) ) / 4.
     +     + ( X43 * A(3)*THETA(4) - Y43 * A(1)*THETA(4) ) / 4.
      BE(2) = BE(2)
     +     + ( X42 * A(2)*THETA(1) - Y42 * A(3)*THETA(1) ) / 4.
     +     + ( X32 * A(2)*THETA(2) - Y32 * A(3)*THETA(2) ) / 4.
     +     + ( X43 * A(2)*THETA(4) - Y43 * A(3)*THETA(4) ) / 4.
      BE(3) = BE(3)
     +     + ( - X41 * A(3)*THETA(1) + Y41 * A(1)*THETA(1) ) / 4.
     +     + ( - X31 * A(3)*THETA(2) + Y31 * A(1)*THETA(2) ) / 4.
     +     + (   X43 * A(3)*THETA(3) - Y43 * A(1)*THETA(3) ) / 4.
      BE(4) = BE(4)
     +     + ( - X41 * A(2)*THETA(1) + Y41 * A(3)*THETA(1) ) / 4.
     +     + ( - X31 * A(2)*THETA(2) + Y31 * A(3)*THETA(2) ) / 4.
     +     + (   X43 * A(2)*THETA(3) - Y43 * A(3)*THETA(3) ) / 4.
      BE(5) = BE(5)
     +     + ( X21 * A(3)*THETA(2) - Y21 * A(1)*THETA(2) ) / 4.
     +     + ( - X42 * A(3)*THETA(3) + Y42 * A(1)*THETA(3) ) / 4.
     +     + ( - X41 * A(3)*THETA(4) + Y41 * A(1)*THETA(4) ) / 4.
      BE(6) = BE(6)
     +     + ( X21 * A(2)*THETA(2) - Y21 * A(3)*THETA(2) ) / 4.
     +     + ( - X42 * A(2)*THETA(3) + Y42 * A(3)*THETA(3) ) / 4.
     +     + ( - X41 * A(2)*THETA(4) + Y41 * A(3)*THETA(4) ) / 4.
      BE(7) = BE(7)
     +     + ( X21 * A(3)*THETA(1) - Y21 * A(1)*THETA(1) ) / 4.
     +     + ( X32 * A(3)*THETA(3) - Y32 * A(1)*THETA(3) ) / 4.
     +     + ( X31 * A(3)*THETA(4) - Y31 * A(1)*THETA(4) ) / 4.
      BE(8) = BE(8)
     +     + ( X21 * A(2)*THETA(1) - Y21 * A(3)*THETA(1) ) / 4.
     +     + ( X32 * A(2)*THETA(3) - Y32 * A(3)*THETA(3) ) / 4.
     +     + ( X31 * A(2)*THETA(4) - Y31 * A(3)*THETA(4) ) / 4.
      END
