1 : C PB663 SOURCE CHAT 05/01/13 02:10:51 5004 2 : SUBROUTINE PB663(X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NGG,NPG,NOM2) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************ 6 : C 7 : C CALCULE LES FONCTIONS DE FORME D'UN : Iso-Q2 (iso P1/P0 nc) PR18 8 : C 9 : C P4 equation du plan zeta=0 10 : C P5 equation du plan zeta=1 11 : C ^ zeta etc 12 : C | 14 ^ eta |X3 3 13 : C | / \ / |Y3 o 14 : C | / \ / |\ 15 : C | 15_____\13 / | \ 16 : C | / \ /\ / P3 | \ P2 17 : C | / \ / \ / | \ 18 : C |/_____\/____\ / | \ 19 : C 10 11 12 / | \ 20 : C | / 1 o______o 2 21 : C | 9 / |X1 |X2 22 : C | / \ / |Y1 P1 |Y2 23 : C | / \ / 24 : C | 18_____17 25 : C | / \ /\ XI abcisse intersection P2 et ksi 26 : C | / \ / \ YI ordonnée intersection P2 et eta 27 : C |/_____\/____\ 28 : C 7 16 8 29 : C | / fct de formes N1=P2*P5/F1C 30 : C | 5 tel que N1=1 au pt 1 31 : C | / \ 32 : C | / \ 33 : C | 6_____\4 34 : C | / \ /\ 35 : C | / \ / \ 36 : C |/_____\/____\ ____________________>ksi 37 : C 1 2 3 38 : C 39 : C 40 : C************************************************************************ 41 : 42 : CHARACTER*4 NOM2 43 : REAL*8 X(NPG),Y(NPG),Z(NPG) 44 : DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) 45 : DIMENSION FM(MP,NPG),GM(ND,MP,NPG) 46 : DIMENSION U(5),H(5),XA(3),XB(3),XC(3),XD(3) 47 : SAVE XA,XB,XC,XD 48 : DATA XA/3*0.25D0/,XB/0.75D0,0.25D0,0.75D0/ 49 : DATA XC/2*0.75D0,0.25D0/,XD/0.25D0,2*0.75D0/ 50 : 51 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 52 : C 53 : 54 : R2=SQRT(2.D0) 55 : R22=R2*0.5D0 56 : 57 : DO 1 K=1,8 58 : C 59 : IF(K.EQ.1)THEN 60 : X1=0.D0 61 : Y1=0.D0 62 : X2=R22 63 : Y2=0.D0 64 : X3=0.D0 65 : Y3=R22 66 : Z1=0.D0 67 : Z2=0.5D0 68 : XI=R22 69 : YI=R22 70 : ELSEIF(K.EQ.2)THEN 71 : X1=R22 72 : Y1=0.D0 73 : X2=R2 74 : Y2=0.D0 75 : X3=R22 76 : Y3=R22 77 : Z1=0.D0 78 : Z2=0.5D0 79 : XI=R2 80 : YI=R2 81 : ELSEIF(K.EQ.3)THEN 82 : X1=0.D0 83 : Y1=R22 84 : X2=R22 85 : Y2=R22 86 : X3=0.D0 87 : Y3=R2 88 : Z1=0.D0 89 : Z2=0.5D0 90 : XI=R2 91 : YI=R2 92 : ELSEIF(K.EQ.4)THEN 93 : X1=R22 94 : Y1=R22 95 : X2=0.D0 96 : Y2=R22 97 : X3=R22 98 : Y3=0.D0 99 : Z1=0.D0 100 : Z2=0.5D0 101 : XI=R22 102 : YI=R22 103 : ELSEIF(K.EQ.5)THEN 104 : X1=0.D0 105 : Y1=0.D0 106 : X2=R22 107 : Y2=0.D0 108 : X3=0.D0 109 : Y3=R22 110 : Z1=0.5D0 111 : Z2=1.D0 112 : XI=R22 113 : YI=R22 114 : ELSEIF(K.EQ.6)THEN 115 : X1=R22 116 : Y1=0.D0 117 : X2=R2 118 : Y2=0.D0 119 : X3=R22 120 : Y3=R22 121 : Z1=0.5D0 122 : Z2=1.D0 123 : XI=R2 124 : YI=R2 125 : ELSEIF(K.EQ.7)THEN 126 : X1=0.D0 127 : Y1=R22 128 : X2=R22 129 : Y2=R22 130 : X3=0.D0 131 : Y3=R2 132 : Z1=0.5D0 133 : Z2=1.D0 134 : XI=R2 135 : YI=R2 136 : ELSEIF(K.EQ.8)THEN 137 : X1=R22 138 : Y1=R22 139 : X2=0.D0 140 : Y2=R22 141 : X3=R22 142 : Y3=0.D0 143 : Z1=0.5D0 144 : Z2=1.D0 145 : XI=R22 146 : YI=R22 147 : ENDIF 148 : 149 : K0=(K-1)*NGG 149 : CALL KHPRM(X(K0+1),Y(K0+1),Z(K0+1),PG(K0+1),NGG, 150 : &X1,Y1,X2,Y3,Z1,Z2) 151 : C write(6,*)' KHPRM NGG=',ngg 152 : C write(6,1002)(x(ii),ii=1,ngg) 153 : C write(6,1002)(y(ii),ii=1,ngg) 154 : C write(6,1002)(z(ii),ii=1,ngg) 155 : 156 : P1C=Y3-Y1 157 : P2C=X1/XI+Y1/YI-1.D0 158 : P3C=X2-X1 159 : P4C=Z2-Z1 160 : P5C=Z1-Z2 161 : C write(6,*)'PiC ' 162 : C write(6,1002)P1C,P2C,P3C,P4C,P5C 163 : 164 : DO 1 L=1,NGG 165 : P1=Y(L+K0)-Y1 166 : P2=X(L+K0)/XI+Y(L+K0)/YI-1.D0 167 : P3=X(L+K0)-X1 168 : P4=Z(L+K0)-Z1 169 : P5=Z(L+K0)-Z2 170 : 171 : F1C=P2C*P5C 172 : F1=P2*P5/F1C 173 : 174 : F2C=P3C*P5C 175 : F2=P3*P5/F2C 176 : 177 : F3C=P1C*P5C 178 : F3=P1*P5/F3C 179 : 180 : F4C=P2C*P4C 181 : F4=P2*P4/F4C 182 : 183 : F5C=P3C*P4C 184 : F5=P3*P4/F5C 185 : 186 : F6C=P1C*P4C 187 : F6=P1*P4/F6C 188 : 189 : C write(6,*)'FiC' 190 : C write(6,1002)F1C,F2C,F3C,F4C,F5C,F6C 191 : 192 : GX1=P5/XI/F1C 193 : GY1=P5/YI/F1C 194 : GZ1=P2/F1C 195 : 196 : GX2=P5/F2C 197 : GY2=0.D0 198 : GZ2=P3/F2C 199 : 200 : GX3=0.D0 201 : GY3=P5/F3C 202 : GZ3=P1/F3C 203 : 204 : GX4=P4/XI/F4C 205 : GY4=P4/YI/F4C 206 : GZ4=P2/F4C 207 : 208 : GX5=P4/F5C 209 : GY5=0.D0 210 : GZ5=P3/F5C 211 : 212 : GX6=0.D0 213 : GY6=P4/F6C 214 : GZ6=P1/F6C 215 : 216 : C write(6,1002)F1,f2,f3,f4,f5,f6 217 : C 218 : C write(6,1002)gx1,gx2,gx3,gx4,gx5,gx6 219 : C write(6,1002)gy1,gy2,gy3,gy4,gy5,gy6 220 : C write(6,1002)gz1,gz2,gz3,gz4,gz5,gz6 221 : 222 : LL=K0+L 223 : 224 : IF(K.EQ.1)THEN 225 : 226 : FN( 1,LL)=F1 227 : FN( 2,LL)=F2 228 : FN( 6,LL)=F3 229 : FN( 7,LL)=F4 230 : FN(16,LL)=F5 231 : FN(18,LL)=F6 232 : 233 : GR(1, 1,LL)=GX1 234 : GR(1, 2,LL)=GX2 235 : GR(1, 6,LL)=GX3 236 : GR(1, 7,LL)=GX4 237 : GR(1,16,LL)=GX5 238 : GR(1,18,LL)=GX6 239 : 240 : GR(2, 1,LL)=GY1 241 : GR(2, 2,LL)=GY2 242 : GR(2, 6,LL)=GY3 243 : GR(2, 7,LL)=GY4 244 : GR(2,16,LL)=GY5 245 : GR(2,18,LL)=GY6 246 : 247 : GR(3, 1,LL)=GZ1 248 : GR(3, 2,LL)=GZ2 249 : GR(3, 6,LL)=GZ3 250 : GR(3, 7,LL)=GZ4 251 : GR(3,16,LL)=GZ5 252 : GR(3,18,LL)=GZ6 253 : 254 : ELSEIF(K.EQ.2)THEN 255 : 256 : FN( 2,LL)=F1 257 : FN( 3,LL)=F2 258 : FN( 4,LL)=F3 259 : FN(16,LL)=F4 260 : FN( 8,LL)=F5 261 : FN(17,LL)=F6 262 : 263 : GR(1, 2,LL)=GX1 264 : GR(1, 3,LL)=GX2 265 : GR(1, 4,LL)=GX3 266 : GR(1,16,LL)=GX4 267 : GR(1, 8,LL)=GX5 268 : GR(1,17,LL)=GX6 269 : 270 : GR(2, 2,LL)=GY1 271 : GR(2, 3,LL)=GY2 272 : GR(2, 4,LL)=GY3 273 : GR(2,16,LL)=GY4 274 : GR(2, 8,LL)=GY5 275 : GR(2,17,LL)=GY6 276 : 277 : GR(3, 2,LL)=GZ1 278 : GR(3, 3,LL)=GZ2 279 : GR(3, 4,LL)=GZ3 280 : GR(3,16,LL)=GZ4 281 : GR(3, 8,LL)=GZ5 282 : GR(3,17,LL)=GZ6 283 : 284 : ELSEIF(K.EQ.3)THEN 285 : 286 : FN( 6,LL)=F1 287 : FN( 4,LL)=F2 288 : FN( 5,LL)=F3 289 : FN(18,LL)=F4 290 : FN(17,LL)=F5 291 : FN( 9,LL)=F6 292 : 293 : GR(1, 6,LL)=GX1 294 : GR(1, 4,LL)=GX2 295 : GR(1, 5,LL)=GX3 296 : GR(1,18,LL)=GX4 297 : GR(1,17,LL)=GX5 298 : GR(1, 9,LL)=GX6 299 : 300 : GR(2, 6,LL)=GY1 301 : GR(2, 4,LL)=GY2 302 : GR(2, 5,LL)=GY3 303 : GR(2,18,LL)=GY4 304 : GR(2,17,LL)=GY5 305 : GR(2, 9,LL)=GY6 306 : 307 : GR(3, 6,LL)=GZ1 308 : GR(3, 4,LL)=GZ2 309 : GR(3, 5,LL)=GZ3 310 : GR(3,18,LL)=GZ4 311 : GR(3,17,LL)=GZ5 312 : GR(3, 9,LL)=GZ6 313 : 314 : ELSEIF(K.EQ.4)THEN 315 : 316 : FN( 4,LL)=F1 317 : FN( 6,LL)=F2 318 : FN( 2,LL)=F3 319 : FN(17,LL)=F4 320 : FN(18,LL)=F5 321 : FN(16,LL)=F6 322 : 323 : GR(1, 4,LL)=GX1 324 : GR(1, 6,LL)=GX2 325 : GR(1, 2,LL)=GX3 326 : GR(1,17,LL)=GX4 327 : GR(1,18,LL)=GX5 328 : GR(1,16,LL)=GX6 329 : 330 : GR(2, 4,LL)=GY1 331 : GR(2, 6,LL)=GY2 332 : GR(2, 2,LL)=GY3 333 : GR(2,17,LL)=GY4 334 : GR(2,18,LL)=GY5 335 : GR(2,16,LL)=GY6 336 : 337 : GR(3, 4,LL)=GZ1 338 : GR(3, 6,LL)=GZ2 339 : GR(3, 2,LL)=GZ3 340 : GR(3,17,LL)=GZ4 341 : GR(3,18,LL)=GZ5 342 : GR(3,16,LL)=GZ6 343 : 344 : ELSEIF(K.EQ.5)THEN 345 : 346 : FN( 7,LL)=F1 347 : FN(16,LL)=F2 348 : FN(18,LL)=F3 349 : FN(10,LL)=F4 350 : FN(11,LL)=F5 351 : FN(15,LL)=F6 352 : 353 : GR(1, 7,LL)=GX1 354 : GR(1,16,LL)=GX2 355 : GR(1,18,LL)=GX3 356 : GR(1,10,LL)=GX4 357 : GR(1,11,LL)=GX5 358 : GR(1,15,LL)=GX6 359 : 360 : GR(2, 7,LL)=GY1 361 : GR(2,16,LL)=GY2 362 : GR(2,18,LL)=GY3 363 : GR(2,10,LL)=GY4 364 : GR(2,11,LL)=GY5 365 : GR(2,15,LL)=GY6 366 : 367 : GR(3, 7,LL)=GZ1 368 : GR(3,16,LL)=GZ2 369 : GR(3,18,LL)=GZ3 370 : GR(3,10,LL)=GZ4 371 : GR(3,11,LL)=GZ5 372 : GR(3,15,LL)=GZ6 373 : 374 : ELSEIF(K.EQ.6)THEN 375 : 376 : FN(16,LL)=F1 377 : FN( 8,LL)=F2 378 : FN(17,LL)=F3 379 : FN(11,LL)=F4 380 : FN(12,LL)=F5 381 : FN(13,LL)=F6 382 : 383 : GR(1,16,LL)=GX1 384 : GR(1, 8,LL)=GX2 385 : GR(1,17,LL)=GX3 386 : GR(1,11,LL)=GX4 387 : GR(1,12,LL)=GX5 388 : GR(1,13,LL)=GX6 389 : 390 : GR(2,16,LL)=GY1 391 : GR(2, 8,LL)=GY2 392 : GR(2,17,LL)=GY3 393 : GR(2,11,LL)=GY4 394 : GR(2,12,LL)=GY5 395 : GR(2,13,LL)=GY6 396 : 397 : GR(3,16,LL)=GZ1 398 : GR(3, 8,LL)=GZ2 399 : GR(3,17,LL)=GZ3 400 : GR(3,11,LL)=GZ4 401 : GR(3,12,LL)=GZ5 402 : GR(3,13,LL)=GZ6 403 : 404 : ELSEIF(K.EQ.7)THEN 405 : 406 : FN(18,LL)=F1 407 : FN(17,LL)=F2 408 : FN( 9,LL)=F3 409 : FN(15,LL)=F4 410 : FN(13,LL)=F5 411 : FN(14,LL)=F6 412 : 413 : GR(1,18,LL)=GX1 414 : GR(1,17,LL)=GX2 415 : GR(1, 9,LL)=GX3 416 : GR(1,15,LL)=GX4 417 : GR(1,13,LL)=GX5 418 : GR(1,14,LL)=GX6 419 : 420 : GR(2,18,LL)=GY1 421 : GR(2,17,LL)=GY2 422 : GR(2, 9,LL)=GY3 423 : GR(2,15,LL)=GY4 424 : GR(2,13,LL)=GY5 425 : GR(2,14,LL)=GY6 426 : 427 : GR(3,18,LL)=GZ1 428 : GR(3,17,LL)=GZ2 429 : GR(3, 9,LL)=GZ3 430 : GR(3,15,LL)=GZ4 431 : GR(3,13,LL)=GZ5 432 : GR(3,14,LL)=GZ6 433 : 434 : ELSEIF(K.EQ.8)THEN 435 : 436 : FN(17,LL)=F1 437 : FN(18,LL)=F2 438 : FN(16,LL)=F3 439 : FN(13,LL)=F4 440 : FN(15,LL)=F5 441 : FN(11,LL)=F6 442 : 443 : GR(1,17,LL)=GX1 444 : GR(1,18,LL)=GX2 445 : GR(1,16,LL)=GX3 446 : GR(1,13,LL)=GX4 447 : GR(1,15,LL)=GX5 448 : GR(1,11,LL)=GX6 449 : 450 : GR(2,17,LL)=GY1 451 : GR(2,18,LL)=GY2 452 : GR(2,16,LL)=GY3 453 : GR(2,13,LL)=GY4 454 : GR(2,15,LL)=GY5 455 : GR(2,11,LL)=GY6 456 : 457 : GR(3,17,LL)=GZ1 458 : GR(3,18,LL)=GZ2 459 : GR(3,16,LL)=GZ3 460 : GR(3,13,LL)=GZ4 461 : GR(3,15,LL)=GZ5 462 : GR(3,11,LL)=GZ6 463 : 464 : ENDIF 465 : 466 : 1 CONTINUE 467 : 468 : IF(NOM2.EQ.'MCP0')THEN 469 : DO 2 L=1,NPG 470 : FM(1,L)=1.D0 471 : GM(1,1,L)=0.D0 472 : GM(2,1,L)=0.D0 473 : GM(3,1,L)=0.D0 474 : 2 CONTINUE 475 : ELSEIF(NOM2.EQ.'MCP1')THEN 476 : 477 : DO 3 LL=1,(2*NGG) 478 : FM(1,LL)=0.25D0 479 : FM(2,LL+2*NGG)=0.25D0 480 : FM(3,LL+4*NGG)=0.25D0 481 : FM(4,LL+6*NGG)=0.25D0 482 : 3 CONTINUE 482 : CALL INITD(GM,(12*NPG),0.D0) 483 : ELSEIF(NOM2.EQ.'MCF1')THEN 484 : DO 4 L=1,NPG 485 : FM(1,L)=(X(L)+Y(L)-R2)*(Z(L)-1.D0)/R2 486 : FM(2,L)=-X(L)*(Z(L)-1.D0)/R2 487 : FM(3,L)=-Y(L)*(Z(L)-1.D0)/R2 488 : FM(4,L)=-(X(L)+Y(L)-R2)*Z(L)/R2 489 : FM(5,L)=X(L)*Z(L)/R2 490 : FM(6,L)=Y(L)*Z(L)/R2 491 : 492 : GM(1,1,L)=(Z(L)-1.D0)/R2 493 : GM(2,1,L)=(Z(L)-1.D0)/R2 494 : GM(3,1,L)=(X(L)+Y(L)-R2)/R2 495 : GM(1,2,L)=-(Z(L)-1.D0)/R2 496 : GM(2,2,L)=0.D0 497 : GM(3,2,L)=-X(L)/R2 498 : GM(1,3,L)=0.D0 499 : GM(2,3,L)=-(Z(L)-1.D0)/R2 500 : GM(3,3,L)=-Y(L)/R2 501 : GM(1,4,L)=-Z(L)/R2 502 : GM(2,4,L)=-Z(L)/R2 503 : GM(3,4,L)=-(X(L)+Y(L)-R2)/R2 504 : GM(1,5,L)=Z(L)/R2 505 : GM(2,5,L)=0.D0 506 : GM(3,5,L)=X(L)/R2 507 : GM(1,6,L)=0.D0 508 : GM(2,6,L)=Z(L)/R2 509 : GM(3,6,L)=Y(L)/R2 510 : C 511 : 4 CONTINUE 512 : ENDIF 513 : C write(6,*)' FM' 514 : C write(6,1002)((FM(II,JJ),II=1,MP),JJ=1,NPG) 515 : 516 : C write(6,*)' VERIF ,PG=' 517 : C write(6,1002)(pg(ii),ii=1,npg) 518 : C do 75 ll=1,npg 519 : C write(6,*)' VERIF ,fn,gr ll=',ll 520 : C write(6,1002)(fn(ii,ll),ii=1,np) 521 : C write(6,1002)(gr(1,ii,ll),ii=1,np) 522 : C write(6,1002)(gr(2,ii,ll),ii=1,np) 523 : C write(6,1002)(gr(3,ii,ll),ii=1,np) 524 : C75 continue 525 : 526 : C write(6,*)' VERIF ,NPG,NP,NGG=',NPG,NP,NGG 527 : C UPG=0.D0 528 : C do 72 L=1,NPG 529 : C UPG=UPG+PG(L) 530 : C UF=0.D0 531 : C UG1=0.D0 532 : C UG2=0.D0 533 : C UG3=0.D0 534 : C DO 71 I=1,NP 535 : C UF=UF+FN(I,L) 536 : C UG1=UG1+GR(1,I,L) 537 : C UG2=UG2+GR(2,I,L) 538 : C UG3=UG3+GR(3,I,L) 539 : C71 CONTINUE 540 : C WRITE(6,*)' VERIF L=',L,UF,UG1,UG2,UG3 541 : C72 CONTINUE 542 : 543 : C WRITE(6,*)' VERIF PG=',UPG 544 : C WRITE(6,101) 545 : 546 : C IF(MP.EQ.1)THEN 547 : C DO 4 L=1,NPG 548 : C FM(1,L)=1.D0 549 : C GM(1,1,L)=0.D0 550 : C GM(2,1,L)=0.D0 551 : C GM(3,1,L)=0.D0 552 : C4 CONTINUE 553 : C ELSEIF(MP.EQ.4)THEN 554 : C XX(1)=X(L) 555 : C XX(2)=Y(L) 556 : C XX(3)=Z(L) 557 : C FM(1,L)=EQPL3P(XX,XB,XC,XD)/ 558 : C & EQPL3P(XA,XB,XC,XD) 559 : C FM(2,L)=EQPL3P(XX,XA,XC,XD)/ 560 : C & EQPL3P(XB,XA,XC,XD) 561 : C FM(3,L)=EQPL3P(XX,XA,XB,XD)/ 562 : C & EQPL3P(XC,XA,XB,XD) 563 : C 564 : C FM(4,L)=EQPL3P(XX,XA,XB,XC)/ 565 : C & EQPL3P(XD,XA,XB,XC) 565 : C CALL INITD(GM,(12*NPG),0.D0) 566 : C3 CONTINUE 567 : C ENDIF 568 : 569 : C WRITE(6,101) 570 : 571 : RETURN 572 : 1002 FORMAT(10(1X,1PD11.4)) 573 : 1001 FORMAT(20(1X,I5)) 574 : 101 FORMAT(1X,'... SUBPB663 ... FM,GM,FN,GR ',9(10H..........)/) 575 : C 576 : END 577 : 578 : 579 : 580 : 581 :
© Cast3M 2003 - All rights reserved.
Disclaimer