i7buc_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/.
!||====================================================================
!||    i7buc_vox1             ../starter/source/interfaces/inter3d1/i7buc_vox1.F
!||--- called by ------------------------------------------------------
!||    inint3                 ../starter/source/interfaces/inter3d1/inint3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||    i7trivox1              ../starter/source/interfaces/inter3d1/i7trivox1.F
!||    margin_reduction       ../starter/source/interfaces/inter3d1/margin.F90
!||--- uses       -----------------------------------------------------
!||    front_mod              ../starter/share/modules1/front_mod.F
!||    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 i7buc_vox1(
     1   X     ,IRECT,NSV    ,BUMULT,NSEG    ,
     2   NMN   ,NRTM ,MWA    ,NSN   ,
     3   GAP  ,XYZM   ,NOINT ,I_STOK  ,
     4   DIST  ,TZINF,MSR     ,   
     5   STF   ,STFN ,MULTIMP,ISTF  ,IDDLEVEL,
     6   ITAB  ,GAP_S,GAP_M  ,IGAP  ,GAPMIN  ,
     7   GAPMAX,INACTI,GAP_S_L,GAP_M_L,I_MEM ,
     8   NCONT ,ICURV ,BGAPSMX,ID, TITR,
     9   DRAD  ,INTERCEP,NIN,
     1   IREMNODE,FLAGREMNODE,KREMNODE,REMNODE,
     2   DGAPLOAD,npari,ipari,intbuf_tab,IS_USED_WITH_LAW151)
C============================================================================
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE tri7box
      USE front_mod
      USE names_and_titles_mod , ONLY : nchartitle
      use intbufdef_mod 
      use margin_reduction_mod , only : margin_reduction
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "vect07_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, NOINT,I_STOK,MULTIMP,ISTF,IGAP,
     .        INACTI,I_MEM,NIN,IREMNODE,FLAGREMNODE
      INTEGER IRECT(4,*),NSV(*),NSEG(*),MWA(*)
      INTEGER MSR(*),IDDLEVEL
      INTEGER ITAB(*),NCONT,ICURV,KREMNODE(*),REMNODE(*)
      my_real
     .   STF(*),STFN(*),X(3,*),XYZM(6,2),GAP_S(*),GAP_M(*),
     .   DIST,BUMULT,GAP,TZINF,GAPMIN,GAPMAX,
     .   GAP_S_L(*),GAP_M_L(*),BGAPSMX, DRAD
      my_real , INTENT(IN) :: DGAPLOAD
      INTEGER ID
      CHARACTER(LEN=NCHARTITLE) :: TITR
      TYPE(INTERSURFP) :: INTERCEP(3,NINTER)
      integer, intent(in) :: npari
      integer, dimension(npari), intent(inout) :: ipari !< interface data
      type(intbuf_struct_), intent(inout) :: intbuf_tab !< interface data
      LOGICAL,INTENT(IN) :: IS_USED_WITH_LAW151
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NRTM_L
      INTEGER, DIMENSION(:), ALLOCATABLE :: INDEX
      INTEGER NBX,NBY,NBZ
      INTEGER (KIND=8) :: NBX8,NBY8,NBZ8,RES8,LVOXEL8
      INTEGER I, J, K, I_ADD, L, LOC_PROC, N, ISZNSNR,
     .        n1, n2, n3, n4, ncontact,i_bid,i_stok_old,
     .        ix1,iy1,iz1,ix2,iy2,iz2,ix,iy,iz
      my_real
     .         marge, aaa,tzinf_st,marge_st 
      my_real
     .   dx1,dy1,dz1,
     .   dx3,dy3,dz3,
     .   dx4,dy4,dz4,
     .   dx6,dy6,dz6,
     .   dd1,dd2,dd3,dd4,dd,dd0,xmin,ymin,zmin,
     .   xmax_m,ymax_m,zmax_m,xmin_m,ymin_m,zmin_m,
     .   xmax_s,ymax_s,zmax_s,xmin_s,ymin_s,zmin_s,
     .   xmax,ymax,zmax,xxx,yyy,zzz,
     .   xminb, yminb, zminb, xmaxb, ymaxb, zmaxb,
     .   mean_x, mean_y, mean_z, dev_x, dev_y, dev_z,
     .   gapv(mvsiz),c_max,
     .   tstart,tstop,
     .   xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,
     .   xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4
        my_real, DIMENSION(:), ALLOCATABLE :: curv_max
      INTEGER,  DIMENSION(:),ALLOCATABLE :: IIX,IIY,IIZ,LOCAL_NEXT_NOD
     
      LOGICAL :: TYPE18
      INTEGER (KIND=8) :: IONE,IHUNDRED !< Integer constants in INTEGER 8 for comparisions
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      IONE=1        !< Constant init
      ihundred=100  !< Constant init
 
      ALLOCATE( index(nrtm) )
      ALLOCATE( curv_max(nrtm) )
      ALLOCATE(local_next_nod(nsn))
      ALLOCATE(iix(nsn))
      ALLOCATE(iiy(nsn))
      ALLOCATE(iiz(nsn))
      type18=.false.
      IF(inacti==7)type18=.true.
 
      mwa(1:numnod+numfakenodigeo) = 0
      ncontact = multimp * ncont
      c_max = zero
      IF(icurv/=0)THEN
        DO i=1,nrtm
          xxx=max(x(1,irect(1,i)),x(1,irect(2,i)),
     .            x(1,irect(3,i)),x(1,irect(4,i)))
     .       -min(x(1,irect(1,i)),x(1,irect(2,i)),
     .            x(1,irect(3,i)),x(1,irect(4,i)))
          yyy=max(x(2,irect(1,i)),x(2,irect(2,i)),
     .            x(2,irect(3,i)),x(2,irect(4,i)))
     .       -min(x(2,irect(1,i)),x(2,irect(2,i)),
     .          x(2,irect(3,i)),x(2,irect(4,i)))
          zzz=max(x(3,irect(1,i)),x(3,irect(2,i)),
     .            x(3,irect(3,i)),x(3,irect(4,i)))
     .       -min(x(3,irect(1,i)),x(3,irect(2,i)),
     .            x(3,irect(3,i)),x(3,irect(4,i)))
          curv_max(i) = half * max(xxx,yyy,zzz)
          c_max = max(c_max,curv_max(i))
        ENDDO
      ELSE
        DO i=1,nrtm
          curv_max(i)=zero
        ENDDO
      ENDIF
C-----
      dd=zero
      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)
  10  CONTINUE
      dd0=dd/nrtm/four
      IF(nrtm >0 .AND. nrtm <= 3 .AND. .not. is_used_with_law151) THEN
        call margin_reduction(x,numnod,irect,nrtm,nsv,nsn,drad,gap,dgapload,bumult,stfn,dd0)
      ENDIF
      dd =dd0
C      DD = MAX(DD0,ONEP251*GAP)
C      DD = MAX(DD ,ONEP251*DRAD)
C remove this additional test which breaks igap performance
C in case this factor is needed, then needs to compute mean gap instead
C 
C Marge: Bumult Bumult Engine Marge in Input (adjustable for performance)
      marge = bumult*dd
      tzinf = marge + max(gap+dgapload,drad)
C marge_st: margin independent of bumult as input necessary to find the initial penalties in a parith/on manner
      marge_st = bmul0*dd
C 1st passage with margin x2 to find more candidates (amounts to taking bmul0 = 0.4 instead of 0.2)
C      IF(IDDLEVEL==0) MARGE_ST = 2*MARGE_ST
      IF(iddlevel==0) marge_st = marge
      tzinf_st = marge_st + max(gap+dgapload,drad)
C     SET TO ZERO TO PERFORM COMPLETE SEARCH CYCLE 0 ENGINE
      dist = zero     
C     TEST TO BY-PASS SECOND SEARCH
C     IF(IDDLEVEL==1 .AND. MARGE >= MARGE_ST) GOTO 999 !
C----- DOMAIN BOUNDARIES
      xmax_m=-ep30
      ymax_m=-ep30
      zmax_m=-ep30 
      xmin_m=ep30
      ymin_m=ep30
      zmin_m=ep30
 100  CONTINUE     
      i_stok = 0
      i_mem = 0
C Loop over the domain to improve sorting performance      
      DO loc_proc=1,nspmd
        nrtm_l=0
        DO i=1,nrtm
          IF(intercep(1,nin)%P(i)==loc_proc)THEN
            nrtm_l=nrtm_l+1
            index(nrtm_l)=i
          END IF
        END DO
C
c        print *,NOINT,'***************start loop:',
c     .  LOC_PROC,NSPMD,NRTM,NRTM_L
        IF(nrtm_l == 0)cycle
        mean_x=zero
        mean_y=zero
        mean_z=zero
        DO k=1,nrtm_l
          i = index(k)
          j=irect(1,i) 
          xmax_m= max(xmax_m,x(1,j))
          ymax_m= max(ymax_m,x(2,j))
          zmax_m= max(zmax_m,x(3,j))
          xmin_m= min(xmin_m,x(1,j))
          ymin_m= min(ymin_m,x(2,j))
          zmin_m= min(zmin_m,x(3,j))
          mean_x=mean_x+x(1,j)
          mean_y=mean_y+x(2,j)
          mean_z=mean_z+x(3,j)
          j=irect(2,i) 
          xmax_m= max(xmax_m,x(1,j))
          ymax_m= max(ymax_m,x(2,j))
          zmax_m= max(zmax_m,x(3,j))
          xmin_m= min(xmin_m,x(1,j))
          ymin_m= min(ymin_m,x(2,j))
          zmin_m= min(zmin_m,x(3,j))
          mean_x=mean_x+x(1,j)
          mean_y=mean_y+x(2,j)
          mean_z=mean_z+x(3,j)
          j=irect(3,i) 
          xmax_m= max(xmax_m,x(1,j))
          ymax_m= max(ymax_m,x(2,j))
          zmax_m= max(zmax_m,x(3,j))
          xmin_m= min(xmin_m,x(1,j))
          ymin_m= min(ymin_m,x(2,j))
          zmin_m= min(zmin_m,x(3,j))
          mean_x=mean_x+x(1,j)
          mean_y=mean_y+x(2,j)
          mean_z=mean_z+x(3,j)
          j=irect(4,i) 
          xmax_m= max(xmax_m,x(1,j))
          ymax_m= max(ymax_m,x(2,j))
          zmax_m= max(zmax_m,x(3,j))
          xmin_m= min(xmin_m,x(1,j))
          ymin_m= min(ymin_m,x(2,j))
          zmin_m= min(zmin_m,x(3,j))
          mean_x=mean_x+x(1,j)
          mean_y=mean_y+x(2,j)
          mean_z=mean_z+x(3,j)
        END DO
 
C MIN/MAX DOMAINE        
        xmin=xmin_m-tzinf_st
        ymin=ymin_m-tzinf_st
        zmin=zmin_m-tzinf_st
        xmax=xmax_m+tzinf_st
        ymax=ymax_m+tzinf_st
        zmax=zmax_m+tzinf_st
c        print*,'bornes min:',XMIN,YMIN,ZMIN
c        print*,'bornes max:',XMAX,YMAX,ZMAX
 
        mean_x=mean_x/max((4*nrtm_l),1)
        mean_y=mean_y/max((4*nrtm_l),1)
        mean_z=mean_z/max((4*nrtm_l),1)
C     STD DEVIATION & CRVOXEL
        dev_x=zero
        dev_y=zero
        dev_z=zero
        crvoxel(0:lrvoxel,0:lrvoxel)=0
        nbx = lrvoxel
        nby = lrvoxel
        nbz = lrvoxel
        DO k=1,nrtm_l
          i = index(k)
          n1 = irect(1,i)
          n2 = irect(2,i)
          n3 = irect(3,i)
          n4 = irect(4,i)
          xx1=x(1,n1)
          xx2=x(1,n2)
          xx3=x(1,n3)
          xx4=x(1,n4)
          xmaxe=max(xx1,xx2,xx3,xx4)
          xmine=min(xx1,xx2,xx3,xx4)
          dev_x=dev_x+(xx1-mean_x)**2+(xx2-mean_x)**2
     .              +(xx3-mean_x)**2+(xx4-mean_x)**2
          yy1=x(2,n1)
          yy2=x(2,n2)
          yy3=x(2,n3)
          yy4=x(2,n4)
          ymaxe=max(yy1,yy2,yy3,yy4)
          ymine=min(yy1,yy2,yy3,yy4)
          dev_y=dev_y+(yy1-mean_y)**2+(yy2-mean_y)**2
     .               +(yy3-mean_y)**2+(yy4-mean_y)**2
          zz1=x(3,n1)
          zz2=x(3,n2)
          zz3=x(3,n3)
          zz4=x(3,n4)
          zmaxe=max(zz1,zz2,zz3,zz4)
          zmine=min(zz1,zz2,zz3,zz4)
          dev_z=dev_z+(zz1-mean_z)**2+(zz2-mean_z)**2
     .               +(zz3-mean_z)**2+(zz4-mean_z)**2
 
C        index of voxels occupied by the facet
 
          ix1=int(nbx*(xmine-tzinf_st-xmin)/(xmax-xmin))
          iy1=int(nby*(ymine-tzinf_st-ymin)/(ymax-ymin))
          iz1=int(nbz*(zmine-tzinf_st-zmin)/(zmax-zmin))
          ix1=max(0,min(nbx,ix1))
          iy1=max(0,min(nby,iy1))
          iz1=max(0,min(nbz,iz1))
          ix2=int(nbx*(xmaxe+tzinf_st-xmin)/(xmax-xmin))
          iy2=int(nby*(ymaxe+tzinf_st-ymin)/(ymax-ymin))
          iz2=int(nbz*(zmaxe+tzinf_st-zmin)/(zmax-zmin))
          ix2=max(0,min(nbx,ix2))
          iy2=max(0,min(nby,iy2))
          iz2=max(0,min(nbz,iz2))
 
          DO iz = iz1, iz2
            DO iy = iy1, iy2
              DO ix = ix1, ix2
                crvoxel(iy,iz)=ibset(crvoxel(iy,iz),ix)
              END DO
            END DO
          END DO
          
        END DO
        dev_x=sqrt(dev_x/max(4*nrtm_l,1))
        dev_y=sqrt(dev_y/max(4*nrtm_l,1))
        dev_z=sqrt(dev_z/max(4*nrtm_l,1))
C we keep all in +/- 2 sigma => 95% of the population for normal law distribution 
        xminb=max(mean_x-2*dev_x,xmin)
        yminb=max(mean_y-2*dev_y,ymin)
        zminb=max(mean_z-2*dev_z,zmin)
        xmaxb=min(mean_x+2*dev_x,xmax)
        ymaxb=min(mean_y+2*dev_y,ymax)
        zmaxb=min(mean_z+2*dev_z,zmax)
C 2D special case test
        IF(abs(xminb-xmaxb) < em10)THEN
          xminb=xmin
          xmaxb=xmax
        END IF
        IF(abs(yminb-ymaxb) < em10)THEN
          yminb=ymin
          ymaxb=ymax
        END IF
        IF(abs(zminb-zmaxb) < em10)THEN
          zminb=zmin
          zmaxb=zmax
        END IF        
 
        xyzm(1,1) = xmin
        xyzm(2,1) = ymin
        xyzm(3,1) = zmin
        xyzm(4,1) = xmax
        xyzm(5,1) = ymax
        xyzm(6,1) = zmax
        xyzm(1,2) = xminb
        xyzm(2,2) = yminb
        xyzm(3,2) = zminb
        xyzm(4,2) = xmaxb
        xyzm(5,2) = ymaxb
        xyzm(6,2) = zmaxb
C
        aaa = sqrt(nmn /
     .           ((xmaxb-xminb)*(ymaxb-yminb)
     .           +(ymaxb-yminb)*(zmaxb-zminb)
     .           +(zmaxb-zminb)*(xmaxb-xminb)))
        aaa = 0.75*aaa
 
        nbx = nint(aaa*(xmaxb-xminb))
        nby = nint(aaa*(ymaxb-yminb))
        nbz = nint(aaa*(zmaxb-zminb))
        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))
        ENDIF
 
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------------------
  200 CONTINUE
C------------------
!$OMP PARALLEL        
        CALL i7trivox1(
     1   nsn      ,i_mem   ,irect    ,x       ,stf     ,
     2   stfn     ,xyzm    ,nsv      ,
     3   ncontact ,noint   ,tzinf_st ,gap_s_l ,gap_m_l ,
     4   voxel1   ,nbx     ,nby      ,nbz     ,nrtm_l  ,
     5   igap     ,gap     ,gap_s    ,gap_m   ,gapmin  ,
     6   gapmax   ,marge_st,curv_max ,bgapsmx ,istf    ,
     7   i_stok   ,nin,
     8   id       ,titr    ,drad     ,index   ,
     9   iremnode,flagremnode,kremnode,remnode,
     1   dgapload,ipari,intbuf_tab,
     2   iix,iiy,iiz,local_next_nod,nrtm,is_used_with_law151 )
!$OMP END PARALLEL
C------------------
 
        IF (i_mem == 2)THEN
          RETURN
        ENDIF
C     I_MEM = 1 ==> PAS ASSEZ DE MEMOIRE PILE
C     I_MEM = 2 ==> PAS ASSEZ DE MEMOIRE CANDIDATES
        IF(i_mem==1)THEN
          i_mem = 0
          GO TO 100
        ELSE IF(i_mem==2) THEN
          marge_st = three_over_4*marge_st
          tzinf_st = marge_st + max(gap,drad)
          i_mem = 0
          IF(marge_st<em03) THEN
            CALL ancmsg(msgid=83,
     .                msgtype=msgerror,
     .                anmode=aninfo,
     .                i1=id,
     .                c1=titr)
          ENDIF
          GO TO 100
        ENDIF
C---------------------------------
      END DO ! end loop over NSPMD domains
      
      IF(.NOT.type18)THEN
        IF(nsn/=0)THEN
         WRITE(iout,*)' POSSIBLE IMPACT NUMBER:',i_stok,' (<=',
     .   1+(i_stok-1)/nsn,'*NSN)'
 
C
        ELSE
         CALL ancmsg(msgid=552,
     .               msgtype=msgwarning,
     .               anmode=aninfo_blind_2,
     .               i1=id,
     .               c1=titr)
        ENDIF
      endif!(.NOT.TYPE18)
C
      DEALLOCATE( index )
      DEALLOCATE( curv_max )
      DEALLOCATE(local_next_nod)
      DEALLOCATE(iix)
      DEALLOCATE(iiy)
      DEALLOCATE(iiz)
      RETURN
      END