i23dst3.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/.
!||====================================================================
!||    i23dst3   ../starter/source/interfaces/inter3d1/i23dst3.F
!||--- called by ------------------------------------------------------
!||    inint3    ../starter/source/interfaces/inter3d1/inint3.f
!||====================================================================
      SUBROUTINE i23dst3(
     1                  JLT    ,CAND_N ,CAND_E ,IRECT  ,NSV    ,
     2                  GAP_S  ,X      ,MSR    ,PENE   ,IFPEN  ,
     3                  IGAP   ,GAP    ,GAPMAX ,GAPMIN ,GAPV   ,
     4                  GAP_M  )
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JLT,  CAND_N(*), CAND_E(*), IRECT(4,*),
     .        NSV(*), IFPEN(*), IGAP, MSR(*)
      my_real
     .     X(3,*), PENE(*), GAP_S(*), 
     .     gap, gapmin, gapmax, gapv(mvsiz), gap_m(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, L, INPROJ(4,MVSIZ)
      INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ), 
     .        NSVG(MVSIZ)
      my_real
     .     x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
     .     y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
     .     z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
     .     xi(mvsiz), yi(mvsiz), zi(mvsiz), 
     .     lb1(mvsiz), lc1(mvsiz), p1(mvsiz), 
     .     lb2(mvsiz), lc2(mvsiz), p2(mvsiz), 
     .     lb3(mvsiz), lc3(mvsiz), p3(mvsiz), 
     .     lb4(mvsiz), lc4(mvsiz), p4(mvsiz),
     .     la1, la2, la3, la4,
     .     nx1(mvsiz), ny1(mvsiz), nz1(mvsiz),
     .     nx2(mvsiz), ny2(mvsiz), nz2(mvsiz),
     .     nx3(mvsiz), ny3(mvsiz), nz3(mvsiz),
     .     nx4(mvsiz), ny4(mvsiz), nz4(mvsiz),
     .     al1(mvsiz),  al2(mvsiz),  al3(mvsiz),  al4(mvsiz),  
     .     x01(mvsiz),  x02(mvsiz),  x03(mvsiz),  x04(mvsiz),
     .     y01(mvsiz),  y02(mvsiz),  y03(mvsiz),  y04(mvsiz),
     .     z01(mvsiz),  z02(mvsiz),  z03(mvsiz),  z04(mvsiz),
     .     xi1(mvsiz),  xi2(mvsiz),  xi3(mvsiz),  xi4(mvsiz), 
     .     yi1(mvsiz),  yi2(mvsiz),  yi3(mvsiz),  yi4(mvsiz), 
     .     zi1(mvsiz),  zi2(mvsiz),  zi3(mvsiz),  zi4(mvsiz), 
     .     x0(mvsiz),  y0(mvsiz),  z0(mvsiz),
     .     xi0v(mvsiz),  yi0v(mvsiz),  zi0v(mvsiz),
     .     la, hla, h0, h00,
     .     hlb1(mvsiz), hlc1(mvsiz),
     .     hlb2(mvsiz), hlc2(mvsiz),
     .     hlb3(mvsiz), hlc3(mvsiz),
     .     hlb4(mvsiz), hlc4(mvsiz),
     .     xp(mvsiz), yp(mvsiz), zp(mvsiz),
     .     xp1(mvsiz), yp1(mvsiz), zp1(mvsiz),
     .     xp2(mvsiz), yp2(mvsiz), zp2(mvsiz),
     .     xp3(mvsiz), yp3(mvsiz), zp3(mvsiz),
     .     xp4(mvsiz), yp4(mvsiz), zp4(mvsiz),
     .     n11(mvsiz), n21(mvsiz), n31(mvsiz),
     .     n12(mvsiz), n22(mvsiz), n32(mvsiz),
     .     n13(mvsiz), n23(mvsiz), n33(mvsiz),
     .     n14(mvsiz), n24(mvsiz), n34(mvsiz),
     .     nx(mvsiz),  ny(mvsiz),  nz(mvsiz),
     .     xg(mvsiz),  yg(mvsiz),  zg(mvsiz),
     .     h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
     .     nvers(4,mvsiz), firstimp(4,mvsiz),
     .     xn1(mvsiz),yn1(mvsiz),zn1(mvsiz),
     .     xn2(mvsiz),yn2(mvsiz),zn2(mvsiz),
     .     xn3(mvsiz),yn3(mvsiz),zn3(mvsiz),
     .     xn4(mvsiz),yn4(mvsiz),zn4(mvsiz),
     .     ssc(mvsiz),ttc(mvsiz),
     .     a12(mvsiz),a23(mvsiz),a34(mvsiz),a41(mvsiz),
     .     b12(mvsiz),b23(mvsiz),b34(mvsiz),b41(mvsiz),
     .     ab1(mvsiz),ab2(mvsiz),
     .     n1(mvsiz),n2(mvsiz),n3(mvsiz),
     .     xx1(mvsiz),xx2(mvsiz),xx3(mvsiz),xx4(mvsiz),
     .     yy1(mvsiz),yy2(mvsiz),yy3(mvsiz),yy4(mvsiz),
     .     zz1(mvsiz),zz2(mvsiz),zz3(mvsiz),zz4(mvsiz),
     .     area(mvsiz),alp(mvsiz),var
      my_real
     .     s2,nn,ps,
     .     x12,x23,x34,x41,sx1,sx2,sx3,sx4,sx0,
     .     y12,y23,y34,y41,sy1,sy2,sy3,sy4,sy0,
     .     z12,z23,z34,z41,sz1,sz2,sz3,sz4,sz0,
     .     gap2(mvsiz), aaa, side, dist,
     .     ll, h
C--------------------------------------------------------
      IF(igap==0)THEN
        DO i=1,jlt
          gapv(i)=gap
        END DO
      ELSE
        DO i=1,jlt
          gapv(i)=gap_s(cand_n(i))+gap_m(cand_e(i))
          gapv(i)=min(gapv(i),gapmax)
          gapv(i)=max(gapmin,gapv(i))
        END DO
      END IF
C--------------------------------------------------------
C
       DO i=1,jlt
         nsvg(i)= nsv(cand_n(i))
         l      = cand_e(i)
C
         ix1(i)=irect(1,l)
         ix2(i)=irect(2,l)
         ix3(i)=irect(3,l)
         ix4(i)=irect(4,l)
       END DO
C
       DO i=1,jlt
         xi(i)=x(1,nsvg(i))
         yi(i)=x(2,nsvg(i))
         zi(i)=x(3,nsvg(i))
         x1(i)=x(1,ix1(i))
         y1(i)=x(2,ix1(i))
         z1(i)=x(3,ix1(i))
         x2(i)=x(1,ix2(i))
         y2(i)=x(2,ix2(i))
         z2(i)=x(3,ix2(i))
         x3(i)=x(1,ix3(i))
         y3(i)=x(2,ix3(i))
         z3(i)=x(3,ix3(i))
         x4(i)=x(1,ix4(i))
         y4(i)=x(2,ix4(i))
         z4(i)=x(3,ix4(i))
C
         pene(i)=zero
       END DO
C--------------------------------------------------------
C  CASE OF MIXED PACKAGES
C--------------------------------------------------------
       DO i=1,jlt
        IF(ix3(i)/=ix4(i))THEN
         x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
         y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
         z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i)) 
        ELSE
         x0(i) = x3(i)
         y0(i) = y3(i)
         z0(i) = z3(i)
        ENDIF
       ENDDO
C
       DO i=1,jlt
C
        x01(i) = x1(i) - x0(i)
        y01(i) = y1(i) - y0(i)
        z01(i) = z1(i) - z0(i)
C
        x02(i) = x2(i) - x0(i)
        y02(i) = y2(i) - y0(i)
        z02(i) = z2(i) - z0(i)
C
        x03(i) = x3(i) - x0(i)
        y03(i) = y3(i) - y0(i)
        z03(i) = z3(i) - z0(i)
C
        x04(i) = x4(i) - x0(i)
        y04(i) = y4(i) - y0(i)
        z04(i) = z4(i) - z0(i)
C
        xi0v(i) = x0(i) - xi(i)
        yi0v(i) = y0(i) - yi(i)
        zi0v(i) = z0(i) - zi(i)
C
        xi1(i) = x1(i) - xi(i)
        yi1(i) = y1(i) - yi(i)
        zi1(i) = z1(i) - zi(i)
C
        xi2(i) = x2(i) - xi(i)
        yi2(i) = y2(i) - yi(i)
        zi2(i) = z2(i) - zi(i)
C
        xi3(i) = x3(i) - xi(i)
        yi3(i) = y3(i) - yi(i)
        zi3(i) = z3(i) - zi(i)
C
        xi4(i) = x4(i) - xi(i)
        yi4(i) = y4(i) - yi(i)
        zi4(i) = z4(i) - zi(i)
C
        sx1 = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
        sy1 = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
        sz1 = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
C
        sx2 = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
        sy2 = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
        sz2 = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
C
        sx0 = y01(i)*z02(i) - z01(i)*y02(i)
        sy0 = z01(i)*x02(i) - x01(i)*z02(i)
        sz0 = x01(i)*y02(i) - y01(i)*x02(i)
        s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
C
        nn=sqrt(s2)
        n11(i)=sx0*nn
        n21(i)=sy0*nn
        n31(i)=sz0*nn
C
        area(i)=nn
C
        lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
        lc1(i) =  (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
C
        sx3 = yi0v(i)*zi3(i) - zi0v(i)*yi3(i)
        sy3 = zi0v(i)*xi3(i) - xi0v(i)*zi3(i)
        sz3 = xi0v(i)*yi3(i) - yi0v(i)*xi3(i)
C
        sx0 = y02(i)*z03(i) - z02(i)*y03(i)
        sy0 = z02(i)*x03(i) - x02(i)*z03(i)
        sz0 = x02(i)*y03(i) - y02(i)*x03(i)
        s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
C
        nn=sqrt(s2)
        n12(i)=sx0*nn
        n22(i)=sy0*nn
        n32(i)=sz0*nn
C
        area(i)=area(i)+nn
C
        lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
        lc2(i) =  (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
C
        sx4 = yi0v(i)*zi4(i) - zi0v(i)*yi4(i)
        sy4 = zi0v(i)*xi4(i) - xi0v(i)*zi4(i)
        sz4 = xi0v(i)*yi4(i) - yi0v(i)*xi4(i)
C
        sx0 = y03(i)*z04(i) - z03(i)*y04(i)
        sy0 = z03(i)*x04(i) - x03(i)*z04(i)
        sz0 = x03(i)*y04(i) - y03(i)*x04(i)
        s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
C
        nn=sqrt(s2)
        n13(i)=sx0*nn
        n23(i)=sy0*nn
        n33(i)=sz0*nn
C
        area(i)=area(i)+nn
C
        lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
        lc3(i) =  (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
C
        sx0 = y04(i)*z01(i) - z04(i)*y01(i)
        sy0 = z04(i)*x01(i) - x04(i)*z01(i)
        sz0 = x04(i)*y01(i) - y04(i)*x01(i)
        s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
C
        nn=sqrt(s2)
        n14(i)=sx0*nn
        n24(i)=sy0*nn
        n34(i)=sz0*nn
C
        area(i)=area(i)+nn
        area(i)=half*area(i)
C
        lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
        lc4(i) =  (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
       ENDDO
C---------------------------------------------------------
C
       DO i=1,jlt
        aaa    = one/max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
        hlc1(i)= lc1(i)*abs(lc1(i))*aaa
        hlb4(i)= lb4(i)*abs(lb4(i))*aaa
        al1(i) = -(xi0v(i)*x01(i)+yi0v(i)*y01(i)+zi0v(i)*z01(i))*aaa
        al1(i) = max(zero,min(one,al1(i)))
        aaa    = one/max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
        hlc2(i)= lc2(i)*abs(lc2(i))*aaa
        hlb1(i)= lb1(i)*abs(lb1(i))*aaa
        al2(i) = -(xi0v(i)*x02(i)+yi0v(i)*y02(i)+zi0v(i)*z02(i))*aaa
        al2(i) = max(zero,min(one,al2(i)))
        aaa    = one/max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
        hlc3(i)= lc3(i)*abs(lc3(i))*aaa
        hlb2(i)= lb2(i)*abs(lb2(i))*aaa
        al3(i) = -(xi0v(i)*x03(i)+yi0v(i)*y03(i)+zi0v(i)*z03(i))*aaa
        al3(i) = max(zero,min(one,al3(i)))
        aaa    = one/max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
        hlc4(i)= lc4(i)*abs(lc4(i))*aaa
        hlb3(i)= lb3(i)*abs(lb3(i))*aaa
        al4(i) = -(xi0v(i)*x04(i)+yi0v(i)*y04(i)+zi0v(i)*z04(i))*aaa
        al4(i) = max(zero,min(one,al4(i)))
       ENDDO
C
       DO i=1,jlt
        x12 = x2(i) - x1(i)
        y12 = y2(i) - y1(i)
        z12 = z2(i) - z1(i)
        la = one - lb1(i) - lc1(i)
C     HLA, HLB1, HLC1 required for triangle sharp angle
        aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
        hla= la*abs(la) * aaa
        inproj(1,i)=0
        IF(la<zero.AND.
     +     hla<=hlb1(i).AND.hla<=hlc1(i))THEN
         lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
         lb1(i) = max(zero,min(one,lb1(i)))
         lc1(i) = one - lb1(i)
         inproj(1,i)=1
        ELSEIF(lb1(i)<zero.AND.
     +         hlb1(i)<=hlc1(i).AND.hlb1(i)<=hla)THEN
         lb1(i) = zero
         lc1(i) = al2(i)
         inproj(1,i)=1
        ELSEIF(lc1(i)<zero.AND.
     +         hlc1(i)<=hla.AND.hlc1(i)<=hlb1(i))THEN
         lc1(i) = zero
         lb1(i) = al1(i)
         inproj(1,i)=1
        ENDIF
       ENDDO
C
       DO i=1,jlt
        IF(ix3(i)==ix4(i))cycle
        x23 = x3(i) - x2(i)
        y23 = y3(i) - y2(i)
        z23 = z3(i) - z2(i)
        la = one - lb2(i) - lc2(i)
C     HLA, HLB1, HLC1 required for triangle sharp angle
        aaa = one / max(em20,x23*x23+y23*y23+z23*z23)
        hla= la*abs(la) * aaa
        inproj(2,i)=0
          IF(la<zero.AND.
     +     hla<=hlb2(i).AND.hla<=hlc2(i))THEN
            lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
            lb2(i) = max(zero,min(one,lb2(i)))
            lc2(i) = one - lb2(i)
          inproj(2,i)=1
          ELSEIF(lb2(i)<zero.AND.
     +          hlb2(i)<=hlc2(i).AND.hlb2(i)<=hla)THEN
         lb2(i) = zero
         lc2(i) = al3(i)
         inproj(2,i)=1
        ELSEIF(lc2(i)<zero.AND.
     +         hlc2(i)<=hla.AND.hlc2(i)<=hlb2(i))THEN
         lc2(i) = zero
         lb2(i) = al2(i)
         inproj(2,i)=1
        ENDIF
       ENDDO
C
       DO i=1,jlt
        IF(ix3(i)==ix4(i))cycle
        x34 = x4(i) - x3(i)
        y34 = y4(i) - y3(i)
        z34 = z4(i) - z3(i)
        la = one - lb3(i) - lc3(i)
C     HLA, HLB1, HLC1 required for triangle sharp angle
        aaa = one / max(em20,x34*x34+y34*y34+z34*z34)
        hla= la*abs(la) * aaa
        inproj(3,i)=0
        IF(la<zero.AND.
     +     hla<=hlb3(i).AND.hla<=hlc3(i))THEN
         lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
         lb3(i) = max(zero,min(one,lb3(i)))
         lc3(i) = one - lb3(i)
         inproj(3,i)=1
        ELSEIF(lb3(i)<zero.AND.
     +         hlb3(i)<=hlc3(i).AND.hlb3(i)<=hla)THEN
         lb3(i) = zero
         lc3(i) = al4(i)
         inproj(3,i)=1
        ELSEIF(lc3(i)<zero.AND.
     +         hlc3(i)<=hla.AND.hlc3(i)<=hlb3(i))THEN
         lc3(i) = zero
         lb3(i) = al3(i)
         inproj(3,i)=1
        ENDIF
       ENDDO
C
       DO i=1,jlt
        IF(ix3(i)==ix4(i))cycle
        x41 = x1(i) - x4(i)
        y41 = y1(i) - y4(i)
        z41 = z1(i) - z4(i)
        la = one - lb4(i) - lc4(i)
C     HLA, HLB1, HLC1 required for triangle sharp angle
        aaa = one / max(em20,x41*x41+y41*y41+z41*z41)
        hla= la*abs(la) * aaa
        inproj(4,i)=0
        IF(la<zero.AND.
     +     hla<=hlb4(i).AND.hla<=hlc4(i))THEN
         lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
         lb4(i) = max(zero,min(one,lb4(i)))
         lc4(i) = one - lb4(i)
         inproj(4,i)=1
        ELSEIF(lb4(i)<zero.AND.
     +         hlb4(i)<=hlc4(i).AND.hlb4(i)<=hla)THEN
         lb4(i) = zero
         lc4(i) = al1(i)
         inproj(4,i)=1
        ELSEIF(lc4(i)<zero.AND.
     +         hlc4(i)<=hla.AND.hlc4(i)<=hlb4(i))THEN
         lc4(i) = zero
         lb4(i) = al4(i)
         inproj(4,i)=1
        ENDIF
       ENDDO
C---------------------------------------------------------
       DO i=1,jlt
C
         gap2(i)=gapv(i)*gapv(i)
C
         la1 = one - lb1(i) - lc1(i)
         xp1(i) = lb1(i)*x1(i) + lc1(i)*x2(i) + la1*x0(i) 
         yp1(i) = lb1(i)*y1(i) + lc1(i)*y2(i) + la1*y0(i) 
         zp1(i) = lb1(i)*z1(i) + lc1(i)*z2(i) + la1*z0(i) 
C
         nx1(i) = xi(i)-xp1(i)
         ny1(i) = yi(i)-yp1(i)
         nz1(i) = zi(i)-zp1(i)
         p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
C
         la2 = one - lb2(i) - lc2(i)
         xp2(i) = lb2(i)*x2(i) + lc2(i)*x3(i) + la2*x0(i) 
         yp2(i) = lb2(i)*y2(i) + lc2(i)*y3(i) + la2*y0(i) 
         zp2(i) = lb2(i)*z2(i) + lc2(i)*z3(i) + la2*z0(i) 
C
         nx2(i) = xi(i)-xp2(i)
         ny2(i) = yi(i)-yp2(i)
         nz2(i) = zi(i)-zp2(i)
         p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
C
         la3 = one - lb3(i) - lc3(i)
         xp3(i) = lb3(i)*x3(i) + lc3(i)*x4(i) + la3*x0(i) 
         yp3(i) = lb3(i)*y3(i) + lc3(i)*y4(i) + la3*y0(i) 
         zp3(i) = lb3(i)*z3(i) + lc3(i)*z4(i) + la3*z0(i) 
C
         nx3(i) = xi(i)-xp3(i)
         ny3(i) = yi(i)-yp3(i)
         nz3(i) = zi(i)-zp3(i)
         p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
C
         la4 = one - lb4(i) - lc4(i)
         xp4(i) = lb4(i)*x4(i) + lc4(i)*x1(i) + la4*x0(i) 
         yp4(i) = lb4(i)*y4(i) + lc4(i)*y1(i) + la4*y0(i) 
         zp4(i) = lb4(i)*z4(i) + lc4(i)*z1(i) + la4*z0(i) 
C
         nx4(i) = xi(i)-xp4(i)
         ny4(i) = yi(i)-yp4(i)
         nz4(i) = zi(i)-zp4(i)
         p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
C
       ENDDO
C
      DO i=1,jlt
      END DO
C
      DO i=1,jlt
C
       IF(ix3(i)/=ix4(i))THEN
C
        h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
        ll=p1(i)-h*h
        IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
          p1(i)=zero
        ELSE
          p1(i)=max(zero,gapv(i)-abs(h))
        END IF
C
        h =nx2(i)*n12(i) + ny2(i)*n22(i) +nz2(i)*n32(i)
        ll=p2(i)-h*h
        IF(inproj(2,i)/=0.AND.ll >= gap2(i))THEN
          p2(i)=zero
        ELSE
          p2(i)=max(zero,gapv(i)-abs(h))
        END IF
C
        h =nx3(i)*n13(i) + ny3(i)*n23(i) +nz3(i)*n33(i)
        ll=p3(i)-h*h
        IF(inproj(3,i)/=0.AND.ll >= gap2(i))THEN
          p3(i)=zero
        ELSE
          p3(i)=max(zero,gapv(i)-abs(h))
        END IF
C
        h =nx4(i)*n14(i) + ny4(i)*n24(i) +nz4(i)*n34(i)
        ll=p4(i)-h*h
        IF(inproj(4,i)/=0.AND.ll >= gap2(i))THEN
          p4(i)=zero
        ELSE
          p4(i)=max(zero,gapv(i)-abs(h))
        END IF
C
        pene(i) = max(p1(i),p2(i),p3(i),p4(i))
C
        n1(i) = p1(i)*n11(i)+p2(i)*n12(i)+p3(i)*n13(i)+p4(i)*n14(i)
        n2(i) = p1(i)*n21(i)+p2(i)*n22(i)+p3(i)*n23(i)+p4(i)*n24(i)
        n3(i) = p1(i)*n31(i)+p2(i)*n32(i)+p3(i)*n33(i)+p4(i)*n34(i)
        nn=sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i))
        nn=one/max(em20,nn)
        n1(i)=n1(i)*nn
        n2(i)=n2(i)*nn
        n3(i)=n3(i)*nn
C
        la1 = one - lb1(i) - lc1(i)
        la2 = one - lb2(i) - lc2(i)
        la3 = one - lb3(i) - lc3(i)
        la4 = one - lb4(i) - lc4(i) 
C
        h0    = fourth * 
     .        (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
        h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
        h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
        h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
        h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
        h00    = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
        h1(i) = h1(i) * h00
        h2(i) = h2(i) * h00
        h3(i) = h3(i) * h00
        h4(i) = h4(i) * h00
C
       ELSE ! IF(IX3(I)/=IX4(I))THEN
C
        h1(i) = lb1(i)
        h2(i) = lc1(i)
        h3(i) = one - lb1(i) - lc1(i)
        h4(i) = zero
C
        n1(i) = n11(i)
        n2(i) = n21(i)
        n3(i) = n31(i)
C
        h =nx1(i)*n11(i) + ny1(i)*n21(i) +nz1(i)*n31(i)
        ll=p1(i)-h*h
        IF(inproj(1,i)/=0.AND.ll >= gap2(i))THEN
          pene(i)=zero
        ELSE
          pene(i)=max(zero,gapv(i)-abs(h))
        END IF
       END IF
      ENDDO
C
      DO 410 i=1,jlt
      xp(i)=h1(i)*x1(i)+h2(i)*x2(i)+h3(i)*x3(i)+h4(i)*x4(i)
      yp(i)=h1(i)*y1(i)+h2(i)*y2(i)+h3(i)*y3(i)+h4(i)*y4(i)
      zp(i)=h1(i)*z1(i)+h2(i)*z2(i)+h3(i)*z3(i)+h4(i)*z4(i)
 410  CONTINUE
C---------------------------------------------------------
       DO i=1,jlt
         nx1(i) = xi(i)-xp(i)
         ny1(i) = yi(i)-yp(i)
         nz1(i) = zi(i)-zp(i)
       ENDDO
C---------------------
       DO i=1,jlt
C
         nvers(1,i)    = zero
         firstimp(1,i) = zero
C
         IF(pene(i)==zero)cycle
C--------------------------------------------------------
         side=n1(i)*nx1(i)+n2(i)*ny1(i)+n3(i)*nz1(i)
C
         IF(ifpen(i)==0)THEN
            firstimp(1,i) = sign(one,side) 
            IF(firstimp(1,i) < zero)THEN
              n1(i)  = -n1(i)
              n2(i)  = -n2(i)
              n3(i)  = -n3(i)
            END IF
            nvers(1,i) = one
C 1st impact below gap (sorting security)
         ELSE ! IF(IFPEN(I)==0.OR.TT==ZERO)THEN
            IF(ifpen(i) < 0)THEN
              n1(i)  = -n1(i)
              n2(i)  = -n2(i)
              n3(i)  = -n3(i)
C             SIDE   = -SIDE
            END IF
            nvers(1,i)= sign(one,side*ifpen(i))
         END IF
C---------------------
         IF(nvers(1,i)==zero)THEN
          pene(i)=zero
          cycle
         END IF
C
C beware of crossing the shell => 1E20...,
C Take normal = normal at triangle
         nn=one/
     .      max(em20,sqrt(n1(i)*n1(i)+n2(i)*n2(i)+n3(i)*n3(i)))
         n1(i)=n1(i)*nn
         n2(i)=n2(i)*nn
         n3(i)=n3(i)*nn
C         
         nx(i)=n1(i)
         ny(i)=n2(i)
         nz(i)=n3(i)
C         
       END DO
C---------------------
       DO i=1,jlt
         IF(pene(i)/=zero)THEN
           IF(firstimp(1,i) > zero)THEN
             ifpen(i) = 1
           ELSEIF(firstimp(1,i) < zero)THEN
             ifpen(i) =-1
           END IF
         END IF
       ENDDO
C---------------------
      RETURN
      END