xini28.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/.
!||====================================================================
!||    xini28        ../starter/source/elements/xelem/xini28.F
!||--- called by ------------------------------------------------------
!||    xinit3        ../starter/source/elements/xelem/xinit3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../starter/source/output/message/message.F
!||    get_u_func    ../starter/source/user_interface/uaccess.F
!||    get_u_geo     ../starter/source/user_interface/uaccess.F
!||    get_u_mat     ../starter/source/user_interface/uaccess.F
!||    get_u_mid     ../starter/source/user_interface/uaccess.F
!||    get_u_mnu     ../starter/source/user_interface/uaccess.F
!||    get_u_pid     ../starter/source/user_interface/uaccess.F
!||    get_u_pnu     ../starter/source/user_interface/uaccess.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE xini28(NX     ,NAX1D ,NAX2D ,NAX3D  ,XEL   ,
     2                  VEL    ,VREL  ,IOUT  ,IPROP  ,IMAT  ,
     3                  IX    ,IDS    ,MASS  ,XINER  ,STIFM ,
     4                  STIFR ,VISCM  ,VISCR ,UVAR   ,NUVAR ,
     5                  UVARN ,NUVARN ,DTE )
      USE message_mod
C-------------------------------------------------------------------------
C     This subroutine initialize a multipurpose element
C     when element uses property TYPE28==NSTRAND.
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 IMAT     |  1      | I | R | MATERIAL NUMBER
C----------+---------+---+---+--------------------------------------------
C NX       |    1    | I | R | NUMBER OF NODES
C----------+---------+---+---+--------------------------------------------
C XEL      | 3*NX    | F | R | NODES COORDINATES
C VEL      | 3*NX    | F | R | NODES VELOCITIES
C VREL     | 3*NX    | F | R | NODES ROTATIONAL VELOCITIES
C----------+---------+---+---+--------------------------------------------
C NAX1D    |    1    | I | W | NUMBER OF EDGES TO BE DRAWN INTO ANIM
C NAX2D    |    1    | I | W | NUMBER OF FACETS TO BE DRAWN INTO ANIM
C NAX3D    |    1    | I | W | NUMBER OF SOLIDS TO BE DRAWN INTO ANIM
C----------+---------+---+---+--------------------------------------------
C IX       |   NX    | I | R | ELEMENT CONNECTIVITY
C                            | IX(J) (1<=J<=NX) : NODE J USER ID
C IDS      |    1    | I | R | ELEMENT USER IDENTIFIER
C----------+---------+---+---+--------------------------------------------
C MASS     |   NX    | F | W | NODAL MASS
C XINER    |   NX    | F | W | NODAL INERTIA (SPHERICAL)
C STIFM    |   NX    | F | W | NODAL STIFNESS (TIME STEP)
C STIFR    |   NX    | F | W | NODAL ROTATION STIFNESS (TIME STEP)
C VISCM    |   NX    | F | W | NODAL VISCOSITY (TIME STEP)
C VISCR    |   NX    | F | W | NODAL ROTATION VISCOSITY (TIME STEP)
C----------+---------+---+---+--------------------------------------------
C UVAR     |NUVAR    | F | W | USER ELEMENT VARIABLES
C NUVAR    |    1    | I | R | NUMBER OF USER ELEMENT VARIABLES
C UVARN    |NUVARN*NX| F | W | USER ELEMENT VARIABLES
C NUVARN   |    1    | I | R | NUMBER OF USER ELEMENT VARIABLES PER NODE.
C----------+---------+---+---+--------------------------------------------
C DTE      |    1    | F | W | MAXIMUM ELEMENT TIME STEP TO BE USED BY RADIOSS.
C----------+---------+---+---+--------------------------------------------
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,KPROP)
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) NUMBERS. 
C         SEE LECG28 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 NUMBERS(function 
C         referred by users materials).
C         SEE LECM28 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-------------------------------------------------------------------------
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-------------------------------------------------------------------------
C INTEGER MID = 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 PID = 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-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   A n a l y s e   M o d u l e
C-----------------------------------------------
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,NUVAR,NUVARN,IPROP,IMAT,
     .        NX ,NAX1D ,NAX2D ,NAX3D , IX(NX), IDS,
     .        GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
     .        KFUNC,KMAT,KPROP
      my_real
     .        xel(3,nx),vel(3,nx),vrel(3,nx),
     .        mass(nx) ,xiner(nx) ,stifm(nx) ,
     .        stifr(nx),viscm(nx) ,viscr(nx) ,uvar(nuvar) ,
     .        uvarn(nuvarn*nx), dte,
     .        get_u_mat,get_u_geo,get_u_func,ffac
      EXTERNAL get_u_pnu,get_u_mnu,get_u_mat,get_u_geo,get_u_pid,
     .         get_u_mid,get_u_func
      parameter(kfunc=29)
      parameter(kmat=31)
      parameter(kprop=33)
C=======================================================================
C
C     EXAMPLE : NSTRAND.
C
C=======================================================================
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real 
     .       ms,xk,xc,epstot,f,dfdx,rho,stif,deps,g,dgdx,l0,
     .       lprev, lnext,
     .       xm, xkm, xcm, fact, xn, dtc, dtk
      INTEGER I,K,NB1,NB2,NB3,MB1,MB2,MB3,MB4,MB5,IFUNCT,IFV
C-----------------------------------------------
C       POLYLINE TO DRAW INTO ANIM.
        NAX1D=nx-1
        nax2d=0
        nax3d=0
C-----------------------------------------------
C       initial total length UVAR(NB1:NB1)
        nb1=1
C       previous elongation  UVAR(NB2:NB2)
        nb2=nb1+1
C       previous force       UVAR(NB3:NB3)
        nb3=nb2+1
C       initial nodes masses UVARN(MB1:MB1+NX-1)
        mb1=1
C       forces into strands UVARN(MB2:MB2+NX-1)
C       using UVARN(MB2:MB2+NX-2) only.
        mb2=mb1+nx
C       strands initial length UVARN(MB3:MB3+NX-1)
C      using UVARN(MB3:MB3+NX-2) only.
        mb3=mb2+nx
C       strands elongations UVARN(MB4:MB4+NX-1)
C       using UVARN(MB4:MB4+NX-2) only.
        mb4=mb3+nx
C       strands internal energy UVARN(MB5:MB5+NX-1)
C       using UVARN(MB5:MB5+NX-2) only.
        mb5=mb4+nx
C-----------------------------------------------
C       ELEMENT CHECK
        IF(nx<3)THEN
C         WRITE(IOUT,*)
C     . ' ** ERROR NSTRAND : LESS THAN 3 NODES, ELEMENT=',           
C     .    IDS
           CALL ancmsg(msgid=381,
     .                 msgtype=msgerror,
     .                 anmode=aninfo,
     .                 i1=ids)
C         IERR=IERR+1
C         GOTO 999
        ENDIF
C-------
        rho   =get_u_geo(3,iprop)
C-------
C       POLYLINE LENGTH & NODES MASS.
C-------
        lprev=
     .       sqrt((xel(1,2)-xel(1,1))*(xel(1,2)-xel(1,1))
     .           +(xel(2,2)-xel(2,1))*(xel(2,2)-xel(2,1))
     .           +(xel(3,2)-xel(3,1))*(xel(3,2)-xel(3,1)))
        uvar(nb1) =lprev
        uvarn(mb3)=lprev
C
        mass(1)  =half*rho*lprev
        IF (lprev<=em15) THEN
C         WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL STRAND LENGTH, ELEMENT=',
C     .      IDS
C        IERR=IERR+1
C        GOTO 999
            CALL ancmsg(msgid=382,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ids)
        ENDIF
        IF(mass(1)<=em15)THEN
C         WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL NODAL MASS, ELEMENT=',
C     .      IDS
C        IERR=IERR+1
C        GOTO 999
            CALL ancmsg(msgid=383,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ids)
        ENDIF
        DO k=2,nx-1
          lnext=
     .       sqrt((xel(1,k+1)-xel(1,k))*(xel(1,k+1)-xel(1,k))
     .           +(xel(2,k+1)-xel(2,k))*(xel(2,k+1)-xel(2,k))
     .           +(xel(3,k+1)-xel(3,k))*(xel(3,k+1)-xel(3,k)))
          IF (lnext<=em15) THEN
C           WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL STRAND LENGTH, ELEMENT=',
C     .      IDS
C          IERR=IERR+1
C           GOTO 999
            CALL ancmsg(msgid=382,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ids)
          ENDIF
          mass(k)       = half*rho*(lprev+lnext)
          uvarn(mb3+k-1)=lnext
          uvar(nb1)     =uvar(nb1)+lnext
          IF(mass(k)<=em15)THEN
C           WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL NODAL MASS, ELEMENT=',
C     .      IDS
C          IERR=IERR+1
C           GOTO 999
            CALL ancmsg(msgid=383,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ids)
          ENDIF
          lprev=lnext
        ENDDO
        mass(nx)  = half*rho*lprev
        IF(mass(nx)<=em15)THEN
C         WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL NODAL MASS, ELEMENT=',
C     .      IDS
C        IERR=IERR+1
C         GOTO 999
            CALL ancmsg(msgid=383,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ids)
        ENDIF
C------------------------------------------
        xk  =get_u_geo(4,iprop)
        dfdx=zero
        ifunct=get_u_pnu(1,iprop,kfunc)
C
        xc  =get_u_geo(5,iprop)
        ifv=get_u_pnu(2,iprop,kfunc)
        ffac=get_u_geo(12,iprop)
C
        IF (ifunct==0.AND.ifv==0) THEN
C------------------------------------------
C        LINEAR ELASTIC
         l0  =uvar(nb1)
         stif=xk
         f   =zero
        ELSEIF (ifunct==0.AND.ifv/=0) THEN
C------------------------------------------
C        G(Deps) only.
         l0  =uvar(nb1)
         stif=zero
         f   =one
        ELSE
C------------------------------------------
C        NON LINEAR ELASTIC
         l0    =uvar(nb1)
         epstot=zero
         f     =get_u_func(ifunct,epstot,dfdx)
         stif  =ffac*dfdx
        ENDIF
C------------------------------------------
        deps=zero
        dgdx=zero
        g   =one
        IF (ifv/=0) THEN
         g=get_u_func(ifv,deps,dgdx)
         stif=stif*g
        ENDIF
C-------
C        CHECK NSTRAND STIFNESS AND DAMPING.
         IF(      stif/uvar(nb1)<=em15
     .      .AND.(f*dgdx+xc)/uvar(nb1)<=em15)THEN
C          WRITE(IOUT,*)
C     .' ** ERROR NSTRAND : NULL NODAL STIFFNESS & DAMPING, ELEMENT=',
C     .     IDS
C         IERR=IERR+1
C          GOTO 999
           CALL ancmsg(msgid=384,
     .                 msgtype=msgerror,
     .                 anmode=aninfo,
     .                 i1=ids)
         ENDIF
C-------
C       RETURN ELEMENT TIME STEP.
        xn=nx
        dte  = ep20
        DO k=1,nx-1
          xm   = rho*uvarn(mb3+k-1)
C         rigidite tangente
C tous les brins ont la meme longueur correspond a :
C         XKM  = STIF*(XN-1.)/L0
          xkm  = stif/uvarn(mb3+k-1)
C tous les brins ont la meme longueur correspond a :
C         XCM  = (F*DGDX+XC)*(XN-1.)/L0
          xcm  = (f*dgdx+xc)/uvarn(mb3+k-1)
          IF(xcm+xkm<em15)xm  =one
          xkm= max(em15,xkm)
C           DTK=(SQRT(0.25*XCM*XCM+XM*XKM)-0.5*XCM)/MAX(EM15,XKM)
C           DTC=0.5 * 2.*XM/MAX(EM15,XCM)
          dtk=(sqrt(xcm*xcm+xm*xkm)-xcm)/max(em15,xkm)
          dtc=xm/max(em15,xcm)
          IF (dtk==zero) THEN
C          XKM==0.
           dtk=dtc
          ELSE
           dtk=min(dtk,dtc)
          ENDIF
          dte=min(dte,dtk)
        ENDDO
C-------
C        FOR NODAL TIME STEP COMPUTATION.
         viscm(1) =(f*dgdx+xc)/uvarn(mb3)
         stifm(1) =stif/uvarn(mb3)
         DO k=2,nx-1
          fact     =one/uvarn(mb3+k-2) + one/uvarn(mb3+k-1)
          stifm(k) =stif*fact
          viscm(k) =(f*dgdx+xc)*fact
         ENDDO
         viscm(nx) =(f*dgdx+xc)/uvarn(mb3+nx-2)
         stifm(nx) =stif/uvarn(mb3+nx-2)
C-------
        DO k=1,nx
         xiner(k) = zero
         stifr(k) = zero
         viscr(k) = zero
 
        ENDDO
C-------
         DO k=1,nx
          uvarn(mb1+k-1)=mass(k)
         ENDDO
C-------
      WRITE(iout,1000) ids,l0,rho*l0,stif/l0
 1000 FORMAT(' NSTRAND ELEMENT CHECKING :',/,
     .       ' ------------------------  ',/,
     .       ' ELEMENT IDENTIFIER . . . .',i8/,
     .       ' TOTAL LENGTH . . . . . . .',e12.4/,
     .       ' MASS . . . . . . . . . . .',e12.4/,
     .       ' INITIAL GLOBAL STIFFNESS .',e12.4//)
C-------
      RETURN
      END