rgwalp.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/.
!||====================================================================
!||    rgwalp        ../engine/source/constraints/general/rwall/rgwalp.F
!||--- called by ------------------------------------------------------
!||    rgwal0        ../engine/source/constraints/general/rwall/rgwal0.F
!||    rgwal0_imp    ../engine/source/constraints/general/rwall/rgwal0.F
!||--- calls      -----------------------------------------------------
!||    my_barrier    ../engine/source/system/machine.F
!||    sum_6_float   ../engine/source/system/parit.F
!||====================================================================
      SUBROUTINE rgwalp(X    ,A    ,V    ,RWL  ,NSW   ,
     1                  NSN  ,ITIED,MSR  ,MS   ,WEIGHT,
     2                  ICONT,FRWL6,IMP_S,NT_RW,IDDL  ,
     3                  IKC  ,NDOF ,NODNX_SMS ,WEIGHT_MD,WFEXT,WFEXT_MD)
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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com06_c.inc"
#include      "com08_c.inc"
#include      "scr05_c.inc"
#include      "scr11_c.inc"
#include      "impl1_c.inc"
#include      "sms_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NSN, ITIED, MSR,ICONT,IMP_S,NT_RW
      INTEGER NSW(*), WEIGHT(*), IDDL(*),IKC(*),NDOF(*), NODNX_SMS(*), WEIGHT_MD(*)
      my_real X(*), A(*), V(*), RWL(*), MS(*)
      DOUBLE PRECISION FRWL6(7,6)
      DOUBLE PRECISION,INTENT(INOUT) :: WFEXT,WFEXT_MD
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER M3, M2, M1, I, N, N3, N2, N1, K, J,JJ,NINDEX, INDEX(NSN)
      my_real
     .   XWL, YWL, ZWL, VXW, VYW, VZW, TFXTN, FACT,
     .   tfxt, xl1, yl1, zl1, xl2, yl2, zl2, sx12, sy12, sz12, s12,
     .   vx, vy, vz, xc, yc, zc, dp, xcp, ycp, zcp,
     .   sx1m, sy1m, sz1m, ps, sm1, sxm2, sym2, szm2, sm2, dv, da, dvt,
     .   fnxn, fnyn, fnzn, fnxt, fnyt, fnzt, fndfn, ftdft, fric, fric2,
     .   fcoe,ax,
     .   xwl0, ywl0, zwl0, dp0, xc0, yc0, zc0, tol, vn, vnold, 
     .   dp0dt, dvx, dvy, dvz, prec, xprec,
     .   f1(nsn), f2(nsn), f3(nsn), f4(nsn), f5(nsn), f6(nsn), f7(nsn),
     .   tfxt2, tfxtn2, wewe2
      DOUBLE PRECISION 
     .   FRWL6_L(7,6)
C-----------------------------------------------
      prec=em07
      ax = one
      icont=0
      CALL my_barrier
!       We need an OMP barrier here because each thread initializes 
!       the variables : without barrier ICONT can be set to 1 by a 
!       thread in the !$OMP DO loop whereas another (late) thread initializes 
!       ICONT to 0
!       ICONT is a shared variable 
      IF(msr==0)THEN
       xwl=rwl(4)
       ywl=rwl(5)
       zwl=rwl(6)
       vxw=zero
       vyw=zero
       vzw=zero
       xwl0=xwl
       ywl0=ywl
       zwl0=zwl
       vnold=zero
       vn   =zero
      ELSE
       m3=3*msr
       m2=m3-1
       m1=m2-1
C   changement formulation : plus d'impasse sur contribution force
       vxw=v(m1)+a(m1)*dt12
       vyw=v(m2)+a(m2)*dt12
       vzw=v(m3)+a(m3)*dt12
       xwl=x(m1)+vxw*dt2
       ywl=x(m2)+vyw*dt2
       zwl=x(m3)+vzw*dt2
       xwl0=x(m1)
       ywl0=x(m2)
       zwl0=x(m3)
       vnold =rwl(4)
       vn    =vxw*rwl(1)+vyw*rwl(2)+vzw*rwl(3)
       rwl(4)=vn
      ENDIF
      tfxt =zero    
      tfxtn=zero
      tfxt2 =zero    
      tfxtn2=zero    
C
      xl1=rwl(7)
      yl1=rwl(8)
      zl1=rwl(9)
      xl2=rwl(10)
      yl2=rwl(11)
      zl2=rwl(12)
      sx12=yl1*zl2-zl1*yl2
      sy12=zl1*xl2-xl1*zl2
      sz12=xl1*yl2-yl1*xl2
      s12=sx12**2+sy12**2+sz12**2
      nindex = 0
C----  
      IF(idtmins==0.AND.idtmins_int==0)THEN
!$OMP DO
       DO i=1,nsn
        n=nsw(i)
        n3=3*n
        n2=n3-1
        n1=n2-1
        IF(n2d==1)THEN
         ax=x(n2)
        ELSE
         ax=one
        ENDIF   
C
        xc0=x(n1)-xwl0
        yc0=x(n2)-ywl0
        zc0=x(n3)-zwl0
        dp0=xc0*rwl(1)+yc0*rwl(2)+zc0*rwl(3)
C
        vx=v(n1)+a(n1)*dt12
        vy=v(n2)+a(n2)*dt12
        vz=v(n3)+a(n3)*dt12
        xc=xc0+(vx-vxw)*dt2
        yc=yc0+(vy-vyw)*dt2
        zc=zc0+(vz-vzw)*dt2
        dp=dp0+((vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3))*dt2
C
        tol=two*dt1*max(
     .   abs( (vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3) ),
     .   abs(vn-vnold) )
        IF(iresp==1)THEN
          xprec=prec*max(abs(xwl),abs(ywl),abs(zwl),
     .              abs(x(n1)),abs(x(n2)),abs(x(n3)),
     .              abs(x(n1)-xwl),abs(x(n2)-ywl),abs(x(n3)-zwl))
          tol=max(tol,xprec)
        END IF
        IF(dp>zero.OR.dp0<=-tol) cycle
C
        xcp=xc-dp*rwl(1)
        ycp=yc-dp*rwl(2)
        zcp=zc-dp*rwl(3)
C
        sx1m=yl1*zcp-zl1*ycp
        sy1m=zl1*xcp-xl1*zcp
        sz1m=xl1*ycp-yl1*xcp
        ps=sx12*sx1m+sy12*sy1m+sz12*sz1m
C
        IF(ps<zero) cycle
C
        sm1=sx1m**2+sy1m**2+sz1m**2
C
        IF(sm1>s12) cycle
C
        sxm2=ycp*zl2-zcp*yl2
        sym2=zcp*xl2-xcp*zl2
        szm2=xcp*yl2-ycp*xl2
        ps=sx12*sxm2+sy12*sym2+sz12*szm2
C
        IF(ps<zero) cycle
C
        sm2=sxm2**2+sym2**2+szm2**2
C
        IF(sm2>s12) cycle
C
        icont=1
        IF((vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3)>zero
     .     .AND.dp0>zero)cycle
        nindex = nindex + 1
        index(nindex) = i
       END DO
!$OMP END DO       
      ELSE
!$OMP DO
       DO i=1,nsn
        n=nsw(i)
        IF(nodnx_sms(n)/=0)cycle
        n3=3*n
        n2=n3-1
        n1=n2-1
        IF(n2d==1)THEN
         ax=x(n2)
        ELSE
         ax=one
        ENDIF   
C
        xc0=x(n1)-xwl0
        yc0=x(n2)-ywl0
        zc0=x(n3)-zwl0
        dp0=xc0*rwl(1)+yc0*rwl(2)+zc0*rwl(3)
C
        vx=v(n1)+a(n1)*dt12
        vy=v(n2)+a(n2)*dt12
        vz=v(n3)+a(n3)*dt12
        xc=xc0+(vx-vxw)*dt2
        yc=yc0+(vy-vyw)*dt2
        zc=zc0+(vz-vzw)*dt2
        dp=dp0+((vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3))*dt2
C
        tol=two*dt1*max(
     .   abs( (vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3) ),
     .   abs(vn-vnold) )
        IF(iresp==1)THEN
          xprec=prec*max(abs(xwl),abs(ywl),abs(zwl),
     .              abs(x(n1)),abs(x(n2)),abs(x(n3)),
     .              abs(x(n1)-xwl),abs(x(n2)-ywl),abs(x(n3)-zwl))
          tol=max(tol,xprec)
        END IF
        IF(dp>zero.OR.dp0<=-tol)cycle
C
        xcp=xc-dp*rwl(1)
        ycp=yc-dp*rwl(2)
        zcp=zc-dp*rwl(3)
C
        sx1m=yl1*zcp-zl1*ycp
        sy1m=zl1*xcp-xl1*zcp
        sz1m=xl1*ycp-yl1*xcp
        ps=sx12*sx1m+sy12*sy1m+sz12*sz1m
C
        IF(ps<zero) cycle
C
        sm1=sx1m**2+sy1m**2+sz1m**2
C
        IF(sm1>s12) cycle
C
        sxm2=ycp*zl2-zcp*yl2
        sym2=zcp*xl2-xcp*zl2
        szm2=xcp*yl2-ycp*xl2
        ps=sx12*sxm2+sy12*sym2+sz12*szm2
C
        IF(ps<zero) cycle
C
        sm2=sxm2**2+sym2**2+szm2**2
C
        IF(sm2>s12) cycle
C
        icont=1
        IF((vx-vxw)*rwl(1)+(vy-vyw)*rwl(2)+(vz-vzw)*rwl(3)>zero
     .     .AND.dp0>zero)cycle
        nindex = nindex + 1
        index(nindex) = i
       END DO
!$OMP END DO
      END IF
C----  
      fact=one/dt12
      IF(itied/=0)THEN
       DO 40 j = 1, nindex
       i = index(j)
       n = nsw(i)
       n3=3*n
       n2=n3-1
       n1=n2-1
C
       xc0=x(n1)-xwl0
       yc0=x(n2)-ywl0
       zc0=x(n3)-zwl0
       dp0  =xc0*rwl(1)+yc0*rwl(2)+zc0*rwl(3)
       dp0dt=-min(dp0,zero)/dt2
       dvx  =dp0dt*rwl(1)
       dvy  =dp0dt*rwl(2)
       dvz  =dp0dt*rwl(3)
C
       dv=(v(n1)-vxw)*rwl(1)+(v(n2)-vyw)*rwl(2)+(v(n3)-vzw)*rwl(3)
       da=a(n1)*rwl(1)+a(n2)*rwl(2)+a(n3)*rwl(3)
       dvt=dv+da*dt12-dp0dt
       fnxn=dvt*rwl(1)*ms(n)
       fnyn=dvt*rwl(2)*ms(n)
       fnzn=dvt*rwl(3)*ms(n)
       tfxtn = tfxtn - weight_md(n)*fact*
     . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)
       wewe2 = (1-weight_md(n))*weight(n)     
       tfxtn2 = tfxtn2 - wewe2*fact*
     . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)     
       f1(j) = fnxn*weight_md(n)
       f2(j) = fnyn*weight_md(n)*ax
       f3(j) = fnzn*weight_md(n)*ax
       f4(j) = ms(n)*weight_md(n)
       fnxt=((v(n1)-vxw)+a(n1)*dt12-dvx)*ms(n)-fnxn
       fnyt=((v(n2)-vyw)+a(n2)*dt12-dvy)*ms(n)-fnyn
       fnzt=((v(n3)-vzw)+a(n3)*dt12-dvz)*ms(n)-fnzn
       fndfn=fnxn**2+fnyn**2+fnzn**2
       ftdft=fnxt**2+fnyt**2+fnzt**2
       fric=rwl(13)
       fric2=fric**2
       IF(ftdft<=fric2*fndfn.OR.itied==1) THEN
C        POINT SECND TIED
         a(n1)=zero
         a(n2)=zero
         a(n3)=zero 
         v(n1)=vxw+dvx
         v(n2)=vyw+dvy
         v(n3)=vzw+dvz
         IF (imp_s==1)THEN
              IF(ndof(n)>0) THEN
                jj=iddl(n)+1
                IF (ikc(jj)==0)ikc(jj)=3
                IF (ikc(jj+1)==0)ikc(jj+1)=3
                IF (ikc(jj+2)==0)ikc(jj+2)=3
              ENDIF
         ENDIF
       ELSE
C        POINT SECND SLIDING
         fcoe=fric*sqrt(fndfn/max(em20,ftdft))
         fnxt=fcoe*fnxt
         fnyt=fcoe*fnyt
         fnzt=fcoe*fnzt
         a(n1)=a(n1)-(da*rwl(1)+fnxt/(dt12*ms(n)))
         a(n2)=a(n2)-(da*rwl(2)+fnyt/(dt12*ms(n)))
         a(n3)=a(n3)-(da*rwl(3)+fnzt/(dt12*ms(n)))
         v(n1)=v(n1)-dv*rwl(1)+dvx
         v(n2)=v(n2)-dv*rwl(2)+dvy
         v(n3)=v(n3)-dv*rwl(3)+dvz
          IF (imp_s==1)THEN
            IF (ndof(n)>0) THEN
              v(n1)=-dv
              a(n1)=rwl(1)
              a(n2)=rwl(2)
              a(n3)=rwl(3)    
              jj=iddl(n)+1
              IF (ikc(jj)==0)ikc(jj)=10
            ENDIF
          ENDIF
         tfxt=tfxt-
     .      ((v(n1)-vxw)*fnxt+(v(n2)-vyw)*fnyt+(v(n3)-vzw)*fnzt)
     .       *weight_md(n)
         wewe2 = (1-weight_md(n))*weight(n)
         tfxt2=tfxt2-
     .      ((v(n1)-vxw)*fnxt+(v(n2)-vyw)*fnyt+(v(n3)-vzw)*fnzt)
     .       *wewe2         
       ENDIF
       f5(j) = fnxt*weight_md(n)
       f6(j) = fnyt*weight_md(n)*ax
       f7(j) = fnzt*weight_md(n)*ax
 40    CONTINUE
      ELSE
       DO 60 j = 1, nindex
       i = index(j)
       n = nsw(i)
       n3=3*n
       n2=n3-1
       n1=n2-1
       xc0=x(n1)-xwl0
       yc0=x(n2)-ywl0
       zc0=x(n3)-zwl0
       dp0  =xc0*rwl(1)+yc0*rwl(2)+zc0*rwl(3)
       dp0dt=-min(dp0,zero)/dt2
       dvx  =dp0dt*rwl(1)
       dvy  =dp0dt*rwl(2)
       dvz  =dp0dt*rwl(3)
       dv=(v(n1)-vxw)*rwl(1)+(v(n2)-vyw)*rwl(2)+(v(n3)-vzw)*rwl(3)
       da=a(n1)*rwl(1)+a(n2)*rwl(2)+a(n3)*rwl(3)
       dvt=dv+da*dt12-dp0dt
       fnxn=dvt*rwl(1)*ms(n)
       fnyn=dvt*rwl(2)*ms(n)
       fnzn=dvt*rwl(3)*ms(n)
       tfxtn = tfxtn - weight_md(n)*fact*
     . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)
       wewe2 = (1-weight_md(n))*weight(n)     
       tfxtn2 = tfxtn2 - wewe2*fact*
     . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)     
       f1(j) = fnxn*weight_md(n)
       f2(j) = fnyn*weight_md(n)*ax
       f3(j) = fnzn*weight_md(n)*ax
       f4(j) = ms(n)*weight_md(n)
       f5(j) = zero
       f6(j) = zero
       f7(j) = zero
       a(n1)=a(n1)-da*rwl(1)
       a(n2)=a(n2)-da*rwl(2)
       a(n3)=a(n3)-da*rwl(3)
       v(n1)=v(n1)-dv*rwl(1)+dvx
       v(n2)=v(n2)-dv*rwl(2)+dvy
       v(n3)=v(n3)-dv*rwl(3)+dvz
        IF(imp_s==1) v(n1) = -dv+dp0dt
 60    CONTINUE
      ENDIF
C
      IF(imconv==1)THEN          
          tfxt=tfxt+half*dt1*tfxtn
          tfxt2=tfxt2+half*dt1*tfxtn2
!$omp atomic
          wfext=wfext+tfxt
          wfext_md=wfext_md+tfxt2
 
          DO k = 1, 6
            frwl6_l(1,k) = zero
            frwl6_l(2,k) = zero
            frwl6_l(3,k) = zero
            frwl6_l(4,k) = zero
            frwl6_l(5,k) = zero
            frwl6_l(6,k) = zero
            frwl6_l(7,k) = zero
          END DO
          CALL sum_6_float(1, nindex, f1, frwl6_l(1,1), 7)
          CALL sum_6_float(1, nindex, f2, frwl6_l(2,1), 7)
          CALL sum_6_float(1, nindex, f3, frwl6_l(3,1), 7)
          CALL sum_6_float(1, nindex, f4, frwl6_l(4,1), 7)
          CALL sum_6_float(1, nindex, f5, frwl6_l(5,1), 7)
          CALL sum_6_float(1, nindex, f6, frwl6_l(6,1), 7)
          CALL sum_6_float(1, nindex, f7, frwl6_l(7,1), 7)
          
#include "lockon.inc"
          DO k = 1, 6
            frwl6(1,k) = frwl6(1,k)+frwl6_l(1,k)
            frwl6(2,k) = frwl6(2,k)+frwl6_l(2,k)
            frwl6(3,k) = frwl6(3,k)+frwl6_l(3,k)
            frwl6(4,k) = frwl6(4,k)+frwl6_l(4,k)
            frwl6(5,k) = frwl6(5,k)+frwl6_l(5,k)
            frwl6(6,k) = frwl6(6,k)+frwl6_l(6,k)
            frwl6(7,k) = frwl6(7,k)+frwl6_l(7,k)
          END DO
#include "lockoff.inc"
      ENDIF
      
      IF (imp_s==1) THEN
        IF(itied==0)THEN
          DO j=1,nindex
           i = index(j)
           n=nsw(i)
           IF (ndof(n)>0) THEN
           n3=3*n
           n2=n3-1
           n1=n2-1
           a(n1)=rwl(1)
           a(n2)=rwl(2)
           a(n3)=rwl(3)    
           jj=iddl(n)+1
           IF (ikc(jj)==0)ikc(jj)=10
           ENDIF 
          ENDDO
        ENDIF
        DO j=1,nindex
          i = index(j)
          n=nsw(i)
          IF (ndof(n)>0) THEN
C to be uncommented the day it is parallel in implicit
!$OMP ATOMIC
            nt_rw = nt_rw + 1
          END IF
        ENDDO
      ENDIF
C
      RETURN
      END