C COQ3KC    SOURCE    KK2000    05/05/25    21:15:12     5088           
      SUBROUTINE COQ3KC(XE,RHO,XEL,BPSS,BB,REL,VROT)
c=======================================================================
c
c      calcule la matrice de raideur centrifuge 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      rel(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,*),REL(18,*),VROT(3)
      DIMENSION REWOLP(18,18),REWO(18,18),ROME(6,6)
      DIMENSION XEL(3,3),BPSS(3,3),BB(3)
      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/0/
C
c
      DO 10 IN=1,6
        DO 10 IM=1,6
         ROME(IN,IM) = 0.D0
   10 CONTINUE
c
      ROME(1,1) = (-1.)*((VROT(2)**2) + (VROT(3)**2))
      ROME(2,2) = (-1.)*((VROT(1)**2) + (VROT(3)**2))
      ROME(3,3) = (-1.)*((VROT(1)**2) + (VROT(2)**2))
      ROME(1,2) = VROT(1)*VROT(2)
      ROME(1,3) = VROT(1)*VROT(3)
      ROME(2,3) = VROT(2)*VROT(3)
      ROME(2,1) = ROME(1,2) 
      ROME(3,1) = ROME(1,3) 
      ROME(3,2) = ROME(2,3) 
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
c     changement de repere
c
      CALL TRANSK(REWO,BPSS,18,3,IZERO)
C
C On bouscule la matrice de raideur
C
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)-6+IDD2,(6*IM)-6+IDD2)*ROME(IDD1,IDD2) 
            REL((6*IN)-6+IDD1,(6*IM)-3+IDD2)=
     . REWO((6*IN)-6+IDD2,(6*IM)-6+IDD2)*ROME(IDD1,IDD2) 
            REL((6*IN)-3+IDD1,(6*IM)-3+IDD2)=
     . REWO((6*IN)-6+IDD2,(6*IM)-6+IDD2)*ROME(IDD1,IDD2) 
C 
  404      CONTINUE    
C 
  403     CONTINUE    
C 
  402   CONTINUE                                                         
C 
  401  CONTINUE   
C 
c
      RETURN
      END