1 : C DEMAIT SOURCE PV 07/10/08 21:15:00 5896 2 : C|-------------------------------------------------------------------| 3 : C| | 4 : C| PROGRAMME PRINCIPAL | 5 : C| MAIN DE DEMETER | 6 : C| | 7 : C|-------------------------------------------------------------------| 8 : C 9 : SUBROUTINE DEMAIT(IDCP,NPTBAS) 10 : C 11 : C 12 : IMPLICIT INTEGER(I-N) 13 : IMPLICIT REAL*8(a-h,o-z) 14 : -INC CCOPTIO 15 : -INC TDEMAIT 16 : SEGMENT IDCP(NPTCOM) 17 : REAL*8 XDEM 18 : LOGICAL REPONS,FACET,INTER,PINT,FERME,VAL,IN,DROIT, 19 : # SOLHEX,SOLPRI,SOLPYR,SOLTET,DIAGO,PLAN 20 : real*8 xval,xo(3),xa(3),xb(3),xc(3),epsi,epsj,xmesu,ymesu 21 : character*4 mcle(8) 22 : data mcle /'TCRI','CFAC','CDIS','TETR','EXPC','FINA','DIAC', 23 : > 'EXPF'/ 24 : ymesu=0 25 : NPTORI=NPTMAX 26 : NFCORI=NFCMAX 27 : * flags type d'operation 28 : ipass=1 29 : tcrit=3. 30 : cfacei=6.0 31 : cfacet=cfacei 32 : cdist=0.125 33 : tetrl=2.75 34 : expcom=1.00 35 : diacri=0.93 36 : diacrd=diacri 37 : diacre=diacri 38 : expfac=sqrt(3.)/2. 39 : faccri=16 40 : volcri=0.01 41 : 10 continue 41 : call lirmot(mcle,8,imot,0) 42 : if (imot.eq.1) then 42 : call lirree(xval,1,iret) 43 : tcrit=xval 44 : endif 45 : if (imot.eq.2) then 45 : call lirree(xval,1,iret) 46 : cfacei=xval 47 : endif 48 : if (imot.eq.3) then 48 : call lirree(xval,1,iret) 49 : cdist=xval 50 : endif 51 : if (imot.eq.4) then 51 : call lirree(xval,1,iret) 52 : tetrl=xval 53 : endif 54 : if (imot.eq.5) then 54 : call lirree(xval,1,iret) 55 : expcom=xval 56 : endif 57 : if (imot.eq.7) then 57 : call lirree(xval,1,iret) 58 : diacre=xval 59 : endif 60 : if (imot.eq.8) then 60 : call lirree(xval,1,iret) 61 : expfac=xval 62 : endif 63 : if (imot.ne.0) goto 10 64 : C 65 : C INITIALISATION DU TABLEAU DES VOLUMES 66 : C ------------------------------------- 67 : NFTOT=IFUT(/1) 68 : NVTOT=IVOL(/2) 69 : DO 130 J=1,NVTOT 70 : DO 120 I=1,9 71 : IVOL(I,J)=0 72 : 120 CONTINUE 73 : 130 CONTINUE 74 : C 75 : C 76 : C CONSTRUCTION DU TABLEAU NPF ( POINTS-FACETTES ) 77 : C ----------------------------------------------- 78 : DO 150 J=1,40 79 : DO 150 I=1,NPTMAX 80 : NPF(J,I)=0 81 : 150 CONTINUE 82 : DO 141 J=1,NFCMAX 83 : DO 140 K=1,4 84 : IP=NFC(K,J) 85 : IF (IP.EQ.0) GOTO 140 86 : L=-NPF(40,IP)+1 86 : IF (L.LE.0) CALL ERREUR(126) 87 : IF (IERR.NE.0) RETURN 88 : NPF(40,IP)=-L 89 : NPF(L,IP)=J 90 : 140 CONTINUE 91 : 141 CONTINUE 92 : DO 145 I=1,NPTMAX 93 : NPF(40,I)=MAX(0,NPF(40,I)) 94 : 145 CONTINUE 95 : C 96 : C RECHERCHE DE LA TAILLE MOYENNE DE MAILLE ASSOCIEE A 97 : C CHAQUE POINT ( 4EME COMPOSANTE DU POINT ) 98 : DO 190 I=1,NPTMAX 99 : DD=0. 100 : KK=0 101 : DO 170 J=1,40 102 : IF (NPF(J,I).EQ.0) GOTO 180 103 : if=npf(j,i) 104 : nc=4 105 : if (nfc(4,if).eq.0) nc=3 106 : jp=nfc(nc,if) 107 : do 175 ic=1,nc 108 : ip=nfc(ic,if) 109 : XX=(XYZ(1,IP)-XYZ(1,JP))**2 110 : YY=(XYZ(2,IP)-XYZ(2,JP))**2 111 : ZZ=(XYZ(3,IP)-XYZ(3,JP))**2 112 : DD=DD+SQRT(XX+YY+ZZ) 113 : KK=KK+1 114 : jp=ip 115 : 175 continue 116 : * IP=ISUCC(NPF(J,I),I) 117 : * XX=(XYZ(1,I)-XYZ(1,IP))**2 118 : * YY=(XYZ(2,I)-XYZ(2,IP))**2 119 : * ZZ=(XYZ(3,I)-XYZ(3,IP))**2 120 : * DD=DD+SQRT(XX+YY+ZZ) 121 : * KK=KK+1 122 : 170 CONTINUE 123 : 180 XYZ(4,I)=DD/KK 124 : 190 CONTINUE 125 : * regularisation locale 126 : * DO 182 I=1,NPTMAX 127 : * DD= XYZ(4,I) 128 : * KK=1 129 : * DO 184 J=1,40 130 : * IF (NPF(J,I).EQ.0) GOTO 186 131 : * IP=ISUCC(NPF(J,I),I) 132 : * DD=DD+XYZ(4,IP) 133 : * KK=KK+1 134 : *84 CONTINUE 135 : *86 XYZ(4,I)=DD/KK 136 : *82 CONTINUE 137 : * taille moyenne generale 138 : xmoy=0 139 : DO 181 I=1,NPTMAX 140 : XMOY=XMOY+LOG(XYZ(4,I)) 141 : 181 CONTINUE 142 : XMOYG=EXP(XMOY/NPTMAX) 143 : WRITE (6,*) ' TAILLE MOYENNE VISEE ',XMOYG 144 : C 145 : C LE MAILLAGE DE LA SURFACE EST-IL FERME ? 146 : C ---------------------------------------- 147 : REPONS=FERME(KKK) 148 : IF (REPONS) GOTO 210 148 : CALL ERREUR(127) 149 : RETURN 150 : 210 continue 151 : xmesu=0 152 : do 100 if=1,nfcmax 153 : xmesu=xmesu+vol(1,nfc(1,if),nfc(2,if),nfc(3,if)) 154 : if (nfc(4,if).ne.0) then 155 : xmesu=xmesu+vol(1,nfc(1,if),nfc(3,if),nfc(4,if)) 156 : ipass=2 157 : endif 158 : 100 continue 159 : xmesu=xmesu/(-6.) 160 : WRITE (6,*) ' volume de la piece ',xmesu 161 : C 162 : C NFACET : NOMBRE DE FACETTES DU MAILLAGE DE SURFACE 163 : C -------------------------------------------------- 164 : NFACET=NFCMAX 165 : NPTCOM=NPTMAX 166 : NVOL=0 167 : NFCPRE=0 168 : NFCTRT=0 169 : NARET=0 170 : NPTOT=XYZ(/2) 171 : NTTRAV=NFCMAX 172 : SEGINI TRAV 173 : YMESU=0 174 : NVOLY=0 175 : NPTDEB=NPTMAX 176 : NPTDIS=1 177 : ICRTS=0 178 : 220 CONTINUE 179 : do 222 jvol=nvoly+1,nvol 180 : if (ivol(9,jvol).eq.25) then 181 : ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), 182 : > ivol(3,jvol),ivol(4,jvol))/6. 183 : endif 184 : if (ivol(9,jvol).eq.35) then 185 : ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), 186 : > ivol(3,jvol),ivol(5,jvol))/6. 187 : > +vol(ivol(1,jvol),ivol(3,jvol), 188 : > ivol(4,jvol),ivol(5,jvol))/6. 189 : endif 190 : if (ivol(9,jvol).eq.30) then 191 : ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), 192 : > ivol(3,jvol),ivol(4,jvol))/6. 193 : > +vol(ivol(2,jvol),ivol(3,jvol), 194 : > ivol(4,jvol),ivol(5,jvol))/6. 195 : > +vol(ivol(3,jvol),ivol(5,jvol), 196 : > ivol(6,jvol),ivol(4,jvol))/6. 197 : endif 198 : if (ivol(9,jvol).eq.20) then 199 : ymesu=ymesu-vol(ivol(1,jvol),ivol(3,jvol), 200 : > ivol(6,jvol),ivol(8,jvol))/6. 201 : > -vol(ivol(5,jvol),ivol(6,jvol), 202 : > ivol(8,jvol),ivol(1,jvol))/6. 203 : > -vol(ivol(2,jvol),ivol(6,jvol), 204 : > ivol(1,jvol),ivol(3,jvol))/6. 205 : ymesu=ymesu-vol(ivol(7,jvol),ivol(8,jvol), 206 : > ivol(6,jvol),ivol(3,jvol))/6. 207 : > -vol(ivol(4,jvol),ivol(1,jvol), 208 : > ivol(8,jvol),ivol(3,jvol))/6. 209 : endif 210 : 222 continue 211 : nvoly=nvol 212 : if (ymesu.gt.xmesu*1.01) goto 340 213 : if (ierr.ne.0) goto 340 214 : * AJUSTEMENT EVENTUEL DES DIMENSIONS DES TABLEAUX 215 : IF (NFCMAX+250.GE.NFTOT) THEN 216 : NFTOT=NFTOT+300 217 : SEGADJ NFC,NFV,IFUT,IFAT 218 : ENDIF 219 : IF (NVOL+10.GE.NVTOT) THEN 220 : NVTOT=NVTOT+50 221 : SEGADJ IVOL 222 : ENDIF 223 : IF (NPTMAX+50.GE.NPTOT) THEN 224 : NPTOT=NPTOT+100 225 : SEGADJ NPF,XYZ 226 : ENDIF 227 : * nouvelle methode de calcul de la taille locale 228 : DO 221 I=NPTDEB+1,NPTMAX 228 : call vcrit(i) 229 : 221 CONTINUE 230 : NPTDEB=NPTMAX 231 : IGAGNE=0 232 : IF (NFACET.EQ.0) GOTO 370 233 : NFCMA=NFCMAX 234 : C 235 : C 236 : C RECHERCHE DES DIEDRES 237 : C --------------------- 238 : * FAIRE ICI LE MENAGE DANS NARET(de temps en temps) 239 : IF (ICRTS.GE.100) THEN 240 : NPTDIS=NPTMAX 241 : ICRTS=0 242 : IVA=0 243 : DO 285 I=1,NARET 244 : if (IIGARD(I).LE.0) goto 285 245 : IF1=IF1GAR(I) 246 : IF (IFAT(IF1).EQ.0) GOTO 285 247 : IF2=IF2GAR(I) 248 : IF (IFAT(IF2).EQ.0) GOTO 285 249 : IVA=IVA+1 250 : II=IIGARD(I) 251 : NPTDIS=MIN(NPTDIS,II) 252 : IIGARD(IVA)=II 253 : JJ=JJGARD(I) 254 : NPTDIS=MIN(NPTDIS,JJ) 255 : JJGARD(IVA)=JJ 256 : IF1=IF1GAR(I) 257 : IF1GAR(IVA)=IF1 258 : IF2=IF2GAR(I) 259 : IF2GAR(IVA)=IF2 260 : ANGAR(IVA)=ANGAR(I) 261 : IF (IVA.NE.1) THEN 262 : ANGMA(IVA)=MAX(ANGAR(IVA),ANGMA(IVA-1)) 263 : ELSE 264 : ANGMA(1)=ANGAR(1) 265 : ENDIF 266 : 285 CONTINUE 267 : * write (6,*) ' demait retassement effectue ',naret,iva 268 : NARET=IVA 269 : ENDIF 270 : DO 290 IF1=NFCTRT+1,NFCMAX 271 : IF (IFAT(IF1).EQ.0) GOTO 290 272 : NBD=4 273 : IF (NFC(4,IF1).EQ.0) NBD=3 274 : DO 292 IC=1,NBD 275 : IC1=IC-1 276 : IF (IC1.EQ.0) IC1=NBD 277 : II=NFC(IC1,IF1) 278 : JJ=NFC(IC,IF1) 279 : DO 294 I=1,40 280 : IF2=NPF(I,II) 281 : IF (IF2.EQ.0) GOTO 292 282 : IF (IF2.GE.IF1) GOTO 294 283 : IF (IPRED(IF2,II).NE.JJ) GOTO 294 284 : C 285 : C COMMENT SONT LES ELEMENTS DU DIEDRE ? 286 : C ------------------------------------- 287 : C 288 : ANGLL=TETA(IF1,IF2,II,JJ) 289 : C 290 : * pour gagner du temps 291 : if (angll.le.-1d4) goto 294 292 : * if (if1.le.nfcori.or.if2.le.nfcori) angll=angll+1d6 293 : NARET=NARET+1 294 : IF (NARET.GT.NTTRAV) THEN 295 : NTTRAV=NARET+20 296 : SEGADJ TRAV 297 : ENDIF 298 : IIGARD(NARET)=II 299 : JJGARD(NARET)=JJ 300 : IF1GAR(NARET)=IF1 301 : IF2GAR(NARET)=IF2 302 : ANGAR(NARET)=ANGLL 303 : IF (NARET.NE.1) THEN 304 : ANGMA(NARET)=MAX(ANGAR(NARET),ANGMA(NARET-1)) 305 : ELSE 306 : ANGMA(1)=ANGAR(1) 307 : ENDIF 308 : 294 CONTINUE 309 : 292 CONTINUE 310 : 290 CONTINUE 311 : NFCTRT=NFCMAX 312 : * 313 : * ON COMMENCE PAR L'ANGLE LE PLUS FERME 314 : * 315 : 315 CONTINUE 316 : ANGMAX=-1.E30 317 : IOK=0 318 : DO 310 I=NARET,1,-1 319 : IF(ANGMAX.GE.ANGMA(I)) GOTO 311 320 : IF(IIGARD(I).LE.0) GOTO 310 321 : IF(ANGMAX.GE.ANGAR(I)) GOTO 310 322 : IOK=I 323 : * if (angar(i).gt.0d0) ANGMAX=ANGAR(I)*0.99999D0 324 : * if (angar(i).lt.0d0) ANGMAX=ANGAR(I)*1.00001D0 325 : angmax=angar(i)-1d-6 326 : 327 : 310 CONTINUE 328 : 311 CONTINUE 329 : IF (IOK.EQ.0) GOTO 320 330 : II=IIGARD(IOK) 331 : JJ=JJGARD(IOK) 332 : IF1=IF1GAR(IOK) 333 : IF2=IF2GAR(IOK) 334 : IIGARD(IOK)=-II 335 : ICRTS=ICRTS+1 336 : IF (IFAT(IF1).EQ.0) GOTO 315 337 : IF (IFAT(IF2).EQ.0) GOTO 315 338 : IGAGNE=0 339 : * write (6,*) 'traitement diedre ',ii,jj,if1,if2,angmax 340 : idiac=5 341 : 313 continue 342 : * on essaie d'abord de faire les hexaedres 343 : if (ipass.eq.2) then 344 : IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0) 344 : # CALL hexa(II,JJ,IF1,IF2,IGAGNE) 345 : GOTO 312 346 : else 347 : IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).EQ.0) 347 : # CALL CONS33(II,JJ,IF1,IF2,IGAGNE,0) 348 : IF (IGAGNE.EQ.1) GOTO 312 349 : IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).NE.0) 349 : # CALL CONS34(II,JJ,IF1,IF2,IGAGNE) 350 : IF (IGAGNE.EQ.1) GOTO 312 351 : IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).EQ.0) 351 : # CALL CONS34(JJ,II,IF2,IF1,IGAGNE) 352 : IF (IGAGNE.EQ.1) GOTO 312 353 : IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0) 353 : # CALL CONS44(II,JJ,IF1,IF2,IGAGNE) 354 : IF (IGAGNE.EQ.1) GOTO 312 355 : endif 356 : * write (6,*) ' demait relachement de diacrd' 357 : diacrd=1-0.3*(1-diacrd) 358 : diacri=diacrd 359 : cfacet=cfacet*1.5 360 : faccri=faccri*2 361 : cdist=0.085 362 : tetrl=tetrl*1.3 363 : idiac=idiac-1 364 : if (idiac.ne.0) goto 313 365 : diacrd=diacre 366 : diacri=diacrd 367 : cfacet=cfacei 368 : faccri=16 369 : cdist=0.125 370 : tetrl=2.75 371 : write (6,*) ' demait echec traitement diedre',ii,jj,if1,if2, 372 : > angmax 373 : goto 315 374 : 312 CONTINUE 375 : diacrd=diacre 376 : diacri=diacrd 377 : cfacet=cfacei 378 : faccri=16 379 : cdist=0.125 380 : tetrl=2.75 381 : if (ipass.eq.-2) then 382 : ipass=-1 383 : * cfacet=6 384 : * cdist=0.125 385 : * tetrl=2.75 386 : * faccri=16 387 : * volcri=0.2 388 : * diacri=0.90 389 : * diacrd=0.92 390 : endif 391 : GOTO 220 392 : 320 CONTINUE 393 : IF (ipass.eq.2) then 394 : ipass=1 395 : write (6,*) ' fin generation de cube ' 396 : * on remet les compteurs a zero 397 : nfctrt=0 398 : naret =0 399 : goto 220 400 : endif 401 : IF (ipass.eq.1) then 402 : ipass=0 403 : write (6,*) ' strategie finale 1' 404 : cfacet=16 405 : cdist=0.085 406 : tetrl=9 407 : diacri=0.95 408 : diacrd=0.95 409 : faccri=64 410 : volcri=0.01 411 : * on remet les compteurs a zero 412 : nfctrt=0 413 : naret =0 414 : goto 220 415 : endif 416 : IF (ipass.eq.0) then 417 : ipass=-1 418 : write (6,*) ' strategie finale 2' 419 : cfacet=100 420 : cdist=0.050 421 : tetrl=9 422 : diacri=0.99 423 : diacrd=0.99 424 : faccri=64 425 : volcri=0.01 426 : * on remet les compteurs a zero 427 : nfctrt=0 428 : naret =0 429 : goto 220 430 : endif 431 : IF (ipass.eq.-1) then 432 : ipass=-2 433 : write (6,*) ' strategie finale 3' 434 : cfacet=1000 435 : cdist=0.005 436 : tetrl=9 437 : diacri=0.999 438 : diacrd=0.999 439 : faccri=64 440 : volcri=0.01 441 : * on remet les compteurs a zero 442 : nfctrt=0 443 : naret =0 444 : goto 220 445 : endif 446 : 340 continue 447 : DO 444 I=1,NFCMAX 448 : IF (IFAT(I).EQ.0) GOTO 444 449 : WRITE (6,*) ' FACETTE RESTANTE IF ',I,NFC(1,I),NFC(2,I), 450 : * NFC(3,I),NFC(4,I) 451 : 444 CONTINUE 452 : 4001 CONTINUE 452 : IF (NFACET.NE.0) CALL ERREUR(27) 453 : 370 CONTINUE 454 : WRITE (6,*) ' volume du maillage ',ymesu 455 : write (6,*) ' nb facette ',nfacet 456 : SEGSUP TRAV 456 : CALL OPTVOL 457 : RETURN 458 : C FIN DU PROGRAMME PRINCIPAL 459 : END 460 : 461 : 462 : 463 : 464 : 465 : 466 :
© Cast3M 2003 - All rights reserved.
Disclaimer