rgwalc.F File Reference

#include "implicit_f.inc"
#include "comlock.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr11_c.inc"
#include "impl1_c.inc"
#include "sms_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	rgwalc (x, a, v, rwl, nsw, nsn, itied, msr, ms, weight, icont, frwl6, imp_s, nt_rw, iddl, ikc, ndof, nodnx_sms, weight_md, wfext, wfext_md)

Function/Subroutine Documentation

rgwalc()

subroutine rgwalc	(		x,
			a,
			v,
			rwl,
		integer, dimension(*)	nsw,
		integer	nsn,
		integer	itied,
		integer	msr,
			ms,
		integer, dimension(*)	weight,
		integer	icont,
		double precision, dimension(7,6)	frwl6,
		integer	imp_s,
		integer	nt_rw,
		integer, dimension(*)	iddl,
		integer, dimension(*)	ikc,
		integer, dimension(*)	ndof,
		integer, dimension(*)	nodnx_sms,
		integer, dimension(*)	weight_md,
		double precision, intent(inout)	wfext,
		double precision, intent(inout)	wfext_md )

Definition at line 32 of file rgwalc.F.

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      "com06_c.inc"
#include      "com08_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,J,JJ,K,NINDEX,
     .        INDEX(NSN)
C     REAL
      my_real
     .   ra2, xwl, ywl, zwl, vxw, vyw, vzw, wfxtn, fact,
     .   wfxt, vx, vy, vz, ux, uy, uz, xc, yc, zc, dd1, dd,
     .   dp, xt, yt, zt, xx, xn, yn, zn, dv, da, dvt, fnxn, fnyn, fnzn,
     .   fnxt, fnyt, fnzt, fndfn, ftdft, fric, fric2, fcoe,
     .   xwlo, ywlo, zwlo,
     .   f1(nsn), f2(nsn), f3(nsn), f4(nsn), f5(nsn), f6(nsn), f7(nsn),
     .   wfxt2, wfxtn2, wewe2
      DOUBLE PRECISION 
     .   FRWL6_L(7,6)
C-----------------------------------------------
      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 
      ra2=(half*rwl(7))**2
      IF(msr==0)THEN
       xwlo=rwl(4)
       ywlo=rwl(5)
       zwlo=rwl(6)
       xwl=rwl(4)
       ywl=rwl(5)
       zwl=rwl(6)
       vxw=zero
       vyw=zero
       vzw=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
       xwlo=x(m1)
       ywlo=x(m2)
       zwlo=x(m3)
       xwl=x(m1)+vxw*dt2
       ywl=x(m2)+vyw*dt2
       zwl=x(m3)+vzw*dt2
      ENDIF
      wfxt =zero
      wfxtn=zero
      wfxt2 =zero
      wfxtn2=zero
      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
        vx=v(n1)+a(n1)*dt12
        vy=v(n2)+a(n2)*dt12
        vz=v(n3)+a(n3)*dt12
        ux=x(n1)+vx*dt2
        uy=x(n2)+vy*dt2
        uz=x(n3)+vz*dt2
        xc=ux-xwl
        yc=uy-ywl
        zc=uz-zwl
        dd1=xc**2+yc**2+zc**2
        dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
        dp=dd1-dd**2
        IF(dp <= ra2)THEN
          icont=1
          nindex = nindex+1
          index(nindex) = i
        END IF
       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
        vx=v(n1)+a(n1)*dt12
        vy=v(n2)+a(n2)*dt12
        vz=v(n3)+a(n3)*dt12
        ux=x(n1)+vx*dt2
        uy=x(n2)+vy*dt2
        uz=x(n3)+vz*dt2
        xc=ux-xwl
        yc=uy-ywl
        zc=uz-zwl
        dd1=xc**2+yc**2+zc**2
        dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
        dp=dd1-dd**2
        IF(dp <= ra2)THEN
          icont=1
          nindex = nindex+1
          index(nindex) = i
        END IF
       END DO
!$OMP END DO
      END IF
C
      fact=one/dt12
      DO j = 1,nindex
       i = index(j)
       n=nsw(i)
       n3=3*n
       n2=n3-1
       n1=n2-1
       xc=x(n1)-xwlo
       yc=x(n2)-ywlo
       zc=x(n3)-zwlo
       dd1=xc**2+yc**2+zc**2
       dd =xc*rwl(1)+yc*rwl(2)+zc*rwl(3)
       dp=dd1-dd**2
       xt=dd*rwl(1)
       yt=dd*rwl(2)
       zt=dd*rwl(3)
       xx=sqrt(dp)
       xn=(xc-xt)/xx
       yn=(yc-yt)/xx
       zn=(zc-zt)/xx
       dv=(v(n1)-vxw)*xn+(v(n2)-vyw)*yn+(v(n3)-vzw)*zn
       da=a(n1)*xn+a(n2)*yn+a(n3)*zn
       dvt=dv+da*dt12
       fnxn=dvt*xn*ms(n)
       fnyn=dvt*yn*ms(n)
       fnzn=dvt*zn*ms(n)
       wfxtn = wfxtn - weight_md(n)*fact*
     . ((v(n1)-vxw)*fnxn+(v(n2)-vyw)*fnyn+(v(n3)-vzw)*fnzn)
       wewe2 = (1-weight_md(n))*weight(n)
       wfxtn2 = wfxtn2 - 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)
       f3(j) = fnzn*weight_md(n)
       f4(j) = ms(n)*weight_md(n)
       IF(itied/=0)THEN
        fnxt=((v(n1)-vxw)+a(n1)*dt12)*ms(n)-fnxn
        fnyt=((v(n2)-vyw)+a(n2)*dt12)*ms(n)-fnyn
        fnzt=((v(n3)-vzw)+a(n3)*dt12)*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
          v(n2)=vyw
          v(n3)=vzw
          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/ftdft)
          fnxt=fcoe*fnxt
          fnyt=fcoe*fnyt
          fnzt=fcoe*fnzt
          a(n1)=a(n1)-(da*xn+fnxt/(dt12*ms(n)))
          a(n2)=a(n2)-(da*yn+fnyt/(dt12*ms(n)))
          a(n3)=a(n3)-(da*zn+fnzt/(dt12*ms(n)))
          v(n1)=v(n1)-dv*xn
          v(n2)=v(n2)-dv*yn
          v(n3)=v(n3)-dv*zn
          IF (imp_s==1) THEN
            IF(ndof(n)>0) THEN
              v(n1)=-dv
              a(n1)=xn
              a(n2)=yn
              a(n3)=zn    
              jj=iddl(n)+1
              IF (ikc(jj)==0)ikc(jj)=10
            ENDIF 
          ENDIF
          wfxt=wfxt-
     .        ((v(n1)-vxw)*fnxt+(v(n2)-vyw)*fnyt+(v(n3)-vzw)*fnzt)
     .        *weight_md(n)
          wewe2 = (1-weight_md(n))*weight(n)
          wfxt2=wfxt2-
     .        ((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)
        f7(j) = fnzt*weight_md(n)
       ELSE
        f5(j) = zero
        f6(j) = zero
        f7(j) = zero
        a(n1)=a(n1)-da*xn
        a(n2)=a(n2)-da*yn
        a(n3)=a(n3)-da*zn
        v(n1)=v(n1)-dv*xn
        v(n2)=v(n2)-dv*yn
        v(n3)=v(n3)-dv*zn
        IF (imp_s==1) THEN
          IF(ndof(n)>0) THEN
            v(n1)=-dv
            a(n1)=xn
            a(n2)=yn
            a(n3)=zn    
            jj=iddl(n)+1
            IF (ikc(jj)==0)ikc(jj)=10
          ENDIF
        ENDIF 
       ENDIF
      ENDDO 
C
      IF(imconv==1)THEN
        wfxt=wfxt+half*dt1*wfxtn
        wfxt2=wfxt2+half*dt1*wfxtn2
!$OMP ATOMIC
        wfext=wfext+wfxt
        wfext_md=wfext_md+wfxt2
        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
        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