freimpl.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/.
!||====================================================================
!||    freimpl                ../engine/source/input/freimpl.F
!||--- called by ------------------------------------------------------
!||    freform                ../engine/source/input/freform.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../engine/source/output/message/message.F
!||    arret                  ../engine/source/system/arret.F
!||    order_dtf              ../engine/source/input/freimpl.F
!||    wriusc2                ../engine/source/input/wriusc2.F
!||--- uses       -----------------------------------------------------
!||    imp_dyna               ../engine/share/modules/impbufdef_mod.F
!||    imp_kbcs               ../engine/share/modules/impbufdef_mod.F
!||    imp_pcg_proj           ../engine/share/modules/impbufdef_mod.F
!||    imp_spbrm              ../engine/share/modules/impbufdef_mod.F
!||    message_mod            ../engine/share/message_module/message_mod.F
!||    names_and_titles_mod   ../common_source/modules/names_and_titles_mod.F
!||====================================================================
      SUBROUTINE freimpl(IKAD,KEY0,KIMPL) 
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE imp_dyna
      USE imp_kbcs
      USE imp_pcg_proj
      USE imp_spbrm
      USE message_mod
      USE names_and_titles_mod, ONLY: ncharline100
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
C-----------------------------------------------
      INTEGER IKAD(0:*),KIMPL
      CHARACTER KEY0(*)*5
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "impl1_c.inc"
#include      "impl2_c.inc"
#include      "parit_c.inc"
#include      "com01_c.inc"
#include      "com06_c.inc"
#include      "buckcom.inc"
#include      "scr06_c.inc"
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      INTEGER NVAR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, NBC, K, IKEY,IM,J,NJ,KK
      CHARACTER  TITLE*72, KEY2*5, KEY3*5, KEY4*5
      CHARACTER(LEN=NCHARLINE100)::CARTE
C----------------------------------------
      ikey=kimpl
      impl_s=0
      idyna=0
      iline=0
      isprb=0
      isolv=0
      insolv=0
      idtc=0
      im=0
      ikg=1
      kz_tol=zero
      sk_int=zero
      d_tol=zero
      lprint=0
      nprint=0
      impdeb=0
      solvnfo=0
      prstifmat=0
      prstifmat_tol=zero
      prstifmat_nc=1
      prstifmat_it=0 
      impmv=1
      isigini=0
      ilintf=0
      iprec = 0
      l_lim = 0
      itol  = 0
      l_tol =zero
      dt_imp = zero
      dt_min = zero
      dt_max = zero
      imp_rby=0
      imp_int=0
      isprn = 1
C      INTP_C = 0
C -----after debugging on int24 spmd, change defaut to INTP_C=1 (INTP_C=0 suppressed good for maintenance)
      intp_c = 1
      l_bfgs = 0
C      IRREF = 0
      irref = 1
      iqstat = 0
      ibuckl = 0
      iscau = 0
      imp_lr=0
      ikproj=0
      ismdisp = 0
      IF(ikad(ikey)/=ikad(ikey+1))THEN
       k=0
       impl_s=1
       ncinp=1
       n_pat = 1
       imp_chk = 0
       imp_int7 = 0
       ittoff = 0
       scal_dtq = one
       idy_damp=0
       iautspc = 1
       itrmax = 0
       msg_lvl = 0
       b_order =0
       b_mcore =0
       irefi = 0
       iline_s = 0
       nls_lim = 0
       ls_tol = zero
       ndiver = 0
       ikt = 0
       ndtfix = 0
       ikpres = 1
       n_tolu=zero
       n_tolf=zero
       n_tole=zero
       ncy_max = 0
       rf_min = zero
       rf_max = zero
       ipupd = 0
       tol_div = zero
       m_vs = 0
       ipro_s0=0
       iikgoff = 1
       m_msg = 0
       m_order =0
       m_ocore =0
       irig_m = 0
 1160  READ(iusc1,rec=ikad(ikey)+k,fmt='(7X,A,1X,A,1X,A,25X,I10)',err=9990)key2,key3,key4,nbc
       k=k+1
C----------------------------
C      Dynamic implicit
C----------------------------
       IF(key2(1:4)=='DYNA')THEN
        IF (idyna==0) idyna=1
        IF(key3(1:4)=='DAMP')THEN
         idy_damp=1
         CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
         READ(iusc2,*) dampa_imp,dampb_imp
        ELSE IF(key3(1:3)=='FSI')THEN
         WRITE(6,*) "ERROR: /IMPL/DYNA/FSI IS A DEPRECATED FEATURE"
          GOTO 9990
        ELSE
         READ(key3,'(I2)')im
         idyna=max(idyna,im)
         IF(idyna==1)THEN
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)hht_a
         ELSEIF(idyna==2)THEN
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)newm_a,newm_b
         ELSE 
          hht_a =-em20
         ENDIF
        ENDIF
C----------------------------
C      Implicit linear
C----------------------------
       ELSEIF(key2(1:4)=='LINE')THEN
        iline=1
        IF(key3(1:5)=='INTER') THEN
         READ(key4,'(I5)')ilintf
         ilintf = max(2,ilintf)
        ELSEIF(key3(1:5)=='SCAUC') THEN
         iscau = 1
        ENDIF
       ELSEIF(key2(1:5)=='MONVO')THEN
        IF(key3(1:3)=='OFF')impmv=0
       ELSEIF(key2(1:5)=='SPRIN')THEN
         IF(key3(1:4)=='NONL')THEN
          isprn = 1 
         ELSEIF(key3(1:4)=='LINE')THEN
          isprn = 0 
         ELSE
          GOTO 9990
         ENDIF
       ELSEIF(key2(1:5)=='PREPA')THEN
        READ(key3,'(I2)')n_pat
       ELSEIF(key2(1:5)=='PROJV')THEN
        READ(key3,'(I2)') m_vs
       ELSEIF(key2(1:5)=='PROSI')THEN
        READ(key3,'(I2)') ipro_s0
C----------------------------
C      Implicit check
C----------------------------
       ELSEIF(key2(1:5)=='CHECK')THEN
        imp_chk = 1
C----------------------------
C      Implicit quasi-static
C----------------------------
       ELSEIF(key2(1:5)=='QSTAT')THEN
        iqstat = 1
        IF(key3(1:5)=='DTSCA')THEN
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)scal_dtq
        ELSEIF(key3(1:5)=='MRIGM')THEN
          irig_m = 1
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*,err=520,END=520)E_REF(1),E_REF(2),E_REF(3)
          IF (e_ref(1)>0.AND.e_ref(2)>0.AND.e_ref(3)>0) irig_m = 2
 520    CONTINUE
        ELSE
          READ(key3,'(I2)')im
          iqstat=max(iqstat,im)
        ENDIF
C----------------------------
C      spring-back
C----------------------------
       ELSEIF(key2(1:4)=='SPRB')THEN
        isprb=1
C----------------------------
C      print-out
C----------------------------
       ELSEIF(key2=='PRINT')THEN
         IF(key3(1:4)=='LINE')THEN
          READ(key4,'(I5)')lprint
         ELSEIF(key3(1:4)=='NONL')THEN
          READ(key4,'(I5)')nprint
         ELSEIF(key3(1:4)=='STIF')THEN
          prstifmat = 1    
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)prstifmat_tol,prstifmat_nc,prstifmat_it          
         ELSE
          GOTO 9990
         ENDIF
C----------------------------
C      Linear SOLVER
C----------------------------
       ELSEIF(key2(1:4)=='SOLV')THEN
        READ(key3,'(I2)')isolv
        CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
        READ(iusc2,*)iprec,l_lim,itol,l_tol
        IF (isolv==3) imumpsd=l_lim
C----------------------------
C      BCS parameters
C----------------------------
       ELSEIF(key2(1:4)=='SBCS')THEN
        IF(key3(1:5)=='MSGLV')THEN
         READ(key4,'(I2)')msg_lvl
        ELSEIF(key3(1:5)=='ORDER')THEN
         READ(key4,'(I2)')b_order
C-------0 default 1 :MMD  2 :metis       
        ELSEIF(key3(1:5)=='outco')THEN
             B_MCORE=1
        ELSE
          GOTO 9990
        ENDIF
C----------------------------
C      MUMPS parameters
C----------------------------
       ELSEIF(KEY2(1:5)=='mumps')THEN
        IF(KEY3(1:5)=='msglv')THEN
         READ(key4,'(I2)')m_msg
        ELSEIF(key3(1:5)=='ORDER')THEN
         IF(key4(1:5)=='METIS')THEN
          m_order = 5
         ELSEIF(key4(1:4)=='PORD')THEN
          m_order = 4
         END IF
C-------0 default 1 :MMD  2 :metis       
        ELSEIF(key3(1:5)=='OUTCO')THEN
             m_ocore=1
        ELSEIF(key3(1:5)=='AUTOC')THEN
             m_ocore=-1
        ELSE
          GOTO 9990
        ENDIF
C----------------------------
C      Nonlinear SOLVER
C----------------------------
       ELSEIF(key2(1:4)=='NONL')THEN
        IF(key3(1:5)=='KTANG')THEN
         ikt = 1
        ELSEIF(key3(1:5)=='KTFUL')THEN
         ikt = 2
        ELSEIF(key3(1:5)=='KTFU8')THEN
         ikt = 3
        ELSEIF(key3(1:5)=='KTCON')THEN
         ikt = 4
        ELSEIF(key3(1:5)=='PITER')THEN
         READ(key4,'(I5)') ipupd
        ELSEIF(key3(1:5)=='SMDIS')THEN
         ismdisp = 1
        ELSEIF(key3(1:5)=='SOLVI')THEN
         solvnfo = 1
        ELSE
         READ(key3,'(I2)')insolv
         CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
         READ(iusc2,'(a)')title
         READ(title,*)n_lim,nitol,n_tol
         IF (nitol>10) THEN
          SELECT CASE (nitol)
           CASE(12)                                            
            READ(title,*)n_lim,nitol,n_tole,n_tolf
           CASE(13)                                   
            READ(title,*)n_lim,nitol,n_tole,n_tolu
           CASE(23)                                   
            READ(title,*)n_lim,nitol,n_tolf,n_tolu
           CASE(123)                                            
            READ(title,*)n_lim,nitol,n_tole,n_tolf,n_tolu
          END SELECT      
         ENDIF !(NITOL>10)
         IF(nitol==1.AND.irref==1) irref = 0
        ENDIF
       ELSEIF(key2(1:5)=='SINIT')THEN
        isigini=1
       ELSEIF(key2(1:5)=='LBFGS')THEN
          READ(key3,'(I5)') l_bfgs
C----------------------------
C      Step Control
C----------------------------
       ELSEIF(key2=='DTINI')THEN
         CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
         READ(iusc2,*)dt_imp
       ELSEIF(key2(1:2)=='DT')THEN
         IF(key3(1:4)=='STOP')THEN
           CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
           READ(iusc2,*)dt_min,dt_max
C----------------------------
C-------------fix point for time step----
C----------------------------
         ELSEIF(key3(1:4)=='FIXP')THEN
          kk =k
          DO i=1,nbc
           READ(iusc1,rec=ikad(ikey)+kk,fmt='(A)',err=9990)carte
           CALL wriusc2(ikad(ikey)+kk,1,key0(ikey))
           nj = nvar(carte)
           IF ((ndtfix+nj)>100) THEN
            nj = 100-ndtfix
            WRITE(istdo,*) 
     .       ' ** WARNING ** : MAXIMUM 100 FIX POINTS PERMITTED '
           ENDIF
           READ(iusc2,*,err=9990,END=9990)(DTIMPF(NDTFIX+J),J=1,NJ)
           kk=kk+1
           ndtfix = ndtfix + nj
          ENDDO
          CALL order_dtf(ndtfix,dtimpf)
         ELSE
           READ(key3,'(I2)')im
           IF (idtc>0.AND.im>0) GOTO 9990
           idtc=im
           IF(im==1)THEN
            CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
            READ(iusc2,*)nl_dtp,scal_dtp,nl_dtn,scal_dtn
           ELSEIF(im==2)THEN
            CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
            READ(iusc2,*)nl_dtp,alen0,nl_dtn,scal_dtn,scal_dtp
           ELSEIF(im==3)THEN
            CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
            READ(iusc2,*)nl_dtp,alen0,nl_dtn,scal_dtn,scal_dtp,ial_m,
     .                   scal_riks
           ELSE
            GOTO 9990
           ENDIF
         ENDIF
C----------------------------
C      NCYCLE stop
C----------------------------
       ELSEIF(key2=='NCYCL')THEN
         IF(key3(1:4)=='STOP')THEN
           CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
           READ(iusc2,*)ncy_max
         ELSE
           GOTO 9990
         ENDIF
C----------------------------
C      interface Control
C----------------------------
       ELSEIF(key2(1:5)=='INTER')THEN
         IF(key3(1:5)=='TTOFF')THEN
          ittoff = 1
         ELSEIF(key3(1:5)=='SINT7')THEN
          READ(key4,'(I2)')imp_int7
C-----0 nonlinear, 1: linear 2: constant---
          imp_int7= min(2,imp_int7)
C---------will be suppressed in the 14.0.210          
         ELSEIF(key3(1:5)=='KNONL')THEN
C-----0 nonlinear, 1: linear ----
          READ(key4,'(I2)')im
          intp_c = -im -1 
C         ELSEIF(KEY3(1:5)=='KCOMP'.AND.INTP_C==0)THEN
         ELSEIF(key3(1:5)=='KCOMP')THEN
c          INTP_C = 1
C-----hide option to activate KG in int24 (/IMPLICIT should be defined in Starter)        
         ELSEIF(key3(1:4)=='KGON')THEN
          iikgoff = 0
         ELSE
          GOTO 9990
         ENDIF
C----------------------------
C      R_ref options
C----------------------------
       ELSEIF(key2(1:4)=='RREF')THEN
         irref = 2
         IF(key3(1:3)=='OFF') THEN
          irref = 0
         ELSEIF(key3(1:5)=='INTER')THEN
C-----0 agressive, 1: moyen ----2--faible  3  non--4 non sauf 1er---
          READ(key4,'(I2)')im
          irefi = im
         ELSEIF(key3(1:5)=='LIMIT')THEN
           CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
           READ(iusc2,*)rf_min,rf_max
         ENDIF
C----------------------------
C      divergence criteria
C----------------------------
       ELSEIF(key2(1:5)=='DIVER')THEN
         IF(key3(1:3)=='TOL')THEN
           CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
           READ(iusc2,*)tol_div
C-----num. of diver---
         ELSE
          READ(key3,'(I2)')im
          ndiver = im
          IF (ndiver ==0) ndiver=-1
         END IF
C----------------------------
C      Geometrical stifness
C----------------------------
       ELSEIF(key2(1:5)=='GSTIF')THEN
        IF(key3(1:3)=='OFF')ikg=0
C----------------------------
C      Geometrical stifness
C----------------------------
       ELSEIF(key2(1:5)=='PSTIF')THEN
        IF(key3(1:3)=='OFF') ikpres=0
C----------------------------
C      buckling analysis
C----------------------------
       ELSEIF(key2=='BUCKL')THEN
        READ(key3,'(I2)')ibuckl
        IF (ibuckl==0) THEN
         WRITE(istdo,*) ' ** ERROR ** : KEYWORD /IMPL/BUCKL OBSOLETE ',
     .                  'USING /IMPL/BUCKL/1 OR /IMPL/BUCKL/2'
         GOTO 9990
        ENDIF
        ibuckl = ibuckl-1
        CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
        READ(iusc2,*) emin_b, emax_b, nbuck, msgl_b, maxset_b, shift_b
        IF (shift_b==zero) shift_b=em02
        shftbuck = shift_b
        IF (maxset_b==0) maxset_b=8
        bniter=300
        bincv=4
        bmaxncv=16
c        BMAXNCV=MAX(BINCV,BMAXNCV)
        bipri =msgl_b
        bisolv=1
C
       ELSEIF(key2(1:5)=='AUTOS')THEN
         IF(key3(1:3)=='OFF')THEN
           iautspc=0
         ELSEIF(key3(1:3)=='ALL')THEN
           iautspc=2
         ENDIF
C----------------------------
C      line_search option
C--------0=3, 1:energy 2: force --3:AUTO (old)------------------
       ELSEIF(key2(1:5)=='LSEAR')THEN
        IF(key3(1:3)=='OFF')THEN
         iline_s = 100
        ELSE
         READ(key3,'(I2)')iline_s
         CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
         READ(iusc2,*)nls_lim,ls_tol
        ENDIF
C----------------------------
C      projection for warped shell elements
C--------0=no proj but keep explicit part 1:doing -1 :no proj, neither for explicit---
       ELSEIF(key2(1:5)=='SHPOF')THEN
           ikproj=-1
C-----------become default-after-----
       ELSEIF(key2(1:5)=='SHPON')THEN
           ikproj=1
C----------------------------
C      OLD CONTROL OPTIONS
C----------------------------
       ELSEIF(key2(1:5)=='CONTR')THEN
         IF(key3(1:2)=='DT')THEN
          IF(key4(1:4)=='STOP')THEN
           CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
           READ(iusc2,*)dt_min,dt_max
          ELSE
           READ(key4,'(I2)')im
           idtc=im
           IF(im==1)THEN
            CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
            READ(iusc2,*)nl_dtp,scal_dtp,nl_dtn,scal_dtn
           ELSEIF(im==2)THEN
            CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
            READ(iusc2,*)nl_dtp,alen0,nl_dtn,scal_dtn,scal_dtp
           ENDIF
          ENDIF
         ELSEIF(key3(1:4)=='SHEL')THEN
C----------------------------
C         Fictif stifness of Mzz for shell
C----------------------------
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)kz_tol
         ELSEIF(key3(1:5)=='INTER')THEN
C----------------------------
C         stifness factor for interface
C----------------------------
          CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
          READ(iusc2,*)sk_int
         ENDIF
C----------------------------
C      hide options
C----------------------------
       ELSEIF(key2(1:5)=='PRTOL')THEN
        CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
        READ(iusc2,*)d_tol
       ELSEIF(key2(1:4)=='NEXP')THEN
        READ(key3,'(I5)')nexp
       ELSEIF(key2=='DEBUG')THEN
        impdeb=1
        CALL wriusc2(ikad(ikey)+k,1,key0(ikey))
        READ(iusc2,*)ndeb0,ndeb1
         IF(ndeb0/=0)ndeb0 = ndeb0 + 1
         ndeb1=max(ndeb0,ndeb1+1)
       ELSEIF(key2(1:3)=='DEL')THEN
         IF(key3(1:5)=='RBODY')THEN
           imp_rby=1
         ELSEIF(key3(1:5)=='INTER')THEN
           imp_int=1
         ENDIF
       ELSEIF(key2(1:5)=='ITRBY')THEN
C-------max iter for secnd dis calculation with finite rotation---       
        READ(key3,'(I3)')itrmax
       ELSEIF(key2(1:4)=='LRIG')THEN
        imp_lr = 1
       ELSE
         GOTO 9990
       ENDIF
       k=k+nbc
       IF(ikad(ikey)+k/=ikad(ikey+1))GO TO 1160
       IF (iparit/=0) THEN
        iparit=0
        ikg=ikg+5
       ENDIF
      ENDIF
C
      RETURN
C
 9990 CONTINUE
      CALL ancmsg(msgid=73,anmode=aninfo,
     .            c1=key0(ikey))
      CALL arret(0)
      END
!||====================================================================
!||    order_dtf   ../engine/source/input/freimpl.F
!||--- called by ------------------------------------------------------
!||    freimpl     ../engine/source/input/freimpl.F
!||====================================================================
      SUBROUTINE order_dtf(N,RC)
C----6---------------------------------------------------------------7---------8
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
C-----------------------------------------------
      INTEGER N 
      my_real
     .        rc(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,II,NN
      my_real
     .        s(n),smin
C
      IF (n==0) RETURN
C-----en ordre croisante-----
      nn =0
      DO i =1,n
       IF (rc(i)>zero) THEN
        nn = nn +1
        s(nn)= rc(i)
       ENDIF
      ENDDO
      n= nn
      DO i =1,n
       smin=s(i)
       ii=i
       DO j =i+1,n
        IF (s(j)<smin) THEN
         ii=j
         smin = s(j)
        ENDIF
       ENDDO
       IF (ii/=i) THEN
        smin =s(i)
        s(i)=s(ii)
        s(ii)=smin
       ENDIF
       rc(i) = s(i)
      ENDDO
C----6---------------------------------------------------------------7---------8
      RETURN
      END