Numérotation des lignes :

tetra
C TETRA     SOURCE    JC220346  16/11/29    21:15:38     9221           C---------------------------------------------------------------------|C                                                                     |       SUBROUTINE TETRA(II,JJ,IF1,IF2,IGAGNE,IBOUT)C                                                                     |C      CETTE SUBROUTINE TENTE DE CREER UN TETRAEDRE A PARTIR          |C      DES 2 TRIANGLES IF1 ET IF2.                                    |C                                                                     |C      -1- RECHERCHE DES FACETTES IF3 ET IF4                          |C      -2- CREATION DES FACETTES INEXISTENTES                         |C      -3- MISE A JOUR DU MAILLAGE DE SURFACE                         |C      -4- TEST DU VOLUME ELEMENTAIRE CREE                            |C      -5- MEMORISATION DU VOLUME CREE                                |C                                                                     |C      - IGAGNE=1 EN CAS DE SUCCES                                    |C      - IGAGNE=0 EN CAS D'ECHEC                                      |C                                                                     |C---------------------------------------------------------------------|C      IMPLICIT INTEGER(I-N)      IMPLICIT REAL*8(A-H,O-Z)-INC TDEMAIT -INC PPARAM-INC CCOPTIO       LOGICAL REPONS,FACET,SOLTET,DIAGO,IN2,VERDIV,IN       dimension ipts(10),jpts(10),qualti(10),qualtj(10)C       inouv=0       idepl=0       if (nfc(4,if1).ne.0) return       if (nfc(4,if2).ne.0) returnC*   on se refuse a mailler un tetraedre trop ouvert*       if (TETA(IF1,IF2,II,JJ).lt.-100) then*        write (6,*) 'if1 if2 ii jj TETA ',IF1,IF2,TETA(IF1,IF2,II,JJ)*        goto 155*       endif       nfcini=nfcmax       nptini=nptmax       ip1=0       ip2=0       N3=0       N4=0       ICTF=0       ICTV=0       ipas=0       IP=IPRED(IF1,II)       JP=ISUCC(IF2,II)*       if (TETA(IF1,IF2,II,JJ).lt.-100D0) goto 250CC      -1- RECHERCHE DES FACETTES IF3 ET IF4C      -------------------------------------C      RECHERCHE D'UNE 3EME FACETTE IF3C      --------------------------------C       IF3=IFACE3(IP,II,JP)         IF (IF3.NE.0) THEN*  si if3 dans le mauvais sens rien a faire       if (isucc(if3,ii).ne.ip) then        IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if3 '        return       endif*  SI IF3 QUADRANGULAIRE RENVOYER SUR PYRA1           IF (NFC(4,IF3).NE.0) THEN              CALL PYRA1(IP,II,IF1,IF3,IGAGNE,IBOUT)              RETURN           ENDIF           N3=IF3*         WRITE(6,1010)IF31010      FORMAT(' ** IF3=',I3)         ENDIFCCCC      RECHERCHE D'UNE 4EME FACETTE IF4C      --------------------------------C       IF4=IFACE3(IP,JJ,JP)         IF (IF4.NE.0) THEN*  si if4 dans le mauvais sens rien a faire       if (isucc(if4,jj).ne.jp) then        IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if4 '        return       endif           IF (NFC(4,IF4).NE.0) THEN*  SI IF4 QUADRANGULAIRE RENVOYER SUR PYRA1              CALL PYRA1(JJ,IP,IF1,IF4,IGAGNE,IBOUT)              RETURN           ENDIF           N4=IF4*          WRITE(6,1020)IF41020      FORMAT(' ** IF4=',I3)         ENDIFCC      -2- CREATION DES FACETTES INEXISTANTESC      --------------------------------------       IF (N3.EQ.0) THENCC      CREATION DE LA FACETTE IF3C      --------------------------       NFCMAX=NFCMAX+1       IF3=NFCMAX       ICTF=ICTF+1C       NFC(1,IF3)=II       NFC(2,IF3)=JP       NFC(3,IF3)=IP       NFC(4,IF3)=0C       ENDIFC       IF (N4.EQ.0) THENCC      CREATION DE LA FACETTE IF4C      --------------------------       NFCMAX=NFCMAX+1       IF4=NFCMAX       ICTF=ICTF+1C       NFC(1,IF4)=JJ       NFC(2,IF4)=IP       NFC(3,IF4)=JP       NFC(4,IF4)=0C       ENDIFCC  VERIFICATION VALIDITE ARETE CREE       IF (n3.EQ.0.AND.n4.EQ.0) THEN       IF (DIAGO(IP,JP,diacri)) then        iff=ifdiag        lpts=0        nfcmax=nfcini        if (iff.ne.0) then        if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)     >     then            ipt1=isucc(iff,ip)            if (ipt1.ne.ii.and.ipt1.ne.jj) then              lpts=lpts+1              ipts(lpts)=ipt1            endif            ipt1=ipred(iff,ip)            if (ipt1.ne.ii.and.ipt1.ne.jj) then              lpts=lpts+1              ipts(lpts)=ipt1            endif         endif        if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)     >     then            ipt1=isucc(iff,jp)            if (ipt1.ne.ii.and.ipt1.ne.jj) then              lpts=lpts+1              ipts(lpts)=ipt1            endif            ipt1=ipred(iff,jp)            if (ipt1.ne.ii.and.ipt1.ne.jj) then              lpts=lpts+1              ipts(lpts)=ipt1            endif         endif         lpts1=0         do 700 lpt=1,lpts          do 701 lpt1=1,lpts1            if (ipts(lpt).eq.jpts(lpt1)) goto 700 701      continue          qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),     >                    qualt(jj,ii,jp,ipts(lpt)))          if (qualti(lpt).le.0.001d0) goto 700          lpts1=lpts1+1          jpts(lpts1)=ipts(lpt)          qualtj(lpts1)=qualti(lpt) 700     continue         lpts=lpts1         lpts1=0         do 702 lpt=1,lpts          do 703 l=1,lpts1           if (qualtj(lpt).gt.qualti(l)) goto 704 703      continue          l=lpts1+1 704      continue          lpts1=lpts1+1          do 705 l1=lpts1,l+1,-1            qualti(l1)=qualti(l1-1)            ipts(l1)=ipts(l1-1) 705      continue            qualti(l)=qualtj(lpt)            ipts(l)=jpts(lpt) 702     continue         do 706 lpt=1,lpts*          write (6,*) 'tetra appel tetra2 type 1 avec ',*    >     ipts(lpt),iff,qualti(lpt)            call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))            if (igagne.eq.1) return 706     continue         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp         goto 350        endif*        write (6,*) ' tetra echec 1 lpts',lpts         return         ENDIF*  test in2*       IF (.NOT.IN2(ip,jp,nfcini)) THEN*         write (6,*) 'tetra test in2 echoue avec ',ip,jp*         nfcmax=nfcini*         nptmax=nptini*         return*       ENDIF       ENDIFCC      -3- MISE A JOUR DU MAILLAGE DE SURFACEC      --------------------------------------       CALL REPSUB(IF1)       CALL REPSUB(IF2)       CALL REPSUB(IF3)       CALL REPSUB(IF4)CC      -4- TEST DU VOLUME ELEMENTAIRE CREEC      -----------------------------------       IF (N3.EQ.0) THEN         IF (.NOT.FACET(IF3))     then*          write(6,*) ' tetra facet if3 invalide'           GOTO 150         endif       ENDIF       IF (N4.EQ.0) THEN         IF (.NOT.FACET(IF4))    then*          write(6,*) ' tetra facet if4 invalide'           GOTO 150         endif       ENDIF**  verification longueur de l'arete a creer*         DNORM=(XYZ(1,IP)-XYZ(1,JP))**2     #        +(XYZ(2,IP)-XYZ(2,JP))**2     #        +(XYZ(3,IP)-XYZ(3,JP))**2         DTEST=tetrl*XYZ(4,IP)*XYZ(4,JP)*        write (6,*) ' tetra dnorm dtest ',dnorm,dtest         IF (DNORM.GT.DTEST) THEN           CALL REPSUB(IF4)           CALL REPSUB(IF3)           CALL REPSUB(IF2)           CALL REPSUB(IF1)           NFCMAX=NFCini           GOTO 250         ENDIF*  verification validite du tetraedre       IF (.NOT.SOLTET(IF1,IF2,IF3,IF4,IPTT)) GOTO 151CC      LE VOLUME CREE EST VALIDE.C      --------------------------C      -5- MEMORISATION DU VOLUME CREEC      -------------------------------*       write (6,*) ' tetra creation tetraedre '       NVOL=NVOL+1       IF (NFV(1,IF1).EQ.0) NFV(1,IF1)=NVOL       IF (NFV(1,IF1).NE.NVOL) NFV(2,IF1)=NVOL       IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL       IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL       IF (NFV(1,IF3).EQ.0) NFV(1,IF3)=NVOL       IF (NFV(1,IF3).NE.NVOL) NFV(2,IF3)=NVOL       IF (NFV(1,IF4).EQ.0) NFV(1,IF4)=NVOL       IF (NFV(1,IF4).NE.NVOL) NFV(2,IF4)=NVOL       IVOL(9,NVOL)=25       DO 130 I=1,3          IVOL(I,NVOL)=NFC(I,IF1)130    CONTINUE       IVOL(4,NVOL)=JPC*      WRITE(6,1070)NVOL,(IVOL(I,NVOL),I=1,9)*1070   FORMAT(I3,4X,14I4)       do 961 npfi=40,1,-1        if (npf(npfi,ii).ne.0) goto 962 961   continue 962   continue       do 963 npfj=40,1,-1        if (npf(npfj,jj).ne.0) goto 964 963   continue 964   continue       qual=qualt(ii,jj,ip,jp)       if (iimpi.eq.1) write (6,1100) nvol,nfcmax,nfacet,qual,     >   (ivol(i,nvol),i=1,4) 1100  format (' TETRA  ',i5,2i6,f8.4,4i6)C*      DO 140 J=1,NPTMAX*          WRITE(6,1080)J,(NPF(I,J),I=1,40)*1080      FORMAT(I4,4X,40I3)* 140    CONTINUEC       IGAGNE=1C       RETURNC350    CONTINUE*  IP JP est invalide (diago)*  IP JP II JJ ont 4 facettes*  on va essayer de mettre un point au centre de*  gravite de ii jj ip jp et des 2 noeuds supplementaires*  des facettes touchant IP JP II JJ       IF3=NOISIN(IP,II,IF1)       if (nfc(4,if3).ne.0) return       IP1=ISUCC(IF3,IP)       IF5=NOISIN(II,JP,IF2)       if (nfc(4,if5).ne.0) return       IP3=ISUCC(IF5,II)       IF4=NOISIN(JP,JJ,IF2)       if (nfc(4,if4).ne.0) return       IP2=ISUCC(IF4,JP)       IF6=NOISIN(JJ,IP,IF1)       if (nfc(4,if6).ne.0) return       IP4=ISUCC(IF6,JJ)*      write (6,*) ' II JJ IP JP IP1 IP2 IP3 IP4',*    >   II,JJ,IP,JP,IP1,IP2,IP3,IP4       NFCMAX=NFCini       idepl=3       inouv=1       DO 351 I=1,4          XYZ(I,NPTMAX+1)=(XYZ(I,II)+XYZ(I,JJ)+XYZ(I,IP)+     >                   XYZ(I,JP)+(XYZ(I,IP1)+XYZ(I,IP3))/2+     >                   (XYZ(I,IP2)+XYZ(I,IP4))/2)/6 351   CONTINUE*      call vcrit(nptmax+1) 356   CONTINUE*  verif de validite de la position du point       iechec=0       v1=0.       v2=0.       v3=0.       v4=0.       v5=0.       v6=0.       v7=0.       v8=0.       v9=0.       v10=0.       v11=0.       v12=0.       v1=vol(nptmax+1,ii,ip,ip1)       if (v1.gt.0.) then       iechec=1*       write (6,*) ' tetra vol 1 positif ',v1       endif       v2=vol(nptmax+1,ii,ip1,jp)       if (v2.gt.0.) then       iechec=1*       write (6,*) ' tetra vol 2 positif ',v2       endif       if (ip1.ne.ip3) then        v3=vol(nptmax+1,ii,ip,ip3)        if (v3.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 3 positif ',v3        endif        v4=vol(nptmax+1,ii,ip3,jp)        if (v4.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 4 positif ',v4        endif       endif       v5=vol(nptmax+1,ip2,jj,jp)       if (v5.gt.0.) then       iechec=1*       write (6,*) ' tetra vol 5 positif ',v5       endif       v6=vol(nptmax+1,ip2,ip,jj)       if (v6.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 6 positif ',v6        endif        if (ip2.ne.ip4) then        v7=vol(nptmax+1,ip4,jj,jp)        if (v7.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 7 positif ',v7        endif        v8=vol(nptmax+1,ip4,ip,jj)        if (v8.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 8 positif ',v8        endif       endif        if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2) then        v9=vol(nptmax+1,ip2,ip1,ip)        if (v9.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 9 positif ',v9        endif        v10=vol(nptmax+1,ip1,ip2,jp)        if (v10.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 10 positif ',v10        endif       endif        v11=vol(nptmax+1,ii,jj,ip)        if (v11.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 11 positif ',v11        endif        v12=vol(nptmax+1,jj,ii,jp)        if (v12.gt.0.) then       iechec=1*        write (6,*) ' tetra vol 12 positif ',v12        endif       vtot=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12       if (iechec.eq.1) then*  deplacement du point        if (idepl.ne.0) then         xyz(1,nptmax+1)=1.10*xyz(1,nptmax+1)-0.10*(xyz(1,ii)+xyz(1,jj))/2         xyz(2,nptmax+1)=1.10*xyz(2,nptmax+1)-0.10*(xyz(2,ii)+xyz(2,jj))/2         xyz(3,nptmax+1)=1.10*xyz(3,nptmax+1)-0.10*(xyz(3,ii)+xyz(3,jj))/2         xyz(4,nptmax+1)=1.10*xyz(4,nptmax+1)-0.10*(xyz(4,ii)+xyz(4,jj))/2*        call vcrit(nptmax+1)*        write (6,*) ' tetra 1 deplacement du point ',nptmax+1         idepl=idepl-1         goto 356        endif*        goto 250*           write (6,*) ' tetra echec 3 pt plus deplacable'         return       endif*       write (6,*) ' tetra v* ',vtot,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10       if (v1.gt.vtot/1000) return       if (v2.gt.vtot/1000) return       if (ip1.ne.ip3.and.v3.gt.vtot/1000) return       if (ip1.ne.ip3.and.v4.gt.vtot/1000) return       if (v5.gt.vtot/1000) return       if (v6.gt.vtot/1000) return       if (ip2.ne.ip4.and.v7.gt.vtot/1000) return       if (ip2.ne.ip4.and.v8.gt.vtot/1000) return       if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.     >      v9.gt.vtot/1000) return       if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.     >      v10.gt.vtot/1000) return       if (v11.gt.vtot/1000) return       if (v12.gt.vtot/1000) return       NPTMAX=NPTMAX+1       GOTO 352250    CONTINUE*   IL FAUT FAIRE DEUX ELEMENTS       idepl=3       inouv=1       NFCMAX=NFCini       NPTMAX=NPTMAX+1*  CREATION POINT       xyz(4,nptmax)=(xyz(4,ii)+xyz(4,jj))/2*  deplacement du point pour l'eloigner de ii jj       xn1=(xyz(2,jj)-xyz(2,ii))*(xyz(3,ip)-xyz(3,ii))-     >     (xyz(3,jj)-xyz(3,ii))*(xyz(2,ip)-xyz(2,ii))       yn1=(xyz(3,jj)-xyz(3,ii))*(xyz(1,ip)-xyz(1,ii))-     >     (xyz(1,jj)-xyz(1,ii))*(xyz(3,ip)-xyz(3,ii))       zn1=(xyz(1,jj)-xyz(1,ii))*(xyz(2,ip)-xyz(2,ii))-     >     (xyz(2,jj)-xyz(2,ii))*(xyz(1,ip)-xyz(1,ii))       sn1=sqrt(xn1**2+yn1**2+zn1**2)       xn1=xn1/sn1       yn1=yn1/sn1       zn1=zn1/sn1       xn2=(xyz(2,jp)-xyz(2,ii))*(xyz(3,jj)-xyz(3,ii))-     >     (xyz(3,jp)-xyz(3,ii))*(xyz(2,jj)-xyz(2,ii))       yn2=(xyz(3,jp)-xyz(3,ii))*(xyz(1,jj)-xyz(1,ii))-     >     (xyz(1,jp)-xyz(1,ii))*(xyz(3,jj)-xyz(3,ii))       zn2=(xyz(1,jp)-xyz(1,ii))*(xyz(2,jj)-xyz(2,ii))-     >     (xyz(2,jp)-xyz(2,ii))*(xyz(1,jj)-xyz(1,ii))       sn2=sqrt(xn2**2+yn2**2+zn2**2)       xn2=xn2/sn2       yn2=yn2/sn2       zn2=zn2/sn2*       xn=(xn1+xn2)/2*       yn=(yn1+yn2)/2*       zn=(zn1+zn2)/2       xn=(xn1+xn2)/2       yn=(yn1+yn2)/2       zn=(zn1+zn2)/2       sn=sqrt(xn**2+yn**2+zn**2)       xn=xn/sn       yn=yn/sn       zn=zn/sn*      xmil=(xyz(1,ii)+xyz(1,jj))/2*      ymil=(xyz(2,ii)+xyz(2,jj))/2*      zmil=(xyz(3,ii)+xyz(3,jj))/2       xv=xyz(1,jj)-xyz(1,ii)       yv=xyz(2,jj)-xyz(2,ii)       zv=xyz(3,jj)-xyz(3,ii)       sv=sqrt(xv**2+yv**2+zv**2)       xv=xv/sv       yv=yv/sv       zv=zv/sv       xli=xv*(xyz(1,ip)-xyz(1,ii))+yv*(xyz(2,ip)-xyz(2,ii))+     >     zv*(xyz(3,ip)-xyz(3,ii))       xlj=xv*(xyz(1,jp)-xyz(1,ii))+yv*(xyz(2,jp)-xyz(2,ii))+     >     zv*(xyz(3,jp)-xyz(3,ii))*      xl=(xli+xlj+2*sv+2*0)/6       xl=0.5*sv*      xl=(xli+xlj)/2       xmil=xyz(1,ii)+xl*xv       ymil=xyz(2,ii)+xl*yv       zmil=xyz(3,ii)+xl*zv       expf = xyz(4,nptmax)       xyz(1,nptmax)=xmil+xn*expf*expfac       xyz(2,nptmax)=ymil+yn*expf*expfac       xyz(3,nptmax)=zmil+zn*expf*expfac*      call vcrit(nptmax)       IP1=0       IP2=0       IP3=0       IP4=0       goto 352 355   CONTINUE*  deplacement du point         xyz(1,nptmax)=0.70*xyz(1,nptmax)+0.30*(xyz(1,ii)+xyz(1,jj))/2         xyz(2,nptmax)=0.70*xyz(2,nptmax)+0.30*(xyz(2,ii)+xyz(2,jj))/2         xyz(3,nptmax)=0.70*xyz(3,nptmax)+0.30*(xyz(3,ii)+xyz(3,jj))/2         xyz(4,nptmax)=0.70*xyz(4,nptmax)+0.30*(xyz(4,ii)+xyz(4,jj))/2*        call vcrit(nptmax)*        write (6,*) ' tetra 2 deplacement du point ',nptmax         idepl=idepl-1 352   CONTINUE       IPTT=NPTMAX       CALL DIST(nptmax,nptaux,GL,IOK,ii,jj,ip,jp,ip1,ip2,ip3,ip4,0,0)       IF (IOK.EQ.0) THEN*        WRITE (6,*) ' tetra DIST ',nptaux         if (nptaux.eq.0) then           if (idepl.ne.0) goto 355           NPTMAX=nptini           NFCMAX=NFCini*           write (6,*) ' tetra echec 4 dist et pt non deplacable'           return         endif          if (idepl.ne.0) goto 355         NPTMAX=nptini         NFCMAX=NFCini         iptt=nptaux         idepl=0         inouv=0*        return       ELSEIF (gl.lt.xyz(4,nptmax)/4) then*        WRITE (6,*) ' tetra GL '          if (idepl.ne.0) goto 355         NPTMAX=nptini         NFCMAX=NFCini*           write (6,*) ' tetra echec 5 gl et pt non deplacable'         return       ENDIF 354   continue*  verification diago avant l'appel a tetra2        iff=0        lpts=0        if (diago(iptt,ii,diacrd)) iff=ifdiag        if (diago(iptt,jj,diacrd)) iff=ifdiag        if (diago(iptt,ip,diacrd)) iff=ifdiag        if (diago(iptt,jp,diacrd)) iff=ifdiag*        write (6,*) ' tetra iff ',iff        if (iff.ne.0) then           nptmaa=nptmax           nptmax=nptini           if (nfc(4,iff).ne.0) return        if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)     >     then            ipt1=isucc(iff,ip)            if (ipt1.ne.ii.and.ipt1.ne.jj.and.     >          ipt1.ne.jp.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif            ipt1=ipred(iff,ip)            if (ipt1.ne.ii.and.ipt1.ne.jj.and.     >          ipt1.ne.jp.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif         endif        if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)     >     then            ipt1=isucc(iff,jp)            if (ipt1.ne.ii.and.ipt1.ne.jj.and.     >          ipt1.ne.ip.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif            ipt1=ipred(iff,jp)            if (ipt1.ne.ii.and.ipt1.ne.jj.and.     >          ipt1.ne.ip.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif         endif        if (nfc(1,iff).eq.ii.or.nfc(2,iff).eq.ii.or.nfc(3,iff).eq.ii)     >     then            ipt1=isucc(iff,ii)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.jj.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif            ipt1=ipred(iff,ii)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.jj.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif         endif        if (nfc(1,iff).eq.jj.or.nfc(2,iff).eq.jj.or.nfc(3,iff).eq.jj)     >     then            ipt1=isucc(iff,jj)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.ii.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif            ipt1=ipred(iff,jj)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.ii.and.ipt1.ne.iptt) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif         endif        if (nfc(1,iff).eq.iptt.or.nfc(2,iff).eq.iptt.or.     >                            nfc(3,iff).eq.iptt)  then            ipt1=isucc(iff,iptt)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.ii.and.ipt1.ne.jj) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif            ipt1=ipred(iff,iptt)            if (ipt1.ne.ip.and.ipt1.ne.jp.and.     >          ipt1.ne.ii.and.ipt1.ne.jj) then                  lpts=lpts+1                  ipts(lpts)=ipt1            endif         endif         lpts1=0         lpts2=lpts         do 710 lpt=1,lpts          do 711 lpt1=1,lpts1            if (ipts(lpt).eq.jpts(lpt1)) goto 710 711      continue          qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),     >                    qualt(jj,ii,jp,ipts(lpt)))          if (qualti(lpt).le.0.001d0) goto 710          lpts1=lpts1+1          jpts(lpts1)=ipts(lpt)          qualtj(lpts1)=qualti(lpt) 710     continue         lpts=lpts1         lpts1=0         do 712 lpt=1,lpts          do 713 l=1,lpts1           if (qualtj(lpt).gt.qualti(l)) goto 714 713      continue          l=lpts1+1 714      continue          lpts1=lpts1+1          do 715 l1=lpts1,l+1,-1            qualti(l1)=qualti(l1-1)            ipts(l1)=ipts(l1-1) 715      continue            qualti(l)=qualtj(lpt)            ipts(l)=jpts(lpt) 712     continue         do 716 lpt=1,lpts*          write (6,*) 'tetra appel tetra2 type 2 avec ',*    >     ipts(lpt),iff,qualti(lpt)            call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))            if (igagne.eq.1) return 716     continue*           write (6,*) ' tetra echec 6 lpts2 lpts',lpts2,lpts,ip,jp*         return          nptmax=nptmaa        endif*  test supplementaire      if (inouv.eq.1) then       xyz(1,nptmax+1)=1.90*xyz(1,nptmax)-0.90*(xyz(1,ii)+xyz(1,jj))/2       xyz(2,nptmax+1)=1.90*xyz(2,nptmax)-0.90*(xyz(2,ii)+xyz(2,jj))/2       xyz(3,nptmax+1)=1.90*xyz(3,nptmax)-0.90*(xyz(3,ii)+xyz(3,jj))/2       xyz(4,nptmax+1)=1.90*xyz(4,nptmax)-0.90*(xyz(4,ii)+xyz(4,jj))/2*      call vcrit(nptmax+1)       IF (.NOT.IN2(ii,nptmax+1,nfcini)) THEN*  deplacement du point        if (idepl.ne.0) goto 355        NPTMAX=nptini        NFCMAX=NFCini*       write (6,*) ' tetra 2eme test in2 echoue ',ii,jj,iptt        RETURN       ENDIF      endif            quala=qualt(ii,jj,ip,iptt)            qualb=qualt(jj,ii,jp,iptt)*           write (6,*) ' tetra quala qualb ',quala,qualb       if (quala.gt.0.001.and.qualb.gt.0.001)     >       call tetra2(ii,jj,if1,if2,igagne,iptt)       if (igagne.eq.1) return*  deplacement du point        if (idepl.ne.0) goto 355       nptmax=nptini       idjf=1       goto 152150    CONTINUE       idjf=0CC      LE VOLUME CREE EST INVALIDE: IL FAUT DONC DETRUIRE LES FACETTESC      CREEES. -------------------------------------------------------       CALL REPSUB(IF4)       CALL REPSUB(IF3)       CALL REPSUB(IF2)       CALL REPSUB(IF1)C       NFCMAX=NFCini 152   CONTINUEC*   on essaye en cas d'intersection de facette d'utiliser cette information*      call facetx(k1,k2,k3,ifr)       k1=iirfac       k2=j1rfac       k3=j2rfac       if (k1.gt.nptmax) k1=0       if (k2.gt.nptmax) k2=0       if (k3.gt.nptmax) k3=0       ifr=ifrfac*       write (6,*) ' tetra retour de facetx ',k1,k2,k3,ifr       if (k1.eq.0) then*          write (6,*) ' tetra echec 7 facetx => 0 '          if (ifdiag.eq.0) return          k1=nfc(1,ifdiag)          k2=nfc(2,ifdiag)          k3=nfc(3,ifdiag)       endif       if (ip.eq.k1.or.jp.eq.k1) then         if (ifr.eq.if3.or.ifr.eq.0) then         iptt=0         if (jj.eq.k2) iptt=k3         if (jj.eq.k3) iptt=k2         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif         if (ifr.eq.if4.or.ifr.eq.0) then         iptt=0         if (ii.eq.k2) iptt=k3         if (ii.eq.k3) iptt=k2         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif       endif       if (ip.eq.k2.or.jp.eq.k2) then         if (ifr.eq.if3.or.ifr.eq.0) then         iptt=0         if (jj.eq.k1) iptt=k3         if (jj.eq.k3) iptt=k1         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif         if (ifr.eq.if4.or.ifr.eq.0) then         iptt=0         if (ii.eq.k1) iptt=k3         if (ii.eq.k3) iptt=k1         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif       endif       if (ip.eq.k3.or.jp.eq.k3) then         if (ifr.eq.if3.or.ifr.eq.0) then         iptt=0         if (jj.eq.k1) iptt=k2         if (jj.eq.k2) iptt=k1         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif         if (ifr.eq.if4.or.ifr.eq.0) then         iptt=0         if (ii.eq.k1) iptt=k2         if (ii.eq.k2) iptt=k1         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         endif       endif         if (idjf.eq.1) then*           write (6,*) ' tetra echec 8 apres facetx',k1,k2,k3,ifr            return         endif         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp          goto 350       RETURN151    CONTINUE       CALL REPSUB(IF4)       CALL REPSUB(IF3)       CALL REPSUB(IF2)       CALL REPSUB(IF1)       NFCMAX=NFCini         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)         if (igagne.eq.1) return         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp          goto 350       return 155   CONTINUE*       write (6,*) ' tetra probleme avec nouvelle arrete '*      CALL COM33(II,JJ,IF1,IF2,IGAGNE)C*       write (6,*) ' tetra echec 9'       RETURN       END     

© Cast3M 2003 - Tous droits réservés.
Mentions légales