Sources Esope | Cast3M

kpcoq4
C KPCOQ4    SOURCE    PV        11/03/07    21:17:18     6885
      SUBROUTINE KPCOQ4(XX, XP, XKP, IANT)
C
C     Procedure de calcul de la matrice Kppour un element COQ4
C     Entrees : XX(4, 3) : REAL*8 : Coordonnees des noeuds
C               XP : REAL : Pression
C               IANT : INTEGER : 1 si calcul asymétrique, 0 sinon
C     Sortie :  XKP(24, 24) : REAL*8 : Matrice Kp elementaire
C
C     (D'apres "Design variations of nonlinear elastic structures
C                subjected to follower forces"
C                 M.J. Poldneff, I.S. Rai, J.S. Arora)
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XX(3, 4), XKP(24, 24), SKRO(3, 3, 3),
     1    XN(4, 9), XNB1(9), XNB2(9), XNB3(9), XNB4(9), XNB(9),
     2    DNB1(9), DNB2(9), DN(2, 4, 9), DX(2, 3, 9), XDUM(24, 24),
     3    XTIN(9), A(9), B(9), C(9), D(9), E(9), F(9)
      DATA XNB1/0.5D0,0.D0,-0.5D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0/
      DATA XNB2/0.5D0,0.D0,0.5D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0/
      DATA XNB3/0.5D0,-0.5D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0/
      DATA XNB4/0.5D0,0.5,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0/
C
C     Initialisation du symbole epsilon
C
      DO 10 I = 1, 3
         DO 10 J = 1, 3
            DO 10 K = 1, 3
   10          SKRO(I, J, K) = 0.D0
      SKRO(1, 2, 3) = 1.D0
      SKRO(1, 3, 2) = -1.D0
      SKRO(2, 3, 1) = 1.D0
      SKRO(2, 1, 3) = -1.D0
      SKRO(3, 1, 2) = 1.D0
      SKRO(3, 2, 1) = -1.D0
C
C     Calcul de la pression
C
      P = XP
*      DO 20 I = 1, 4
*   20    P = P + XP(I)
*      P = P/4.
*      WRITE (*,*) 'Pression : ', P
C
C     Fonctions de forme et derivees
C
C     Les coefficients sont ranges comme suit :
C        indice : 1 2  3    4    5    6      7       8        9
C        terme :  1 T2 T1 T1*T2 T2^2 T1^2 T1*T2^2 T1^2*T2 T1^2*T2^2
C
      CALL MULTP2(XNB1, XNB3, XNB)
      CALL DERIP2(XNB, 1, DNB1)
      CALL DERIP2(XNB, 2, DNB2)
      DO 31 I = 1, 9
         XN(1, I) = XNB(I)
         DN(1, 1, I) = DNB1(I)
         DN(2, 1, I) = DNB2(I)
   31 CONTINUE
      CALL MULTP2(XNB2, XNB3, XNB)
      CALL DERIP2(XNB, 1, DNB1)
      CALL DERIP2(XNB, 2, DNB2)
      DO 32 I = 1, 9
         XN(2, I) = XNB(I)
         DN(1, 2, I) = DNB1(I)
         DN(2, 2, I) = DNB2(I)
   32 CONTINUE
      CALL MULTP2(XNB2, XNB4, XNB)
      CALL DERIP2(XNB, 1, DNB1)
      CALL DERIP2(XNB, 2, DNB2)
      DO 33 I = 1, 9
         XN(3, I) = XNB(I)
         DN(1, 3, I) = DNB1(I)
         DN(2, 3, I) = DNB2(I)
   33 CONTINUE
      CALL MULTP2(XNB1, XNB4, XNB)
      CALL DERIP2(XNB, 1, DNB1)
      CALL DERIP2(XNB, 2, DNB2)
      DO 34 I = 1, 9
         XN(4, I) = XNB(I)
         DN(1, 4, I) = DNB1(I)
         DN(2, 4, I) = DNB2(I)
   34 CONTINUE
C
C     Vecteurs tangents a la coque dans la configuration initiale
C
C     Initialisation
      DO 40 I = 1, 2
C          Boucle sur les parametres
         DO 40 J = 1, 3
C             Boucle sur les composantes
            DO 40 K = 1, 9
C                Boucle sur les coefficients des polynomes
               DX(I, J, K) = 0.D0
   40 CONTINUE
C     Calcul
      DO 50 I = 1, 2
C          Boucle sur les parametres
         DO 50 J = 1, 3
C             Boucle sur les composantes
            DO 50 K = 1, 9
C                Boucle sur les coefficients des polynomes
               DO 50 L = 1, 4
C                   Boucle sur les noeuds
            DX(I, J, K) = DX(I, J, K) + XX(J, L)*DN(I, L, K)
   50 CONTINUE
C
C
C     Calcul des termes de la matrice Kp
C
C     Initialisation
      DO 55 I = 1, 24
         DO 55 J = 1, 24
            XDUM(I,J) = 0.D0
   55       XKP(I, J) = 0. D0
C     Calcul
      DO 60 II = 1, 3
         DO 60 IL = 1, 4
            DO 60 IS = 1, 3
               DO 60 IT = 1, 4
                  XRES = 0.D0
                  DO 70 J = 1, 3
                     IFLAG = 0
                     DO 71 K = 1, 9
                        A(K) = XN(IL, K)
                        D(K) = 0.D0
   71                CONTINUE
                     IF (SKRO(II, J, IS) .NE. 0.D0) THEN
                        DO 81 K = 1, 9
                           B(K) = DN(2, IT, K)
                           C(K) = DX(1, J, K)
   81                   CONTINUE
                        CALL MULTP2(B, C, E)
                        DO 72 K = 1, 9
   72                      D(K) = SKRO(II, J, IS)*E(K)
                        IFLAG = 1
                     ENDIF
                     IF (SKRO(II, IS, J) .NE. 0.D0) THEN
                        DO 82 K = 1, 9
                           B(K) = DN(1, IT, K)
                           C(K) = DX(2, J, K)
   82                   CONTINUE
                        CALL MULTP2(B, C, E)
                        IF (IFLAG .EQ. 1) THEN
                           DO 73 K = 1, 9
   73                         D(K) = D(K) + SKRO(II, IS, J)*E(K)
                        ELSE
                           DO 74 K = 1, 9
   74                         D(K) = SKRO(II, IS, J)*E(K)
                        ENDIF
                        IFLAG = 1
                     ENDIF
                     IF (IFLAG .NE. 0) THEN
                        CALL MULTP2(A, D, XTIN)
                        CALL INTGP2(XTIN, TING)
                        XRES = XRES + TING
                     ENDIF
   70             CONTINUE
                  XRES = XRES*P
   60             XDUM(6*(IL-1) + II, 6*(IT-1) + IS) = XRES
      IF (IANT .EQ. 0) THEN
         DO 90 II = 1, 24
            DO 90 IJ = II, 24
               XKP(II, IJ) = (XDUM(II, IJ) + XDUM(IJ, II))*0.5D0
   90          XKP(IJ, II) = XKP(II, IJ)
      ELSE
         DO 91 II = 1, 24
            DO 91 IJ = 1, 24
   91          XKP(II, IJ) = XDUM(II, IJ)
      ENDIF
      RETURN
      END