C DEMAIT SOURCE PV 07/10/08 21:15:00 5896 C|-------------------------------------------------------------------| C| | C| PROGRAMME PRINCIPAL | C| MAIN DE DEMETER | C| | C|-------------------------------------------------------------------| C SUBROUTINE DEMAIT(IDCP,NPTBAS) C C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(a-h,o-z) -INC CCOPTIO -INC TDEMAIT SEGMENT IDCP(NPTCOM) REAL*8 XDEM LOGICAL REPONS,FACET,INTER,PINT,FERME,VAL,IN,DROIT, # SOLHEX,SOLPRI,SOLPYR,SOLTET,DIAGO,PLAN real*8 xval,xo(3),xa(3),xb(3),xc(3),epsi,epsj,xmesu,ymesu character*4 mcle(8) data mcle /'TCRI','CFAC','CDIS','TETR','EXPC','FINA','DIAC', > 'EXPF'/ ymesu=0 NPTORI=NPTMAX NFCORI=NFCMAX * flags type d'operation ipass=1 tcrit=3. cfacei=6.0 cfacet=cfacei cdist=0.125 tetrl=2.75 expcom=1.00 diacri=0.93 diacrd=diacri diacre=diacri expfac=sqrt(3.)/2. faccri=16 volcri=0.01 10 continue call lirmot(mcle,8,imot,0) if (imot.eq.1) then call lirree(xval,1,iret) tcrit=xval endif if (imot.eq.2) then call lirree(xval,1,iret) cfacei=xval endif if (imot.eq.3) then call lirree(xval,1,iret) cdist=xval endif if (imot.eq.4) then call lirree(xval,1,iret) tetrl=xval endif if (imot.eq.5) then call lirree(xval,1,iret) expcom=xval endif if (imot.eq.7) then call lirree(xval,1,iret) diacre=xval endif if (imot.eq.8) then call lirree(xval,1,iret) expfac=xval endif if (imot.ne.0) goto 10 C C INITIALISATION DU TABLEAU DES VOLUMES C ------------------------------------- NFTOT=IFUT(/1) NVTOT=IVOL(/2) DO 130 J=1,NVTOT DO 120 I=1,9 IVOL(I,J)=0 120 CONTINUE 130 CONTINUE C C C CONSTRUCTION DU TABLEAU NPF ( POINTS-FACETTES ) C ----------------------------------------------- DO 150 J=1,40 DO 150 I=1,NPTMAX NPF(J,I)=0 150 CONTINUE DO 141 J=1,NFCMAX DO 140 K=1,4 IP=NFC(K,J) IF (IP.EQ.0) GOTO 140 L=-NPF(40,IP)+1 IF (L.LE.0) CALL ERREUR(126) IF (IERR.NE.0) RETURN NPF(40,IP)=-L NPF(L,IP)=J 140 CONTINUE 141 CONTINUE DO 145 I=1,NPTMAX NPF(40,I)=MAX(0,NPF(40,I)) 145 CONTINUE C C RECHERCHE DE LA TAILLE MOYENNE DE MAILLE ASSOCIEE A C CHAQUE POINT ( 4EME COMPOSANTE DU POINT ) DO 190 I=1,NPTMAX DD=0. KK=0 DO 170 J=1,40 IF (NPF(J,I).EQ.0) GOTO 180 if=npf(j,i) nc=4 if (nfc(4,if).eq.0) nc=3 jp=nfc(nc,if) do 175 ic=1,nc ip=nfc(ic,if) XX=(XYZ(1,IP)-XYZ(1,JP))**2 YY=(XYZ(2,IP)-XYZ(2,JP))**2 ZZ=(XYZ(3,IP)-XYZ(3,JP))**2 DD=DD+SQRT(XX+YY+ZZ) KK=KK+1 jp=ip 175 continue * IP=ISUCC(NPF(J,I),I) * XX=(XYZ(1,I)-XYZ(1,IP))**2 * YY=(XYZ(2,I)-XYZ(2,IP))**2 * ZZ=(XYZ(3,I)-XYZ(3,IP))**2 * DD=DD+SQRT(XX+YY+ZZ) * KK=KK+1 170 CONTINUE 180 XYZ(4,I)=DD/KK 190 CONTINUE * regularisation locale * DO 182 I=1,NPTMAX * DD= XYZ(4,I) * KK=1 * DO 184 J=1,40 * IF (NPF(J,I).EQ.0) GOTO 186 * IP=ISUCC(NPF(J,I),I) * DD=DD+XYZ(4,IP) * KK=KK+1 *84 CONTINUE *86 XYZ(4,I)=DD/KK *82 CONTINUE * taille moyenne generale xmoy=0 DO 181 I=1,NPTMAX XMOY=XMOY+LOG(XYZ(4,I)) 181 CONTINUE XMOYG=EXP(XMOY/NPTMAX) WRITE (6,*) ' TAILLE MOYENNE VISEE ',XMOYG C C LE MAILLAGE DE LA SURFACE EST-IL FERME ? C ---------------------------------------- REPONS=FERME(KKK) IF (REPONS) GOTO 210 CALL ERREUR(127) RETURN 210 continue xmesu=0 do 100 if=1,nfcmax xmesu=xmesu+vol(1,nfc(1,if),nfc(2,if),nfc(3,if)) if (nfc(4,if).ne.0) then xmesu=xmesu+vol(1,nfc(1,if),nfc(3,if),nfc(4,if)) ipass=2 endif 100 continue xmesu=xmesu/(-6.) WRITE (6,*) ' volume de la piece ',xmesu C C NFACET : NOMBRE DE FACETTES DU MAILLAGE DE SURFACE C -------------------------------------------------- NFACET=NFCMAX NPTCOM=NPTMAX NVOL=0 NFCPRE=0 NFCTRT=0 NARET=0 NPTOT=XYZ(/2) NTTRAV=NFCMAX SEGINI TRAV YMESU=0 NVOLY=0 NPTDEB=NPTMAX NPTDIS=1 ICRTS=0 220 CONTINUE do 222 jvol=nvoly+1,nvol if (ivol(9,jvol).eq.25) then ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), > ivol(3,jvol),ivol(4,jvol))/6. endif if (ivol(9,jvol).eq.35) then ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), > ivol(3,jvol),ivol(5,jvol))/6. > +vol(ivol(1,jvol),ivol(3,jvol), > ivol(4,jvol),ivol(5,jvol))/6. endif if (ivol(9,jvol).eq.30) then ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol), > ivol(3,jvol),ivol(4,jvol))/6. > +vol(ivol(2,jvol),ivol(3,jvol), > ivol(4,jvol),ivol(5,jvol))/6. > +vol(ivol(3,jvol),ivol(5,jvol), > ivol(6,jvol),ivol(4,jvol))/6. endif if (ivol(9,jvol).eq.20) then ymesu=ymesu-vol(ivol(1,jvol),ivol(3,jvol), > ivol(6,jvol),ivol(8,jvol))/6. > -vol(ivol(5,jvol),ivol(6,jvol), > ivol(8,jvol),ivol(1,jvol))/6. > -vol(ivol(2,jvol),ivol(6,jvol), > ivol(1,jvol),ivol(3,jvol))/6. ymesu=ymesu-vol(ivol(7,jvol),ivol(8,jvol), > ivol(6,jvol),ivol(3,jvol))/6. > -vol(ivol(4,jvol),ivol(1,jvol), > ivol(8,jvol),ivol(3,jvol))/6. endif 222 continue nvoly=nvol if (ymesu.gt.xmesu*1.01) goto 340 if (ierr.ne.0) goto 340 * AJUSTEMENT EVENTUEL DES DIMENSIONS DES TABLEAUX IF (NFCMAX+250.GE.NFTOT) THEN NFTOT=NFTOT+300 SEGADJ NFC,NFV,IFUT,IFAT ENDIF IF (NVOL+10.GE.NVTOT) THEN NVTOT=NVTOT+50 SEGADJ IVOL ENDIF IF (NPTMAX+50.GE.NPTOT) THEN NPTOT=NPTOT+100 SEGADJ NPF,XYZ ENDIF * nouvelle methode de calcul de la taille locale DO 221 I=NPTDEB+1,NPTMAX call vcrit(i) 221 CONTINUE NPTDEB=NPTMAX IGAGNE=0 IF (NFACET.EQ.0) GOTO 370 NFCMA=NFCMAX C C C RECHERCHE DES DIEDRES C --------------------- * FAIRE ICI LE MENAGE DANS NARET(de temps en temps) IF (ICRTS.GE.100) THEN NPTDIS=NPTMAX ICRTS=0 IVA=0 DO 285 I=1,NARET if (IIGARD(I).LE.0) goto 285 IF1=IF1GAR(I) IF (IFAT(IF1).EQ.0) GOTO 285 IF2=IF2GAR(I) IF (IFAT(IF2).EQ.0) GOTO 285 IVA=IVA+1 II=IIGARD(I) NPTDIS=MIN(NPTDIS,II) IIGARD(IVA)=II JJ=JJGARD(I) NPTDIS=MIN(NPTDIS,JJ) JJGARD(IVA)=JJ IF1=IF1GAR(I) IF1GAR(IVA)=IF1 IF2=IF2GAR(I) IF2GAR(IVA)=IF2 ANGAR(IVA)=ANGAR(I) IF (IVA.NE.1) THEN ANGMA(IVA)=MAX(ANGAR(IVA),ANGMA(IVA-1)) ELSE ANGMA(1)=ANGAR(1) ENDIF 285 CONTINUE * write (6,*) ' demait retassement effectue ',naret,iva NARET=IVA ENDIF DO 290 IF1=NFCTRT+1,NFCMAX IF (IFAT(IF1).EQ.0) GOTO 290 NBD=4 IF (NFC(4,IF1).EQ.0) NBD=3 DO 292 IC=1,NBD IC1=IC-1 IF (IC1.EQ.0) IC1=NBD II=NFC(IC1,IF1) JJ=NFC(IC,IF1) DO 294 I=1,40 IF2=NPF(I,II) IF (IF2.EQ.0) GOTO 292 IF (IF2.GE.IF1) GOTO 294 IF (IPRED(IF2,II).NE.JJ) GOTO 294 C C COMMENT SONT LES ELEMENTS DU DIEDRE ? C ------------------------------------- C ANGLL=TETA(IF1,IF2,II,JJ) C * pour gagner du temps if (angll.le.-1d4) goto 294 * if (if1.le.nfcori.or.if2.le.nfcori) angll=angll+1d6 NARET=NARET+1 IF (NARET.GT.NTTRAV) THEN NTTRAV=NARET+20 SEGADJ TRAV ENDIF IIGARD(NARET)=II JJGARD(NARET)=JJ IF1GAR(NARET)=IF1 IF2GAR(NARET)=IF2 ANGAR(NARET)=ANGLL IF (NARET.NE.1) THEN ANGMA(NARET)=MAX(ANGAR(NARET),ANGMA(NARET-1)) ELSE ANGMA(1)=ANGAR(1) ENDIF 294 CONTINUE 292 CONTINUE 290 CONTINUE NFCTRT=NFCMAX * * ON COMMENCE PAR L'ANGLE LE PLUS FERME * 315 CONTINUE ANGMAX=-1.E30 IOK=0 DO 310 I=NARET,1,-1 IF(ANGMAX.GE.ANGMA(I)) GOTO 311 IF(IIGARD(I).LE.0) GOTO 310 IF(ANGMAX.GE.ANGAR(I)) GOTO 310 IOK=I * if (angar(i).gt.0d0) ANGMAX=ANGAR(I)*0.99999D0 * if (angar(i).lt.0d0) ANGMAX=ANGAR(I)*1.00001D0 angmax=angar(i)-1d-6 310 CONTINUE 311 CONTINUE IF (IOK.EQ.0) GOTO 320 II=IIGARD(IOK) JJ=JJGARD(IOK) IF1=IF1GAR(IOK) IF2=IF2GAR(IOK) IIGARD(IOK)=-II ICRTS=ICRTS+1 IF (IFAT(IF1).EQ.0) GOTO 315 IF (IFAT(IF2).EQ.0) GOTO 315 IGAGNE=0 * write (6,*) 'traitement diedre ',ii,jj,if1,if2,angmax idiac=5 313 continue * on essaie d'abord de faire les hexaedres if (ipass.eq.2) then IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0) # CALL hexa(II,JJ,IF1,IF2,IGAGNE) GOTO 312 else IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).EQ.0) # CALL CONS33(II,JJ,IF1,IF2,IGAGNE,0) IF (IGAGNE.EQ.1) GOTO 312 IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).NE.0) # CALL CONS34(II,JJ,IF1,IF2,IGAGNE) IF (IGAGNE.EQ.1) GOTO 312 IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).EQ.0) # CALL CONS34(JJ,II,IF2,IF1,IGAGNE) IF (IGAGNE.EQ.1) GOTO 312 IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0) # CALL CONS44(II,JJ,IF1,IF2,IGAGNE) IF (IGAGNE.EQ.1) GOTO 312 endif * write (6,*) ' demait relachement de diacrd' diacrd=1-0.3*(1-diacrd) diacri=diacrd cfacet=cfacet*1.5 faccri=faccri*2 cdist=0.085 tetrl=tetrl*1.3 idiac=idiac-1 if (idiac.ne.0) goto 313 diacrd=diacre diacri=diacrd cfacet=cfacei faccri=16 cdist=0.125 tetrl=2.75 write (6,*) ' demait echec traitement diedre',ii,jj,if1,if2, > angmax goto 315 312 CONTINUE diacrd=diacre diacri=diacrd cfacet=cfacei faccri=16 cdist=0.125 tetrl=2.75 if (ipass.eq.-2) then ipass=-1 * cfacet=6 * cdist=0.125 * tetrl=2.75 * faccri=16 * volcri=0.2 * diacri=0.90 * diacrd=0.92 endif GOTO 220 320 CONTINUE IF (ipass.eq.2) then ipass=1 write (6,*) ' fin generation de cube ' * on remet les compteurs a zero nfctrt=0 naret =0 goto 220 endif IF (ipass.eq.1) then ipass=0 write (6,*) ' strategie finale 1' cfacet=16 cdist=0.085 tetrl=9 diacri=0.95 diacrd=0.95 faccri=64 volcri=0.01 * on remet les compteurs a zero nfctrt=0 naret =0 goto 220 endif IF (ipass.eq.0) then ipass=-1 write (6,*) ' strategie finale 2' cfacet=100 cdist=0.050 tetrl=9 diacri=0.99 diacrd=0.99 faccri=64 volcri=0.01 * on remet les compteurs a zero nfctrt=0 naret =0 goto 220 endif IF (ipass.eq.-1) then ipass=-2 write (6,*) ' strategie finale 3' cfacet=1000 cdist=0.005 tetrl=9 diacri=0.999 diacrd=0.999 faccri=64 volcri=0.01 * on remet les compteurs a zero nfctrt=0 naret =0 goto 220 endif 340 continue DO 444 I=1,NFCMAX IF (IFAT(I).EQ.0) GOTO 444 WRITE (6,*) ' FACETTE RESTANTE IF ',I,NFC(1,I),NFC(2,I), * NFC(3,I),NFC(4,I) 444 CONTINUE 4001 CONTINUE IF (NFACET.NE.0) CALL ERREUR(27) 370 CONTINUE WRITE (6,*) ' volume du maillage ',ymesu write (6,*) ' nb facette ',nfacet SEGSUP TRAV CALL OPTVOL RETURN C FIN DU PROGRAMME PRINCIPAL END