i24pen3.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/.
C
!||====================================================================
!||    i24pen3                ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- called by ------------------------------------------------------
!||    inint3                 ../starter/source/interfaces/inter3d1/inint3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||    i24fic_getn            ../starter/source/interfaces/inter3d1/i24surfi.F
!||    i24penmax              ../starter/source/interfaces/inter3d1/i24pen3.F
!||    iconnet                ../starter/source/interfaces/inter3d1/i24pen3.F
!||    ini_st3                ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- uses       -----------------------------------------------------
!||    format_mod             ../starter/share/modules1/format_mod.F90
!||    message_mod            ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE i24pen3(X    ,IRECT ,GAPV ,CAND_E,CAND_N,
     2                   NSV  ,INACTI,ITAB ,TAG   ,IWPENE,
     3                   NSN  ,IRTLM ,MSEGTYP   ,IWPENE0 ,
     4                   PMIN ,GAP_N ,MVOISN    ,IXS     ,
     5                   IXS10,IXS16 ,IXS20,PENMAX,PENMIN,
     6                   ID,TITR ,ILEV ,PEN_OLD,KNOD2ELS,
     7                   NOD2ELS,IPARTNS,IPEN0 ,ICONT_I ,
     8                   XFIC   ,NRTM ,IRTSE   ,IS2SE)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
#ifndef HYPERMESH_LIB
      USE message_mod
#endif
      USE names_and_titles_mod , ONLY : nchartitle
      USE format_mod , ONLY : fmt_i_3f
C
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      "vect07_c.inc"
#include      "scr05_c.inc"
#include      "com04_c.inc"
C-----------------------------------------------
      INTEGER IWPENE,TAG(*),INACTI,NSV(*),NSN,MSEGTYP(*),IWPENE0,
     .        MVOISN(4,*),ILEV,KNOD2ELS(*),NOD2ELS(*),IPARTNS(*),NRTM
C----  IWPENE0 : nb of tiny initial pene;  IWPENE nb of initial pene
      my_real GAPV(*)
      INTEGER IRECT(4,*), ITAB(*),CAND_E(*),CAND_N(*),IRTLM(2,*)
      INTEGER IXS(NIXS,*),IXS10(6,*), IXS16(8,*), IXS20(12,*),ICONT_I(*),
     .        IRTSE(*) ,IS2SE(*)
      my_real x(3,*),pmin(*),gap_n(12,*),penmax,penmin,pen_old(5,nsn),xfic(3,*)
C--------GAP_N(1,*) stock temporarily characteristic length for Pen_max
      INTEGER ID,IPEN0
      CHARACTER(LEN=NCHARTITLE) :: TITR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER II, I, J, K, L, JJ, NJ, IER,NS,IC,I0,IELIM,NI,ICONN,ip,NS1,
     .        IDEL,NN1,NN2,IE
C     REAL
      my_real
     .   pen, alp,xx(4),yy(4),zz(4),ssc,ttc,dist,dist0,
     .   xi,yi,zi,xc,yc,zc,nn(3),tol,pen0,dpen,norm,maxpen
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
C---  MVOISN,IRTLM(2,*) is used temporarily for Pen_ini MVOISN(1,*) -> MTYPE(solid),MVOISN(2,*) -> E_id
C-----MVOISN(3,*) -> part_id, PEN_OLD(1-3,*)->SECONDARY nodal normal, reput PEN_OLD(1-2,*)=0, 3,)=5,)after
C-----to be consistent with engine for PENMIN
      IF (iresp==1.AND.penmin<=em06) penmin = two*em06
          tol = penmin
      alp = two*em06
      IF (iresp==1) alp = two*em05
      DO i=lft,llt
        l  = cand_e(i)
        ni = cand_n(i)         
        ns = nsv(ni)
        IF (ns >numnod) THEN
         ns1 = ns -numnod
         xi=xfic(1,ns1)
         yi=xfic(2,ns1)
         zi=xfic(3,ns1)
        ELSE
         xi=x(1,ns)
         yi=x(2,ns)
         zi=x(3,ns)
        END IF
        DO jj=1,4
         nj=irect(jj,l)
         xx(jj)=x(1,nj)
         yy(jj)=x(2,nj)
         zz(jj)=x(3,nj)
        END DO
C--------
        CALL ini_st3(xx,yy,zz,xi,yi,zi,nn,ssc,ttc,ier,alp,
     2               xc,yc,zc)
        IF(ier==-1)THEN
#ifndef HYPERMESH_LIB
         CALL ancmsg(msgid=85,
     .               msgtype=msgerror,
     .               anmode=aninfo,
     .               i1=id,
     .               c1=titr,
     .               i2=itab(ns),
     .               i3=l,
     .               i4=l,
     .               i5=itab(irect(1,l)),
     .               i6=itab(irect(2,l)),
     .               i7=itab(irect(3,l)),
     .               i8=itab(irect(4,l)))
#endif
C
        ELSE IF(ier==1.AND.(msegtyp(l)/=0.AND.msegtyp(l)<=nrtm))THEN
C                             shells except coating shells
C---------outside        
c         IF(IPRI>=1)WRITE(IOUT,FMT=FMT_6I)ITAB(NS),L,
c     .                    (ITAB(IRECT(JJ,L)),JJ=1,4)
        ELSE
C------initial penetration case 1) |PEN|<TOL : on, 2) Inacti<0 on,3)Inacti=3,4
C-------warnning out for 2),3),4)
C --------MSEGTYP /=0 ->shell ---
          pen0=nn(1)*(xi-xc)+nn(2)*(yi-yc)+nn(3)*(zi-zc)
          IF(ier==1) THEN
           dist = sqrt((xi-xc)*(xi-xc)+(yi-yc)*(yi-yc)+(zi-zc)*(zi-zc))
          ELSE
           dist = abs(pen0)
          END IF
C--------for exception of elimination---           
          idel = 1
C----------coating shell is like solid now----           
          IF (msegtyp(l)/=0.AND.msegtyp(l)<=nrtm) THEN
           pen=gapv(i)-abs(pen0)
           IF (pen > penmax ) idel = 0
C-debug sandwish shell : avoid elimination with high thickness outside
           IF (pen > zero) dist = abs(gapv(i)-pen0)
C----------give up the the wrong one (normal direction, and)           
           IF (pen0 < zero .OR. pen > penmax) pen=-abs(pen)-tol
C----------distance after shifted           
          ELSE
           pen=gapv(i)-pen0
C------------used only for eliminating wrong contact w/ smaller distance            
           IF(ier==1) pen=-abs(pen)-tol
           IF (pen > zero .OR. abs(pen) < tol) THEN
            maxpen = gap_n(1,l)
            IF (inacti /= 0) maxpen = penmax
            CALL i24penmax(pen,maxpen ,mvoisn(1,l),mvoisn(2,l),
     +                    ns  ,ixs,   ixs10,     ixs16, ixs20   ,
     +                    ielim)
            iconn = 0
            IF (ns>numnod) THEN
              CALL i24fic_getn(ns1     ,irtse   ,is2se   ,ie    ,nn1     ,
     4                       nn2     )
             CALL iconnet(irect(1,l),ixs   ,knod2els,nod2els,
     .                    ixs10 ,ixs16 ,ixs20 ,nn1  ,iconn   )
             IF (iconn == 0)
     .       CALL iconnet(irect(1,l),ixs   ,knod2els,nod2els,
     .                    ixs10 ,ixs16 ,ixs20 ,nn2  ,iconn   )
            ELSE
             CALL iconnet(irect(1,l),ixs   ,knod2els,nod2els,
     .                    ixs10 ,ixs16 ,ixs20 ,ns   ,iconn   )
            END IF
             IF ((ielim+iconn) > 0) pen = -abs(pen)-tol
             IF (pen < zero ) idel = 0
           END IF
C------Elimine the impact take into account to SECONDARY nodal normal          
           IF (inacti/=0.AND.(pen > zero .OR. abs(pen) < tol).AND.ilev/=3) THEN
            norm = nn(1)*pen_old(1,ni)+nn(2)*pen_old(2,ni)
     +             +nn(3)*pen_old(3,ni)
             IF (norm >= zero) THEN
              pen = -abs(pen)-tol
c             print *,'impact pair eliminated due to N_SECONDARY'
              idel = 0
             END IF
           END IF
          END IF !(MSEGTYP(L)/=0 ) THEN
C------Elimine the impact between the same part
          IF (ipen0==0) THEN           
           IF (inacti/=0.AND.(pen > zero .OR. abs(pen) < tol)) THEN
            IF (ipartns(ni) == mvoisn(3,l)) THEN
             pen = -abs(pen)-tol
            END IF
           END IF
          END IF !(IPEN0==0) THEN
C------don't take into account auto-impact case for elimination          
          IF (ipartns(ni) == mvoisn(3,l)) idel = 0
C--------exception for SECONDARY shell (test TV)----          
          IF (gapv(i)>zero.AND.(msegtyp(l)==0.OR.msegtyp(l)>nrtm))idel=0
C--------PMIN() has been changed from PENE to dist excepting (INACTI ==0,1)        
C--------PMIN() = -dist for  ABS(PEN) < TOL .OR. PEN<ZERO       
C------------ cas 1  this part is removed in Engine at T=0 for consisting
          IF(abs(pen) < tol .OR. (pen<zero.AND.idel>0)) THEN
C---------only used to calculate Dist_min and to eliminate wrong contact (higher)            
           IF (tag(ns)==0) THEN
            pmin(ni)=-dist
            tag(ns)=ni
           ELSE
            i0=tag(ns)
            pen0=pmin(i0)
            IF (dist <abs(pen0)) THEN
C----------only update dist_min            
              pmin(ni)=-dist
              tag(ns)=ni
             IF (pen0 > zero) THEN
C----------elimine wrong contact            
              irtlm(1,i0)=0
              irtlm(2,i0)=0
              pen_old(5,i0)=zero
             END IF
            END IF
           END IF
          ELSEIF(pen > penmax) THEN
C----------warning w/o treatment
#ifndef HYPERMESH_LIB
            WRITE(iout,1200)pen
#endif
          ELSEIF(pen > zero) THEN
C------Warning anyway-------------
           IF (tag(ns)==0) iwpene=iwpene+1
C------------PMIN has been changed from PENE to dist excepting (INACTI ==0,1)        
           IF(inacti ==0 .OR. inacti ==1) THEN
C------------use IRTLM(2,NI)-> ICONT_I<0 for initial penetration
            IF (tag(ns)>0) THEN
             i0=tag(ns)
             pen0=pmin(i0)
C----------exclude case of PMIN(I0)<0 : -dist            
             IF (pen < pen0) THEN
              icont_i(ni)=-l
              pmin(ni)=pen
              tag(ns) = ni
#ifdef HYPERMESH_LIB
             pen_old(1:3,ni) = nn(1:3)
#endif
             ENDIF             
            ELSE
             icont_i(ni)=-l
             pmin(ni)=pen
             tag(ns) = ni
#ifdef HYPERMESH_LIB
             pen_old(1:3,ni) = nn(1:3)
#endif
            END IF
           ELSEIF(inacti ==-1) THEN
C------------multi-cont-> single by overwriting with min_pene
            IF (tag(ns)>0) THEN
             i0=tag(ns)
             pen0=pmin(i0)
             dist0 = abs(pmin(i0))
             IF (dist < dist0) THEN
              irtlm(1,ni)=l
              irtlm(2,ni)=1
              pmin(ni)=dist
              pen_old(5,ni)=pen
              tag(ns) = ni
#ifdef HYPERMESH_LIB
              pen_old(1:3,ni) = nn(1:3)
#endif
             ENDIF             
            ELSE
             irtlm(1,ni)=l
             irtlm(2,ni)=1
             pmin(ni)=dist
             pen_old(5,ni)=pen
             tag(ns) = ni
#ifdef HYPERMESH_LIB
             pen_old(1:3,ni) = nn(1:3)
#endif
            END IF
C--------hide option------           
           ELSEIF(inacti ==3 ) THEN
            IF (ilev ==3) THEN
             dpen = pen + tol
            ELSE
             dpen = half*(pen + tol)
            END IF
C-------change SECONDARY node
            IF (tag(ns)==0) THEN
              irtlm(1,ni)=l
              irtlm(2,ni)=1
              iwpene=iwpene+1
              tag(ns)=ni
#ifndef HYPERMESH_LIB
             WRITE(iout,1000)pen
#endif
             IF (ns >numnod) THEN
               ns1 = ns -numnod
               xfic(1,ns1) = xi + dpen*nn(1)
               xfic(2,ns1) = yi + dpen*nn(2)
               xfic(3,ns1) = zi + dpen*nn(3)
#ifndef HYPERMESH_LIB
              WRITE(iout,fmt=fmt_i_3f)(itab(numnod)+ns1),xfic(1,ns1),xfic(2,ns1),xfic(3,ns1)
#endif
             ELSE
              x(1,ns) = xi + dpen*nn(1)
              x(2,ns) = yi + dpen*nn(2)
              x(3,ns) = zi + dpen*nn(3)
#ifndef HYPERMESH_LIB
              WRITE(iout,fmt=fmt_i_3f)itab(ns),x(1,ns),x(2,ns),x(3,ns)
#endif
             END IF !(NS >NUMNOD) THEN
            END IF
           ELSEIF(inacti ==5) THEN
C------------multi-cont-> single by overwriting with min_pene         
            IF (tag(ns)>0) THEN
             i0=tag(ns)
             pen0=pen_old(5,i0)
             dist0 = abs(pmin(i0))
             IF (dist < dist0) THEN
              irtlm(1,ni)=l
              irtlm(2,ni)=1
              pen_old(5,ni)=pen
              pmin(ni)=dist
              tag(ns) = ni
#ifdef HYPERMESH_LIB
              pen_old(1:3,ni) = nn(1:3)
#endif
             ENDIF             
            ELSE
             irtlm(1,ni)=l
             irtlm(2,ni)=1
             pen_old(5,ni)=pen
             pmin(ni)=dist
             tag(ns) = ni
#ifdef HYPERMESH_LIB
             pen_old(1:3,ni) = nn(1:3)
#endif
            END IF
           END IF !IF(INACTI ==0 .OR. INACTI ==1) THEN
          END IF !(pen > zero) THEN
        END IF !(IER==-1)THEN
      END DO !I=LFT,LLT              
C      
      RETURN
 1000 FORMAT(2x,'** INITIAL PENETRATION =',1pg20.13,
     . ' CHANGE COORDINATES OF SECONDARY NODE TO:')
 1100 FORMAT(2x,'** INITIAL PENETRATION =',1pg20.13,
     . ' CHANGE COORDINATES OF MAIN NODE TO:')
 1200 FORMAT(2x,'** TOO HIGH INITIAL PENETRATION=, WILL BE IGNORED',
     .       1pg20.13)
C
      END
!||====================================================================
!||    ini_st3   ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- called by ------------------------------------------------------
!||    i24pen3   ../starter/source/interfaces/inter3d1/i24pen3.F
!||====================================================================
      SUBROUTINE ini_st3(XX,YY,ZZ,XI,YI,ZI,NN,SSC,TTC,IER,ALP,
     1                   XC,YC,ZC)
C
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 IER
      my_real
     .   XX(4),YY(4),ZZ(4),NN(3), SSC, TTC, ALP,XI,YI,ZI,XC,YC,ZC
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real
     .   H(4), X0, Y0, Z0, XL1, XL2, XL3, XL4, YY1, YY2, YY3, YY4,
     .   zz1, zz2, zz3, zz4, xi1, xi2, xi3, xi4, yi1, yi2, yi3, yi4,
     .   zi1, zi2, zi3, zi4, xn1, yn1, zn1, xn2, yn2, zn2, xn3, yn3,
     .   zn3, xn4, yn4, zn4, an, area, a12, a23, a34, a41, b12, b23,
     .   b34, b41, ab1, ab2, tp, tm, sp, sm, x1,x2,x3,x4,
     .   y1,y2,y3,y4,z1,z2,z3,z4,n1,n2,n3,la,lb,lc,lbs,lcs,tt1,ss1
C
      x1=xx(1)
      x2=xx(2)
      x3=xx(3)
      x4=xx(4)
      y1=yy(1)
      y2=yy(2)
      y3=yy(3)
      y4=yy(4)
      z1=zz(1)
      z2=zz(2)
      z3=zz(3)
      z4=zz(4)
C
      x0 = fourth*(x1+x2+x3+x4)
      y0 = fourth*(y1+y2+y3+y4)
      z0 = fourth*(z1+z2+z3+z4)
C
      xl1 = x1-x0
      xl2 = x2-x0
      xl3 = x3-x0
      xl4 = x4-x0
      yy1 = y1-y0
      yy2 = y2-y0
      yy3 = y3-y0
      yy4 = y4-y0
      zz1 = z1-z0
      zz2 = z2-z0
      zz3 = z3-z0
      zz4 = z4-z0
C
      xi1 = x1-xi
      xi2 = x2-xi
      xi3 = x3-xi
      xi4 = x4-xi
      yi1 = y1-yi
      yi2 = y2-yi
      yi3 = y3-yi
      yi4 = y4-yi
      zi1 = z1-zi
      zi2 = z2-zi
      zi3 = z3-zi
      zi4 = z4-zi
C
      xn1 = yy1*zz2 - yy2*zz1
      yn1 = zz1*xl2 - zz2*xl1
      zn1 = xl1*yy2 - xl2*yy1
      n1=xn1
      n2=yn1
      n3=zn1
C
      xn2 = yy2*zz3 - yy3*zz2
      yn2 = zz2*xl3 - zz3*xl2
      zn2 = xl2*yy3 - xl3*yy2
      n1=n1+xn2
      n2=n2+yn2
      n3=n3+zn2
C
      xn3 = yy3*zz4 - yy4*zz3
      yn3 = zz3*xl4 - zz4*xl3
      zn3 = xl3*yy4 - xl4*yy3
      n1=n1+xn3
      n2=n2+yn3
      n3=n3+zn3
C
      xn4 = yy4*zz1 - yy1*zz4
      yn4 = zz4*xl1 - zz1*xl4
      zn4 = xl4*yy1 - xl1*yy4
      n1=n1+xn4
      n2=n2+yn4
      n3=n3+zn4
C
      an= max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      n1=n1/an
      n2=n2/an
      n3=n3/an
      nn(1)=n1
      nn(2)=n2
      nn(3)=n3
      IF(an<=em19) THEN
        ier=-1
        RETURN
      ENDIF
      area=half*an
C
      a12=(n1*xn1+n2*yn1+n3*zn1)
      a23=(n1*xn2+n2*yn2+n3*zn2)
      a34=(n1*xn3+n2*yn3+n3*zn3)
      a41=(n1*xn4+n2*yn4+n3*zn4)
C
      xn1 = yi1*zi2 - yi2*zi1
      yn1 = zi1*xi2 - zi2*xi1
      zn1 = xi1*yi2 - xi2*yi1
      b12=(n1*xn1+n2*yn1+n3*zn1)
C
      xn2 = yi2*zi3 - yi3*zi2
      yn2 = zi2*xi3 - zi3*xi2
      zn2 = xi2*yi3 - xi3*yi2
      b23=(n1*xn2+n2*yn2+n3*zn2)
C
      xn3 = yi3*zi4 - yi4*zi3
      yn3 = zi3*xi4 - zi4*xi3
      zn3 = xi3*yi4 - xi4*yi3
      b34=(n1*xn3+n2*yn3+n3*zn3)
C
      xn4 = yi4*zi1 - yi1*zi4
      yn4 = zi4*xi1 - zi1*xi4
      zn4 = xi4*yi1 - xi1*yi4
      b41=(n1*xn4+n2*yn4+n3*zn4)
C
      ab1=a23*b41
      ab2=b23*a41
C
      IF(abs(ab1+ab2)/area>em10)THEN
       ssc=(ab1-ab2)/(ab1+ab2)
      ELSE
       ssc=zero
      ENDIF
      IF(abs(a34/area)>em10)THEN
       ab1=b12*a34
       ab2=b34*a12
       IF(abs(ab1+ab2)/area>em10)THEN
        ttc=(ab1-ab2)/(ab1+ab2)
       ELSE
        ttc=zero
       END IF
      ELSE
       ttc=b12/a12-one
       IF(b23<=zero.AND.b41<=zero)THEN
         IF(-b23/a12<=alp.AND.-b41/a12<=alp)ssc=zero
       ELSEIF(b23<=zero)THEN
         IF(-b23/a12<=alp) THEN
          ssc=one
         ELSE
          ssc=two
         END IF
       ELSEIF(b41<=zero)THEN
         IF(-b41/a12<=alp) THEN
          ssc=-one
         ELSE
          ssc=-two
         END IF
       ENDIF
      ENDIF
C-------------out of seg
      IF(abs(ssc)>one+alp.OR.abs(ttc)>one+alp) THEN
        ier=1
C------case tria  re-compute      
        IF (a34==zero.AND.ttc< one) THEN
             lb=fourth*(one - ttc)*(one - ssc)
             lc=fourth*(one - ttc)*(one + ssc)
             la = one - lb - lc
             IF (la>=zero) THEN
               lb= min(one,max(zero,lb))
               lc= min(one,max(zero,lc))
             ELSEIF(lc>lb.AND.lc >= one) THEN
               lc = one
               lb = zero
             ELSEIF(lb >= one) THEN
               lc = zero
               lb = one
             ELSE             
               lbs = half*(one+lb-lc)
               lcs = half*(one-lb+lc)
               lb= min(one,max(zero,lbs))
               lc= min(one,max(zero,lcs))
             ENDIF
             ssc= (lc-lb)/(lc+lb)
             ttc= one - two*lb - two*lc 
        END IF
        IF(abs(ssc)>one)ssc=ssc/abs(ssc)
        IF(abs(ttc)>one)ttc=ttc/abs(ttc)
      ELSE
        ier=0
        IF(abs(ssc)>one)ssc=ssc/abs(ssc)
        IF(abs(ttc)>one)ttc=ttc/abs(ttc)
      ENDIF
C
      tp=fourth*(one+ttc)
      tm=fourth*(one-ttc)
C
      sp=one+ssc
      sm=one-ssc
      h(1)=tm*sm
      h(2)=tm*sp
      h(3)=tp*sp
      h(4)=tp*sm
C
      xc =h(1)*x1+h(2)*x2+h(3)*x3+h(4)*x4
      yc =h(1)*y1+h(2)*y2+h(3)*y3+h(4)*y4
      zc =h(1)*z1+h(2)*z2+h(3)*z3+h(4)*z4
      RETURN
      END
!||====================================================================
!||    i24penmax   ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- called by ------------------------------------------------------
!||    i24pen3     ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- calls      -----------------------------------------------------
!||    intab       ../starter/source/interfaces/inter3d1/i24tools.F
!||====================================================================
      SUBROUTINE i24penmax(PEN  ,PENMAX,ETYP  ,EL    ,NS   ,
     +                     IXS   ,IXS10 ,IXS16 , IXS20  ,IELIM )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "com04_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IXS(NIXS,*),IXS10(6,*), IXS16(8,*), IXS20(12,*)
      INTEGER ETYP  ,EL    ,NS,IELIM
C     REAL
      my_real
     .   PEN  ,PENMAX
C-----------------------------------------------
C   External function
C-----------------------------------------------
      LOGICAL INTAB
      EXTERNAL intab
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
C     REAL
      INTEGER  EL2
      my_real
     .   S
C----add commun ID--at end--------------------------
C--------eliminier some initial penetrations---------
C------automatic for self_contact
       ielim=0
       SELECT CASE (etyp)
         CASE(1)
           IF (intab(8,ixs(2,el),ns)) ielim=1
         CASE(10)
           el2=el-numels8
           IF (intab(8,ixs(2,el),ns).OR.intab(6,ixs10(1,el2),ns))
     +        ielim=1
         CASE(20)
           el2=el-numels8-numels10
           IF (intab(8,ixs(2,el),ns).OR.intab(12,ixs20(1,el2),ns))
     +        ielim=1
         CASE(16)
           el2=el-numels8-numels10-numels20
           IF (intab(8,ixs(2,el),ns).OR.intab(8,ixs16(1,el2),ns))
     +        ielim=1
       END SELECT
C-------
       IF (pen >= penmax ) ielim = 1
 
      RETURN
      END
!||====================================================================
!||    iconnet   ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- called by ------------------------------------------------------
!||    i24pen3   ../starter/source/interfaces/inter3d1/i24pen3.F
!||====================================================================
      SUBROUTINE iconnet(IRECT ,IXS   ,KNOD2ELS,NOD2ELS,
     .                   IXS10 ,IXS16 ,IXS20 ,NS  ,ICONN   )
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-----------------------------------------------
C     REAL
      INTEGER IRECT(4), IXS(NIXS,*), KNOD2ELS(*), NOD2ELS(*), 
     .        IXS10(6,*), IXS16(8,*), IXS20(12,*),ICONN,NS
C     REAL
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER N, JJ, II, K, NN, KK, IC, IAD
C     REAL
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
C
      iconn = 0
      IF(numels==0) RETURN
       DO 230 iad=knod2els(ns)+1,knod2els(ns+1)
        n = nod2els(iad)
        IF(n <= numels8)THEN
          DO jj=1,4
            ii=irect(jj)
            DO k=1,8
              IF(ixs(k+1,n)==ii) iconn = 1
            ENDDO
          ENDDO
        ELSEIF(n <= numels8+numels10)THEN
          DO jj=1,4
            ii=irect(jj)
            DO k=1,8
              IF(ixs(k+1,n)==ii) iconn = 1
            ENDDO
            DO k=1,6
              IF(ixs10(k,n-numels8)==ii) iconn = 1
            ENDDO
          ENDDO
        ELSEIF(n <= numels8+numels10+numels20)THEN
          DO jj=1,4
            ii=irect(jj)
            DO k=1,8
              IF(ixs(k+1,n)==ii) iconn = 1
            ENDDO
            DO k=1,12
              IF(ixs20(k,n-numels8-numels10)==ii) iconn = 1
            ENDDO
          ENDDO
        ELSEIF(n <= numels8+numels10+numels20+numels16)THEN
          DO jj=1,4
            ii=irect(jj)
            DO k=1,8
              IF(ixs(k+1,n)==ii) iconn = 1
            ENDDO
            DO k=1,8
              IF(ixs16(k,n-numels8-numels10-numels20)==ii) iconn = 1
            ENDDO
          ENDDO
        END IF
  230  CONTINUE
      RETURN
      END
C-----------------------------------------------
!||====================================================================
!||    i24cand   ../starter/source/interfaces/inter3d1/i24pen3.F
!||--- called by ------------------------------------------------------
!||    inint3    ../starter/source/interfaces/inter3d1/inint3.F
!||====================================================================
      SUBROUTINE i24cand(CAND_E,CAND_N,NSN  ,IRTLM ,II_STOK ,
     *                    MSEGTYP)
C
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-----------------------------------------------
      INTEGER CAND_E(*),CAND_N(*),NSN,IRTLM(2,*),II_STOK,
     *                    msegtyp(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER E, I,ISH
     .        
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
            II_STOK = 0
            DO i=1,nsn
             e = irtlm(1,i)
             IF (e >0) THEN
              ii_stok =ii_stok + 1
              cand_n(ii_stok) = i
              cand_e(ii_stok) = e
 
              ish = abs(msegtyp(e))
              IF (ish > 0)THEN
                 ii_stok =ii_stok + 1
                 cand_n(ii_stok) = i
                 cand_e(ii_stok) = ish
              ENDIF
 
             END IF
            END DO          
C      
      RETURN
      END