C ELLA21 SOURCE CB215821 25/04/24 21:15:08 12248 SUBROUTINE ELLA21 (CARACT,INP,ZS,COOR,GAMA,ZC1,XPI) C IMPLICIT INTEGER(I-N) REAL*8 CARACT,COOR,GAMA,XPI,CE,CNU,CRO,CRIN,CREX, & CKC,CAM,CETA,CROF,CSON,CSE,CSEF,CC,CI,CIP, & COEF,X1,Y1,Z1,X2,Y2,Z2,XL,XI1,XI2,XI3, & GX,GY,GZ,GG,DELTA,DET,XJ1,XJ2,XJ3, & XK1,XK2,XK3 C COMPLEX*16 ZA,ZA1,ZC1,ZD,ZE,ZER,ZR, & ZI,ZS,ZCE,ZCP,ZCMU,ZCR,ZALFA2,ZALFAX,ZALFAR, & ZALFA4,ZALFAF,ZZ3,ZALFAA, & ZA56,ZB56,ZA89,ZB89,ZA1112,ZB1112, & ZG010,ZG110,ZG020,ZG120,ZG030,ZG130,ZG230, & ZG330,ZG430,ZG040,ZG140, & ZG011,ZG111,ZG021,ZG121,ZG031,ZG131,ZG231,ZG331, & ZG431,ZG041,ZG141 C INTEGER INP C C OPERATEUR ELFE LAPLACE ACOU C C CALCUL POUR LA POUTRE Num INP DE LA MATRICE CI TEL QUE CI * U = SM C C PARAMETRES : C C CARACT : TABLEAU DES CARACTERISTIQUE DES POUTRES (10,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 XPI : VALEUR PRECISE DE PI C C SORTIES : C C ZC1 : MATRICE 14*28 CALCULEE EN DIMENSIONNEL DANS LE REPERE C GLOBLAL C C C AUTEURS : SAINT-DIZIER ET GORCY C DATE : 24 JANVIER 1991 C C DIMENSION CARACT(10,*),COOR(3,*),GAMA(3,*) DIMENSION ZC1(14,28) DIMENSION ICHANG(28) DIMENSION ZR(3,3),ZD(14),ZE(14,14),ZA(14,28),ZA1(14,28) DIMENSION ZER(28,28) C C -------------- ACQUISITION DES CARACTERISTIQUES DE LA POUTRE C --------------------------------------------- C CE = CARACT( 1,INP) CNU = CARACT( 2,INP) CRO = CARACT( 3,INP) CRIN = CARACT( 4,INP) CREX = CARACT( 5,INP) CKC = CARACT( 6,INP) CAM = CARACT( 7,INP) CETA = CARACT( 8,INP) CROF = CARACT( 9,INP) CSON = CARACT(10,INP) C ZI = CMPLX(0.D0 , 1.D0) C CSE = XPI * ((CREX*CREX)-(CRIN*CRIN)) CSEF = XPI * (CRIN*CRIN) CC = XPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 2.D0 CI = XPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 4.D0 CIP = 2.D0 * CI C COEF = 1.D0+( 2.D0*CROF*CSON*CSON*CRIN / (CE*(CREX-CRIN))) CSON = CSON / SQRT(COEF) 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*CI*ZS*ZS*ZS*ZS/ZCMU/CKC * + CAM*CRO*CI*ZS*ZS*ZS/ZCMU/CKC/CSE * + ((CRO*CSE)+(CROF*CSEF))*ZS*ZS * + CAM*ZS ZALFA4= ZALFA4/(ZCE*CI) ZALFA2= SQRT(ZALFA4) ZALFAF= SQRT(ZALFA2) C ZZ3 = CRO*CI*(CMPLX(1.D0)+ZCE/(ZCMU*CKC))*ZS*ZS * + CAM*ZCE*CI*ZS/(ZCMU*CKC*CSE) ZZ3 = ZZ3 / (ZCE*CI*ZALFA2) C ZALFAA = ZS / CSON 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.1.D-12) THEN XJ1 = -XI2 XJ2 = 0.D0 XJ3 = 0.D0 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 END IF 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) = CMPLX(XJ2*XK3 - XJ3*XK2) ZR(1,2) = CMPLX(XJ3*XK1 - XJ1*XK3) ZR(1,3) = CMPLX(XJ1*XK2 - XJ2*XK1) ZR(2,1) = CMPLX(XI3*XK2 - XI2*XK3) ZR(2,2) = CMPLX(XI1*XK3 - XI3*XK1) ZR(2,3) = CMPLX(XI2*XK1 - XI1*XK2) ZR(3,1) = CMPLX(XI2*XJ3 - XI3*XJ2) ZR(3,2) = CMPLX(XI3*XJ1 - XI1*XJ3) ZR(3,3) = CMPLX(XI1*XJ2 - XI2*XJ1) 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*ZALFAF ZD( 6) = XL**2/ZCP*ZALFAF ZD( 7) = XL**2/ZCP*ZALFAX ZD( 8) = XL**2/ZCP*(ZALFAF**3) ZD( 9) = XL**2/ZCP*(ZALFAF**3) ZD(10) = XL/ZCR*ZALFAR ZD(11) = XL**2/ZCP*(ZALFAF**2) ZD(12) = XL**2/ZCP*(ZALFAF**2) ZD(13) = CMPLX(XL*CROF*CSON,0.D0) ZD(14) = XL*CROF*CSON*ZALFAA C C ------------------------------- MATRICE DES EFFORTS ZE C ---------------------- DO 20 I = 1 , 14 DO 21 J = 1 , 14 ZE(I,J) = CMPLX(0.D0,0.D0) 21 CONTINUE 20 CONTINUE C ZA56 = CMPLX(1.D0,0.D0) ZB56 = CMPLX(1.D0) / (ZCMU*CKC*CSE) C ZA89 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*ZCMU*CKC*CKC*CSE*CSE) * - CMPLX(1.D0)/(ZCE*CI) ZB89 = CRO*ZS*ZS*(CMPLX(1.D0)/(ZCMU*CKC)+CMPLX(1.D0)/ZCE) * + CAM*ZS/(ZCMU*CKC*CSE) C ZA1112= CMPLX(1.D0)/(ZCE*CI) ZB1112= (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*CKC*CSE) C ZE( 1, 1) = CMPLX(1.D0,0.D0) ZE( 2, 2) = CMPLX(1.D0,0.D0) ZE( 3, 3) = CMPLX(1.D0,0.D0) ZE( 4, 4) = CMPLX(1.D0,0.D0) ZE( 5, 5) = -ZA56 C ZE( 5, 9) = ZB56 ZE( 6, 6) = ZA56 ZE( 6, 8) = ZB56 ZE( 7, 7) = CMPLX(1.D0)/(ZCE*CSE) ZE( 8, 6) = ZB89 ZE( 8, 8) = ZA89 ZE( 9, 5) = -ZB89 C ZE( 9, 9) = ZA89 ZE(10,10) = CMPLX(1.D0)/(ZCMU*CC) C ZE(11, 3) = ZB1112 ZE(11,11) = -ZA1112 C ZE(12, 2) = ZB1112 ZE(12,12) = ZA1112 ZE(13,13) = CMPLX(1.D0,0.D0) ZE(14,14) = -ZS/CSEF C C -- CALCUL DES FONCTIONS DE GREEN C CALL ELLA31(0.D0,ZALFAX,ZALFAR,ZALFAF,ZALFAA, * ZG010,ZG110,ZG020,ZG120, * ZG030,ZG130,ZG230,ZG330,ZG430,ZZ3, * ZG040,ZG140) C CALL ELLA31(XL,ZALFAX,ZALFAR,ZALFAF,ZALFAA, * ZG011,ZG111,ZG021,ZG121, * ZG031,ZG131,ZG231,ZG331,ZG431,ZZ3, * ZG041,ZG141) C C -- REMPLISSAGE DE LA MATRICE A LOCAL DANS L'ORDRE INITIAL C DO 50 J = 1,28 DO 51 I = 1,14 ZA(I,J)= CMPLX(0.D0,0.D0) 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) =-ZG330+ZZ3*ZG130+CMPLX(1.D0) ZA( 9,18) =-ZG230+ZZ3*ZG030 ZA( 9,19) =-ZG130 ZA( 9,20) =-ZG030 ZA( 9,21) =-ZG331+ZZ3*ZG131 ZA( 9,22) = ZG231-ZZ3*ZG031 ZA( 9,23) =-ZG131 ZA( 9,24) = ZG031 ZA(10,17) =-ZG331+ZZ3*ZG131 ZA(10,18) =-ZG231+ZZ3*ZG031 ZA(10,19) =-ZG131 ZA(10,20) =-ZG031 ZA(10,21) =-ZG330+ZZ3*ZG130+CMPLX(1.D0) ZA(10,22) = ZG230-ZZ3*ZG030 ZA(10,23) =-ZG130 ZA(10,24) = ZG030 ZA(11,17) =-ZG430+ZZ3*ZG230 ZA(11,18) =-ZG330+ZZ3*ZG130+CMPLX(1.D0) ZA(11,19) =-ZG230 ZA(11,20) =-ZG130 ZA(11,21) = ZG431-ZZ3*ZG231 ZA(11,22) =-ZG331+ZZ3*ZG131 ZA(11,23) = ZG231 ZA(11,24) =-ZG131 ZA(12,17) =-ZG431+ZZ3*ZG231 ZA(12,18) =-ZG331+ZZ3*ZG131 ZA(12,19) =-ZG231 ZA(12,20) =-ZG131 ZA(12,21) = ZG430-ZZ3*ZG230 ZA(12,22) =-ZG330+ZZ3*ZG130+CMPLX(1.D0) ZA(12,23) = ZG230 ZA(12,24) =-ZG130 C ZA(13,25) = ZG140+CMPLX(1.D0) ZA(13,26) = ZG040 ZA(13,27) = ZG141 ZA(13,28) =-ZG041 ZA(14,25) = ZG141 ZA(14,26) = ZG041 ZA(14,27) = ZG140+CMPLX(1.D0) ZA(14,28) =-ZG040 C C -- TRANSFORMATION DE CETTE MATRICE DANS L'ORDRE U R F M P Q C C - ICHANG : MATRICE DE TRANSFORMATION DES COLONNES C ICHANG( 1) = 1 ICHANG( 2) = 7 ICHANG( 3) = 15 ICHANG( 4) = 21 ICHANG( 5) = 4 ICHANG( 6) = 10 ICHANG( 7) = 18 ICHANG( 8) = 24 ICHANG( 9) = 2 ICHANG(10) = 6 ICHANG(11) = 12 ICHANG(12) = 8 ICHANG(13) = 16 ICHANG(14) = 20 ICHANG(15) = 26 ICHANG(16) = 22 ICHANG(17) = 3 ICHANG(18) = 5 ICHANG(19) = 11 ICHANG(20) = 9 ICHANG(21) = 17 ICHANG(22) = 19 ICHANG(23) = 25 ICHANG(24) = 23 ICHANG(25) = 13 ICHANG(26) = 14 ICHANG(27) = 27 ICHANG(28) = 28 C DO 10 I = 1,28 K = ICHANG(I) DO 11 J = 1,14 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 , 28 DO 25 J = 1 , 28 ZER(I,J)= CMPLX(0.D0,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 ZER(13,13) = ZE(13,13) ZER(14,14) = ZE(14,14) C DO 34 L = 1,2 DO 35 I = 1,12 DO 36 J = 1,3 DO 37 K = 1,3 ZER(I+14,(L-1)*3+J+14) = ZER(I+14,(L-1)*3+J+14) * + ZE(I,(L-1)*3+K)*ZR(K,J) 37 CONTINUE 36 CONTINUE 35 CONTINUE 34 CONTINUE C DO 54 L = 3,4 DO 55 I = 1,12 DO 56 J = 1,3 DO 57 K = 1,3 ZER(I+14,(L-1)*3+J+14) = ZER(I+14,(L-1)*3+J+14) * - ZE(I,(L-1)*3+K)*ZR(K,J) 57 CONTINUE 56 CONTINUE 55 CONTINUE 54 CONTINUE C ZER(27,27) = ZE(13,13) ZER(28,28) = ZE(14,14) C DO 38 I = 1 , 14 DO 39 J = 1 , 28 ZER(I ,J) = ZER(I ,J)/ZD(I ) ZER(I+14,J) = ZER(I+14,J)/ZD(I ) 39 CONTINUE 38 CONTINUE C C --------------------------------- CALCUL DE ZA1*ZER C ----------------- DO 60 J = 1,28 DO 61 I = 1,14 ZC1(I,J)= CMPLX(0.D0,0.D0) 61 CONTINUE 60 CONTINUE C DO 40 I = 1,14 DO 41 J = 1,28 DO 42 K = 1,28 ZC1(I,J) = ZC1(I,J) + ZA1(I,K)*ZER(K,J) 42 CONTINUE 41 CONTINUE 40 CONTINUE C END