C COQ3GY    SOURCE    CHAT      05/01/12    22:21:10     5004           
      SUBROUTINE COQ3GY(XE,RHO,XEL,BPSS,BB,REL,VROT)
c=======================================================================
c
c       calcule la matrice de couplage gyroscopique des elements de coque
c       a 3 noeud
c
c   entree
c      xe(3,3)=coodonnees de l element
c      rho    =masse volumique * epaisseur
c      vrot (3) =  vecteur rotation
c   travail
c      xel(3,3)   =coordonnees locales de l element
c      bpss(3,3)  =matrice de passage
c      bb(9)      =stocke les fonctions de formes de flexion
c   sorties
c      rewo(18,18)=matrice dans le repere local puis global
c
c      Didier COMBESCURE mars 2003
c
c=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XE(3,*),XEL(3,*),BPSS(3,*),REL(18,*),BB(*),VROT(*)
      DIMENSION REWOLP(18,18),REWO(18,18),ROME(6,6)
      DIMENSION XX(3),YY(3)
      DATA XX/0.5D0,0.5D0,0.0D0/
      DATA YY/0.0D0,0.5D0,0.5D0/
      DATA UNSIX,UN12/.166666666666666666D0,.833333333333333333D-1/
      DATA UNTIER/.33333333333333333D0/
      DATA IZERO,IUN/0,1/
C
C
      DO 10 IN=1,6
        DO 10 IM=1,6
         ROME(IN,IM) = 0.D0
   10 CONTINUE
C
      ROME(1,1) = 0.
      ROME(2,2) = 0.
      ROME(3,3) = 0.
      ROME(1,2) = (-1.)*VROT(3)
      ROME(1,3) = VROT(2)
      ROME(2,3) = (-1.)*VROT(1)
      ROME(2,1) = VROT(3)
      ROME(3,1) = (-1.)*VROT(2)
      ROME(3,2) = VROT(1)
C
c     matrice de passage
c
      CALL VPAST(XE,BPSS)
c
c     coordonnees locales
c
      CALL VCORLC(XE,XEL,BPSS)
c
c     mise a 0 de la matrice de masse
c
      CALL ZERO(REWO,18,18)
      X21=XEL(1,2)-XEL(1,1)
      Y31=XEL(2,3)-XEL(2,1)
      SURF=X21*Y31*.5D0*RHO
c
c     termes diagonaux  de la membrane
c
      DO 100 IA=1,13,6
          REWO(IA  ,IA  )=UNSIX
          REWO(IA+1,IA+1)=UNSIX
          DO 110 IB=1,13,6
              IF(IB.EQ.IA) GOTO 110
              REWO(IA  ,IB  )=UN12
              REWO(IA+1,IB+1)=UN12
  110     CONTINUE
  100 CONTINUE
c
c     partie flexion   integration numerique
c
c     ia numero du point de gauss
c     ib numero d un noeud
c     ic numero d un noeud
c
      DO 200 IA=1,3
          CALL MFDKT(XX(IA),YY(IA),XEL,BB)
          DO 210 IB=1,3
              IBDERI=(IB-1)*6+2
              IBDEBB=(IB-1)*3
              DO 210 IC=1,3
                  ICDERI=(IC-1)*6+2
                  ICDEBB=(IC-1)*3
                  DO 220 ID=1,3
                   DO 220 IE=1,3
                    REWO(IBDERI+ID,ICDERI+IE)=REWO(IBDERI+ID,ICDERI+IE)
     1               +BB(IBDEBB+ID)*BB(ICDEBB+IE)*UNTIER
  220             CONTINUE
  210     CONTINUE
  200 CONTINUE
c
c     multiplication par rho*surf
c
      DO 300 IA=1,18
          DO 300 IB=1,18
              REWO(IA,IB)=REWO(IA,IB)*SURF
  300 CONTINUE
c   
C       WRITE(*,*) 'REWO avant'    
       DO 517 I=1,18
C         WRITE(*,*)  REWO(I,I)
  517  CONTINUE       
c
c     changement de repere
c
      CALL TRANSK(REWO,BPSS,18,3,IUN)
C
C On bouscule la matrice de masse
C
      CALL ZERO(REL,18,18)
C
      DO   401 IN=1,3
C
        DO   402 IM=1,3
C
          DO   403 IDD1=1,3
C
            DO   404 IDD2=1,3
C
            REL((6*IN)-6+IDD1,(6*IM)-6+IDD2)=
     . REWO((6*IN)-6+IDD2,(6*IM)-6+IDD2)*ROME(IDD1,IDD2) 
            REL((6*IN)-3+IDD1,(6*IM)-6+IDD2)=
     . REWO((6*IN)-3+IDD2,(6*IM)-6+IDD2)*ROME(IDD1,IDD2) 
            REL((6*IN)-6+IDD1,(6*IM)-3+IDD2)=
     . REWO((6*IN)-6+IDD2,(6*IM)-3+IDD2)*ROME(IDD1,IDD2) 
            REL((6*IN)-3+IDD1,(6*IM)-3+IDD2)=
     . REWO((6*IN)-3+IDD2,(6*IM)-3+IDD2)*ROME(IDD1,IDD2) 
C 
  404      CONTINUE    
C 
  403     CONTINUE    
C 
  402   CONTINUE                                                         
C 
  401  CONTINUE   
C 
C       WRITE(*,*) 'REWO'    
       DO 518 I=1,6
C         WRITE(*,*)  REWO(I,I)
  518  CONTINUE
C       WRITE(*,*) 'REL'    
       DO 519 I=1,3
       DO 519 J=1,3
C         WRITE(*,*)  REL((I + 3) ,(J + 3) )
C         WRITE(*,*)  REL((J + 3) ,(I + 3) )
  519  CONTINUE
c
      RETURN
      END