C COQ8KS    SOURCE    CHAT      05/01/12    22:22:05     5004
      SUBROUTINE COQ8KS(AKS,XX,SHPCOQ,SHPELE,NBPGAU,POIGAU,
     1 DZEGAU,ES,SIGGAU,NBNN,LRE,LRI,RI)
C==============================================================
C    MATRICE DES CONTRAINTES INITIALES  COQ8 ET COQ6
C ENTREES
C    XX(3,NBNN)  = COORDONNEES DE L ELELEMENT
C    SHPCOQ(6,NBNN,NBPGAU) = FONCTIONS DE FORME POINTS DE CONTRAINTES
C    SHPELE(6,NBNN,NBNN) = FONCTIONS DE FORME NOEUDS
C    NBPGAU              = NOMBRE DE POINTS  DE GAUSS
C    POIGAU(NBPGAU)      = POIDS D INTEGRATION
C    DZEGAU(NBPGAU)      = TROISIEME COORDONNEE REDUITE
C    ES(NBNN)            = EPAISSEURS
C    SIGGAU(5*NBPGAU)    = LES CONTRAINTES
C    NBNN                = NOMBRE DE NOEUDS
C    LRE                 = NOMBRE DE D.D.L. ( NBNN*6)
C    LRI                 = NOMBRE DE D.D.L. REDUITS ( NBNN*5)
C TABLEAU DE TRAVAIL
C    RI(LRI,LRI)         = TABLEAU DE TRAVAIL
C SORTIES
C    AKS(LRE,LRE)        = MATRICE DES CONTRAINTES INITIALES
C
C=================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(XZER=0.D0,UN=1.D0,UNDEMI=.5D0)
      DIMENSION AKS(LRE,*),XX(3,*),SHPCOQ(6,NBNN,*),SHPELE(6,NBNN,*)
      DIMENSION POIGAU(*),DZEGAU(*)
      DIMENSION SIGGAU(*),ES(*),RI(LRI,*)
      DIMENSION H(8),P(2,8),TH(8),XJ(3,3),XJI(3,3)
      DIMENSION XJIT(3,3),TET(3,3),TT(9),TE(3,3),SIG(3,3),SIGT(3,3)
      DIMENSION SI(3,3),A(3,8),V2V1(3,2),V2V1T(2,3),RIP(2,2)
      DIMENSION VR(3,2),AKSP(3,3),TXR(3,3,8)
      DIMENSION AMA1(3,3),AMA2(3,3),AMA3(3,3)
C
C
      EQUIVALENCE(TET(1,1),TT(1))
C
C        OBJET  MATRICE DE RIGIDITE GEOMETRIQUE ELEMENTAIRE
C
      DO 5 I=1,NBNN
  5   TH(I)=ES(I)
      DO 10 I=1,LRE
      DO 10 J=1,LRE
 10   AKS(I,J)=XZER
      CALL CQ8LOC(XX,NBNN,SHPELE,TXR,IRR)
      DO 50 NGAU=1,NBPGAU
      DO 20 I=1,NBNN
      H(I)  =SHPCOQ(1,I,NGAU)
      P(1,I)=SHPCOQ(2,I,NGAU)
      P(2,I)=SHPCOQ(3,I,NGAU)
 20   CONTINUE
      DO 11 I=1,LRI
      DO 11 J=1,LRI
      RI(I,J)=XZER
 11   CONTINUE
      E3=DZEGAU(NGAU)
      WT=POIGAU(NGAU)
      CALL COQ8JC(NGAU,NBNN,E3,XX,TH,TXR,SHPCOQ,XJ,DET,IRR)
      FACT=WT*DET
C
C       INVERSION DE LA MATRICE JACOBIENNE
C
      DUM= UN/DET
      XJI(1,1)=DUM*( XJ(2,2)*XJ(3,3)-XJ(2,3)*XJ(3,2))
      XJI(2,1)=DUM*(-XJ(2,1)*XJ(3,3)+XJ(2,3)*XJ(3,1))
      XJI(3,1)=DUM*( XJ(2,1)*XJ(3,2)-XJ(2,2)*XJ(3,1))
      XJI(1,2)=DUM*(-XJ(1,2)*XJ(3,3)+XJ(1,3)*XJ(3,2))
      XJI(2,2)=DUM*( XJ(1,1)*XJ(3,3)-XJ(1,3)*XJ(3,1))
      XJI(3,2)=DUM*(-XJ(1,1)*XJ(3,2)+XJ(1,2)*XJ(3,1))
      XJI(1,3)=DUM*( XJ(1,2)*XJ(2,3)-XJ(1,3)*XJ(2,2))
      XJI(2,3)=DUM*(-XJ(1,1)*XJ(2,3)+XJ(1,3)*XJ(2,1))
      XJI(3,3)=DUM*( XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1))
C
C       CALCUL DE J-1 TRANSPOSEE
C
      DO 30 I=1,3
      DO 30 J=1,3
      XJIT(I,J)=XJI(J,I)
 30   CONTINUE
C
C  CALCUL DE TET MATRICE DE LA ROTATION LOCALE AU POINT DE GAUSS
C
      DO 31 I=1,3
      DO 31 J=1,2
      TET(I,J)=XJ(J,I)
 31   CONTINUE
      CALL CROSS2(TT(1),TT(4),TT(7),IERR)
      CALL CROSS2(TT(7),TT(1),TT(4),IERR)
      CALL CROSS2(TT(4),TT(7),TT(1),IERR)
C
C       TE EST LA TRANSPOSEE DE TET
C
      DO 32 I=1,3
      DO 32 J=1,3
      TE(I,J)=TET(J,I)
 32    CONTINUE
C
C       ON RECUPERE LES CONTRAINTES D UN POINT DE GAUSS
C
      SIG(1,1)=SIGGAU((NGAU-1)*5+1)
      SIG(2,2)=SIGGAU((NGAU-1)*5+2)
      SIG(1,2)=SIGGAU((NGAU-1)*5+3)
      SIG(1,3)=SIGGAU((NGAU-1)*5+4)
      SIG(2,3)=SIGGAU((NGAU-1)*5+5)
      SIG(2,1)=SIG(1,2)
      SIG(3,2)=SIG(2,3)
      SIG(3,1)=SIG(1,3)
      SIG(3,3)=XZER
C
C       CALCUL DE SIGT=TET*SIG*TE
C
      CALL MULMAT(AMA1,XJIT,TET,3,3,3)
      CALL MULMAT(AMA2,AMA1,SIG,3,3,3)
      CALL MULMAT(AMA3,AMA2,TE,3,3,3)
      CALL MULMAT(SIGT,AMA3,XJI,3,3,3)
C
C       CALCUL DES TERMES RELATIFS AU NOEUD I,J
C
      DO 40 I=1,NBNN
      DO 100 J=I,NBNN
      A(1,I)=P(1,I)*E3*TH(I)*UNDEMI
      A(2,I)=P(2,I)*E3*TH(I)*UNDEMI
      A(1,J)=P(1,J)*E3*TH(J)*UNDEMI
      A(2,J)=P(2,J)*E3*TH(J)*UNDEMI
      A(3,I)=H(I)*TH(I)*UNDEMI
      A(3,J)=H(J)*TH(J)*UNDEMI
C
      DO 116 K=1,3
      IK=5*(I-1)+K
      JK=5*(J-1)+K
      DO 200 M=1,2
      DO 200 N=1,2
      RI(IK,JK)=RI(IK,JK)+SIGT(M,N)*P(M,I)*P(N,J)
 200  CONTINUE
116   CONTINUE
C
      DO 117 K=1,3
      IS=-1
      DO 117 L=4,5
      LL=L-3
C
C
      IK=5*(I-1)+K
      JL=5*(J-1)+L
      IS=-IS
      DO 210 M=1,2
      DO 210 N=1,3
      RI(IK,JL)=RI(IK,JL)+SIGT(M,N)*P(M,I)*TXR(K,LL,J)*A(N,J)*IS
 210  CONTINUE
117   CONTINUE
      IS=-1
C
      DO 119 K=4,5
      IS=-IS
      DO 119 L=1,3
      IK=5*(I-1)+K
      JL=5*(J-1)+L
      KK=K-3
      DO 230 N=1,2
      DO 230 M=1,3
      RI(IK,JL)=RI(IK,JL)+SIGT(M,N)*P(N,J)*TXR(L,KK,I)*A(M,I)*IS
 230  CONTINUE
119   CONTINUE
      IS=1
C
      DO 121 K=4,5
      IS=-IS
      DO 121 L=4,5
      IS=-IS
      KK=K-3
      LL=L-3
      IK=5*(I-1)+K
      JL=5*(J-1)+L
      COE1=XZER
      COE2=XZER
      DO 240 M=1,3
      DO 240 N=1,3
      COE1=COE1+SIGT(M,N)*A(M,I)*A(N,J)
 240  CONTINUE
      DO 241 MM=1,3
      COE2=COE2+TXR(MM,KK,I)*TXR(MM,LL,J)
 241  CONTINUE
       RI(IK,JL)=COE1*COE2*IS
121   CONTINUE
C
C       CHANGEMENT DE VARIABLE ROTATIONS LOCALES-GLOBALES
C                  ET INTEGRATION
C
      DO 130 K=1,3
      DO 130 L=1,3
      IK6=6*(I-1)+K
      JL6=6*(J-1)+L
      IK5=5*(I-1)+K
      JL5=5*(J-1)+L
      AKS(IK6,JL6)=AKS(IK6,JL6)+FACT*RI(IK5,JL5)
130   CONTINUE
C
100   CONTINUE
 40   CONTINUE
 50   CONTINUE
C
      DO 150 I=1,LRE
      DO 150 J=1,I
      AKS(I,J)=AKS(J,I)
150   CONTINUE
C
      RETURN
      END