rmatacce.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/.
!||====================================================================
!||    rmatacce   ../engine/source/materials/mat/mat013/rmatacce.F
!||--- called by ------------------------------------------------------
!||    resol      ../engine/source/engine/resol.F
!||====================================================================
      SUBROUTINE rmatacce(RBYM ,ARBYM ,ARRBYM ,VRBYM,VRRBYM,
     2                    IRBYM,LNRBYM,X      ,A    ,AR    ,
     3                    V    ,VR    ,KIND   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IRBYM(NIRBYM,*),LNRBYM(*), KIND(*)
C     REAL
      my_real
     .   x(3,*), a(3,*), ar(3,*), vr(3,*),v(3,*),
     .   vrrbym(3,*), arbym(3,*),arrbym(3,*),vrbym(3,*),
     .   rbym(nfrbym,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J, N, M, NSL, NSLT
C     REAL
      my_real rx, ry, rz, msx, msy, msz, fsx, fsy, fsz,
     .        xg,yg,zg,usdt,v1x2,v2x1,v2x3,v3x2,v3x1,v1x3,
     .        vx2,vx1,vx3,vg(3)
 
C======================================================================|
       usdt = one /dt12
 
!$OMP DO
 
      DO m=1,nrbym
C     
        nsl = irbym(2,m)
        xg = rbym(2,m)
        yg = rbym(3,m)
        zg = rbym(4,m)
C       
        nslt= kind(m)-1
C
        vg(1) = vrrbym(1,m) + arrbym(1,m)*dt12
        vg(2) = vrrbym(2,m) + arrbym(2,m)*dt12
        vg(3) = vrrbym(3,m) + arrbym(3,m)*dt12                
        DO i=1,nsl
          n = lnrbym(nslt + i)
C          
          rx = x(1,n) - xg
          ry = x(2,n) - yg
          rz = x(3,n) - zg
C
          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)*ry
          v2x1 = vg(2)*rx
          v2x3 = vg(2)*rz
          v3x2 = vg(3)*ry
          v3x1 = vg(3)*rx
          v1x3 = vg(1)*rz
C
          vx1 = v2x3 - v3x2
          vx2 = v3x1 - v1x3
          vx3 = v1x2 - v2x1
C          
          a(1,n)= arbym(1,m) + usdt*(
     .     vrbym(1,m) + vx1+half*dt2*(vg(2)*vx3-vg(3)*vx2)-v(1,n) )
          a(2,n)= arbym(2,m) + usdt*(
     .      vrbym(2,m) + vx2+half*dt2*(vg(3)*vx1-vg(1)*vx3)-v(2,n))
          a(3,n)= arbym(3,m) + usdt*(
     .      vrbym(3,m)+vx3+half*dt2*(vg(1)*vx2-vg(2)*vx1)-v(3,n))
C compute rotational velocity
C
        ENDDO
        
          vrbym(1,m) = vrbym(1,m) + arbym(1,m)*dt12
          vrbym(2,m) = vrbym(2,m) + arbym(2,m)*dt12
          vrbym(3,m) = vrbym(3,m) + arbym(3,m)*dt12
C       
        vrrbym(1,m) = vg(1)
        vrrbym(2,m) = vg(2)
        vrrbym(3,m) = vg(3)      
      ENDDO
 
!$OMP END DO
 
C---
      RETURN
      END