Sources Esope | Cast3M

ellp21
C ELLP21    SOURCE    KK2000    14/04/09    21:15:19     8027
        SUBROUTINE ELLP21 (CARACT,INP,ZS,COOR,GAMA,ICORR,ZC1)
C
      IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8     (A-H,O-Y)
        IMPLICIT COMPLEX*16 (Z)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   OPERATEUR ELFE LAPLACE POUTRE
C
C   CALCUL POUR LA POUTRE N INP DE LA MATRICE CI TEL QUE  CI * U = SM
C
C   PARAMETRES :
C
C   CARACT : TABLEAU DES CARACTERISTIQUE DES POUTRES (9,NP)
C   INP    : NUMERO DE LA POUTRE
C   ZS     : VALEUR DE S = S0 + I W
C   COOR   : TABLEAU DES COORDONNEES DES NOEUDS (3,2*NP)
C   GAMA   : TABLEAU DES VECTEUR DEFINISSANT LE PLAN LOCALE OXY (3,2*NP)
C   ICORR  : TABLEAU DE CORRESP. DE NOUEDS REELS &lt;--> NOEUDS FICTIFS
C            (2*NP)
C   SORTIES :
C
C   ZC1    : MATRICE 12*24 CALCULEE EN DIMENSIONNEL DANS LE REPERE
C            GLOBLAL
C
C
C   AUTEUR : SAINT-DIZIER
C   DATE   : 30 NOVEMBRE 1989 (VERSION DU 03 AVRIL 1990)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DIMENSION CARACT(12,*),COOR(3,*),GAMA(3,*),ICORR(*)
        DIMENSION ZC1(12,24)
        DIMENSION ICHANG(24)
        DIMENSION ZR(3,3),ZD(12),ZE(12,12),ZA(12,24),ZA1(12,24)
        DIMENSION ZER(24,24)
C
C -------------- ACQUISITION DES CARACTERISTIQUES DE LA POUTRE
C                ---------------------------------------------
        CE    = CARACT( 1,INP)
        CNU   = CARACT( 2,INP)
        CRO   = CARACT( 3,INP)
        CSE   = CARACT( 4,INP)
        CC    = CARACT( 5,INP)
        CIP   = CARACT( 6,INP)
        CIY   = CARACT( 7,INP)
        CIZ   = CARACT( 8,INP)
        CKCY  = CARACT( 9,INP)
        CKCZ  = CARACT(10,INP)
        CAM   = CARACT(11,INP)
        CETA  = CARACT(12,INP)
C
        ZI = (0.D0 , 1.D0)
C
        ZCE   = CE*(CMPLX(1.D0)+CETA*ZI)
        ZCP   = SQRT(ZCE/CRO)
        ZCMU  = ZCE/(CMPLX(2.D0)*(CMPLX(1.D0)+CNU))
        ZCR   = SQRT(ZCMU/CRO)
C
        ZALFA2= ZS*ZS*CRO*CSE + CAM*ZS
        ZALFA2= ZALFA2 / (ZCE*CSE)
        ZALFAX= SQRT(ZALFA2)
C
        ZALFA2= ZS*ZS*CRO*CIP + CAM*ZS
        ZALFA2= ZALFA2 / (ZCMU*CC)
        ZALFAR= SQRT(ZALFA2)
C
        ZALFA4= CRO*CRO*CIZ*ZS*ZS*ZS*ZS/ZCMU/CKCY
     *        + CAM*CRO*CIZ*ZS*ZS*ZS/ZCMU/CKCY/CSE
     *        + CRO*CSE*ZS*ZS
     *        + CAM*ZS
        ZALFA4= ZALFA4/(ZCE*CIZ)
        ZALFA2= SQRT(ZALFA4)
        ZALFAY= SQRT(ZALFA2)
C
        ZZ3   = CRO*CIZ*(CMPLX(1.D0)+ZCE/(ZCMU*CKCY))*ZS*ZS
     *        + CAM*ZCE*CIZ*ZS/(ZCMU*CKCY*CSE)
        ZZ3   = ZZ3 / (ZCE*CIZ*ZALFA2)
C
        ZALFA4= CRO*CRO*CIY*ZS*ZS*ZS*ZS/ZCMU/CKCZ
     *        + CAM*CRO*CIY*ZS*ZS*ZS/ZCMU/CKCZ/CSE
     *        + CRO*CSE*ZS*ZS
     *        + CAM*ZS
        ZALFA4= ZALFA4/(ZCE*CIY)
        ZALFA2= SQRT(ZALFA4)
        ZALFAZ= SQRT(ZALFA2)
C
        ZZ4   = CRO*CIY*(CMPLX(1.D0)+ZCE/(ZCMU*CKCZ))*ZS*ZS
     *        + CAM*ZCE*CIY*ZS/(ZCMU*CKCZ*CSE)
        ZZ4   = ZZ4 / (ZCE*CIY*ZALFA2)
C
        N1  = 2*INP-1
        N2  = 2*INP
        X1  = COOR (1,N1)
        Y1  = COOR (2,N1)
        Z1  = COOR (3,N1)
        X2  = COOR (1,N2)
        Y2  = COOR (2,N2)
        Z2  = COOR (3,N2)
C
        XL = SQRT((X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2)
C
C ------------------------------ VECTEUR UNITAIRE OX REPERE LOCALE
C                                ---------------------------------
        XI1 = (X2-X1)/XL
        XI2 = (Y2-Y1)/XL
        XI3 = (Z2-Z1)/XL
C
C ------------------------- VECTEUR UNITAIRE OY REPERE LOCALE
C                           ---------------------------------
        GX = GAMA(1,INP)
        GY = GAMA(2,INP)
        GZ = GAMA(3,INP)
        GG = SQRT(GX*GX + GY*GY + GZ*GZ)
        GX = GX/GG
        GY = GY/GG
        GZ = GZ/GG
C
        DELTA =  SQRT (1.D0 - (XI1*GX + XI2*GY + XI3*GZ)**2)
C
        DET=-(GX*XI3-GZ*XI1)**2-(GY*XI1-GX*XI2)**2-(GY*XI3-GZ*XI2)**2
C
        IF (ABS(DET).LT.1E-12) THEN
           XJ1 = -XI2
           XJ2 = 0.
           XJ3 = 0.
        ELSE
           XJ1 = (XI2*(GY*XI1-GX*XI2)-XI3*(GX*XI3-GZ*XI1))*DELTA/DET
           XJ2 = (XI3*(GZ*XI2-GY*XI3)-XI1*(GY*XI1-GX*XI2))*DELTA/DET
           XJ3 = (XI1*(GX*XI3-GZ*XI1)-XI2*(GZ*XI2-GY*XI3))*DELTA/DET
        ENDIF
C
C ---------------------------- VECTEUR UNITAIRE OZ REPERE LOCALE
C                              ---------------------------------
        XK1 = XI2*XJ3 - XI3*XJ2
        XK2 = XI3*XJ1 - XI1*XJ3
        XK3 = XI1*XJ2 - XI2*XJ1
C
        ZR(1,1) = XJ2*XK3 - XJ3*XK2
        ZR(1,2) = XJ3*XK1 - XJ1*XK3
        ZR(1,3) = XJ1*XK2 - XJ2*XK1
        ZR(2,1) = XI3*XK2 - XI2*XK3
        ZR(2,2) = XI1*XK3 - XI3*XK1
        ZR(2,3) = XI2*XK1 - XI1*XK2
        ZR(3,1) = XI2*XJ3 - XI3*XJ2
        ZR(3,2) = XI3*XJ1 - XI1*XJ3
        ZR(3,3) = XI1*XJ2 - XI2*XJ1
C
C
C -----------------------------  MATRICE DIMENSIONNEMENT ZD
C                                --------------------------
        ZD( 1) = XL**2/ZCP
        ZD( 2) = XL**2/ZCP
        ZD( 3) = XL**2/ZCP
        ZD( 4) = XL/ZCR
        ZD( 5) = XL**2/ZCP*ZALFAZ
        ZD( 6) = XL**2/ZCP*ZALFAY
        ZD( 7) = XL**2/ZCP*ZALFAX
        ZD( 8) = XL**2/ZCP*(ZALFAY**3)
        ZD( 9) = XL**2/ZCP*(ZALFAZ**3)
        ZD(10) = XL/ZCR*ZALFAR
        ZD(11) = XL**2/ZCP*(ZALFAZ**2)
        ZD(12) = XL**2/ZCP*(ZALFAY**2)
C
C ------------------------------- MATRICE DES EFFORTS ZE
C                                 ----------------------
        DO 20 I = 1 , 12
            DO 21 J = 1 , 12
                ZE(I,J) = CMPLX(0.D0)
21          CONTINUE
20      CONTINUE
C
C        ZEPSY = CRO*CIZ*ZS*ZS/(ZCMU*CKCY*CSE)
C        ZEPSZ = CRO*CIY*ZS*ZS/(ZCMU*CKCZ*CSE)
C
        ZA5 = CMPLX(1.D0)
        ZA6 = CMPLX(1.D0)
        ZB6 = CMPLX(1.D0) / (ZCMU*CKCY*CSE)
        ZB5 = CMPLX(1.D0) / (ZCMU*CKCZ*CSE)
C
        ZA8 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*ZCMU*CKCY*CKCY*CSE*CSE)
     *      -  CMPLX(1.D0)/(ZCE*CIZ)
        ZA9 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*ZCMU*CKCZ*CKCZ*CSE*CSE)
     *      -  CMPLX(1.D0)/(ZCE*CIY)
        ZB8 = CRO*ZS*ZS*(CMPLX(1.D0)/(ZCMU*CKCY)+CMPLX(1.D0)/ZCE)
     *      + CAM*ZS/(ZCMU*CKCY*CSE)
        ZB9 = CRO*ZS*ZS*(CMPLX(1.D0)/(ZCMU*CKCZ)+CMPLX(1.D0)/ZCE)
     *      + CAM*ZS/(ZCMU*CKCZ*CSE)
C
        ZA12= CMPLX(1.D0)/(ZCE*CIZ)
        ZA11= CMPLX(1.D0)/(ZCE*CIY)
        ZB12= (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*CKCY*CSE)
        ZB11= (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*CKCZ*CSE)
C
C
        ZE( 1, 1) = CMPLX(1.D0)
        ZE( 2, 2) = CMPLX(1.D0)
        ZE( 3, 3) = CMPLX(1.D0)
        ZE( 4, 4) = CMPLX(1.D0)
        ZE( 5, 5) = ZA5
        ZE( 5, 9) = ZB5
        ZE( 6, 6) = ZA6
        ZE( 6, 8) = ZB6
        ZE( 7, 7) = CMPLX(1.D0)/(ZCE*CSE)
        ZE( 8, 6) = ZB8
        ZE( 8, 8) = ZA8
        ZE( 9, 5) = ZB9
        ZE( 9, 9) = ZA9
        ZE(10,10) =-CMPLX(1.D0)/(ZCMU*CC)
        ZE(11, 3) = ZB11
        ZE(11,11) = ZA11
        ZE(12, 2) = ZB12
        ZE(12,12) = ZA12
C
C -- CALCUL DES FONCTIONS DE GREEN
C
        CALL ELLP31(0.D0,INP,ZS,ZALFAX,ZALFAR,ZALFAY,ZALFAZ,
     *              ZG010,ZG110,ZG020,ZG120,
     *              ZG030,ZG130,ZG230,ZG330,ZG430,ZZ3,
     *              ZG040,ZG140,ZG240,ZG340,ZG440,ZZ4)
C
        CALL ELLP31(XL,INP,ZS,ZALFAX,ZALFAR,ZALFAY,ZALFAZ,
     *              ZG011,ZG111,ZG021,ZG121,
     *              ZG031,ZG131,ZG231,ZG331,ZG431,ZZ3,
     *              ZG041,ZG141,ZG241,ZG341,ZG441,ZZ4)
C
C -- REMPLISSAGE DE LA MATRICE A LOCAL DANS L'ORDRE INITIAL
C
        DO 50 J = 1,24
            DO 51 I = 1,12
                ZA(I,J)= 0.
51          CONTINUE
50      CONTINUE
C
        ZA( 1, 1) = ZG110+CMPLX(1.D0)
        ZA( 1, 2) = ZG010
        ZA( 1, 3) = ZG111
        ZA( 1, 4) =-ZG011
        ZA( 2, 1) = ZG111
        ZA( 2, 2) = ZG011
        ZA( 2, 3) = ZG110+CMPLX(1.D0)
        ZA( 2, 4) =-ZG010
C
        ZA( 3, 5) = ZG120+CMPLX(1.D0)
        ZA( 3, 6) = ZG020
        ZA( 3, 7) = ZG121
        ZA( 3, 8) =-ZG021
        ZA( 4, 5) = ZG121
        ZA( 4, 6) = ZG021
        ZA( 4, 7) = ZG120+CMPLX(1.D0)
        ZA( 4, 8) =-ZG020
C
        ZA( 5, 9) =-ZG330+ZZ3*ZG130+CMPLX(1.D0)
        ZA( 5,10) =-ZG230+ZZ3*ZG030
        ZA( 5,11) =-ZG130
        ZA( 5,12) =-ZG030
        ZA( 5,13) =-ZG331+ZZ3*ZG131
        ZA( 5,14) = ZG231-ZZ3*ZG031
        ZA( 5,15) =-ZG131
        ZA( 5,16) = ZG031
        ZA( 6, 9) =-ZG331+ZZ3*ZG131
        ZA( 6,10) =-ZG231+ZZ3*ZG031
        ZA( 6,11) =-ZG131
        ZA( 6,12) =-ZG031
        ZA( 6,13) =-ZG330+ZZ3*ZG130+CMPLX(1.D0)
        ZA( 6,14) = ZG230-ZZ3*ZG030
        ZA( 6,15) =-ZG130
        ZA( 6,16) = ZG030
        ZA( 7, 9) =-ZG430+ZZ3*ZG230
        ZA( 7,10) =-ZG330+ZZ3*ZG130+CMPLX(1.D0)
        ZA( 7,11) =-ZG230
        ZA( 7,12) =-ZG130
        ZA( 7,13) = ZG431-ZZ3*ZG231
        ZA( 7,14) =-ZG331+ZZ3*ZG131
        ZA( 7,15) = ZG231
        ZA( 7,16) =-ZG131
        ZA( 8, 9) =-ZG431+ZZ3*ZG231
        ZA( 8,10) =-ZG331+ZZ3*ZG131
        ZA( 8,11) =-ZG231
        ZA( 8,12) =-ZG131
        ZA( 8,13) = ZG430-ZZ3*ZG230
        ZA( 8,14) =-ZG330+ZZ3*ZG130+CMPLX(1.D0)
        ZA( 8,15) = ZG230
        ZA( 8,16) =-ZG130
C
        ZA( 9,17) =-ZG340+ZZ4*ZG140+CMPLX(1.D0)
        ZA( 9,18) =-ZG240+ZZ4*ZG040
        ZA( 9,19) =-ZG140
        ZA( 9,20) =-ZG040
        ZA( 9,21) =-ZG341+ZZ4*ZG141
        ZA( 9,22) = ZG241-ZZ4*ZG041
        ZA( 9,23) =-ZG141
        ZA( 9,24) = ZG041
        ZA(10,17) =-ZG341+ZZ4*ZG141
        ZA(10,18) =-ZG241+ZZ4*ZG041
        ZA(10,19) =-ZG141
        ZA(10,20) =-ZG041
        ZA(10,21) =-ZG340+ZZ4*ZG140+CMPLX(1.D0)
        ZA(10,22) = ZG240-ZZ4*ZG040
        ZA(10,23) =-ZG140
        ZA(10,24) = ZG040
        ZA(11,17) =-ZG440+ZZ4*ZG240
        ZA(11,18) =-ZG340+ZZ4*ZG140+CMPLX(1.D0)
        ZA(11,19) =-ZG240
        ZA(11,20) =-ZG140
        ZA(11,21) = ZG441-ZZ4*ZG241
        ZA(11,22) =-ZG341+ZZ4*ZG141
        ZA(11,23) = ZG241
        ZA(11,24) =-ZG141
        ZA(12,17) =-ZG441+ZZ4*ZG241
        ZA(12,18) =-ZG341+ZZ4*ZG141
        ZA(12,19) =-ZG241
        ZA(12,20) =-ZG141
        ZA(12,21) = ZG440-ZZ4*ZG240
        ZA(12,22) =-ZG340+ZZ4*ZG140+CMPLX(1.D0)
        ZA(12,23) = ZG240
        ZA(12,24) =-ZG140
C
C -- TRANSFORMATION DE CETTE MATRICE DANS L'ORDRE U R F M
C
C - ICHANG : MATRICE DE TRANSFORMATION DES COLONNES
C
        ICHANG( 1) = 1
        ICHANG( 2) = 7
        ICHANG( 3) = 13
        ICHANG( 4) = 19
        ICHANG( 5) = 4
        ICHANG( 6) = 10
        ICHANG( 7) = 16
        ICHANG( 8) = 22
        ICHANG( 9) = 2
        ICHANG(10) = 6
        ICHANG(11) = 12
        ICHANG(12) = 8
        ICHANG(13) = 14
        ICHANG(14) = 18
        ICHANG(15) = 24
        ICHANG(16) = 20
        ICHANG(17) = 3
        ICHANG(18) = 5
        ICHANG(19) = 11
        ICHANG(20) = 9
        ICHANG(21) = 15
        ICHANG(22) = 17
        ICHANG(23) = 23
        ICHANG(24) = 21
C
        DO 10 I = 1,24
            K = ICHANG(I)
            DO 11 J = 1,12
                ZA1(J,K) = ZA(J,I)
11          CONTINUE
10      CONTINUE
C
C
C -- CACUL DE ZCI = ZA1*ZR (DANS CET ORDRE CAR LE PRODUIT N'EST
C -- PAS COMMUTATIF
C
C
C - CALCUL DE  ZR :
C
        DO 26 I = 1 , 24
            DO 25 J = 1 , 24
               ZER(I,J)= 0.D0
25          CONTINUE
26      CONTINUE
C
        DO 30 L = 1,4
            DO 31 I = 1,12
                DO 32 J = 1,3
                    DO 33 K = 1,3
        ZER(I,(L-1)*3+J) = ZER(I,(L-1)*3+J) + ZE(I,(L-1)*3+K)*ZR(K,J)
33                  CONTINUE
32              CONTINUE
31          CONTINUE
30      CONTINUE
C
C
        DO 34 L = 1,2
            DO 35 I = 1,12
                DO 36 J = 1,3
                    DO 37 K = 1,3
                        ZER(I+12,(L-1)*3+J+12) = ZER(I+12,(L-1)*3+J+12)
     *                  + ZE(I,(L-1)*3+K)*ZR(K,J)
37                  CONTINUE
36              CONTINUE
35          CONTINUE
34      CONTINUE
C
C
        DO 54 L = 3,4
            DO 55 I = 1,12
                DO 56 J = 1,3
                    DO 57 K = 1,3
                        ZER(I+12,(L-1)*3+J+12) = ZER(I+12,(L-1)*3+J+12)
     *                  - ZE(I,(L-1)*3+K)*ZR(K,J)
57                  CONTINUE
56              CONTINUE
55          CONTINUE
54      CONTINUE
C
C
        DO 38 I = 1 , 12
            DO 39 J = 1 , 24
                ZER(I   ,J) = ZER(I   ,J)/ZD(I  )
                ZER(I+12,J) = ZER(I+12,J)/ZD(I  )
39          CONTINUE
38      CONTINUE
C
C --------------------------------- CALCUL DE ZA1*ZER
C                                   -----------------
        DO 60 J = 1,24
            DO 61 I = 1,12
                ZC1(I,J)= (0.D0,0.D0)
61          CONTINUE
60      CONTINUE
C
        DO 40 I = 1,12
            DO 41 J = 1,24
                DO 42 K = 1,24
                    ZC1(I,J) = ZC1(I,J) + ZA1(I,K)*ZER(K,J)
42              CONTINUE
41         CONTINUE
40      CONTINUE
C
        END