1 : C ELLA21 SOURCE CHAT 05/01/12 23:34:09 5004 2 : 3 : SUBROUTINE ELLA21 (CARACT,INP,ZS,COOR,GAMA,ZC1,XPI) 4 : C 5 : IMPLICIT INTEGER(I-N) 6 : REAL*8 CARACT,COOR,GAMA,XPI,CE,CNU,CRO,CRIN,CREX, 7 : & CKC,CAM,CETA,CROF,CSON,CSE,CSEF,CC,CI,CIP, 8 : & COEF,X1,Y1,Z1,X2,Y2,Z2,XL,XI1,XI2,XI3, 9 : & GX,GY,GZ,GG,DELTA,DET,XJ1,XJ2,XJ3, 10 : & XK1,XK2,XK3 11 : C 12 : COMPLEX*16 ZA,ZA1,ZC1,ZD,ZE,ZER,ZR, 13 : & ZI,ZS,ZCE,ZCP,ZCMU,ZCR,ZALFA2,ZALFAX,ZALFAR, 14 : & ZALFA4,ZALFAF,ZZ3,ZALFAA, 15 : & ZA56,ZB56,ZA89,ZB89,ZA1112,ZB1112, 16 : & ZG010,ZG110,ZG020,ZG120,ZG030,ZG130,ZG230, 17 : & ZG330,ZG430,ZG040,ZG140, 18 : & ZG011,ZG111,ZG021,ZG121,ZG031,ZG131,ZG231,ZG331, 19 : & ZG431,ZG041,ZG141 20 : C 21 : INTEGER INP 22 : C 23 : 24 : 25 : C OPERATEUR ELFE LAPLACE ACOU 26 : C 27 : C CALCUL POUR LA POUTRE Num INP DE LA MATRICE CI TEL QUE CI * U = SM 28 : C 29 : C PARAMETRES : 30 : C 31 : C CARACT : TABLEAU DES CARACTERISTIQUE DES POUTRES (10,NP) 32 : C INP : NUMERO DE LA POUTRE 33 : C ZS : VALEUR DE S = S0 + I W 34 : C COOR : TABLEAU DES COORDONNEES DES NOEUDS (3,2*NP) 35 : C GAMA : TABLEAU DES VECTEUR DEFINISSANT LE PLAN LOCALE OXY (3,2*NP) 36 : C XPI : VALEUR PRECISE DE PI 37 : C 38 : C SORTIES : 39 : C 40 : C ZC1 : MATRICE 14*28 CALCULEE EN DIMENSIONNEL DANS LE REPERE 41 : C GLOBLAL 42 : C 43 : C 44 : C AUTEURS : SAINT-DIZIER ET GORCY 45 : C DATE : 24 JANVIER 1991 46 : C 47 : C 48 : DIMENSION CARACT(10,*),COOR(3,*),GAMA(3,*) 49 : DIMENSION ZC1(14,28) 50 : DIMENSION ICHANG(28) 51 : DIMENSION ZR(3,3),ZD(14),ZE(14,14),ZA(14,28),ZA1(14,28) 52 : DIMENSION ZER(28,28) 53 : C 54 : C -------------- ACQUISITION DES CARACTERISTIQUES DE LA POUTRE 55 : C --------------------------------------------- 56 : C 57 : CE = CARACT( 1,INP) 58 : CNU = CARACT( 2,INP) 59 : CRO = CARACT( 3,INP) 60 : CRIN = CARACT( 4,INP) 61 : CREX = CARACT( 5,INP) 62 : CKC = CARACT( 6,INP) 63 : CAM = CARACT( 7,INP) 64 : CETA = CARACT( 8,INP) 65 : CROF = CARACT( 9,INP) 66 : CSON = CARACT(10,INP) 67 : C 68 : ZI = CMPLX(0.D0 , 1.D0) 69 : C 70 : CSE = XPI * ((CREX*CREX)-(CRIN*CRIN)) 71 : CSEF = XPI * (CRIN*CRIN) 72 : CC = XPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 2.D0 73 : CI = XPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 4.D0 74 : 75 : CIP = 2.D0 * CI 76 : C 77 : COEF = 1.D0+( 2.D0*CROF*CSON*CSON*CRIN / (CE*(CREX-CRIN))) 78 : CSON = CSON / SQRT(COEF) 79 : C 80 : ZCE = CE*(1.D0+CETA*ZI) 81 : ZCP = SQRT(ZCE/CRO) 82 : ZCMU = ZCE/(2.D0*(1.D0+CNU)) 83 : ZCR = SQRT(ZCMU/CRO) 84 : C 85 : ZALFA2= ZS*ZS*CRO*CSE + CAM*ZS 86 : ZALFA2= ZALFA2 / (ZCE*CSE) 87 : ZALFAX= SQRT(ZALFA2) 88 : C 89 : ZALFA2= ZS*ZS*CRO*CIP + CAM*ZS 90 : ZALFA2= ZALFA2 / (ZCMU*CC) 91 : ZALFAR= SQRT(ZALFA2) 92 : C 93 : ZALFA4= CRO*CRO*CI*ZS*ZS*ZS*ZS/ZCMU/CKC 94 : * + CAM*CRO*CI*ZS*ZS*ZS/ZCMU/CKC/CSE 95 : * + ((CRO*CSE)+(CROF*CSEF))*ZS*ZS 96 : * + CAM*ZS 97 : ZALFA4= ZALFA4/(ZCE*CI) 98 : ZALFA2= SQRT(ZALFA4) 99 : ZALFAF= SQRT(ZALFA2) 100 : C 101 : ZZ3 = CRO*CI*(1.D0+ZCE/(ZCMU*CKC))*ZS*ZS 102 : * + CAM*ZCE*CI*ZS/(ZCMU*CKC*CSE) 103 : ZZ3 = ZZ3 / (ZCE*CI*ZALFA2) 104 : C 105 : ZALFAA = ZS / CSON 106 : C 107 : N1 = 2*INP-1 108 : N2 = 2*INP 109 : X1 = COOR (1,N1) 110 : Y1 = COOR (2,N1) 111 : Z1 = COOR (3,N1) 112 : X2 = COOR (1,N2) 113 : Y2 = COOR (2,N2) 114 : Z2 = COOR (3,N2) 115 : C 116 : XL = SQRT((X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2) 117 : C 118 : C ------------------------------ VECTEUR UNITAIRE OX REPERE LOCALE 119 : C --------------------------------- 120 : XI1 = (X2-X1)/XL 121 : XI2 = (Y2-Y1)/XL 122 : XI3 = (Z2-Z1)/XL 123 : C 124 : C ------------------------- VECTEUR UNITAIRE OY REPERE LOCALE 125 : C --------------------------------- 126 : GX = GAMA(1,INP) 127 : GY = GAMA(2,INP) 128 : GZ = GAMA(3,INP) 129 : GG = SQRT(GX*GX + GY*GY + GZ*GZ) 130 : GX = GX/GG 131 : GY = GY/GG 132 : GZ = GZ/GG 133 : C 134 : DELTA = SQRT (1.D0 - (XI1*GX + XI2*GY + XI3*GZ)**2) 135 : C 136 : DET=-(GX*XI3-GZ*XI1)**2-(GY*XI1-GX*XI2)**2-(GY*XI3-GZ*XI2)**2 137 : C 138 : IF (ABS(DET).LT.1.D-12) THEN 139 : XJ1 = -XI2 140 : XJ2 = 0.D0 141 : XJ3 = 0.D0 142 : ELSE 143 : XJ1 = (XI2*(GY*XI1-GX*XI2)-XI3*(GX*XI3-GZ*XI1))*DELTA/DET 144 : XJ2 = (XI3*(GZ*XI2-GY*XI3)-XI1*(GY*XI1-GX*XI2))*DELTA/DET 145 : XJ3 = (XI1*(GX*XI3-GZ*XI1)-XI2*(GZ*XI2-GY*XI3))*DELTA/DET 146 : END IF 147 : C 148 : C ---------------------------- VECTEUR UNITAIRE OZ REPERE LOCALE 149 : C --------------------------------- 150 : XK1 = XI2*XJ3 - XI3*XJ2 151 : XK2 = XI3*XJ1 - XI1*XJ3 152 : XK3 = XI1*XJ2 - XI2*XJ1 153 : C 154 : ZR(1,1) = CMPLX(XJ2*XK3 - XJ3*XK2) 155 : ZR(1,2) = CMPLX(XJ3*XK1 - XJ1*XK3) 156 : ZR(1,3) = CMPLX(XJ1*XK2 - XJ2*XK1) 157 : ZR(2,1) = CMPLX(XI3*XK2 - XI2*XK3) 158 : ZR(2,2) = CMPLX(XI1*XK3 - XI3*XK1) 159 : ZR(2,3) = CMPLX(XI2*XK1 - XI1*XK2) 160 : ZR(3,1) = CMPLX(XI2*XJ3 - XI3*XJ2) 161 : ZR(3,2) = CMPLX(XI3*XJ1 - XI1*XJ3) 162 : ZR(3,3) = CMPLX(XI1*XJ2 - XI2*XJ1) 163 : C 164 : C ----------------------------- MATRICE DIMENSIONNEMENT ZD 165 : C -------------------------- 166 : ZD( 1) = XL**2/ZCP 167 : ZD( 2) = XL**2/ZCP 168 : ZD( 3) = XL**2/ZCP 169 : ZD( 4) = XL/ZCR 170 : ZD( 5) = XL**2/ZCP*ZALFAF 171 : ZD( 6) = XL**2/ZCP*ZALFAF 172 : ZD( 7) = XL**2/ZCP*ZALFAX 173 : ZD( 8) = XL**2/ZCP*(ZALFAF**3) 174 : ZD( 9) = XL**2/ZCP*(ZALFAF**3) 175 : ZD(10) = XL/ZCR*ZALFAR 176 : ZD(11) = XL**2/ZCP*(ZALFAF**2) 177 : ZD(12) = XL**2/ZCP*(ZALFAF**2) 178 : ZD(13) = CMPLX(XL*CROF*CSON,0.D0) 179 : ZD(14) = XL*CROF*CSON*ZALFAA 180 : C 181 : C ------------------------------- MATRICE DES EFFORTS ZE 182 : C ---------------------- 183 : DO 20 I = 1 , 14 184 : DO 21 J = 1 , 14 185 : ZE(I,J) = CMPLX(0.D0,0.D0) 186 : 21 CONTINUE 187 : 20 CONTINUE 188 : C 189 : ZA56 = CMPLX(1.D0,0.D0) 190 : ZB56 = 1.D0 / (ZCMU*CKC*CSE) 191 : C 192 : ZA89 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*ZCMU*CKC*CKC*CSE*CSE) 193 : * - 1.D0/(ZCE*CI) 194 : ZB89 = CRO*ZS*ZS*(1.D0/(ZCMU*CKC)+1.D0/ZCE) 195 : * + CAM*ZS/(ZCMU*CKC*CSE) 196 : C 197 : ZA1112= 1.D0/(ZCE*CI) 198 : ZB1112= (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*CKC*CSE) 199 : C 200 : ZE( 1, 1) = CMPLX(1.D0,0.D0) 201 : ZE( 2, 2) = CMPLX(1.D0,0.D0) 202 : ZE( 3, 3) = CMPLX(1.D0,0.D0) 203 : ZE( 4, 4) = CMPLX(1.D0,0.D0) 204 : ZE( 5, 5) = -ZA56 205 : C 206 : ZE( 5, 9) = ZB56 207 : ZE( 6, 6) = ZA56 208 : ZE( 6, 8) = ZB56 209 : ZE( 7, 7) = 1.D0/(ZCE*CSE) 210 : ZE( 8, 6) = ZB89 211 : ZE( 8, 8) = ZA89 212 : ZE( 9, 5) = -ZB89 213 : C 214 : ZE( 9, 9) = ZA89 215 : ZE(10,10) = 1.D0/(ZCMU*CC) 216 : C 217 : ZE(11, 3) = ZB1112 218 : ZE(11,11) = -ZA1112 219 : C 220 : ZE(12, 2) = ZB1112 221 : ZE(12,12) = ZA1112 222 : ZE(13,13) = CMPLX(1.D0,0.D0) 223 : ZE(14,14) = -ZS/CSEF 224 : C 225 : C -- CALCUL DES FONCTIONS DE GREEN 226 : C 226 : CALL ELLA31(0.D0,ZALFAX,ZALFAR,ZALFAF,ZALFAA, 227 : * ZG010,ZG110,ZG020,ZG120, 228 : * ZG030,ZG130,ZG230,ZG330,ZG430,ZZ3, 229 : * ZG040,ZG140) 230 : C 230 : CALL ELLA31(XL,ZALFAX,ZALFAR,ZALFAF,ZALFAA, 231 : * ZG011,ZG111,ZG021,ZG121, 232 : * ZG031,ZG131,ZG231,ZG331,ZG431,ZZ3, 233 : * ZG041,ZG141) 234 : C 235 : C -- REMPLISSAGE DE LA MATRICE A LOCAL DANS L'ORDRE INITIAL 236 : C 237 : DO 50 J = 1,28 238 : DO 51 I = 1,14 239 : ZA(I,J)= CMPLX(0.D0,0.D0) 240 : 51 CONTINUE 241 : 50 CONTINUE 242 : C 243 : ZA( 1, 1) = ZG110+1.D0 244 : ZA( 1, 2) = ZG010 245 : ZA( 1, 3) = ZG111 246 : ZA( 1, 4) =-ZG011 247 : ZA( 2, 1) = ZG111 248 : ZA( 2, 2) = ZG011 249 : ZA( 2, 3) = ZG110+1.D0 250 : ZA( 2, 4) =-ZG010 251 : C 252 : ZA( 3, 5) = ZG120+1.D0 253 : ZA( 3, 6) = ZG020 254 : ZA( 3, 7) = ZG121 255 : ZA( 3, 8) =-ZG021 256 : ZA( 4, 5) = ZG121 257 : ZA( 4, 6) = ZG021 258 : ZA( 4, 7) = ZG120+1.D0 259 : ZA( 4, 8) =-ZG020 260 : C 261 : ZA( 5, 9) =-ZG330+ZZ3*ZG130+1.D0 262 : ZA( 5,10) =-ZG230+ZZ3*ZG030 263 : ZA( 5,11) =-ZG130 264 : ZA( 5,12) =-ZG030 265 : ZA( 5,13) =-ZG331+ZZ3*ZG131 266 : ZA( 5,14) = ZG231-ZZ3*ZG031 267 : ZA( 5,15) =-ZG131 268 : ZA( 5,16) = ZG031 269 : ZA( 6, 9) =-ZG331+ZZ3*ZG131 270 : ZA( 6,10) =-ZG231+ZZ3*ZG031 271 : ZA( 6,11) =-ZG131 272 : ZA( 6,12) =-ZG031 273 : ZA( 6,13) =-ZG330+ZZ3*ZG130+1.D0 274 : ZA( 6,14) = ZG230-ZZ3*ZG030 275 : ZA( 6,15) =-ZG130 276 : ZA( 6,16) = ZG030 277 : ZA( 7, 9) =-ZG430+ZZ3*ZG230 278 : ZA( 7,10) =-ZG330+ZZ3*ZG130+1.D0 279 : ZA( 7,11) =-ZG230 280 : ZA( 7,12) =-ZG130 281 : ZA( 7,13) = ZG431-ZZ3*ZG231 282 : ZA( 7,14) =-ZG331+ZZ3*ZG131 283 : ZA( 7,15) = ZG231 284 : ZA( 7,16) =-ZG131 285 : ZA( 8, 9) =-ZG431+ZZ3*ZG231 286 : ZA( 8,10) =-ZG331+ZZ3*ZG131 287 : ZA( 8,11) =-ZG231 288 : ZA( 8,12) =-ZG131 289 : ZA( 8,13) = ZG430-ZZ3*ZG230 290 : ZA( 8,14) =-ZG330+ZZ3*ZG130+1.D0 291 : ZA( 8,15) = ZG230 292 : ZA( 8,16) =-ZG130 293 : C 294 : ZA( 9,17) =-ZG330+ZZ3*ZG130+1.D0 295 : ZA( 9,18) =-ZG230+ZZ3*ZG030 296 : ZA( 9,19) =-ZG130 297 : ZA( 9,20) =-ZG030 298 : ZA( 9,21) =-ZG331+ZZ3*ZG131 299 : ZA( 9,22) = ZG231-ZZ3*ZG031 300 : ZA( 9,23) =-ZG131 301 : ZA( 9,24) = ZG031 302 : ZA(10,17) =-ZG331+ZZ3*ZG131 303 : ZA(10,18) =-ZG231+ZZ3*ZG031 304 : ZA(10,19) =-ZG131 305 : ZA(10,20) =-ZG031 306 : ZA(10,21) =-ZG330+ZZ3*ZG130+1.D0 307 : ZA(10,22) = ZG230-ZZ3*ZG030 308 : ZA(10,23) =-ZG130 309 : ZA(10,24) = ZG030 310 : ZA(11,17) =-ZG430+ZZ3*ZG230 311 : ZA(11,18) =-ZG330+ZZ3*ZG130+1.D0 312 : ZA(11,19) =-ZG230 313 : ZA(11,20) =-ZG130 314 : ZA(11,21) = ZG431-ZZ3*ZG231 315 : ZA(11,22) =-ZG331+ZZ3*ZG131 316 : ZA(11,23) = ZG231 317 : ZA(11,24) =-ZG131 318 : ZA(12,17) =-ZG431+ZZ3*ZG231 319 : ZA(12,18) =-ZG331+ZZ3*ZG131 320 : ZA(12,19) =-ZG231 321 : ZA(12,20) =-ZG131 322 : ZA(12,21) = ZG430-ZZ3*ZG230 323 : ZA(12,22) =-ZG330+ZZ3*ZG130+1.D0 324 : ZA(12,23) = ZG230 325 : ZA(12,24) =-ZG130 326 : C 327 : ZA(13,25) = ZG140+1.D0 328 : ZA(13,26) = ZG040 329 : ZA(13,27) = ZG141 330 : ZA(13,28) =-ZG041 331 : ZA(14,25) = ZG141 332 : ZA(14,26) = ZG041 333 : ZA(14,27) = ZG140+1.D0 334 : ZA(14,28) =-ZG040 335 : C 336 : C -- TRANSFORMATION DE CETTE MATRICE DANS L'ORDRE U R F M P Q 337 : C 338 : C - ICHANG : MATRICE DE TRANSFORMATION DES COLONNES 339 : C 340 : ICHANG( 1) = 1 341 : ICHANG( 2) = 7 342 : ICHANG( 3) = 15 343 : ICHANG( 4) = 21 344 : ICHANG( 5) = 4 345 : ICHANG( 6) = 10 346 : ICHANG( 7) = 18 347 : ICHANG( 8) = 24 348 : ICHANG( 9) = 2 349 : ICHANG(10) = 6 350 : ICHANG(11) = 12 351 : ICHANG(12) = 8 352 : ICHANG(13) = 16 353 : ICHANG(14) = 20 354 : ICHANG(15) = 26 355 : ICHANG(16) = 22 356 : ICHANG(17) = 3 357 : ICHANG(18) = 5 358 : ICHANG(19) = 11 359 : ICHANG(20) = 9 360 : ICHANG(21) = 17 361 : ICHANG(22) = 19 362 : ICHANG(23) = 25 363 : ICHANG(24) = 23 364 : ICHANG(25) = 13 365 : ICHANG(26) = 14 366 : ICHANG(27) = 27 367 : ICHANG(28) = 28 368 : C 369 : DO 10 I = 1,28 370 : K = ICHANG(I) 371 : DO 11 J = 1,14 372 : ZA1(J,K) = ZA(J,I) 373 : 11 CONTINUE 374 : 10 CONTINUE 375 : C 376 : C 377 : C -- CACUL DE ZCI = ZA1*ZR (DANS CET ORDRE CAR LE PRODUIT N'EST 378 : C -- PAS COMMUTATIF 379 : C 380 : C 381 : C - CALCUL DE ZR : 382 : C 383 : DO 26 I = 1 , 28 384 : DO 25 J = 1 , 28 385 : ZER(I,J)= CMPLX(0.D0,0.D0) 386 : 25 CONTINUE 387 : 26 CONTINUE 388 : C 389 : DO 30 L = 1,4 390 : DO 31 I = 1,12 391 : DO 32 J = 1,3 392 : DO 33 K = 1,3 393 : ZER(I,(L-1)*3+J) = ZER(I,(L-1)*3+J) + ZE(I,(L-1)*3+K)*ZR(K,J) 394 : 33 CONTINUE 395 : 32 CONTINUE 396 : 31 CONTINUE 397 : 30 CONTINUE 398 : C 399 : ZER(13,13) = ZE(13,13) 400 : ZER(14,14) = ZE(14,14) 401 : C 402 : DO 34 L = 1,2 403 : DO 35 I = 1,12 404 : DO 36 J = 1,3 405 : DO 37 K = 1,3 406 : ZER(I+14,(L-1)*3+J+14) = ZER(I+14,(L-1)*3+J+14) 407 : * + ZE(I,(L-1)*3+K)*ZR(K,J) 408 : 37 CONTINUE 409 : 36 CONTINUE 410 : 35 CONTINUE 411 : 34 CONTINUE 412 : C 413 : DO 54 L = 3,4 414 : DO 55 I = 1,12 415 : DO 56 J = 1,3 416 : DO 57 K = 1,3 417 : ZER(I+14,(L-1)*3+J+14) = ZER(I+14,(L-1)*3+J+14) 418 : * - ZE(I,(L-1)*3+K)*ZR(K,J) 419 : 57 CONTINUE 420 : 56 CONTINUE 421 : 55 CONTINUE 422 : 54 CONTINUE 423 : C 424 : ZER(27,27) = ZE(13,13) 425 : ZER(28,28) = ZE(14,14) 426 : C 427 : DO 38 I = 1 , 14 428 : DO 39 J = 1 , 28 429 : ZER(I ,J) = ZER(I ,J)/ZD(I ) 430 : ZER(I+14,J) = ZER(I+14,J)/ZD(I ) 431 : 39 CONTINUE 432 : 38 CONTINUE 433 : C 434 : C --------------------------------- CALCUL DE ZA1*ZER 435 : C ----------------- 436 : DO 60 J = 1,28 437 : DO 61 I = 1,14 438 : ZC1(I,J)= CMPLX(0.D0,0.D0) 439 : 61 CONTINUE 440 : 60 CONTINUE 441 : C 442 : DO 40 I = 1,14 443 : DO 41 J = 1,28 444 : DO 42 K = 1,28 445 : ZC1(I,J) = ZC1(I,J) + ZA1(I,K)*ZER(K,J) 446 : 42 CONTINUE 447 : 41 CONTINUE 448 : 40 CONTINUE 449 : C 450 : END 451 : 452 :
© Cast3M 2003 - All rights reserved.
Disclaimer