scfint_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/.
!||====================================================================
!||    scfint_reg       ../engine/source/elements/thickshell/solidec/scfint_reg.F
!||--- called by ------------------------------------------------------
!||    scforc3          ../engine/source/elements/thickshell/solidec/scforc3.f
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod     ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    nlocal_reg_mod   ../common_source/modules/nlocal_reg_mod.F
!||====================================================================
      SUBROUTINE scfint_reg(
     1   NLOC_DMG ,VAR_REG  ,NEL      ,OFF      ,
     2   VOL      ,NC1      ,NC2      ,NC3      ,
     3   NC4      ,NC5      ,NC6      ,NC7      ,
     4   NC8      ,PX1      ,PX2      ,PX3      ,
     5   PX4      ,PY1      ,PY2      ,PY3      ,
     6   PY4      ,PZ1      ,PZ2      ,PZ3      ,
     7   PZ4      ,IMAT     ,ITASK    ,DT2T     ,
     8   VOL0     ,NFT      ,NLAY     ,WS       ,
     9   AS       ,AREA     ,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-----------------------------------------------
      TYPE(nlocal_str_), INTENT(INOUT), TARGET :: NLOC_DMG !< Non-local data structure
      my_real, DIMENSION(NEL,NLAY), INTENT(INOUT) :: VAR_REG !< Variable to regularize
      INTEGER, INTENT(IN)    :: NEL               !< Number of elements in the group
      my_real, DIMENSION(NEL), INTENT(IN) :: OFF  !< Flag for element deletion
      my_real, DIMENSION(NEL), INTENT(IN) :: VOL  !< Current volume of the element
      INTEGER, DIMENSION(NEL) :: NC1              !< Node 1 of the elements
      INTEGER, DIMENSION(NEL) :: NC2              !< Node 2 of the elements
      INTEGER, DIMENSION(NEL) :: NC3              !< Node 3 of the elements
      INTEGER, DIMENSION(NEL) :: NC4              !< Node 4 of the elements
      INTEGER, DIMENSION(NEL) :: NC5              !< Node 5 of the elements
      INTEGER, DIMENSION(NEL) :: NC6              !< Node 6 of the elements
      INTEGER, DIMENSION(NEL) :: NC7              !< Node 7 of the elements
      INTEGER, DIMENSION(NEL) :: NC8              !< Node 8 of the elements
      my_real, DIMENSION(NEL), INTENT(IN) :: px1  !< Derivative of the first node shape function along x
      my_real, DIMENSION(NEL), INTENT(IN) :: px2  !< Derivative of the second node shape function along x
      my_real, DIMENSION(NEL), INTENT(IN) :: px3  !< Derivative of the third node shape function along x
      my_real, DIMENSION(NEL), INTENT(IN) :: px4  !< Derivative of the fourth node shape function along x
      my_real, DIMENSION(NEL), INTENT(IN) :: py1  !< Derivative of the first node shape function along y
      my_real, DIMENSION(NEL), INTENT(IN) :: py2  !< Derivative of the second node shape function along y
      my_real, DIMENSION(NEL), INTENT(IN) :: py3  !< Derivative of the third node shape function along y
      my_real, DIMENSION(NEL), INTENT(IN) :: py4  !< Derivative of the fourth node shape function along y
      my_real, DIMENSION(NEL), INTENT(IN) :: pz1  !< Derivative of the first node shape function along z
      my_real, DIMENSION(NEL), INTENT(IN) :: pz2  !< Derivative of the second node shape function along z
      my_real, DIMENSION(NEL), INTENT(IN) :: pz3  !< Derivative of the third node shape function along z
      my_real, DIMENSION(NEL), INTENT(IN) :: pz4  !< Derivative of the fourth node shape function along z
      INTEGER, INTENT(IN)    :: IMAT              !< Material internal Id
      INTEGER, INTENT(IN)    :: ITASK             !< Number of the thread
      my_real, INTENT(INOUT) :: dt2t              !< Element critical timestep
      my_real, DIMENSION(NEL), INTENT(IN) :: vol0 !< Initial volume of the element
      INTEGER, INTENT(IN)    :: NFT               !< Address of the first element of the group
      INTEGER, INTENT(IN)    :: NLAY              !< Number of layers in the elements
      my_real, DIMENSION(9,9), INTENT(IN) :: ws   !< Weight of integration points through thickness
      my_real, DIMENSION(9,9), INTENT(IN) :: as   !< Position of integration points through thickness
      my_real, DIMENSION(NEL), INTENT(IN) :: area !< Area of the thickshell element
      TYPE(buf_nlocts_), INTENT(INOUT), TARGET :: BUFNLTS !< Non-local specific buffer for thickshell thickness
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, II, J, K, N1, N2, N3, N4, N5, N6, N7, N8, 
     . l_nloc,nddl,ndnod
      INTEGER :: NODA_DT
      INTEGER, DIMENSION(:), ALLOCATABLE   ::
     .   pos1,pos2,pos3,pos4,pos5,pos6,pos7,pos8
      my_real 
     .   l2,ntn_unl,ntn_vnl,xi,ntvar,a,dtnl,le_max,
     .   b1,b2,b3,b4,b5,b6,b7,b8,zeta,sspnl,maxstif,
     .   bth1,bth2,nth1,nth2,dt2p,dtnod,k1,k2,k12,
     .   dtnl_th,ntn
      my_real, DIMENSION(:) ,ALLOCATABLE   :: 
     .   btb11,btb12,btb13,btb14,btb22,btb23,btb24,
     .   btb33,btb34,btb44,lc,thk,lthk
      my_real, DIMENSION(:,:) ,ALLOCATABLE ::
     .   f1,f2,f3,f4,f5,f6,f7,f8,stifnlth,dtn,
     .   sti1,sti2,sti3,sti4,sti5,sti6,sti7,sti8
      my_real, POINTER, DIMENSION(:) :: 
     . vnl,fnl,unl,mass,mass0,vnl0
      my_real, POINTER, DIMENSION(:,:) :: 
     .   massth,unlth,vnlth,fnlth
      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=======================================================================
      noda_dt = nodadt
c 
      ! Recover non-local parameters
      l2     = nloc_dmg%LEN(imat)**2 !< Squared internal length 
      xi     = nloc_dmg%DAMP(imat)   !< Damping parameter
      l_nloc = nloc_dmg%L_NLOC       !< Size of non-local tables
      zeta   = nloc_dmg%DENS(imat)   !< Density parameter
      sspnl  = nloc_dmg%SSPNL(imat)  !< Non-local "Sound speed"
      le_max = nloc_dmg%LE_MAX(imat) !< Maximal length of convergence
c
      ! Allocation of variable table
      ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb22(nel),
     .         btb23(nel),btb24(nel),btb33(nel),btb34(nel),btb44(nel),
     .         f1(nel,nlay),f2(nel,nlay),f3(nel,nlay),f4(nel,nlay),
     .         f5(nel,nlay),f6(nel,nlay),f7(nel,nlay),f8(nel,nlay),
     .         pos1(nel),pos2(nel),pos3(nel),pos4(nel),pos5(nel),
     .         pos6(nel),pos7(nel),pos8(nel),lc(nel),thk(nel),lthk(nel))
      ! Factor for non-local forces
      ntn = eight*eight
      ! Local variables initialization
      lc(1:nel) = zero              
      ! For nodal timestep
      IF (noda_dt > 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),sti7(nel,nlay),sti8(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,mass0)
        ALLOCATE(sti1(1,1),sti2(1,1),sti3(1,1),sti4(1,1),
     .          sti5(1,1),sti6(1,1),sti7(1,1),sti8(1,1))
      ENDIF
      ! Current timestep non-local velocities
      vnl  => nloc_dmg%VNL(1:l_nloc)
      ! Previous timestep non-local velocities
      vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
      ! Non-local displacements
      unl  => nloc_dmg%UNL(1:l_nloc)
c
      !--------------------------------------------------------------------------------
      ! Computation of the position of the non-local d.o.fs and the BtB matrix product
      !--------------------------------------------------------------------------------
      ! Loop over elements
      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))
        n7 = nloc_dmg%IDXI(nc7(i))
        n8 = nloc_dmg%IDXI(nc8(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)
        pos7(i) = nloc_dmg%POSI(n7)
        pos8(i) = nloc_dmg%POSI(n8)
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
      !-----------------------------------------------------------------------
      ! 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 (noda_dt > 0) THEN 
          ALLOCATE(stifnlth(nel,nddl+1))
          ALLOCATE(dtn(nel,nddl+1))
        ELSE
          ALLOCATE(dtn(1,1))
          dtn = ep20
          ALLOCATE(stifnlth(1,1))
          stifnlth = 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 (noda_dt > 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 (noda_dt > 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 (noda_dt > 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
              dtnod = dtmin1(11)*(sqrt(csta))
            ENDIF
          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 layers
      DO k = 1,nlay
c
        ! Loop over elements
        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) + unl(pos7(i)+k-1) + unl(pos8(i)+k-1)) / ntn
c        
            ! Computing the product 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) + vnl(pos7(i)+k-1) + vnl(pos8(i)+k-1)) / ntn
            IF (noda_dt > 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)), 
     .                      sqrt(mass(pos7(i)+k-1)/mass0(pos7(i)+k-1)), 
     .                      sqrt(mass(pos8(i)+k-1)/mass0(pos8(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) - btb13(i)*unl(pos5(i)+k-1)
     .                - btb14(i)*unl(pos6(i)+k-1) - btb11(i)*unl(pos7(i)+k-1) - btb12(i)*unl(pos8(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) - btb23(i)*unl(pos5(i)+k-1)
     .                - btb24(i)*unl(pos6(i)+k-1) - btb12(i)*unl(pos7(i)+k-1) - btb22(i)*unl(pos8(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) - btb33(i)*unl(pos5(i)+k-1)
     .                - btb34(i)*unl(pos6(i)+k-1) - btb13(i)*unl(pos7(i)+k-1) - btb23(i)*unl(pos8(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) - btb34(i)*unl(pos5(i)+k-1)
     .                - btb44(i)*unl(pos6(i)+k-1) - btb14(i)*unl(pos7(i)+k-1) - btb24(i)*unl(pos8(i)+k-1))
c       
            b5 = 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) + btb33(i)*unl(pos5(i)+k-1)
     .                + btb34(i)*unl(pos6(i)+k-1) + btb13(i)*unl(pos7(i)+k-1) + btb23(i)*unl(pos8(i)+k-1))
c        
            b6 = 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) + btb34(i)*unl(pos5(i)+k-1)
     .                + btb44(i)*unl(pos6(i)+k-1) + btb14(i)*unl(pos7(i)+k-1) + btb24(i)*unl(pos8(i)+k-1))
c        
            b7 = 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) + btb13(i)*unl(pos5(i)+k-1)
     .                + btb14(i)*unl(pos6(i)+k-1) + btb11(i)*unl(pos7(i)+k-1) + btb12(i)*unl(pos8(i)+k-1))
c        
            b8 = 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) + btb23(i)*unl(pos5(i)+k-1)
     .                + btb24(i)*unl(pos6(i)+k-1) + btb12(i)*unl(pos7(i)+k-1) + btb22(i)*unl(pos8(i)+k-1))       
c
            ! Multiplication by the volume of the element (and damping parameter XI)  
            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_8* 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
            f7(i,k) = a + b7
            f8(i,k) = a + b8
c
            ! Computing nodal equivalent stiffness
            IF (noda_dt > 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*btb13(i) + one/ntn) + abs(-l2*btb14(i) + one/ntn) +
     .                     abs(-l2*btb11(i) + one/ntn) + abs(-l2*btb12(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*btb23(i) + one/ntn) + abs(-l2*btb24(i) + one/ntn) +
     .                     abs(-l2*btb12(i) + one/ntn) + abs(-l2*btb22(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*btb33(i) + one/ntn) + abs(-l2*btb34(i) + one/ntn) +
     .                     abs(-l2*btb13(i) + one/ntn) + abs(-l2*btb23(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*btb34(i) + one/ntn) + abs(-l2*btb44(i) + one/ntn) +
     .                     abs(-l2*btb14(i) + one/ntn) + abs(-l2*btb24(i) + one/ntn))*vol(i)*ws(k,nlay)*half
              sti5(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*btb33(i)  + one/ntn) + abs( l2*btb34(i) + one/ntn) +
     .                     abs(l2*btb13(i)  + one/ntn) + abs(l2*btb23(i)  + one/ntn))*vol(i)*ws(k,nlay)*half
              sti6(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*btb34(i)  + one/ntn) + abs( l2*btb44(i) + one/ntn)  +
     .                     abs(l2*btb14(i)  + one/ntn) + abs(l2*btb24(i)  + one/ntn))*vol(i)*ws(k,nlay)*half
              sti7(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*btb13(i)  + one/ntn) + abs( l2*btb14(i) + one/ntn)  +
     .                     abs(l2*btb11(i)  + one/ntn) + abs(l2*btb12(i)  + one/ntn))*vol(i)*ws(k,nlay)*half    
              sti8(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*btb23(i)  + one/ntn) + abs( l2*btb24(i) + one/ntn)  +
     .                     abs(l2*btb12(i)  + one/ntn) + abs(l2*btb22(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 (noda_dt > 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))*(three/four)*(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))*(three/four)*(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))*(three/four)*(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))*(three/four)*(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))*(three/four)*(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))*(three/four)*(lc(i)**2)
              f7(i,k) = sqrt(mass(pos7(i)+k-1)/mass0(pos7(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos7(i)+k-1)+vnl0(pos7(i)+k-1))*(three/four)*(lc(i)**2)
              f8(i,k) = sqrt(mass(pos8(i)+k-1)/mass0(pos8(i)+k-1))*zeta*sspnl*half*
     .                       (vnl(pos8(i)+k-1)+vnl0(pos8(i)+k-1))*(three/four)*(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
              sti7(i,k) = em20
              sti8(i,k) = em20
            ELSE
              ! Non-local absorbing forces
              f1(i,k) = zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*(three/four)*(lc(i)**2)
              f2(i,k) = zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*(three/four)*(lc(i)**2)
              f3(i,k) = zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*(three/four)*(lc(i)**2)
              f4(i,k) = zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*(three/four)*(lc(i)**2)
              f5(i,k) = zeta*sspnl*half*(vnl(pos5(i)+k-1)+vnl0(pos5(i)+k-1))*(three/four)*(lc(i)**2)
              f6(i,k) = zeta*sspnl*half*(vnl(pos6(i)+k-1)+vnl0(pos6(i)+k-1))*(three/four)*(lc(i)**2)
              f7(i,k) = zeta*sspnl*half*(vnl(pos7(i)+k-1)+vnl0(pos7(i)+k-1))*(three/four)*(lc(i)**2)
              f8(i,k) = zeta*sspnl*half*(vnl(pos8(i)+k-1)+vnl0(pos8(i)+k-1))*(three/four)*(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
c       IF (NODA_DT > 0) STIFNL => NLOC_DMG%STIFNL(1:L_NLOC,ITASK+1) ! Non-local equivalent nodal stiffness
c        
        ! Loop over non-local degrees of freedom
        DO k=1,nlay
          ! Loop over elements
          DO i=1,nel
            ! 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)
            fnl(pos7(i)+k-1) = fnl(pos7(i)+k-1) - f7(i,k)
            fnl(pos8(i)+k-1) = fnl(pos8(i)+k-1) - f8(i,k)     
            IF (noda_dt > 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),sti7(i,k),sti8(i,k))
              ! Computing nodal stiffness
              nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos1(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos2(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos3(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos4(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos5(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos6(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos7(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos7(i)+k-1,itask+1) + maxstif
              nloc_dmg%STIFNL(pos8(i)+k-1,itask+1) = nloc_dmg%STIFNL(pos8(i)+k-1,itask+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 (noda_dt > 0) THEN
              maxstif = max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j),
     .                      sti5(i,j),sti6(i,j),sti7(i,j),sti8(i,j))
            ENDIF
c            
            k = nloc_dmg%IADS(1,ii)
            nloc_dmg%FSKY(k,j) = -f1(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(2,ii)
            nloc_dmg%FSKY(k,j) = -f2(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(3,ii)
            nloc_dmg%FSKY(k,j) = -f3(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(4,ii)
            nloc_dmg%FSKY(k,j) = -f4(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(5,ii)
            nloc_dmg%FSKY(k,j) = -f5(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(6,ii)
            nloc_dmg%FSKY(k,j) = -f6(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(7,ii)
            nloc_dmg%FSKY(k,j) = -f7(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
            k = nloc_dmg%IADS(8,ii)
            nloc_dmg%FSKY(k,j) = -f8(i,j)
            IF (noda_dt > 0) nloc_dmg%STSKY(k,j) = maxstif
c
          ENDDO
        ENDDO
      ENDIF
c      
      !-----------------------------------------------------------------------
      ! Computing non-local timestep
      !-----------------------------------------------------------------------
      IF (noda_dt == 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 ((l2>zero).AND.(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
            ! Keeping 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(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(pos5))     DEALLOCATE(pos5)
      IF (ALLOCATED(pos6))     DEALLOCATE(pos6) 
      IF (ALLOCATED(pos7))     DEALLOCATE(pos7)
      IF (ALLOCATED(pos8))     DEALLOCATE(pos8)
      IF (ALLOCATED(f1))       DEALLOCATE(f1)
      IF (ALLOCATED(f2))       DEALLOCATE(f2) 
      IF (ALLOCATED(f3))       DEALLOCATE(f3)
      IF (ALLOCATED(f4))       DEALLOCATE(f4) 
      IF (ALLOCATED(f5))       DEALLOCATE(f5)
      IF (ALLOCATED(f6))       DEALLOCATE(f6) 
      IF (ALLOCATED(f7))       DEALLOCATE(f7)
      IF (ALLOCATED(f8))       DEALLOCATE(f8)
      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(sti7))     DEALLOCATE(sti7)
      IF (ALLOCATED(sti8))     DEALLOCATE(sti8)
      IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
      IF (ALLOCATED(dtn))      DEALLOCATE(dtn)
      IF (ALLOCATED(lc))       DEALLOCATE(lc)
      IF (ALLOCATED(thk))      DEALLOCATE(thk)
      IF (ALLOCATED(lthk))     DEALLOCATE(lthk)
c-----------
      END