Sources Esope | Cast3M

trplsi
C TRPLSI    SOURCE    CHAT      06/03/29    21:36:47     5360
      SUBROUTINE TRPLSI(X,IPOLYG,NCP,
     >            IND,ICARD,JT,IFD,
     >            IFG,IT,ITM,IPERE,
     >            ITPOLY,ITRIA,
     >            ITRMIN,
     >            QTRIA,QMIN,NPMAX,NCMAX,
C     >            ITRPOL,FCRMIN,QTMIN,iarr)  modif TC
     >            ITRPOL,QTMIN,iarr)
C     *****************************************************************
C     OBJET : CALCULE LA TRIANGULATION D'UN POLYGONE SIMPLE
C             QUI MAXIMISE LA VALEUR MINIMUM D'UN CRITERE DONNE
C
C     EN ENTREE :
C        X     : COORDONNEES DES POINTS DU POLYGONE
C        IPOLYG : NUMERO DES NOEUDS DU POLYGONE
C        NCP   : NOMBRE DE POINT DU POLYGONE
C        FCRMIN: FONCTION RENVOYANT LA VALEUR DU CRITERE
C                FUNCTION REAL FCRMIN(P1,P2,P3)
C                REAL*8 P1(*),P2(*),P3(*)
C                OU P1,P2,P3 SONT LES COORDONNEES DES POINTS
C                DU TRIANGLE
C        QTMIN : VALEUR MINIMUM DU CRITERE
C
C     LES TABLEAUX DES POLYGONES
C     --------------------------
C     IND(IPI)  DONNE L'ADRESSE DU POLYGONE "IPI" DANS "ITPOLY"
C     ICARD(IPI) DONNE LE NOMBRE DE COTES DU POLYGONE "IPI"
C     JT(IPI)   DONNE L'INDICE DES TRIANGLES DEJA TESTES
C     IFG(IPI),IFD(IPI)  DONNENT LES ADRESSES DES POLYGONES
C                    FILS DROIT ET FILS GAUCHE
C     QTRIA(IPI) DONNE LA QUALITE DE LA TRIANGULATION COURANTE
C     QMIN(IPI)  DONNE LA QUALITE MINI DEJA ATTEINTE POUR UNE
C               TRIANGULATION DE IPI
C     IT(IPI)    L'ADRESSE DE LA TRIANGULATION COURANTE DE IPI
C               LE NB DE TRIANGLES = ICARD(IPI) - 2
C     ITM(IPI)   L'ADRESSE DE LA TRIANGULATION MINI DE IPI
C     ----------------------------
C     ITPOLY : TABLEAU OU SONT DECRITS LES POLYGONES
C              LA LISTE DES INDICES DES NOEUDS
C     ITRIA   : TABLEAU OU SONT DECRITES LES TRIANGULATIONS
C     --------------------------------------------------------
C     NCMAX : LE NOMBRE MAXIMUM DE COTE DU POLYGONE
C     NPMAX : LE NOMBRE MAXIMUM DE POLYGONES EMPILES
C             POUR UN POLY DE N COTE = 2*(N-2) POLYGONES
C             DONT UN FILS PEUT AVOIR N-1 COTE.
C             DANS LE PIRE CAS NPMAX = (NCMAX-2)!
C             (NCMAX-2)!= 6.4E15
C             1000 > 7!
C      INTEGER   IND(NPMAX),ICARD(NPMAX),JT(NPMAX),IFD(NPMAX)
C      INTEGER   IFG(NPMAX),IT(NPMAX),ITM(NPMAX),IPERE(NPMAX)
C      INTEGER   ITPOLY(NPMAX*NCMAX),ITRIA((NCMAX-2)*3)
C      INTEGER   ITRMIN((NCMAX-2)*3*NPMAX)
C      REAL*8      QTRIA(NPMAX),QMIN(NPMAX)
C
C     EN SORTIE :
C        ITRPOL: TRIANGULATION RESULANTE
C                 ITRPOL((I-1)*3+1) PREMIER NOEUD DU TRIANGLE I
C                 ITRPOL((I-1)*3+2) DEUXIEME NOEUD DU TRIANGLE I
C                 ITRPOL((I-1)*3+3) TROISIEME NOEUD DU TRIANGLE I
C        QTMIN  : VALEUR MINIMUM DE FCRMIN SUR ITRPOL
C        iarr   : 0 SI TOUT EST OK
C                -1 SI QTMIN N'A PAS PU ETRE ATTEINT
C                -2 SI UN PROBLEME DE TAILLE MEMOIRE
C     REMARQUE :
C     ---------
C     POLYGONE DE 5 COTE  => 5 TRIANGULATIONS POSSIBLES
C     POLYGONE DE 10 COTE => 1430 TRIANGULATIONS POSSIBLES
C     POLYGONE DE 14 COTE => 208012 TRIANGULATIONS POSSIBLES
C     POLYGONE DE 20 COTE => 477638700 TRIANGULATIONS POSSIBLES
C
C     DANS LE CAS GENERAL LE RESULTAT EST TOUT A FAIT ACCEPTABLE :
C     POUR UN POLYGONE 14 COTES, CORRESPONDANT A 14 A 28 TRIANGLES
C     INTERSECTANTS
C     SUIVANT LE CAS : 0.08 SECONDES SUR HP700
C     AVEC 1939 POLYGONES TESTES, 54 POLYGONES EMPILES
C     POUR 208012 TRIANGULATIONS POSSIBLES.
C
C     *****************************************************************
      IMPLICIT INTEGER(I-N)
      REAL*8    X(*),QTMIN
      INTEGER IPOLYG(*),NCP,ITRPOL(*),iarr
C      REAL*8    FCRMIN
C      EXTERNAL FCRMIN
      REAL*8 TRRILF
      EXTERNAL TRRILF
      INTEGER   NPMAX,NCMAX
C
C     ---- LES TABLEAUX DE TRAVAIL ----
C
      INTEGER   IND(*),ICARD(*),JT(*),IFD(*)
      INTEGER   IFG(*),IT(*),ITM(*),IPERE(*)
      INTEGER   ITPOLY(*),ITRIA(*)
      INTEGER   ITRMIN(*)
      REAL*8      QTRIA(*),QMIN(*)
C
C     ---- VARIABLES  LOCALES ----
C
      INTEGER IDIMC
      PARAMETER (IDIMC = 2)
      INTEGER NBP,II,I,J,K,N,ITP,IPI,LIBTL,LIBTR,LIBPL
      INTEGER LPLMAX,NPTEST
C
      IF( NCP .GT. NCMAX )GO TO 888
C
      NPTEST = 0
      LPLMAX = 0
      iarr = 0
      NBP = 1
      LIBPL = 2
      IND(1)   = 1
      ICARD(1)  = NCP
      JT(1)    = 1
      IFD(1)    = -1
      IFG(1)    = -1
      QTRIA(1) = 0.0D0
      QMIN(1)  = -1.0D0
      IT(1)  = 1
      ITM(1) = 1
      IPERE(1) = 0
      DO 5 I=1,NCP
        ITPOLY(I) = IPOLYG(I)
    5 CONTINUE
      ITPOLY(NCP+1) = IPOLYG(1)
      LIBTL = NCP+1
      LIBTR = (NCP-2)*3
C
C     --- 1. BOUCLE SUR LES POLYGONES ---
C     -----------------------------------
   10 IPI  = NBP
      ITP = IND(IPI)
      N   = ICARD(IPI)
      J   = JT(IPI)
      I   = IT(IPI)
C     ------- POUR LE DEBUG ----------
      LPLMAX = MAX(LIBPL,LPLMAX)
      NPTEST = NPTEST + 1
      IF((IFD(IPI).NE.-1).AND.(IFG(IPI).NE.-1))THEN
C
C     --- ON A 1 TRIANGULATION DU POLY IPI ---
C
        IF(IFD(IPI).NE.0)
     >    QTRIA(IPI)=MIN(QMIN(IFD(IPI)),QTRIA(IPI))
        IF(IFG(IPI).NE.0)
     >    QTRIA(IPI)=MIN(QMIN(IFG(IPI)),QTRIA(IPI))
C
C     --- ON A TROUVER UNE MEILLEURE TRIANGULATION ---
C
        IF(QTRIA(IPI) .GT. QMIN(IPI) )THEN
          K = ITM(IPI)-1
          DO 15 II=1,3
            ITRMIN(K+II) = ITRIA((I-1)*3+II)
   15     CONTINUE
          IF(IFG(IPI).NE.0)THEN
          K = ITM(IPI) + 2
          DO 16 II=1,(ICARD(IFG(IPI))-2)*3
           ITRMIN(K+II)=ITRMIN(ITM(IFG(IPI))-1+II)
   16     CONTINUE
          ENDIF
C
          IF(IFD(IPI).NE.0)THEN
          IF(IFG(IPI).NE.0)THEN
            K = ITM(IPI) + 2 + (ICARD(IFG(IPI))-2)*3
          ELSE
            K = ITM(IPI) + 2
          ENDIF
          DO 17 II=1,(ICARD(IFD(IPI))-2)*3
           ITRMIN(K+II)=ITRMIN(ITM(IFD(IPI))-1+II)
   17     CONTINUE
          ENDIF
          QMIN(IPI) = QTRIA(IPI)
C     --- POUR LE DEBUG ----
C      PRINT *,'---------------------------------------'
C      PRINT *,' TRIANGULATION RETENUE POUR LE POLYGON ',IPI
C      PRINT *,' PERE DE ',IFG(IPI),' ET ',IFD(IPI)
C      PRINT *,' ',(ITRMIN(ITM(IPI)-1+II),II=1,((N-2)*3))
C      PRINT *,' QMIN ',QMIN(IPI)
C      PRINT *,'---------------------------------------'
        ENDIF
      ENDIF
C
C     --- 2. BOUCLE SUR LES TRIANGULATIONS D'UN POLYGONE
C     --------------------------------------------------
   20 J=J+1
      IF( J.GE.N )THEN
C
C     --- ON A TOUTES LES TRIANGULATIONS DU POLY IPI ---
C
        IF( IPERE(NBP).EQ.0 )THEN
C         --- ON EST A LA RACINE : ON A FINI ---
          DO 25 II=1,(ICARD(1)-2)*3
            ITRPOL(II) = ITRMIN(ITM(1)-1+II)
   25     CONTINUE
          QTMIN = QMIN(1)
C         --- POUR LE DEBUG -------------------------------
C          PRINT *,' ',(ITRPOL(II),II=1,((N-2)*3))
C          PRINT *,' NOMBRE DE POLYS TESTES : ',NPTEST
C          PRINT *,' NOMBRE MAX DE POLYS EMPILES : ',LPLMAX
          GO TO 999
        ENDIF
C       --- ON REMONTE SUR LE FILS GAUCHE ---
        IF((IFD(IPERE(IPI)).EQ.IPI).AND.
     >     (IFG(IPERE(IPI)).NE.0))THEN
           NBP  = IFG(IPERE(IPI))
           IF( LIBPL.GT.(IPI+1))THEN
C            --- ON LIBERE LES FILS DU FILS DROIT ---
             LIBPL = IPI+1
             LIBTR = ITM(IPI+1)
C             PRINT *,'LIBPL = ',LIBPL
C             PRINT *,'LIBTR = ',LIBTR
C             PRINT *,'LIBTL = ',LIBTL
           ENDIF
        ELSE
           NBP = IPERE(IPI)
           IF( LIBPL.GT.(IPI+2))THEN
C            --- ON LIBERE LES FILS DU FILS GAUCHE---
             LIBPL = IPI+2
             LIBTR = ITM(IPI+2)
             LIBTL = IND(IPI+2)
C             PRINT *,'LIBPL = ',LIBPL
C             PRINT *,'LIBTR = ',LIBTR
C             PRINT *,'LIBTL = ',LIBTL
           ENDIF
        ENDIF
        GO TO 10
      ENDIF
C     --- ON CALCULE TOUTES LES TRIANGULATIONS CONTENANT
C     --- LE TRIANGLE ITP,ITP+1,ITP+J
C     --------------------------------------------------
      QTRIA(IPI) = TRRILF(X((ITPOLY(ITP  )-1)*IDIMC+1),
     >                   X((ITPOLY(ITP+1)-1)*IDIMC+1),
     >                   X((ITPOLY(ITP+J)-1)*IDIMC+1))
C     --------- POUR LE DEBUG -------------
C      PRINT *,'POLYGON ',IPI,' ',IND(IPI),' ',ICARD(IPI),
C     >  ' TRIANGLE ',ITPOLY(ITP),' ',(ITPOLY(ITP+1)),
C     >  ' ',(ITPOLY(ITP+J)),' RIL = ',QTRIA(IPI)
      IF( QTRIA(IPI) .LE. QMIN(IPI) ) GO TO 20
      ITRIA((I-1)*3+1) = ITPOLY(ITP)
      ITRIA((I-1)*3+2) = ITPOLY(ITP+1)
      ITRIA((I-1)*3+3) = ITPOLY(ITP+J)
C
C     --- FILS GAUCHE ---
C
      IF( (N-J+1) .GT. 2 )THEN
C     --- ON DESCEND SUR LE FILS GAUCHE ---
        IF( LIBPL .EQ. NPMAX )GO TO 888
        NBP = LIBPL
        IND(NBP)  = LIBTL
        ICARD(NBP) = N-J+1
        IFG(NBP)   = -1
        IFD(NBP)   = -1
        IPERE(NBP) = IPI
        JT(NBP)   = 1
        IT(NBP)   = I+1+J-2
        ITM(NBP)  = LIBTR
        QMIN(NBP) = QMIN(IPI)
        QTRIA(NBP)= 0.0D0
        IFG(IPI) = NBP
C       ---
        IF( ((NPMAX*NCMAX)-LIBTL) .LT. (N-J+1) )GO TO 888
        DO 30 II=0,(N-(J+1))
         ITPOLY(LIBTL+II) = ITPOLY(ITP+J+II)
   30   CONTINUE
        ITPOLY(LIBTL+N-J) = ITPOLY(ITP)
        LIBTL =  LIBTL + ICARD(NBP)
        IF(((NPMAX*NCMAX)-LIBTR).LT.((ICARD(NBP)-2)*3))GO TO 888
        LIBTR =  LIBTR + ((ICARD(NBP)-2)*3)
        LIBPL = LIBPL + 1
      ELSE
        IFG(IPI) = 0
      ENDIF
C
C     --- FILS DROIT ---
C
      IF( J .GE. 3 )THEN
C     --- ON DESCEND SUR LE FILS DROIT ---
        IF( LIBPL .EQ. NPMAX )GO TO 888
        NBP = LIBPL
        IND(NBP) = ITP+1
        ICARD(NBP) = J
        IFG(NBP)   = -1
        IFD(NBP)   = -1
        IPERE(NBP) = IPI
        JT(NBP)   = 1
        IT(NBP)   = I+1
        ITM(NBP)  = LIBTR
        QMIN(NBP) = QMIN(IPI)
        QTRIA(NBP)= 0.0
        IFD(IPI)  = NBP
        IF(((NPMAX*NCMAX)-LIBTR).LT.((ICARD(NBP)-2)*3))GO TO 888
        LIBTR =  LIBTR + ((ICARD(NBP)-2)*3)
        LIBPL = LIBPL + 1
      ELSE
        IFD(IPI) = 0
      ENDIF
C
      JT(IPI) = J
C
C     --- ON TRAITE DANS L'ORDRE : IFD, IFG, IPERE
C
      GOTO 10
C
  888 iarr = -2
C
  999 END