1 : C CRCOU2 SOURCE PV 09/03/12 21:18:36 6325 2 : C CALCUL D'UN OBJET COUPE VIRTUEL (CAS DE DEFORMEES) 3 : C 4 : C CRCOU2 REMPLACE MELEME PAR L'OBJET COUPE 5 : C MCOUP ET LE DERNIER COMPOSANT DE MELEME REMPLACE DECRIVENT 6 : C L'INTERSECTION 7 : C IL FAUDRA METTRE UNE INDICATION AUX TRIANGLES SE TROUVANT DANS LE PLA 8 : C DE COUPE POUR REUSSIR LES PARTIES VUES PARTIES CACHEES 9 : SUBROUTINE CRCOU2(JOEIL,ICOUP1,ICOUP2,ICOUP3,MELEME,MCOUP,VCHCA, 10 : # XCORD,ICPR,MELEM2,MCOU2,ITE,IVU,mchamor,isect) 11 : IMPLICIT INTEGER(I-N) 12 : -INC CCGEOME 13 : -INC CCOPTIO 14 : -INC SMCOORD 15 : -INC SMELEME 16 : -INC SMCHAML 17 : SEGMENT XCORD(IDIM,ITE) 18 : REAL*8 XN,YN,ZN,SIGA,SIGB,SIGC,SIGRAV 19 : SEGMENT IVU(0) 20 : SEGMENT VCHCA(0) 21 : SEGMENT ITR(0) 22 : SEGMENT ITRT(0) 23 : SEGMENT XTR(0) 24 : SEGMENT XTRT(0) 25 : SEGMENT IQUAT(0) 26 : SEGMENT XQUAT(0) 27 : SEGMENT ISEGM(0) 28 : SEGMENT XSEGM(NBPOIN) 29 : SEGMENT MCOUP(0) 30 : SEGMENT ICPR(0) 31 : SEGMENT IJGRAV(2,100) 32 : DIMENSION JTT(9),ISEG(9) 33 : LOGICAL VOLUM 34 : n2ptel=0 35 : n2el=0 36 : mcham=mchamor 37 : MEDSAU=0 38 : SEGINI ITR,ISEGM,MCOUP,IJGRAV,ITRT,IQUAT 39 : SEGACT XCORD,MCOORD 40 : NBPOIN=XCOOR(/1)/4 41 : if (mcham.ne.0) segini xtr,xtrt,xquat,xsegm 42 : ITE=XCORD(/2) 43 : IREF=(JOEIL-1)*4 44 : XO=XCOOR(IREF+1) 45 : YO=XCOOR(IREF+2) 46 : ZO=XCOOR(IREF+3) 47 : IREF=(ICOUP1-1)*4 48 : XCT=XCOOR(IREF+1) 49 : YCT=XCOOR(IREF+2) 50 : ZCT=XCOOR(IREF+3) 51 : IREF=(ICOUP2-1)*4 52 : XC2=XCOOR(IREF+1) 53 : YC2=XCOOR(IREF+2) 54 : ZC2=XCOOR(IREF+3) 55 : IREF=(ICOUP3-1)*4 56 : XC3=XCOOR(IREF+1) 57 : YC3=XCOOR(IREF+2) 58 : ZC3=XCOOR(IREF+3) 59 : C NORMALE AU PLAN : 60 : XV=XC2-XCT 61 : YV=YC2-YCT 62 : ZV=ZC2-ZCT 63 : XW=XC3-XCT 64 : YW=YC3-YCT 65 : ZW=ZC3-ZCT 66 : XN=YV*ZW-ZV*YW 67 : YN=ZV*XW-XV*ZW 68 : ZN=XV*YW-YV*XW 69 : RN=XN**2+YN**2+ZN**2 69 : IF (RN.EQ.0.) CALL ERREUR(21) 70 : IF (IERR.NE.0) RETURN 71 : RN=SQRT(RN) 72 : XN=XN/RN 73 : YN=YN/RN 74 : ZN=ZN/RN 75 : SIG=(XO-XCT)*XN+(YO-YCT)*YN+(ZO-ZCT)*ZN 76 : IF (SIG.LE.0.) THEN 77 : XN=-XN 78 : YN=-YN 79 : ZN=-ZN 80 : ENDIF 81 : * DANS IVU ON NOTE LES POINTS DU MAUVAIS COTE DU PLAN DE COUPE 82 : DO 60 I=1,XCORD(/2) 83 : XP=XCORD(1,I) 84 : YP=XCORD(2,I) 85 : ZP=XCORD(3,I) 86 : SIG=(XP-XCT)*XN+(YP-YCT)*YN+(ZP-ZCT)*ZN 87 : IF (SIG.GT.0.) IVU(I)=0 88 : 60 CONTINUE 89 : 5000 CONTINUE 90 : SEGACT MELEME 91 : NBSOU=LISOUS(/1) 92 : NBSOUS=MAX(2,NBSOU+1) 93 : NBREF=0 94 : if (mcham.ne.0) nbref=nbsous 95 : NBNN=0 96 : NBELEM=0 97 : SEGINI IPT8 98 : IPT1=MELEME 99 : ISU=0 100 : * prob optimiseur 101 : melva1=meleme 102 : melva2=meleme 103 : DO 10 ISOUS=1,MAX(1,NBSOU) 104 : IF (NBSOU.NE.0) THEN 105 : IPT1=LISOUS(ISOUS) 106 : if (mcham.ne.0) then 107 : melva1=lisref(isous) 108 : segact melva1 109 : lva1=melva1.velche(/1) 110 : endif 111 : SEGACT IPT1 112 : ENDIF 113 : NBREF=0 114 : NBSOUS=0 115 : NBNN=IPT1.NUM(/1) 116 : NBELEM=IPT1.NUM(/2) 117 : SEGINI IPT2 118 : IPT2.ITYPEL=IPT1.ITYPEL 119 : if (mcham.ne.0) then 120 : n1ptel=nbnn 121 : n1el=nbelem 122 : segini melva2 123 : endif 124 : JELN=0 125 : DO 20 JEL=1,NBELEM 126 : DO 25 INOEU=1,NBNN 127 : IPT=IPT1.NUM(INOEU,JEL) 128 : IF (IVU(ICPR(IPT)).EQ.0) GOTO 30 129 : 25 CONTINUE 130 : if (isect.ne.0) goto 20 131 : JELN=JELN+1 132 : * write (6,*) ' mcham melva1 melva2 meleme ',mcham,melva1,melva2, 133 : * > meleme 134 : DO INOEU=1,NBNN 135 : IPT2.NUM(INOEU,JELN)=IPT1.NUM(INOEU,JEL) 136 : ENDDO 137 : if (mcham.ne.0) then 138 : DO INOEU=1,NBNN 139 : melva2.velche(INOEU,JELN)=melva1.velche(min(INOEU,lva1),JEL) 140 : ENDDO 141 : ENDIF 142 : IPT2.ICOLOR(JELN)=IPT1.ICOLOR(JEL) 143 : GOTO 20 144 : 30 CONTINUE 145 : VOLUM=KSURF(IPT1.ITYPEL).NE.IPT1.ITYPEL 146 : * IF (VOLUM) THEN 147 : IPT=IPT1.NUM(1,JEL) 148 : XGRAV=XCORD(1,ICPR(IPT)) 149 : YGRAV=XCORD(2,ICPR(IPT)) 150 : ZGRAV=XCORD(3,ICPR(IPT)) 151 : IF (VCHCA.NE.0) VGRAV=VCHCA(IPT) 152 : if (mcham.ne.0) chgrav=melva1.velche(1,JEL) 153 : DO 21 INOEU=2,NBNN 154 : IPT=IPT1.NUM(INOEU,JEL) 155 : XGRAV=XGRAV+XCORD(1,ICPR(IPT)) 156 : YGRAV=YGRAV+XCORD(2,ICPR(IPT)) 157 : ZGRAV=ZGRAV+XCORD(3,ICPR(IPT)) 158 : IF (VCHCA.NE.0) VGRAV=VGRAV+VCHCA(IPT) 159 : if (mcham.ne.0) chgrav=chgrav+melva1.velche(min(INOEU,lva1),JEL) 160 : 21 CONTINUE 161 : XGRAV=XGRAV/NBNN 162 : YGRAV=YGRAV/NBNN 163 : ZGRAV=ZGRAV/NBNN 164 : C ON PREND COMME POINT DU PLAN LA PROJECTION DU CENTRE GRAVITE 165 : SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN 166 : XCT=XGRAV-SIGRAV*XN 167 : YCT=YGRAV-SIGRAV*YN 168 : ZCT=ZGRAV-SIGRAV*ZN 169 : SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN 170 : IF (ABS(SIGRAV).LE.1E-30) SIGRAV=1E-30 171 : NJGRAV=0 172 : IF (VCHCA.NE.0) VGRAV=VGRAV/NBNN 173 : if (mcham.ne.0) CHGRAV=CHGRAV/NBNN 174 : * ENDIF 175 : C DECOMPOSITION EN FACE PUIS EN TRIANGLE 176 : NBFAC=LTEL(1,IPT1.ITYPEL) 177 : IAD=LTEL(2,IPT1.ITYPEL)-1 178 : IF (NBFAC.EQ.0) GOTO 20 179 : DO 161 IFAC=1,NBFAC 180 : ITYP=LDEL(1,IAD+IFAC) 181 : JAD=LDEL(2,IAD+IFAC)-1 182 : NPFAC=KDFAC(1,ITYP) 183 : IDEP=KDFAC(2,ITYP) 184 : IFEP=IDEP+3*(KDFAC(3,ITYP)-1) 185 : DO 160 ITRIAN=IDEP,IFEP,3 186 : JFA=JAD+KFAC(ITRIAN) 187 : JFB=JAD+KFAC(ITRIAN+1) 188 : JFC=JAD+KFAC(ITRIAN+2) 189 : ISEG(1)=0 190 : ISEG(2)=0 191 : ISEG(3)=0 192 : ISEG(4)=0 193 : ISEG1=0 194 : ISEG2=0 195 : ISEG3=0 196 : IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.1) ISEG1=1 197 : IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.NPFAC-1) ISEG1=1 198 : IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.1) ISEG2=1 199 : IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.NPFAC-1) ISEG2=1 200 : IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.1) ISEG3=1 201 : IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.NPFAC-1) ISEG3=1 202 : IAFA=LFAC(JFA) 203 : IBFA=LFAC(JFB) 204 : ICFA=LFAC(JFC) 205 : IA=IPT1.NUM(IAFA,JEL) 206 : IB=IPT1.NUM(IBFA,JEL) 207 : IC=IPT1.NUM(ICFA,JEL) 208 : XA=XCORD(1,ICPR(IA)) 209 : YA=XCORD(2,ICPR(IA)) 210 : ZA=XCORD(3,ICPR(IA)) 211 : XB=XCORD(1,ICPR(IB)) 212 : YB=XCORD(2,ICPR(IB)) 213 : ZB=XCORD(3,ICPR(IB)) 214 : XC=XCORD(1,ICPR(IC)) 215 : YC=XCORD(2,ICPR(IC)) 216 : ZC=XCORD(3,ICPR(IC)) 217 : C TESTONS SI LES POINTS SONT DU BON COTE 218 : SIGA=(XA-XCT)*XN+(YA-YCT)*YN+(ZA-ZCT)*ZN 219 : IF (ABS(SIGA).LE.1E-30) SIGA=1E-30 220 : SIGB=(XB-XCT)*XN+(YB-YCT)*YN+(ZB-ZCT)*ZN 221 : IF (ABS(SIGB).LE.1E-30) SIGB=1E-30 222 : SIGC=(XC-XCT)*XN+(YC-YCT)*YN+(ZC-ZCT)*ZN 223 : IF (ABS(SIGC).LE.1E-30) SIGC=1E-30 224 : if (isect.eq.1) goto 350 225 : IF (SIGA.GT.0..OR.SIGB.GT.0..OR.SIGC.GT.0.) GOTO 103 226 : IF (VOLUM) THEN 227 : IA1=MIN(IA,IB,IC) 228 : IB1=MAX(IA,IB,IC) 229 : IC1=IA+IB+IC 230 : DO 102 IDJF=ITRT(/1)-4,1,-5 231 : ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2)) 232 : ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2)) 233 : ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2) 234 : IF (IA1.NE.ITR1) GOTO 102 235 : IF (IB1.NE.ITR2) GOTO 102 236 : IF (IC1.NE.ITR3) GOTO 102 237 : * WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 1 ',IA,IB,IC 238 : C LE TRIANGLE EST DECLARE NON TRACABLE 239 : ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16 240 : GOTO 301 241 : 102 CONTINUE 242 : ENDIF 243 : * WRITE (6,*) ' NOUVEAU TRIANGLE ',IA,IB,IC 244 : ITRT(**)=IA 245 : ITRT(**)=IB 246 : ITRT(**)=IC 247 : ITRT(**)=IPT1.ICOLOR(JEL) 248 : ITRT(**)=ISEG1+2*ISEG2+4*ISEG3 249 : if (mcham.ne.0) then 250 : xtrt(**)=melva1.velche(min(lva1,iafa),jel) 251 : xtrt(**)=melva1.velche(min(lva1,ibfa),jel) 252 : xtrt(**)=melva1.velche(min(lva1,icfa),jel) 253 : endif 254 : 301 CONTINUE 255 : GOTO 101 256 : 103 CONTINUE 257 : IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0.) GOTO 101 258 : IPOS=0 259 : IF (SIGA.LE.0.) THEN 260 : IPOS=IPOS+1 261 : JTT(IPOS)=IA 262 : ISEG(IPOS)=ISEG1 263 : if (mcham.ne.0) xsegm(ia)= 264 : > melva1.velche(min(lva1,iafa),jel) 265 : ENDIF 266 : IF (SIGA*SIGB.LT.0..AND.(.NOT.(SIGA.LE.0..AND.SIGB.LE.0.))) THEN 267 : IA1=MIN(IA,IB) 268 : IB1=MAX(IA,IB) 269 : DO 104 NSEGM=1,ISEGM(/1),3 270 : IF (IA1.NE.ISEGM(NSEGM)) GOTO 104 271 : IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 104 272 : ISUP=ISEGM(NSEGM+2) 273 : if (mcham.ne.0) 274 : > xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigb- 275 : > melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga) 276 : * WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP 277 : GOTO 105 278 : 104 CONTINUE 279 : ITE=ITE+1 280 : SEGADJ XCORD 281 : XCORD(1,ITE)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA) 282 : XCORD(2,ITE)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA) 283 : XCORD(3,ITE)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA) 284 : IF (VCHCA.NE.0) 285 : # VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA) 286 : ICPR(**)=ITE 287 : ISUP=ICPR(/1) 288 : ISEGM(**)=IA1 289 : ISEGM(**)=IB1 290 : ISEGM(**)=ISUP 291 : if (mcham.ne.0) 292 : > xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb- 293 : > melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga) 294 : 105 CONTINUE 295 : IPOS=IPOS+1 296 : JTT(IPOS)=ISUP 297 : IF (SIGB.LE.0.) THEN 298 : ISEG(IPOS)=ISEG1 299 : ELSE 300 : ISEG(IPOS)=1 301 : ENDIF 302 : ENDIF 303 : IF (SIGB.LE.0.) THEN 304 : IPOS=IPOS+1 305 : JTT(IPOS)=IB 306 : ISEG(IPOS)=ISEG2 307 : if (mcham.ne.0) xsegm(ib)= 308 : > melva1.velche(min(lva1,ibfa),jel) 309 : ENDIF 310 : IF (SIGB*SIGC.LT.0..AND.(.NOT.(SIGB.LE.0..AND.SIGC.LE.0.))) THEN 311 : IB1=MIN(IB,IC) 312 : IC1=MAX(IB,IC) 313 : DO 106 NSEGM=1,ISEGM(/1),3 314 : IF (IB1.NE.ISEGM(NSEGM)) GOTO 106 315 : IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 106 316 : ISUP=ISEGM(NSEGM+2) 317 : if (mcham.ne.0) 318 : > xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigc- 319 : > melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb) 320 : * WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP 321 : GOTO 107 322 : 106 CONTINUE 323 : ITE=ITE+1 324 : SEGADJ XCORD 325 : XCORD(1,ITE)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB) 326 : XCORD(2,ITE)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB) 327 : XCORD(3,ITE)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB) 328 : IF (VCHCA.NE.0) 329 : # VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB) 330 : ICPR(**)=ITE 331 : ISUP=ICPR(/1) 332 : ISEGM(**)=IB1 333 : ISEGM(**)=IC1 334 : ISEGM(**)=ISUP 335 : if (mcham.ne.0) 336 : > xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigc- 337 : > melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb) 338 : 107 CONTINUE 339 : IPOS=IPOS+1 340 : JTT(IPOS)=ISUP 341 : IF (SIGC.LE.0.) THEN 342 : ISEG(IPOS)=ISEG2 343 : ELSE 344 : ISEG(IPOS)=1 345 : ENDIF 346 : ENDIF 347 : IF (SIGC.LE.0.) THEN 348 : IPOS=IPOS+1 349 : JTT(IPOS)=IC 350 : ISEG(IPOS)=ISEG3 351 : if (mcham.ne.0) xsegm(ic)= 352 : > melva1.velche(min(lva1,icfa),jel) 353 : ENDIF 354 : IF (SIGC*SIGA.LT.0..AND.(.NOT.(SIGC.LE.0..AND.SIGA.LE.0.))) THEN 355 : IC1=MIN(IC,IA) 356 : IA1=MAX(IC,IA) 357 : DO 108 NSEGM=1,ISEGM(/1),3 358 : IF (IC1.NE.ISEGM(NSEGM)) GOTO 108 359 : IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 108 360 : ISUP=ISEGM(NSEGM+2) 361 : if (mcham.ne.0) 362 : > xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*siga- 363 : > melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc) 364 : * WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP 365 : GOTO 109 366 : 108 CONTINUE 367 : ITE=ITE+1 368 : SEGADJ XCORD 369 : XCORD(1,ITE)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC) 370 : XCORD(2,ITE)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC) 371 : XCORD(3,ITE)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC) 372 : IF (VCHCA.NE.0) 373 : # VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC) 374 : ICPR(**)=ITE 375 : ISUP=ICPR(/1) 376 : ISEGM(**)=IC1 377 : ISEGM(**)=IA1 378 : ISEGM(**)=ISUP 379 : if (mcham.ne.0) 380 : > xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*siga- 381 : > melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc) 382 : 109 CONTINUE 383 : IPOS=IPOS+1 384 : JTT(IPOS)=ISUP 385 : IF (SIGA.LE.0.) THEN 386 : ISEG(IPOS)=ISEG3 387 : ELSE 388 : ISEG(IPOS)=1 389 : ENDIF 390 : ENDIF 391 : IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS 392 : IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS 393 : IF (IPOS.LT.3) GOTO 101 394 : IF (IPOS.EQ.3) THEN 395 : IF (VOLUM) THEN 396 : IA1=MIN(JTT(1),JTT(2),JTT(3)) 397 : IB1=MAX(JTT(1),JTT(2),JTT(3)) 398 : IC1=JTT(1)+JTT(2)+JTT(3) 399 : DO 110 IDJF=ITRT(/1)-4,1,-5 400 : ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2)) 401 : ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2)) 402 : ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2) 403 : IF (IA1.NE.ITR1) GOTO 110 404 : IF (IB1.NE.ITR2) GOTO 110 405 : IF (IC1.NE.ITR3) GOTO 110 406 : * WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 2 ',JTT(1),JTT(2),JTT(3) 407 : C LE TRIANGLE EST DECLARE NON TRACABLE 408 : ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16 409 : GOTO 303 410 : 110 CONTINUE 411 : ENDIF 412 : * WRITE (6,*) ' NOUVEAU TRIANGLE ',JTT(1),JTT(2),JTT(3) 413 : ITRT(**)=JTT(1) 414 : ITRT(**)=JTT(2) 415 : ITRT(**)=JTT(3) 416 : ITRT(**)=IPT1.ICOLOR(JEL) 417 : ITRT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3) 418 : if (mcham.ne.0) then 419 : xtrt(**)=xsegm(jtt(1)) 420 : xtrt(**)=xsegm(jtt(2)) 421 : xtrt(**)=xsegm(jtt(3)) 422 : endif 423 : ENDIF 424 : 303 CONTINUE 425 : IF (IPOS.EQ.4) THEN 426 : IF (VOLUM) THEN 427 : IA1=MIN(JTT(1),JTT(2),JTT(3),JTT(4)) 428 : IB1=MAX(JTT(1),JTT(2),JTT(3),JTT(4)) 429 : IC1=JTT(1)+JTT(2)+JTT(3)+JTT(4) 430 : ID1=JTT(1)*JTT(2)*JTT(3)*JTT(4) 431 : DO 111 IDJF=IQUAT(/1)-5,1,-6 432 : ITR1=MIN(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3)) 433 : ITR2=MAX(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3)) 434 : ITR3=IQUAT(IDJF)+IQUAT(IDJF+1)+IQUAT(IDJF+2)+IQUAT(IDJF+3) 435 : ITR4=IQUAT(IDJF)*IQUAT(IDJF+1)*IQUAT(IDJF+2)*IQUAT(IDJF+3) 436 : IF (IA1.NE.ITR1) GOTO 111 437 : IF (IB1.NE.ITR2) GOTO 111 438 : IF (IC1.NE.ITR3) GOTO 111 439 : IF (ID1.NE.ITR4) GOTO 111 440 : * WRITE (6,*) ' QUADRILA EN DOUBLE ',JTT(1),JTT(2),JTT(3),JTT(4) 441 : C LE QUADRILATERE EST DECLARE NON TRACABLE 442 : IQUAT(IDJF+5)=MOD(IQUAT(IDJF+5),16)+16 443 : GOTO 302 444 : 111 CONTINUE 445 : ENDIF 446 : * WRITE (6,*) ' NOUVEAU QUADRILATERE ',JTT(1),JTT(2),JTT(3),JTT(4) 447 : IQUAT(**)=JTT(1) 448 : IQUAT(**)=JTT(2) 449 : IQUAT(**)=JTT(3) 450 : IQUAT(**)=JTT(4) 451 : IQUAT(**)=IPT1.ICOLOR(JEL) 452 : IQUAT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8*ISEG(4) 453 : if (mcham.ne.0) then 454 : xquat(**)=xsegm(jtt(1)) 455 : xquat(**)=xsegm(jtt(2)) 456 : xquat(**)=xsegm(jtt(3)) 457 : xquat(**)=xsegm(jtt(4)) 458 : endif 459 : ENDIF 460 : 302 CONTINUE 461 : GOTO 101 462 : 101 CONTINUE 463 : 350 CONTINUE 464 : IF (.NOT.VOLUM) GOTO 160 465 : C ETUDE POUR LE CAS DES VOLUMES DU TETRAEDRE S'APPUYANT SUR LE 466 : C TRIANGLE ET LE CENTRE DE GRAVITE 467 : C TETRAEDRE 468 : IF (SIGA.LT.0..AND.SIGB.LT.0..AND.SIGC.LT.0..AND.SIGRAV 469 : # .LT.0.) GOTO 160 470 : IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0..AND.SIGRAV 471 : # .GT.0.) GOTO 160 472 : IPOS=0 473 : ISEG(1)=0 474 : ISEG(2)=0 475 : ISEG(3)=0 476 : ISEG(4)=0 477 : IF (SIGRAV*SIGA.LE.0.) THEN 478 : IF (SIGA.EQ.0.) THEN 479 : ISUP=IA 480 : if (mcham.ne.0) 481 : > xsegm(isup)=melva1.velche(min(lva1,iafa),jel) 482 : GOTO 359 483 : ENDIF 484 : DO 358 NSEGM=1,NJGRAV 485 : IF (IA.NE.IJGRAV(1,NSEGM)) GOTO 358 486 : ISUP=IJGRAV(2,NSEGM) 487 : if (mcham.ne.0) 488 : > xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigrav- 489 : > chgrav*siga)/(sigrav-siga) 490 : GOTO 359 491 : 358 CONTINUE 492 : ITE=ITE+1 493 : SEGADJ XCORD 494 : XCORD(1,ITE)=(XA*SIGRAV-XGRAV*SIGA)/(SIGRAV-SIGA) 495 : XCORD(2,ITE)=(YA*SIGRAV-YGRAV*SIGA)/(SIGRAV-SIGA) 496 : XCORD(3,ITE)=(ZA*SIGRAV-ZGRAV*SIGA)/(SIGRAV-SIGA) 497 : IF (VCHCA.NE.0) 498 : # VCHCA(**)=(VCHCA(IA)*SIGRAV-VGRAV*SIGA)/(SIGRAV-SIGA) 499 : ICPR(**)=ITE 500 : ISUP=ICPR(/1) 501 : if (mcham.ne.0) 502 : > xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigrav- 503 : > chgrav*siga)/(sigrav-siga) 504 : NJGRAV=NJGRAV+1 505 : IJGRAV(1,NJGRAV)=IA 506 : IJGRAV(2,NJGRAV)=ISUP 507 : 359 CONTINUE 508 : IPOS=IPOS+1 509 : JTT(IPOS)=ISUP 510 : ENDIF 511 : IF (SIGA*SIGB.LT.0..AND. 512 : # (.NOT.(SIGRAV*SIGA.LE.0..AND.SIGRAV*SIGB.LE.0.))) THEN 513 : IA1=MIN(IA,IB) 514 : IB1=MAX(IA,IB) 515 : DO 404 NSEGM=1,ISEGM(/1),3 516 : IF (IA1.NE.ISEGM(NSEGM)) GOTO 404 517 : IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 404 518 : ISUP=ISEGM(NSEGM+2) 519 : * WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP 520 : if (mcham.ne.0) 521 : > xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigb- 522 : > melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga) 523 : GOTO 405 524 : 404 CONTINUE 525 : ITE=ITE+1 526 : SEGADJ XCORD 527 : XCORD(1,ITE)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA) 528 : XCORD(2,ITE)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA) 529 : XCORD(3,ITE)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA) 530 : IF (VCHCA.NE.0) 531 : # VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA) 532 : ICPR(**)=ITE 533 : ISUP=ICPR(/1) 534 : ISEGM(**)=IA1 535 : ISEGM(**)=IB1 536 : ISEGM(**)=ISUP 537 : if (mcham.ne.0) 538 : > xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb- 539 : > melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga) 540 : 405 CONTINUE 541 : IPOS=IPOS+1 542 : JTT(IPOS)=ISUP 543 : IF (SIGRAV*SIGB.GT.0.) THEN 544 : ISEG(IPOS)=1 545 : ENDIF 546 : ENDIF 547 : IF (SIGRAV*SIGB.LE.0.) THEN 548 : IF (SIGB.EQ.0.) THEN 549 : ISUP=IB 550 : if (mcham.ne.0) 551 : > xsegm(isup)=melva1.velche(min(lva1,ibfa),jel) 552 : GOTO 369 553 : ENDIF 554 : DO 368 NSEGM=1,NJGRAV 555 : IF (IB.NE.IJGRAV(1,NSEGM)) GOTO 368 556 : ISUP=IJGRAV(2,NSEGM) 557 : GOTO 369 558 : 368 CONTINUE 559 : ITE=ITE+1 560 : SEGADJ XCORD 561 : XCORD(1,ITE)=(XB*SIGRAV-XGRAV*SIGB)/(SIGRAV-SIGB) 562 : XCORD(2,ITE)=(YB*SIGRAV-YGRAV*SIGB)/(SIGRAV-SIGB) 563 : XCORD(3,ITE)=(ZB*SIGRAV-ZGRAV*SIGB)/(SIGRAV-SIGB) 564 : IF (VCHCA.NE.0) 565 : # VCHCA(**)=(VCHCA(IB)*SIGRAV-VGRAV*SIGB)/(SIGRAV-SIGB) 566 : ICPR(**)=ITE 567 : ISUP=ICPR(/1) 568 : if (mcham.ne.0) 569 : > xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigrav- 570 : > chgrav*sigb)/(sigrav-sigb) 571 : NJGRAV=NJGRAV+1 572 : IJGRAV(1,NJGRAV)=IB 573 : IJGRAV(2,NJGRAV)=ISUP 574 : 369 CONTINUE 575 : IPOS=IPOS+1 576 : JTT(IPOS)=ISUP 577 : if (mcham.ne.0) 578 : > xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigrav- 579 : > chgrav*sigb)/(sigrav-sigb) 580 : ENDIF 581 : IF (SIGB*SIGC.LT.0..AND. 582 : # (.NOT.(SIGRAV*SIGB.LE.0..AND.SIGRAV*SIGC.LE.0.))) THEN 583 : IB1=MIN(IB,IC) 584 : IC1=MAX(IB,IC) 585 : DO 406 NSEGM=1,ISEGM(/1),3 586 : IF (IB1.NE.ISEGM(NSEGM)) GOTO 406 587 : IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 406 588 : ISUP=ISEGM(NSEGM+2) 589 : * WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP 590 : if (mcham.ne.0) 591 : > xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigc- 592 : > melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb) 593 : GOTO 407 594 : 406 CONTINUE 595 : ITE=ITE+1 596 : SEGADJ XCORD 597 : XCORD(1,ITE)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB) 598 : XCORD(2,ITE)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB) 599 : XCORD(3,ITE)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB) 600 : IF (VCHCA.NE.0) 601 : # VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB) 602 : ICPR(**)=ITE 603 : ISUP=ICPR(/1) 604 : ISEGM(**)=IB1 605 : ISEGM(**)=IC1 606 : ISEGM(**)=ISUP 607 : if (mcham.ne.0) 608 : > xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigc- 609 : > melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb) 610 : 407 CONTINUE 611 : IPOS=IPOS+1 612 : JTT(IPOS)=ISUP 613 : IF (SIGRAV*SIGC.GT.0.) THEN 614 : ISEG(IPOS)=1 615 : ENDIF 616 : ENDIF 617 : IF (SIGRAV*SIGC.LE.0.) THEN 618 : IF (SIGC.EQ.0.) THEN 619 : ISUP=IC 620 : if (mcham.ne.0) 621 : > xsegm(isup)=melva1.velche(min(lva1,icfa),jel) 622 : GOTO 379 623 : ENDIF 624 : DO 378 NSEGM=1,NJGRAV 625 : IF (IC.NE.IJGRAV(1,NSEGM)) GOTO 378 626 : ISUP=IJGRAV(2,NSEGM) 627 : if (mcham.ne.0) 628 : > xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*sigrav- 629 : > chgrav*sigc)/(sigrav-sigc) 630 : GOTO 379 631 : 378 CONTINUE 632 : ITE=ITE+1 633 : SEGADJ XCORD 634 : XCORD(1,ITE)=(XC*SIGRAV-XGRAV*SIGC)/(SIGRAV-SIGC) 635 : XCORD(2,ITE)=(YC*SIGRAV-YGRAV*SIGC)/(SIGRAV-SIGC) 636 : XCORD(3,ITE)=(ZC*SIGRAV-ZGRAV*SIGC)/(SIGRAV-SIGC) 637 : IF (VCHCA.NE.0) 638 : # VCHCA(**)=(VCHCA(IC)*SIGRAV-VGRAV*SIGC)/(SIGRAV-SIGC) 639 : ICPR(**)=ITE 640 : ISUP=ICPR(/1) 641 : if (mcham.ne.0) 642 : > xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*sigrav- 643 : > chgrav*sigc)/(sigrav-sigc) 644 : NJGRAV=NJGRAV+1 645 : IJGRAV(1,NJGRAV)=IC 646 : IJGRAV(2,NJGRAV)=ISUP 647 : 379 CONTINUE 648 : IPOS=IPOS+1 649 : JTT(IPOS)=ISUP 650 : ENDIF 651 : IF (SIGC*SIGA.LT.0..AND. 652 : # (.NOT.(SIGRAV*SIGC.LE.0..AND.SIGRAV*SIGA.LE.0.))) THEN 653 : IC1=MIN(IC,IA) 654 : IA1=MAX(IC,IA) 655 : DO 408 NSEGM=1,ISEGM(/1),3 656 : IF (IC1.NE.ISEGM(NSEGM)) GOTO 408 657 : IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 408 658 : ISUP=ISEGM(NSEGM+2) 659 : * WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP 660 : if (mcham.ne.0) 661 : > xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*siga- 662 : > melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc) 663 : GOTO 409 664 : 408 CONTINUE 665 : ITE=ITE+1 666 : SEGADJ XCORD 667 : XCORD(1,ITE)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC) 668 : XCORD(2,ITE)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC) 669 : XCORD(3,ITE)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC) 670 : IF (VCHCA.NE.0) 671 : # VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC) 672 : ICPR(**)=ITE 673 : ISUP=ICPR(/1) 674 : ISEGM(**)=IC1 675 : ISEGM(**)=IA1 676 : ISEGM(**)=ISUP 677 : if (mcham.ne.0) 678 : > xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*siga- 679 : > melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc) 680 : 409 CONTINUE 681 : IPOS=IPOS+1 682 : JTT(IPOS)=ISUP 683 : IF (SIGRAV*SIGA.GT.0.) THEN 684 : ISEG(IPOS)=1 685 : ENDIF 686 : ENDIF 687 : IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS 688 : IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS 689 : IF (IPOS.LT.3) GOTO 160 690 : ITR(**)=JTT(1) 691 : ITR(**)=JTT(2) 692 : ITR(**)=JTT(3) 693 : ITR(**)=IPT1.ICOLOR(JEL) 694 : if (mcham.ne.0) then 695 : xtr(**)=xsegm(jtt(1)) 696 : xtr(**)=xsegm(jtt(2)) 697 : xtr(**)=xsegm(jtt(3)) 698 : endif 699 : * WRITE (6,*) ' TRIANGLE (CAS 1) DU PLAN DE COUPE ' 700 : * WRITE (6,*) JTT(1),JTT(2),JTT(3) 701 : C TRIANGLE DU PLAN DE COUPE 702 : IF (IPOS.EQ.3) THEN 703 : MCOUP(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8 704 : ELSE 705 : MCOUP(**)=ISEG(1)+2*ISEG(2)+4*0 +8 706 : ENDIF 707 : 403 CONTINUE 708 : IF (IPOS.GE.4) THEN 709 : ITR(**)=JTT(1) 710 : ITR(**)=JTT(3) 711 : ITR(**)=JTT(4) 712 : ITR(**)=IPT1.ICOLOR(JEL) 713 : if (mcham.ne.0) then 714 : xtr(**)=xsegm(jtt(1)) 715 : xtr(**)=xsegm(jtt(3)) 716 : xtr(**)=xsegm(jtt(4)) 717 : endif 718 : * WRITE (6,*) ' TRIANGLE (CAS 2) DU PLAN DE COUPE ' 719 : * WRITE (6,*) JTT(1),JTT(3),JTT(4) 720 : C TRIANGLE DU PLAN DE COUPE 721 : MCOUP(**)=0 +2*ISEG(3)+4*ISEG(4)+8 722 : ENDIF 723 : 402 CONTINUE 724 : 160 CONTINUE 725 : 161 CONTINUE 726 : 20 CONTINUE 727 : NBELEM=JELN 728 : SEGADJ IPT2 729 : ISU=ISU+1 730 : IPT8.LISOUS(ISU)=IPT2 731 : if (mcham.ne.0) then 732 : n1el=nbelem 733 : segadj melva2 734 : ipt8.lisref(isu)=melva2 735 : endif 736 : SEGDES IPT1 737 : 10 CONTINUE 738 : IF (NBSOU.NE.0) SEGDES MELEME 739 : 12 CONTINUE 740 : NBSOUS=ISU+1 741 : * IF (ITR(/1).EQ.0) NBSOUS=ISU 742 : NBELEM=0 743 : NBNN=0 744 : NBREF=0 745 : if (mcham.ne.0) nbref=nbsous 746 : SEGADJ IPT8 747 : MELEME=IPT8 748 : * IF (ITR(/1).EQ.0) RETURN 749 : * WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES DANS LE PLAN DE COUPE', 750 : * # ITR(/1)/4 751 : * WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES HORS LE PLAN DE COUPE', 752 : * # ITRT(/1)/5 753 : * WRITE(6,*) ' NOMBRE DE QUADRILA ENGENDREES HORS LE PLAN DE COUPE', 754 : * # IQUAT(/1)/6 755 : NBELEM=ITR(/1)/4+ITRT(/1)/5+(IQUAT(/1)/6)*2 756 : NBNN=3 757 : NBREF=0 758 : NBSOUS=0 759 : SEGINI IPT4 760 : IPT4.ITYPEL=4 761 : if (mcham.ne.0) then 762 : n1el=nbelem 763 : n1ptel=nbnn 764 : segini melva4 765 : endif 766 : ITRL1=ITR(/1)/4 767 : DO 200 IEL=1,ITRL1 768 : DO NBNN=1,3 769 : IPT4.NUM(NBNN,IEL)=ITR((IEL-1)*4+NBNN) 770 : ENDDO 771 : if (mcham.ne.0) then 772 : DO NBNN=1,3 773 : melva4.velche(nbnn,iel)=xtr((IEL-1)*3+nbnn) 774 : ENDDO 775 : ENDIF 776 : IPT4.ICOLOR(IEL)=ITR(IEL*4) 777 : 200 CONTINUE 778 : ITRL2=ITRT(/1)/5 779 : DO 202 IEL=1,ITRL2 780 : DO 203 NBNN=1,3 781 : IPT4.NUM(NBNN,IEL+ITRL1)=ITRT((IEL-1)*5+NBNN) 782 : if (mcham.ne.0) 783 : > melva4.velche(nbnn,iel+itrl1)=xtrt((IEL-1)*3+nbnn) 784 : 203 CONTINUE 785 : IPT4.ICOLOR(IEL+ITRL1)=ITRT(IEL*5-1) 786 : MCOUP(**)=ITRT(IEL*5) 787 : 202 CONTINUE 788 : ITRL3=(IQUAT(/1)/6) 789 : DO 204 IEL=1,ITRL3 790 : IEL1=2*IEL-1 791 : DO 205 NBNN=1,3 792 : IPT4.NUM(NBNN,IEL1+ITRL1+ITRL2)=IQUAT((IEL-1)*6+NBNN) 793 : if (mcham.ne.0) 794 : > melva4.velche(nbnn,iel1+itrl1+itrl2)=xquat((IEL-1)*4+nbnn) 795 : 205 CONTINUE 796 : IPT4.ICOLOR(IEL1+ITRL1+ITRL2)=IQUAT(IEL*6-1) 797 : MCOUP(**)=MOD(IQUAT(IEL*6),4)+16*(IQUAT(IEL*6)/16) 798 : IEL2=2*IEL 799 : IPT4.NUM(1,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+1) 800 : IPT4.NUM(2,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+3) 801 : IPT4.NUM(3,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+4) 802 : IPT4.ICOLOR(IEL2+ITRL1+ITRL2)=IQUAT(IEL*6-1) 803 : MCOUP(**)=0+MOD(IQUAT(IEL*6)/4,2)*2+MOD(IQUAT(IEL*6)/8,2)*4 804 : # +16*(IQUAT(IEL*6)/16) 805 : if (mcham.ne.0) then 806 : melva4.velche(1,iel2+itrl1+itrl2)=xquat((IEL-1)*4+1) 807 : melva4.velche(2,iel2+itrl1+itrl2)=xquat((IEL-1)*4+3) 808 : melva4.velche(3,iel2+itrl1+itrl2)=xquat((IEL-1)*4+4) 809 : endif 810 : 204 CONTINUE 811 : * WRITE (6,*) ' LISTE DES ELEMENTS RAJOUTES ' 812 : * DO 9834 J=1,NBELEM 813 : * WRITE (6,*) (IPT4.NUM(I,J),I=1,3) 814 : *9834 CONTINUE 815 : LISOUS(ISU+1)=IPT4 816 : if (mcham.ne.0) lisref(ISU+1)=melva4 817 : IF (MELEM2.NE.0) THEN 818 : IF (MEDSAU.EQ.0) THEN 819 : C ON COUPE MELEM2 820 : mcham=0 821 : MEDSAU=MELEME 822 : MEDCOU=MCOUP 823 : MELEME=MELEM2 824 : SEGSUP ITR,ITRT,IQUAT 825 : SEGINI MCOUP,ITR,ITRT,IQUAT 826 : GOTO 5000 827 : ELSE 828 : MELEM2=MELEME 829 : MELEME=MEDSAU 830 : MCOU2=MCOUP 831 : MCOUP=MEDCOU 832 : ENDIF 833 : ENDIF 834 : SEGSUP ITR,ISEGM,IJGRAV,ITRT,IQUAT 835 : if (mchamor.ne.0) segsup xtr,xtrt,xquat,xsegm 836 : END 837 : 838 : 839 : 840 : 841 : 842 : 843 : 844 : 845 : 846 : 847 : 848 : 849 :
© Cast3M 2003 - All rights reserved.
Disclaimer