ruser46.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/.
!||====================================================================
!||    ruser46      ../engine/source/elements/spring/ruser46.F
!||--- called by ------------------------------------------------------
!||    rforc3       ../engine/source/elements/spring/rforc3.F
!||--- calls      -----------------------------------------------------
!||    get_u_func   ../engine/source/user_interface/ufunc.F
!||    get_u_mid    ../engine/source/user_interface/upidmid.F
!||    get_u_mnu    ../engine/source/user_interface/upidmid.F
!||    get_u_pid    ../engine/source/user_interface/upidmid.F
!||    get_u_pnu    ../engine/source/user_interface/upidmid.F
!||    get_u_time   ../engine/source/user_interface/uaccess.F
!||====================================================================
      SUBROUTINE ruser46(NEL,IOUT   ,IPROP  ,UVAR   ,NUVAR  ,
     2             FX      ,FY      ,FZ     ,XMOM   ,YMOM   ,
     3             ZMOM    ,E       ,OFF    ,STIFM  ,STIFR  ,
     4             VISCM   ,VISCR   ,MASS   ,XINER  ,DT     ,
     5             XL      ,VX      ,RY1    ,RZ1    ,RX     ,
     6             RY2     ,RZ2     ,FR_WAVE)
C-------------------------------------------------------------------------
C     This subroutine compute springs forces and moments.
C-------------------------------------------------------------------------
C----------+---------+---+---+--------------------------------------------
C VAR      | SIZE    |TYP| RW| DEFINITION
C----------+---------+---+---+--------------------------------------------
C IOUT     |  1      | I | R | OUTPUT FILE UNIT (L00 file)
C IPROP    |  1      | I | R | PROPERTY NUMBER
C----------+---------+---+---+--------------------------------------------
C XL       |   NEL   | F | R | ELEMENT LENGTH
C----------+---------+---+---+--------------------------------------------
C UVAR     |NUVAR*NEL| F |R/W| USER ELEMENT VARIABLES
C NUVAR    |  1      | I | R | NUMBER OF USER ELEMENT VARIABLES
C----------+---------+---+---+--------------------------------------------
C The user routine does not need to return elements mass in vector MASS : 
C this vector is not used by RADIOSS since version 4.1E.
C The mass which is used by RADIOSS for time step (and output) is the 
C initial mass which was returned by user routine RINI29 into starter.
C-------------------------------------------------------------------------
C FUNCTION 
C-------------------------------------------------------------------------
C INTEGER II = GET_U_PNU(I,IP,KK)
C         IFUNCI = GET_U_PNU(I,IP,KFUNC)
C         IPROPI = GET_U_PNU(I,IP,KFUNC)
C         IMATI  = GET_U_PNU(I,IP,KMAT)
C         I     :     VARIABLE INDEX(1 for first variable,...)
C         IP    :     PROPERTY NUMBER
C         KK    :     PARAMETER KFUNC,KMAT,KPROP
C         THIS FUNCTION RETURN THE USER STORED FUNCTION(IF KK=KFUNC), 
C         MATERIAL(IF KK=KMAT) OR PROPERTY(IF KK=KPROP) NUMBER. 
C         SEE LECG29 FOR CORRESPONDING ID STORAGE.
C-------------------------------------------------------------------------
C INTEGER IFUNCI = GET_U_MNU(I,IM,KFUNC)
C         I     :     VARIABLE INDEX(1 for first function)
C         IM    :     MATERIAL NUMBER
C         KFUNC :     ONLY FUNCTION ARE YET AVAILABLE.
C         THIS FUNCTION RETURN THE USER STORED FUNCTION NUMBER(function 
C         referred by users materials).
C         SEE LECM29 FOR CORRESPONDING ID STORAGE.
C-------------------------------------------------------------------------
C my_real PARAMI = GET_U_GEO(I,IP)
C         I     :     PARAMETER INDEX(1 for first parameter,...)
C         IP    :     PROPERTY NUMBER
C         THIS FUNCTION RETURN THE USER GEOMETRY PARAMETERS 
C         NOTE: IF(IP==IPROP) UPARAG(I) == GET_U_GEO(I,IPROP)
C         see lecg30 for storage
C-------------------------------------------------------------------------
C my_real PARAMI = GET_U_MAT(I,IM)
C         I     :     PARAMETER INDEX(1 for first parameter,...)
C         IM    :     MATERIAL NUMBER
C         THIS FUNCTION RETURN THE USER MATERIAL PARAMETERS 
C         NOTE: GET_U_MAT(0,IMAT) RETURN THE DENSITY
C         see lecm29,30,31 for storage
C-------------------------------------------------------------------------
C INTEGER PID = GET_U_PID(IP)
C         IP    :     PROPERTY NUMBER
C         THIS FUNCTION RETURN THE USER PROPERTY ID CORRESPONDING TO
C         USER PROPERTY NUMBER IP. 
C-------------------------------------------------------------------------
C INTEGER MID = GET_U_MID(IM)
C         IM   :     MATERIAL NUMBER
C         THIS FUNCTION RETURN THE USER MATERIAL ID CORRESPONDING TO
C         USER MATERIAL NUMBER IM. 
C-------------------------------------------------------------------------
C my_real Y = GET_U_FUNC(IFUNC,X,DXDY)
C         IFUNC :     function number obtained by 
C                     IFUNC = GET_U_MNU(I,IM,KFUNC) or IFUNC = GET_U_PNU(I,IP,KFUNC)
C         X     :     X value
C         DXDY  :     slope dX/dY
C         THIS FUNCTION RETURN Y(X)
C-------------------------------------------------------------------------
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C----------------------------------------------------------
C   D u m m y   A r g u m e n t s   a n d   F u n c t i o n
C----------------------------------------------------------
      INTEGER IOUT,NEL,NUVAR,IPROP,
     .        GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
     .        KFUNC,KMAT,KPROP
      my_real
     .   UVAR(NUVAR,*),DT ,
     .   fx(*), fy(*), fz(*), e(*), vx(*),mass(*) ,xiner(*),
     .   ry1(*), rz1(*), off(*), xmom(*), ymom(*),
     .   zmom(*), rx(*), ry2(*), rz2(*),xl(*),
     .   stifm(*) ,stifr(*) , viscm(*) ,viscr(*) ,fr_wave(*) ,
     .   get_u_mat, get_u_geo, get_u_func, get_u_time
      EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
     .         get_u_mat,get_u_geo, get_u_func, get_u_time
      parameter(kfunc=29)
!      PARAMETER (KMAT=31)
!      PARAMETER (KPROP=33)
C=======================================================================
C
C     EXAMPLE 2 : Elastoplastic truss defined with 1 user property.
C                 
C               
C=======================================================================
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,IFUNC1,IFUNC2,IFUNC3,IFUNC4,IFUNC5,EPSI
      my_real
     .        ELASTIF,X,DXDY,XLIM1,XLIM2,FMX,FMN,DX,DAMP,
     .        FUX1,FUX2,FUX3,FUX4,DEBUG,SLOPE
      my_real
     .        tt,newlen,deltalen,deltalendot,scalet,scalex,scalev,scalef
C-----------------------------------------------
C
      elastif= get_u_geo(2,iprop)
      xlim1  = get_u_geo(3,iprop) !vmax
      xlim2  = get_u_geo(4,iprop) !fmax
      damp   = get_u_geo(6,iprop)
      scalet = get_u_geo(8,iprop)
      scalex = get_u_geo(9,iprop)
      scalev = get_u_geo(10,iprop)
      scalef = get_u_geo(11,iprop)
 
      epsi   = int(get_u_geo(7,iprop))
      ifunc1= get_u_pnu(1,iprop,kfunc) 
      ifunc2= get_u_pnu(2,iprop,kfunc) 
      ifunc3= get_u_pnu(3,iprop,kfunc) 
      ifunc4= get_u_pnu(4,iprop,kfunc) 
      ifunc5= get_u_pnu(5,iprop,kfunc)
C 
      tt= get_u_time()
      tt = tt / scalet
C 
      deltalendot = zero
C
C      DEBUG = GET_U_FUNC(IFUNC5,TT,DXDY)
C Loop over elements
      DO i=1,nel
          dx = dt * vx(i)
          x = uvar(1,i) + dx
          uvar(1,i) = x 
          newlen = uvar(3,i) + uvar(1,i)!UVAR(3,I)=initial length of the spring at time=0
C
          IF (epsi == 0) THEN
             deltalen   = newlen/uvar(3,i) - one
                IF (dt == zero) THEN
                   deltalen = zero
                ELSE
                   deltalendot = newlen/uvar(3,i)*vx(i)/xlim1
                ENDIF
          ELSE
             deltalen = x
             deltalendot = vx(i)
          ENDIF
C
C Active force vs. time
          fux1 = get_u_func(ifunc1,tt,dxdy)
C Active force vs. elongation
          IF(epsi/= 0)deltalen = deltalen / scalex 
          fux2 = get_u_func(ifunc2,deltalen,dxdy)
C Active force vs. velocity
C          FUX3 = GET_U_FUNC(IFUNC3,VX(I),DXDY)
          deltalendot = deltalendot / scalev 
          fux3 = get_u_func(ifunc3,deltalendot,dxdy)
C Passive force vs. elongation
          fux4 = get_u_func(ifunc4,deltalen,dxdy)
          fux4 = fux4 * scalef
C Check if velocity > limit1 ...if YES, then use XLIM1 = max_vel
C Sum FX + damping*velocity
          IF (vx(i)>xlim1)THEN
             fx(i) = xlim2 * fux1 * fux2 * fux3 + fux4 + damp*xlim1
          ELSEIF (vx(i)<-xlim1)THEN
             fx(i) = xlim2 * fux1 * fux2 * fux3 + fux4 - damp*xlim1
          ELSE
             fx(i) = xlim2 * fux1 * fux2 * fux3 + fux4 + damp*vx(i)
          ENDIF
C
C        ***       TIMESTEP       ***
C=======================================================================
C       TIME STEP
C=======================================================================
C          STIFM(I) = ELASTIF / XL(I)
          IF(deltalen /= zero)THEN
            slope = abs(fx(i)-uvar(4,i))/abs(dx)
C            WRITE(IOUT,*) 'Slope: FR_WAVE = ZERO'
          ELSE
            slope = elastif
          ENDIF
C
          
C
          stifm(i) = slope
          stifr(i) = zero
          viscm(i) = zero
          viscr(i) = zero
          xiner(i) = zero
          uvar(4,i)= fx(i)
      ENDDO
C-------------------------------
      RETURN
      END