rgbodv.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/.
!||====================================================================
!||    rgbodv                ../engine/source/constraints/general/rbody/rgbodv.F
!||--- called by ------------------------------------------------------
!||    rbyvit                ../engine/source/constraints/general/rbody/rbyvit.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                ../engine/source/output/message/message.F
!||    arret                 ../engine/source/system/arret.F
!||    velrot_explicit       ../engine/source/constraints/general/rbody/velrot_explicit.F90
!||--- uses       -----------------------------------------------------
!||    message_mod           ../engine/share/message_module/message_mod.F
!||    velrot_explicit_mod   ../engine/source/constraints/general/rbody/velrot_explicit.F90
!||====================================================================
      SUBROUTINE rgbodv(V     ,VR   ,X     ,RBY ,NOD  ,
     2                  NBY   ,SKEW ,ISKEW ,FS  ,ITAB ,
     3                  WEIGHT,A    ,AR    ,MS  ,IN   ,
     4                  ISENS ,ID   ,IFAIL ,FNY ,EXPN ,
     5                  FTY   ,EXPT ,CRIT  ,NODREAC,FTHREAC ,
     6                  FREAC )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE velrot_explicit_mod, ONLY : velrot_explicit
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NOD(*), NBY(*), ISKEW(*),ITAB(*), WEIGHT(*),ID, NODREAC(*)
      INTEGER, INTENT(IN) :: IFAIL
C     REAL
      my_real
     .   V(3,*), VR(3,*), X(3,*), RBY(*), SKEW(LSKEW,*), FS(*),
     .   a(3,*),ar(3,*),in(*),ms(*),freac(6,*)
      my_real, INTENT(IN) :: 
     .   fny,fty,expn,expt
      my_real, INTENT(INOUT) :: 
     .   crit, fthreac(6,*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "param_c.inc"
#include      "parit_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER M, NSN, ICODR, ISK, I, N,ISENS, J, L,NSN_G
C     REAL
      my_real
     .   XDUM(3), VG(3), VI(3),
     .   V1X2, V2X1, V2X3, V3X2, V3X1, V1X3,USDT,DT05,
     .   vx1, vx2, vx3, mas,am1,am2,am3,vxx1,vxx2,vxx3,vy1
      my_real
     .   ux, uy, uz, nn, fn, ft,lsm(3),vs(3),as(3)
      my_real, DIMENSION(:,:), ALLOCATABLE :: 
     .   r, rr
C-----------------------------------------------
      m    =nby(1)
C optimisation spmd
      IF (m<0) RETURN
      nsn  =nby(2)
      nsn_g = nby(19)
C en spmd il ne faut sommer qu'une fois la vitesse du centre d'inertie
      fs(7)=fs(7)+vr(1,m)*dt2*weight(m)
      fs(8)=fs(8)+vr(2,m)*dt2*weight(m)
      fs(9)=fs(9)+vr(3,m)*dt2*weight(m)
C
      IF(dt2*dt2*(vr(1,m)**2+vr(2,m)**2+vr(3,m)**2)>1.0.AND.nsn_g>2)THEN 
         CALL ancmsg(msgid=110,anmode=aninfo,
     .            i1=id)
         CALL arret(2)
         RETURN
      ENDIF
C
      IF(n2d==0) THEN
        vg(1)=vr(1,m)+ar(1,m)*dt12
        vg(2)=vr(2,m)+ar(2,m)*dt12
        vg(3)=vr(3,m)+ar(3,m)*dt12
C
        usdt = one/dt12
C
        am1 = a(1,m)
        am2 = a(2,m)
        am3 = a(3,m)
       IF (nsn_g<=2) THEN 
          DO i=1,nsn
            n = nod(i)
            ar(1:3,n)= (vg(1:3)-vr(1:3,n)) * usdt
            a(1:3,n)=  a(1:3,m)+(v(1:3,m)-v(1:3,n))*usdt
C
            lsm(1:3) = x(1:3,n)-x(1:3,m)
            CALL velrot_explicit(vg, lsm,vs,dt12)
C  
            vx1=vs(1)
            vx2=vs(2)
            vx3=vs(3)
 
            a(1,n)= a(1,n) +
     .         (vx1+half*dt2*(vg(2)*vx3-vg(3)*vx2))*usdt
            a(2,n)= a(2,n) +
     .         (vx2+half*dt2*(vg(3)*vx1-vg(1)*vx3))*usdt
            a(3,n)= a(3,n) +
     .         (vx3+half*dt2*(vg(1)*vx2-vg(2)*vx1))*usdt
          ENDDO
       ELSE
          DO i=1,nsn
            n = nod(i)
            ar(1,n)= (vg(1)-vr(1,n)) * usdt
            ar(2,n)= (vg(2)-vr(2,n)) * usdt
            ar(3,n)= (vg(3)-vr(3,n)) * usdt
C
            v1x2=vg(1)*(x(2,n)-x(2,m))
            v2x1=vg(2)*(x(1,n)-x(1,m))
            v2x3=vg(2)*(x(3,n)-x(3,m))
            v3x2=vg(3)*(x(2,n)-x(2,m))
            v3x1=vg(3)*(x(1,n)-x(1,m))
            v1x3=vg(1)*(x(3,n)-x(3,m))
C  
            vx1=v2x3-v3x2
            vx2=v3x1-v1x3
            vx3=v1x2-v2x1
 
            a(1,n)= am1 +
     .        (v(1,m)+vx1+half*dt2*(vg(2)*vx3-vg(3)*vx2)-v(1,n))*usdt
            a(2,n)= am2 +
     .        (v(2,m)+vx2+half*dt2*(vg(3)*vx1-vg(1)*vx3)-v(2,n))*usdt
            a(3,n)= am3 +
     .        (v(3,m)+vx3+half*dt2*(vg(1)*vx2-vg(2)*vx1)-v(3,n))*usdt
          ENDDO
       END IF ! (NSN<=2) THEN 
      ELSEIF(n2d ==1) THEN
        vg(1)=zero
        vg(2)=zero
        vg(3)=vr(3,m)+ar(3,m)*dt12
C
        usdt = one/dt12
C
        am1 = zero
        am2 = a(2,m)
        am3 = a(3,m)
         DO i=1,nsn
           n = nod(i)
           ar(1,n)= zero
           ar(2,n)= zero
           ar(3,n)= (vg(3)-vr(3,n)) * usdt
C  
           vx1= vg(3)*(x(1,n)-x(1,m))
           vy1=-vg(3)*(x(2,n)-x(2,m))
 
           a(1,n)= am1 + (v(1,m)+vy1-half*dt2*vg(3)*vx1-v(1,n))*usdt
           a(2,n)= am2 + (v(2,m)+vx1+half*dt2*vg(3)*vy1-v(2,n))*usdt
           a(3,n)= am3 + (v(3,m)-v(3,n))*usdt
         ENDDO
      ELSEIF(n2d ==2) THEN
        vg(1)=vr(1,m)+ar(1,m)*dt12
        vg(2)=zero
        vg(3)=zero
C
        usdt = one/dt12
C
        am1 = zero
        am2 = a(2,m)
        am3 = a(3,m)
         DO i=1,nsn
           n = nod(i)
           ar(1,n)= (vg(1)-vr(1,n)) * usdt
           ar(2,n)= zero
           ar(3,n)= zero
C  
           vx1=zero
           vx2=-vg(1)*(x(3,n)-x(3,m))
           vx3=vg(1)*(x(2,n)-x(2,m))
           vxx1=zero
           vxx2=-vg(1)*vg(1)*(x(2,n)-x(2,m))
           vxx3=-vg(1)*vg(1)*(x(3,n)-x(3,m))
           a(1,n)= zero
           a(2,n)= am2 +
     .       (v(2,m)+vx2+half*dt2*vxx2-v(2,n))*usdt
           a(3,n)= am3 +
     .       (v(3,m)+vx3+half*dt2*vxx3-v(3,n))*usdt
         ENDDO
      ENDIF
C-----------------------------------------------
      IF (ireac == 0 .AND. ifail /= 1) RETURN
C-----------------------------------------------
C     Reaction forces and/or failure
C       Note : At this stage, Freac contains Fint+Fext+Fcont (Freac has not been finalized yet)
C-----------------------------------------------
      ALLOCATE(r(3,nsn))
      ALLOCATE(rr(3,nsn))
C
      DO i=1,nsn
        n = nod(i)
C
C       Reaction force is also valid with zero mass or inertia
        r(1,i)  =  ms(n)*a(1,n)  - freac(1,n)
        r(2,i)  =  ms(n)*a(2,n)  - freac(2,n)
        r(3,i)  =  ms(n)*a(3,n)  - freac(3,n)
C
        rr(1,i) =  in(n)*ar(1,n) - freac(4,n)
        rr(2,i) =  in(n)*ar(2,n) - freac(5,n)
        rr(3,i) =  in(n)*ar(3,n) - freac(6,n)
C
      ENDDO
C
      IF(ireac/=0)THEN
        DO i=1,nsn
          n = nod(i)
          l=nodreac(n)
          IF(l/=0)THEN
            fthreac(1,l)=fthreac(1,l)+r(1,i)*dt12*weight(n) ! no other kin condition can apply to a secondary node 
            fthreac(2,l)=fthreac(2,l)+r(2,i)*dt12*weight(n)
            fthreac(3,l)=fthreac(3,l)+r(3,i)*dt12*weight(n)
            fthreac(4,l)=fthreac(4,l)+rr(1,i)*dt12*weight(n) ! no other kin condition can apply to a secondary node 
            fthreac(5,l)=fthreac(5,l)+rr(2,i)*dt12*weight(n)
            fthreac(6,l)=fthreac(6,l)+rr(3,i)*dt12*weight(n)
          END IF
        END DO
      END IF
C
      IF(ifail==1)THEN
        DO i=1,nsn
          n = nod(i)
          ux    = x(1,n)-x(1,m)
          uy    = x(2,n)-x(2,m)
          uz    = x(3,n)-x(3,m)
          nn    = one/max(em20,sqrt(ux*ux+uy*uy+uz*uz))
          ux    = ux*nn
          uy    = uy*nn
          uz    = uz*nn
C
          fn    = r(1,i)*ux+r(2,i)*uy+r(3,i)*uz
          r(1,i) = r(1,i)-fn*ux
          r(2,i) = r(2,i)-fn*uy
          r(3,i) = r(3,i)-fn*uz
C     
          fn    = abs(min(fn,zero)) ! Tension <=> Reaction force FN < 0
          ft    = sqrt(r(1,i)*r(1,i)+r(2,i)*r(2,i)+r(3,i)*r(3,i))
C
          crit  = max(crit,exp(expn*log(max(em20,fn/fny)))+exp(expt*log(max(em20,ft/fty)))) ! Max over secondary nodes and over time
        ENDDO
      END IF
C
      DEALLOCATE(r,rr)
C-----------------------------------------------
      RETURN
      END