Sources Esope | Cast3M

hbmini
C HBMINI    SOURCE    CB215821  23/01/25    21:15:22     11573          
*
      SUBROUTINE HBMINI(KTKAM,KTQ,KTFEX,KTPAS,KTLIAA,KTEMP,KTLIAB,KTPHI,
     &                  KCPR,KOCLFA,KOCLB1,KPARNUM,NHBM,NFFT,CHECK,ITER)
*    
C=======================================================================
C=                            HBMINI                                   =
C=                           --------                                  =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  A partir d'une approximation initiale, converge vers une solution  =
C=  des equations d'equilibre pour une frequence donnee.               =
C=   -->iterations hybrides NR - gradient conjuge                      =
C=                                                                     =          
C=======================================================================
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL CHECK
      CHARACTER*8 FLAG
*        
-INC PPARAM
-INC CCOPTIO
*-INC TMDYNC.INC
************************** debut TMDYNC.INC ****************************
 
*    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
*    TODO   : a extraire dans un include des que stabilise
*
*     Segment des variables generalisees:
*     -----------------------------------
      SEGMENT MTQ
         REAL*8 Q1(NT1)
         REAL*8 OMEG,XPARA
         REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
         REAL*8 dX(NT1), dw, dv
      ENDSEGMENT
*     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
*            Q1 = {q_0 q_c1 q_s1 ... q_sh}
*            avec q_i vecteur de dimension n ou n=nombre de modes
*     OMEG : frequence fondamentale de l'approximation
*     XPARA: parametre de continuation (par defaut la frequence) 
*            \in [PARINI,PARFIN]
*     RX   : matrice jacobienne = ZZ + dFnl/dX
*     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
*     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
*            lineaires et constantes
*     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
*
*
*     Segment contenant les matrices XK, XASM et XM:
*     ---------------------------------------------
      SEGMENT MTKAM
         REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
         REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
*         REAL*8 GAMFIN(NPC2,nl1)
      ENDSEGMENT
*     XK,XASM et XM : matrices de raideur, amortissement et masse
*     GAM et IGAM : matrices pour la FFT et son inverse
*     GAMFIN      : 
*
*     Segment des deformees modales:
*     ------------------------------
*     (idem DYNE)
      SEGMENT MTPHI
         INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
         INTEGER IAROTA(NSB)
         REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
      ENDSEGMENT
*
*     Segment descriptif des liaisons en base A:
*     ------------------------------------------
*     (idem DYNE)
      SEGMENT MTLIAA
         INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
         REAL*8  XPALA(NLIAA,NXPALA)
      ENDSEGMENT
*
*     Segment descriptif des liaisons en base B:
*     ------------------------------------------
*     (idem DYNE)
      SEGMENT MTLIAB
         INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
         REAL*8  XPALB(NLIAB,NXPALB)
         REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
      ENDSEGMENT
*
*     Segment representant les chargements exterieurs:
*     -----------------------------------------------
      SEGMENT MTFEX
         REAL*8  FEXA(NT1)
         REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
         INTEGER BAL
      ENDSEGMENT
*     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
*            Fourier et exprimes en base A
*     FEXPSM: chargement/deplacement statique lie aux modes negliges
*             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
*             compatibilite avec calcul des Fnl.
*     BAL  : indique s'il s'agit d'un chargement de type balourd
*            (cad proportionnel a OMEG**2)
*
*     Segment "local" pour DEVLFA:
*     ----------------------------
      SEGMENT LOCLFA
         REAL*8 FTEST(NA1,4)
      ENDSEGMENT
*      
*     Segment "local" pour DEVLB1:
*     ----------------------------
      SEGMENT LOCLB1
         REAL*8 FTEST2(NPLB,6)
      ENDSEGMENT
*      
*     Segment contenant les variables au cours d un pas de temps:
*     ----------------------------------------------------------
      SEGMENT MTPAS
         REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
         REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
         REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
         REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
         REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
         REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
      ENDSEGMENT
*     FTOTA/B/BA : forces sur base A, B et B projetees sur A
*     XPTB       : deplacement du point d'une liaison en base B
*     XVALA/B    : grandeurs de la liaison en base A/B a stocker
*     FEXB       : forces exterieures en base B (a priori uniquement 
*                  pour les moments appliques aux rotations rigides ?)
*     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
*                  maillage de contact dans certaines liaisons)
*     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
*                       base A/B (= contributions a dFnl/dX)
*
*
*     Segment des points de reference des modes (base A):
*     --------------------------------------------------
      SEGMENT MPREF
         INTEGER IPOREF(NPREF)
      ENDSEGMENT
*
*     Segment des points en base B:
*     -----------------------------
      SEGMENT NCPR(NBPTS)
*     NCRP(#global) = #local dans XPTB (1er indice)      
*
*     Segment des parametres numeriques pour la continuation: 
*     ------------------------------------------------------
      SEGMENT PARNUM
         CHARACTER*4 TYPS
         REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
         REAL*8   PARINI,PARFIN
         INTEGER  ITERMAX,NBPAS
         LOGICAL  JANAL
      ENDSEGMENT
*
*     Segment des resultats:
*     ---------------------
      SEGMENT  PSORT
        REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
        REAL*8 VSAVE(NPAS)
        LOGICAL ZSAVE(NPAS)
        CHARACTER*2 TYPBIF(NBIFU)
        REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
        REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
        INTEGER CBIF
      ENDSEGMENT
*     QSAVE(i,j) = Q harmonique i au pas j
*     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
*     ZSAVE(j)   = stabilite au j-eme pas
*     LSAVE(1,j) : partie reelle de l'exposant de Floquet
*     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
*     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
*     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
*           WBIF2 : partie imaginaire de l'exposant de Floquet
*     QPSIR,QPSII : vecteur propre au point de bifurcation
 
*     Segment des tableaux de travail:
*     -------------------------------
      SEGMENT  MTEMP
         REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
         REAL*8 T02(NT1+2), TP2(NT1+2)
         INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
         REAL*8 res
         REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
         REAL*8 QOLD(NT1),OMEGOLD
         REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
         REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
      ENDSEGMENT
*     Jacobiennes augmentees 
*     Ja : [ RX  Rw ; dX dw]
*     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
 
*      SEGMENT  NNNN
*         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
*      ENDSEGMENT
 
*************************** fin TMDYNC.INC *****************************
 
*-----------------------------------------------------------------------
*       Le contenu des segments est donne en entree a NEWT: 
*
*       Matrices XK,XASM,XM,GAM,IGAM      
        MTKAM=KTKAM
*        
*       Variables generalisees: coefficients de Fourier et frequence
        MTQ=KTQ
*
*       Forces externes
        MTFEX=KTFEX
*        
*       Variables sur un pas de temps
        MTPAS=KTPAS 
*        
*       Liaisons sur base A
        MTLIAA = KTLIAA
*       Liaisons sur base B
        MTLIAB = KTLIAB
*
*       Truc local base A
        LOCLFA = KOCLFA
*       Truc local base B
        LOCLB1 = KOCLB1
*        
*       Deformees modales
        MTPHI = KTPHI 
*       
*       Parametres numeriques
        PARNUM = KPARNUM
* 
        MTEMP=KTEMP
*       Nombre de modes      
        NDDL = XM(/1)   
        NDDL2 = XM(/2)  
*       Nombre d'inconnues
        NT = Q1(/1) 
 
*----- Iterations de Newton Raphson  -----------------------------------
 
*      -systeme NON-Autonome (= force)
       IF (TYPS.EQ.'FORC') THEN
        CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
     &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'IF',ITER)
 
*      -systeme AUTOnome (sans forces exterieures)
       ELSEIF (TYPS.EQ.'AUTO') THEN
        CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
     &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'IA',ITER)
 
*      -systeme Hamiltonien (modes non lineaires) 
       ELSEIF (TYPS.EQ.'NNM') THEN
        CALL HBMNEWT(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTEMP,PARNUM,
     &       MTLIAA,MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,'IN',ITER)
       ELSE
         WRITE(IOIMP,*) 'DYNC : Type TYPS de probleme inconnu !'
         CALL ERREUR(5)
         RETURN
 
       ENDIF
 
 
*     Sortie
      KTQ = MTQ
      KTEMP=MTEMP
 
      END