C ZONAG     SOURCE    OF166741  25/07/28    21:15:15     12336          

C====================================================================
C   localisation de points dans des elements  d une zone elementaire
C     zonage  et construction  d un tableau de correspondance
C
C Entrees :
C    ip1 pointeur sur le maillage elementaire  massif
C    ip2 pointeur sur le maillage  poi1
C    iexx pointeur sur un segment contenant les points deja vus
C         (eviter comptage des points aux frontieres de 2 sous zones
C          et le traitement des points qui sont des noeuds de ip1)
C    ITR : COMPTEUR MIS A JOUR DANS ZONAG
C Sorties :
C    iccoun  pointeur sur le segment ICOUNT contenant le tableau de
C    correspondance suivant
C        IEINT(1,N) numero du point du meleme IP2
C        IEINT(2,N) numero de l element du meleme IP1 le contenant
C        IEINT(3,N) numero de la souszone ISZ correspondant a IP1
C        MAXPZ nombre de points de IP2 contenus dans la s_zone (N max)
C====================================================================
      SUBROUTINE ZONAG(ISZ,IP1,IP2,ICCOUN,IEXX,ITR)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC CCGEOME

-INC SMCOORD
-INC SMELEME

      DATA EPSI/1.D-5/

      SEGMENT iexclu
        INTEGER IDEJVU(nbpti)
      ENDSEGMENT

      SEGMENT icount
        INTEGER NBPTOT
        INTEGER MAXPZ(nbsz)
        INTEGER IEINT(3,NODES)
        REAL*8 QQQ(3,NODES)
      ENDSEGMENT

      SEGMENT inomil(0)

      segment wrk4
        real*8 qsi(3),xel(3,nbnn),shpp(6,nbnn)
      endsegment
      segment wrkl
        real*8 shpl(6,nbnnl),xell(3,nbnnl)
      endsegment

      SEGMENT ISEG1
        REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL)
      ENDSEGMENT
      SEGMENT ISEG3
        INTEGER NIZO(NZO+1), NUMZO(NZO), IDEJ(NZO)
      ENDSEGMENT
      SEGMENT ISEG5
        INTEGER NNMEL(ILON)
      ENDSEGMENT

      DIMENSION XPU(3)

      IDIMP1 = IDIM + 1
c*d       write(ioimp,*) 'ZONAGE',isz,ip1,ip2

      IPT2 = IP2
C*      SEGACT IPT2 <- Actif en E/S et non modifie
      NODES = IPT2.NUM(/2)

      IPT1 = IP1
C*      SEGACT IPT1 <- Actif en E/S et non modifie
      NBNN1 = IPT1.NUM(/1)
      NBEL1 = IPT1.NUM(/2)
      KELE1 = IPT1.ITYPEL

c*d      write(ioimp,*) 'ZONAGE',nbnn1,nbel1,kele1, nodes

      icount = ICCOUN
C*      SEGACT,icount*MOD <- Actif en E/S en etat modifiable

      iexclu = IEXX
C*      SEGACT,iexclu*MOD <- Actif en E/S en etat modifiable

      NBEL = NBEL1
      SEGINI ISEG1

      XZO = 0.D0
      YZO = 0.D0
      ZZO = 0.D0
      XZA = XGRAND
      YZA = XGRAND
      ZZA = XGRAND
      XTOMI = +XGRAND
      XTOMA = -XGRAND
      YTOMI = +XGRAND
      YTOMA = -XGRAND
      ZTOMI = +XGRAND
      ZTOMA = -XGRAND
      D_X = 0.D0
      D_Y = 0.D0
      D_Z = 0.D0
      DO I1 = 1, NBEL1
        XMI = +XGRAND
        XMA = -XGRAND
        YMI = +XGRAND
        YMA = -XGRAND
        ZMI = +XGRAND
        ZMA = -XGRAND
        DO I2 = 1, NBNN1
          IB = IPT1.NUM(I2,I1)
c*          iexclu.IDEJVU(IB) = 2 <- Deja fait dans ACCRO.eso
          IA = (IB-1)*IDIMP1
          XMI = MIN(XMI,XCOOR(IA+1))
          XMA = MAX(XMA,XCOOR(IA+1))
          YMI = MIN(YMI,XCOOR(IA+2))
          YMA = MAX(YMA,XCOOR(IA+2))
          IF (IDIM.EQ.3) THEN
            ZMI = MIN(ZMI,XCOOR(IA+3))
            ZMA = MAX(ZMA,XCOOR(IA+3))
          ENDIF
        ENDDO
        XLIM(1,I1) = XMI
        XLIM(2,I1) = XMA
        YLIM(1,I1) = YMI
        YLIM(2,I1) = YMA
        D_X = XMA - XMI
        XZO = MAX(XZO,D_X)
        XZA = MIN(XZA,D_X)
        XTOMI = MIN(XMI,XTOMI)
        XTOMA = MAX(XMA,XTOMA)
        D_Y = YMA - YMI
        YZO = MAX(YZO,D_Y)
        YZA = MIN(YZA,D_Y)
        YTOMI = MIN(YMI,YTOMI)
        YTOMA = MAX(YMA,YTOMA)
        IF (IDIM.EQ.3) THEN
          ZLIM(1,I1)=ZMI
          ZLIM(2,I1)=ZMA
          D_Z = ZMA - ZMI
          ZZO = MAX(ZZO,D_Z)
          ZZA = MIN(ZZA,D_Z)
          ZTOMI = MIN(ZMI,ZTOMI)
          ZTOMA = MAX(ZMA,ZTOMA)
        ENDIF
      ENDDO
      D_X = 1.D-4 * (XTOMA - XTOMI)
      XPR = MIN(XZO*1.D-2, D_X*0.5D0)
      XTOMI = XTOMI - D_X
      XTOMA = XTOMA + D_X
      D_X = XTOMA - XTOMI
      XZA = MIN(XZA*0.97D0, D_X)
      NXZO = D_X / XZA + 1
      XZO = XZA

      D_Y = 1.D-4 * (YTOMA - YTOMI)
      YPR = MIN(YZO*1.D-2, D_Y*0.5D0)
      YTOMI = YTOMI - D_Y
      YTOMA = YTOMA + D_Y
      D_Y = YTOMA - YTOMI
      YZA = MIN(YZA*0.97D0, D_Y)
      NYZO = D_Y / YZA + 1
      YZO = YZA
c*d      WRITE(6,FMT='('' NXZO NYZO'',2I7)') NXZO,NYZO
      NZZO = 1
      IF (IDIM.EQ.3) THEN
        D_Z = (ZTOMA - ZTOMI) * 1.D-4
        ZPR = MIN(ZZO*1.D-2, 0.5D0*D_Z)
        ZTOMI = ZTOMI - D_Z
        ZTOMA = ZTOMA + D_Z
        D_Z = ZTOMA - ZTOMI
        ZZA = MIN(ZZA*0.97D0, D_Z)
        IF (ZZA.EQ.0.D0) THEN
          NZZO = 1
        ELSE
          NZZO = D_Z / ZZA + 1
        ENDIF
        ZZO = ZZA
c*d      WRITE(6,FMT='('' zz0,zzA,ztomi,ztoma'',4e12.5)')
c*d     $ xzo,xza,ztomi,ztoma
      ENDIF
c*d      WRITE(6,*)  XTOMI, XTOMA
c*d      WRITE(6,*)  YTOMI ,YTOMA
c*d      WRITE(6,*)  ZTOMI ,ZTOMA
      IF (D_X.LT.EPSI) THEN
        XTOMA = XTOMA + EPSI/2
        XTOMI = XTOMI - EPSI/2
        XZO = EPSI
        D_X = XTOMA - XTOMI
      ENDIF
      IF (D_Y.LT.EPSI) THEN
        YTOMA = YTOMA + EPSI/2
        YTOMI = YTOMI - EPSI/2
        YZO = EPSI
        D_Y = YTOMA - YTOMI
      ENDIF
      IF (IDIM.EQ.3) THEN
        IF (D_Z.LT.EPSI) THEN
          ZTOMA = ZTOMA + EPSI/2
          ZTOMI = ZTOMI - EPSI/2
          ZZO = EPSI
          D_Z = ZTOMA - ZTOMI
        ENDIF
      ENDIF
C_______________________________
      NXDEP = MIN(NXZO,10)
      NYDEP = MIN(NYZO,10)
      IF (IDIM.EQ.2) THEN
        r_z = FLOAT(NXZO)*FLOAT(NYZO)
        IF (r_z.GT.100000.) THEN
          XY=SQRT(r_z)/90
          NXZO=MAX(INT(NXZO/XY),NXDEP)
          NYZO=MAX(INT(NYZO/XY),NYDEP)
          r_z = FLOAT(NXZO)*FLOAT(NYZO)
          IF (r_z.GT.100000.) THEN
            XY=SQRT(r_z)/60
            NXZO=MAX(INT(NXZO/XY),NXDEP)
            NYZO=MAX(INT(NYZO/XY),NYDEP)
          ENDIF
          XZO= D_X /NXZO
          YZO= D_Y /NYZO
          NXZO = D_X /XZO +1
          NYZO = D_Y /YZO +1
        ENDIF
c*d        WRITE(6,FMT='(''  XZO NXZO YZO NYZO '' , E12.5,I5,E12.5,I5)')
c*d     $  XZO ,NXZO, YZO, NYZO
      ELSE
         NZDEP=MIN(NZZO,10)
         r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO)
         IF (r_z .GT. 250000.) THEN
            XYZ = r_z**0.33333333D0/25
            NXZO=MAX(INT(NXZO/XYZ),NXDEP)
            NYZO=MAX(INT(NYZO/XYZ),NYDEP)
            NZZO=MAX(INT(NZZO/XYZ),NZDEP)
            r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO)
            IF (r_z.GT.250000.) THEN
               XYZ = r_z**0.3333333333D0/25
               NXZO=MAX(INT(NXZO/XYZ),NXDEP)
               NYZO=MAX(INT(NYZO/XYZ),NYDEP)
               NZZO=MAX(INT(NZZO/XYZ),NZDEP)
c*d               IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO  ''
c*d     $              ,4I7) ') NXZO,NYZO,NZZO
               r_z = FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO)
               IF (r_z.GT.250000.) THEN
                 XYZ = r_z**0.3333333333D0/25
                  NXZO=MAX(INT(NXZO/XYZ),NXDEP)
                  NYZO=MAX(INT(NYZO/XYZ),NYDEP)
                  NZZO=MAX(INT(NZZO/XYZ),NZDEP)
               ENDIF
            ENDIF
            XZO=D_X / NXZO
            YZO=D_Y / NYZO
            ZZO=D_Z / NZZO
            NXZO= D_X / XZO +1
            NYZO= D_Y / YZO +1
            NZZO= D_Z / ZZO +1
        ENDIF
      ENDIF
*
*  ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE
*  POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENTS TOUCHENT
*  CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES
*  PAR CHAQUE ELEMENT ET LEUR NOMBRE
      NXYZO = NXZO*NYZO
      NZO   = NXYZO*NZZO
c*d      WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO  ''
c*d     $,4I7) ') NZO,NXZO,NYZO,NZZO

      SEGINI ISEG3

      DO I1=1,NBEL1
        NIZ1X = INT( (XLIM(1,I1)-XTOMI-XPR) / XZO ) + 1
        NIZ1Y = INT( (YLIM(1,I1)-YTOMI-YPR) / YZO ) + 1
        NIZ2X = INT( (XLIM(2,I1)-XTOMI+XPR) / XZO ) + 1
        NIZ2Y = INT( (YLIM(2,I1)-YTOMI+YPR) / YZO ) + 1
        IF (IDIM.EQ.3) THEN
          NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1
          NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1
          DO L3=NIZ1Z,NIZ2Z
            DO L1=NIZ1Y,NIZ2Y
              DO L2=NIZ1X,NIZ2X
                NIZA = L2 + (L1-1) * NXZO + (L3-1)*NXYZO
                NUMZO(NIZA) = NUMZO(NIZA) +1
              ENDDO
            ENDDO
          ENDDO
        ELSE
          DO L1=NIZ1Y,NIZ2Y
             DO L2=NIZ1X,NIZ2X
                NIZA = L2 + ( L1-1) * NXZO
                NUMZO(NIZA) = NUMZO(NIZA) +1
             ENDDO
          ENDDO
        ENDIF
      ENDDO
*
*  CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
*  ELEMENTS CONCERNES PAR UNE ZONE
*
      ILON=0
      NIZO(1)=1
      DO L1=1,NZO
        NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
        ILON=ILON+ NUMZO(L1)
      ENDDO
C      WRITE(6,FMT='(''  ILON '',I5)') ILON
C      WRITE(6,109) (KKK,NUMZO(KKK),(NELZO(KI,KKK),KI=1,4),KKK=1,NBEL)
C  109 FORMAT(I6,I5,4I5)
C      WRITE(6,110)( NIZO(KI),KI=1,NZO+1)
C 110  FORMAT(16I5)

      SEGINI ISEG5

      DO I1=1,NBEL
        NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
        NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
        NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
        NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
        IF(IDIM.EQ.3) THEN
          NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1
          NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1
          DO 251 L3=NIZ1Z,NIZ2Z
            DO 2511 L1=NIZ1Y,NIZ2Y
              DO 25111 L2=NIZ1X,NIZ2X
                NIZA = L2 + (L1-1) * NXZO + (L3-1)*NXYZO
                IAD=NIZO(NIZA)+IDEJ(NIZA)
                NNMEL(IAD)=I1
                IDEJ(NIZA)=IDEJ(NIZA)+1
25111         CONTINUE
 2511       CONTINUE
 251      CONTINUE
        ELSE
          DO 252 L1=NIZ1Y,NIZ2Y
            DO 2521 L2=NIZ1X,NIZ2X
              NIZA = L2 + ( L1-1) * NXZO
              IAD=NIZO(NIZA)+IDEJ(NIZA)
             NNMEL(IAD)=I1
              IDEJ(NIZA)=IDEJ(NIZA)+1
 2521       CONTINUE
 252      CONTINUE
        ENDIF
      ENDDO

      nbnn = NBNN1

      segini wrk4
      segini inomil
C element quaf ? si oui kell=numero lineaire correspondant
C                si non kell=0
      CALL NQF2NL(KELE1,kell)
      IF (IERR.NE.0) RETURN
      IF (KELL.NE.0) THEN
        nsom  = nbsom(KELE1)
        nbnnl = NBNNE(kell)
        IF (nsom.NE.nbnnl) THEN
          call erreur(5)
          return
        endif
        idx = nspos(KELE1)-1
        segini wrkl
        do i=1,nbnnl
          inomil(**)=i
        enddo
      ELSE
        kd = kdegre(KELE1)
C element quadratiques
        if (kd.eq.3) then
          nso = nbsom(KELE1)
          idx = nspos(KELE1)-1
          do 762  i=1,nbnn
            do 763  j=1,nso
              iso = ibsom(idx+j)
              if (i.eq.iso) goto 762
 763        continue
            inomil(**)= i
 762      continue
C  elements  lineaires
        else if (kd.eq.2) then
          do i=1,nbnn
            inomil(**)=i
          enddo
        endif
      ENDIF
c*d       write(6,*) (inomil(i),i=1,inomil(/1))

C*      IA1 = 0
C*      IF (KELE1.EQ.14 .OR. KELE1.EQ.15) IA1 = 1
C*      IF (KELE1.EQ.16 .OR. KELE1.EQ.17) IA1 = 7

C   on va tourner sur les elements poi1

      DO 2 K = 1, NODES
        IP = IPT2.NUM(1,K)
C      write(ioimp,*) 'on cherche l element contenant le noeud' ,ip
        IF (iexclu.IDEJVU(IP).NE.0) GOTO 2

        IREFP = (IP-1) * IDIMP1
        XPU(1) = XCOOR(IREFP+1)
        IF (XPU(1).LT.XTOMI.OR.XPU(1).GT.XTOMA) GOTO 2
        XPU(2) = XCOOR(IREFP+2)
        IF (XPU(2).LT.YTOMI.OR.XPU(2).GT.YTOMA) GOTO 2
        XPU(3) = 0.D0
        IF (IDIM.EQ.3) THEN
          XPU(3) = XCOOR(IREFP+3)
          IF (XPU(3).LT.ZTOMI.OR.XPU(3).GT.ZTOMA) GOTO 2
        ENDIF
        INDZO=INT((XPU(1)-XTOMI)/XZO)+ 1 +INT((XPU(2)-YTOMI)/YZO)*NXZO
        IF (IDIM.EQ.3) INDZO=INDZO+INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO
        IDEB = NIZO(INDZO)
        IFIN = NIZO(INDZO+1)-1
C        write(6,fmt='('' ideb ifin'',2i5)') ideb,ifin ,indzo
        IF (IDEB.GT.IFIN) GO TO 2
C
        ITROU = 1
C-------------------------------------------------
        DO 11 KK = IDEB,IFIN
          jkk = NNMEL(KK)

*       write(6,*) '        '
*       write(6,*) 'recherche point ',ip, 'ds elem ',jkk,'itrou',ITROU

C Dans le cas "ACCRO", on cherche tout d'abord si le noeud 
C appartient aux deux maillages car si oui on ne le stocke pas dans 
C les noeuds a traiter ensuite (pas de rigidite a construire)
C Deja fait auparavant !
c*          DO jl = 1, NBNN1
c*            IF (IP.EQ.IPT1.NUM(jl,jkk)) THEN
c*              ITROU = 2
c*              GOTO 19
c*            ENDIF
c*          ENDDO

          CALL DOXE(xcoor,IDIM,nbnn,ipt1.num,jkk,xel)
* cas quaf
          IF (kell.NE.0) THEN
            DO innl = 1, nbnnl
              jl = ibsom(idx+innl)
              DO iid = 1, IDIM
                xell(iid,innl) = xel(iid,jl)
              ENDDO
            ENDDO
            qsi(1) = 0.D0
            qsi(2) = 0.D0
            qsi(3) = 0.D0
            CALL qsijs(xell,kell ,nbnnl,IDIM,xpu,SHPL,qsi,iret)
            CALL qsijs(xel ,KELE1,nbnn ,IDIM,xpu,SHPP,qsi,iret)
          ELSE
            qsi(1) = 0.D0
            qsi(2) = 0.D0
            qsi(3) = 0.D0
            CALL qsijs(xel,KELE1,nbnn,IDIM,xpu,SHPP,qsi,iret)
          ENDIF
          IF (IRET.NE.0) GOTO 11

c*d        write(6,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn)
          if (kell.ne.0) then
            do i=1,inomil(/1)
              ilp= inomil(i)
              if (shpl(1,ilp).lt.0.D0) goto 11
            enddo
          else
            do i=1,inomil(/1)
              ilp= inomil(i)
              if (shpp(1,ilp).lt.0.D0) goto 11
            enddo
          endif
          GOTO 100
 11     CONTINUE
C  on n a rien trouve pour ce point
        ITROU = 0
 100    CONTINUE

 19     CONTINUE
        IF (ITROU.EQ.1) THEN
C    on a trouve  l element contenant le point K  attention
          iexclu.IDEJVU(ip) = 1
          itr = itr + 1
          icount.IEINT(1,itr) = IP
          icount.IEINT(2,itr) = JKK
          icount.IEINT(3,itr) = ISZ
          icount.QQQ(1,itr) = qsi(1)
          icount.QQQ(2,itr) = qsi(2)
          icount.QQQ(3,itr) = qsi(3)
*       write(6,*) 'trouve  point' ,IP, 'Element' ,JKK
C      write(6,2375)ip,j,(xpu(i1),i1=1,IDIM),(qsi(i2),i2=1,IDIM)
2375    format(2i4,6e12.5)
C*        ELSE IF (ITROU.EQ.2) THEN
C*          iexclu.IDEJVU(ip) = 2
C*C     write(6,*) 'trouve  point' ,IP, 'noeud Element' ,jll,JKK
        ENDIF

 2    CONTINUE

      icount.MAXPZ(ISZ) = ITR  - ICOUNT.NBPTOT
      icount.NBPTOT =  ITR

      segsup wrk4,inomil
      IF (KELL.NE.0) SEGSUP wrkl
      SEGSUP,ISEG1,ISEG5,ISEG3
C
C    le tableau des correspondances est termine
C
C*      SEGDES,IPT1,IPT2 <- Actifs en E/S et non modifies
C*      SEGDES,icount,iexclu <- Actifs en E/S et modifies

      RETURN
      END