Télécharger la plaquette Cast3M


fr en

Télécharger tetra.eso

Retour à la liste

Numérotation des lignes : 

tetra

  1. C TETRA     SOURCE    JC220346  16/11/29    21:15:38     9221           
  2. C---------------------------------------------------------------------|
  3. C                                                                     |
  4.        SUBROUTINE TETRA(II,JJ,IF1,IF2,IGAGNE,IBOUT)
  5. C                                                                     |
  6. C      CETTE SUBROUTINE TENTE DE CREER UN TETRAEDRE A PARTIR          |
  7. C      DES 2 TRIANGLES IF1 ET IF2.                                    |
  8. C                                                                     |
  9. C      -1- RECHERCHE DES FACETTES IF3 ET IF4                          |
  10. C      -2- CREATION DES FACETTES INEXISTENTES                         |
  11. C      -3- MISE A JOUR DU MAILLAGE DE SURFACE                         |
  12. C      -4- TEST DU VOLUME ELEMENTAIRE CREE                            |
  13. C      -5- MEMORISATION DU VOLUME CREE                                |
  14. C                                                                     |
  15. C      - IGAGNE=1 EN CAS DE SUCCES                                    |
  16. C      - IGAGNE=0 EN CAS D'ECHEC                                      |
  17. C                                                                     |
  18. C---------------------------------------------------------------------|
  19. C
  20.       IMPLICIT INTEGER(I-N)
  21.       IMPLICIT REAL*8(A-H,O-Z)
  22. -INC TDEMAIT
  23.  
  24. -INC PPARAM
  25. -INC CCOPTIO
  26.        LOGICAL REPONS,FACET,SOLTET,DIAGO,IN2,VERDIV,IN
  27.        dimension ipts(10),jpts(10),qualti(10),qualtj(10)
  28. C
  29.        inouv=0
  30.        idepl=0
  31.        if (nfc(4,if1).ne.0) return
  32.        if (nfc(4,if2).ne.0) return
  33. C
  34. *   on se refuse a mailler un tetraedre trop ouvert
  35. *       if (TETA(IF1,IF2,II,JJ).lt.-100) then
  36. *        write (6,*) 'if1 if2 ii jj TETA ',IF1,IF2,TETA(IF1,IF2,II,JJ)
  37. *        goto 155
  38. *       endif
  39.        nfcini=nfcmax
  40.        nptini=nptmax
  41.        ip1=0
  42.        ip2=0
  43.        N3=0
  44.        N4=0
  45.        ICTF=0
  46.        ICTV=0
  47.        ipas=0
  48.        IP=IPRED(IF1,II)
  49.        JP=ISUCC(IF2,II)
  50. *
  51.        if (TETA(IF1,IF2,II,JJ).lt.-100D0) goto 250
  52. C
  53. C      -1- RECHERCHE DES FACETTES IF3 ET IF4
  54. C      -------------------------------------
  55. C      RECHERCHE D'UNE 3EME FACETTE IF3
  56. C      --------------------------------
  57. C
  58.        IF3=IFACE3(IP,II,JP)
  59.          IF (IF3.NE.0) THEN
  60. *  si if3 dans le mauvais sens rien a faire
  61.        if (isucc(if3,ii).ne.ip) then
  62.         IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if3 '
  63.         return
  64.        endif
  65. *  SI IF3 QUADRANGULAIRE RENVOYER SUR PYRA1
  66.            IF (NFC(4,IF3).NE.0) THEN
  67.               CALL PYRA1(IP,II,IF1,IF3,IGAGNE,IBOUT)
  68.               RETURN
  69.            ENDIF
  70.            N3=IF3
  71. *         WRITE(6,1010)IF3
  72. 1010      FORMAT(' ** IF3=',I3)
  73.          ENDIF
  74. C
  75. C
  76. C
  77. C      RECHERCHE D'UNE 4EME FACETTE IF4
  78. C      --------------------------------
  79. C
  80.        IF4=IFACE3(IP,JJ,JP)
  81.          IF (IF4.NE.0) THEN
  82. *  si if4 dans le mauvais sens rien a faire
  83.        if (isucc(if4,jj).ne.jp) then
  84.         IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if4 '
  85.         return
  86.        endif
  87.            IF (NFC(4,IF4).NE.0) THEN
  88. *  SI IF4 QUADRANGULAIRE RENVOYER SUR PYRA1
  89.               CALL PYRA1(JJ,IP,IF1,IF4,IGAGNE,IBOUT)
  90.               RETURN
  91.            ENDIF
  92.            N4=IF4
  93. *          WRITE(6,1020)IF4
  94. 1020      FORMAT(' ** IF4=',I3)
  95.          ENDIF
  96. C
  97. C      -2- CREATION DES FACETTES INEXISTANTES
  98. C      --------------------------------------
  99.        IF (N3.EQ.0) THEN
  100. C
  101. C      CREATION DE LA FACETTE IF3
  102. C      --------------------------
  103.        NFCMAX=NFCMAX+1
  104.        IF3=NFCMAX
  105.        ICTF=ICTF+1
  106. C
  107.        NFC(1,IF3)=II
  108.        NFC(2,IF3)=JP
  109.        NFC(3,IF3)=IP
  110.        NFC(4,IF3)=0
  111. C
  112.        ENDIF
  113. C
  114.        IF (N4.EQ.0) THEN
  115. C
  116. C      CREATION DE LA FACETTE IF4
  117. C      --------------------------
  118.        NFCMAX=NFCMAX+1
  119.        IF4=NFCMAX
  120.        ICTF=ICTF+1
  121. C
  122.        NFC(1,IF4)=JJ
  123.        NFC(2,IF4)=IP
  124.        NFC(3,IF4)=JP
  125.        NFC(4,IF4)=0
  126. C
  127.        ENDIF
  128. C
  129. C  VERIFICATION VALIDITE ARETE CREE
  130.        IF (n3.EQ.0.AND.n4.EQ.0) THEN
  131.        IF (DIAGO(IP,JP,diacri)) then
  132.         iff=ifdiag
  133.         lpts=0
  134.         nfcmax=nfcini
  135.         if (iff.ne.0) then
  136.         if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)
  137.      >     then
  138.             ipt1=isucc(iff,ip)
  139.             if (ipt1.ne.ii.and.ipt1.ne.jj) then
  140.               lpts=lpts+1
  141.               ipts(lpts)=ipt1
  142.             endif
  143.             ipt1=ipred(iff,ip)
  144.             if (ipt1.ne.ii.and.ipt1.ne.jj) then
  145.               lpts=lpts+1
  146.               ipts(lpts)=ipt1
  147.             endif
  148.          endif
  149.         if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)
  150.      >     then
  151.             ipt1=isucc(iff,jp)
  152.             if (ipt1.ne.ii.and.ipt1.ne.jj) then
  153.               lpts=lpts+1
  154.               ipts(lpts)=ipt1
  155.             endif
  156.             ipt1=ipred(iff,jp)
  157.             if (ipt1.ne.ii.and.ipt1.ne.jj) then
  158.               lpts=lpts+1
  159.               ipts(lpts)=ipt1
  160.             endif
  161.          endif
  162.          lpts1=0
  163.          do 700 lpt=1,lpts
  164.           do 701 lpt1=1,lpts1
  165.             if (ipts(lpt).eq.jpts(lpt1)) goto 700
  166.  701      continue
  167.           qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),
  168.      >                    qualt(jj,ii,jp,ipts(lpt)))
  169.           if (qualti(lpt).le.0.001d0) goto 700
  170.           lpts1=lpts1+1
  171.           jpts(lpts1)=ipts(lpt)
  172.           qualtj(lpts1)=qualti(lpt)
  173.  700     continue
  174.          lpts=lpts1
  175.          lpts1=0
  176.          do 702 lpt=1,lpts
  177.           do 703 l=1,lpts1
  178.            if (qualtj(lpt).gt.qualti(l)) goto 704
  179.  703      continue
  180.           l=lpts1+1
  181.  704      continue
  182.           lpts1=lpts1+1
  183.           do 705 l1=lpts1,l+1,-1
  184.             qualti(l1)=qualti(l1-1)
  185.             ipts(l1)=ipts(l1-1)
  186.  705      continue
  187.             qualti(l)=qualtj(lpt)
  188.             ipts(l)=jpts(lpt)
  189.  702     continue
  190.          do 706 lpt=1,lpts
  191. *          write (6,*) 'tetra appel tetra2 type 1 avec ',
  192. *    >     ipts(lpt),iff,qualti(lpt)
  193.             call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))
  194.             if (igagne.eq.1) return
  195.  706     continue
  196.          if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
  197. *        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
  198.          goto 350
  199.         endif
  200. *        write (6,*) ' tetra echec 1 lpts',lpts
  201.          return
  202.          ENDIF
  203. *  test in2
  204. *       IF (.NOT.IN2(ip,jp,nfcini)) THEN
  205. *         write (6,*) 'tetra test in2 echoue avec ',ip,jp
  206. *         nfcmax=nfcini
  207. *         nptmax=nptini
  208. *         return
  209. *       ENDIF
  210.        ENDIF
  211. C
  212. C      -3- MISE A JOUR DU MAILLAGE DE SURFACE
  213. C      --------------------------------------
  214.        CALL REPSUB(IF1)
  215.        CALL REPSUB(IF2)
  216.        CALL REPSUB(IF3)
  217.        CALL REPSUB(IF4)
  218. C
  219. C      -4- TEST DU VOLUME ELEMENTAIRE CREE
  220. C      -----------------------------------
  221.        IF (N3.EQ.0) THEN
  222.          IF (.NOT.FACET(IF3))     then
  223. *          write(6,*) ' tetra facet if3 invalide'
  224.            GOTO 150
  225.          endif
  226.        ENDIF
  227.        IF (N4.EQ.0) THEN
  228.          IF (.NOT.FACET(IF4))    then
  229. *          write(6,*) ' tetra facet if4 invalide'
  230.            GOTO 150
  231.          endif
  232.        ENDIF
  233. *
  234. *  verification longueur de l'arete a creer
  235. *
  236.          DNORM=(XYZ(1,IP)-XYZ(1,JP))**2
  237.      #        +(XYZ(2,IP)-XYZ(2,JP))**2
  238.      #        +(XYZ(3,IP)-XYZ(3,JP))**2
  239.          DTEST=tetrl*XYZ(4,IP)*XYZ(4,JP)
  240. *        write (6,*) ' tetra dnorm dtest ',dnorm,dtest
  241.          IF (DNORM.GT.DTEST) THEN
  242.            CALL REPSUB(IF4)
  243.            CALL REPSUB(IF3)
  244.            CALL REPSUB(IF2)
  245.            CALL REPSUB(IF1)
  246.            NFCMAX=NFCini
  247.            GOTO 250
  248.          ENDIF
  249. *  verification validite du tetraedre
  250.        IF (.NOT.SOLTET(IF1,IF2,IF3,IF4,IPTT)) GOTO 151
  251. C
  252. C      LE VOLUME CREE EST VALIDE.
  253. C      --------------------------
  254. C      -5- MEMORISATION DU VOLUME CREE
  255. C      -------------------------------
  256. *       write (6,*) ' tetra creation tetraedre '
  257.        NVOL=NVOL+1
  258.        IF (NFV(1,IF1).EQ.0) NFV(1,IF1)=NVOL
  259.        IF (NFV(1,IF1).NE.NVOL) NFV(2,IF1)=NVOL
  260.        IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL
  261.        IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL
  262.        IF (NFV(1,IF3).EQ.0) NFV(1,IF3)=NVOL
  263.        IF (NFV(1,IF3).NE.NVOL) NFV(2,IF3)=NVOL
  264.        IF (NFV(1,IF4).EQ.0) NFV(1,IF4)=NVOL
  265.        IF (NFV(1,IF4).NE.NVOL) NFV(2,IF4)=NVOL
  266.        IVOL(9,NVOL)=25
  267.        DO 130 I=1,3
  268.           IVOL(I,NVOL)=NFC(I,IF1)
  269. 130    CONTINUE
  270.        IVOL(4,NVOL)=JP
  271. C
  272. *      WRITE(6,1070)NVOL,(IVOL(I,NVOL),I=1,9)
  273. *1070   FORMAT(I3,4X,14I4)
  274.        do 961 npfi=40,1,-1
  275.         if (npf(npfi,ii).ne.0) goto 962
  276.  961   continue
  277.  962   continue
  278.        do 963 npfj=40,1,-1
  279.         if (npf(npfj,jj).ne.0) goto 964
  280.  963   continue
  281.  964   continue
  282.        qual=qualt(ii,jj,ip,jp)
  283.        if (iimpi.eq.1) write (6,1100) nvol,nfcmax,nfacet,qual,
  284.      >   (ivol(i,nvol),i=1,4)
  285.  1100  format (' TETRA  ',i5,2i6,f8.4,4i6)
  286. C
  287. *      DO 140 J=1,NPTMAX
  288. *          WRITE(6,1080)J,(NPF(I,J),I=1,40)
  289. *1080      FORMAT(I4,4X,40I3)
  290. * 140    CONTINUE
  291. C
  292.        IGAGNE=1
  293. C
  294.        RETURN
  295. C
  296. 350    CONTINUE
  297. *  IP JP est invalide (diago)
  298. *  IP JP II JJ ont 4 facettes
  299. *  on va essayer de mettre un point au centre de
  300. *  gravite de ii jj ip jp et des 2 noeuds supplementaires
  301. *  des facettes touchant IP JP II JJ
  302.        IF3=NOISIN(IP,II,IF1)
  303.        if (nfc(4,if3).ne.0) return
  304.        IP1=ISUCC(IF3,IP)
  305.        IF5=NOISIN(II,JP,IF2)
  306.        if (nfc(4,if5).ne.0) return
  307.        IP3=ISUCC(IF5,II)
  308.        IF4=NOISIN(JP,JJ,IF2)
  309.        if (nfc(4,if4).ne.0) return
  310.        IP2=ISUCC(IF4,JP)
  311.        IF6=NOISIN(JJ,IP,IF1)
  312.        if (nfc(4,if6).ne.0) return
  313.        IP4=ISUCC(IF6,JJ)
  314. *      write (6,*) ' II JJ IP JP IP1 IP2 IP3 IP4',
  315. *    >   II,JJ,IP,JP,IP1,IP2,IP3,IP4
  316.        NFCMAX=NFCini
  317.        idepl=3
  318.        inouv=1
  319.        DO 351 I=1,4
  320.           XYZ(I,NPTMAX+1)=(XYZ(I,II)+XYZ(I,JJ)+XYZ(I,IP)+
  321.      >                   XYZ(I,JP)+(XYZ(I,IP1)+XYZ(I,IP3))/2+
  322.      >                   (XYZ(I,IP2)+XYZ(I,IP4))/2)/6
  323.  351   CONTINUE
  324. *      call vcrit(nptmax+1)
  325.  356   CONTINUE
  326. *  verif de validite de la position du point
  327.        iechec=0
  328.        v1=0.
  329.        v2=0.
  330.        v3=0.
  331.        v4=0.
  332.        v5=0.
  333.        v6=0.
  334.        v7=0.
  335.        v8=0.
  336.        v9=0.
  337.        v10=0.
  338.        v11=0.
  339.        v12=0.
  340.        v1=vol(nptmax+1,ii,ip,ip1)
  341.        if (v1.gt.0.) then
  342.        iechec=1
  343. *       write (6,*) ' tetra vol 1 positif ',v1
  344.        endif
  345.        v2=vol(nptmax+1,ii,ip1,jp)
  346.        if (v2.gt.0.) then
  347.        iechec=1
  348. *       write (6,*) ' tetra vol 2 positif ',v2
  349.        endif
  350.        if (ip1.ne.ip3) then
  351.         v3=vol(nptmax+1,ii,ip,ip3)
  352.         if (v3.gt.0.) then
  353.        iechec=1
  354. *        write (6,*) ' tetra vol 3 positif ',v3
  355.         endif
  356.         v4=vol(nptmax+1,ii,ip3,jp)
  357.         if (v4.gt.0.) then
  358.        iechec=1
  359. *        write (6,*) ' tetra vol 4 positif ',v4
  360.         endif
  361.        endif
  362.        v5=vol(nptmax+1,ip2,jj,jp)
  363.        if (v5.gt.0.) then
  364.        iechec=1
  365. *       write (6,*) ' tetra vol 5 positif ',v5
  366.        endif
  367.        v6=vol(nptmax+1,ip2,ip,jj)
  368.        if (v6.gt.0.) then
  369.        iechec=1
  370. *        write (6,*) ' tetra vol 6 positif ',v6
  371.         endif
  372.         if (ip2.ne.ip4) then
  373.         v7=vol(nptmax+1,ip4,jj,jp)
  374.         if (v7.gt.0.) then
  375.        iechec=1
  376. *        write (6,*) ' tetra vol 7 positif ',v7
  377.         endif
  378.         v8=vol(nptmax+1,ip4,ip,jj)
  379.         if (v8.gt.0.) then
  380.        iechec=1
  381. *        write (6,*) ' tetra vol 8 positif ',v8
  382.         endif
  383.        endif
  384.         if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2) then
  385.         v9=vol(nptmax+1,ip2,ip1,ip)
  386.         if (v9.gt.0.) then
  387.        iechec=1
  388. *        write (6,*) ' tetra vol 9 positif ',v9
  389.         endif
  390.         v10=vol(nptmax+1,ip1,ip2,jp)
  391.         if (v10.gt.0.) then
  392.        iechec=1
  393. *        write (6,*) ' tetra vol 10 positif ',v10
  394.         endif
  395.        endif
  396.         v11=vol(nptmax+1,ii,jj,ip)
  397.         if (v11.gt.0.) then
  398.        iechec=1
  399. *        write (6,*) ' tetra vol 11 positif ',v11
  400.         endif
  401.         v12=vol(nptmax+1,jj,ii,jp)
  402.         if (v12.gt.0.) then
  403.        iechec=1
  404. *        write (6,*) ' tetra vol 12 positif ',v12
  405.         endif
  406.        vtot=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12
  407.        if (iechec.eq.1) then
  408. *  deplacement du point
  409.         if (idepl.ne.0) then
  410.          xyz(1,nptmax+1)=1.10*xyz(1,nptmax+1)-0.10*(xyz(1,ii)+xyz(1,jj))/2
  411.          xyz(2,nptmax+1)=1.10*xyz(2,nptmax+1)-0.10*(xyz(2,ii)+xyz(2,jj))/2
  412.          xyz(3,nptmax+1)=1.10*xyz(3,nptmax+1)-0.10*(xyz(3,ii)+xyz(3,jj))/2
  413.          xyz(4,nptmax+1)=1.10*xyz(4,nptmax+1)-0.10*(xyz(4,ii)+xyz(4,jj))/2
  414. *        call vcrit(nptmax+1)
  415. *        write (6,*) ' tetra 1 deplacement du point ',nptmax+1
  416.          idepl=idepl-1
  417.          goto 356
  418.         endif
  419. *        goto 250
  420. *           write (6,*) ' tetra echec 3 pt plus deplacable'
  421.          return
  422.        endif
  423. *       write (6,*) ' tetra v* ',vtot,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10
  424.        if (v1.gt.vtot/1000) return
  425.        if (v2.gt.vtot/1000) return
  426.        if (ip1.ne.ip3.and.v3.gt.vtot/1000) return
  427.        if (ip1.ne.ip3.and.v4.gt.vtot/1000) return
  428.        if (v5.gt.vtot/1000) return
  429.        if (v6.gt.vtot/1000) return
  430.        if (ip2.ne.ip4.and.v7.gt.vtot/1000) return
  431.        if (ip2.ne.ip4.and.v8.gt.vtot/1000) return
  432.        if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.
  433.      >      v9.gt.vtot/1000) return
  434.        if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.
  435.      >      v10.gt.vtot/1000) return
  436.        if (v11.gt.vtot/1000) return
  437.        if (v12.gt.vtot/1000) return
  438.        NPTMAX=NPTMAX+1
  439.        GOTO 352
  440. 250    CONTINUE
  441. *   IL FAUT FAIRE DEUX ELEMENTS
  442.        idepl=3
  443.        inouv=1
  444.        NFCMAX=NFCini
  445.        NPTMAX=NPTMAX+1
  446. *  CREATION POINT
  447.        xyz(4,nptmax)=(xyz(4,ii)+xyz(4,jj))/2
  448. *  deplacement du point pour l'eloigner de ii jj
  449.        xn1=(xyz(2,jj)-xyz(2,ii))*(xyz(3,ip)-xyz(3,ii))-
  450.      >     (xyz(3,jj)-xyz(3,ii))*(xyz(2,ip)-xyz(2,ii))
  451.        yn1=(xyz(3,jj)-xyz(3,ii))*(xyz(1,ip)-xyz(1,ii))-
  452.      >     (xyz(1,jj)-xyz(1,ii))*(xyz(3,ip)-xyz(3,ii))
  453.        zn1=(xyz(1,jj)-xyz(1,ii))*(xyz(2,ip)-xyz(2,ii))-
  454.      >     (xyz(2,jj)-xyz(2,ii))*(xyz(1,ip)-xyz(1,ii))
  455.        sn1=sqrt(xn1**2+yn1**2+zn1**2)
  456.        xn1=xn1/sn1
  457.        yn1=yn1/sn1
  458.        zn1=zn1/sn1
  459.        xn2=(xyz(2,jp)-xyz(2,ii))*(xyz(3,jj)-xyz(3,ii))-
  460.      >     (xyz(3,jp)-xyz(3,ii))*(xyz(2,jj)-xyz(2,ii))
  461.        yn2=(xyz(3,jp)-xyz(3,ii))*(xyz(1,jj)-xyz(1,ii))-
  462.      >     (xyz(1,jp)-xyz(1,ii))*(xyz(3,jj)-xyz(3,ii))
  463.        zn2=(xyz(1,jp)-xyz(1,ii))*(xyz(2,jj)-xyz(2,ii))-
  464.      >     (xyz(2,jp)-xyz(2,ii))*(xyz(1,jj)-xyz(1,ii))
  465.        sn2=sqrt(xn2**2+yn2**2+zn2**2)
  466.        xn2=xn2/sn2
  467.        yn2=yn2/sn2
  468.        zn2=zn2/sn2
  469. *       xn=(xn1+xn2)/2
  470. *       yn=(yn1+yn2)/2
  471. *       zn=(zn1+zn2)/2
  472.        xn=(xn1+xn2)/2
  473.        yn=(yn1+yn2)/2
  474.        zn=(zn1+zn2)/2
  475.        sn=sqrt(xn**2+yn**2+zn**2)
  476.        xn=xn/sn
  477.        yn=yn/sn
  478.        zn=zn/sn
  479. *      xmil=(xyz(1,ii)+xyz(1,jj))/2
  480. *      ymil=(xyz(2,ii)+xyz(2,jj))/2
  481. *      zmil=(xyz(3,ii)+xyz(3,jj))/2
  482.        xv=xyz(1,jj)-xyz(1,ii)
  483.        yv=xyz(2,jj)-xyz(2,ii)
  484.        zv=xyz(3,jj)-xyz(3,ii)
  485.        sv=sqrt(xv**2+yv**2+zv**2)
  486.        xv=xv/sv
  487.        yv=yv/sv
  488.        zv=zv/sv
  489.        xli=xv*(xyz(1,ip)-xyz(1,ii))+yv*(xyz(2,ip)-xyz(2,ii))+
  490.      >     zv*(xyz(3,ip)-xyz(3,ii))
  491.        xlj=xv*(xyz(1,jp)-xyz(1,ii))+yv*(xyz(2,jp)-xyz(2,ii))+
  492.      >     zv*(xyz(3,jp)-xyz(3,ii))
  493. *      xl=(xli+xlj+2*sv+2*0)/6
  494.        xl=0.5*sv
  495. *      xl=(xli+xlj)/2
  496.        xmil=xyz(1,ii)+xl*xv
  497.        ymil=xyz(2,ii)+xl*yv
  498.        zmil=xyz(3,ii)+xl*zv
  499.        expf = xyz(4,nptmax)
  500.        xyz(1,nptmax)=xmil+xn*expf*expfac
  501.        xyz(2,nptmax)=ymil+yn*expf*expfac
  502.        xyz(3,nptmax)=zmil+zn*expf*expfac
  503. *      call vcrit(nptmax)
  504.        IP1=0
  505.        IP2=0
  506.        IP3=0
  507.        IP4=0
  508.        goto 352
  509.  355   CONTINUE
  510. *  deplacement du point
  511.          xyz(1,nptmax)=0.70*xyz(1,nptmax)+0.30*(xyz(1,ii)+xyz(1,jj))/2
  512.          xyz(2,nptmax)=0.70*xyz(2,nptmax)+0.30*(xyz(2,ii)+xyz(2,jj))/2
  513.          xyz(3,nptmax)=0.70*xyz(3,nptmax)+0.30*(xyz(3,ii)+xyz(3,jj))/2
  514.          xyz(4,nptmax)=0.70*xyz(4,nptmax)+0.30*(xyz(4,ii)+xyz(4,jj))/2
  515. *        call vcrit(nptmax)
  516. *        write (6,*) ' tetra 2 deplacement du point ',nptmax
  517.          idepl=idepl-1
  518.  352   CONTINUE
  519.        IPTT=NPTMAX
  520.        CALL DIST(nptmax,nptaux,GL,IOK,ii,jj,ip,jp,ip1,ip2,ip3,ip4,0,0)
  521.        IF (IOK.EQ.0) THEN
  522. *        WRITE (6,*) ' tetra DIST ',nptaux
  523.          if (nptaux.eq.0) then
  524.            if (idepl.ne.0) goto 355
  525.            NPTMAX=nptini
  526.            NFCMAX=NFCini
  527. *           write (6,*) ' tetra echec 4 dist et pt non deplacable'
  528.            return
  529.          endif
  530.           if (idepl.ne.0) goto 355
  531.          NPTMAX=nptini
  532.          NFCMAX=NFCini
  533.          iptt=nptaux
  534.          idepl=0
  535.          inouv=0
  536. *        return
  537.        ELSEIF (gl.lt.xyz(4,nptmax)/4) then
  538. *        WRITE (6,*) ' tetra GL '
  539.           if (idepl.ne.0) goto 355
  540.          NPTMAX=nptini
  541.          NFCMAX=NFCini
  542. *           write (6,*) ' tetra echec 5 gl et pt non deplacable'
  543.          return
  544.        ENDIF
  545.  354   continue
  546. *  verification diago avant l'appel a tetra2
  547.         iff=0
  548.         lpts=0
  549.         if (diago(iptt,ii,diacrd)) iff=ifdiag
  550.         if (diago(iptt,jj,diacrd)) iff=ifdiag
  551.         if (diago(iptt,ip,diacrd)) iff=ifdiag
  552.         if (diago(iptt,jp,diacrd)) iff=ifdiag
  553. *        write (6,*) ' tetra iff ',iff
  554.         if (iff.ne.0) then
  555.            nptmaa=nptmax
  556.            nptmax=nptini
  557.            if (nfc(4,iff).ne.0) return
  558.         if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)
  559.      >     then
  560.             ipt1=isucc(iff,ip)
  561.             if (ipt1.ne.ii.and.ipt1.ne.jj.and.
  562.      >          ipt1.ne.jp.and.ipt1.ne.iptt) then
  563.                   lpts=lpts+1
  564.                   ipts(lpts)=ipt1
  565.             endif
  566.             ipt1=ipred(iff,ip)
  567.             if (ipt1.ne.ii.and.ipt1.ne.jj.and.
  568.      >          ipt1.ne.jp.and.ipt1.ne.iptt) then
  569.                   lpts=lpts+1
  570.                   ipts(lpts)=ipt1
  571.             endif
  572.          endif
  573.         if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)
  574.      >     then
  575.             ipt1=isucc(iff,jp)
  576.             if (ipt1.ne.ii.and.ipt1.ne.jj.and.
  577.      >          ipt1.ne.ip.and.ipt1.ne.iptt) then
  578.                   lpts=lpts+1
  579.                   ipts(lpts)=ipt1
  580.             endif
  581.             ipt1=ipred(iff,jp)
  582.             if (ipt1.ne.ii.and.ipt1.ne.jj.and.
  583.      >          ipt1.ne.ip.and.ipt1.ne.iptt) then
  584.                   lpts=lpts+1
  585.                   ipts(lpts)=ipt1
  586.             endif
  587.          endif
  588.         if (nfc(1,iff).eq.ii.or.nfc(2,iff).eq.ii.or.nfc(3,iff).eq.ii)
  589.      >     then
  590.             ipt1=isucc(iff,ii)
  591.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  592.      >          ipt1.ne.jj.and.ipt1.ne.iptt) then
  593.                   lpts=lpts+1
  594.                   ipts(lpts)=ipt1
  595.             endif
  596.             ipt1=ipred(iff,ii)
  597.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  598.      >          ipt1.ne.jj.and.ipt1.ne.iptt) then
  599.                   lpts=lpts+1
  600.                   ipts(lpts)=ipt1
  601.             endif
  602.          endif
  603.         if (nfc(1,iff).eq.jj.or.nfc(2,iff).eq.jj.or.nfc(3,iff).eq.jj)
  604.      >     then
  605.             ipt1=isucc(iff,jj)
  606.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  607.      >          ipt1.ne.ii.and.ipt1.ne.iptt) then
  608.                   lpts=lpts+1
  609.                   ipts(lpts)=ipt1
  610.             endif
  611.             ipt1=ipred(iff,jj)
  612.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  613.      >          ipt1.ne.ii.and.ipt1.ne.iptt) then
  614.                   lpts=lpts+1
  615.                   ipts(lpts)=ipt1
  616.             endif
  617.          endif
  618.         if (nfc(1,iff).eq.iptt.or.nfc(2,iff).eq.iptt.or.
  619.      >                            nfc(3,iff).eq.iptt)  then
  620.             ipt1=isucc(iff,iptt)
  621.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  622.      >          ipt1.ne.ii.and.ipt1.ne.jj) then
  623.                   lpts=lpts+1
  624.                   ipts(lpts)=ipt1
  625.             endif
  626.             ipt1=ipred(iff,iptt)
  627.             if (ipt1.ne.ip.and.ipt1.ne.jp.and.
  628.      >          ipt1.ne.ii.and.ipt1.ne.jj) then
  629.                   lpts=lpts+1
  630.                   ipts(lpts)=ipt1
  631.             endif
  632.          endif
  633.          lpts1=0
  634.          lpts2=lpts
  635.          do 710 lpt=1,lpts
  636.           do 711 lpt1=1,lpts1
  637.             if (ipts(lpt).eq.jpts(lpt1)) goto 710
  638.  711      continue
  639.           qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),
  640.      >                    qualt(jj,ii,jp,ipts(lpt)))
  641.           if (qualti(lpt).le.0.001d0) goto 710
  642.           lpts1=lpts1+1
  643.           jpts(lpts1)=ipts(lpt)
  644.           qualtj(lpts1)=qualti(lpt)
  645.  710     continue
  646.          lpts=lpts1
  647.          lpts1=0
  648.          do 712 lpt=1,lpts
  649.           do 713 l=1,lpts1
  650.            if (qualtj(lpt).gt.qualti(l)) goto 714
  651.  713      continue
  652.           l=lpts1+1
  653.  714      continue
  654.           lpts1=lpts1+1
  655.           do 715 l1=lpts1,l+1,-1
  656.             qualti(l1)=qualti(l1-1)
  657.             ipts(l1)=ipts(l1-1)
  658.  715      continue
  659.             qualti(l)=qualtj(lpt)
  660.             ipts(l)=jpts(lpt)
  661.  712     continue
  662.          do 716 lpt=1,lpts
  663. *          write (6,*) 'tetra appel tetra2 type 2 avec ',
  664. *    >     ipts(lpt),iff,qualti(lpt)
  665.             call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))
  666.             if (igagne.eq.1) return
  667.  716     continue
  668. *           write (6,*) ' tetra echec 6 lpts2 lpts',lpts2,lpts,ip,jp
  669. *         return
  670.           nptmax=nptmaa
  671.         endif
  672. *  test supplementaire
  673.       if (inouv.eq.1) then
  674.        xyz(1,nptmax+1)=1.90*xyz(1,nptmax)-0.90*(xyz(1,ii)+xyz(1,jj))/2
  675.        xyz(2,nptmax+1)=1.90*xyz(2,nptmax)-0.90*(xyz(2,ii)+xyz(2,jj))/2
  676.        xyz(3,nptmax+1)=1.90*xyz(3,nptmax)-0.90*(xyz(3,ii)+xyz(3,jj))/2
  677.        xyz(4,nptmax+1)=1.90*xyz(4,nptmax)-0.90*(xyz(4,ii)+xyz(4,jj))/2
  678. *      call vcrit(nptmax+1)
  679.        IF (.NOT.IN2(ii,nptmax+1,nfcini)) THEN
  680. *  deplacement du point
  681.         if (idepl.ne.0) goto 355
  682.         NPTMAX=nptini
  683.         NFCMAX=NFCini
  684. *       write (6,*) ' tetra 2eme test in2 echoue ',ii,jj,iptt
  685.         RETURN
  686.        ENDIF
  687.       endif
  688.             quala=qualt(ii,jj,ip,iptt)
  689.             qualb=qualt(jj,ii,jp,iptt)
  690. *           write (6,*) ' tetra quala qualb ',quala,qualb
  691.        if (quala.gt.0.001.and.qualb.gt.0.001)
  692.      >       call tetra2(ii,jj,if1,if2,igagne,iptt)
  693.        if (igagne.eq.1) return
  694. *  deplacement du point
  695.         if (idepl.ne.0) goto 355
  696.        nptmax=nptini
  697.        idjf=1
  698.        goto 152
  699. 150    CONTINUE
  700.        idjf=0
  701. C
  702. C      LE VOLUME CREE EST INVALIDE: IL FAUT DONC DETRUIRE LES FACETTES
  703. C      CREEES. -------------------------------------------------------
  704.        CALL REPSUB(IF4)
  705.        CALL REPSUB(IF3)
  706.        CALL REPSUB(IF2)
  707.        CALL REPSUB(IF1)
  708. C
  709.        NFCMAX=NFCini
  710.  152   CONTINUE
  711. C
  712. *   on essaye en cas d'intersection de facette d'utiliser cette information
  713. *      call facetx(k1,k2,k3,ifr)
  714.        k1=iirfac
  715.        k2=j1rfac
  716.        k3=j2rfac
  717.        if (k1.gt.nptmax) k1=0
  718.        if (k2.gt.nptmax) k2=0
  719.        if (k3.gt.nptmax) k3=0
  720.        ifr=ifrfac
  721. *       write (6,*) ' tetra retour de facetx ',k1,k2,k3,ifr
  722.        if (k1.eq.0) then
  723. *          write (6,*) ' tetra echec 7 facetx => 0 '
  724.           if (ifdiag.eq.0) return
  725.           k1=nfc(1,ifdiag)
  726.           k2=nfc(2,ifdiag)
  727.           k3=nfc(3,ifdiag)
  728.        endif
  729.        if (ip.eq.k1.or.jp.eq.k1) then
  730.          if (ifr.eq.if3.or.ifr.eq.0) then
  731.          iptt=0
  732.          if (jj.eq.k2) iptt=k3
  733.          if (jj.eq.k3) iptt=k2
  734.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  735.          if (igagne.eq.1) return
  736.          endif
  737.          if (ifr.eq.if4.or.ifr.eq.0) then
  738.          iptt=0
  739.          if (ii.eq.k2) iptt=k3
  740.          if (ii.eq.k3) iptt=k2
  741.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  742.          if (igagne.eq.1) return
  743.          endif
  744.        endif
  745.        if (ip.eq.k2.or.jp.eq.k2) then
  746.          if (ifr.eq.if3.or.ifr.eq.0) then
  747.          iptt=0
  748.          if (jj.eq.k1) iptt=k3
  749.          if (jj.eq.k3) iptt=k1
  750.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  751.          if (igagne.eq.1) return
  752.          endif
  753.          if (ifr.eq.if4.or.ifr.eq.0) then
  754.          iptt=0
  755.          if (ii.eq.k1) iptt=k3
  756.          if (ii.eq.k3) iptt=k1
  757.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  758.          if (igagne.eq.1) return
  759.          endif
  760.        endif
  761.        if (ip.eq.k3.or.jp.eq.k3) then
  762.          if (ifr.eq.if3.or.ifr.eq.0) then
  763.          iptt=0
  764.          if (jj.eq.k1) iptt=k2
  765.          if (jj.eq.k2) iptt=k1
  766.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  767.          if (igagne.eq.1) return
  768.          endif
  769.          if (ifr.eq.if4.or.ifr.eq.0) then
  770.          iptt=0
  771.          if (ii.eq.k1) iptt=k2
  772.          if (ii.eq.k2) iptt=k1
  773.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  774.          if (igagne.eq.1) return
  775.          endif
  776.        endif
  777.          if (idjf.eq.1) then
  778. *           write (6,*) ' tetra echec 8 apres facetx',k1,k2,k3,ifr
  779.             return
  780.          endif
  781.          if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
  782. *        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
  783.           goto 350
  784.        RETURN
  785. 151    CONTINUE
  786.        CALL REPSUB(IF4)
  787.        CALL REPSUB(IF3)
  788.        CALL REPSUB(IF2)
  789.        CALL REPSUB(IF1)
  790.        NFCMAX=NFCini
  791.          if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
  792.          if (igagne.eq.1) return
  793.          if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
  794. *        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
  795.           goto 350
  796.        return
  797.  155   CONTINUE
  798. *       write (6,*) ' tetra probleme avec nouvelle arrete '
  799. *      CALL COM33(II,JJ,IF1,IF2,IGAGNE)
  800. C
  801. *       write (6,*) ' tetra echec 9'
  802.        RETURN
  803.        END
  804.  
  805.  
  806.  
  807.  
  808.  

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