i7lagm.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/.
!||====================================================================
!||    i7lagm         ../engine/source/interfaces/int07/i7lagm.F
!||--- called by ------------------------------------------------------
!||    i7main_lmult   ../engine/source/interfaces/int07/i7main_lmult.F
!||--- calls      -----------------------------------------------------
!||    ancmsg         ../engine/source/output/message/message.F
!||    arret          ../engine/source/system/arret.F
!||--- uses       -----------------------------------------------------
!||    message_mod    ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE i7lagm(LLL     ,JLL     ,SLL    ,XLL   ,IADLL  ,
     2                  N_MUL_MX,ITASK   ,NINT   ,NKMAX , 
     3                  JLT     ,A       ,V      ,ITAG  ,XTAG   ,
     4                  GAP     ,NOINT   ,STFN   ,ITAB  ,CN_LOC ,
     5                  NX1     ,NX2     ,NX3    ,NX4   ,NY1    ,
     6                  NY2     ,NY3     ,NY4    ,NZ1   ,NZ2    ,
     7                  NZ3     ,NZ4     ,LB1    ,LB2   ,LB3    ,
     8                  LB4     ,LC1     ,LC2    ,LC3   ,LC4    ,
     9                  P1      ,P2      ,P3     ,P4    ,
     A                  IX1     ,IX2     ,IX3    ,IX4   ,NSVG   ,
     B                  GAPV    ,NEWFRONT,IBAG   ,ICONTACT,STIF ,
     C                  COMNTAG ,IADM    )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
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      "units_c.inc"
#include      "com08_c.inc"
      COMMON /lagglob/n_mult
      INTEGER N_MULT
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER N_MUL_MX,ITASK,ITIED,NINT,NKMAX ,
     .  LLL(*),JLL(*),SLL(*),IADLL(*),COMNTAG(*)
C     REAL
      my_real
     .  V(3,*),XLL(*),A(3,*),XTAG(*)
      INTEGER JLT, IBAG  ,NOINT,NEWFRONT, IADM
      INTEGER ITAB(*),ICONTACT(*),ITAG(*)
      INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
     .        NSVG(MVSIZ), CN_LOC(MVSIZ)
      my_real
     .   GAP, STFN(*)
      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)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,K,IK,IE,IS,IC,NK,III(MVSIZ,17),LLT,NFT,LE,FIRST,LAST,
     .        I16,IAD,LL
      my_real
     .   AA,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX
      INTEGER IG 
      my_real
     .   NX(MVSIZ), NY(MVSIZ), NZ(MVSIZ), PENE(MVSIZ),
     .   H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
     .   VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ), 
     .   H0, LA1, LA2, LA3, LA4,D1,D2,D3,D4,A1,A2,A3,A4
C-----------------------------------------------
C
C
C       | M | Lt| | a    | M ao
C       |---+---| |    = |
C       | L | 0 | | la   | bo
C
C        [M] a + [L]t la = [M] ao
C        [L] a = bo
C
C         a = -[M]-1[L]t la + ao
C        [L][M]-1[L]t la  = [L] ao - bo
C
C on pose:
C        [H] = [L][M]-1[L]t
C        b  = [L] ao - bo
C
C        [H] la = b
C
C        a = ao - [M]-1[L]t la
C-----------------------------------------------
C
C        la              : LAMBDA(N_MULT)
C        ao              : A(NUMNOD)
C        L               : XLL(NK,N_MULT)
C        M               : MAS(NUMNOD)
C        [L][M]-1[L]t la : HLA(N_MULT)
C        [L] ao - b      : B(N_MULT)
C        [M]-1[L]t la    : LTLA(NUMNOD)
C
C        N_MULT : nombre de contact 
C        NK : nombre de noeud pour un contact (8+1,16+1,8+8,16+16)
C
C        IC : numero du contact (1,N_MULT)
C        IK : numero de noeud local a un contact (1,NK)
C        I  : numero global du noeud (1,NUMNOD)
C
C        IADLL(N_MULT)        : IAD = IADLL(IC)
C        LLL(N_MULT*(17,51))  : I  = LLL(IAD+1,2...IADNEXT-1)
C-----------------------------------------------
C  evaluation de b:
C
C         Vs = Somme(Ni Vi)
C         Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai)
C         Somme(dt Ni Ai) - dt As =  Vs_ -Somme(Ni Vi_)
C         [L] = dt {N1,N2,..,N15,-1}
C         bo = [L] a = -[L]/dt v_
C         b  = [L] ao - bo
C         b  = [L] ao + [L]/dt v_ = [L] (v_ + ao dt)/dt
C-----------------------------------------------
C                 b  = [L] vo+/dt   +   vout
C-----------------------------------------------
 
C--------------------------------------------------------
C  SEUL CAS RESTANT : PAQUETS MIXTES
C--------------------------------------------------------
 
      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)
        nx(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i) 
        ny(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i) 
        nz(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)
         h0    = 1./max(em20,h1(i) + h2(i) + h3(i) + h4(i))
         h1(i) = h1(i) * h0
         h2(i) = h2(i) * h0
         h3(i) = h3(i) * h0
         h4(i) = h4(i) * h0
C
        ELSE
         pene(i) = p1(i)
         nx(i) = nx1(i)
         ny(i) = ny1(i)
         nz(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      DO I=1,JLT
cC       correction hourglass
c         H0 = -.25*(H1(I) - H2(I) + H3(I) - H4(I))
c         H0 = MIN(H0,H2(I),H4(I))
c         H0 = MAX(H0,-H1(I),-H3(I))
c         IF(IX3(I)==IX4(I))H0 = 0.0
c         H1(I) = H1(I) + H0
c         H2(I) = H2(I) - H0
c         H3(I) = H3(I) + H0
c         H4(I) = H4(I) - H0
c      ENDDO
C
C---------------------
C
      DO i=1,jlt
       IF( (gapv(i)-pene(i))/gapv(i) <em10 .AND.  stif(i)>zero) THEN
         stif(i) = zero
         newfront = -1
#include "lockon.inc"
            stfn(cn_loc(i)) = -abs(stfn(cn_loc(i)))
 
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,I8,A)')' NODE ',itab(nsvg(i)),
     .                   ' DE-ACTIVATED FROM INTERFACE'                 
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout ,'(A,I8,A)')' NODE ',itab(nsvg(i)),
     .                   ' DE-ACTIVATED FROM INTERFACE'   
#include "lockoff.inc"
       ENDIF
      ENDDO
C
      DO i=1,jlt
        ig=nsvg(i)
        vx(i) = v(1,ig)+dt12*a(1,ig)
     .                  - h1(i)*(v(1,ix1(i))+dt12*a(1,ix1(i)))
     .                  - h2(i)*(v(1,ix2(i))+dt12*a(1,ix2(i)))
     .                  - h3(i)*(v(1,ix3(i))+dt12*a(1,ix3(i)))
     .                  - h4(i)*(v(1,ix4(i))+dt12*a(1,ix4(i)))
        vy(i) = v(2,ig)+dt12*a(2,ig)
     .                  - h1(i)*(v(2,ix1(i))+dt12*a(2,ix1(i)))
     .                  - h2(i)*(v(2,ix2(i))+dt12*a(2,ix2(i)))
     .                  - h3(i)*(v(2,ix3(i))+dt12*a(2,ix3(i)))
     .                  - h4(i)*(v(2,ix4(i))+dt12*a(2,ix4(i)))
        vz(i) = v(3,ig)+dt12*a(3,ig)  
     .                  - h1(i)*(v(3,ix1(i))+dt12*a(3,ix1(i)))
     .                  - h2(i)*(v(3,ix2(i))+dt12*a(3,ix2(i)))
     .                  - h3(i)*(v(3,ix3(i))+dt12*a(3,ix3(i)))
     .                  - h4(i)*(v(3,ix4(i))+dt12*a(3,ix4(i)))
        vn(i) = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
#include "lockon.inc"
       IF(stif(i)/=zero.AND.pene(i)>zero.AND.vn(i)<xtag(ig))THEN
            aa = one/sqrt(nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i))
            nx(i) = nx(i)*aa
            ny(i) = ny(i)*aa
            nz(i) = nz(i)*aa
            IF(itag(nsvg(i))==0)then
             n_mult = n_mult+1
             itag(nsvg(i)) = n_mult
             xtag(nsvg(i)) = vn(i)
             IF(n_mult > n_mul_mx)THEN
#include "lockoff.inc"
               CALL ancmsg(msgid=95,anmode=aninfo)
               CALL arret(2)
             ENDIF
             iadll(n_mult+1)=iadll(n_mult) + 15
             IF(iadll(n_mult+1)-1 > nkmax)THEN
#include "lockoff.inc"
               CALL ancmsg(msgid=96,anmode=aninfo,
     .                     i1=iadll(n_mult+1)-1,
     .                     i2=nkmax)
               CALL arret(2)
             ENDIF
             iad = iadll(n_mult) - 1
            else
             xtag(nsvg(i)) = vn(i)
             iad = iadll(itag(nsvg(i))) - 1
             ll = lll(iad+1)
             comntag(ll)= comntag(ll) - 1
             ll = lll(iad+2)
             comntag(ll)= comntag(ll) - 1
             ll = lll(iad+3)
             comntag(ll)= comntag(ll) - 1
             ll = lll(iad+4)
             comntag(ll)= comntag(ll) - 1
             ll = lll(iad+5)
             comntag(ll)= comntag(ll) - 1
            ENDIF
C
            lll(iad+1)    = ix1(i)
            jll(iad+1)    = 1
            sll(iad+1)    = 0
            xll(iad+1)    = nx(i)*h1(i)
C---
            lll(iad+2)    = ix2(i)
            jll(iad+2)    = 1
            sll(iad+2)    = 0
            xll(iad+2)    = nx(i)*h2(i)
C---
            lll(iad+3)    = ix3(i)
            jll(iad+3)    = 1
            sll(iad+3)    = 0
            xll(iad+3)    = nx(i)*h3(i)
C---
            lll(iad+4)    = ix4(i)
            jll(iad+4)    = 1
            sll(iad+4)    = 0
            xll(iad+4)    = nx(i)*h4(i)
C---
            lll(iad+5)    = nsvg(i)
            jll(iad+5)    = 1
            sll(iad+5)    = nint
            xll(iad+5)    = -nx(i)
C-----------------------------
            lll(iad+6)    = ix1(i)
            jll(iad+6)    = 2
            sll(iad+6)    = 0
            xll(iad+6)    = ny(i)*h1(i)
C---
            lll(iad+7)    = ix2(i)
            jll(iad+7)    = 2
            sll(iad+7)    = 0
            xll(iad+7)    = ny(i)*h2(i)
C---
            lll(iad+8)    = ix3(i)
            jll(iad+8)    = 2
            sll(iad+8)    = 0
            xll(iad+8)    = ny(i)*h3(i)
C---
            lll(iad+9)    = ix4(i)
            jll(iad+9)    = 2
            sll(iad+9)    = 0
            xll(iad+9)    = ny(i)*h4(i)
C---
            lll(iad+10)    = nsvg(i)
            jll(iad+10)    = 2
            sll(iad+10)    = nint
            xll(iad+10)    = -ny(i)
C------------------------------------
            lll(iad+11)    = ix1(i)
            jll(iad+11)    = 3
            sll(iad+11)    = 0
            xll(iad+11)    = nz(i)*h1(i)
C---
            lll(iad+12)    = ix2(i)
            jll(iad+12)    = 3
            sll(iad+12)    = 0
            xll(iad+12)    = nz(i)*h2(i)
C---
            lll(iad+13)    = ix3(i)
            jll(iad+13)    = 3
            sll(iad+13)    = 0
            xll(iad+13)    = nz(i)*h3(i)
C---
            lll(iad+14)    = ix4(i)
            jll(iad+14)    = 3
            sll(iad+14)    = 0
            xll(iad+14)    = nz(i)*h4(i)
C---
            lll(iad+15)    = nsvg(i)
            jll(iad+15)    = 3
            sll(iad+15)    = nint
            xll(iad+15)    = -nz(i)
C--------------------------------------
            ll = ix1(i)
            comntag(ll) = comntag(ll) + 1
            ll = ix2(i)
            comntag(ll) = comntag(ll) + 1
            ll = ix3(i)
            comntag(ll) = comntag(ll) + 1
            ll = ix4(i)
            comntag(ll) = comntag(ll) + 1
            ll = nsvg(i)
            comntag(ll) = comntag(ll) + 1
C
        ENDIF
#include "lockoff.inc"
      ENDDO
C
      IF(ibag/=0.OR.iadm/=0)THEN
       DO i=1,jlt
        IF(pene(i)/=zero)THEN
          icontact(nsvg(i))=1
          icontact(ix1(i))=1
          icontact(ix2(i))=1
          icontact(ix3(i))=1
          icontact(ix4(i))=1
        ENDIF
       ENDDO
      ENDIF
C
C-----------------------------------------------
      RETURN
      END