i7keg3.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/.
!||====================================================================
!||    i7keg3     ../engine/source/interfaces/int07/i7keg3.F
!||--- called by ------------------------------------------------------
!||    i7ke3      ../engine/source/interfaces/int07/i7ke3.F
!||--- uses       -----------------------------------------------------
!||    imp_intm   ../engine/share/modules/imp_intm.F
!||====================================================================
      SUBROUTINE i7keg3(JLT    ,A      ,V      ,MS    ,FRIC   ,
     1                  NX1    ,NX2    ,NX3    ,NX4   ,NY1    ,
     2                  NY2    ,NY3    ,NY4    ,NZ1   ,NZ2    ,
     3                  NZ3    ,NZ4    ,LB1    ,LB2   ,LB3    ,
     4                  LB4    ,LC1    ,LC2    ,LC3   ,LC4    ,
     5                  P1     ,P2     ,P3     ,P4    ,NIN    ,
     6                  IX1    ,IX2    ,IX3    ,IX4   ,NSVG   ,
     7                  GAPV   ,INACTI ,CAND_P ,INDEX ,STIGLO ,
     8                  STIF   ,VXI    ,VYI    ,VZI   ,MSI    ,
     9                  KI11   ,KI12   ,KJ11  ,KJ12   ,KK11   ,
     A                  KK12   ,KL11   ,KL12  ,OFF    ,SCALK  ,
     B                  LREM   ,ICURV  ,X      ,
     2                  X1     ,X2     ,X3     ,X4     ,Y1     ,
     3                  Y2     ,Y3     ,Y4     ,Z1     ,Z2     ,
     4                  Z3     ,Z4     ,XI     ,YI     ,ZI     )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE imp_intm
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      "scr05_c.inc"
Ctmp+1
#include      "com01_c.inc"
#include      "impl1_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JLT, INACTI,NIN,LREM
      INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
     .        NSVG(MVSIZ), INDEX(MVSIZ),ICURV(3)
      my_real
     .   A(3,*), MS(*), V(3,*),
     .   STIGLO,CAND_P(*),FRIC,OFF(*),SCALK,
     .   VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ)
      my_real
     .     NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
     .     NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
     .     NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
     .     LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
     .     LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
     .    P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ), STIF(MVSIZ),
     .    gapv(mvsiz),ki11(3,3,mvsiz),kj11(3,3,mvsiz),
     .    kk11(3,3,mvsiz),kl11(3,3,mvsiz),ki12(3,3,mvsiz),
     .    kj12(3,3,mvsiz),kk12(3,3,mvsiz),kl12(3,3,mvsiz)
      my_real
     .     x1(mvsiz),y1(mvsiz),z1(mvsiz),
     .     x2(mvsiz),y2(mvsiz),z2(mvsiz),
     .     x3(mvsiz),y3(mvsiz),z3(mvsiz),
     .     x4(mvsiz),y4(mvsiz),z4(mvsiz),
     .     xi(mvsiz),yi(mvsiz),zi(mvsiz),x(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J1, J, K,IG,ISF,NN,NS,JLTF,NE
      my_real
     .   N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
     .   H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
     .   VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ), 
     .   s2,fac,facf, h0, la1, la2, la3, la4,fact(mvsiz),
     .   d1,d2,d3,d4,a1,a2,a3,a4,kn(4,mvsiz),q(3,3,mvsiz)
      my_real
     .   prec,q11,q12,q13,q22,q23,q33,h00,vtx,vty,vtz,vt,
     .   kt1,kt2,kt3,kt4,q1,q2
      INTEGER NA1,NA2
      my_real
     .   A0X,A0Y,A0Z,RX,RY,RZ,
     .   ANX,ANY,ANZ,AAN,AAX,AAY,AAZ ,RR,RS,AAA ,TM,TS
C-----------------------------------------------
C  dev021
      IF (iresp==1) THEN
           prec = fiveem4
      ELSE
           prec = em10
      ENDIF
C--------------------------------------------------------
C  CAS DES PAQUETS MIXTES
C--------------------------------------------------------
      IF (imp_int7==3) THEN
       DO i=1,jlt
        d1 = sqrt(p1(i))
        p1(i) = fourth*gapv(i) 
        d2 = sqrt(p2(i))
        p2(i) = fourth*gapv(i)
        d3 = sqrt(p3(i))
        p3(i) = fourth*gapv(i)
        d4 = sqrt(p4(i))
        p4(i) = fourth*gapv(i)
        a1 = p1(i)/max(em20,d1)
        a2 = p2(i)/max(em20,d2)
        a3 = p3(i)/max(em20,d3)
        a4 = p4(i)/max(em20,d4)
        n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i) 
        n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i) 
        n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i) 
       ENDDO
      ELSE
       DO i=1,jlt
C
        d1 = sqrt(p1(i))
        p1(i) = max(zero, gapv(i) - d1)
C
        d2 = sqrt(p2(i))
        p2(i) = max(zero, gapv(i) - d2)
C
        d3 = sqrt(p3(i))
        p3(i) = max(zero, gapv(i) - d3)
C
        d4 = sqrt(p4(i))
        p4(i) = max(zero, gapv(i) - d4)
C
        a1 = p1(i)/max(em20,d1)
        a2 = p2(i)/max(em20,d2)
        a3 = p3(i)/max(em20,d3)
        a4 = p4(i)/max(em20,d4)
        n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i) 
        n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i) 
        n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i) 
       ENDDO
      ENDIF !(IMP_INT7==3) 
C
       DO i=1,jlt
        IF(ix3(i)/=ix4(i))THEN
         pene(i) = max(p1(i),p2(i),p3(i),p4(i))
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
         pene(i) = p1(i)
         n1(i) = nx1(i)
         n2(i) = ny1(i)
         n3(i) = nz1(i)
         h1(i) = lb1(i)
         h2(i) = lc1(i)
         h3(i) = one - lb1(i) - lc1(i)
         h4(i) = zero
        ENDIF
       ENDDO
C---------------------
C     COURBURE FIXE
C---------------------
      IF(icurv(1)==1)THEN
C       spherique (que concave pour le moment)
        na1 = icurv(2)
        DO i=1,jlt
          rr = 1.e30
          a0x = x(1,na1)
          a0y = x(2,na1)
          a0z = x(3,na1)
C
          rx = x1(i)-a0x
          ry = y1(i)-a0y
          rz = z1(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          rx = x2(i)-a0x
          ry = y2(i)-a0y
          rz = z2(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          rx = x3(i)-a0x
          ry = y3(i)-a0y
          rz = z3(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          IF(ix3(i)/=ix4(i))THEN
            rx = x4(i)-a0x
            ry = y4(i)-a0y
            rz = z4(i)-a0z
            rr = min(rr , rx*rx + ry*ry + rz*rz)
          ENDIF
          rx = xi(i)-a0x
          ry = yi(i)-a0y
          rz = zi(i)-a0z
          rs = sqrt(rx*rx + ry*ry + rz*rz)
          rr = sqrt(rr)
          IF(rs-rr+gapv(i)<0.)THEN
            stif(i) = 0.
            pene(i) = 0.
          ELSEIF(rs-rr+gapv(i)<pene(i))THEN
            pene(i) = rs-rr+gapv(i)
          ENDIF
          n1(i) = -rx 
          n2(i) = -ry 
          n3(i) = -rz 
        ENDDO
      ELSEIF(icurv(1)==2)THEN
C       cylindrique (que concave pour le moment)
        na1 = icurv(2)
        na2 = icurv(3)
        DO i=1,jlt
          rr = 1.e30
          a0x = x(1,na1)
          a0y = x(2,na1)
          a0z = x(3,na1)
          anx = x(1,na2)-a0x
          any = x(2,na2)-a0y
          anz = x(3,na2)-a0z
          aan = 1. / (anx*anx + any*any + anz*anz) 
C
          aax = x1(i)-a0x
          aay = y1(i)-a0y
          aaz = z1(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          aax = x2(i)-a0x
          aay = y2(i)-a0y
          aaz = z2(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          aax = x3(i)-a0x
          aay = y3(i)-a0y
          aaz = z3(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          IF(ix3(i)/=ix4(i))THEN
            aax = x4(i)-a0x
            aay = y4(i)-a0y
            aaz = z4(i)-a0z
            aaa = (aax*anx + aay*any + aaz*anz) * aan
            rx = aax - aaa * anx 
            ry = aay - aaa * any 
            rz = aaz - aaa * anz 
            rr = min(rr , rx*rx + ry*ry + rz*rz)
          ENDIF
          aax = xi(i)-a0x
          aay = yi(i)-a0y
          aaz = zi(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rs = sqrt(rx*rx + ry*ry + rz*rz)
          rr = sqrt(rr)
          IF(rs-rr+gapv(i)<0.)THEN
            stif(i) = 0.
            pene(i) = 0.
          ELSEIF(rs-rr+gapv(i)<pene(i))THEN
            pene(i) = rs-rr+gapv(i)
            n1(i) = -rx 
            n2(i) = -ry 
            n3(i) = -rz 
          ENDIF
        ENDDO
      ELSEIF(icurv(1) == 3)THEN
c        print *,'ICURV=3 NOT AVAILABLE WITH IMPLICIT OPTION'
      ENDIF
C
      DO i=1,jlt
         s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
         n1(i) = n1(i)*s2
         n2(i) = n2(i)*s2
         n3(i) = n3(i)*s2
      ENDDO
C 
      DO 500 i=1,jlt
        vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
     .                 - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
        vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
     .                 - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
        vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
     .                 - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
        vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
C       correction hourglass
         h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
         h0 = min(h0,h2(i),h4(i))
         h0 = max(h0,-h1(i),-h3(i))
         IF(ix3(i)==ix4(i))h0 = zero
         h1(i) = h1(i) + h0
         h2(i) = h2(i) - h0
         h3(i) = h3(i) + h0
         h4(i) = h4(i) - h0
 500  CONTINUE
C---------------------
C      PENE INITIALE
C---------------------
      IF(inacti==5)THEN
         DO i=1,jlt
           pene(i)=max(zero,pene(i)-cand_p(index(i)))
           IF( pene(i)==zero )  stif(i) = zero
           gapv(i)=gapv(i)-cand_p(index(i))
         ENDDO
      ELSE IF(inacti==6)THEN
         DO i=1,jlt
           cand_p(index(i))=min(cand_p(index(i)),
     .        ( (one-fiveem2)*cand_p(index(i))
     .          +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i))))  )
           pene(i)=max(zero,pene(i)-cand_p(index(i)))
           IF( pene(i)==zero )  stif(i) = zero
           gapv(i)=gapv(i)-cand_p(index(i))
         ENDDO
      ENDIF
C-------------------------------------------
c      print *,'IMP_INT7=',IMP_INT7,jlt
      IF(imp_int7>=2)THEN
      DO i=1,jlt
       IF(stiglo<=zero)THEN
        stif(i) = half*stif(i) 
       ELSEIF(stif(i)/=zero)THEN
        stif(i) = stiglo 
       ENDIF
      ENDDO
      ELSEIF(imp_int7==1)THEN
      DO i=1,jlt
       fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
       IF( (gapv(i)-pene(i))/gapv(i) <prec .AND. 
     .     stif(i)>0. ) stif(i) = 0.
       IF(stiglo<=zero)THEN
        stif(i) = half*stif(i) * fac
       ELSEIF(stif(i)/=zero)THEN
        stif(i) = stiglo * fac
       ENDIF
      ENDDO
      ELSE
C
      DO i=1,jlt
       fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
       IF( (gapv(i)-pene(i))/gapv(i) <prec .AND. 
     .     stif(i)>0. ) stif(i) = 0.
       IF(stiglo<=zero)THEN
        stif(i) = half*stif(i) * fac
       ELSEIF(stif(i)/=zero)THEN
        stif(i) = stiglo * fac
       ENDIF
      ENDDO
      DO i=1,jlt
          stif(i) = stif(i) * gapv(i) / 
     .            max((gapv(i) - pene(i)),em10) 
      ENDDO
      ENDIF
C---------------------------------
C    ----sans frottement d'abord--- 
      DO i=1,jlt
        vtx = vx(i) -vn(i)*n1(i)
        vty = vy(i) -vn(i)*n2(i)
        vtz = vz(i) -vn(i)*n3(i)
        vt  = vtx*vtx+vty*vty+vtz*vtz
        IF (vt>em20) THEN
         s2=one/sqrt(vt)
         q(1,1,i)=vtx*s2
         q(1,2,i)=vty*s2
         q(1,3,i)=vtz*s2
         q(3,1,i)=n1(i)
         q(3,2,i)=n2(i)
         q(3,3,i)=n3(i)
         q(2,1,i)=q(3,2,i)*q(1,3,i)-q(3,3,i)*q(1,2,i)
         q(2,2,i)=q(3,3,i)*q(1,1,i)-q(3,1,i)*q(1,3,i)
         q(2,3,i)=q(3,1,i)*q(1,2,i)-q(3,2,i)*q(1,1,i)
         fact(i)=fric
        ELSE
         fact(i)=zero
        ENDIF
      ENDDO
      IF (scalk<0) THEN
       isf=1
      ELSE
       isf=0
      ENDIF
      facf=abs(scalk)
      IF (isf==1) THEN
       DO i=1,jlt
        IF (vn(i)>zero) THEN
         fac=stif(i)/facf
        ELSEIF (vn(i)<zero) THEN
         fac=stif(i)*facf
        ELSE
         fac=stif(i)
        ENDIF
        kn(1,i)=fac*h1(i)
        kn(2,i)=fac*h2(i)
        kn(3,i)=fac*h3(i)
        kn(4,i)=fac*h4(i)
        fact(i)=fac*fact(i)
       ENDDO
      ELSE
       DO i=1,jlt
        fac=stif(i)*facf
        kn(1,i)=fac*h1(i)
        kn(2,i)=fac*h2(i)
        kn(3,i)=fac*h3(i)
        kn(4,i)=fac*h4(i)
        fact(i)=fac*fact(i)
       ENDDO
      ENDIF
      DO i=1,jlt
       q11=n1(i)*n1(i)
       q12=n1(i)*n2(i)
       q13=n1(i)*n3(i)
       q22=n2(i)*n2(i)
       q23=n2(i)*n3(i)
       q33=n3(i)*n3(i)
       ki11(1,1,i)=kn(1,i)*q11
       ki11(1,2,i)=kn(1,i)*q12
       ki11(1,3,i)=kn(1,i)*q13
       ki11(2,2,i)=kn(1,i)*q22
       ki11(2,3,i)=kn(1,i)*q23
       ki11(3,3,i)=kn(1,i)*q33
       kj11(1,1,i)=kn(2,i)*q11
       kj11(1,2,i)=kn(2,i)*q12
       kj11(1,3,i)=kn(2,i)*q13
       kj11(2,2,i)=kn(2,i)*q22
       kj11(2,3,i)=kn(2,i)*q23
       kj11(3,3,i)=kn(2,i)*q33
       kk11(1,1,i)=kn(3,i)*q11
       kk11(1,2,i)=kn(3,i)*q12
       kk11(1,3,i)=kn(3,i)*q13
       kk11(2,2,i)=kn(3,i)*q22
       kk11(2,3,i)=kn(3,i)*q23
       kk11(3,3,i)=kn(3,i)*q33
       kl11(1,1,i)=kn(4,i)*q11
       kl11(1,2,i)=kn(4,i)*q12
       kl11(1,3,i)=kn(4,i)*q13
       kl11(2,2,i)=kn(4,i)*q22
       kl11(2,3,i)=kn(4,i)*q23
       kl11(3,3,i)=kn(4,i)*q33
      ENDDO
C    ----avec frottement --- 
       DO j=1,3 
        DO k=j,3 
         DO i=1,jlt
          IF (fact(i)>zero) THEN
           q1 =q(1,j,i)*q(1,k,i)
           q2 =q(2,j,i)*q(2,k,i)
           fac=fact(i)*(q1+q2)
           kt1=fac*h1(i)
           ki11(j,k,i)=ki11(j,k,i)+kt1
           kt2=fac*h2(i)
           kj11(j,k,i)=kj11(j,k,i)+kt2
           kt3=fac*h3(i)
           kk11(j,k,i)=kk11(j,k,i)+kt3
           kt4=fac*h4(i)
           kl11(j,k,i)=kl11(j,k,i)+kt4
          ENDIF 
         ENDDO
        ENDDO
       ENDDO
C
       DO j=1,3 
        DO k=j,3 
         DO i=1,jlt
          ki12(j,k,i)=-ki11(j,k,i)
          kj12(j,k,i)=-kj11(j,k,i)
          kk12(j,k,i)=-kk11(j,k,i)
          kl12(j,k,i)=-kl11(j,k,i)
         ENDDO
        ENDDO
       ENDDO
       DO j=1,3 
        DO k=j+1,3 
         DO i=1,jlt
          ki12(k,j,i)=-ki11(j,k,i)
          kj12(k,j,i)=-kj11(j,k,i)
          kk12(k,j,i)=-kk11(j,k,i)
          kl12(k,j,i)=-kl11(j,k,i)
         ENDDO
        ENDDO
       ENDDO
C
       DO i=1,jlt
        off(i)=one
       ENDDO
C
       IF (intp_d <= 0 .AND. nspmd > 1) THEN
        jltf = 0
        DO i=1,jlt
         IF(nsvg(i)<0) THEN
            nn=-nsvg(i)
          jltf = jltf + 1
            ne=shf_int(nin) + jltf +lrem
            ns=ind_int(nin)%P(nn)
            stifs(ne)=stif(i)
            h_e(1,ne)=h1(i)
            h_e(2,ne)=h2(i)
            h_e(3,ne)=h3(i)
            h_e(4,ne)=h4(i)
            n_e(1,ne)=n1(i)
            n_e(2,ne)=n2(i)
            n_e(3,ne)=n3(i)
           ENDIF
        ENDDO
       ENDIF 
C
       RETURN
      END
!||====================================================================
!||    i7frf3     ../engine/source/interfaces/int07/i7keg3.F
!||--- called by ------------------------------------------------------
!||    i7forcf3   ../engine/source/interfaces/int07/i7ke3.F
!||--- uses       -----------------------------------------------------
!||    imp_intm   ../engine/share/modules/imp_intm.F
!||====================================================================
      SUBROUTINE i7frf3(JLT    ,A      ,V      ,MS     ,FRIC   ,
     1                  N1     ,N2     ,N3     ,H1     ,H2     ,
     2                  H3     ,H4     ,IX1    ,IX2    ,IX3    ,
     2                  IX4    ,INDEX  ,VXI    ,VYI    ,VZI    ,   
     3                  MSI    ,DXI    ,DYI    ,DZI    ,STIF   ,
     4                  NIN    ,D      ,SCALK  )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE imp_intm
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, INACTI,NIN
      INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
     .        INDEX(MVSIZ)
      my_real
     .   A(3,*), MS(*), V(3,*),D(3,*),
     .   FRIC,SCALK,DXI(MVSIZ),DYI(MVSIZ),DZI(MVSIZ),
     .   H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
     .   vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
      my_real
     .     n1(mvsiz), n2(mvsiz), n3(mvsiz),stif(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J1, J, K,IG,ISF,NN,NS,NI
      my_real
     .   VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ), 
     .   DX(MVSIZ), DY(MVSIZ), DZ(MVSIZ), DN(MVSIZ), 
     .   DNI(MVSIZ),DT(MVSIZ), DTI(MVSIZ), 
     .   S2,FACN(MVSIZ),FACF, FACT(MVSIZ)
      my_real
     .   fx,fy,fz,fn,ft,fni,fti,vtx,vty,vtz,vt,
     .   t1(mvsiz), t2(mvsiz), t3(mvsiz),
     .   kt1,kt2,kt3,kt4,q1,q2,h0
 
C-----------------------------------------------
C
      DO i=1,jlt
        vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
     .                 - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
        vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
     .                 - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
        vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
     .                 - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
        vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
        dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
     .                 - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
        dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
     .                 - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
        dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
     .                 - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
        dn(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
        dni(i) = n1(i)*dxi(i) + n2(i)*dyi(i) + n3(i)*dzi(i)
      ENDDO
C-------------------------------------------
      DO i=1,jlt
        vtx = vx(i) -vn(i)*n1(i)
        vty = vy(i) -vn(i)*n2(i)
        vtz = vz(i) -vn(i)*n3(i)
        vt  = vtx*vtx+vty*vty+vtz*vtz
        IF (vt>em20) THEN
         s2=one/sqrt(vt)
         t1(i)=vtx*s2
         t2(i)=vty*s2
         t3(i)=vtz*s2
         fact(i)=fric
        ELSE
         fact(i)=zero
         t1(i)=one
         t2(i)=zero
         t3(i)=zero
        ENDIF
      ENDDO
      DO i=1,jlt
       dt(i) = t1(i)*dx(i) + t2(i)*dy(i) + t3(i)*dz(i)
       dti(i) = t1(i)*dxi(i) + t2(i)*dyi(i) + t3(i)*dzi(i)
      ENDDO
      IF (scalk<0) THEN
       isf=1
      ELSE
       isf=0
      ENDIF
      facf=abs(scalk)
      IF (isf==1) THEN
       DO i=1,jlt
        IF (vn(i)>zero) THEN
         facn(i)=stif(i)*facf
        ELSEIF (vn(i)<zero) THEN
         facn(i)=stif(i)/facf
        ELSE
         facn(i)=stif(i)
        ENDIF
        fact(i)=facn(i)*fact(i)
       ENDDO
      ELSE
       DO i=1,jlt
        facn(i)=stif(i)*facf
        fact(i)=facn(i)*fact(i)
       ENDDO
      ENDIF
C--------partie NML-------
      DO i=1,jlt
        fn = -facn(i)*dni(i)
        fx=fn*n1(i)
        fy=fn*n2(i)
        fz=fn*n3(i)
        IF (fact(i)/=zero) THEN
         ft = -fact(i)*dti(i)
         fx = fx + ft*t1(i)
         fy = fy + ft*t2(i)
         fz = fz + ft*t3(i)
        ENDIF
        a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
        a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
        a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
        a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
        a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
        a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
        a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
        a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
        a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
        a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
        a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
        a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
      ENDDO
C--------partie NSL-------
      DO i=1,jlt
        fni = facn(i)*dn(i)
        ni = index(i)
        ffi(1,ni)=ffi(1,ni)+fni*n1(i)
        ffi(2,ni)=ffi(2,ni)+fni*n2(i)
        ffi(3,ni)=ffi(3,ni)+fni*n3(i)
        IF (fact(i)/=zero) THEN
         fti = fact(i)*dt(i)
         ffi(1,ni)=ffi(1,ni)+fti*t1(i)
         ffi(2,ni)=ffi(2,ni)+fti*t2(i)
         ffi(3,ni)=ffi(3,ni)+fti*t3(i)
        ENDIF
      ENDDO
C
       RETURN
      END
!||====================================================================
!||    i7kfor3    ../engine/source/interfaces/int07/i7keg3.F
!||--- called by ------------------------------------------------------
!||    i7fku3     ../engine/source/interfaces/int07/i7ke3.F
!||--- uses       -----------------------------------------------------
!||    imp_intm   ../engine/share/modules/imp_intm.F
!||====================================================================
      SUBROUTINE i7kfor3(JLT    ,A      ,V      ,MS    ,FRIC   ,
     1                  NX1    ,NX2    ,NX3    ,NX4   ,NY1    ,
     2                  NY2    ,NY3    ,NY4    ,NZ1   ,NZ2    ,
     3                  NZ3    ,NZ4    ,LB1    ,LB2   ,LB3    ,
     4                  LB4    ,LC1    ,LC2    ,LC3   ,LC4    ,
     5                  P1     ,P2     ,P3     ,P4    ,NIN    ,
     6                  IX1    ,IX2    ,IX3    ,IX4   ,NSVG   ,
     7                  GAPV   ,INACTI ,CAND_P ,INDEX ,STIGLO ,
     8                  STIF   ,VXI    ,VYI    ,VZI   ,MSI    ,
     F                  X      ,IRECT  ,CE_LOC ,MFROT ,IFQ    ,
     G                  FROT_P ,CAND_FX,CAND_FY,CAND_FZ,ALPHA0,
     H                  DXI    ,DYI    ,DZI    ,D     ,SCALK  ,
     2                  X1     ,X2     ,X3     ,X4     ,Y1     ,
     3                  Y2     ,Y3     ,Y4     ,Z1     ,Z2     ,
     4                  Z3     ,Z4     ,XI     ,YI     ,ZI     ,
     5                  ICURV  )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE imp_intm
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.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      "com08_c.inc"
#include      "scr05_c.inc"
#include      "impl1_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JLT,INACTI,NIN,MFROT, IFQ, NOINT,IRECT(4,*)
      INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
     .        CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ)
      my_real
     .   STIGLO,CAND_P(*),FROT_P(*), X(3,*),A(3,*), MS(*), V(3,*), 
     .   CAND_FX(*),CAND_FY(*),CAND_FZ(*),ALPHA0,FRIC,SCALK,ICURV(3)
      my_real
     .     NX1(MVSIZ), NX2(MVSIZ), NX3(MVSIZ), NX4(MVSIZ),
     .     NY1(MVSIZ), NY2(MVSIZ), NY3(MVSIZ), NY4(MVSIZ),
     .     NZ1(MVSIZ), NZ2(MVSIZ), NZ3(MVSIZ), NZ4(MVSIZ),
     .     LB1(MVSIZ), LB2(MVSIZ), LB3(MVSIZ), LB4(MVSIZ),
     .     LC1(MVSIZ), LC2(MVSIZ), LC3(MVSIZ), LC4(MVSIZ),
     .     P1(MVSIZ), P2(MVSIZ), P3(MVSIZ), P4(MVSIZ), STIF(MVSIZ),
     .     GAPV(MVSIZ),VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ),
     .     DXI(MVSIZ),DYI(MVSIZ),DZI(MVSIZ),D(3,*)
      my_real
     .     x1(mvsiz),y1(mvsiz),z1(mvsiz),
     .     x2(mvsiz),y2(mvsiz),z2(mvsiz),
     .     x3(mvsiz),y3(mvsiz),z3(mvsiz),
     .     x4(mvsiz),y4(mvsiz),z4(mvsiz),
     .     xi(mvsiz),yi(mvsiz),zi(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,NS
      my_real
     .   FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ), FNI(MVSIZ),
     .   FXT(MVSIZ),FYT(MVSIZ),FZT(MVSIZ),STIF0(MVSIZ),
     .   N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PENE(MVSIZ),
     .   H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
     .   VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ),XMU(MVSIZ), 
     .   DX(MVSIZ), DY(MVSIZ), DZ(MVSIZ), DN0(MVSIZ),
     .    
     .   VNX, VNY, VNZ, AA, CRIT,S2,DIST,RDIST,
     .   v2, fm2, fac,ff,alphi,alpha,beta,
     .   fx, fy, fz, f2, mas2, m2sk, dti,ft,fn,fmax,ftn,
     .   h0, la1, la2, la3, la4,d1,d2,d3,d4,a1,a2,a3,a4, ff0,
     .   vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,t1,t2,t3,dxt,
     .   area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
     .   anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa ,gap2 ,pene2,prec 
      INTEGER NA1,NA2
C-----------------------------------------------
      IF (IRESP==1) then
           prec = fiveem4
      ELSE
           prec = em10
      ENDIF
C---------------------
C      PENE INITIALE
C---------------------
       IF(inacti==5.OR.inacti==6)THEN
        DO i=1,jlt
          IF(stif(i)==zero)THEN
            cand_p(index(i))=0
          ENDIF
        ENDDO
       ENDIF
C--------------------------------------------------------
C  actualise stif
C--------------------------------------------------------
       DO i=1,jlt
        gap2=gapv(i)*gapv(i)
C
        d1 = max(zero, gap2 - p1(i))
        d2 = max(zero, gap2 - p2(i))
        d3 = max(zero, gap2 - p3(i))
        d4 = max(zero, gap2 - p4(i))
        pene2 = max(d1,d2,d3,d4)
        IF (pene2<=zero) stif(i) = zero
       ENDDO
C--------------------------------------------------------
C  CAS DES PAQUETS MIXTES
C--------------------------------------------------------
       DO i=1,jlt
        d1 = sqrt(p1(i))
        p1(i) = max(zero, gapv(i) - d1)
C
        d2 = sqrt(p2(i))
        p2(i) = max(zero, gapv(i) - d2)
C
        d3 = sqrt(p3(i))
        p3(i) = max(zero, gapv(i) - d3)
C
        d4 = sqrt(p4(i))
        p4(i) = max(zero, gapv(i) - d4)
C
        a1 = p1(i)/max(em20,d1)
        a2 = p2(i)/max(em20,d2)
        a3 = p3(i)/max(em20,d3)
        a4 = p4(i)/max(em20,d4)
        n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i) 
        n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i) 
        n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i) 
       ENDDO
C
       DO i=1,jlt
        IF(ix3(i)/=ix4(i))THEN
         pene(i) = max(p1(i),p2(i),p3(i),p4(i))
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
         pene(i) = p1(i)
         n1(i) = nx1(i)
         n2(i) = ny1(i)
         n3(i) = nz1(i)
         h1(i) = lb1(i)
         h2(i) = lc1(i)
         h3(i) = one - lb1(i) - lc1(i)
         h4(i) = zero
        ENDIF
       ENDDO
C---------------------
C     COURBURE FIXE
C---------------------
      IF(icurv(1)==1)THEN
C       spherique (que concave pour le moment)
        na1 = icurv(2)
        DO i=1,jlt
          rr = 1.e30
          a0x = x(1,na1)
          a0y = x(2,na1)
          a0z = x(3,na1)
C
          rx = x1(i)-a0x
          ry = y1(i)-a0y
          rz = z1(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          rx = x2(i)-a0x
          ry = y2(i)-a0y
          rz = z2(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          rx = x3(i)-a0x
          ry = y3(i)-a0y
          rz = z3(i)-a0z
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          IF(ix3(i)/=ix4(i))THEN
            rx = x4(i)-a0x
            ry = y4(i)-a0y
            rz = z4(i)-a0z
            rr = min(rr , rx*rx + ry*ry + rz*rz)
          ENDIF
          rx = xi(i)-a0x
          ry = yi(i)-a0y
          rz = zi(i)-a0z
          rs = sqrt(rx*rx + ry*ry + rz*rz)
          rr = sqrt(rr)
          IF(rs-rr+gapv(i)<0.)THEN
            stif(i) = 0.
            pene(i) = 0.
          ELSEIF(rs-rr+gapv(i)<pene(i))THEN
            pene(i) = rs-rr+gapv(i)
          ENDIF
          n1(i) = -rx 
          n2(i) = -ry 
          n3(i) = -rz 
        ENDDO
      ELSEIF(icurv(1)==2)THEN
C       cylindrique (que concave pour le moment)
        na1 = icurv(2)
        na2 = icurv(3)
        DO i=1,jlt
          rr = 1.e30
          a0x = x(1,na1)
          a0y = x(2,na1)
          a0z = x(3,na1)
          anx = x(1,na2)-a0x
          any = x(2,na2)-a0y
          anz = x(3,na2)-a0z
          aan = 1. / (anx*anx + any*any + anz*anz) 
C
          aax = x1(i)-a0x
          aay = y1(i)-a0y
          aaz = z1(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          aax = x2(i)-a0x
          aay = y2(i)-a0y
          aaz = z2(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          aax = x3(i)-a0x
          aay = y3(i)-a0y
          aaz = z3(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rr = min(rr , rx*rx + ry*ry + rz*rz)
          IF(ix3(i)/=ix4(i))THEN
            aax = x4(i)-a0x
            aay = y4(i)-a0y
            aaz = z4(i)-a0z
            aaa = (aax*anx + aay*any + aaz*anz) * aan
            rx = aax - aaa * anx 
            ry = aay - aaa * any 
            rz = aaz - aaa * anz 
            rr = min(rr , rx*rx + ry*ry + rz*rz)
          ENDIF
          aax = xi(i)-a0x
          aay = yi(i)-a0y
          aaz = zi(i)-a0z
          aaa = (aax*anx + aay*any + aaz*anz) * aan
          rx = aax - aaa * anx 
          ry = aay - aaa * any 
          rz = aaz - aaa * anz 
          rs = sqrt(rx*rx + ry*ry + rz*rz)
          rr = sqrt(rr)
          IF(rs-rr+gapv(i)<0.)THEN
            stif(i) = 0.
            pene(i) = 0.
          ELSEIF(rs-rr+gapv(i)<pene(i))THEN
            pene(i) = rs-rr+gapv(i)
            n1(i) = -rx 
            n2(i) = -ry 
            n3(i) = -rz 
          ENDIF
        ENDDO
      ELSEIF(icurv(1) == 3)THEN
c        print *,'ICURV=3 NOT AVAILABLE WITH IMPLICIT OPTION'
      ENDIF
C
       DO i=1,jlt
         s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
         n1(i) = n1(i)*s2
         n2(i) = n2(i)*s2
         n3(i) = n3(i)*s2
       ENDDO
C
       DO i=1,jlt
        dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
     .                 - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
        dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
     .                 - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
        dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
     .                 - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
        dn0(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
       ENDDO
C
      DO 500 i=1,jlt
        vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
     .                 - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
        vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
     .                 - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
        vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
     .                 - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
        vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
C       correction hourglass
         h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
         h0 = min(h0,h2(i),h4(i))
         h0 = max(h0,-h1(i),-h3(i))
         IF(ix3(i)==ix4(i))h0 = zero
         h1(i) = h1(i) + h0
         h2(i) = h2(i) - h0
         h3(i) = h3(i) + h0
         h4(i) = h4(i) - h0
 500  CONTINUE
C---------------------
C      PENE INITIALE
C---------------------
      IF(inacti==5)THEN
         DO i=1,jlt
C REDUCTION DE LA PENE INITIALE
           cand_p(index(i))=min(cand_p(index(i)),
     .            ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i))  )
C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU GAP
           pene(i)=max(zero,pene(i)-cand_p(index(i)))
           IF( pene(i)==zero )  stif(i) = zero
           gapv(i)=gapv(i)-cand_p(index(i))
         ENDDO
      ELSE IF(inacti==6)THEN
         DO i=1,jlt
C REDUCTION DE LA PENE INITIALE
C           CAND_P(INDEX(I))=MIN(CAND_P(INDEX(I)),PENE(I))
           cand_p(index(i))=min(cand_p(index(i)),
     .        ( (one-fiveem2)*cand_p(index(i))
     .          +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i))))  )
C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU GAP
           pene(i)=max(zero,pene(i)-cand_p(index(i)))
           IF( pene(i)==zero )  stif(i) = zero
           gapv(i)=gapv(i)-cand_p(index(i))
         ENDDO
      ENDIF
C       
       DO i=1,jlt
        stif0(i) = stif(i)
       ENDDO
C-------------------------------------------
      IF(imp_int7>=2)THEN
       DO i=1,jlt
        IF(stiglo<=zero)THEN
         stif(i) = half*stif(i) 
        ELSEIF(stif(i)/=zero)THEN
         stif(i) = stiglo 
        ENDIF
       ENDDO
      ELSEIF(imp_int7==1)THEN
       DO i=1,jlt
        fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
        IF( (gapv(i)-pene(i))/gapv(i) <prec .AND. 
     .     stif(i)>0. ) stif(i) = 0.
        IF(stiglo<=zero)THEN
         stif(i) = half*stif(i) * fac
        ELSEIF(stif(i)/=zero)THEN
         stif(i) = stiglo * fac
        ENDIF
       ENDDO
      ELSE
       DO i=1,jlt
        fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
        IF( (gapv(i)-pene(i))/gapv(i) <prec .AND. 
     .     stif(i)>0. ) stif(i) = 0.
        IF(stiglo<=zero)THEN
         stif(i) = half*stif(i) * fac
        ELSEIF(stif(i)/=zero)THEN
         stif(i) = stiglo * fac
        ENDIF
       ENDDO
       DO i=1,jlt
          stif(i) = stif(i) * gapv(i) / 
     .            max((gapv(i) - pene(i)),em10) 
       ENDDO
      ENDIF ! (IMP_INT7>=2)
C 
       fac =   abs(scalk)   
       DO i=1,jlt
        stif(i)=stif(i)*fac
        fni(i)= -stif(i) * dn0(i) 
        fxi(i)=n1(i)*fni(i)
        fyi(i)=n2(i)*fni(i)
        fzi(i)=n3(i)*fni(i)
       ENDDO
C------------------
C    TANGENT FORCE CALCULATION
C------------------
C---------------------------------
C       NEW FRICTION MODELS
C---------------------------------
        IF (mfrot==0) THEN
C---      Coulomb friction
          DO i=1,jlt
            xmu(i) = fric
          ENDDO
        ELSEIF (mfrot==1) THEN
C---      Viscous friction
          DO i=1,jlt
            ie=ce_loc(i)
            aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
            v2 = (vx(i) - n1(i)*aa)**2 
     .         + (vy(i) - n2(i)*aa)**2 
     .         + (vz(i) - n3(i)*aa)**2
            vv  = sqrt(max(em30,v2))
            ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
            ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
            az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
            ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
            ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
            az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
            ax  = ay1*az2 - az1*ay2
            ay  = az1*ax2 - ax1*az2
            az  = ax1*ay2 - ay1*ax2
            area = half*sqrt(ax*ax+ay*ay+az*az)
            p = -fni(i)/area
            xmu(i) = fric + (frot_p(1) + frot_p(4)*p ) * p 
     .             +(frot_p(2) + frot_p(3)*p) * vv + frot_p(5)*v2
          ENDDO
        ELSEIF(mfrot==2)THEN
C---        Loi Darmstad
          DO i=1,jlt
            ie=ce_loc(i)
            aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
            v2 = (vx(i) - n1(i)*aa)**2 
     .         + (vy(i) - n2(i)*aa)**2 
     .         + (vz(i) - n3(i)*aa)**2
            vv = sqrt(max(em30,v2))
            ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
            ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
            az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
            ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
            ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
            az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
            ax  = ay1*az2 - az1*ay2
            ay  = az1*ax2 - ax1*az2
            az  = ax1*ay2 - ay1*ax2
                area = half*sqrt(ax*ax+ay*ay+az*az)
            p = -fni(i)/area
            xmu(i) = frot_p(1)*exp(frot_p(2)*vv)*p*p
     .             + frot_p(3)*exp(frot_p(4)*vv)*p
     .             + frot_p(5)*exp(frot_p(6)*vv)
          ENDDO
        ELSEIF (mfrot==3) THEN
C---        Loi Renard
          DO i=1,jlt
            aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
            v2 = (vx(i) - n1(i)*aa)**2 
     .         + (vy(i) - n2(i)*aa)**2 
     .         + (vz(i) - n3(i)*aa)**2
            vv = sqrt(max(em30,v2))
            IF(vv>=0.AND.vv<=frot_p(5)) THEN
              dmu = frot_p(3)-frot_p(1)
              vv1 = vv / frot_p(5)
              xmu(i) = frot_p(1)+ dmu*vv1*(two-vv1)
            ELSEIF(vv>frot_p(5).AND.vv<frot_p(6)) THEN
              dmu = frot_p(4)-frot_p(3) 
              vv1 = (vv - frot_p(5))/(frot_p(6)-frot_p(5))
              xmu(i) = frot_p(3)+ dmu * (three-two*vv1)*vv1**2
            ELSE
              dmu = frot_p(2)-frot_p(4)
              vv2 = (vv - frot_p(6))**2
              xmu(i) = frot_p(2) - dmu / (one + dmu*vv2)
            ENDIF
          ENDDO
        ENDIF
C
C---------------------------------
C       INCREMENTAL (STIFFNESS) FORMULATION
C---------------------------------
        IF (ifq>=10.AND.scalk<zero) THEN
         IF (ifq==13) THEN
            alpha = max(one,alpha0*dt12)
         ELSE
            alpha = alpha0
         ENDIF
          DO i=1,jlt
            fx = stif0(i)*dx(i)
            fy = stif0(i)*dy(i)
            fz = stif0(i)*dz(i)
            fx = cand_fx(index(i)) + alpha*fx
            fy = cand_fy(index(i)) + alpha*fy
            fz = cand_fz(index(i)) + alpha*fz
            ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
            fx = fx - ftn*n1(i)
            fy = fy - ftn*n2(i)
            fz = fz - ftn*n3(i)
            ft = fx*fx + fy*fy + fz*fz
            ft = max(ft,em30)
            fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
            beta = min(one,xmu(i)*sqrt(fn/ft))
            fxt(i) = fx * beta
            fyt(i) = fy * beta
            fzt(i) = fz * beta
C-------      total force
            fxi(i)=fxi(i) + fxt(i)
            fyi(i)=fyi(i) + fyt(i)
            fzi(i)=fzi(i) + fzt(i)
          ENDDO
        ELSE 
C---------------------------------
C       TOTAL (VISCOUS) FORMULATION / NO FRICTION FILTERING 
C---------------------------------
         IF (fric>0) THEN
          DO i=1,jlt
            vnx = n1(i)*dn0(i)
            vny = n2(i)*dn0(i)
            vnz = n3(i)*dn0(i)
            vx(i) = dx(i) - vnx
            vy(i) = dy(i) - vny
            vz(i) = dz(i) - vnz
            v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
            dxt = sqrt(v2)
            aa = dxt/max(em30,v2)
            t1 = vx(i)*aa
            t2 = vy(i)*aa
            t3 = vz(i)*aa
            ftn = -xmu(i)*stif(i) * dxt
            fxt(i) = ftn * t1
            fyt(i) = ftn * t2
            fzt(i) = ftn * t3
C-------      total force
            fxi(i)=fxi(i) + fxt(i)
            fyi(i)=fyi(i) + fyt(i)
            fzi(i)=fzi(i) + fzt(i)
          ENDDO
         ENDIF
        ENDIF 
C      
C--------main part-------
c
       DO i=1,jlt
        fx=fxi(i)
        fy=fyi(i)
        fz=fzi(i)
        a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
        a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
        a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
        a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
        a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
        a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
        a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
        a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
        a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
        a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
        a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
        a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
       ENDDO
C--------secnd part-------
 
       DO i=1,jlt
         ig=nsvg(i)
         IF(ig>0)THEN
          a(1,ig)=a(1,ig)-fxi(i)
          a(2,ig)=a(2,ig)-fyi(i)
          a(3,ig)=a(3,ig)-fzi(i)
           ELSE
            nn=-ig
            ns=ind_int(nin)%P(nn)
            ffi(1,ns)=ffi(1,ns)-fxi(i)
            ffi(2,ns)=ffi(2,ns)-fyi(i)
            ffi(3,ns)=ffi(3,ns)-fzi(i)
         ENDIF
       ENDDO
C
      RETURN
      END