1 : C BIOSAV SOURCE CHAT 05/11/29 21:15:02 5153 2 : SUBROUTINE BIOSAV 3 : *---------------------------------------------------------- 4 : * 5 : * CALCUL DES CHAMPS DE BIOT ET SAVART 6 : * 7 : *---------------------------------------------------------- 8 : IMPLICIT INTEGER(I-N) 9 : IMPLICIT REAL*8(A-H,O-Z) 10 : -INC CCOPTIO 11 : -INC CCREEL 12 : -INC SMELEME 13 : -INC SMCOORD 14 : -INC TMTRAV 15 : -INC SMCHPOI 16 : INTEGER ENT 17 : DIMENSION C(3),PN1(3),RESP(3) 18 : DIMENSION VX(3),VZ(3),VY(3),VI(3),PP3(3) 19 : DIMENSION ZF1(1),ZALF(1),ZDENS(1),ZBET(1) 20 : DIMENSION P1(3),P2(3),P3(3),ROTA(3,3) 21 : SEGMENT ICPR(XCOOR(/1)/(IDIM+1)) 22 : SEGMENT SANGLE 23 : REAL*8 TETM(NT4) 24 : REAL*8 TETI(NT4) 25 : REAL*8 DTE(NT4) 26 : ENDSEGMENT 27 : * 28 : CHARACTER*4 MCLE(4),CHSECT(1),CHOIX(2) 29 : LOGICAL INDU,POTE 30 : * 31 : DATA EPSIL3/1.D-12/,EPSIL2/1.D-8/ 32 : * 33 : DATA MCLE/'CERC','ARC ','BARR','FIL '/ 34 : DATA CHSECT/'TRAP'/ 35 : DATA CHOIX/'INDU','POTE'/ 36 : * 36 : CALL LIROBJ('CHPOINT ',MCHPO1,0,IRETOU) 37 : IF (IERR.NE.0) RETURN 38 : IF (MCHPO1.GT.0) THEN 38 : CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) 39 : IF (IERR.NE.0) RETURN 40 : 41 : * OPTION CALCUL PAR INTEGRALE ELLIPTIQUE (PHR 2002) 42 : 42 : CALL BIOKAM (MCHPO1,IPT1,ENT) 43 : IF (IERR.NE.0) RETURN 44 : MCHPOI = ENT 45 : GOTO 9990 46 : ENDIF 47 : * CHOIX DE L'OPTION DE CALCUL 48 : * (INDUCTION PAR DEFAUT) 49 : * 50 : INDU=.TRUE. 51 : POTE=.FALSE. 51 : CALL LIRMOT(CHOIX,2,IMC,0) 52 : IF (IERR.NE.0) RETURN 53 : IF(IMC.EQ.2) THEN 54 : POTE=.TRUE. 55 : INDU=.FALSE. 56 : ENDIF 57 : * 58 : * LECTURE D'UN MAILLAGE 59 : * 59 : CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU) 60 : IF (IERR.NE.0) RETURN 61 : * 62 : * LECTURE D'UN MOT CLE 63 : * 63 : CALL LIRMOT(MCLE,4,IMC,1) 64 : IF (IERR.NE.0) RETURN 65 : * 66 : * IMC = 1 ON A LU CERCLE 67 : * IMC = 2 ON A LU ARC 68 : * IMC = 3 ON A LU BARRE 69 : * IMC = 4 ON A LU FIL 70 : * 71 : * 72 : * LECTURE DE 2 ou 3 POINTS 73 : * 73 : CALL LIROBJ('POINT ',IP1,1,IRETOU) 73 : CALL LIROBJ('POINT ',IP2,1,IRETOU) 74 : IP3=0 75 : IF(IMC.LT.4) THEN 75 : CALL LIROBJ('POINT ',IP3,1,IRETOU) 76 : ELSE 77 : C CAS DU FIL RECTILIGNE 78 : GO TO 10 79 : ENDIF 80 : IF (IERR.NE.0) RETURN 81 : * 82 : * LECTURE DE 4 OU 5 FLOTTANTS SELON LES CAS 83 : * 83 : CALL LIRREE(F1,1,IRETOU) 83 : CALL LIRREE(F2,1,IRETOU) 84 : F3=0.D0 84 : IF(IMC.LT.3) CALL LIRREE(F3,1,IRETOU) 85 : * 86 : * LECTURE DU TYPE DE SECTION 87 : * 88 : ISECT=0 88 : CALL LIRMOT(CHSECT,1,ISECT,0) 89 : IF(ISECT.EQ.1) THEN 89 : CALL LIRREE(PENT1,1,IRETOU) 89 : CALL LIRREE(PENT2,1,IRETOU) 90 : IF(IERR.NE.0) RETURN 91 : ELSE 92 : ISECT=3 93 : ENDIF 94 : 10 CONTINUE 94 : CALL LIRREE(DENS,1,IRETOU) 94 : CALL LIRREE(XIU,1,IRETOU) 95 : C ON VIENT DE LIRE MU0 DANS L UNITE UTILISATEUR 96 : RUMKSA= 4.D0 * XPI * 1.D-7 97 : RAP=(XIU / RUMKSA) 98 : IU= 1 99 : IF (IERR.NE.0) RETURN 100 : ZF1(1)=F1 101 : ZDENS(1)=DENS 102 : * 103 : * EXEMPLE DE RECUPERATION DES COORDONNEES DU POINT P1 104 : * 105 : IREF=(IDIM+1)*IP1 106 : P1(1)=XCOOR(IREF-IDIM) 107 : P1(2)=XCOOR(IREF-IDIM+1) 108 : P1(3)=0.D0 109 : IF (IDIM.EQ.3) P1(3)=XCOOR(IREF-IDIM+2) 110 : IREF=(IDIM+1)*IP2 111 : P2(1)=XCOOR(IREF-IDIM) 112 : P2(2)=XCOOR(IREF-IDIM+1) 113 : P2(3)=0.D0 114 : IF (IDIM.EQ.3) P2(3)=XCOOR(IREF-IDIM+2) 115 : IF(IP3.NE.0) THEN 116 : IREF=(IDIM+1)*IP3 117 : P3(1)=XCOOR(IREF-IDIM) 118 : P3(2)=XCOOR(IREF-IDIM+1) 119 : P3(3)=0.D0 120 : IF (IDIM.EQ.3) P3(3)=XCOOR(IREF-IDIM+2) 121 : ENDIF 122 : * 123 : * ----------------------- 124 : * | LES CALCULS .... | 125 : * ----------------------- 126 : * 127 : * ON COMMENCE PAR EXPLORER LE MAILLAGE 128 : * 129 : SEGACT MELEME 129 : IF (ITYPEL.NE.1.OR.LISOUS(/1).NE.0) CALL CHANGE(MELEME,1) 130 : * 131 : * ON VERIFIE QUE CHAQUE POINT N'APPARAIT QU'UNE SEULE FOIS 132 : * 133 : NNNOE=NUM(/2) 134 : SEGINI ICPR 135 : ICON=0 136 : DO 12 I=1,NNNOE 137 : IKI=NUM(1,I) 138 : IF(ICPR(IKI).NE.0) GO TO 12 139 : ICON=ICON+1 140 : ICPR(IKI)=ICON 141 : 12 CONTINUE 142 : SEGSUP ICPR 143 : * 144 : IF(ICON.NE.NNNOE) THEN 144 : CALL ERREUR(75) 145 : SEGDES MELEME 146 : RETURN 147 : ENDIF 148 : * 149 : * ON INITIALISE LES NOMS DES COMPOSANTES ET LES HARMONIQUES 150 : * POUR L'INSTANT LIMITE AU 3D ET 2D PLAN 151 : * 152 : * 153 : NNIN=IDIM 154 : SEGINI MTRAV 155 : NHRM=NIFOUR 156 : IF(NIFOUR.EQ.1) NHRM=IFOUR 157 : INCO(1)='BX ' 158 : IF(POTE) INCO(1)='AX ' 159 : INCO(2)='BY ' 160 : IF(POTE) INCO(2)='AY ' 161 : NHAR(1)=NHRM 162 : NHAR(2)=NHRM 163 : IF(IDIM.EQ.3) THEN 164 : INCO(3)='BZ ' 165 : IF(POTE) INCO(3)='AZ ' 166 : NHAR(3)=NHRM 167 : ENDIF 168 : C 169 : GOTO (100,100,200,300),IMC 170 : 100 CONTINUE 171 : C BOBINE CIRCULAIRE 172 : 173 : C ETABLISSEMENT DU REPERE LOCAL 174 : C AXE DES Z (P2 - P1 ) * (P3 - P1) 175 : DO 801 I=1,3 176 : VX(I)=P2(I)- P1(I) 177 : VI(I)=P3(I)- P1(I) 178 : 801 CONTINUE 178 : CALL PROVEC(VX,VI,VZ) 178 : CALL PROVEC(VZ,VX,VY) 178 : CALL ROTREP (ROTA,VX,VY,VZ) 179 : IF(IERR.NE.0) THEN 180 : SEGDES MELEME 181 : SEGSUP MTRAV 182 : RETURN 183 : ENDIF 184 : NB=1 185 : IF(F1.EQ.0.D0) RETURN 186 : ALF=F2/F1 187 : BET=F3/(2.D0*F1) 188 : ZALF(1)=ALF 189 : ZBET(1)=BET 190 : IC= 3 191 : DO 802 J=1,3 192 : PP3(J)= P3(J) - P1(J) 193 : 802 CONTINUE 194 : RAYO2=SDOT(3,VX,1,VX,1) 195 : PROD=SDOT(3,VX,1,VI,1) 196 : RAYO2=RAYO2+EPSIL3 197 : TM=ACOS(PROD/RAYO2) 198 : GOTO 350 199 : 200 CONTINUE 200 : C CAS DES BARRES ON VA TOURNER / OX POUR AMENER OY // P1P2 201 : GY=0.D0 202 : DO 803 I=1,3 203 : VX(I)=P2(I)- P1(I) 204 : VI(I)=P3(I)- P1(I) 205 : GY=GY +( VX(I)*VX(I)) 206 : C(I)= (P2(I)+ P1(I))/2.D0 207 : 803 CONTINUE 208 : GY= SQRT( GY) 208 : CALL PROVEC(VX,VI,VZ) 208 : CALL PROVEC(VZ,VX,VY) 208 : CALL ROTREP (ROTA,VX,VY,VZ) 209 : IF(IERR.NE.0) THEN 210 : SEGDES MELEME 211 : SEGSUP MTRAV 212 : RETURN 213 : ENDIF 214 : DO 804 I= 1,3 215 : VZ(I)=0.D0 216 : DO 8041 J= 1,3 217 : VZ(I)= VZ(I) + ROTA(I,J) * C(J) 218 : 8041 CONTINUE 219 : 804 CONTINUE 219 : * CALL ROTATI (C,1,1,1) 220 : GO TO 350 221 : C 222 : C CAS DU FIL 223 : 300 CONTINUE 224 : C 225 : 350 CONTINUE 226 : UPSI =1.D-9 227 : * 228 : * ON BOUCLE SUR LES NOEUDS 229 : * 230 : NT4=8000 231 : SEGINI SANGLE 232 : DO 30 IPT=1,NNNOE 233 : L=NUM(1,IPT) 234 : IREF= (IDIM+1)*L 235 : * L EST LE NUMERO DU NOEUD DANS LE TABLEAU DES COORDONNEES 236 : DO 805 J=1,3 237 : PN1(J)=XCOOR(IREF-IDIM-1+J) 238 : 805 CONTINUE 239 : * 240 : GOTO (400,400,500,600),IMC 241 : 400 CONTINUE 242 : * CAS DES CERCLES ET DES ARC 243 : C MESURE DU POINT DANS LE REPERE LOCAL 244 : DO 806 J=1,3 245 : PN1(J)=PN1(J) - P1(J) 246 : 806 CONTINUE 247 : DO 807 I= 1,3 248 : VZ(I)=0.D0 249 : VY(I)=0.D0 250 : DO 8071 J= 1,3 251 : VY(I)= VY(I) + ROTA(I,J) * PN1(J) 252 : VZ(I)= VZ(I) + ROTA(I,J) * PP3(J) 253 : 8071 CONTINUE 254 : 807 CONTINUE 254 : * CALL ROTATI (PP3,1,1,1) 254 : * CALL ROTATI (PN1,1,1,1) 255 : C ON CHANGE DE RERERENTIEL 256 : C MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS 257 : 258 : ZP= VY(3) 259 : RP= SQRT ( VY(1)**2 + VY(2)**2 ) 260 : IF(RP.EQ.0.D0) THEN 261 : TMIN=0.D0 262 : COST=1.D0 263 : SINT=0.D0 264 : RP=EPSIL2 265 : ELSE 266 : TMIN=ACOS(VY(1)/(RP+EPSIL3)) 267 : COST=COS(TMIN) 268 : SINT=VY(2)/(RP+EPSIL3) 269 : ENDIF 270 : DO 808 J= 1,3 271 : C(J)=0.D0 272 : 808 CONTINUE 273 : C 274 : IF(BET.EQ.0.D0.AND.ALF.EQ.1.D0) THEN 275 : C CAS DU FIL CIRCULAIRE 276 : RAYON=F1 277 : COUR=DENS 278 : C ON CONSTRUIT LE REPERE LOCAL 279 : sigt=1.d0 280 : if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT) 281 : TMIN=SIGT*(TMIN-XPI)+XPI 282 : IF(IMC.EQ.1) TM=2.D0*XPI 283 : TMAX=TMIN+TM 283 : CALL TEGRAL(SANGLE) 284 : IF(INDU) THEN 284 : CALL BFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ) 285 : ENDIF 286 : IF(POTE) THEN 286 : CALL AFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ) 287 : ENDIF 288 : RESP(1)=BX*COST-BY*SINT 289 : RESP(2)=BX*SINT+BY*COST 290 : RESP(3)=BZ 291 : ELSE IF(BET.EQ.0.D0.AND.ISECT.EQ.3) THEN 292 : C CAS DE LA PLAQUE CIRCULAIRE 293 : RINT=F1 294 : REXT=F2 295 : DR=REXT-RINT 296 : COUR=DENS*DR 297 : C ON CONSTRUIT LE REPERE LOCAL 298 : sigt=1.d0 299 : if( abs(sint).gt.1.e-6) SIGT=SIGN(1.D0,-SINT) 300 : TMIN=SIGT*(TMIN-XPI)+XPI 301 : IF(IMC.EQ.1) TM=2.D0*XPI 302 : TMAX=TMIN+TM 302 : CALL TEGRAL(SANGLE) 303 : IF(INDU) THEN 303 : CALL BPLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP, 304 : & BX,BY,BZ) 305 : ENDIF 306 : IF(POTE) THEN 306 : CALL APLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP, 307 : & BX,BY,BZ) 308 : ENDIF 309 : RESP(1)=BX*COST-BY*SINT 310 : RESP(2)=BX*SINT+BY*COST 311 : RESP(3)=BZ 312 : ELSE IF(ALF.EQ.1D0) THEN 313 : C CAS DU CYLINDRE 314 : RAYON=F1 315 : H=F3 316 : COUR=DENS*H 317 : C ON CONSTRUIT LE REPERE LOCAL 318 : sigt=1.d0 319 : if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT) 320 : TMIN=SIGT*(TMIN-XPI)+XPI 321 : IF(IMC.EQ.1) TM=2.D0*XPI 322 : TMAX=TMIN+TM 322 : CALL TEGRAL(SANGLE) 323 : IF(INDU) THEN 323 : CALL BCYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP, 324 : & BX,BY,BZ) 325 : ENDIF 326 : IF(POTE) THEN 326 : CALL ACYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP, 327 : & BX,BY,BZ) 328 : ENDIF 329 : RESP(1)=BX*COST-BY*SINT 330 : RESP(2)=BX*SINT+BY*COST 331 : RESP(3)=BZ 332 : ELSE 333 : IF (IMC.EQ.1.AND.ISECT.EQ.3) THEN 334 : IF(INDU) THEN 334 : CALL CHAMBO(IU,IC,NB,C(3),ZF1,ZALF,ZBET,ZDENS,ZP,RP,BZ 335 : $ ,BR,BP) 336 : IF (RP.LE.UPSI) THEN 337 : RESP(1)=0.D0 338 : RESP(2)=0.D0 339 : ELSE 340 : RESP(1)=BR * VY(1) / RP 341 : RESP(2)=BR * VY(2) / RP 342 : END IF 343 : RESP(3)=BZ 344 : GO TO 700 345 : ENDIF 346 : IF(POTE) THEN 347 : RINT=F1 348 : REXT=F2 349 : HCEN=F3 350 : WIDTH=F2-F1 351 : AIRE=F3*WIDTH 352 : COUR=DENS*AIRE 353 : DELTA=0.D0 354 : C ON CONSTRUIT LE REPERE LOCAL 355 : SIGT=1.D0 356 : if(ABS(SINT).gt.1.e-8)SIGT=SIGN(1.D0,-SINT) 357 : TMIN=SIGT*(TMIN-XPI)+XPI 358 : TM=2.D0*XPI 359 : TMAX=TMIN+TM 360 : PENT1=0.D0 361 : PENT2=0.D0 361 : CALL TEGRAL(SANGLE) 361 : CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT, 362 : * PENT1,PENT2,TMIN,TMAX, 363 : * RP,ZP,BX,BY,BZ) 364 : RESP(1)=BX*COST-BY*SINT 365 : RESP(2)=BX*SINT+BY*COST 366 : RESP(3)=BZ 367 : GO TO 700 368 : ENDIF 369 : ENDIF 370 : IF(ISECT.EQ.3) THEN 371 : C CAS DES SECTIONS RECTANGULAIRES 372 : C ON CALCULE POUR UN ARC ---> ANGLES 373 : A1=0.D0 374 : IF (ABS(VZ(1)).LT.UPSI) THEN 375 : A2=XPI/2.D0 376 : ELSE 377 : AA= VZ(2)/VZ(1) 378 : A2=ATAN ( AA ) 379 : IF ( A2.LT.0.D0) A2 = XPI + A2 380 : END IF 381 : IF(INDU) THEN 381 : CALL CHAMAR(IU,IC,C,A1,A2,F1,ALF,BET,DENS, 382 : * VY(1),VY(2),VY(3),BX,BY,BZ) 383 : RESP(1)=BX 384 : RESP(2)=BY 385 : RESP(3)=BZ 386 : ENDIF 387 : IF(POTE) THEN 388 : RINT=F1 389 : REXT=F2 390 : HCEN=F3 391 : WIDTH=F2-F1 392 : AIRE=F3*WIDTH 393 : COUR=DENS*AIRE 394 : PENT1=0.D0 395 : PENT2=0.D0 396 : C ON CONSTRUIT LE REPERE LOCAL 397 : sigt=1.d0 398 : if( abs(sint).gt.1.e-6)SIGT=SIGN(1.d0,-SINT) 399 : TMIN=SIGT*(TMIN-XPI)+XPI 400 : TMAX=TMIN+TM 400 : CALL TEGRAL(SANGLE) 400 : CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT, 401 : * PENT1,PENT2,TMIN,TMAX, 402 : * RP,ZP,BX,BY,BZ) 403 : RESP(1)=BX*COST-BY*SINT 404 : RESP(2)=BX*SINT+BY*COST 405 : RESP(3)=BZ 406 : ENDIF 407 : ELSEIF(ISECT.EQ.1) THEN 408 : C CAS DES SECTIONS TRAPEZOIDALES 409 : RINT=F1 410 : REXT=F2 411 : HCEN=F3 412 : WIDTH=F2-F1 413 : AIRE=F3*WIDTH 414 : COUR=DENS*AIRE 415 : DELTA=0.5D0*WIDTH*(PENT2-PENT1) 416 : HINT=HCEN-DELTA 417 : HEXT=HCEN+DELTA 418 : C ON CONSTRUIT LE REPERE LOCAL 419 : SIGT=1.d0 420 : if(abs(sint).gt.1.e-06)SIGT=SIGN(1.D0,-SINT) 421 : TMIN=SIGT*(TMIN-XPI)+XPI 422 : IF(IMC.EQ.1) TM=2.D0*XPI 423 : TMAX=TMIN+TM 423 : CALL TEGRAL(SANGLE) 424 : IF(HCEN.EQ.0.D0) THEN 425 : C CAS DU TRONC DE CONE 426 : TETA=ATAN(PENT1) 427 : COUR=DENS*WIDTH/COS(TETA) 428 : IF(INDU) THEN 428 : CALL BTRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX, 429 : & RP,ZP,BX,BY,BZ) 430 : ENDIF 431 : IF(POTE) THEN 431 : CALL ATRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX, 432 : & RP,ZP,BX,BY,BZ) 433 : ENDIF 434 : ELSE 435 : IF(INDU) THEN 435 : CALL ARCTRA(SANGLE,COUR,HINT,HEXT,RINT,REXT, 436 : * PENT1,PENT2,TMIN,TMAX, 437 : * RP,ZP,BX,BY,BZ) 438 : ENDIF 439 : IF(POTE) THEN 439 : CALL AARCTR(SANGLE,COUR,HINT,HEXT,RINT,REXT, 440 : * PENT1,PENT2,TMIN,TMAX, 441 : * RP,ZP,BX,BY,BZ) 442 : ENDIF 443 : ENDIF 444 : RESP(1)=BX*COST-BY*SINT 445 : RESP(2)=BX*SINT+BY*COST 446 : RESP(3)=BZ 447 : ENDIF 448 : END IF 449 : IF(IERR.NE.0) THEN 450 : SEGDES MELEME 451 : SEGSUP MTRAV 452 : RETURN 453 : ENDIF 454 : C 455 : GOTO 700 456 : * 457 : 500 CONTINUE 458 : C CAS DES BARRES ON CHANGE DE RERERENTIEL 459 : DO 809 I= 1,3 460 : VY(I)=0.D0 461 : DO 8091 J= 1,3 462 : VY(I)= VY(I) + ROTA(I,J) * PN1(J) 463 : 8091 CONTINUE 464 : 809 CONTINUE 464 : * CALL ROTATI (PN1,1,1,1) 465 : C MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS 466 : IF(ISECT.EQ.3) THEN 467 : C CAS DES SECTIONS RETANGULAIRES 468 : IF(F2.EQ.0.D0) THEN 469 : C CAS DE LA PLAQUE RECTANGULAIRE 470 : DO 813 I=1,3 471 : VY(I)=0.D0 472 : DO 8131 J=1,3 473 : VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) 474 : 8131 CONTINUE 475 : 813 CONTINUE 476 : DY=F1 477 : REXT=0.5D0*DY 478 : RINT=-REXT 479 : COUR=DENS*DY 480 : YMIN=-VY(1) 481 : YMAX=YMIN+GY 482 : XP=-VY(2) 483 : ZP=VY(3) 484 : IF(INDU) THEN 484 : CALL BPLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP 485 : $ ) 486 : ENDIF 487 : IF(POTE) THEN 487 : CALL APLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP 488 : $ ) 489 : ENDIF 490 : RESP(1)=BYP 491 : RESP(2)=-BXP 492 : RESP(3)= BZP 493 : ELSE 494 : IF(INDU) THEN 494 : CALL CHABAR(IU,VZ,GY,F1,F2,DENS,VY(1),VY(2),VY(3), 495 : * BXP,BYP,BZP) 496 : RESP(1)=BXP 497 : RESP(2)=BYP 498 : RESP(3)=BZP 499 : ENDIF 500 : IF(POTE) THEN 501 : DO 814 I=1,3 502 : VY(I)=0.D0 503 : DO 8141 J=1,3 504 : VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) 505 : 8141 CONTINUE 506 : 814 CONTINUE 507 : DY=F1 508 : DZ=F2 509 : REXT=0.5D0*DY 510 : RINT=-REXT 511 : AIRE=DY*DZ 512 : COUR=DENS*AIRE 513 : PENT1=0.D0 514 : PENT2=0.D0 515 : YMIN=-VY(1) 516 : YMAX=YMIN+GY 517 : XP=-VY(2) 518 : ZP=VY(3) 518 : CALL ASEGTR(COUR,DZ,DZ,RINT,REXT,PENT1,PENT2, 519 : * YMIN,YMAX,XP,ZP,BXP,BYP,BZP) 520 : RESP(1)=BYP 521 : RESP(2)=-BXP 522 : RESP(3)= BZP 523 : ENDIF 524 : ENDIF 525 : ELSE IF(ISECT.EQ.1) THEN 526 : C CAS DES SECTIONS TRAPEZOIDALES 527 : DO 812 I=1,3 528 : VY(I)=0.D0 529 : DO 8121 J=1,3 530 : VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J)) 531 : 8121 CONTINUE 532 : 812 CONTINUE 533 : DY=F1 534 : DZ=F2 535 : REXT=0.5D0*DY 536 : RINT=-REXT 537 : AIRE=DY*DZ 538 : COUR=DENS*AIRE 539 : DELTA=0.5D0*DY*(PENT2-PENT1) 540 : HINT=DZ-DELTA 541 : HEXT=DZ+DELTA 542 : YMIN=-VY(1) 543 : YMAX=YMIN+GY 544 : XP=-VY(2) 545 : ZP=VY(3) 546 : IF(INDU) THEN 546 : CALL SEGTRA(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2, 547 : * YMIN,YMAX,XP,ZP,BXP,BYP,BZP) 548 : ENDIF 549 : IF(POTE) THEN 549 : CALL ASEGTR(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2, 550 : * YMIN,YMAX,XP,ZP,BXP,BYP,BZP) 551 : ENDIF 552 : RESP(1)=BYP 553 : RESP(2)=-BXP 554 : RESP(3)= BZP 555 : ENDIF 556 : GO TO 700 557 : C 558 : C CAS DU FIL 559 : 600 CONTINUE 560 : GY=0.D0 561 : DO 601 I=1,3 562 : VY(I)=P2(I)-P1(I) 563 : VI(I)=PN1(I)-P1(I) 564 : GY=GY+(VY(I)*VY(I)) 565 : 601 CONTINUE 566 : GY=SQRT(GY) 566 : CALL PROVEC(VI,VY,VZ) 566 : CALL PROVEC(VY,VZ,VX) 566 : CALL ROTREP (ROTA,VX,VY,VZ) 567 : IF(IERR.NE.0) THEN 568 : SEGDES MELEME 569 : SEGSUP MTRAV 570 : RETURN 571 : ENDIF 572 : XP=0.D0 573 : YMIN=0.D0 574 : DO 602 J= 1,3 575 : XP=XP+ROTA(1,J)*VI(J) 576 : YMIN=YMIN-ROTA(2,J)*VI(J) 577 : 602 CONTINUE 578 : YMAX=YMIN+GY 579 : COUR=DENS 580 : IF(INDU) THEN 580 : CALL BFILDR(COUR,YMIN,YMAX,XP,BZP) 581 : RESP(1)=0.D0 582 : RESP(2)=0.D0 583 : RESP(3)=BZP 584 : ENDIF 585 : IF(POTE) THEN 585 : CALL AFILDR(COUR,YMIN,YMAX,XP,AYP) 586 : RESP(1)=0.D0 587 : RESP(2)=AYP 588 : RESP(3)=0.D0 589 : ENDIF 590 : CC 591 : 700 CONTINUE 592 : * 593 : C RETOUR DU CHAMP DANS LE REPERE GENERAL 594 : 595 : DO 810 I= 1,3 596 : VY(I)=0.D0 597 : DO 8101 J= 1,3 598 : VY(I)= VY(I) + ROTA(J,I) * RESP(J) * RAP 599 : 8101 CONTINUE 600 : 810 CONTINUE 600 : * CALL RET_REF (RESP) 601 : * REMPLISSAGE DU SEGMENT DE TRAVAIL MTRAV 602 : * 603 : IGEO(IPT)=L 604 : DO 811 J= 1,3 605 : BB(J,IPT) = VY(J) 606 : IBIN(J,IPT) = 1 607 : 811 CONTINUE 608 : C FIN BOUCLE SUR LES POINTS 609 : * 610 : 30 CONTINUE 611 : SEGSUP SANGLE 612 : SEGDES,MELEME 613 : * 613 : CALL CRECHP(MTRAV,MCHPOI) 614 : SEGSUP MTRAV 615 : * 616 : * ATTRIBUTION D'UNE NATURE DISCRETE AU CHAMP QUI SORT 617 : * 618 : 619 : 9990 SEGACT MCHPOI*MOD 620 : JATTRI(1) = 2 621 : SEGDES MCHPOI 622 : * 623 : * 623 : CALL ECROBJ('CHPOINT ',MCHPOI) 624 : RETURN 625 : END 626 : 627 : 628 : 629 : 630 : 631 : 632 : 633 : 634 : 635 : 636 : 637 : 638 : 639 : 640 : 641 : 642 :
© Cast3M 2003 - All rights reserved.
Disclaimer