i11remlin.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/.
!||====================================================================
!||    i11remline    ../starter/source/interfaces/inter3d1/i11remlin.F
!||--- called by ------------------------------------------------------
!||    inint3        ../starter/source/interfaces/inter3d1/inint3.F
!||--- calls      -----------------------------------------------------
!||    i11pene_lin   ../starter/source/interfaces/inter3d1/i11remlin.f
!||    origin        ../starter/source/model/remesh/build_admesh.F
!||====================================================================
      SUBROUTINE i11remline(
     1         X,NRTM,IRECTM,NRTS,IRECTS,
     2         NUMNOD,GAP_S ,GAP_M, GAPMIN,IGAP,
     3         KREMNODE,REMNODE,GAP,DRAD,NREMNODE,
     4         I_START,I_MEM_REM,INOD2LIN,TAGSECND,NOD2LIN,
     5         DGAPLOAD,GAP_S_L,GAP_M_L)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NRTM,NRTS, NUMNOD,IGAP, NREMNODE
      INTEGER IRECTM(2,*),IRECTS(2,*),KREMNODE(*),REMNODE(*),I_START,I_MEM_REM
      INTEGER INOD2LIN(NUMNOD+1),TAGSECND(NUMNOD),NOD2LIN(2*NRTM)
      my_real
     .        X(3,*),GAP_S(*),GAP_M(*),GAP,DRAD,GAPMIN
      my_real , INTENT(IN) :: dgapload
      my_real , INTENT(IN) :: gap_s_l(nrts), gap_m_l(nrtm)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,K,LIN,ILIN,LEVEL,CPT,NBLIN,LIN1,L,CPT1,N,NBLIN_MAX,CPT_TOTAL
      INTEGER, dimension(:),ALLOCATABLE :: ITAG,LISTLIN,LISTLINTMP,LISTLINTOTAL
      INTEGER ::  IM1,IM2
      INTEGER, DIMENSION(:),ALLOCATABLE ::
     .        KNOD2LIN,TAGNOD,ORIGIN
      my_real
     .        dmax,new_dist,pene,i11pene_lin,xl,gapv
      my_real, DIMENSION(:),ALLOCATABLE ::
     .        dist1
C-----------------------------------------------
c Build inverse connectivity for segments - only at first pass (I_START=1)
C-----------------------------------------------
      ALLOCATE(itag(nrtm))
      ALLOCATE(listlin(nrtm))
      ALLOCATE(listlintmp(nrtm))
      ALLOCATE(listlintotal(nrtm))
C
      cpt = 0
      IF (i_start ==1) THEN
C
        ALLOCATE(knod2lin(numnod+1))
C
        kremnode(1) = 1
        nod2lin(1:2*nrtm)    = 0
        knod2lin(1:numnod+1) = 0
        inod2lin(1:numnod+1) = 0
        tagsecnd(1:numnod)    = 0
C
        DO i=1,nrts
          tagsecnd(irects(1,i)) = 1
          tagsecnd(irects(2,i)) = 1
        ENDDO
C
        DO i=1,nrtm
          DO j=1,2
            IF( tagsecnd(irectm(j,i)) == 1 ) cpt = cpt + 1
          ENDDO
        ENDDO
C
        IF (cpt == 0) THEN
C--
          DO i=1,nrtm
            kremnode(i+1) = 0
          ENDDO
C
        ELSE
C-----------------------------------------------
C  Definition of node to segment connections
C-----------------------------------------------
C
          DO i=1,nrtm
            cpt = 0
            DO k=1,2
              IF(tagsecnd(irectm(k,i)) == 1) cpt = cpt + 1
            END DO
            IF (cpt /= 0 ) THEN
              DO k=1,2
                n = irectm(k,i)
                knod2lin(n) = knod2lin(n) + 1
              END DO
            ENDIF
          END DO
C
          inod2lin(1) = 1
          DO i=1,numnod
            inod2lin(i+1) = inod2lin(i) + knod2lin(i)
          END DO
          knod2lin(1:numnod+1) = inod2lin(1:numnod+1)
C
          DO i=1,nrtm
            cpt = 0
            DO k=1,2
              IF(tagsecnd(irectm(k,i)) == 1) cpt = cpt + 1
            END DO
            IF (cpt /= 0) THEN
              DO k=1,2
                n = irectm(k,i)
                nod2lin(knod2lin(n)) = i
                knod2lin(n) = knod2lin(n) + 1
              END DO
            ENDIF
          END DO
C
          DEALLOCATE(knod2lin)
C
        ENDIF
C
      ENDIF
C
      IF (((i_start==1).AND.(cpt > 0)).OR.(i_start>1)) THEN
C
C-----------------------------------------------
C  Searching Algorithm Connected nodes : D < SQRT(2.) * GAP
C-----------------------------------------------
C
        ALLOCATE(tagnod(numnod),origin(numnod),dist1(numnod))
        tagnod(1:numnod) = 0
        origin(1:numnod)     = 0
        dist1(1:numnod) = ep30
        itag(1:nrtm) = 0
        listlin(1:nrtm) = 0
        listlintmp(1:nrtm)=0
        listlintotal(1:nrtm) = 0
        cpt_total = 0
C
        dmax = sqrt(two) * max(gap+dgapload,drad)
C
        DO i=i_start,nrtm
 
          level = 1
          lin = i
C          IF ((ITAB(IRECTM(1,LIN))/=30151).OR.(ITAB(IRECTM(2,LIN))/=30197)) CYCLE
          itag(lin) = level
          listlin(1)=lin
          nblin=1
          nblin_max=1
          cpt = 0
          cpt_total = 0
          xl = (x(1,irectm(1,i))-x(1,irectm(2,i)))**2+(x(2,irectm(1,i))-x(2,irectm(2,i)))**2+(x(3,irectm(1,i))-x(3,irectm(2,i)))**2
          xl = sqrt(xl)
C
          DO j=1,2
            tagnod(irectm(j,lin)) = 1
            dist1(irectm(j,lin)) = zero
          ENDDO
C
          DO WHILE (nblin/=0)
C
            level = level+1
            cpt = 0
            DO ilin=1,nblin
              lin=listlin(ilin)
              tagnod(irectm(1:2,lin))=2
C
C              ESTA = (DIST1(IRECTM(2,LIN))*DIST1(IRECTM(2,LIN))-DIST1(IRECTM(1,LIN))*DIST1(IRECTM(1,LIN))-XL*XL)
C     .               /(TWO*XL*DIST1(IRECTM(1,LIN)))
C              DIST_AXIS = DIST1(IRECTM(1,LIN))*SQRT(ONE-ESTA*ESTA)
C
              pene = zero
              IF ((dist1(irectm(1,lin)) > dmax).AND.(dist1(irectm(2,lin)) > dmax).AND.(level>2)) THEN
                pene = i11pene_lin(x,irectm(1,lin),irectm(2,lin),irectm(1,i),irectm(2,i),dmax)
              ENDIF
C
              IF ((level <= 2).OR.(dist1(irectm(1,lin)) <= dmax).OR.(dist1(irectm(2,lin)) <= dmax).OR.(pene > zero)) THEN
              DO j=1,2
                DO k=inod2lin(irectm(j,lin)),inod2lin(irectm(j,lin)+1)-1
                  lin1 = nod2lin(k)
                  IF( (itag(lin1) == 0  .OR. itag(lin1) == level)) THEN
                    IF(itag(lin1) == 0)THEN
                      cpt = cpt + 1
                      listlintmp(cpt)=lin1
                    ENDIF
                    itag(lin1)=level
                    DO l=1,2
 
                      IF ((tagsecnd(irectm(l,lin1))== 1).AND.(origin(irectm(l,lin1)) /= irectm(j,lin))
     .                     .AND.((irectm(l,lin1)) /= irectm(j,lin)).AND.(tagnod(irectm(l,lin1)) /= 2)) THEN
C
                        new_dist=dist1(irectm(j,lin))+
     .                  sqrt((x(1,irectm(l,lin1))-x(1,irectm(j,lin)))**2 + 
     .                       (x(2,irectm(l,lin1)) - x(2,irectm(j,lin)))**2 +
     .                       (x(3,irectm(l,lin1)) - x(3,irectm(j,lin)))**2 )
C
                        IF (new_dist < dist1(irectm(l,lin1))) THEN
                          dist1(irectm(l,lin1)) = new_dist
                        ENDIF
C
                        IF(tagnod(irectm(l,lin1))==0) THEN
                          tagnod(irectm(l,lin1)) = 1
                        ENDIF                     
C
                      ENDIF
                    ENDDO 
                  ENDIF                  
                ENDDO   
              ENDDO
              ENDIF
C
              tagnod(irectm(1:2,lin))=1
            ENDDO
C
            nblin = cpt
C 
            nblin_max = max(nblin_max,nblin)
            IF(nblin ==0)EXIT
            DO j=1,cpt
              listlin(j)=listlintmp(j)
              listlintmp(j) = 0
              listlintotal(j+cpt_total) = listlin(j)
            ENDDO
            cpt_total = cpt_total + cpt
C
C----------------
          ENDDO
C
CC END DO WHILE
C
C-- Check memory for data storage
C
          i_start = i
          IF (kremnode(i)+cpt_total > nremnode) THEN
C--         Not enough memory - upgrade_remnode
            i_mem_rem = 1
            EXIT
          ENDIF
C
          cpt1 = 0
          im1 = irectm(1,i)
          im2 = irectm(2,i)
C
          IF (igap == 0) THEN 
            DO l=1,cpt_total
              lin = listlintotal(l)
              IF ((im1 /= irectm(1,lin)).AND.(im1 /= irectm(2,lin))
     .           .AND.(im2 /= irectm(1,lin)).AND.(im2 /= irectm(2,lin))) THEN
C--- lines with common nodes with main lines are already removed - no need to store them in remnode
                IF ((dist1(irectm(1,lin)) <= dmax).OR.(dist1(irectm(2,lin)) <= dmax)) THEN
                  remnode(kremnode(i)+cpt1) = lin
                  cpt1 = cpt1 + 1
                ELSE
                  pene = i11pene_lin(x,irectm(1,lin),irectm(2,lin),irectm(1,i),irectm(2,i),sqrt(two)*max(gap+dgapload,drad))
                  IF (pene > 0) THEN
                    remnode(kremnode(i)+cpt1) = lin
                    cpt1 = cpt1 + 1
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
            kremnode(i+1) = kremnode(i) + cpt1
          ELSE 
            DO l=1,cpt_total
              lin = listlintotal(l)
              IF ((im1 /= irectm(1,lin)).AND.(im1 /= irectm(2,lin))
     .           .AND.(im2 /= irectm(1,lin)).AND.(im2 /= irectm(2,lin))) THEN
C--- lines with common nodes with main lines are already removed - no need to store them in remnode
                gapv = gap_s(lin)+gap_m(i)
                IF(igap == 3) gapv=min(gap_s_l(lin)+gap_m_l(i),gapv)
                gapv = sqrt(two)*max(drad,gapmin,gapv+dgapload)
                IF ((dist1(irectm(1,lin)) <= gapv).OR.(dist1(irectm(2,lin)) <= gapv)) THEN
                  remnode(kremnode(i)+cpt1) = lin
                  cpt1 = cpt1 + 1
                ELSE
                  pene = i11pene_lin(x,irectm(1,lin),irectm(2,lin),irectm(1,i),irectm(2,i),gapv)
                  IF (pene > 0) THEN
                    remnode(kremnode(i)+cpt1) = lin
                    cpt1 = cpt1 + 1
                  ENDIF
                ENDIF
              ENDIF
            ENDDO
            kremnode(i+1) = kremnode(i) + cpt1
          ENDIF
C
C-----------------------------------------------
C  Clean of used arrays
C-----------------------------------------------
C
          dist1(irectm(1,i)) = ep30
          dist1(irectm(2,i)) = ep30
          origin(irectm(1,i)) = 0
          origin(irectm(2,i)) = 0
          tagnod(irectm(1,i)) = 0
          tagnod(irectm(2,i)) = 0
          itag(i) = 0
C
          DO l=1,cpt_total
            lin = listlintotal(l)
            itag(lin) = 0
            listlintotal(l) = 0
            tagnod(irectm(1,lin)) = 0
            tagnod(irectm(2,lin)) = 0
            dist1(irectm(1,lin)) = ep30
            dist1(irectm(2,lin)) = ep30
            origin(irectm(1,lin)) = 0
            origin(irectm(2,lin)) = 0
          ENDDO
          listlintmp(1:nblin_max)=0
          listlin(1:nblin_max)=0
C
        ENDDO
CC END DO NRTM
        DEALLOCATE(dist1,tagnod,origin)
C
      ELSE
        i_start = nrtm  ! avoid infinite loop later
      ENDIF
 
      DEALLOCATE(itag)
      DEALLOCATE(listlin)
      DEALLOCATE(listlintmp)
      DEALLOCATE(listlintotal)
C
      RETURN
      END
C
C=======================================================================
!||====================================================================
!||    i11pene_lin   ../starter/source/interfaces/inter3d1/i11remlin.F
!||--- called by ------------------------------------------------------
!||    i11remline    ../starter/source/interfaces/inter3d1/i11remlin.F
!||    i25remline    ../starter/source/interfaces/int25/i25remlin.F
!||====================================================================
      my_real FUNCTION i11pene_lin(X,N1,N2,M1,M2,GAP)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER n1,n2,m1,m2
      my_real
     .     x(3,*),gap
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real
     .     xs12,ys12,zs12,xm12,ym12,zm12,xa,xb,
     .     xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
     .     xx,yy,zz,als,alm,det,h1s,h2s,h1m,h2m,nx,ny,nz,gap2
C-----------------------------------------------
C
C      COMPTUTATION OF PENE*PENE FOR PENETRATION CHECK
C
       gap2 = gap**2
C
       xs12 = x(1,n2)-x(1,n1)
       ys12 = x(2,n2)-x(2,n1)
       zs12 = x(3,n2)-x(3,n1)
       xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
       xm12 = x(1,m2)-x(1,m1)
       ym12 = x(2,m2)-x(2,m1)
       zm12 = x(3,m2)-x(3,m1)
       xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
       xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
       xs2m2 = x(1,m2)-x(1,n2)
       ys2m2 = x(2,m2)-x(2,n2)
       zs2m2 = x(3,m2)-x(3,n2)
       xa =  xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
       xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2 
       det = xm2*xs2 - xsm*xsm
       det = max(em20,det)
C
       h1m = (xa*xsm-xb*xs2) / det
C
       xs2 = max(xs2,em20)
       xm2 = max(xm2,em20)
       h1m=min(one,max(zero,h1m))
       h1s = -(xa + h1m*xsm) / xs2
       h1s=min(one,max(zero,h1s))
       h1m = -(xb + h1s*xsm) / xm2
       h1m=min(one,max(zero,h1m))
C
       h2s = one - h1s
       h2m = one - h1m
C !!!!!!!!!!!!!!!!!!!!!!!
C  PENE = GAP^2 - DIST^2 UTILISE POUR TESTER SI NON NUL
C!!!!!!!!!!!!!!!!!!!!!!!!
       nx =  h1s*x(1,n1) + h2s*x(1,n2)
     .      - h1m*x(1,m1) - h2m*x(1,m2)
       ny =  h1s*x(2,n1) + h2s*x(2,n2)
     .      - h1m*x(2,m1) - h2m*x(2,m2)
       nz =  h1s*x(3,n1) + h2s*x(3,n2)
     .      - h1m*x(3,m1) - h2m*x(3,m2)
       i11pene_lin = gap2 - nx*nx - ny*ny - nz*nz
       i11pene_lin = max(zero,i11pene_lin)
C
       RETURN
       END