s6cfint_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/.
!||====================================================================
!||    s6cfint_reg      ../engine/source/elements/thickshell/solide6c/s6cfint_reg.F
!||--- called by ------------------------------------------------------
!||    s6cforc3         ../engine/source/elements/thickshell/solide6c/s6cforc3.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod     ../common_source/modules/mat_elem/elbufdef_mod.f90
!||    nlocal_reg_mod   ../common_source/modules/nlocal_reg_mod.F
!||====================================================================
      SUBROUTINE s6cfint_reg(
     1   NLOC_DMG,VAR_REG, NEL,     OFF,
     2   VOL,     NC1,     NC2,     NC3,
     3   NC4,     NC5,     NC6,     PX1,
     4   PX2,     PX3,     PX4,     PY1,
     5   PY2,     PY3,     PY4,     PZ1,
     6   PZ2,     PZ3,     PZ4,     IMAT,
     7   ITASK,   DT2T,    VOL0,    NFT,
     8   NLAY,    WS,      AS,      AREA,
     9   BUFNLTS)
C-----------------------------------------------
C   M o d u l e s
C----------------------------------------------- 
      USE nlocal_reg_mod
      USE elbufdef_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"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: 
     .   NFT,NLAY,NEL,IMAT,ITASK 
      INTEGER, INTENT(IN), DIMENSION(NEL) :: 
     .   NC1,NC2,NC3,NC4,NC5,NC6
      my_real, INTENT(INOUT) ::
     .   DT2T
      my_real, DIMENSION(9,9), INTENT(IN) :: 
     .   WS,AS
      my_real, DIMENSION(NEL,NLAY), INTENT(INOUT) ::
     .   var_reg
      my_real, DIMENSION(NEL), INTENT(IN) :: 
     .   vol,off,vol0,px1,px2,px3,px4,area,
     .   py1,py2,py3,py4,pz1,pz2,pz3,pz4
      TYPE(nlocal_str_), INTENT(INOUT), TARGET :: NLOC_DMG 
      TYPE(BUF_NLOCTS_), INTENT(INOUT), TARGET :: BUFNLTS
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,J,K,NNOD,N1,N2,N3,N4,N5,N6,
     .        l_nloc,nddl,ndnod
      my_real 
     . l2,ntn,ntn_unl,ntn_vnl,xi,ntvar,a,dtnl,le_max,
     . b1,b2,b3,b4,b5,b6,zeta,sspnl,maxstif,
     . bth1,bth2,nth1,nth2,dt2p,dtnod,k1,k2,k12,
     . dtnl_th
      my_real, DIMENSION(NEL,NLAY) :: 
     . f1,f2,f3,f4,f5,f6
      my_real, DIMENSION(NEL) ::
     . lc,thk,lthk,pxx1,pxx2,pxx3,
     . pxx4,pxx5,pxx6,pyy1,pyy2,pyy3,
     . pyy4,pyy5,pyy6,pzz1,pzz2,pzz3,
     . pzz4,pzz5,pzz6
      my_real, DIMENSION(:) ,ALLOCATABLE   :: 
     . btb11,btb12,btb13,btb14,btb15,btb16,
     . btb22,btb23,btb24,btb25,btb26,btb33,
     . btb34,btb35,btb36,btb44,btb45,btb46,
     . btb55,btb56,btb66
      my_real, DIMENSION(:,:) ,ALLOCATABLE ::
     . sti1,sti2,sti3,sti4,sti5,sti6
      INTEGER, DIMENSION(:), ALLOCATABLE   ::
     . POS1,POS2,POS3,POS4,POS5,POS6
      my_real, POINTER, DIMENSION(:) :: 
     . vnl,fnl,unl,stifnl,mass,mass0,vnl0
      my_real, POINTER, DIMENSION(:,:)     :: 
     . massth,unlth,vnlth,fnlth
      my_real, DIMENSION(:,:), ALLOCATABLE :: 
     . stifnlth,dtn
      my_real, target :: nothing(1)
      ! Safety coefficient for non-local stability vs mechanical stability
      ! (it has been slightly increased vs nloc_dmg_init.F)
      my_real, PARAMETER :: csta  = 40.0d0
      ! Coefficient for non-local stability to take into account damping
      my_real, PARAMETER :: cdamp = 0.7d0
      my_real
     . zs(10,9)      
      ! Position of nodes in the thickshell thickness
      DATA zs / 
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,
     2 -1.              ,0.               ,1.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,
     3 -1.              ,-.549193338482966,0.549193338482966,
     3 1.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     3 0.               ,
     4 -1.              ,-.600558677589454,0.               ,
     4 0.600558677589454,1.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     4 0.               ,
     5 -1.              ,-.812359691877328,-.264578928334038,
     5 0.264578928334038,0.812359691877328,1.               ,
     5 0.               ,0.               ,0.               ,
     5 0.               ,
     6 -1.              ,-.796839450334708,-.449914286274731,
     6 0.               ,0.449914286274731,0.796839450334708,
     6 1.               ,0.               ,0.               ,
     6 0.               ,
     7 -1.              ,-.898215824685518,-.584846546513270,
     7 -.226843756241524,0.226843756241524,0.584846546513270,
     7 0.898215824685518,1.               ,0.               ,
     7 0.               ,
     8 -1.              ,-.878478166955581,-.661099443664978,
     8 -.354483526205989,0.               ,0.354483526205989,
     8 0.661099443664978,0.878478166955581,1.               ,
     8 0.               ,
     9 -1.              ,-.936320479015252,-.735741735638020,
     9 -.491001129763160,-.157505717044458,0.157505717044458,
     9 0.491001129763160,0.735741735638020,0.936320479015252,
     9 1.               /
C=======================================================================
      stifnl => nothing
      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    = six*six
      lc(1:nel) = zero
      ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb15(nel),
     .         btb16(nel),btb22(nel),btb23(nel),btb24(nel),btb25(nel),
     .         btb26(nel),btb33(nel),btb34(nel),btb35(nel),btb36(nel),
     .         btb44(nel),btb45(nel),btb46(nel),btb55(nel),btb56(nel),
     .         btb66(nel),pos1(nel) ,pos2(nel) ,pos3(nel) ,pos4(nel) ,
     .         pos5(nel) ,pos6(nel) )
      ! For nodal timestep
      IF (nodadt > 0) THEN
        ! Non-local nodal stifness
        ALLOCATE(sti1(nel,nlay),sti2(nel,nlay),sti3(nel,nlay),
     .           sti4(nel,nlay),sti5(nel,nlay),sti6(nel,nlay))
        ! Non-local mass
        mass =>  nloc_dmg%MASS(1:l_nloc)
        ! Initial non-local mass
        mass0 => nloc_dmg%MASS0(1:l_nloc)
      ELSE
       NULLIFY(mass)
       NULLIFY(mass0)
       ALLOCATE(sti1(1,1),sti2(1,1),sti3(1,1),
     .           sti4(1,1),sti5(1,1),sti6(1,1))
      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 brick 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))
        n5 = nloc_dmg%IDXI(nc5(i))
        n6 = nloc_dmg%IDXI(nc6(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)
        pos5(i) = nloc_dmg%POSI(n5)
        pos6(i) = nloc_dmg%POSI(n6)
c
        ! Computation of derivatives of shape functions
        pxx1(i) = px1(i)-px4(i)
        pyy1(i) = py1(i)-py4(i)
        pzz1(i) = pz1(i)-pz4(i)
c
        pxx2(i) = px2(i)-px4(i)
        pyy2(i) = py2(i)-py4(i)
        pzz2(i) = pz2(i)-pz4(i)
c
        pxx3(i) = px3(i)-px4(i)
        pyy3(i) = py3(i)-py4(i)
        pzz3(i) = pz3(i)-pz4(i)
c
        pxx4(i) = px1(i)+px4(i)
        pyy4(i) = py1(i)+py4(i)
        pzz4(i) = pz1(i)+pz4(i)
c
        pxx5(i) = px2(i)+px4(i)
        pyy5(i) = py2(i)+py4(i)
        pzz5(i) = pz2(i)+pz4(i)
c
        pxx6(i) = px3(i)+px4(i)
        pyy6(i) = py3(i)+py4(i)
        pzz6(i) = pz3(i)+pz4(i)
c 
        ! Computation of the product BtxB 
        btb11(i) = pxx1(i)**2 + pyy1(i)**2 + pzz1(i)**2
        btb12(i) = pxx1(i)*pxx2(i) + pyy1(i)*pyy2(i) + pzz1(i)*pzz2(i)
        btb13(i) = pxx1(i)*pxx3(i) + pyy1(i)*pyy3(i) + pzz1(i)*pzz3(i)
        btb14(i) = pxx1(i)*pxx4(i) + pyy1(i)*pyy4(i) + pzz1(i)*pzz4(i)
        btb15(i) = pxx1(i)*pxx5(i) + pyy1(i)*pyy5(i) + pzz1(i)*pzz5(i)
        btb16(i) = pxx1(i)*pxx6(i) + pyy1(i)*pyy6(i) + pzz1(i)*pzz6(i)
c 
        btb22(i) = pxx2(i)**2 + pyy2(i)**2 + pzz2(i)**2
        btb23(i) = pxx2(i)*pxx3(i) + pyy2(i)*pyy3(i) + pzz2(i)*pzz3(i)
        btb24(i) = pxx2(i)*pxx4(i) + pyy2(i)*pyy4(i) + pzz2(i)*pzz4(i)
        btb25(i) = pxx2(i)*pxx5(i) + pyy2(i)*pyy5(i) + pzz2(i)*pzz5(i)
        btb26(i) = pxx2(i)*pxx6(i) + pyy2(i)*pyy6(i) + pzz2(i)*pzz6(i)
c
        btb33(i) = pxx3(i)**2 + pyy3(i)**2 + pzz3(i)**2
        btb34(i) = pxx3(i)*pxx4(i) + pyy3(i)*pyy4(i) + pzz3(i)*pzz4(i)
        btb35(i) = pxx3(i)*pxx5(i) + pyy3(i)*pyy5(i) + pzz3(i)*pzz5(i)
        btb36(i) = pxx3(i)*pxx6(i) + pyy3(i)*pyy6(i) + pzz3(i)*pzz6(i)
c
        btb44(i) = pxx4(i)**2 + pyy4(i)**2 + pzz4(i)**2
        btb45(i) = pxx4(i)*pxx5(i) + pyy4(i)*pyy5(i) + pzz4(i)*pzz5(i)
        btb46(i) = pxx4(i)*pxx6(i) + pyy4(i)*pyy6(i) + pzz4(i)*pzz6(i)
c
        btb55(i) = pxx5(i)**2 + pyy5(i)**2 + pzz5(i)**2
        btb56(i) = pxx5(i)*pxx6(i) + pyy5(i)*pyy6(i) + pzz5(i)*pzz6(i)
c
        btb66(i) = pxx6(i)**2 + pyy6(i)**2 + pzz6(i)**2
c        
      ENDDO
c
      !-----------------------------------------------------------------------
      ! Pre-treatment non-local regularization in the thickshell thickness
      !-----------------------------------------------------------------------
      IF ((l2>zero).AND.(nlay > 1)) THEN 
c
        ! Compute thickshell thickness
        DO i = 1,nel
          thk(i)  = vol(i)/area(i)
          lthk(i) = (zs(nlay+1,nlay)-zs(nlay,nlay))*thk(i)*half
        ENDDO
c
        ! Allocation of the velocities predictor
        nddl = nlay
        IF (nodadt > 0) THEN 
          ALLOCATE(stifnlth(nel,nddl+1))
          ALLOCATE(dtn(nel,nddl+1))
        ELSE
          ALLOCATE(dtn(1,1))
          ALLOCATE(stifnlth(1,1))
          dtn(1,1) = ep20
          stifnlth(1,1) = ep20
        ENDIF
        ndnod = nddl+1
c 
        ! Pointing the non-local values in the thickness of the corresponding element
        massth => bufnlts%MASSTH(1:nel,1:ndnod)
        unlth  => bufnlts%UNLTH(1:nel ,1:ndnod)
        vnlth  => bufnlts%VNLTH(1:nel ,1:ndnod)
        fnlth  => bufnlts%FNLTH(1:nel ,1:ndnod)    
c
        DO k = 1,ndnod
          DO i = 1,nel
            ! Resetting non-local forces
            fnlth(i,k) = zero
            ! Resetting non-local nodal stiffness
            IF (nodadt > 0) THEN
              stifnlth(i,k) = em20
            ENDIF
          ENDDO
        ENDDO
c
        ! Computation of non-local forces in the shell thickness
        DO k = 1, nddl
c        
          ! Computation of shape functions value
          nth1 = (as(k,nddl)   - zs(k+1,nddl)) / 
     .           (zs(k,nddl)   - zs(k+1,nddl))
          nth2 = (as(k,nddl)   - zs(k,nddl))   / 
     .           (zs(k+1,nddl) - zs(k,nddl))
c          
          ! Loop over elements
          DO i = 1,nel
c
            ! Computation of B-matrix values
            bth1 = (one/(zs(k,nddl)   - zs(k+1,nddl)))*(two/thk(i))
            bth2 = (one/(zs(k+1,nddl) - zs(k,nddl)))*(two/thk(i))   
c         
            ! Computation of the non-local K matrix
            k1   = l2*(bth1**2)  + nth1**2
            k12  = l2*(bth1*bth2)+ (nth1*nth2)
            k2   = l2*(bth2**2)  + nth2**2
c
            ! Computation of the non-local forces
            fnlth(i,k)   = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1) 
     .                                + xi*((nth1**2)*vnlth(i,k) 
     .                                + (nth1*nth2)*vnlth(i,k+1))
     .                                - (nth1*var_reg(i,k)))*half*ws(k,nddl)*vol(i)  
            fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i,k) + k2*unlth(i,k+1)
     .                                + xi*(nth1*nth2*vnlth(i,k) 
     .                                + (nth2**2)*vnlth(i,k+1))
     .                                - nth2*var_reg(i,k))*half*ws(k,nddl)*vol(i)  
c
            ! Computation of non-local nodal stiffness
            IF (nodadt > 0) THEN 
              stifnlth(i,k)   = stifnlth(i,k)   + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)
              stifnlth(i,k+1) = stifnlth(i,k+1) + max(abs(k1)+abs(k12),abs(k12)+abs(k2))*half*ws(k,nddl)*vol(i)            
            ENDIF
c
          ENDDO
        ENDDO
c       
        ! Updating non-local mass with /DT/NODA
        IF (nodadt > 0) THEN 
C
          ! Initial computation of the nodal timestep
          dtnod = ep20
          DO k = 1,ndnod
            DO i = 1,nel
              dtn(i,k) = dtfac1(11)*cdamp*sqrt(two*massth(i,k)/max(stifnlth(i,k),em20)) 
              dtnod    = min(dtn(i,k),dtnod)
            ENDDO
          ENDDO
C
          ! /DT/NODA/CSTX - Constant timestep with added mass
          IF ((idtmin(11)==3).OR.(idtmin(11)==4).OR.(idtmin(11)==8)) THEN  
            ! Added mass computation if necessary
            IF (dtnod < dtmin1(11)*(sqrt(csta))) THEN
              DO k = 1,ndnod
                DO i = 1,nel
                  IF (dtn(i,k) < dtmin1(11)) THEN
                    dt2p        = dtmin1(11)/(dtfac1(11)*cdamp)
                    massth(i,k) = max(massth(i,k),csta*half*stifnlth(i,k)*dt2p*dt2p*onep00001)
                  ENDIF
                ENDDO
              ENDDO
            ENDIF
            dtnod = dtmin1(11)*(sqrt(csta))
          ENDIF
C
          ! Classical nodal timestep check
          IF (dtnod < dt2t) THEN
            dt2t = min(dt2t,dtnod)
          ENDIF
        ENDIF
c
        DO k = 1,ndnod
          DO i = 1,nel
            ! Updating the non-local in-thickness velocities   
            vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt12
          ENDDO
        ENDDO
c          
        DO k = 1,ndnod
          DO i = 1,nel
            ! computing the non-local in-thickness cumulated values
            unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt1
          ENDDO
        ENDDO
c
        ! Transfert at the integration point
        DO k = 1, nddl
          !Computation of shape functions value
          nth1 = (as(k,nddl)   - zs(k+1,nddl))/
     .           (zs(k,nddl)   - zs(k+1,nddl))
          nth2 = (as(k,nddl)   - zs(k,nddl))/
     .           (zs(k+1,nddl) - zs(k,nddl))
          ! Loop over elements
          DO i = 1,nel
            !Integration points non-local variables
            var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i,k+1)
          ENDDO  
        ENDDO
      ENDIF
c      
      !-----------------------------------------------------------------------
      ! Computation of non-local forces
      !-----------------------------------------------------------------------
      ! Loop over elements
      DO k = 1,nlay
c
        ! 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)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) + unl(pos4(i)+k-1)
     .              +  unl(pos5(i)+k-1) + unl(pos6(i)+k-1)) / ntn
c        
            ! Computing the product XDAMP*NtN*VNL
            ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) + vnl(pos4(i)+k-1)
     .              +  vnl(pos5(i)+k-1) + vnl(pos6(i)+k-1)) / ntn
            IF (nodadt > 0) THEN 
              ntn_vnl = min(sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1)), 
     .                      sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1)), 
     .                      sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1)), 
     .                      sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1)),
     .                      sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1)), 
     .                      sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1)))*ntn_vnl
            ENDIF
c        
            ! Computation of the product LEN**2 * BtxB
            b1 = l2 * vol(i) * ws(k,nlay) *half * ( btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1) 
     .                + btb13(i)*unl(pos3(i)+k-1) + btb14(i)*unl(pos4(i)+k-1) + btb15(i)*unl(pos5(i)+k-1)
     .                + btb16(i)*unl(pos6(i)+k-1) )
c        
            b2 = l2 * vol(i) * ws(k,nlay) *half * ( btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1) 
     .                + btb23(i)*unl(pos3(i)+k-1) + btb24(i)*unl(pos4(i)+k-1) + btb25(i)*unl(pos5(i)+k-1)
     .                + btb26(i)*unl(pos6(i)+k-1) )
c        
            b3 = l2 * vol(i) * ws(k,nlay) *half * ( btb13(i)*unl(pos1(i)+k-1) + btb23(i)*unl(pos2(i)+k-1) 
     .                + btb33(i)*unl(pos3(i)+k-1) + btb34(i)*unl(pos4(i)+k-1) + btb35(i)*unl(pos5(i)+k-1)
     .                + btb36(i)*unl(pos6(i)+k-1) )
c        
            b4 = l2 * vol(i) * ws(k,nlay) *half * ( btb14(i)*unl(pos1(i)+k-1) + btb24(i)*unl(pos2(i)+k-1) 
     .                + btb34(i)*unl(pos3(i)+k-1) + btb44(i)*unl(pos4(i)+k-1) + btb45(i)*unl(pos5(i)+k-1)
     .                + btb46(i)*unl(pos6(i)+k-1) )
c       
            b5 = l2 * vol(i) * ws(k,nlay) *half * ( btb15(i)*unl(pos1(i)+k-1) + btb25(i)*unl(pos2(i)+k-1) 
     .                + btb35(i)*unl(pos3(i)+k-1) + btb45(i)*unl(pos4(i)+k-1) + btb55(i)*unl(pos5(i)+k-1)
     .                + btb56(i)*unl(pos6(i)+k-1) )
c        
            b6 = l2 * vol(i) * ws(k,nlay) *half * ( btb16(i)*unl(pos1(i)+k-1) + btb26(i)*unl(pos2(i)+k-1) 
     .                + btb36(i)*unl(pos3(i)+k-1) + btb46(i)*unl(pos4(i)+k-1) + btb56(i)*unl(pos5(i)+k-1)
     .                + btb66(i)*unl(pos6(i)+k-1) )
c
            ! Multiplication by the volume of the element    
            ntn_unl = ntn_unl * vol(i) * ws(k,nlay) * half
            ntn_vnl = ntn_vnl * xi * vol(i) * ws(k,nlay) * half
c
            ! Introducing the internal variable to be regularized
            ntvar   = var_reg(i,k)*one_over_6* vol(i) * ws(k,nlay) * half
c
            ! Computing the elementary non-local forces
            a = ntn_unl + ntn_vnl - ntvar
            f1(i,k) = a + b1
            f2(i,k) = a + b2
            f3(i,k) = a + b3
            f4(i,k) = a + b4
            f5(i,k) = a + b5
            f6(i,k) = a + b6
c
            ! Computing nodal equivalent stiffness
            IF (nodadt > 0) THEN 
              sti1(i,k) = (abs(l2*btb11(i) + one/ntn) + abs(l2*btb12(i) + one/ntn) + abs(l2*btb13(i) + one/ntn) +
     .                     abs(l2*btb14(i) + one/ntn) + abs(l2*btb15(i) + one/ntn) + abs(l2*btb16(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
              sti2(i,k) = (abs(l2*btb12(i) + one/ntn) + abs(l2*btb22(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) +
     .                     abs(l2*btb24(i) + one/ntn) + abs(l2*btb25(i) + one/ntn) + abs(l2*btb26(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
              sti3(i,k) = (abs(l2*btb13(i) + one/ntn) + abs(l2*btb23(i) + one/ntn) + abs(l2*btb33(i) + one/ntn) +
     .                     abs(l2*btb34(i) + one/ntn) + abs(l2*btb35(i) + one/ntn) + abs(l2*btb36(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
              sti4(i,k) = (abs(l2*btb14(i) + one/ntn) + abs(l2*btb24(i) + one/ntn) + abs(l2*btb34(i) + one/ntn) +
     .                     abs(l2*btb44(i) + one/ntn) + abs(l2*btb45(i) + one/ntn) + abs(l2*btb46(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
              sti5(i,k) = (abs(l2*btb15(i) + one/ntn) + abs(l2*btb25(i) + one/ntn) + abs(l2*btb35(i) + one/ntn) +
     .                     abs(l2*btb45(i) + one/ntn) + abs(l2*btb55(i) + one/ntn) + abs(l2*btb56(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
              sti6(i,k) = (abs(l2*btb16(i) + one/ntn) + abs(l2*btb26(i) + one/ntn) + abs(l2*btb36(i) + one/ntn) +
     .                     abs(l2*btb46(i) + one/ntn) + abs(l2*btb56(i) + one/ntn) + abs(l2*btb66(i) + one/ntn))
     .                     *vol(i)*ws(k,nlay)*half
            ENDIF
c            
          ! If the element is broken, the non-local wave is absorbed  
          ELSE
c
            ! Initial element characteristic length
            lc(i) = (vol0(i)*ws(k,nlay)*half)**third  
c
            IF (nodadt > 0) THEN
     
              ! Non-local absorbing forces
              f1(i,k) = sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
              f2(i,k) = sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
              f3(i,k) = sqrt(mass(pos3(i)+k-1)/mass0(pos3(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
              f4(i,k) = sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
              f5(i,k) = sqrt(mass(pos5(i)+k-1)/mass0(pos5(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
              f6(i,k) = sqrt(mass(pos6(i)+k-1)/mass0(pos6(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)
              ! Computing nodal equivalent stiffness
              sti1(i,k) = em20
              sti2(i,k) = em20
              sti3(i,k) = em20
              sti4(i,k) = em20
              sti5(i,k) = em20
              sti6(i,k) = em20
            ELSE
              ! Non-local absorbing forces
              f1(i,k) = zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(two*third)*(lc(i)**2)
              f2(i,k) = zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(two*third)*(lc(i)**2)
              f3(i,k) = zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(two*third)*(lc(i)**2)
              f4(i,k) = zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(two*third)*(lc(i)**2)
              f5(i,k) = zeta*sspnl*half*(vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(two*third)*(lc(i)**2)
              f6(i,k) = zeta*sspnl*half*(vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(two*third)*(lc(i)**2)         
            ENDIF
          ENDIF
        ENDDO
      ENDDO
c      
      !-----------------------------------------------------------------------
      ! Assembling of the non-local forces
      !-----------------------------------------------------------------------
c
      ! If PARITH/OFF
      IF (iparit == 0) THEN 
        ! Recovering non-local internal forces
        fnl => nloc_dmg%FNL(1:l_nloc,itask+1)                       ! Non-local forces
        IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
        ! Loop over elements
        DO i=1,nel
          ! Loop over non-local degrees of freedom
# include "vectorize.inc"
          DO k=1,nlay
            ! Assembling the forces in the classic way 
            fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
            fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
            fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
            fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
            fnl(pos5(i)+k-1) = fnl(pos5(i)+k-1) - f5(i,k)
            fnl(pos6(i)+k-1) = fnl(pos6(i)+k-1) - f6(i,k)  
            IF (nodadt > 0) THEN
              ! Spectral radius of stiffness matrix
              maxstif = max(sti1(i,k),sti2(i,k),sti3(i,k),
     .                      sti4(i,k),sti5(i,k),sti6(i,k))
              ! Computing nodal stiffness
              stifnl(pos1(i)+k-1) = stifnl(pos1(i)+k-1) + maxstif
              stifnl(pos2(i)+k-1) = stifnl(pos2(i)+k-1) + maxstif
              stifnl(pos3(i)+k-1) = stifnl(pos3(i)+k-1) + maxstif
              stifnl(pos4(i)+k-1) = stifnl(pos4(i)+k-1) + maxstif
              stifnl(pos5(i)+k-1) = stifnl(pos5(i)+k-1) + maxstif
              stifnl(pos6(i)+k-1) = stifnl(pos6(i)+k-1) + maxstif
            ENDIF
          ENDDO
        ENDDO
c
      ! If PARITH/ON
      ELSE
        ! Loop over additional d.o.fs
        DO j = 1,nlay
c
          ! Loop over elements
          DO i=1,nel
            ii  = i + nft
c
            ! Spectral radius of stiffness matrix
            IF (nodadt > 0) THEN
              maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j),
     .                      sti5(i,j),sti6(i,j))
            ENDIF
c            
            k = nloc_dmg%IADS(1,ii)
            nloc_dmg%FSKY(k,j) = -f1(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(2,ii)
            nloc_dmg%FSKY(k,j) = -f2(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(3,ii)
            nloc_dmg%FSKY(k,j) = -f3(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(5,ii)
            nloc_dmg%FSKY(k,j) = -f4(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(6,ii)
            nloc_dmg%FSKY(k,j) = -f5(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(7,ii)
            nloc_dmg%FSKY(k,j) = -f6(i,j)
            IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
          ENDDO
        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)
            ! Non-local critical time-step in the thickness
            IF (nlay > 1) THEN 
              dtnl_th = (two*(min(lthk(i),le_max))*sqrt(three*zeta))/
     .                sqrt(twelve*l2 + (min(lthk(i),le_max))**2)
            ELSE 
              dtnl_th = ep20
            ENDIF
            ! Retaining the minimal value
            dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl_th,dtfac1(1)*cdamp*dtnl)
          ENDIF
        ENDDO
      ENDIF
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(btb15))    DEALLOCATE(btb15)
      IF (ALLOCATED(btb16))    DEALLOCATE(btb16)
      IF (ALLOCATED(btb22))    DEALLOCATE(btb22)
      IF (ALLOCATED(btb23))    DEALLOCATE(btb23)      
      IF (ALLOCATED(btb24))    DEALLOCATE(btb24)
      IF (ALLOCATED(btb25))    DEALLOCATE(btb25)
      IF (ALLOCATED(btb26))    DEALLOCATE(btb26)
      IF (ALLOCATED(btb33))    DEALLOCATE(btb33)
      IF (ALLOCATED(btb34))    DEALLOCATE(btb34)
      IF (ALLOCATED(btb35))    DEALLOCATE(btb35)
      IF (ALLOCATED(btb36))    DEALLOCATE(btb36)
      IF (ALLOCATED(btb44))    DEALLOCATE(btb44)
      IF (ALLOCATED(btb45))    DEALLOCATE(btb45)
      IF (ALLOCATED(btb46))    DEALLOCATE(btb46)
      IF (ALLOCATED(btb55))    DEALLOCATE(btb55)
      IF (ALLOCATED(btb56))    DEALLOCATE(btb56)
      IF (ALLOCATED(btb66))    DEALLOCATE(btb66)
      IF (ALLOCATED(pos1))     DEALLOCATE(pos1)
      IF (ALLOCATED(pos2))     DEALLOCATE(pos2) 
      IF (ALLOCATED(pos3))     DEALLOCATE(pos3)
      IF (ALLOCATED(pos4))     DEALLOCATE(pos4) 
      IF (ALLOCATED(pos5))     DEALLOCATE(pos5)
      IF (ALLOCATED(pos6))     DEALLOCATE(pos6) 
      IF (ALLOCATED(sti1))     DEALLOCATE(sti1)
      IF (ALLOCATED(sti2))     DEALLOCATE(sti2) 
      IF (ALLOCATED(sti3))     DEALLOCATE(sti3)
      IF (ALLOCATED(sti4))     DEALLOCATE(sti4) 
      IF (ALLOCATED(sti5))     DEALLOCATE(sti5)
      IF (ALLOCATED(sti6))     DEALLOCATE(sti6) 
      IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
c-----------
      END