rgwals.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	rgwals (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

rgwals()

subroutine rgwals	(		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 rgwals.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
      DOUBLE PRECISION,INTENT(INOUT) :: 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, K, JJ,NINDEX, INDEX(NSN)
      my_real
     .   ra2, xwl, ywl, zwl, vxw, vyw, vzw, tfxtn, fact,
     .   tfxt, vx, vy, vz, ux, uy, uz, xc, yc, zc, dp, 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),
     .   tfxt2, tfxtn2, 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
      tfxt =zero    
      tfxtn=zero
      tfxt2 =zero    
      tfxtn2=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
        dp=xc**2+yc**2+zc**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
        dp=xc**2+yc**2+zc**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
       xx=sqrt(xc**2+yc**2+zc**2)
       xn=xc/xx
       yn=yc/xx
       zn=zc/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)
       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)
       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
          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)
        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
          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
         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