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