i25buc_vox1.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    i25buc_vox1            ../starter/source/interfaces/inter3d1/i25buc_vox1.F
!||--- called by ------------------------------------------------------
!||    inint3                 ../starter/source/interfaces/inter3d1/inint3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||    i25trivox1             ../starter/source/interfaces/inter3d1/i25trivox1.F
!||    insol25                ../starter/source/interfaces/inter3d1/i25buc_vox1.F
!||    margin_reduction       ../starter/source/interfaces/inter3d1/margin.F90
!||--- uses       -----------------------------------------------------
!||    margin_reduction_mod   ../starter/source/interfaces/inter3d1/margin.F90
!||    message_mod            ../starter/share/message_module/message_mod.F
!||    tri7box                ../starter/share/modules1/tri7box.F
!||====================================================================
      SUBROUTINE i25buc_vox1(
     1   X     ,IRECT,NSV    ,BUMULT,
     2   NMN   ,NRTM ,NSN   ,INTBUF_TAB  ,
     3   GAP  ,I_STOK  ,
     4   DIST  ,TZINF,MAXBOX ,MINBOX,MSR     ,   
     5   STF   ,STFN ,IDDLEVEL,
     6   GAP_S,GAP_M  ,IGAP  ,GAPMIN  ,
     7   GAPMAX,INACTI,GAP_S_L,GAP_M_L,
     8   MARGE ,ID   ,TITR   ,NBINFLG,MBINFLG,
     9   ILEV  ,MSEGTYP,GAP_N,BGAPSMX,
     A   IPARTS,KNOD2ELS,NOD2ELS,
     B   KREMNODE,REMNODE,
     C   IXS, IXS10, IXS16, IXS20,ICODE,ISKEW ,
     D   DRAD ,DGAPLOAD,NRTMT,FLAG_REMOVED_NODE,
     E   IELEM_M,nin,npari,ipari )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE tri7box
      USE names_and_titles_mod , ONLY : nchartitle
      USE intbufdef_mod 
      use margin_reduction_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "com04_c.inc"
#include      "scr06_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NMN, NRTM, NSN, I_STOK,IGAP,
     .        INACTI,
     .        IPARTS(*), KNOD2ELS(*), NOD2ELS(*)
      INTEGER IRECT(4,*),NSV(*),ICODE(*),ISKEW(*)
      INTEGER MSR(*),IDDLEVEL
      INTEGER NBINFLG(*),MBINFLG(*),ILEV,MSEGTYP(*)
      INTEGER KREMNODE(*),REMNODE(*)
      INTEGER IXS(*), IXS10(*), IXS16(*), IXS20(*)
      LOGICAL, INTENT(in) :: FLAG_REMOVED_NODE !< flag to remove some S node from the list of candidates
      INTEGER, INTENT(IN) :: IELEM_M(2,NRTM)
      integer, intent(in) :: nin !< interface index
      integer, intent(in) :: npari !< first dim of ipari
      integer, dimension(npari), intent(inout) :: ipari !< interface data
C     REAL
      my_real
     .   stf(*),stfn(*),x(3,*),gap_s(*),gap_m(*),
     .   dist,bumult,gap,tzinf,maxbox,minbox,gapmin,gapmax,
     .   gap_s_l(*),gap_m_l(*),marge,gap_n(4,*)
      my_real
     .   bgapsmx
      my_real , INTENT(IN) :: drad, dgapload
      INTEGER ID
      CHARACTER(LEN=NCHARTITLE) :: TITR
 
      INTEGER , INTENT(IN) :: NRTMT
      type(intbuf_struct_), intent(inout) :: intbuf_tab !< interface data 
      integer :: candidate_count
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J, L, N1, N2, N3, N4,N_SOL, ESHIFT
      INTEGER IBID, NELS
      INTEGER,  DIMENSION(:),ALLOCATABLE :: IIX,IIY,IIZ
C     REAL
      my_real
     .   dx1,dy1,dz1,
     .   dx3,dy3,dz3,
     .   dx4,dy4,dz4,
     .   dx6,dy6,dz6,
     .   dd1,dd2,dd3,dd4,dd,dd0,xmin,ymin,zmin,
     .   xmax,ymax,zmax,tzinf0,minbox0,maxbox0,
     .   bid,tzinf_st,marge_st,dd_st,d_max,pensol,d_moy,
     .   xyzm(6),bminma(6),aaa,ledgmax
      my_real,
     .   DIMENSION(:), ALLOCATABLE :: edge_l2,edge_l2_tmp
      INTEGER TAGP(NPART),IAD,N,IE,IL,IP
      INTEGER, DIMENSION(:), ALLOCATABLE :: NPARTNS,IELEM,LPARTNS
      INTEGER NBX,NBY,NBZ
      INTEGER, DIMENSION(:), ALLOCATABLE :: IS_LARGE_NODE,LARGE_NODE, TAGNOD,LOCAL_NEXT_NOD
      INTEGER :: NB_LARGE_NODES
      INTEGER (KIND=8) :: NBX8,NBY8,NBZ8,RES8,LVOXEL8
      INTEGER (KIND=8) :: IONE,IHUNDRED !< Integer constants in INTEGER 8 for comparisions
      my_real :: xmin_base, ymin_base, zmin_base, xmax_base, ymax_base, zmax_base
C-----------------------------------------------
      IONE=1
      ihundred=100
      gapmax=ep30
      gapmin=zero
 
C-------use temporarily GAP_N(1,*)=V/A
C       Not used :: MVOISN(1,*)-> MTYPE(solid),MVOISN(2,*)-> E_id
C
C
C     1-CALCULATION TAILLE DES ZONES INFLUENCES
C     dd is the average element length
C     dd_st is the max element length
      dd=zero
      dd_st=zero
      pensol=ep30
      d_moy = zero
      n_sol = 0
      DO 10 l=1,nrtm
C      ELEMENT CONNECTIVITIES
       n1=irect(1,l)
       n2=irect(2,l)
       n3=irect(3,l)
       n4=irect(4,l)
C      LONGUEUR COTE 1
       dx1=(x(1,n1)-x(1,n2))
       dy1=(x(2,n1)-x(2,n2))
       dz1=(x(3,n1)-x(3,n2))
       dd1=sqrt(dx1**2+dy1**2+dz1**2)
C      LONGUEUR COTE 2
       dx3=(x(1,n1)-x(1,n4))
       dy3=(x(2,n1)-x(2,n4))
       dz3=(x(3,n1)-x(3,n4))
       dd2=sqrt(dx3**2+dy3**2+dz3**2)
C      LONGUEUR COTE 3
       dx4=(x(1,n3)-x(1,n2))
       dy4=(x(2,n3)-x(2,n2))
       dz4=(x(3,n3)-x(3,n2))
       dd3=sqrt(dx4**2+dy4**2+dz4**2)
C      LONGUEUR COTE 4
       dx6=(x(1,n4)-x(1,n3))
       dy6=(x(2,n4)-x(2,n3))
       dz6=(x(3,n4)-x(3,n3))
       dd4=sqrt(dx6**2+dy6**2+dz6**2)
       dd=dd+ (dd1+dd2+dd3+dd4)
C-------only for solid---  and coating shell
       IF (msegtyp(l)==0.OR.msegtyp(l)>nrtmt) THEN
        d_max=max(dd1,dd2,dd3,dd4)
        d_max=min(d_max,gap_n(1,l))
C--------correction of too huge GAP_N(1,L) w/ irregular mesh         
        gap_n(1,l)=d_max
        dd_st=max(dd_st,d_max)
        n_sol = n_sol + 1
        d_moy = d_moy + d_max
       END IF
C       IF (MSEGTYP(L)==0) DD_ST=MAX(DD_ST,DD1,DD2,DD3,DD4)
  10  CONTINUE
C     TAILLE ZINF = .1*TAILLE MOYENNE ELEMENT DE CHAQUE COTE
C     TAILLE BUCKET MIN = TZINF  * BUMULT
      dd0=dd/nrtm/four
      IF(nrtm > 0 .AND. nrtm <= 3) THEN
        call margin_reduction(x,numnod,irect,nrtm,nsv,nsn,drad,gap,dgapload,bumult,stfn,dd0)
! special case for floor and wall with only one element
! the margin may be too large when a wall is modeled with only one element
      END IF
 
      dd=dd0
 
C     TZINF = BUMULT*DD
      marge = bumult*dd
 
C calculation of tzinf based on the margin and not the other way around
      tzinf = marge + max(gap+dgapload,drad)
C marge_st: margin independent of bumult as input to find the same initial penalties (delete or coordinate change cases)
      marge_st = bmul0*dd
 
C if iddlevel = 0, then 1st pass here, with the engine's margin to find the same candidates
      IF(iddlevel==0) marge_st = marge
      tzinf_st = marge_st + max(gap+dgapload,drad)
      bid = four_over_5*dd
      IF (inacti/=7.AND.tzinf>bid) THEN
        ibid = nint(tzinf/dd0)
        ibid =(2*ibid+4)*ibid*2
      ENDIF
C      MAXBOX= ZEP9*(DD - TZINF)
C DD + 2*Tzinfo: Size element increased in influence zone
      maxbox= half*(dd + 2*tzinf)
      minbox= half*maxbox
      tzinf0 = tzinf
      minbox0 = minbox
      maxbox0 = maxbox
C     SET TO ZERO TO PERFORM COMPLETE SEARCH CYCLE 0 ENGINE
      dist = zero     
 
C--------------------------------
C     calculation of the domain boundaries
C--------------------------------
      bminma(1)=-ep30
      bminma(2)=-ep30
      bminma(3)=-ep30
      bminma(4)= ep30
      bminma(5)= ep30
      bminma(6)= ep30
C
      xmin=ep30
      xmax=-ep30
      ymin=ep30
      ymax=-ep30
      zmin=ep30
      zmax=-ep30 
C
      DO 20 i=1,nmn
        j=msr(i) 
        xmin= min(xmin,x(1,j))
        ymin= min(ymin,x(2,j))
        zmin= min(zmin,x(3,j))
        xmax= max(xmax,x(1,j))
        ymax= max(ymax,x(2,j))
        zmax= max(zmax,x(3,j))
 20   CONTINUE
C
      xmin=xmin-tzinf_st
      ymin=ymin-tzinf_st
      zmin=zmin-tzinf_st
      xmax=xmax+tzinf_st
      ymax=ymax+tzinf_st
      zmax=zmax+tzinf_st
C
      DO 25 i=1,nsn
        j=nsv(i)
        xmin= min(xmin,x(1,j))
        ymin= min(ymin,x(2,j))
        zmin= min(zmin,x(3,j))
        xmax= max(xmax,x(1,j))
        ymax= max(ymax,x(2,j))
        zmax= max(zmax,x(3,j))
 25   CONTINUE
C
      bminma(1) = max(bminma(1),xmax)
      bminma(2) = max(bminma(2),ymax)
      bminma(3) = max(bminma(3),zmax)
      bminma(4) = min(bminma(4),xmin)
      bminma(5) = min(bminma(5),ymin)
      bminma(6) = min(bminma(6),zmin)
 
      xyzm(1) = bminma(4)
      xyzm(2) = bminma(5)
      xyzm(3) = bminma(6)
      xyzm(4) = bminma(1)
      xyzm(5) = bminma(2)
      xyzm(6) = bminma(3)
 
      aaa = sqrt(nmn /
     .           ((bminma(1)-bminma(4))*(bminma(2)-bminma(5))
     .           +(bminma(2)-bminma(5))*(bminma(3)-bminma(6))
     .           +(bminma(3)-bminma(6))*(bminma(1)-bminma(4))))
 
      aaa = 0.75*aaa
 
      nbx = nint(aaa*(bminma(1)-bminma(4)))
      nby = nint(aaa*(bminma(2)-bminma(5)))
      nbz = nint(aaa*(bminma(3)-bminma(6)))
      nbx = max(nbx,1)
      nby = max(nby,1)
      nbz = max(nbz,1)
 
      nbx8=nbx
      nby8=nby
      nbz8=nbz
      res8=(nbx8+2)*(nby8+2)*(nbz8+2)
      lvoxel8 = lvoxel      
 
      IF(res8 > lvoxel8) THEN
        aaa = lvoxel
        aaa = aaa/((nbx8+2)*(nby8+2)*(nbz8+2))
        aaa = aaa**(third)
        nbx = int((nbx+2)*aaa)-2
        nby = int((nby+2)*aaa)-2
        nbz = int((nbz+2)*aaa)-2
        nbx = max(nbx,1)
        nby = max(nby,1)
        nbz = max(nbz,1)
      ENDIF
      
      nbx8=nbx
      nby8=nby
      nbz8=nbz
      res8=(nbx8+2)*(nby8+2)*(nbz8+2)
      
      IF(res8 > lvoxel8) THEN
        nbx = min(ihundred,max(nbx8,ione))
        nby = min(ihundred,max(nby8,ione))
        nbz = min(ihundred,max(nbz8,ione))
      END IF  
 
C--------------------------------
C     Initializations case of IDDLEVEL=0 or INACTI/=5 and INACTI/=-1
C--------------------------------
      ALLOCATE(npartns(nsn+1),ielem(nrtm),edge_l2(nsn))
 
      npartns(1:nsn+1) = 0
      ielem(1:nrtm)    = 0
      edge_l2(1:nsn)   = zero 
      ledgmax          = zero
C--------------------------------
 
 
 
       nb_large_nodes = 0
       ALLOCATE(large_node(nsn))
       ALLOCATE(is_large_node(nsn))
       ALLOCATE(tagnod(numnod))
       is_large_node(1:nsn) = 0
       large_node(1:nsn) = 0
       nb_large_nodes = 0
       tagnod(1:numnod) = 0
 
      IF(iddlevel==1)THEN
        ! Looking for initial penetration (vs solids) if and only if IDDLEVEL=1
C            
        DO ie=1,nrtm
           ! looks for solid element supporting the segment (cf initial penetrations vs solids)
           IF(ielem_m(1,ie)<=numels) THEN
              ielem(ie)= ielem_m(1,ie)
           ELSE
             CALL insol25(irect    ,ixs     ,ixs10,ixs16,ixs20,
     .                  knod2els ,nod2els ,nels ,ie   )
             ielem(ie)=nels
           ENDIF
        END DO
C
        IF(inacti==5.OR.inacti==-1)THEN 
 
          ALLOCATE(edge_l2_tmp(numnod))
          edge_l2_tmp(1:numnod)=zero
          n_sol = 0
          DO ie=1,nrtm
            IF(stf(ie)> zero)THEN ! Solids only
              DO il=1,4
                n1=irect(il,ie)
                n2=irect(mod(il,4)+1,ie)
 
                aaa = (x(1,n2)-x(1,n1))*(x(1,n2)-x(1,n1))
     .              + (x(2,n2)-x(2,n1))*(x(2,n2)-x(2,n1))
     .              + (x(3,n2)-x(3,n1))*(x(3,n2)-x(3,n1)) 
                edge_l2_tmp(n1) = max(edge_l2_tmp(n1), aaa )
                edge_l2_tmp(n2) = max(edge_l2_tmp(n2), aaa )
                IF (msegtyp(ie)==0.OR.msegtyp(ie)>nrtmt) THEN ! Solids only
                   IF(tagnod(n1)==0) THEN
                      n_sol= n_sol + 1
                      tagnod(n1) = 1
                   ENDIF
                   IF(tagnod(n2)==0) THEN
                      n_sol= n_sol + 1
                      tagnod(n2) = 1
                   ENDIF
                ENDIF
              END DO
            ENDIF
 
          END DO
C
          
          DO i=1,nsn
            n=nsv(i)
            IF(stfn(i)/=zero) THEN
              edge_l2(i) = half*sqrt(edge_l2_tmp(n))
              IF(tagnod(n)==1) ledgmax=ledgmax+edge_l2(i)!MAX(LEDGMAX,EDGE_L2(I))
            END IF
          END DO 
 
          IF(n_sol > 0) ledgmax=half*ledgmax/n_sol
 
 
          DEALLOCATE(edge_l2_tmp)
 
        ENDIF
 
c        IF(INACTI==5.OR.INACTI==-1)THEN ! Do it for all Inacti ( warnings initial penetrations)
          tagp(1:npart)   =0
C
          DO i=1,nsn
            n=nsv(i)
            DO iad=knod2els(n)+1,knod2els(n+1)
              ie=nod2els(iad)
              ip=iparts(ie)
              IF(tagp(ip)==0)THEN
                npartns(i)=npartns(i)+1
                tagp(ip)  =1
              END IF
            END DO
            DO iad=knod2els(n)+1,knod2els(n+1)
              ie=nod2els(iad)
              ip=iparts(ie)
              tagp(ip)  =0 ! Reset
            END DO
          END DO
C
          DO i=1,nsn
            npartns(i+1) = npartns(i+1) + npartns(i)
          END DO
          DO i=nsn,1,-1
            npartns(i+1) = npartns(i)
          END DO
          npartns(1)=0
C
          ALLOCATE(lpartns(npartns(nsn+1)))
C
          DO i=1,nsn
            n=nsv(i)
            DO iad=knod2els(n)+1,knod2els(n+1)
              ie=nod2els(iad)
              ip=iparts(ie)
              IF(tagp(ip)==0)THEN
                npartns(i)=npartns(i)+1
                lpartns(npartns(i))=ip
                tagp(ip)  =1
              END IF
            END DO
            DO iad=knod2els(n)+1,knod2els(n+1)
              ie=nod2els(iad)
              ip=iparts(ie)
              tagp(ip)  =0 ! Reset
            END DO
          END DO
C
          DO i=nsn,1,-1
            npartns(i+1) = npartns(i)
          END DO
          npartns(1)=0
C
c        ELSE
c          ALLOCATE(LPARTNS(0))
c        END IF
      ELSE
        ALLOCATE(lpartns(0))
      END IF
C
 
C     complete initialization of voxel
      DO i=inivoxel,(nbx+2)*(nby+2)*(nbz+2)
        voxel1(i)=0
      ENDDO
      inivoxel = max(inivoxel,(nbx+2)*(nby+2)*(nbz+2)+1)
C
      eshift=0
C
      i_stok = 0
 
      ALLOCATE(local_next_nod(nsn))
      ALLOCATE(iix(nsn))
      ALLOCATE(iiy(nsn))
      ALLOCATE(iiz(nsn))
C
!$OMP PARALLEL
      CALL i25trivox1(
     1   nsn     ,irect   ,x       ,
     2   stfn    ,xyzm    ,nsv     ,i_stok  ,
     3   eshift  ,bgapsmx ,
     4   voxel1  ,nbx     ,nby     ,nbz     ,nrtm    ,
     5   gap_s   ,gap_m   ,marge_st,
     6   nbinflg ,mbinflg ,ilev    ,msegtyp ,
     7   igap    ,gap_s_l ,gap_m_l ,edge_l2 ,ledgmax ,
     8   kremnode,remnode,
     9   iparts  ,npartns ,lpartns ,ielem   ,icode   ,
     a   iskew   ,drad, is_large_node, large_node, nb_large_nodes,
     b   dgapload,nrtmt,flag_removed_node,
     c   ielem_m,local_next_nod,iix,iiy,iiz,
     d   intbuf_tab,ipari,nin)
!$OMP END PARALLEL
      DEALLOCATE(local_next_nod)
      DEALLOCATE(iix)
      DEALLOCATE(iiy)
      DEALLOCATE(iiz)
 
      DEALLOCATE(edge_l2,npartns,ielem,lpartns)
      DEALLOCATE(is_large_node,large_node,tagnod)
C
      IF(nsn/=0)THEN
       WRITE(iout,*)' POSSIBLE IMPACT NUMBER, NSN:',i_stok,nsn
C
      ELSE
       CALL ancmsg(msgid=552,
     .             msgtype=msgwarning,
     .             anmode=aninfo_blind_2,
     .             i1=id,
     .             c1=titr)
      ENDIF
C
      RETURN
      END
 
!||====================================================================
!||    insol25       ../starter/source/interfaces/inter3d1/i25buc_vox1.F
!||--- called by ------------------------------------------------------
!||    i25buc_vox1   ../starter/source/interfaces/inter3d1/i25buc_vox1.F
!||--- uses       -----------------------------------------------------
!||====================================================================
      SUBROUTINE insol25(IRECT   ,IXS    ,IXS10,IXS16,IXS20,
     .                   KNOD2ELS,NOD2ELS,NEL  ,I    )
      use element_mod , only :nixs
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL, I
      INTEGER IRECT(4,*), IXS(NIXS,*), KNOD2ELS(*), NOD2ELS(*), 
     .        IXS10(6,*), IXS16(8,*), IXS20(12,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER N, JJ, II, K, IC, IAD,
     .        NUSER, NUSERM
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
C
      NEL=0
 
      if(numels==0) RETURN
 
      IF(irect(1,i)>numnod) RETURN
 
      ic=0
      nuserm = -1
      DO 230 iad=knod2els(irect(1,i))+1,knod2els(irect(1,i)+1)
       n = nod2els(iad)
       IF(n <= numels8)THEN
         DO 210 jj=1,4
           ii=irect(jj,i)
           DO k=1,8
             IF(ixs(k+1,n)==ii) GOTO 210
           ENDDO
           GOTO 230
  210    CONTINUE
         ic=ic+1
         nuser = ixs(11,n)
         IF (nuser>nuserm) THEN
           nel = n
           nuserm = nuser
         ENDIF
       ELSEIF(n <= numels8+numels10)THEN
         DO 220 jj=1,4
           ii=irect(jj,i)
           DO k=1,8
             IF(ixs(k+1,n)==ii) GOTO 220
           ENDDO
           DO k=1,6
             IF(ixs10(k,n-numels8)==ii) GOTO 220
           ENDDO
           GOTO 230
  220    CONTINUE
         ic=ic+1
         nuser = ixs(11,n)
         IF (nuser>nuserm) THEN
           nel = n
           nuserm = nuser
         ENDIF
       ELSEIF(n <= numels8+numels10+numels20)THEN
         DO 222 jj=1,4
           ii=irect(jj,i)
           DO k=1,8
             IF(ixs(k+1,n)==ii) GOTO 222
           ENDDO
           DO k=1,12
             IF(ixs20(k,n-numels8-numels10)==ii) GOTO 222
           ENDDO
           GOTO 230
  222    CONTINUE
         ic=ic+1
         nuser = ixs(11,n)
         IF (nuser>nuserm) THEN
           nel = n
           nuserm = nuser
         ENDIF
       ELSEIF(n <= numels8+numels10+numels20+numels16)THEN
         DO 224 jj=1,4
           ii=irect(jj,i)
           DO k=1,8
             IF(ixs(k+1,n)==ii) GOTO 224
           ENDDO
           DO k=1,8
             IF(ixs16(k,n-numels8-numels10-numels20)==ii) GOTO 224
           ENDDO
           GOTO 230
  224    CONTINUE
         ic=ic+1
         nuser = ixs(11,n)
         IF (nuser>nuserm) THEN
           nel = n
           nuserm = nuser
         ENDIF
       ELSE
         GOTO 230
       END IF
  230 CONTINUE
      RETURN
      END