lag_direct.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/.
!||====================================================================
!||    lag_direct    ../engine/source/tools/lagmul/lag_direct.F
!||--- called by ------------------------------------------------------
!||    i2lagm        ../engine/source/tools/lagmul/lag_i2main.F
!||    lag_gjnt      ../engine/source/tools/lagmul/lag_gjnt.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 lag_direct(
     1      IADLL    ,LLL      ,JLL      ,XLL      ,LTSM     ,
     2      V        ,VR       ,A        ,AR       ,MS       ,
     3      IN       ,NC_INI   ,NCL      )
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"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "lagmult.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NC_INI, NCL, IADLL(*), LLL(*), JLL(*)
C     REAL
      my_real
     .   ltsm(6,*),xll(*),ms(*),in(*),v(3,*),vr(3,*),a(3,*),ar(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,K,L,IC,JC,IK
      my_real
     .        HLOC(MXDLEN,MXDLEN),RLOC(MXDLEN),S,HIJ
C======================================================================|
      IF (ncl>mxdlen) THEN
        CALL ancmsg(msgid=111,anmode=aninfo,
     .            i1=ncl)
        CALL arret(2)
      ENDIF
C---  Local H matrix
      DO k=1,ncl
        rloc(k) = zero
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF (j>3) THEN
            ltsm(j,i) = xll(ik)/in(i)
          ELSE
            ltsm(j,i) = xll(ik)/ms(i)
          ENDIF
        ENDDO
        DO l = 1,k
          jc = nc_ini + l
          hij = zero
          DO ik=iadll(jc),iadll(jc+1)-1
            hij = hij + xll(ik)*ltsm(jll(ik),lll(ik))
          ENDDO
          hloc(l,k)  = hij
        ENDDO
        DO ik=iadll(ic),iadll(ic+1)-1
          ltsm(jll(ik),lll(ik)) = zero
        ENDDO
      ENDDO
      DO k = 2,ncl
        DO l = 1,k
          hloc(k,l) = hloc(l,k)
        ENDDO
      ENDDO
C---  second membre
      DO k = 1,ncl
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF (j>3) THEN
            j = j-3
            rloc(k) = rloc(k) + xll(ik)*(vr(j,i)/dt12+ar(j,i))
          ELSE
            rloc(k) = rloc(k) + xll(ik)*(v(j,i)/dt12+a(j,i))
          ENDIF
        ENDDO
      ENDDO
C---  Factorise H (Full Cholesky)
      DO j=1,ncl
              IF (hloc(j,j)<=zero) THEN
          CALL ancmsg(msgid=112,anmode=aninfo,
     .            i1=j)
          CALL arret(2)
              ENDIF
        hloc(j,j) = sqrt(hloc(j,j))
        DO k=1,j-1
          DO i=j+1,ncl
                  hloc(i,j) = hloc(i,j) - hloc(i,k)*hloc(j,k)
          ENDDO
        ENDDO
        DO i=j+1,ncl
          hloc(i,j) = hloc(i,j)/hloc(j,j)
          hloc(i,i) = hloc(i,i) - hloc(i,j)*hloc(i,j)
        ENDDO
      ENDDO
C---  back subst Ly = r,
      DO i=1,ncl
        s = rloc(i)
        DO j=1,i-1
          s = s - hloc(i,j)*rloc(j)
        ENDDO
              rloc(i) = s / hloc(i,i)
      ENDDO
C---  back subst Lz = y
      DO i=ncl,1,-1
        s = rloc(i)
        DO j=i+1,ncl
          s = s - hloc(j,i)*rloc(j)
        ENDDO
              rloc(i) = s / hloc(i,i)
      ENDDO
C---  update accelerations
      DO k=1,ncl
        ic = nc_ini + k
        DO ik=iadll(ic),iadll(ic+1)-1
          i = lll(ik)
          j = jll(ik)
          IF(j>3) THEN
            j = j-3
            ar(j,i) = ar(j,i) - xll(ik)*rloc(k)/in(i)
          ELSE
            a(j,i)  = a(j,i)  - xll(ik)*rloc(k)/ms(i)
          ENDIF
        ENDDO
      ENDDO
C---
      RETURN
      END