C ELLP21 SOURCE CB215821 25/04/24 21:15:11 12248 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 <--> 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 = CMPLX(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)= CMPLX(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