s4fint_reg.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/.
!||====================================================================
!||    s4fint_reg       ../engine/source/elements/solid/solide4/s4fint_reg.F
!||--- called by ------------------------------------------------------
!||    s4forc3          ../engine/source/elements/solid/solide4/s4forc3.F
!||--- uses       -----------------------------------------------------
!||    nlocal_reg_mod   ../common_source/modules/nlocal_reg_mod.F
!||====================================================================
      SUBROUTINE s4fint_reg(
     1   NLOC_DMG,VAR_REG, NEL,     OFF,
     2   VOL,     NC1,     NC2,     NC3,
     3   NC4,     PX1,     PX2,     PX3,
     4   PX4,     PY1,     PY2,     PY3,
     5   PY4,     PZ1,     PZ2,     PZ3,
     6   PZ4,     IMAT,    ITASK,   DT2T,
     7   VOL0,    NFT)
C-----------------------------------------------
C   M o d u l e s
C----------------------------------------------- 
      USE nlocal_reg_mod
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      "parit_c.inc"
#include      "scr02_c.inc"
#include      "scr18_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NFT
      INTEGER ::  NEL, IMAT, ITASK
      INTEGER, DIMENSION(NEL) :: NC1,NC2,NC3,NC4
      my_real, INTENT(INOUT)                     ::
     .  DT2T
      my_real, DIMENSION(NEL), INTENT(IN) :: 
     .   VOL,OFF,VAR_REG,VOL0,
     .   px1,px2,px3,px4,py1,py2,py3,py4,pz1,pz2,pz3,pz4
      TYPE(nlocal_str_), TARGET :: NLOC_DMG 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,K,NNOD,N1,N2,N3,N4,L_NLOC
      my_real 
     . DX, DY, DZ, L2,NTN,NTN_UNL,NTN_VNL,XI,NTVAR,A,
     . b1,b2,b3,b4,zeta,sspnl,dtnl,le_max,maxstif
      my_real, DIMENSION(NEL) :: 
     . f1,f2,f3,f4,lc
      my_real, DIMENSION(:) ,ALLOCATABLE   :: 
     . btb11,btb12,btb13,btb14,btb22,btb23,btb24,
     . btb33,btb34,btb44,sti1,sti2,sti3,sti4
      INTEGER, DIMENSION(:), ALLOCATABLE   ::
     . POS1,POS2,POS3,POS4
      my_real, POINTER, DIMENSION(:) :: 
     . vnl,fnl,unl,stifnl,mass,mass0,vnl0
      ! Coefficient for non-local stability to take into account damping
      my_real, PARAMETER :: cdamp = 0.7d0
c-----------------------------------------------------------------------
c     VAR_REG :  variable a regulariser (local cumulated plastic strain)
c     NTVAR  =  NT * VAR_REG
C=======================================================================
      l2     = nloc_dmg%LEN(imat)**2
      xi     = nloc_dmg%DAMP(imat)
      nnod   = nloc_dmg%NNOD
      l_nloc = nloc_dmg%L_NLOC
      zeta   = nloc_dmg%DENS(imat)
      sspnl  = nloc_dmg%SSPNL(imat)
      le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
      ntn    = four*four
      lc(1:nel) = zero
      ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
     . btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),pos1(nel),
     . pos2(nel),pos3(nel),pos4(nel))
      ! For nodal timestep
      IF (nodadt > 0) THEN
        ! Non-local nodal stifness
        ALLOCATE(sti1(nel),sti2(nel),sti3(nel),sti4(nel))
        ! Non-local mass
        mass =>  nloc_dmg%MASS(1:l_nloc)
        ! Initial non-local mass
        mass0 => nloc_dmg%MASS0(1:l_nloc)
      ENDIF
      vnl  => nloc_dmg%VNL(1:l_nloc)
      vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
      unl  => nloc_dmg%UNL(1:l_nloc)
c
      !-----------------------------------------------------------------------
      ! Computation of the element volume and the BtB matrix product
      !-----------------------------------------------------------------------
      ! Loop over elements
# include "vectorize.inc"
      DO i=1,nel
c
        ! Recovering the nodes of the tetra element
        n1 = nloc_dmg%IDXI(nc1(i))
        n2 = nloc_dmg%IDXI(nc2(i))
        n3 = nloc_dmg%IDXI(nc3(i))
        n4 = nloc_dmg%IDXI(nc4(i))
c        
        ! Recovering the positions of the first d.o.fs of each nodes
        pos1(i) = nloc_dmg%POSI(n1)
        pos2(i) = nloc_dmg%POSI(n2)
        pos3(i) = nloc_dmg%POSI(n3)
        pos4(i) = nloc_dmg%POSI(n4) 
c        
        ! Computation of the product BtxB 
        btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
        btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
        btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
        btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
        btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
        btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
        btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
        btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
        btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
        btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
c
      ENDDO
c      
      !-----------------------------------------------------------------------
      ! Computation of non-local forces
      !-----------------------------------------------------------------------
      ! Loop over elements
# include "vectorize.inc"
      DO i = 1, nel      
c     
        ! If the element is not broken, normal computation
        IF (off(i)/=zero) THEN 
c
          ! Computing the product NtN*UNL
          ntn_unl = (unl(pos1(i)) + unl(pos2(i)) + unl(pos3(i)) + unl(pos4(i))) / ntn
c        
          ! Computing the product XDAMP*NtN*VNL
          ntn_vnl = (vnl(pos1(i)) + vnl(pos2(i)) + vnl(pos3(i)) + vnl(pos4(i))) / ntn
          IF (nodadt > 0) THEN 
            ntn_vnl = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
     .                    sqrt(mass(pos2(i))/mass0(pos2(i))), 
     .                    sqrt(mass(pos3(i))/mass0(pos3(i))), 
     .                    sqrt(mass(pos4(i))/mass0(pos4(i))))*ntn_vnl
          ENDIF
c
          ! Computation of the product LEN**2 * BtxB
          b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i)) 
     .       + btb13(i)*unl(pos3(i)) + btb14(i)*unl(pos4(i)))
c        
          b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i)) 
     .       + btb23(i)*unl(pos3(i)) + btb24(i)*unl(pos4(i)))
c        
          b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i)) 
     .       + btb33(i)*unl(pos3(i)) + btb34(i)*unl(pos4(i)))
c        
          b4 = l2 * vol(i) * ( btb14(i)*unl(pos1(i)) + btb24(i)*unl(pos2(i)) 
     .       + btb34(i)*unl(pos3(i)) + btb44(i)*unl(pos4(i)))  
c
          ! Multiplication by the volume of the element    
          ntn_unl = ntn_unl * vol(i)
          ntn_vnl = ntn_vnl * xi * vol(i)
c
          ! Introducing the internal variable to be regularized
          ntvar   = var_reg(i)*fourth* vol(i)
c
          ! Computing the elementary non-local forces
          a = ntn_unl + ntn_vnl - ntvar
          f1(i) = a + b1
          f2(i) = a + b2
          f3(i) = a + b3
          f4(i) = a + b4
c
          ! Computing nodal equivalent stiffness
          IF (nodadt > 0) THEN 
            sti1(i) = (abs(l2*btb11(i)  + one/ntn) + abs(l2*btb12(i)  + one/ntn) + abs(l2*btb13(i)  + one/ntn) +
     .                 abs(l2*btb14(i)  + one/ntn))* vol(i)
            sti2(i) = (abs(l2*btb12(i)  + one/ntn) + abs(l2*btb22(i)  + one/ntn) + abs(l2*btb23(i)  + one/ntn) +
     .                 abs(l2*btb24(i)  + one/ntn))* vol(i) 
            sti3(i) = (abs(l2*btb13(i)  + one/ntn) + abs(l2*btb23(i)  + one/ntn) + abs(l2*btb33(i)  + one/ntn) +
     .                 abs(l2*btb34(i)  + one/ntn))* vol(i) 
            sti4(i) = (abs(l2*btb14(i)  + one/ntn) + abs(l2*btb24(i)  + one/ntn) + abs(l2*btb34(i)  + one/ntn) +
     .                 abs(l2*btb44(i)  + one/ntn))* vol(i) 
          ENDIF
c
        ! If the element is broken, computation of absorbing forces
        ELSE
c
          ! Initial element characteristic length
          lc(i) = vol0(i)**third  
c
          IF (nodadt > 0) THEN     
            ! Non-local absorbing forces
            f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*zeta*sspnl*half*
     .                          (vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
            f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*zeta*sspnl*half*
     .                          (vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
            f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*zeta*sspnl*half*
     .                          (vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
            f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*zeta*sspnl*half*
     .                          (vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)
            ! Computing nodal equivalent stiffness
            sti1(i) = em20
            sti2(i) = em20
            sti3(i) = em20
            sti4(i) = em20        
          ELSE
            ! Non-local absorbing forces
            f1(i) = zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(sqrt(three)/four)*(lc(i)**2)
            f2(i) = zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(sqrt(three)/four)*(lc(i)**2)
            f3(i) = zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(sqrt(three)/four)*(lc(i)**2)
            f4(i) = zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(sqrt(three)/four)*(lc(i)**2)          
          ENDIF
        ENDIF
      ENDDO
c-----------------------------------------------------------------------
c     Assemblage       
c-----------------------------------------------------------------------
      ! If PARITH/OFF
      IF (iparit == 0) THEN 
        fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
        IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
        DO i=1,nel
          ! Assembling the forces in the classic way 
          fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
          fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
          fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
          fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
          IF (nodadt > 0) THEN
            ! Spectral radius of stiffness matrix
            maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
            ! Computing nodal stiffness
            stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
            stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
            stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
            stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
          ENDIF
        ENDDO
c
      ! If PARITH/ON
      ELSE
c
        DO i=1,nel
          ii  = i + nft
c
          ! Spectral radius of stiffness matrix
          IF (nodadt > 0) THEN
            maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i))
          ENDIF
c
          k = nloc_dmg%IADS(1,ii)
          nloc_dmg%FSKY(k,1) = -f1(i)
          IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
c
          k = nloc_dmg%IADS(3,ii)
          nloc_dmg%FSKY(k,1) = -f2(i)
          IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
c
          k = nloc_dmg%IADS(6,ii)
          nloc_dmg%FSKY(k,1) = -f3(i)
          IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
c
          k = nloc_dmg%IADS(5,ii)
          nloc_dmg%FSKY(k,1) = -f4(i)
          IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
c
        ENDDO
      ENDIF     
c      
      !-----------------------------------------------------------------------
      ! Computing non-local timestep
      !-----------------------------------------------------------------------
      IF (nodadt == 0) THEN
        DO i = 1,nel
          ! If the element is not broken, normal computation
          IF (off(i)/=zero) THEN
            ! Non-local critical time-step in the plane
            dtnl = (two*(min(vol(i)**third,le_max))*sqrt(three*zeta))/
     .                  sqrt(twelve*l2 + (min(vol(i)**third,le_max))**2)
            ! Retaining the minimal value
            dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl)
          ENDIF
        ENDDO
      ENDIF
c
c-----------
      ! Deallocation of tables
      IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
      IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
      IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
      IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
      IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
      IF (ALLOCATED(btb23)) DEALLOCATE(btb23)      
      IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
      IF (ALLOCATED(btb33)) DEALLOCATE(btb33)
      IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
      IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
      IF (ALLOCATED(pos1))  DEALLOCATE(pos1)
      IF (ALLOCATED(pos2))  DEALLOCATE(pos2) 
      IF (ALLOCATED(pos3))  DEALLOCATE(pos3)
      IF (ALLOCATED(pos4))  DEALLOCATE(pos4)
      IF (ALLOCATED(sti1))  DEALLOCATE(sti1)
      IF (ALLOCATED(sti2))  DEALLOCATE(sti2) 
      IF (ALLOCATED(sti3))  DEALLOCATE(sti3)
      IF (ALLOCATED(sti4))  DEALLOCATE(sti4) 
c
      END