1 : C FACAXI SOURCE CHAT 11/03/16 21:21:55 6902 2 : SUBROUTINE FACAXI (MYMOD,INFOEL,NPAX,NGAX,KACHE,INOR,ICHFAC 3 : & ,EXTINC) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : C---------------------------------------------------------------------- 7 : C SP appele par FFORME 8 : C 9 : C Calcul des facteurs de forme en axisymetrique 10 : C entée: 11 : C MYMOD : pointeur de l'objet modèle 12 : C INFOEL : utile pour les coques ou quadratiques 13 : C NPAX : nombre de points d integration (disc.reguliere) 14 : C NGAX : nombre de points de Gauss 15 : C KACHE : 0 si option convexe, sinon option cache 16 : C INOR : 1 si normalisation, 0 sinon 17 : C EXTINC : coefficient d'extinction (si cavite absorbante) 18 : C sortie: 19 : C ICHFAC : pointeur sur l'objet MCHAML resultat 20 : C 21 : C---------------------------------------------------------------------- 22 : C traitement des coques par dedoublement des elements a partir 23 : C de la normale 24 : C -> 25 : C A (inverse) = A - e*n (e=1e-3) 26 : C cas des boucles 1 sur k1 et 2 sur k2 27 : C mais pas de la boucle 3 obstructeurs 28 : C 29 : C bcl face k1 30 : C ** face k1 ** 31 : C bcl face k2 32 : C .. face k2 .. 33 : C bcl 3 obstructeurs 34 : C .. si coq: inverse face k2 .. 35 : C les obstructeurs sont les memes que pour k2 36 : C fin bcl face k2 37 : C 38 : C ** si coq : inverse face k1 ** 39 : C bcl face k2 40 : C .. face k2 .. 41 : C bcl 3 obstructeurs 42 : C .. si coq: inverse face k2 .. 43 : C les obstructeurs sont les memes que pour k2 44 : C fin bcl face k2 45 : C fin bcl face k1 46 : C---------------------------------------------------------------------- 47 : C 48 : LOGICAL ICOQ 49 : C 50 : -INC CCREEL 51 : -INC SMELEME 52 : -INC SMMODEL 53 : -INC CCOPTIO 54 : -INC SMCOORD 55 : POINTEUR MYMOD.MMODEL 56 : 57 : C---------------------------------------------------------------------- 58 : SEGMENT , INFOEL 59 : LOGICAL KCOQ(N1),KQUAD(N1) 60 : ENDSEGMENT 61 : C---------------------------------------------------------------------- 62 : C FACTEURS DE FORME 63 : C NNBEL1 = NOMBRE DE LIGNES + 1 64 : C NBEL2 = NOMBRE DE COLONNES 65 : C LFACT(NNBEL1) POINTE SUR LE TABLEAU DES SURFACES 66 : C 67 : SEGMENT IFACFO 68 : INTEGER LFACT(NNBEL1) 69 : ENDSEGMENT 70 : SEGMENT LFAC 71 : REAL*8 FACT(NBEL2) 72 : ENDSEGMENT 73 : POINTEUR PSUR.LFAC, PLIG.LFAC 74 : POINTEUR PCOL.LFAC 75 : C---------------------------------------------------------------------- 76 : C coordonnees des obstructeurs 77 : SEGMENT SFOBS 78 : REAL*8 OBS(2,NTOBS) 79 : ENDSEGMENT 80 : C---------------------------------------------------------------------- 81 : SEGMENT STRAV 82 : REAL*8 A1(NA,NA),A2(NA,NA),A3(NA,NA),AA2(NA,NA) 83 : REAL*8 U1(NA1),U2(NA1),UU2(NA1) 84 : ENDSEGMENT 85 : C---------------------------------------------------------------------- 86 : C 87 : EPS = 1D-5 88 : 89 : KIMP = IIMPI 90 : NES = IDIM 91 : 92 : C... critere de dedoublement des coques 93 : ECOQ=1.D-3 94 : IF (INFOEL.EQ.0) THEN 95 : ICOQ = .FALSE. 96 : ELSE 97 : ICOQ = .TRUE. 98 : SEGACT INFOEL 99 : ENDIF 100 : C... quadratique 101 : NSPA1=1 102 : NSPA2=1 103 : NSPA3=1 104 : NS=2 105 : 106 : RAD = 0 107 : C 108 : C Calcul du nombre de faces NFACE et du rayon RAD 109 : SEGACT MYMOD 110 : NTYP = MYMOD.KMODEL(/1) 111 : NFACE = 0 112 : DO 10 ITYP=1,NTYP 113 : IMODEL = MYMOD.KMODEL(ITYP) 114 : SEGACT IMODEL 115 : IPT1 = IMAMOD 116 : SEGDES IMODEL 117 : SEGACT IPT1 118 : NEL = IPT1.NUM(/2) 119 : NSGEO = IPT1.NUM(/1) 120 : C Recherche du max sur Ox 121 : DO 5 IEL = 1,NEL 122 : DO 6 IPT = 1,NSGEO 123 : IREF = (IDIM+1)*(IPT1.NUM(IPT,IEL)-1) 124 : VALX = XCOOR(IREF+1) 125 : IF (VALX.GT.RAD) RAD = VALX 126 : 6 CONTINUE 127 : 5 CONTINUE 128 : IF (ICOQ.AND.KCOQ(ITYP)) THEN 129 : NFACE = NFACE + 2 * NEL 130 : ELSE 131 : NFACE = NFACE + NEL 132 : ENDIF 133 : 134 : 10 CONTINUE 135 : 136 : C 137 : 138 : IF (KIMP.GE.3) THEN 139 : WRITE( 6,*) ' DIMENSIONNEMENT : RAD =',RAD 140 : WRITE( 6,*) ' NOMBRE TOTAL DE FACES ',NFACE 141 : ENDIF 142 : C 143 : C>>> INITIALISATION SFOBS et STRAV 144 : C ----------------------------- 145 : IF(KACHE.NE.0) THEN 146 : NTOBS = 2*NFACE 147 : SEGINI SFOBS 148 : ENDIF 149 : NA=2 150 : NA1=3 151 : SEGINI STRAV 152 : 153 : C 154 : C>>> INITIALISATION OBJET FACFOR 155 : C --------------------------- 156 : C 157 : NNBEL1=NFACE+1 158 : NBEL2=NFACE 159 : SEGINI IFACFO 160 : DO 900 LS=1,NNBEL1 161 : SEGINI PLIG 162 : LFACT(LS)=PLIG 163 : SEGACT PLIG*MOD 164 : 900 CONTINUE 165 : PSUR = LFACT(NNBEL1) 166 : SEGACT PSUR*MOD 167 : 168 : C ------------------------------------------------------------- 169 : C>>> CALCUL 170 : 171 : C 172 : C>> BOUCLE FACE 1 173 : C 174 : 175 : NELT1= 0 176 : DO 100 ITYP = 1,NTYP 177 : IMODEL = MYMOD.KMODEL(ITYP) 178 : SEGACT IMODEL 179 : IPT1 = IMAMOD 180 : SEGDES IMODEL 181 : IF(KIMP.GE.4) THEN 182 : WRITE (6,*) 'Maillage :',IPT1 183 : ENDIF 184 : 185 : NSGEO1 = IPT1.NUM(/1) 186 : NSPA1=1 187 : IF(ICOQ.AND.KQUAD(ITYP)) NSPA1=2 188 : NEL1 = IPT1.NUM(/2) 189 : DO 110 I1 = 1,NEL1 190 : 191 : IF (ICOQ.AND.KCOQ(ITYP)) THEN 192 : K1 = (2*I1-1) + NELT1 193 : ELSE 194 : K1 = I1+ NELT1 195 : ENDIF 196 : 197 : C*** traitement de la face K1 *************************************** 198 : 199 : PLIG=LFACT(K1) 200 : C 201 : DO 111 IS = 1,NSGEO1,NSPA1 202 : LS=IS 203 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 204 : IREF = (IDIM+1)*(IPT1.NUM(IS,I1)-1) 205 : IF(KIMP.GE.4) THEN 206 : WRITE (6,*) 'Noeud numéro',IPT1.NUM(IS,I1),'ref :' 207 : $ ,IREF 208 : ENDIF 209 : DO 112 K = 1,NES 210 : C** A1(K,IS) = XCOOR(IREF+K)/RAD 211 : A1(K,LS) = XCOOR(IREF+K)/RAD 212 : 112 CONTINUE 213 : 111 CONTINUE 213 : CALL KNORM2(NES,A1,NS,U1,S1) 214 : S1=XPI*S1*DABS(A1(1,1)+A1(1,2)) 215 : PSUR.FACT(K1)=S1 216 : ZMIN1=DMIN1(A1(2,1),A1(2,2)) 217 : ZMAX1=DMAX1(A1(2,1),A1(2,2)) 218 : RMAX1=DMAX1(A1(1,1),A1(1,2)) 219 : 220 : 221 : C>> BOUCLE FACE 2 222 : C ------------------------------------------------------------- 223 : 224 : NELT2= 0 225 : DO 200 JTYP = 1,NTYP 226 : IMODEL = MYMOD.KMODEL(JTYP) 227 : SEGACT IMODEL 228 : IPT2 = IMAMOD 229 : SEGDES IMODEL 230 : 231 : NSGEO2 = IPT2.NUM(/1) 232 : NSPA2=1 233 : IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2 234 : NEL2 = IPT2.NUM(/2) 235 : DO 210 I2 = 1,NEL2 236 : 237 : IF (ICOQ.AND.KCOQ(JTYP))THEN 238 : K2 = 2*I2-1 + NELT2 239 : ELSE 240 : K2 = I2 + NELT2 241 : ENDIF 242 : 243 : C... face K2 ..................................................... 244 : 245 : KVU=1 246 : FF=0.D0 247 : 248 : C.. UTILISATION DE LA RECIPROCITE 249 : C 250 : IF(K1.GT.K2) THEN 251 : S2=PSUR.FACT(K2) 252 : PCOL=LFACT(K2) 253 : 254 : PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1) 255 : 256 : ELSE 257 : C.. TEST K1=K2 ET VISIBILITE A PRIORI 258 : 259 : C------------------------------------------------------------------ 260 : 261 : C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0 262 : C 263 : IF(K1.EQ.K2) THEN 264 : IF(KIMP.GE.4) THEN 265 : WRITE(6,*) ' A1 ',A1(1,1),A1(2,1) 266 : WRITE(6,*) ' A1 ',A1(1,2),A1(2,2) 267 : WRITE(6,*) ' U1 ',U1(1),U1(2) 268 : ENDIF 269 : IF (DABS(U1(1)).LT.EPS) KVU=0 270 : IF (U1(1).GE.EPS) KVU=0 271 : 272 : IF(KVU.NE.0) THEN 273 : DO 214 K = 1,NES 274 : A2(K,1) = A1(K,1) 275 : A2(K,2) = A1(K,2) 276 : 214 CONTINUE 277 : ENDIF 278 : 279 : C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE 280 : C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z 281 : C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z) 282 : 283 : ELSE 284 : DO 211 IS = 1,NSGEO2,NSPA2 285 : LS=IS 286 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 287 : IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1) 288 : DO 212 K = 1,NES 289 : C** A2(K,IS) = XCOOR(IREF+K)/RAD 290 : A2(K,LS) = XCOOR(IREF+K)/RAD 291 : 212 CONTINUE 292 : 211 CONTINUE 292 : CALL KNORM2(NES,A2,NS,U2,S2) 293 : S2=XPI*S2*DABS(A2(1,1)+A2(1,2)) 294 : 295 : C.. calcul du symetrique 296 : AA2(1,1)= -A2(1,2) 297 : AA2(2,1)= A2(2,2) 298 : AA2(1,2)= -A2(1,1) 299 : AA2(2,2)= A2(2,1) 299 : CALL KNORM2(NES,AA2,NS,UU2,SS2) 300 : 301 : C.. visibilite a priori 301 : CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU) 301 : CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS) 302 : 303 : IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2 304 : $ ,IVU,IVUS 305 : IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0 306 : 307 : ENDIF 308 : C --------------------------------------------------------------- 309 : C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI 310 : IF(KIMP.GE.4.AND.KVU.NE.0) THEN 311 : WRITE(6,*) ' FACES VISIBLES ',K1,K2 312 : ENDIF 313 : 314 : C>> CALCUL 315 : C... option convexe 316 : IF(KACHE.EQ.0) THEN 317 : IF(KVU.NE.0) THEN 318 : 318 : CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD) 319 : IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2 320 : $ ,FF/S1 321 : 322 : ENDIF 323 : C.. option cache 324 : ELSE 325 : IF(KVU.NE.0) THEN 326 : 327 : C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.-------------------------- 328 : 329 : NOBS=0 330 : 331 : C>> BOUCLE FACE 3--------------------------------------------------- 332 : NELT3= 0 333 : DO 300 KTYP = 1,NTYP 334 : IMODEL = MYMOD.KMODEL(KTYP) 335 : SEGACT IMODEL 336 : IPT3 = IMAMOD 337 : SEGDES IMODEL 338 : 339 : NSGEO3 = IPT3.NUM(/1) 340 : NSPA3=1 341 : IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2 342 : NEL3 = IPT3.NUM(/2) 343 : DO 310 I3 = 1,NEL3 344 : K3 = I3+ NELT3 345 : 346 : IF(K3.NE.K1.AND.K3.NE.K2) THEN 347 : 348 : DO 311 IS = 1,NSGEO3,NSPA3 349 : LS=IS 350 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1 351 : $ )/2 352 : IREF = (IDIM+1)*(IPT3.NUM(IS,I3) 353 : $ -1) 354 : DO 312 K = 1,NES 355 : C** A3(K,IS) = XCOOR(IREF+K)/RAD 356 : A3(K,LS) = XCOOR(IREF+K)/RAD 357 : 312 CONTINUE 358 : 311 CONTINUE 359 : 360 : ZMIN2=DMIN1(A2(2,1),A2(2,2)) 361 : ZMAX2=DMAX1(A2(2,1),A2(2,2)) 362 : RMAX2=DMAX1(A2(1,1),A2(1,2)) 363 : ZTMIN=DMIN1(ZMIN1,ZMIN2) 364 : ZTMAX=DMAX1(ZMAX1,ZMAX2) 365 : RMAX =DMAX1(RMAX1,RMAX2) 366 : 367 : IF(DMAX1(A3(2,1),A3(2,2)).LE.ZTMIN) 368 : $ THEN 369 : KOBS=0 370 : ELSEIF(DMIN1(A3(2,1),A3(2,2)).GE 371 : $ .ZTMAX) THEN 372 : KOBS=0 373 : ELSEIF(DMIN1(A3(1,1),A3(1,2)).GE 374 : $ .RMAX) THEN 375 : KOBS=0 376 : ELSE 377 : KOBS=1 378 : OBS(1,2*NOBS+1)=A3(1,1) 379 : OBS(2,2*NOBS+1)=A3(2,1) 380 : OBS(1,2*NOBS+2)=A3(1,2) 381 : OBS(2,2*NOBS+2)=A3(2,2) 382 : NOBS=NOBS+1 383 : 384 : ENDIF 385 : ENDIF 386 : 387 : 310 CONTINUE 388 : C 389 : NELT3=NELT3+NEL3 390 : C 391 : 300 CONTINUE 392 : 393 : C<< FIN RECHERCHE OBSTRUCTEURS-------------------------------------- 394 : 395 : IF(KIMP.GE.4) THEN 396 : WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS ' 397 : $ ,NOBS 398 : ENDIF 399 : 399 : CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX,FF 400 : $ ,KIMP,EXTINC,RAD) 401 : IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF/S1 402 : 403 : ENDIF 404 : 405 : ENDIF 406 : C<< FIN CALCUL 407 : C --------------------------------------------------------------- 408 : 409 : C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG 410 : PLIG.FACT(K2) = FF/S1 411 : 412 : ENDIF 413 : C... fin face K2 .................................................. 414 : 415 : 416 : C... face inverse de K2 (cas des coques) .......................... 417 : 418 : IF (ICOQ.AND.KCOQ(JTYP)) THEN 419 : 420 : K2=K2 + 1 421 : 422 : KVU=1 423 : FF=0.D0 424 : 425 : C UTILISATION DE LA RECIPROCITE 426 : C 427 : IF(K1.GT.K2) THEN 428 : S2=PSUR.FACT(K2) 429 : PCOL=LFACT(K2) 430 : 431 : PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1) 432 : 433 : ELSE 434 : 435 : C>> TEST K1=K2 ET VISIBILITE A PRIORI 436 : 437 : C------------------------------------------------------------------ 438 : 439 : C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0 440 : C 441 : IF(K1.EQ.K2) THEN 442 : IF(KIMP.GE.4) THEN 443 : WRITE(6,*) ' A1 ',A1(1,1),A1(2,1) 444 : WRITE(6,*) ' A1 ',A1(1,2),A1(2,2) 445 : WRITE(6,*) ' U1 ',U1(1),U1(2) 446 : ENDIF 447 : IF (DABS(U1(1)).LT.EPS) KVU=0 448 : IF (U1(1).GE.EPS) KVU=0 449 : 450 : IF(KVU.NE.0) THEN 451 : DO 514 K = 1,NES 452 : A2(K,1) = A1(K,1) 453 : A2(K,2) = A1(K,2) 454 : 514 CONTINUE 455 : ENDIF 456 : 457 : C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE 458 : C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z 459 : C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z) 460 : 461 : 462 : ELSE 463 : 464 : DO 512 K = 1,NES 465 : U2(K)=-U2(K) 466 : XX1= A2(K,1) 467 : A2(K,1) = A2(K,2) + ECOQ*U2(K) 468 : A2(K,2) = XX1 + ECOQ*U2(K) 469 : 512 CONTINUE 470 : 471 : C... calcul du symetrique 472 : AA2(1,1)= -A2(1,2) 473 : AA2(2,1)= A2(2,2) 474 : AA2(1,2)= -A2(1,1) 475 : AA2(2,2)= A2(2,1) 475 : CALL KNORM2(NES,AA2,NS,UU2,SS2) 476 : 477 : C... visibilité a priori 477 : CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU) 477 : CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS) 478 : 479 : IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1 480 : $ ,K2,IVU,IVUS 481 : IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0 482 : 483 : ENDIF 484 : C --------------------------------------------------------------- 485 : C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI 486 : IF(KIMP.GE.4.AND.KVU.NE.0) THEN 487 : WRITE(6,*) ' FACES VISIBLES ',K1,K2 488 : ENDIF 489 : 490 : C>> CALCUL -------------------------------------------------------- 491 : C.. option convexe 492 : IF(KACHE.EQ.0) THEN 493 : IF(KVU.NE.0) THEN 494 : 494 : CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC 495 : $ ,RAD) 496 : IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1 497 : $ ,K2,FF/S1 498 : 499 : ENDIF 500 : C.. option cache 501 : ELSE 502 : IF(KVU.NE.0) THEN 503 : 504 : C>> les obstructeurs potentiels sont les memes que pour la face K2 505 : 506 : IF(KIMP.GE.4) THEN 507 : WRITE(6,*) ' FACES K1 K2 ',K1,K2 508 : $ ,' NOBS ',NOBS 509 : ENDIF 510 : 510 : CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX 511 : $ ,FF,KIMP,EXTINC,RAD) 512 : 513 : 514 : IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF 515 : $ /S1 516 : 517 : ENDIF 518 : 519 : ENDIF 520 : C<< FIN CALCUL ------------------------------------------------------ 521 : 522 : C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG 523 : PLIG.FACT(K2) = FF/S1 524 : 525 : ENDIF 526 : 527 : ENDIF 528 : 529 : C... fin face inverse de K2 (cas des coques) ......................... 530 : 531 : 210 CONTINUE 532 : 533 : IF (ICOQ.AND.KCOQ(JTYP)) THEN 534 : NELT2 = NELT2 + 2 * NEL2 535 : ELSE 536 : NELT2 = NELT2 + NEL2 537 : ENDIF 538 : 539 : 200 CONTINUE 540 : 541 : C ---FIN BOUCLE FACE 2----------------------------------------- 542 : 543 : 544 : 545 : 546 : C*** fin traitement de la face K1 *********************************** 547 : 548 : C*** traitement de la face inverse de K1***************************** 549 : IF (ICOQ.AND.KCOQ(ITYP)) THEN 550 : 551 : K1=K1+1 552 : PLIG=LFACT(K1) 553 : 554 : DO 122 K = 1,NES 555 : U1(K) =-U1(K) 556 : XX1=A1(K,1) 557 : A1(K,1) = A1(K,2) + ECOQ * U1(K) 558 : A1(K,2) = XX1 + ECOQ * U1(K) 559 : 122 CONTINUE 560 : 561 : PSUR.FACT(K1)=S1 562 : 563 : C>> BOUCLE FACE 2 564 : C ------------------------------------------------------------- 565 : 566 : NELT2= 0 567 : DO 2000 JTYP = 1,NTYP 568 : IMODEL = MYMOD.KMODEL(JTYP) 569 : SEGACT IMODEL 570 : IPT2 = IMAMOD 571 : SEGDES IMODEL 572 : 573 : NSGEO2 = IPT2.NUM(/1) 574 : NSPA2=1 575 : IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2 576 : NEL2 = IPT2.NUM(/2) 577 : DO 2100 I2 = 1,NEL2 578 : 579 : IF (ICOQ.AND.KCOQ(JTYP))THEN 580 : K2 = 2*I2-1 + NELT2 581 : ELSE 582 : K2 = I2 + NELT2 583 : ENDIF 584 : 585 : C... face K2 ..................................................... 586 : 587 : KVU=1 588 : FF=0.D0 589 : 590 : C.. UTILISATION DE LA RECIPROCITE 591 : C 592 : IF(K1.GT.K2) THEN 593 : S2=PSUR.FACT(K2) 594 : PCOL=LFACT(K2) 595 : 596 : PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1) 597 : 598 : ELSE 599 : C.. TEST K1=K2 ET VISIBILITE A PRIORI 600 : 601 : C------------------------------------------------------------------ 602 : 603 : C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0 604 : C 605 : IF(K1.EQ.K2) THEN 606 : IF(KIMP.GE.4) THEN 607 : WRITE(6,*) ' A1 ',A1(1,1),A1(2,1) 608 : WRITE(6,*) ' A1 ',A1(1,2),A1(2,2) 609 : WRITE(6,*) ' U1 ',U1(1),U1(2) 610 : ENDIF 611 : IF (DABS(U1(1)).LT.EPS) KVU=0 612 : IF (U1(1).GE.EPS) KVU=0 613 : 614 : IF(KVU.NE.0) THEN 615 : DO 2140 K = 1,NES 616 : A2(K,1) = A1(K,1) 617 : A2(K,2) = A1(K,2) 618 : 2140 CONTINUE 619 : ENDIF 620 : 621 : C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE 622 : C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z 623 : C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z) 624 : 625 : ELSE 626 : DO 2110 IS = 1,NSGEO2,NSPA2 627 : LS=IS 628 : IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2 629 : IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1) 630 : DO 2120 K = 1,NES 631 : C** A2(K,IS) = XCOOR(IREF+K)/RAD 632 : A2(K,LS) = XCOOR(IREF+K)/RAD 633 : 2120 CONTINUE 634 : 2110 CONTINUE 634 : CALL KNORM2(NES,A2,NS,U2,S2) 635 : S2=XPI*S2*DABS(A2(1,1)+A2(1,2)) 636 : 637 : C.. calcul du symetrique 638 : AA2(1,1)= -A2(1,2) 639 : AA2(2,1)= A2(2,2) 640 : AA2(1,2)= -A2(1,1) 641 : AA2(2,2)= A2(2,1) 641 : CALL KNORM2(NES,AA2,NS,UU2,SS2) 642 : 643 : C.. visibilite a priori 643 : CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU) 643 : CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS) 644 : 645 : IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1 646 : $ ,K2,IVU,IVUS 647 : IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0 648 : 649 : ENDIF 650 : C --------------------------------------------------------------- 651 : C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI 652 : IF(KIMP.GE.4.AND.KVU.NE.0) THEN 653 : WRITE(6,*) ' FACES VISIBLES ',K1,K2 654 : ENDIF 655 : 656 : C>> CALCUL 657 : C... option convexe 658 : IF(KACHE.EQ.0) THEN 659 : IF(KVU.NE.0) THEN 660 : 660 : CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC 661 : $ ,RAD) 662 : IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1 663 : $ ,K2,FF/S1 664 : 665 : ENDIF 666 : C.. option cache 667 : ELSE 668 : IF(KVU.NE.0) THEN 669 : 670 : C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.-------------------------- 671 : 672 : NOBS=0 673 : 674 : C>> BOUCLE FACE 3--------------------------------------------------- 675 : NELT3= 0 676 : DO 3000 KTYP = 1,NTYP 677 : IMODEL = MYMOD.KMODEL(KTYP) 678 : SEGACT IMODEL 679 : IPT3 = IMAMOD 680 : SEGDES IMODEL 681 : 682 : NSGEO3 = IPT3.NUM(/1) 683 : NSPA3=1 684 : IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2 685 : NEL3 = IPT3.NUM(/2) 686 : DO 3100 I3 = 1,NEL3 687 : K3 = I3+ NELT3 688 : 689 : IF(K3.NE.K1.AND.K3.NE.K2) THEN 690 : 691 : DO 3110 IS = 1,NSGEO3,NSPA3 692 : LS=IS 693 : IF(ICOQ.AND.KQUAD(ITYP)) LS 694 : $ =(IS+1)/2 695 : IREF = (IDIM+1)*(IPT3.NUM(IS 696 : $ ,I3)-1) 697 : DO 3120 K = 1,NES 698 : C** A3(K,IS) = XCOOR(IREF+K)/RAD 699 : A3(K,LS) = XCOOR(IREF+K) 700 : $ /RAD 701 : 3120 CONTINUE 702 : 3110 CONTINUE 703 : 704 : ZMIN2=DMIN1(A2(2,1),A2(2,2)) 705 : ZMAX2=DMAX1(A2(2,1),A2(2,2)) 706 : RMAX2=DMAX1(A2(1,1),A2(1,2)) 707 : ZTMIN=DMIN1(ZMIN1,ZMIN2) 708 : ZTMAX=DMAX1(ZMAX1,ZMAX2) 709 : RMAX =DMAX1(RMAX1,RMAX2) 710 : 711 : IF(DMAX1(A3(2,1),A3(2,2)).LE 712 : $ .ZTMIN) THEN 713 : KOBS=0 714 : ELSEIF(DMIN1(A3(2,1),A3(2,2)).GE 715 : $ .ZTMAX) THEN 716 : KOBS=0 717 : ELSEIF(DMIN1(A3(1,1),A3(1,2)).GE 718 : $ .RMAX) THEN 719 : KOBS=0 720 : ELSE 721 : KOBS=1 722 : OBS(1,2*NOBS+1)=A3(1,1) 723 : OBS(2,2*NOBS+1)=A3(2,1) 724 : OBS(1,2*NOBS+2)=A3(1,2) 725 : OBS(2,2*NOBS+2)=A3(2,2) 726 : NOBS=NOBS+1 727 : 728 : ENDIF 729 : ENDIF 730 : 731 : 3100 CONTINUE 732 : C 733 : NELT3=NELT3+NEL3 734 : C 735 : 3000 CONTINUE 736 : 737 : C<< FIN RECHERCHE OBSTRUCTEURS-------------------------------------- 738 : 739 : IF(KIMP.GE.4) THEN 740 : WRITE(6,*) ' FACES K1 K2 ',K1,K2 741 : $ ,' NOBS ',NOBS 742 : ENDIF 743 : 743 : CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX 744 : $ ,FF,KIMP,EXTINC,RAD) 745 : IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF 746 : $ /S1 747 : 748 : ENDIF 749 : 750 : ENDIF 751 : C<< FIN CALCUL 752 : C --------------------------------------------------------------- 753 : 754 : C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG 755 : PLIG.FACT(K2) = FF/S1 756 : 757 : ENDIF 758 : C... fin face K2 .................................................. 759 : 760 : 761 : C... face inverse de K2 (cas des coques) .......................... 762 : 763 : IF (ICOQ.AND.KCOQ(JTYP)) THEN 764 : 765 : K2=K2 + 1 766 : 767 : KVU=1 768 : FF=0.D0 769 : 770 : C UTILISATION DE LA RECIPROCITE 771 : C 772 : IF(K1.GT.K2) THEN 773 : S2=PSUR.FACT(K2) 774 : PCOL=LFACT(K2) 775 : 776 : PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1) 777 : 778 : ELSE 779 : C** 780 : 781 : C>> TEST K1=K2 ET VISIBILITE A PRIORI 782 : 783 : C------------------------------------------------------------------ 784 : 785 : C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0 786 : C 787 : IF(K1.EQ.K2) THEN 788 : IF(KIMP.GE.4) THEN 789 : WRITE(6,*) ' A1 ',A1(1,1),A1(2,1) 790 : WRITE(6,*) ' A1 ',A1(1,2),A1(2,2) 791 : WRITE(6,*) ' U1 ',U1(1),U1(2) 792 : ENDIF 793 : IF (DABS(U1(1)).LT.EPS) KVU=0 794 : IF (U1(1).GE.EPS) KVU=0 795 : 796 : IF(KVU.NE.0) THEN 797 : DO 5140 K = 1,NES 798 : A2(K,1) = A1(K,1) 799 : A2(K,2) = A1(K,2) 800 : 5140 CONTINUE 801 : ENDIF 802 : 803 : C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE 804 : C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z 805 : C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z) 806 : 807 : 808 : ELSE 809 : 810 : DO 5120 K = 1,NES 811 : U2(K)=-U2(K) 812 : XX1=A2(K,1) 813 : A2(K,1) = A2(K,2) + ECOQ*U2(K) 814 : A2(K,2) = XX1 + ECOQ*U2(K) 815 : 5120 CONTINUE 816 : 817 : C... calcul du symetrique 818 : AA2(1,1)= -A2(1,2) 819 : AA2(2,1)= A2(2,2) 820 : AA2(1,2)= -A2(1,1) 821 : AA2(2,2)= A2(2,1) 821 : CALL KNORM2(NES,AA2,NS,UU2,SS2) 822 : 823 : C... visibilité a priori 823 : CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU) 823 : CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS) 824 : 825 : IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ' 826 : $ ,K1,K2,IVU,IVUS 827 : IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0 828 : 829 : ENDIF 830 : C --------------------------------------------------------------- 831 : C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI 832 : IF(KIMP.GE.4.AND.KVU.NE.0) THEN 833 : WRITE(6,*) ' FACES VISIBLES ',K1,K2 834 : ENDIF 835 : 836 : C>> CALCUL -------------------------------------------------------- 837 : C.. option convexe 838 : IF(KACHE.EQ.0) THEN 839 : IF(KVU.NE.0) THEN 840 : 840 : CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP 841 : $ ,EXTINC,RAD) 842 : IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ' 843 : $ ,K1,K2,FF/S1 844 : 845 : ENDIF 846 : C.. option cache 847 : ELSE 848 : IF(KVU.NE.0) THEN 849 : 850 : C>> les obstructeurs potentiels sont les memes que pour la face K2 851 : 852 : IF(KIMP.GE.4) THEN 853 : WRITE(6,*) ' FACES K1 K2 ',K1,K2 854 : $ ,' NOBS ',NOBS 855 : ENDIF 856 : 856 : CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX 857 : $ ,NGAX,FF,KIMP,EXTINC,RAD) 858 : 859 : IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2 860 : $ ,FF/S1 861 : 862 : ENDIF 863 : 864 : ENDIF 865 : C<< FIN CALCUL ------------------------------------------------------ 866 : 867 : C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG 868 : PLIG.FACT(K2) = FF/S1 869 : 870 : ENDIF 871 : 872 : ENDIF 873 : 874 : C... fin face inverse de K2 (cas des coques) ......................... 875 : 876 : 2100 CONTINUE 877 : 878 : IF (ICOQ.AND.KCOQ(JTYP)) THEN 879 : NELT2 = NELT2 + 2 * NEL2 880 : ELSE 881 : NELT2 = NELT2 + NEL2 882 : ENDIF 883 : 884 : 2000 CONTINUE 885 : 886 : C ---FIN BOUCLE FACE 2----------------------------------------- 887 : 888 : ENDIF 889 : C*** fin traitement de la face inverse de K1************************* 890 : 891 : 110 CONTINUE 892 : 893 : 894 : IF (ICOQ.AND.KCOQ(ITYP)) THEN 895 : NELT1 = NELT1 + 2 * NEL1 896 : ELSE 897 : NELT1 = NELT1 + NEL1 898 : ENDIF 899 : 900 : 100 CONTINUE 901 : 902 : C ---FIN BOUCLE FACE 1----------------------------------------- 903 : C 904 : 905 : 906 : 907 : C Désactivation des maillages élémentaires 908 : DO 920 ITYP = 1,NTYP 909 : IMODEL = MYMOD.KMODEL(ITYP) 910 : SEGACT IMODEL 911 : IPT1 = IMAMOD 912 : SEGDES IPT1 913 : SEGDES IMODEL 914 : 920 CONTINUE 915 : SEGDES MYMOD 916 : 917 : C>> SURFACES DIMENSIONNEES 918 : 919 : DO 910 LS=1,NFACE 920 : PSUR.FACT(LS)=PSUR.FACT(LS)*RAD*RAD 921 : 910 CONTINUE 922 : C 923 : IF(KACHE.NE.0) SEGSUP SFOBS 924 : C 925 : 926 : C>>> NORMALISATION,CACUL DES BILANS ET IMPRESSION 927 : C -------------------------------------------- 928 : 929 : IF(EXTINC.GT.1D-3) THEN 930 : INOR=0 931 : ENDIF 932 : 932 : CALL KFN(IFACFO,INOR,KIMP) 933 : 934 : C Traduction puis suppression de IFACFO 935 : 936 : IF (KIMP.GE.3) THEN 936 : CALL PRFACF(IFACFO) 937 : ENDIF 938 : 939 : LTITR=1 939 : CALL FFMCHA(MYMOD,INFOEL,IFACFO,ICHFAC,LTITR) 940 : 941 : SEGACT IFACFO 942 : DO 930 NNEL = 1,LFACT(/1) 943 : LFAC = LFACT(NNEL) 944 : SEGSUP LFAC 945 : 930 CONTINUE 946 : SEGSUP IFACFO 947 : 948 : SEGSUP STRAV 949 : 950 : RETURN 951 : END 952 : C/////////// 953 : C mettre apres KPRIOR 954 : * WRITE (6,*) 'K1 =',K1,' et K2 =',K2 955 : * WRITE (6,*) 'NES =',NES,'NS =',NS 956 : * WRITE (6,*) 'A1(*,1) :',(A1(K,1),K=1,IDIM) 957 : * WRITE (6,*) 'A1(*,2) :',(A1(K,2),K=1,IDIM) 958 : * WRITE (6,*) 'U1 :',(U1(K),K=1,IDIM+1) 959 : * WRITE (6,*) 'A2(*,1) :',(A2(K,1),K=1,IDIM) 960 : * WRITE (6,*) 'A2(*,2) :',(A2(K,2),K=1,IDIM) 961 : * WRITE (6,*) 'U2 :',(U2(K),K=1,IDIM+1) 962 : * WRITE (6,*) 'IVU =',IVU 963 : * WRITE (6,*) 'AA2(*,1) :',(AA2(K,1),K=1,IDIM) 964 : * WRITE (6,*) 'AA2(*,2) :',(AA2(K,2),K=1,IDIM) 965 : * WRITE (6,*) 'UU2 :',(UU2(K),K=1,IDIM+1) 966 : * WRITE (6,*) 'IVUS =',IVUS 967 : C/////////// 968 : 969 : 970 : 971 : 972 : 973 : 974 : 975 : 976 : 977 : 978 : 979 : 980 :
© Cast3M 2003 - All rights reserved.
Disclaimer